Effective Roughness and Displaced Mean Flow over Complex Terrain

Via analysis of velocity and stress fields from Reynolds-Averaged Navier–Stokes simulations over diverse complex terrains spanning several continents, in neutral conditions we find displaced areal-mean logarithmic wind speed profiles. The corresponding effective roughness length (z0,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_\text {0,eff}$$\end{document}), friction velocity (u∗,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u _{*\text {,eff}}$$\end{document}), and displacement height (deff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_\text {eff}$$\end{document}) characterise the drag exerted by the terrain. Simulations and spectral analyses reveal that the terrain statistics—and consequently deff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_\text {eff}$$\end{document}, u∗,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u _{*\text {,eff}}$$\end{document} and z0,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_\text {0,eff}$$\end{document}—can change significantly with flow direction, including flow in opposite directions. Previous studies over scaled or simulated fractal surfaces reported z0,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_\text {0,eff}$$\end{document} to depend on the standard deviation of terrain elevation (σh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _h$$\end{document}), but over real terrains we find z0,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_\text {0,eff}$$\end{document} varies with standard deviation of terrain slopes (σΔh/Δx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{\Delta h/\Delta x}$$\end{document}). Terrain spectra show the dominant scales contributing to σΔh/Δx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{\Delta h/\Delta x}$$\end{document} vary from ∼\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document}1–10 km, with power-law behaviour over smaller scales corresponding to fractal terrain used in earlier works. The dependence of z0,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_\text {0,eff}$$\end{document} on σΔh/Δx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{\Delta h/\Delta x}$$\end{document} is consistent with fractal terrain having σΔh/Δx∝σh\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{\Delta h/\Delta x} \propto \sigma _h$$\end{document}, as well as classic theory for individual hills. We obtain relationships for z0,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_\text {0,eff}$$\end{document}, deff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_\text {eff}$$\end{document}, and u∗,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u _{*\text {,eff}}$$\end{document} in terms of σΔh/Δx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _{\Delta h/\Delta x}$$\end{document}, finding that deff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d_\text {eff}$$\end{document} acts as a characteristic length scale within z0,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_\text {0,eff}$$\end{document}. Considering flow in opposite directions, use of upslope statistics did not improve z0,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_\text {0,eff}$$\end{document} predictions; sheltering effects likely require more sophisticated treatment. Our findings impact practical applications and research, including micrometeorological flow, computational fluid dynamics, atmospheric model coupling, and mesoscale and climate modelling. We discuss limitations of the z0,eff\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_\text {0,eff}$$\end{document} formulations developed herein, and provide recommendations for practical use.


Motivation and Outline
In general we are interested in understanding flow in the atmospheric boundary layer (ABL) over hilly terrain, and simulation of such. This includes use of computational fluid dynamics (CFD); because CFD is used both in research and various applications, we have some focus on Reynolds-averaged Navier-Stokes (RANS) solvers, which give mean flow solutions while also modelling turbulence. We also wish to better understand some simpler (one-dimensional) aspects, such as vertical profiles of stress and mean wind speed; typical measurements capture these, and they are also present in parametrisations of larger-scale models such as mesoscale, weather, and even climate models. For various applications one wishes to use mesoscale-model output, but needs to know the degree to which the model does or does not treat variations in terrain elevation; we wish to advance quantification of the unresolved (subgrid) terrain-induced drag, for a given model resolution and terrain. Further, for driving RANS (or large-eddy) simulations over complex terrain-whether via mesoscale output or prescription of inflow conditions-we wish to account for the terrain-induced drag experienced by the flow in a given simulation, as well as minimise adjustment of simulated inflow into the most finely resolved computational domains. Yet another motivation arises when considering applications where a large-scale roughness is required in flow parametrisation; for example, use of the geostrophic drag law in wind energy applications requires a geostrophic-scale roughness length.
More directly, we wish to test the effective roughness concept over actual terrain (elevation maps) spanning a range of complexities from locations around the world. Further, we aim to find not only the degree to which complex terrain increases drag on the flow (via a roughness length), but how the effect of the terrain may systematically vary with height depending on statistical terrain characteristics; this also includes checking where and whether a surfacelayer (logarithmic) type of wind profile exists above the surface.

Outline
First we introduce the concept of effective roughness and review previous work on the topic; this ranges from analytical estimates, based on flow physics for single-scale hills, to empirical expressions for multiscale terrain. Section 2 looks into spectra and statistics of terrain elevation, with particular attention given to terrain-slope statistics. Section 3 presents analysis of mean flow from computational fluid dynamics simulations, including determination of quantities describing flow over complex terrain (effective friction velocity, diplacement height, and roughness) for five different cases. Section 4 compares predictions using previous formulations with effective roughnesses diagnosed from the flow simulations; it further presents new forms for effective roughness and analysis of their predictions. Section 5 distills the utility of the new forms developed herein, and failure of previous forms over multiscale terrain; it also discusses combination of slope-associated roughness with background (aerodynamic) roughness, limitations of our investigation, and recommendations for application of the new forms. We then summarise the main points of this work, and share our outlook on its continuation and future research.

Effective Roughness Concept and Theory
The concept of an effective roughness is not new, particularly if one generalises to include roughness over obstacle arrays, forests, and urban canopies-let alone over complex terrain. Fiedler and Panofsky (1972) first described it as "that roughness length which homogeneous terrain would have in order to produce the correct space-average downward flux of momentum near the ground, with a given wind near the ground." In the current paper we further develop this concept without the constraint of being near the ground, further including the terraininduced displacement of mean flow analogous to that due to forests, as shown the sections below.
To deal with the mean effect of hills, Smith and Carson (1977) suggested an effective roughness relation, for use in single-layer calculations over "variable" (non-uniform) terrain at resolutions coarser than 10 km; h was the average height range in a given area, and L c was the average distance between peaks or ridges. For the relatively flat terrain of the Netherlands Wieringa (1986) used (1), but set h to be the maximum elevation difference within 5×5 km squares, picking L c =5 km, to obtain an 'elevation variability roughness length'. Noting a number of observational studies over hilly terrain and numerical simulations over hills, Wood and Mason (1993) noted that areally averaged flow follows the generic law of the wall; they thus defined effective roughness via the displaced logarithmic profile: where U is the areally averaged wind speed profile, and u 2 * 0,eff is the areally averaged surface force (per area and normalised by density). 1 They proposed expressions for the pressure force due to turbulent flow over hills, employing such to derive an approximate relation for the effective roughness. Adding the pressure-force and surface-stress contributions to total momentum transfer, Wood and Mason (1993) gave: for hills of height h and streamwise scale λ x ; C a was a constant proportional to the product of bulk slope (h/λ x ) and ratio of frontal hill area to domain size ("packing fraction"), with Z m being a pressure scale height. 2 The effective roughness in (3) is non-monotonic in hill slope and packing fraction, though for each it has an increasing behaviour for moderate values (after first dropping at small values); it converges to z 0 for flat terrain. Though Finnigan (2002) called the Wood and Mason (1993) formulation (3) an "ad hoc derivation," he admitted that 1 Note Taylor et al. (1987) also proposed a form like (2) for flow over flat surfaces with heterogeneous roughnesses (but with d = 0). They showed the z 0,eff concept applies, for heights above ∼200 m, with areal averaging of ln z 0 sufficing to get z 0,eff ; the latter corresponds to a geometric mean of z 0 (e.g. Kelly and Jørgensen 2017). 2 In Wood and Mason (1993) the pressure scale height Z m was simply prescribed as the maximum of hill height or middle-layer height h m , with h m = (λ x /4) ln −1/2 (h m /z 0 ) found numerically (Belcher et al. 1993).
We note one can directly find the analytical solution h m = λ x 8W λ 2 x /(8z 2 0 ) −1/2 via Kelly (2021), where W[x] is the Lambert-W function. it performs reasonably well for slopes less than ∼0.3; Hignett and Hopwood (1994) found a crude version of it to estimate z 0,eff within a factor of ∼3 for sites over moderately hilly terrain. 3 Finnigan (2002) further wrote that the terrain-perturbed region of mean flow extends up to z ∼ L (where L is the horizontal lengthscale of the hills)-and that unlike canopies, flow patterns around hills are not commonly affected by their neighbours. However, here we note that (3) does not account for collections of hills of different slopes and sizes, as found in most complex terrain-i.e. over terrain with a broad-band spectrum of elevations; perhaps more importantly, it does not deal with the interacting effects of shears and stresses induced by neighbouring hills, which may become significant where there are steep slopes.
Addressing the multiscale aspect of terrain, Beljaars et al. (2004) built upon the Wood and Mason (1993) formulation, and additionally upon Wood et al. (2001) for stress induced by sinusoidal hills, to derive height-dependent stress due to orography having a broad-band spectrum. However, because their work was aimed at parametrisation of unresolved terraininduced stress in mesoscale and regional models (having resolutions on the order of ∼10 km or coarser), Beljaars et al. (2004) prescribed a fixed form for terrain spectra, and did not include an effective roughness.
Towards characterization of flow over complex multiscale surfaces, particularly at smaller scales, large-eddy simulation (LES) has been used to find the effective roughness over surfaces having elevation spectra that follow a power law; examples are Wan and Porte-Agel (2011) and Anderson and Meneveau (2011). The latter works, along with the review of Flack and Schultz (2010) for engineering flows based on wind-tunnel observations, gave some evidence for a connection between the standard deviation of surface elevation (σ h ) and the effective roughness; this is elucidated more in the next subsection. Recently Finnigan et al. (2020) summarised progress in simulating flow over complex topography, stating that future parametrisations for drag due to unresolved terrain "would use the statistics of the topography to inform the statistics of the space-time structure of the airflow"; this is precisely what the current work has aimed to do, as previously seen in presentation of preliminary results (Kelly et al. 2019a). Here we investigate the connection between areal statistics of terrain elevation and the flow above, particularly considering complex terrain having broad-band elevation spectra.

Previous Attempts to Find Effective Roughness over Multiscale Terrain
Several studies have worked towards finding an effective roughness to account for the (pressure) drag induced by surface height variations occuring at horizontal scales smaller than some threshold; the latter can correspond to a flow model resolution, to allow parametrisation of roughness due to unresolved terrain. One commonality they share is the assumption of σ h as the characteristic length scale for z 0,eff ; this was first proposed by Stone and Dugundji (1965) for broad geophysical application. Flack and Schultz (2010) considered numerous engineering forms for frictional drag due to rough surfaces, and subsequently derived an empirical relationship for roughness over irregular (random) surfaces. Their roughness estimate depended primarily on the standard deviation of surface elevation variations (σ h ), as well as incorporating the skewness Sk h ; it is expressible as: where they reported a proportionality constant equivalent to c fs = 0.148. Through large-eddy simulations over a series of synthesised multiscale topographies, Wan and Porte-Agel (2011) investigated the postulate that the effective roughness length scales with the standard deviation of unresolved surface heights, σ h,sgs . To account for the lower limit of roughness in the case of zero subgrid height variations, they suggested a form expressible as: Examining flows over synthesised terrain whose spectra had roughly the same spectral form (following a power-law with exponent β −2), they found (5) to hold, but with the coefficient C increasing with z 0 . Anderson and Meneveau (2011) also used LES to examine the response of flow over flat rough surfaces whose roughness-element heights were synthesised to possess power-law spectra (i.e. fractal terrain), considering different surfaces over a range of spectral exponents β. Over statistically homogeneous surfaces they found the hydrodynamic roughness length (effective z 0 within the log-law) to follow z 0, = ασ h, , where the corresponds to the filter scale of the LES and α is a coefficient depending on the power-law exponent (β) of the terrain spectrum. To account for flows over intrinsically rough surfaces 4 with (background) roughness z 0 , for "numerical convenience" they suggested an effective roughness of the form: This expression appears to rectify the problem of C not being constant in the form (5) of Wan and Porte-Agel (2011). However, Anderson and Meneveau (2011) did not include a background z 0 in their analysis; this is because they state consideration of cases (LES) only where ασ h| z 0 . But they did find α to be monotonic in β, and showed α was invariant to LES resolution and filter scale ( ) for the fractal surfaces and flows simulated. From Fig. 8 of Anderson and Meneveau (2011), which plots α diagnosed over different surfaces with spectral exponents in the range −3 β −1.4, the dependence of α on β can be inferred as: We note their assumption ασ h| z 0 means that for surfaces with more negative spectral exponents (steeper roll-off of the spectrum), such as their flattest case having β = −3, then σ h /z 0 needs to be at least ∼ 10 6 .
Another form for effective roughness was developed by Kelly (2016), but through a different means, stemming from simplified treatment of flow over roughness changes. To account for flow perturbations induced by roughness changes within a quasi-linearised spectral model, Astrup and Larsen (1999) described friction velocity perturbations in Fourier space using a characteristic length scale of (z 0 2 s ) 1/3 , where s is the reciprocal of streamwise wavenumber. Adapting this form for an efffective roughness to be used in estimation of uncertainty in shearbased vertical extrapolation of mean winds, Kelly (2016) suggested z 0,eff = [z 0 (σ h +z 0 ) 2 ] 1/3 , taking the characteristic length scale as σ h .
We note that the form of Astrup and Larsen (1999) can be seen to follow from Jensen et al. (1984), who described the depth of the turbulence-dominated inner layer of (unseparated) flow over rough hills, where the characteristic hill size (outer scale) L h is used instead of s . Recalling also that both the Jensen et al. (1984) and Astrup and Larsen (1999) forms describe equilibrium-layer depths, it is sensible to adapt such for a statistical or mean description of multiscale terrain effects. Over complex terrain the characteristic length scale due to terrain elevation variations has been presumed to be σ h in the studies named thus far; this was also done for the Kelly (2016) form, which can be written more generally as: with the coefficient c m requiring empirical determination. Looking at the forms (5), (6), and (8) from Wan and Porte-Agel (2011), Anderson and Meneveau (2011), and Kelly (2016), respectively, we see that they can be generalised to the σ h -dependent form the coefficient c can depend on the shape of terrain spectrum, as in Anderson and Meneveau (2011).

Terrain Analysis
Here we start by considering the statistical character of complex terrain, in terms of surface elevation variations. We examine a number of sites, each possessing different geological and geometric character (created via various physical mechanisms over different time scales). These sites include: a hilly area in north-east China with modest slopes, a Portuguese area of sharp orography with steep slopes, a South African area with a range of hills, and two Norwegian areas dominated by mountains and valleys. Again, one aim is to develop robust, universal statistical descriptions of non-uniform terrain and the (mean) flow above it, so we have chosen sites with different spectral and statistical character, as will be seen in the next subsections; however, the number of sites is limited by computational resources needed for the corresponding flow simulations. The elevation maps have resolutions which range from ∼ 20-100 m, according to the spacing of contours/points in their digital formats (however, as we will see below, the actual resolution of the maps can be significantly cruder). The resolution of the terrain data analyzed herein is sampled at a constant 56 m in the streamwise direction parallel to the terrain. To calculate terrain statistics we use linear transects across the domains, which span roughly 40 km.

Spectra of Terrain Elevation and Slopes
Towards statistical characterization of terrain variations in a general sense, we first examine their spectra. One-dimensional power spectra of the terrain elevation h, defined as S hh (k 1 ) ≡ F [h(x)], where F denotes the Fourier transform and k 1 is the wavenumber corresponding to a given x direction (later corresponding also to the mean flow direction), were calculated via fast Fourier transform. This was done in 12 different directions (every 30 • ), along 21 transects in each direction over a lateral span of 4 km (i.e., every 200 m). An average over the 21 lateral spectra was also taken to provide a single 'mean' spectrum in each direction; these are shown in the upper plots of Fig. 1 for the least complex (old hills in north-east China, at left) and most complex (Portuguese mountains, at right) areas considered in this study. Due to symmetry, only six unique spectra are shown in each plot, along with the omnidirectional spectrum (average of the six directional spectra). To mitigate the aliasing effect, which is caused by sampling with a constant spacing x that is less than the resolution in the outer  part of the domain 5 we have applied a fourth-order low-pass Butterworth filter with cutoff wavenumber of k cut = 0.5k Nyquist (i.e., 2π/k cut = 4 x) to the spectra displayed. Note the lowest wavenumber portion of the terrain elevation spectra represents variations at horizontal scales approaching the domain size; thus the spectra and the total variance σ 2 h can be influenced by the choice of domain size. The latter is particularly evident in Fig. 1 due to the highest spectral magnitudes being found at the lowest wavenumbers, which is not uncommon for complex terrain. However, the spectra of terrain slope are generally not affected by the domain size 6 ; these are shown in the bottom of Fig. 1 below S hh (k 1 ). Subsequently the domain size has less impact on statistics of dh/dx than on statistics of h. This is seen by noting that the power spectrum of terrain slope is simply the product of wavenumber squared and the terrain spectrum, i.e., for a given x-direction. 5 The uniform sampling was identical for the terrain elevation and all analyzed flow fields; this was done to deal practically with data from a large number of simulations, which have a horizontal resolution ranging from 20 m in the central area of the domain to more than 70 m at r > 17 km from the center (see Sect. 3 for more detail). 6 Domain size will not affect terrain-slope statistics, unless such a small domain is chosen that the statistics become very uncertain or the domain becomes unrepresentative.
We can see that a power law is generally observable from scales of several hundred metres up to several kilometres; this is shown by the terrain-slope power spectra k 2 S hh (k) plots as well as the terrain spectra in Fig. 1. The plots of S hh further show that a different behaviour occurs at wavenumbers smaller than (scales larger than) the peak of k 2 S hh (k); such a peak (k peak ) implies a characteristic length scale 2π/k peak for the terrain along a given direction. The peak in k 2 S hh (k) corresponds to the change in power-law exponent (spectral slope) of S hh (k) noted by Beljaars et al. (2004); they parametrised a change in spectral exponent to occur at k = 0.003 m −1 , i.e. at a horizontal scale of roughly 2 km. We find that the peak of the one-dimensional terrain-slope spectrum varies around this value from one site and direction to another, generally in the range ∼1-10 km. The Beljaars et al. (2004) prescription for terrain spectra is shown by the cyan line in Fig. 1 for the omnidirectional mean. Its quite practical parametrisation is anchored to the observed value of S hh near k = 0.00035 m −1 , resulting in spectra which basically follow observed omnidirectional spectra within a factor of ∼ 3. However, looking at the terrain-slope spectra one sees that the shape of k 2 S hh (k)-including the peak wavenumber-varies sufficiently from one place or even direction to another; the robust Beljaars et al. (2004) form can (counterintuitively) overpredict slope variance for complex sites and underpredict it for simpler ones by a factor of ∼ 3, due to its anchor point being at a relatively low wavenumber and its prescription of a low-k spectral exponent of −1.9 differing from observations.
The description of smaller-scale terrain features by single power-law spectra is not new. Numerous works simply assume this for the entire landscape in a given area; i.e., for fractal terrain the surface elevation spectrum is taken to follow the simple form: Following Mandelbrot (1989), for fractal terrain its fractal dimension can be described by At any rate, like various studies in previous decades (Anderson and Meneveau 2011) and our earlier work considering various sites around the world for the Global Wind Atlas (Badger et al. 2015), here we find a power law with β ranging between approximately −4.3 and −1.2; for the most complex (Aveiro) area the mean β over all directions was equal to the Beljaars et al. (2004) prescription of β = −2.8 for k > 0.003 m −1 , with less complex areas having smaller (more negative) β. We note this range of β only corresponds to scales smaller than the spectral peak of the slope spectrum (i.e., at k > k peak ). For larger scales (k < k peak ) one sees either a different power law, or yet more complicated behaviour, as notable from the terrain-slope spectra shown in Fig. 1. In other words, terrain is not mono-fractal in nature-or at least not having the same fractal dimension across all scales. 7

Terrain Slopes and Finite Differencing
In Fig. 1 we showed spectra that had an anti-aliasing filter applied to treat the effects of sampling a domain having non-constant horizontal spacing. In practice we can also mitigate the aliasing issue by calculating the spectra of ∂h/∂ x using F ( h/ x), where x is optimally the map resolution (grid spacing). A first-order finite difference in physical space corresponds to combination of low-pass filter sinc 2 (k 1 x/2) and multiplying by k 2 in wavenumber space, i.e., F ( h/ x) = F (∂h/∂ x)sinc 2 (k 1 x/2). For the least complex terrain we studied (north-east China case in left-hand plots of Fig. 1), without the anti-aliasing filter a spurious upturn at high k 1 occurs not only in F (∂h/∂ x), but even in F (h) (not shown); such noise is more easily seen in terrain-slope spectra, because the factor of k 2 amplifies such. 8 Due to this k 2 amplification, map processing artifacts may be substantial enough to also cause an upturn in the spectrum of finite-differenced slope, S h/ x at the finest scales; this is the case for the north-east China region when considering the entire map, due to the subsampling which occurs for the parts of the map at distances increasing beyond ∼ 8 km from the domain center. 9 This is seen in Fig. 2, which displays S h/ x for the same cases shown in Fig. 1, but without an anti-aliasing filter applied. The top plots in Fig. 2 present finite-differenced terrain slope spectra, with the upper-left plot (north-east China region) exhibiting a high-wavenumber upturn.
The spectra of finite-differenced slopes from Aveiro (upper-right) do not suffer from the high-wavenumber upturn; unlike all other regions considered, its finest-resolution area was larger (the central 17×17 km), so that sampling issues only occurred over a small part near the outer-most domain edges.
To further demonstrate terrain-slope spectra via finite-differencing, the bottom plots of Fig. 2 display S h/ x calculated over the finest resolution central 1/3 of the maps (r < 12 km), where there is less undersampling. These plots in Fig. 2 are also made in variance-preserving coordinates, specifically k 1 S h/ x (k 1 ) versus ln k 1 , meaning they directly show the terrain scales which contribute to σ h/ x . For hilly terrain we find the majority of variance in slopes is contributed by variations at horizontal scales ranging from ∼ 0.5 km to 3 km, as seen by the peaks in the bottom plots of Fig. 2 for the least and most complex cases analyzed in this work.
Employment of anti-aliasing filters or h/ x to damp the high-wavenumber noise is defensible from a physical standpoint, as the geological processes which create the terrain on average cause an observed decay in S hh (k) with wavenumber, particularly at the smallest scales k 0.1 m −1 ; erosion processes, such as those driven by the wind, help to wear down sharp edges which would otherwise give contributions to the spectrum at these scales (e.g. Brown and Scholz 1985;Sulebak 1999). The use of low-order finite difference is justifiable, as long as the map resolution is fine enough relative to the variance-containing part of the spectrum; this is evident from Figs. 1, 2.
For robustness and universal applicability, hereafter we primarily consider statistics of h/ x for our analyses, because the anti-aliasing filter needed could depend on both the range of resolutions for a given terrain map and any processing that may have been done on the map. Further, calculation of such statistics are more accessible to readers, and advanced filter techniques are not the focus of this work. As we will see in later sections, the highwavenumber noise, discussed above, does not prevent relation of terrain-slope statistics to large-scale effective roughness-rather, it facilitates the latter.
If a spectral exponent can be estimated for k > k peak , this can be used to extend the spectrum and avoid the effects of processing-induced noise, provided that care is taken to ensure that the exponent is not affected by the filtering caused by finite-differencing; this is pertinent considering the left-hand plot in Fig. 2, for example. For the purposes of this work-notably considering effective roughness as seen from hundreds of metres (or more) from the surface-the smallest-scale effects due to the highest wavenumber slope variations act more locally and do not contribute substantially to the mean flow at such heights. 10

Terrain Asymmetry: Upslopes
In the above presentation of one-dimensional terrain spectra, for the 12 directions considered, only 6 unique spectra were found for each area, as shown in the figures above; this is because the Fourier spectra of elevation are invariant to coordinate reversal. But the geological processes that form terrain, including wind-driven erosion (for example), are not expected to be symmetric: for a one-dimensional transect across the terrain, the upslopes and their statistics differ when viewed from opposite directions.
Indeed one can calculate upslope spectra, e.g., by setting the negative slopes to zero. Doing so, we obtained different spectra in all 12 directions for each terrain map (not shown); the upslope spectra have the same peaks as the slope spectra shown in Fig. 2, but with more variation of spectral shapes. Across all sites and directions we find the standard deviation of upslope in a given direction can be as much as 8/5 times larger than in the opposite direction, with the ratio of mean up/downslopes reaching 4/3. Collecting all cases together we note both these ratios are centered around 1 (with the mean slope ratio following from neglect of domain-mean slope, see Sect. 3.1.2). Considering the slope asymmetry found, and that the upslope portion of hills tend to dictate more of the form drag, upslope terrain statistics may be of use in predicting effective roughness and other complex terrain effects; this is discussed in the following sections.

Terrain Statistics and Universal Statistical Features
One exploitable aspect, which can find use in predictive applications such as ours and potentially be described in a universal way, is the statistics of terrain slopes; most notably, this includes the probability distribution function (PDF) of slopes and subsequently low-order moments. In particular one sees a pattern emerge upon comparing the PDFs of slopes, in contrast with terrain elevation distributions. To show this, the probability density function of h/ x is displayed in Fig. 3 for the Portuguese and Chinese areas analysed, for all 12 directions and across 21 different transects over 4 km (separated by 200 m) for each direction; i.e., each line corresponds to a histogram over 21 transects, covering an area of 4 km×34 km. Only positive slopes are shown, because the twelve directions span 360 • ; for a given direction h/ x < 0 simply appears as h/ x > 0 from the opposite direction.
One can see a collapse of P( h/ x) for the most common slopes; we note that the shape of the PDF in this range depends on the power-law exponent β defined in the highwavenumber (sub-mesoscale) range. The magnitude of what is common (e.g., the slope where the cumulative density function exceeds 0.3) depends on the complexity of the area, with more complex sites having larger slopes occuring more frequently. Again, as with Figs. 1 and 2, only the least and most complex sites are shown; the other three sites have plots which fall between these two, depending on the complexity.
With regards to the PDFs of terrain slopes shown in Fig. 3 and asymmetry of opposing directions, we note that the mean upslope varied from 1/40 to 1/8, while the standard deviation of upslope (σ h + / x ) ranged from roughly 0.035 to 0.21.

Flow Simulations
Simulations were made for 36 directions, i.e. every 10 • , for each site; however, due to the terrain and flow statistics not changing significantly over such a small angular increment, we analyse the flow in every third 10 • sector, as done for the spectra in Sect. 2.1.
Though we present results for five cases, several more had been considered, including the Columbia Gorge region in the north-west contiguous United States and a Pakistani site in the outer Himalayas; however, achieving convergence in all directions for these simulations required excessive resources, so these sites were not included. 11

Flow Model and Mesh
The RANS solver Ellipsys3D (Sørensen 1995) was used to simulate mean flow in 36 directions over each of the various areas considered here, calculating steady-state solutions (as opposed to unsteady solutions from so-called 'URANS'). Ellipsys3D is a parallelised multiblock finite-volume code, employing the two-equation k-turbulence model and the SIMPLE pressure-solving algorithm (Patankar and Spalding 1972) 12 . The circular computational surface mesh, suitable for calculation of any wind direction, was created with the hyperbolic generator HypGrid2d (Sørensen 1998). The three-dimensional hexagonal volume mesh is grown away from the surface (thus vertical faces can be treated) using the HypGrid3d mesh generator, with nearly orthogonal cells having negligible skewness. The surface boundary is treated via shear stress prescribed with a high-Reynolds number assumption, consistent with both k-closure and surface-layer theory (Cavar et al. 2016).
The simulations were conducted with neither buoyancy nor Coriolis force, and were consequently Reynolds-number independent (Chew et al. 2018;van der Laan et al. 2020). So the solutions do not depend on the geostrophic wind or U in (z) or u * in , and only one simulation is needed per direction per map.

Model Boundary Conditions
Bottom conditions: The bottom boundary has a uniform roughness z 0 = 9 cm, essentially equal to the roughness z 0,in characterising the logarithmic inflow wind profile. The height of the bottom surface is dictated by an elevation map for each site/area considered. For manageable calculations the RANS mesh resolution is 20 m only in the innermost 4 × 4 km of the domain, and gradually stretches outward, reaching 25-40 m within a 7 × 7 km extent around domain center. Beyond that, the map is gradually smoothed, and the computational grid spacing likewise increases, with both approaching a resolution of 80-100 m for r 10 km from the center and progressing to ∼ 200 m for r ∼ 20 km. The elevation is smoothed to a uniform height towards the outermost part of the circular domain bottom. The flat outer zone is designed to both ensure the incoming log-law profile from inlet on the upwind side, 11 We note that due to its high complexity and difficulty converging, a map with inner-domain extent of 17×17 km (constant app. 20 m resolution) and total radial extent of 40 km was used in the Aveiro-Viseu case. 12 Here we have used the k-turbulence model, with ABL-appropriate values of its constants (C μ , C 1 , C 2 ) as reported in Bechmann (2016) and used in validation studies of e.g. Troen et al. (2014) and Troen and Hansen (2015). A first-order upwind discretization scheme was used for the k and transport equations, and the QUICK scheme was employed for the discretization of the RANS (velocity) equations (Leonard 1979).

Fig. 4
Depiction of flow domain with terrain, inflow direction, outflow region (shaded), and coordinate system; analysis domain denoted by black circle. Case shown corresponds to Norwegian 'NO52' area and that the flow can settle to a steady-state near the outflow boundary on the downwind side with proper convergence (flow variable residuals less than 10 −6 ). Top boundary conditions: The top boundary is located at a constant height of about 15 km above the mean surface and has a constant wind velocity boundary condition imposed. It is sufficiently high to ensure (a) no damping of any variable is necessary, and (b) no artificial speed-up effects arise due to height differences between the inner domain's hilly surface and the uniformly flat surface in the horizontal outer zone. A constant-velocity condition was also chosen at the top boundary to enhance the numerical stability of the computations, ensuring that only the streamwise velocity component has a constant non-zero value at domain top. Thus we are considering RANS in neutral conditions without a capping inversion, where the domain height z top is much greater than the dominant terrain scale (z top σ h z 0 ) such that there are no 'top-down' effects. Inlet/outlet conditions: The RANS simulations have an inlet condition described by a logarithmic wind profile, covering the incoming part of the cylindrical surface of the grid domain (which has a uniform flat surface). A zero-gradient boundary condition was applied to define the outflow part of the cylindrical domain boundary on the downwind side; the outflow zone width is approximately 45 • , which gives optimal trade-off between the outlet zone needed for the finite-volume RANS solver, and ability to obtain numerically stable convergence. An example flow domain is shown in Fig. 4.

Consideration of Flow-Field Analysis Domain
Here we have avoided the effect of transition from flat inflow far-field to hilly surface (and more broadly the momentum flux footprint aspect), by conservatively limiting the calculation of flow statistics to the downwind two-thirds of each (∼40 km long) simulated domain; even for the least complex terrain cases, such a buffer (wider than 12 km) eliminated the impact of such artifacts on finding the effective roughnesses and (displaced) friction velocities. Also, the peaks of the terrain slope spectra correspond to scales much smaller than the domain extent, so the terrain was statistically well-sampled; the loss at the lowest wavenumbers due to taking flow statistics from the downwind 2/3 of each domain, did not affect the terrain statistics (< 1% loss in standard deviations of slopes). Further, when taken over the lateral extent considered (or just the 4 km width corresponding to the flow variables output for each direction simulated), the terrain statistics in the downwind 2/3 of the domains differ negligibly from the statistics over the entire (unsmoothed portion of the) domain: i.e., despite the progressive smoothing of the outermost parts approaching the domain edges, the terrain statistics were essentially homogeneous. Wieringa (1993) discussed the general notion of a transition sublayer over heterogeneous terrain, whose depth was also known as the blending height (z b ); above this height, horizontal homogeneity of the mean flow is approached. Similarly over forests, a logarithmic profile and effective roughness length can be found for the areal-mean flow above a displacement height d; this was derived for flow over forest canopies (Thom 1971;Finnigan 2002;Sogachev and Kelly 2016) and has become a common approach there. In order for the concept of an effective roughness to be valid, the displaced logarithmic form (2) should apply to the mean wind speed profile; an effective displacement height d eff can be defined as the height above which a log-law describes the temporal and areal mean speed following the surface, depending on the character (statistics) of the terrain. 13 Towards describing the flow, it is useful to consider how the kinematic streamwise momentum flux (negative shear stress) uw or local friction velocity u * ≡ (−uw) 1/2 behaves, particularly the vertical profile of its horizontal-area mean in terrain-following coordinates. Similar to behaviour in the ideal surface layer over a uniform rough surface, the height of maximum u * xy defines d eff ; above this virtual displaced surface, u * xy decreases with z, as in a flat homogeneous ABL. This is shown in Fig. 5, which presents vertical profiles of the terrain-following areal (and implicitly temporal) mean u * xy from the RANS simulations for all directions at the sites considered. Figure 5 illustrates that for different directions in a given area, the character of the terrain affects the height of maximum u * xy (i.e. d eff ), as well as the magnitude of u * xy and shape of its vertical profile. Larger variations in both of these are seen for increasingly complex terrain; this is elucidated further below.

Diagnosing Effective Friction Velocity and Roughness Length
To find the roughness length representative of the terrain's effect on the flow in a given direction over a given area, there are a number of methods one could consider, when invoking the logarithmic profile (2). In addition to considering the flow above d eff , these generally depend on how one finds an effective friction velocity (u * ,eff ).
Starting with the blending-height (z b ) concept, Wieringa (1986) and Mason (1988) proposed finding u * ,eff at z = z b . This was used by Bou-Zeid et al. (2004) to find z 0,eff over a heterogeneous flat surface possessing different z 0 . However, picking u * at just one height can incur significant uncertainty, due to vertical variations in u * xy . A more robust method is to fit the (displaced) log law to the profile of terrain-following areal mean wind speed above the estimated d eff . This is what we do here; it gives both u * ,eff and z 0,eff , with less sensitivity to diagnosed (or estimated) d eff than selecting u * xy from a single height, as shown below. Further, it is sensible given the area-mean friction velocity profiles observed: ∂ u * xy /∂z is For the results and analysis presented here, we calculate vertical profiles of U xy and u * xy in terrain-following coordinates, for each wind direction. In order to examine and mitigate the effects of adjustment of the inflow, we first consider areal means for calculating flow statistics using both the downwind two-thirds of the domain (i.e., excluding the first ∼10 km of the flow) in addition to flow statistics calculated over the entire horizontal domain.
Examining the plane-mean profiles of friction velocity and wind speed, we indeed find a logarithmic profile for U xy above the height of maximum u * xy , following (2); again we take this maximum to be the effective displacement height, d eff . Examples of this logarithmic fit and the corresponding planar-mean profiles of u * are shown in Fig. 6.
Instead of showing dozens of profiles and fits, in Fig. 6 we present the two most extreme and different cases: one site and direction where the terrain is quite complex, such that the downwind-2/3 and full-domain u * ,eff differ the most; and a site where the terrain variations are mildest but with a flatter u * xy profile, which gives largest variation in d eff . The maximum u * xy and the height of its occurence (d eff ) are both higher in the downwind subdomain, due to the full-domain flow still retaining some residual artifacts of the more uniform inflow. One can see that even in these most deviant cases (compared to other flow directions at all sites), the logarithmic profile fits the data well. Also slightly notable is that near the top of the domain there can be a deviation from the terrain-induced mean log law, where the background upwind roughness over smoothed simulated inflow surface is still felt. To mitigate this effect, Fig. 6 Terrain-following areal-averaged profiles of wind speed U (z) and friction velocity u * (z) for two cases, to illustrate diagnosis of effective roughness and displacement height. Blue lines denote full-domain values; red lines denote 'downwind' 2/3 of domain values. Bars (on U ) and points (on u * ) correspond to d, and dashed lines correspond to log-law fits above d for the full and 2/3-domain profiles when fitting (2) above d eff , the top height of the U xy profile being fit was taken to be the maximum of half the domain height and 3d eff ; thus one can see the dotted/dashed lines (fit log-law) in Fig. 6 falling slightly below the mean profiles at the very top of the left-hand plot. Further, to check this fitting, the dimensionless sensitivity of fitted z 0,eff to the upper cutoff height z cut was calculated: d ln z 0,eff /d ln z cut was found to be on average less than 0.5%, and did not exceed 1%.

Prediction of Effective Displacement Height
One might expect simple relations based on σ h -perhaps like previously postulated dependences for z 0,eff (such as Eqs. 5-9)-to suffice for predicting d eff . But σ h depends on the domain size, as shown earlier considering the terrain spectra ( Fig. 1): larger domains include yet more contributions from smaller wavenumbers (larger scales). Further, σ h does not account for the steepness of slopes and contains no information about different steepness occurring in opposite directions. Thus low-order statistics of terrain slope-such as the standard deviation (σ ) of ∂h/∂ x, its finite-difference approximation h/ x, or the upslope finite-difference h + / x ≡ ( h/ x)| h>0 -are expected to be more appropriate for predicting d eff and ultimately z 0,eff . Figure 7 demonstrates this, showing simple predictions for d eff assuming it is proportional to σ h , σ h/ x , or σ h + / x , respectively. In the figure the estimates for d eff are plotted against diagnosed d eff , where the latter is again simply the height of maximum stress (thus maximum u * ) as shown earlier in Fig. 5; the respective proportionality constants relating σ h , σ h/ x , and σ h + / x to d eff are shown above each plot.
The effective diplacement height is somewhat difficult to predict accurately from terrain statistics, though a monotonic and roughly linear relationship is observed in Fig. 7. The figure demonstrates that σ h might only predict a lower-bound for d eff (but it may contain some useful information for simpler sites). However, taking d eff to be simply proportional to σ h/ x or σ h + / x , as: provides modestly good approximation (with an r.m.s. error of ∼40%). While one can see that σ h/ x and σ h + / x are much better than σ h for predicting d eff , they also require prescription of a length scale (r d + or r d ). For the estimations shown in Fig. 7, these were found to be r d + = 1 km and r d = 1.65 km; these values give distributions of d pred /d eff with maxima at 1, as does the coefficient in d eff ≈ 1.6σ h ; the latter approximiation provides a much cruder estimate. The figure also shows that the upslope statistic σ h + / x provides little improvement over σ h/ x . This is likely because calculating the variance of upslopes does not consider the shadowing effect that bigger hills have on the smaller ones immediately downwind. It also hints at the difficulty of defining a unique upslope version of terrain-slope variance; e.g., including shadowing effects would itself introduce an extra length scale, and we leave such to future work. We also note that the slope statistics thus far have been calculated using first-order finite differences, and thus differ from (are smaller than) σ ∂h/∂ x due to the filtering effect at the finest scales. This can be seen considering the previous section on terrain spectra, as in Figs. 1 and 2. We also calculated σ ∂h/∂ x , by integrating the spectra F (dh/dx) = k 2 x F (h) over wavenumber and averaging across all y for each wind direction. But use of the exact slope statistic σ ∂h/∂ x leads to predictions of d eff significantly poorer than those using σ h/ x or σ h + / x . Figure 8 compares predictions using areal σ h/ x and σ ∂h/∂ x , respectively, showing the latter to be less accurate. We postulate that the decreased utility of the exact terrain slope ∂h/∂ x, i.e., the increased error and reduced correlation with d eff , is due to the smallest-scale fluctuations contributing a relatively larger amount of variance (in part due to sampling issues as discussed earlier in Sect. 2.1), while contributing relatively little to the aggregate drag. The smallest features are the most likely to be sheltered by larger hills upwind, and the pressure-blocking effect of larger hills contributes more to the drag. Further, since the maps used for the simulations had 20-m horizontal resolution, the 56 m-resolution output data (which we use to calculate statistics) is not likely to be affected by any processing on the original maps.
Given that expressions simply based on slope statistics such as these also require a length scale (such as r d + or r d in Eq. 12), dimensionally one might expect that a form such as σ h σ h/ x would improve results; in fact, combinations involving σ h were found to degrade estimates of d eff . We also note that fits using combinations of σ h and σ h/ x (including additive and multiplicative forms, including exponents on each σ ) did not give d eff estimates appreciably better than those given in (12). Linear (arithmetic, additive) combinations of h/ x e μ xy and σ e σ h/ x , where e μ and e σ are exponents of order 1 from fits to data, gave only slightly better results than the simple forms in (12).
However, it is not necessary to have a very accurate d eff in order to obtain z 0,eff . That is, fitting planar-mean wind profiles above d eff to the logarithmic form (2) does not significantly affect z 0,eff ; more on this is shown in the following subsections. We also note that using upslope standard deviations, whether σ e σ + h + / x by itself or in combination with h/ x e μ xy , did not improve results; again using a simple form in terms of σ h + / x gave slightly worse results than using σ h/ x , as shown in Fig. 7.

Prediction of Effective Friction Velocity
Associated with the displaced areal-mean logarithmic flow is an effective friction velocity, u * ,eff ; this is diagnosed via log-law fits to the areal-mean wind profile above d eff , as detailed further in the next section. As with the effective displacement d eff , the friction velocity u * ,eff can be predicted using a simple form based on the standard deviation of terrain slope or upslope. The most robust form is found to be a simple linear relationship between u * ,eff /u * in and the standard deviation of slope; for terrain slope and upslope, we find: These give predictions with an r.m.s. error of 8-9%; expression (13), based on σ h/ x , performs slightly better. More complicated empirical relations may be fit, but do not offer much improvement over the above; further, they tend to lack physical justification. However, crudely including the effect of the cross-sectional area through the mean absolute lateral slope μ | h/ y| is possible via: and is seen to give a small amount of improvement for c sμ between 4 and 5; e.g., for c sμ = 4.5 the r.m.s. error is reduced to 6%. 14 This is shown in Fig. 9, along with predictions using (13) and (14). We also note that a pair of relations like (13)-(14), but which add to u * in instead of scaling by it, perform slightly better. That is, the form: and its counterpart u * ,eff = u * in + a u *up σ h + / x reduce the r.m.s. error to 6-7%, if choosing a u * = 1.7 m s −1 and a u *up = 2.8 m s −1 . Extending them by incorporating (1 − c sμ μ | h/ y| )   (15), but with c sμ = 3, reduces the error to 5-6%. Nevertheless, it is not clear what the velocity scales a u * and a u *up are comprised of or represent, and the improvement they appear to offer is small compared to the error itself. We view (13)-(14) as defensible, with more work needed to understand or endorse additive forms like (16) which use a characteristic velocity scale that is independent of u * in , though such independence is logical considering the terrain-affected integral length scales of the turbulent flow.

Performance of Earlier Theories
For context and comparison, we first examine the behavior of the parametrisations which precede this work and have been discussed in the previous sections. Using the simulated flow data and corresponding terrain statistics discussed in Sects. 2-3, we calculated z 0,eff via the form (4) of Flack and Schultz (2010), (6) of Anderson and Meneveau (2011), and (8) of Kelly (2016). For the Anderson and Meneveau (2011) form we have already obtained the high-wavenumber spectral exponent β from fits to the terrain spectra in every direction at every site; then (7) gives α(β) for use in (6). Figure 10 shows z 0,eff predicted by these earlier formulations, versus the effective roughness diagnosed from the simulated flows. The figure shows that none of the extant forms capture the effective roughness, giving predictions which do not generally follow z 0,eff . The Flack and Schultz (2010) parametrisation overpredicts z 0,eff by one order of magnitude, though it does show a slight correlation (upward trend) for z 0,eff beyond 1 m. Similarly, the Kelly (2016) form overpredicts z 0,eff , but not as Fig. 11 Left: z 0,eff predicted using only terrain statistics, via (17)-(18); right: predictions using terrain statistics along with diagnosed d eff , via (19)-(20). Both have γ = γ + = 3. Circles are results using terrain slopes, triangles using terrain-upslopes severely, though it does not increase as much with z 0,eff . The Anderson and Meneveau (2011) form, on the other hand, underpredicts z 0,eff by roughly an order of magnitude, and is not correlated with z 0,eff . In the right-hand plot of Fig. 10 one can also see that the general form (9) can be made to give predictions within one order of magnitude of diagnosed z 0,eff by using a = b = 2 and c = 0.01, but it does not really follow z 0,eff .

New Forms for Effective Roughness Based on Terrain Slope Statistics
Considering all of the data analyzed from the five areas, using only the standard deviation of terrain slope σ h/ x , an optimal fit in logarithmic z 0,eff space was found when using the form: considering r.m.s. upslope σ h + / x , the analogous relationship: emerged. Whether using σ h/ x or σ h + / x , i.e., for both (17) and (18) above, the best fits were found when the exponent value was between 2.5 and 3. For γ = γ + = 2.5, characteristic length scales r s 125 m and r s + 500 m were found; but for γ = γ + = 3, the corresponding scales were found to be r s 325 m and r s + 1450 m. The standard deviation of ln(z 0eff,pred /z 0eff,obs ) was essentially equal for both sets of exponents and coefficients mentioned above, but we note that a better fit across all z 0,eff was exhibited for γ = γ + = 3. The left-hand plot of Fig. 11 shows estimates of the effective roughness length employing (17) and (18), with the optimal exponent γ = γ + = 3. Given the character of the problem, one may expect the length scale in simple formulations such as (17) or (18) to have physical meaning; we note r s is on the order of the diagnosed terrain-induced displacement over all sites and directions, d eff 150 m. Further noting that earlier in (12) we found d eff ∝ σ h/ x , we anticipate z 0,eff to then be proportional to either d eff σ γ −1 h/ x or alternately d eff σ γ + −1 h + / x . Indeed fitting diagnosed z 0,eff and using d eff as the characteristic length scale (instead of r s or r s + ), we find the best results are offered by: and:   (23), using the terrain-driven part of (19), z 0,terrain ∝ d eff σ 2 h/ x (black triangles); blue circles are (19) as in Fig. 11. Right: error metric of predictions, ln(z 0,eff | pred /z 0,eff | obs ), corresponding to triangles in left-hand figure; colours correspond to different terrains as in Fig. 7 where the pressure scale-height is Z p , and the constant c σ is of order 1. By fitting to the diagnosed z 0,eff we find the scale Z p 100 m and c σ 4, but this still gives predictions with 40% larger error than (17); the analagous version of (22) using upslopes gives yet worse predictions. However, taking the scale height in (22) to be proportional to the diagnosed effective displacement height, optimal results using (22) arise for Z p 0.06d eff , with the same r.m.s. log-error as (19); we note that expression (21), which also uses | h/ y| , gives better results. Alternate forms were also examined for the terrain-drag portion of (22), but did not give reasonable predictions (e.g. using other combinations involving | h/ x| and/or means and standard deviations of h/ x, instead of the product σ h/ x | h/ y| ).
Recognizing the physical basis for the summation of roughness-and terrain-induced stresses which gave rise to (3), if the latter is applied with the terrain-induced contributions shown in (17)-(21), we have the generalised expression: where z 0,terrain is equivalent to the right-most terms in each of (17)-(21); i.e., z 0,terrain ≡ (z 0,eff −z 0,in ) for linear combinations of roughness. Here in (23) we have written Z p = a d d eff with a d a constant; this is preferable to attempting to find some Z p applicable over all terrain, which is physically difficult to defend, as the pressure scale-height should depend on the terrain statistics. It is analogous to use of c d d eff in (19) giving substantially better results than the implied all-terrain scale r s in (17). For the form z 0,terrain = 1 3 d eff σ 2 h/ x of (19), fitting (23) to diagnosed z 0,eff we find a d 0.04; the mean of implied pressure scale-heights is then 0.04 d eff 6 m, varying from about 2-17 m over the range of terrains considered here. 15 The resultant z 0,eff are shown in Fig. 13, along with those predicted by (19) alone.
From the left-hand plot in the figure one can see that the two expressions are essentially identical over complex terrain with z 0,eff 1 m; over simpler terrain, the form (23) gives yet better predictions. The right-hand plot in Fig. 13 also shows that using (23) with z 0,terrain from (19) gives predictions of z 0,eff that do not appear to favor complex or simpler terrain. Similar to using (23) with z 0,terrain from (19), if using the z 0,terrain part of (20) in (23), one obtains analogous results, with slightly better predictions by (23) at the smallest z 0,eff compared to (20); however, this is still not an improvement over (19) or using its z 0,terrain in (23). Further, use of z 0,terrain from (21) in (23) gives slightly improved predictions over (21) alone, though negligibly compared to the r.m.s. log error of 0.25 found for (21).
Although the summed-stress form (23) gives better results for smaller z 0,eff when employing z 0,terrain = 1 3 d eff σ 2 h/ x from (19) with a d = 0.04, there is no clear basis for this value of a d . Considering the coefficient C a of (3) is proportional to the ratio of hill cross-sectional area to domain area (Wood and Mason 1993), one might might suppose a d ∝ dh + /dx; but taking a d ∝ h + / x considerably degraded predictions when used in (23).
Earlier we obtained estimates for d eff such as (12), which could also be used with (19) or (20), but doing so does not improve predictions over the simple estimates given by (17)-(18). Further, one could also simply invert the log law, using u * ,eff based on terrain statistics; using (13) and assuming d eff L z , as z → L z one obtains: However, this relation, as well as its analogue using h + / x from (14), result in logprediction errors 2-3 times larger than using (17)-(20); because of this, and since it depends on the domain depth L z , use of (24) is not recommended.

Discussion
Above we found from RANS simulations over large areas (> 40 × 40 km) of complex terrain that the horizontally (areally) averaged flow exhibits a displaced logarithmic wind speed profile, with an effective roughness and displacement height diagnosed from the output of O(100) simulations. Effective roughness formulations based on terrain elevation variance (e.g. Flack and Schultz 2010;Anderson and Meneveau 2011;Kelly 2016) have difficulty predicting z 0,eff within one order of magnitude, over the range of terrain complexities found in nature. This is due primarily to form-drag being connected to the slopes of hills, while σ h is not directly related to slope statistics-unless the terrain elevation (and slope) spectra exhibit a simple power-law dependence. As shown in Sect. 2.1, actual terrain slope spectra exhibit a more complex behavior, with a peak occuring at scales between 1 km and 10 km. For example, the Anderson and Meneveau (2011) expression (6) underpredicts z 0,eff by at least an order of magnitude for z 0,eff greater than about 0.5 m, increasingly so for more complex terrain. This is presumably due to its design originally being for smaller-scale roughness, where a fractal-like surface was (reasonably) assumed to induce the drag. Fractal terrain, i.e. those whose surface-slope spectra have simple power-law behavior, permit direct connection between σ h and σ h/ x ; this kind of idealised terrain allows σ h -based forms for effective roughness. The Flack and Schultz (2010) formulation (4) exhibits the opposite behavior to Anderson and Meneveau (2011), predicting within an order of magnitude for most cases with z 0,eff >∼ 5 m but giving worse results for less complex sites. Incorporation of the skewness of h in Flack and Schultz's model provides a slight correlation in predictions (rising with increasing z 0,eff ); we find Sk h exhibits modest correlation with z 0,eff only for more complex cases having z 0,eff > 1, while it offers no predictive value for less complex terrain.

Combining Surface Roughness with Roughness Induced by Terrain Slopes
The addition of stresses caused by complex terrain at different distances upwind is still somewhat of an open research question, due in part to the difficulty of defining or calculating a momentum-flux footprint. One can expect the extent of upwind terrain that contributes significantly to z 0,eff over a given area to depend on d eff and the terrain-slope statistics themselves, though analysis of such is beyond the scope of the current article. At any rate, for practical use, the excess stress caused by hills (pressure drag) must be combined along with the roughness-induced stress for a given area-or expressed as an effective roughness length which is somehow added to the background z 0 .

Summation of Stresses
A number of works have summed stresses due to heterogeneous surfaces, to obtain an effective roughness. This has included contributions from different roughnesses (i.e. patches) on a flat surface (e.g. Wieringa 1986), or summing contributions due to a collection of identical hills (e.g. Wood and Mason 1993). These provide z 0,eff formulae which all have the same mathematical character as (3) and (22)-(23), i.e. quadratic in 1/ ln(Z p /z 0,i ) where Z p is a pressure-drag scale height and z 0,i represent the different roughness contributions. The form (23) indeed offers good predictions of z 0,eff , when using (17)-(21) for the terrain-induced part of roughness; e.g. z 0,terrain = d eff σ 2 h/ x /3 from (19), as shown earlier in Fig. 13. However, such combination via summation of stresses, as in (23), requires knowledge of the effective pressure scale-height a d d eff ; one might question how 'universal' the value of a d is, as its physical basis is not clear. We note that different values of a d were found, depending on whether we used z 0,terrain (σ h/ x ) from (19), z 0,terrain (σ h + / x ) from (20), or alternately the z 0,terrain (σ h/ x , μ | h/ y| ) from (21) for the terrain-drag contribution to (23). Mason (1988) reported this vertical scale to be ∼ L c /200 for summing drag quadratically as in (23), where L c is the characteristic horizontal scale of variations; here taking L c to be the peak scale from terrain-slope spectra, this and a d = 0.04 would imply L c ∼ 8d eff , though we find d eff to be roughly 1/25th to 1/15th of the peak spectral scale; this difference is reasonable, considering the PDF and spectra of slopes are quite broad (with significant contributions spanning several orders of magnitude of slope and wavenumber, respectively). Although a d is of similar magnitude as | h/ y| , and the constant C a in Wood and Mason's formulation suggests use of such (it is proportional to sillhouette hλ y per area for uniform hills), incorporating first-or second-order statistics of h/ y within the pressure-scale height (via a d ) caused degradation of predictions. 16

Simpler Combinations
The best results for z 0,eff were obtained using the stress-based form (23) to add contributions from background roughness and terrain drag, as shown in the previous section. However, we also find robust results using simple addition: the inflow roughness z 0,in , which corresponds to the profile over flat terrain upwind and the surface roughness, is simply added to the terrain-induced drag contribution as in (17)-(20). Such summation gives the same results as stress-based addition, except for some lower-complexity cases where it underestimates z 0,eff , as shown in Fig. 11. We again note that (19)-(20) offered better predictions compared to (17)-(18); i.e., d eff appears to be the characteristic length scale which acts as a coefficient multiplying the terrain slope variance and carries additional information about z 0,eff , while use of a constant length scale resulted in ln z 0,eff error that was at least 1/3 larger.
Alternatively, formulations like (5) of Wan and Porte-Agel (2011) and (6) of Anderson and Meneveau (2011) follow a quadratic form, summing the squares of roughness contributions; adapting these for use with terrain-slope variance, they can be generally expressed with the form z 2 0,eff = z 2 0,in + z 2 0,sl , where z 0,sl is the terrain slope-induced component. Employing this form to combine z 0,in with the respective z 0,sl parts of (17)-(21), we get results nearly equal to the linear combinations (17)-(21) for more complex terrain, z 0,eff 1 m. Unsurprisingly, however, over less complex terrain this quadratic combination causes yet more substantial underpredictions; this follows by noting that:

Limitations and Interpretation
For clarity in assessing the slope-induced contributions to z 0,eff , we used only a uniform background (surface) roughness in all simulations, z 0,in = 9 cm. Previous studies (see Sect. 1) indicate that other (uniform) values of z 0,in would combine with terrain slope-induced roughnesses in the same way as found here; however, for smaller z 0,in , our results and analysis imply that simple addition of z 0,in to z 0,sl would be equal to the stress-based combination over a wider range of z 0,eff than shown here, i.e. down to lower z 0,eff . For areas where z 0,in is appreciably larger than the 9 cm considered here, the stress-based summation of z 0 would give relatively better results, but such large roughnesses generally occur only in forested or built-up areas, which tend to be inhomogeneous. We do not consider inhomogeneous surface roughnesses, as this is beyond the work of the present study. The high-resolution part of the domains in the simulations had limited extent, with the amount of high-wavenumber terrain-slope variance decreasing with distance beyond 7 km from the center (due to progressive smoothing, see Sect. 3); this could be interpreted as a weak violation of homogeneity of terrain statistics. However, the limited extent of the high-resolution area had negligible effect, for multiple reasons: first, the relative contribution of smallest-scale terrain-slope variance (2π/k <∼ 300 m) was much smaller than that from larger scales; secondly, the use of the downwind two-thirds of the domain for analysis reduced this effect; third, the output data analyzed had a resolution of 56 m, corresponding to a region about 8-10 km from domain center. Reduction of the analysis domain to the downwind 2/3 portion for each wind direction was intentionally done, to reduce the aforementioned issue: conducting tests where we analysed subdomains progressively further downwind and with decreasing areas, we found that the displacement heights first increased and then converged, with u * ,eff and z 0,eff following suit.
The results here correspond to terrain-induced momentum transfer under neutral atmospheric stratification, such that hundreds of metres above the terrain a logarithmic profile is recovered-with an effective roughness that we can predict based on the terrain statistics. In that regard, the simulations and analysis have been capturing 'bottom-up' effects, ignoring anything from above; we remind that the RANS simulations have no capping temperature inversion, but rather a tall domain much larger than the diagnosed terrain-displacement heights. Consequently for shallower ABLs observed to have depths approaching or less than the terrain-induced mean displacement (z i 3d) occurring under neutral conditions, obvi-ously the theory presented here would break down. Stability effects are also neglected, but the effective roughness results and expressions here can be taken as a representative average: most places have nearly neutral conditions in the mean, i.e. the long-term distribution of inverse Obukhov length has a significant peak and is centered around zero (Kelly and Gryning 2010).

Beyond 1h/1x : Consideration of Other Terrain Statistics
We have shown that the variance of streamwise terrain slope is the key quantity for characterizing the effective roughness of areal-mean wind profiles displaced by complex terrain.
In particular, σ h/ x was found to be the optimal predictive terrain statistic, in contrast with σ ∂h/∂ x . Statistics of ∂h/∂ x (calculated in Fourier space) were not as well-correlated with z 0,eff and offered less predictive power; we attribute this to high-wavenumber noise, as well as the more local impact of the smaller scales filtered out by h/ x. For example, analogous relations to (17)-(21), but using σ ∂h/∂ x instead of σ h/ x (with different fitted constants), gave 2-4 times larger error in predicted ln(z 0,eff ).
Streamwise upslope statistics were also considered, since the flow experiences different slopes, and presumably different terrain drag, in opposite directions; we thus expect statistics of h + , notably σ h + / x , to facilitate better prediction of effective roughness. Indeed u * ,eff , d eff , and z 0,eff were all found to differ in opposing directions. The effective roughness differs by a factor as large as 5, while the rectified geometric mean 17 of z 0,eff|ϕ+180 • /z 0,eff|ϕ , over all directions ϕ and sites considered, is a roughly 1.6; i.e., on average z 0,eff in a given direction is different from z 0,eff in the opposite direction by a factor of about 1.6 ±1 . Complementing this, the difference in d eff seen from opposite directions was on average about 40% of d eff from either direction, with the largest difference having d eff being three times larger in one direction than its opposite. However, σ h + / x did not give better predictions of z 0,eff (or d eff ) compared to σ h/ x . The inability of σ h + / x to predict this asymmetry appears to be due to the simple definition of h + , along with the complexity of the flow: when considering upslopes, h + does not account for the sheltering effect of larger hills on smaller hills. In other words, the form drag due to a hill can also depend on the flow both upwind and downwind; a hill in the wake of a significantly taller (but not necessarily steeper) mountain will tend to provide relatively less drag than if it were not sheltered. Introducing a spatially-varying weighting, to incorporate sheltering in the definition of h + , could in principle help to exploit observed slope asymmetry; but doing so may depend on the terrain shape itself, quickly becomes complicated, and is beyond the scope of current work.
Statistics of the lateral slope were also considered here, because form drag for single hills has been shown to relate to frontal area (e.g. Wood and Mason 1993) and thus h/ y. The use of | h/ y| xy as in (21) indeed added small improvements to predictions of z 0,eff , shown earlier in Fig. 12; in contrast, σ h/ y is not seen to carry information about z 0,eff . However, the improvements which can be offered by using a representative frontal area or h/ y are also limited by the sheltering aspect.

Turbulence Model Considerations
As noted by Wood and Mason (1993) and later summarised by Wood et al. (2001), the amount of form drag captured will tend to depend on the turbulence closure, with eddy-viscosity closures giving larger drag and z 0,eff compared to second-order closure. But, as Ayotte et al. (1994) pointed out, the mean flow over hills will only be weakly affected by the turbulence parametrisation. A weakness of commonly-used two-equation RANS closures, such as the k-one used here, is the lack of accounting for the pressure-strain terms (see e.g. Pope 2000) which become significant due to hills. The results of the current work are expected to hold for RANS using typical (eddy-viscosity) closures. Use of more advanced turbulence modelling which includes e.g. pressure-transport or pressure-related anisotropy for RANS (e.g. Reece and Rodi 1975;Speziale et al. 1991;Wallin and Johansson 2000), or alternately using LES, will reduce the total form drag; the latter would subsequently require reduction of the leading coefficient in each of (17)-(21), consistent with the reduction in Wood and Mason (1993) of their shear factor. The equations given in the current work may thus slightly overpredict z 0,eff if compared with LES or RANS with higher-order closure. However, for common stand-alone application or for RANS driven by mesoscale modelling, this will not be an issue: the forms given here correspond to the displaced effective roughness reproduced by RANS with eddy-viscosity closure.

Application and Recommendations
Because the terrain-slope variance is easier to calculate than the corresponding upslope statistic, and because the latter was not found to improve predictions, we recommend use of (17), (19), or (21) to find z 0,eff in applications. Though it is not as accurate as (19) or (21), the form (17), using γ = γ + = 3, serves as a practical estimate of effective roughness in applications where one cannot observe or diagnose d eff a priori. But given the results shown in Sect. 4.3.2, for applications where d eff can be diagnosed following the method given in Sect. 3.3 (such as those involving RANS), we recommend the simple and robust form (19). In such applications, if the map and model resolution is finer than ∼1 km and the lateral slopes can be readily calculated, one may obtain further improvement by using (21).

Summary and Outlook
For brevity, we summarise what we have found and reported, via bullet points below.
• Complex terrain causes a displaced areal-mean logarithmic profile, with effective displacement height d eff , friction velocity u * ,eff , and roughness length z 0,eff .
• Previous works cast z 0,eff in terms of σ h , but σ h does not correlate with z 0,eff unless terrain is mono-fractal.
-Previous methods mis-predict z 0,eff by one order of magnitude.
σ h/ x using first-order finite differences gives robust results. · Incorporating d eff from flow field improves results, where d eff is diagnosed from terrain-following plane-mean profiles u * (z). · Further improvement found by also using lateral terrain slope. -Using up-slope σ h + / x does not improve over σ h/ x ; · this upslope statistic does not account for hill sheltering.
We note that the arguments of Beljaars et al. (2004) for choosing to parametrise terraininduced stress (instead of z 0,eff ) can still be considered, if one is looking to parametrise unresolved terrain drag in mesoscale and climate models having (effective) resolutions of several kilometres or greater. Consistent with this, and because we have seen that terrain slope spectra can differ significantly from the Beljaars et al. (2004) form (see Fig. 1), ongoing work includes extension of (17)-(21) to further account for limited mesoscale resolutions. One can argue that simple analytical adjustment of our z 0,eff for a given model resolution is possible, if one knows σ h/ x along with the high-wavenumber spectral exponent β and peak of the terrain-slope spectrum for a given terrain (map); however, this needs to be further explored and verified.
Our findings and statistical relations for terrain-induced effective roughness have multiple applications involving microscale flow models, not limited to RANS. This includes choice of z 0,eff for inflow profiles to microscale flow models, to reduce the amount and spatial extent of adjustment of the incoming flow as it approaches finely-resolved terrain. Another usage involves the geostrophic drag law (GDL) in applications such as wind energy; use of z 0,eff in the GDL allows observations from an area of one terrain complexity to be used for predictions in another area following the European Wind Atlas method (Troen and Petersen 1989), whether using e.g. linearised flow models or RANS solvers. The forms given here also facilitate better coupling of mesoscale with RANS models: output from the former, which typically do not resolve the drag below scales of several or even tens of kilometres (e.g. for WRF having an effective resolution of 6-8 times its horizontal grid spacing, Pielke 2013), may not resolve the peak part of the terrain-slope spectrum, while typical RANS application resolves the terrain drag as found in this work. Mesoscale model wind speed (and possibly velocity) output could be scaled using the ratio of low-pass filtered terrain slope variance to said variance resolved by a microscale model, per direction; study and validation of such is an ongoing task. Brown and Wood (2003) found for stable conditions that using an effective roughness length, independent of stability, gave reasonable results within one-dimensional models predicting total surface drag. Future work can include use of RANS solvers in stable conditions (see van der Laan et al. 2017van der Laan et al. , 2020 to check the validity of such, and potentially extend our z 0,eff forms for these conditions. Further, RANS with length-scale limited k-turbulence closure (Apsley and Castro 1997) or capping inversion can be used to investigate the development of displaced mean profiles where the ABL depth is limited and stable stratification spreads downward (Kelly et al. 2019b), though we note that RANS solvers generally cannot address hill-or inversion-induced waves. Large-eddy simulation can more effectively facilitate such investigation, though it is currently expensive in terms of computational resources.
Other ongoing and future work involves continued investigation into discoveries made in this study, but beyond the scope of the current article. This includes: a relation observed between the profile of terrain-following areal mean eddy viscosity and shape of the probability density function of terrain slopes, which can be used to inform or improve e.g. single-column models or PBL schemes; relations between terrain-slope PDF shape and the high-wavenumber power-law exponent β; and relating spatial spectra of mean (Reynolds-averaged) velocity and shear stress to terrain spectra at heights approaching (or exceeding) d eff . Growth of d eff over complex terrain, downwind of simpler terrain, is also ripe for investigation; this goes along with further testing of the results here in larger finely-resolved domains, as well as investigation over terrain having different aerodynamic roughnesses.