Atmospheric boundary-layer (ABL) flows over complex terrain have been the focus of active research, given their impact on weather and climate variability. Surface complexity is understood in a broad sense and includes variation in roughness properties, inclination of the underlying surface, presence of heterogeneous forcing mechanisms (e.g., buoyancy, humidity), to name but a few. Most assumptions of classical boundary-layer similarity theory do not hold under such conditions, complicating matters from both a measurement and modeling perspective. Here, a combination of analytical and numerical approaches are used to address two among the most relevant problems: turbulent slope flows, and ABL flows over multi-scale rough surfaces. The first part of the thesis focuses on slope flows: the building blocks of local weather in mountainous regions. To understand the system conceptually, a closed-form analytic solution to the Prandtl slope flow model is first derived, prescribing transfer coefficients in accordance to the O'Brien K-theory model. Profiles are characterized by stark variations in both phase and amplitude of extrema compared to the classic constant-K and a more recent solution, valid within the Wentzel-Kramers-Brillouin theory, shedding new light on this long-standing geophysical problem. In addition, direct numerical simulation is used to study the turbulent structure of anabatic and katabatic flows, and to describe the sensitivity of the solution to variations in the parameter space, within the conceptual framework of the Prandtl model. Variations in the sloping angle from the vertical wall setup are shown to induce a progressive departure of averaged profiles between the two flow regimes, ultimately resulting in stark differences at gentle sloping angles. The thermodynamical mechanisms responsible for sustaining mean and turbulent kinetic energy are used to further distinguish between flow regimes, and to propose a qualitative partition of the boundary layer in slope flows. The DNS setup is additionally adopted to identify coherent structures in katabatic flows over steep slopes. Coherent motions are responsible for the maintenance of turbulence in the ABL, hence their characterization is of fundamental importance toward a better understanding of boundary-layer dynamics. Packets of hairpins are found to connect in the streamwise direction to form large-scale motions (LSMs). For the lower sloping angles that are considered, it is then shown how LSMs further align to form very-large-scale motions (VLSMs). LSMs and VLSMs are found to be dominant contributors to streamwise momentum variance and turbulent momentum transfer in the above-jet regions. Next, drag properties of fractal-like sea ice surface morphologies are examined within the large-eddy simulation framework. The effects of large-scale surface features on wind flow are accounted for by an immersed boundary method. Conversely, the drag forces caused by subgrid-scale features are modeled through a novel dynamic roughness approach, in which the hydrodynamic roughness length parameter is determined using the first-principles based constraint that the total momentum flux (drag) must be independent of the grid-filter scale. This approach leads to accurate flow predictions, and provides an estimate of the otherwise unknown roughness parameter for sea ice surfaces, of use in climate, weather prediction and scalar transport models to evaluate the hydrodynamic roughness length.