Large Eddy Simulation Study of Atmospheric Boundary Layer Flow over an Abrupt Rough-to-Smooth Surface Roughness Transition

The atmospheric boundary layer flow downstream of an abrupt rough-to-smooth surface roughness transition is studied using large eddy simulations (LES) for a range of surface roughness ratios. Standard wall models assume horizontal homogeneity and are inapplicable for heterogeneous surfaces. Two heterogeneous-surface wall models are evaluated, one based on a local application of similarity theory using a twice-filtered velocity field (BZ model) and another based on a local friction-velocity obtained by blending the upstream and downstream profiles (APA model). The wall shear stress and the turbulence intensity (TI) are sensitive to the wall model while the mean streamwise velocity and the total shear stress (TSS) are less sensitive. The APA model is more accurate than the BZ model on comparison to previous experiments. The wall shear stress obtained using the APA wall model is sensitive to the ratio of the equilibrium and the internal boundary layer (IBL) heights, while other statistics are not. The IBL height is insensitive to the turbulent quantity (TSS or TI) on which it is based. Several analytical relations for the IBL height are evaluated using the LES data. Two models are found to be accurate for different roughness ratios while one model is reasonable over the full range investigated. A phenomenological model is developed for the TI downstream of the roughness jump using a weighted average of the upstream and far-downstream profiles. The model yields reasonable predictions for all roughness ratios investigated.


Introduction
The flow in the atmospheric boundary layer (ABL), formed in the lowest part of the atmosphere, critically affects several aspects of human activity ranging from weather, air quality, agriculture to energy extraction from the wind (Stull 1988). The Earth's surface is ubiquitously heterogeneous (Bou-Zeid et al. 2020) and the fluxes of momentum, heat, moisture and other passive scalars imposed on the flow depend on the type of land surface. Changes in land characteristics can be due to changing landscapes (e.g. a transition from a water body to land) or due to changing land use (grassland, forest land and cultivated or fallow land). High-fidelity large-eddy simulations (LES) of the ABL flow over heterogeneously rough surfaces forms the topic of the current study.
We assume that the surface irregularities (typical height k) are much smaller than the boundary-layer height (δ), and the resulting roughness can be considered as 'sub-surface', with its effect on the flow encoded as an aerodynamic roughness height, z 0 . We restrict attention to the geophysical regime (δ/k 80 as per Jimenez 2004) and a configuration wherein the aerodynamic roughness of the bottom surface undergoes an abrupt change from a relatively higher value (z 01 ) to a lower value (z 02 ) in the direction normal to the mean flow. Upstream of the rough-to-smooth (RS) transition, the flow is in 'equilibrium' with the surface of roughness z 01 . The mean and fluctuations of the streamwise velocity and the vertical momentum flux are given by:ū with z denoting the vertical coordinate, u * 1 being the upstream friction velocity, κ ≈ 0.41 being the von-Karman constant and A, B are empirical constants (Stull 1988;Bou-Zeid et al. 2004;Stevens et al. 2018). Downstream of the RS transition, the flow is a function of the streamwise coordinate (x) as well as z. An internal boundary layer (IBL) is formed close to the bottom surface that grows in height along x as the flow increasingly farther away from the ground is affected by the roughness jump. Different features of the flow downstream of a RS transition have been studied by field observations (Bradley 1968), wind tunnel experiments (Antonia and Luxton 1971;Chamorro and Porté-Agel 2009;Li et al. 2021), analytical models (Elliott 1958;Panofsky and Townsend 1964;Plate and Hidy 1967;Taylor 1969;Chamorro and Porté-Agel 2009;Ghaisas 2020;Li et al. 2022) and numerical simulations (Shir 1972;Rao et al. 1974;Abkar and Porté-Agel 2012;Sridhar et al. 2017;Anderson 2020).
The primary challenge in LES of the flow over heterogeneous surfaces is related to the modelling of the shear stresses at the bottom wall. Since the nominal Reynolds numbers based on the friction velocity and the boundary layer height are very high for atmospheric flows, a common practice to enable simulations on reasonably-sized grids is to neglect the viscous terms from the Navier-Stokes equations and to introduce an additional stress at the bottom wall. The Monin-Obukhov Similarity theory (MOST) (Monin and Obukhov 1959) has provided the most commonly-used wall model formulation for LES of ABL flows (Moeng 1984;Khanna and Brasseur 1997;Brasseur and Wei 2010;Xie et al. 2015). However, since these wall models assume the flow conditions to be statistically identical at all horizontal locations, they are inappropriate for simulating flows over heterogeneous surfaces.
Two wall models accounting for heterogeneously rough ground surface that prescribe the wall shear stress in a localized manner have been proposed in the literature. The wall model by Bou-Zeid et al. (2004), denoted as the 'BZ model' hereafter, is based on filtering the velocity field at a scale larger than the LES-filter width. The wall model by Abkar and Porté-Agel (2012), denoted as the 'APA model', was proposed by recasting a slightly modified diagnostic analytical model of Chamorro and Porté-Agel (2009). This wall model introduces a so-called blending function that allows for the mean streamwise velocity profile to vary smoothly from its upstream profile to its profile far downstream of the roughness jump. The APA model requires specification of the ratio of the equilibrium boundary layer (EBL) height to the IBL height, α = δ e /δ i , as an input parameter. Here, δ e and δ i are the EBL and IBL heights, respectively. The APA model has been tested in LES (Abkar and Porté-Agel 2012) with only one value of α = 0.027 and for only one value of the roughness ratio, m = z 01 /z 02 = 83.3. Furthermore, the results of Abkar and Porté-Agel (2012) focused only on the mean streamwise velocity profiles and the surface shear stress. Other quantities of interest such as the turbulence intensity, the vertical momentum flux and the internal boundary layer height evolution have not been studied using different heterogeneous-surface wall models in a systematic manner over a range of m and α values.
Several of the previously-mentioned analytical models require specification of the IBL height, δ i (x), as an input. A number of analytical and/or empirical models for δ i (x) have in turn been proposed in the literature. The Elliott (1958) model is purely empirical but has been widely used as a building block in several further studies. This model as well as the models of Wood (1982) and Jegede and Foken (1999) assume that the IBL height grows as a power-law, δ i ∼ x 0.8 . The Panofsky and Dutton (1984) model proposes an implicit non-linear relation for δ i (x) that relies only on z 02 . A similar implicit relation, but one that involves the ratio m, was proposed by Savelyev and Taylor (2005). It is unclear as to which of these IBL height models is accurate, particularly over a large range of m.
The above-mentioned theoretical studies mainly focused on modelling the mean streamwise velocity behind an abrupt surface roughness transition. The streamwise turbulence intensity plays a key role in determining fatigue loads on passive structures, such as trees or buildings, and engineering systems, such as solar or wind farms. Despite its importance, analytical models for the streamwise turbulence intensity are largely missing.
This paper describes LES of the flow over an abrupt rough-to-smooth surface roughness transition using a high-order numerical framework. Our LES framework is validated by comparing to two different experimental data sets, of Chamorro andPorté-Agel (2009) andLi et al. (2021). Our aim is to assess the performance of the BZ and APA wall models by evaluating a range of turbulent statistics beyond only the mean velocity and surface shear stress. A second aim is to study the sensitivity of the APA model results to the input parameter α = δ e /δ i . We also report simulation results for different roughness ratios, m = z 01 /z 02 . Several previously developed models for the IBL height are evaluated using our LES results. Finally, a phenomenological model is proposed for the turbulence intensity profile downstream of a step change in surface roughness.
The numerical methodology and cases studied are described in Sect. 2. Details of the wall models that are assessed here are given in Sect. 2.2. Results are presented and discussed in Sect. 3 and conclusions are presented in Sect. 4.

LES Methodology
The incompressible, LES-filtered Navier-Stokes (NS) equations that are solved are: whereũ i is the instantaneous resolved velocity in i-direction, t is the time, and x i with i = 1, 2 and 3 are the three Cartesian coordinate directions which can be used interchangeably with x, y and z. The filtered pressure field is denoted byp and ρ is the (constant) density of air.
The forcing term f i drives the flow and is described below. The tensor τ i j is usually comprised of the viscous stresses (−2νS i j , where ν is the viscosity andS i j is the strain-rate tensor) and the subgrid scale (SGS) stresses (τ sgs i j = u i u j −ũ iũ j ). However, the Reynolds number based on the free-stream velocity, the height of the boundary layer and the viscosity of air is in excess of the order 10 5 in atmospheric flows. As a result, the direct effects of viscosity are negligible over most of the domain except for a very thin region close to the bottom surface. The effects of these extremely thin viscous sub-layer and transition layer are modelled through a wall model and the viscous terms are neglected over the entire domain. The tensor τ i j is thus comprised of the SGS stresses and the wall stresses, τ i j = τ sgs i j + τ wm i j . The anisotropic minimum dissipation (AMD) model (Rozema et al. 2015) is used for subgrid scale closure. The trace of the SGS stress tensor is incorporated along with the pressure while the deviatoric part is given by τ with the eddy viscosity ν sgs given by the AMD model. This model has been extensively tested previously for a variety of flow configurations including simulations of atmospheric boundary layers and turbulent channels (Rozema et al. 2015;Abkar et al. 2016;Vreugdenhil and Taylor 2018;Zahiri and Roohi 2019;Ghaisas et al. 2020). The wall stresses, τ wm i j , are prescribed as discussed in detail in the next section.
The above-mentioned equations are solved using the concurrent precursor simulation method (Stevens et al. 2014). This methodology is similar to the one explained in detail in Ghaisas et al. (2020). Two computational domains of sizes (L x × L y × L z ) each are used. The first, 'precursor', domain is driven by an imposed constant pressure gradient, f i = −u 2 * 1 /L z , and has a homogeneous surface roughness, z 01 . A shifted periodic boundary condition (Munters et al. 2016) is applied to the precursor domain to ensure that spurious, infinitely long, streamwise-aligned streaks do not develop in the domain and contaminate the solution. The second, 'main', domain (see Fig. 1) has an upstream portion with aerodynamic roughness z 01 , followed by a transition to a surface with roughness z 02 . The last portion of the main domain is a 'fringe' region, wherein the flow is nudged towards the same upstream conditions as in the precursor domain using the additional forcing term f i . The surface roughness in the fringe region is the same as in the upstream and precursor domains, i.e. z 01 . The top boundary (z = L z ) is imposed with no-penetration (ũ 3 = 0), free-slip (∂ũ 1 /∂z = ∂ũ 2 /∂z = 0) conditions and all the horizontal (x, y) boundaries are periodic. The bottom boundary is a no-penetration (ũ 3 = 0), slip (ũ i = 0 for i = 1, 2) wall. The slip velocity for the horizontal components is not explicitly specified since these components are stored at a The fringe region has the same roughness as the upwind of the step jump; RP stands for rough patch, SP stands for smooth patch vertically staggered location, z = z/2. The slip velocity is governed primarily by a shear stress is applied using a wall model as described below.
The 'PadeOps-igrid' code (Subramaniam et al. 2021) developed over the years is used for the simulations. This code uses Fourier-spectral discretization in the horizontal (x, y) directions, 6th-order staggered compact finite-difference scheme in the vertical (z) direction and a 3rd-order Runge-Kutta method for time advancement. The code is well validated and has been used previously for several problems including rough-wall turbulent channels (Ghate and Lele 2020), stratified and unstratified atmospheric boundary layers (Ghate and Lele 2017), and problems involving wind turbines or farms Ghaisas et al. 2020;Howland et al. 2020).

Wall Models
The effect of viscosity is modelled by introducing the term τ wm i j in the total stress in the Navier-Stokes equations. Since viscous effects are important only close to the wall, τ wm i j is zero at all points in the domain except at the bottom wall, z = 0. Furthermore, only the vertical shear components (i.e. τ wm i3 with i = 1, 2) are non-zero. The standard 'equilibrium' wall model based on the Monin-Obhukhov Similarity Theory (Moeng 1984) first estimates the magnitude of the mean shear stress (− u * 2 ) by inverting the logarithmic law-of-the-wall. The horizontally-averaged streamwise velocity ( ũ 1 ) available during the LES at the first grid point from the wall (i.e. z = z/2) is used in these models to determine the mean shear stress. The mean shear stress is then distributed into its two components, i.e. τ wm 13 and τ wm 23 , with each component being proportional to the corresponding component of the instantaneous local filtered horizontal velocity, Here, .. denotes a horizontal averaging operation. This model is inapplicable for flows over heterogeneously rough surfaces since it is inappropriate to carry out a horizontal average in such flows. A few studies tried to overcome this problem by applying the logarithmic law-ofthe-wall in a strictly local sense (Albertson and Parlange 1999). However, it is easily shown, using the Cauchy-Schwartz inequality, that this leads to an over-prescription of the wall shear stress. Two previously-proposed wall models that try to account for non-equilibrium effects induced by the presence of the surface roughness step are evaluated in this paper. These models are described next. The first model evaluated here was proposed by Bou-Zeid et al. (2004) and is denoted as the 'BZ' wall model in this paper. This strictly local wall model uses a velocity field that is filtered using a width of size 2 , where = ( x y z) 1/3 is the characteristic grid size. Referring to this velocity field filtered at the 2 -scale as ũ i , the wall shear stress is given by: where the local friction velocity is given by assuming that the 2 −filtered horizontal velocity satisfies the law-of-the-wall locally, In this study, z 0 (x, y) is either z 01 or z 02 depending on the roughness of the underlying surface at the horizontal location given by coordinates (x, y). This model reduces the overprescription of the mean wall shear stress that is observed on applying the law-of-the-wall locally because the 2 -filtered velocity field has much smaller fluctuations (Bou-Zeid et al. 2004). This model also does not require a horizontal averaging operation and hence can be applied to heterogeneously rough surfaces. The second wall model evaluated here is that originally proposed by Chamorro and Porté-Agel (2009) and modified by Abkar and Porté-Agel (2012). We denote this as the 'APA' model. This model uses the same formulation as Eq. 7 but it does not rely on the assumption that the local friction velocity satisfies a local logarithmic law-of-the-wall immediately downstream of an abrupt surface roughness transition. Instead, this model explicitly accounts for the fact that the local friction velocity changes along the streamwise direction and gradually approaches its equilibrium value by using a blending function λ(x, z). The local friction velocity is given by: Here, (..) denotes averaging in the spanwise (y) direction. u * 1 is the friction velocity in the upstream region, which is also the friction velocity in the precursor domain since the flow in the upstream region is driven by the flow in the precursor domain. The blending function is modelled as λ(x, z) = ln[z/αδ i (x)]/ ln [1/α], and is evaluated at x = z/2 in Eq. 9. The IBL height is specified using the empirical relation proposed by Elliott (1958), This model is applicable only in the region where the surface roughness has abruptly changed from its upstream value to z 02 because λ < 1 in this region. The value of λ is set to zero once the equilibrium boundary layer crosses the first computational grid point, i.e. once δ e (x) > z/2. Beyond this streamwise location, the model reduces to the BZ model with the surface roughness equal to z 02 . In the entirety of the precursor domain, in the portion of the main domain that is upstream of the abrupt surface roughness transition, and in the fringe portion of the main domain where the flow is nudged towards the same flow field as in the precursor domain, the surface roughness is uniformly z 01 . In these three regions, the BZ wall model with surface roughness equal to z 01 is applied. In the portion of the main domain between the abrupt transition and the fringe, either the BZ model (Eqs. 7 and 8) with roughness z 02 or the APA model (Eqs. 7 and 9) is applied. We refer to these two combinations as 'BZ' model and 'APA' model respectively.

Cases Simulated
A total of twenty-three LES simulations are carried out which are listed in Table 1. The cases are selected so as to cover different reference experiments, roughness ratios (m = z 01 /z 02 ), wall models, grid sizes and different values of the parameter α = δ e /δ i which is an input to the APA model. Figure  First, a set of six simulations (Cases 1-6) for m = 83.3 covers the two wall models (BZ and APA) and three grid sizes comprised of 128 × 32 × 32 (G1), 192 × 64 × 64 (G2) and 240 × 80 × 80 (G3) points per domain. The actual number of computational points in each simulation is twice that mentioned above and in Table 1. Following a grid sensitivity study, we use the finest grid (G3) for all further runs. The second set of simulations (Cases 7-10) covers two additional values of m = 20 and 125 and the two wall models. In all the cases mentioned above that involve the APA wall model, the value of α = 0.027 is used. To study the sensitivity to this input parameter, six additional cases (11-16) using the APA model, with α = 0.054 and 0.108, for three surface roughness ratios (m = 20, 83.3, 125), are carried out. Further, Cases 17-20 are simulations corresponding to the 'Re21ks16' experiment of Li  (2021) using the BZ and APA wall models with different α values. This experimental data set is selected because it has sufficiently large δ/k for it to be of geophysical relevance. Finally, simulations over homogeneously rough surfaces with the tabulated roughness values are conducted. These cases are useful for developing an analytical model for the turbulence intensity described later in this paper. The upstream friction velocity is u * 1 = 0.5473 m/s following Chamorro and Porté-Agel (2009) for Cases 1-16 and 21-23 whereas for Cases 17-20, u * 1 = 1.0191 m/s as measured by Li et al. (2021). All the simulations are carried over 100 non-dimensional time units (normalized using δ ex p m and u * 1 as reference scales). Statistical averaging is performed over the last 60 time units. Averaging is performed in time and along the horizontal (x, y) directions for the precursor domain and over time and the spanwise (y) direction in the main domain.

Validation with Experiments of Chamorro and Porté-Agel (2009) and Effect of Wall Models
The sensitivity of different statistics to the BZ and APA wall models is studied first. Experimental results of Chamorro and Porté-Agel (2009) (2012), is used throughout this section. Following a grid convergence study described in Appendix 1, 240 × 80 × 80 grid points are used per simulation domain for these runs. Figure 3 shows the streamwise evolution of the wall shear stress (τ ) normalized by its upstream value (τ 0 ). Following the surface roughness jump at x = 0, the shear stress applied by the new smooth surface on the flow is smaller compared to that applied by the upstream rough surface, which results in τ/τ 0 < 1. The BZ model under-predicts the experimental results while the APA wall model results agree well with the experimental results. A careful inspection of Fig. 3 shows that the computed shear stress values using the APA model reduce slightly between x/δ = 0.5 and x/δ = 1.5. As discussed later, in Sect. 3.2, this feature is absent for larger values of α. The asymptotic value (for large x) is the same for both models and is in good agreement with the experimental results. Figure 4 compares the mean velocity profiles obtained from the LES using the two wall models at several locations downstream of the roughness jump. It is clear that the LES results adhere to the law-of-the-wall closely above a certain height. This height above which the downstream and upstream profiles are identical is the internal boundary layer height and is discussed in detail later. Below the IBL height, the APA model results are in slightly better agreement with the experimental results than the BZ model results. This is consistent with the under-predictions of the surface shear stresses found in Fig. 3.
Vertical profiles of T SS obtained from different regions of the simulation domains are shown in Fig.5. The T SS has contributions due to the resolved scales, subgrid scales and the wall model i.e. T SS = u w + τ sgs 13 + τ wm 13 . The profile obtained from averaging over the precursor domain (surface roughness z 01 everywhere) agrees very well with the theoretically expected line (Eq. 3). The profile obtained by averaging over the upstream portion of the The profiles of T SS averaged over the spanwise coordinate at different locations downstream of the surface roughness jump (x > 0) are found to be insensitive to the wall model. Close to the bottom wall, the magnitude of the T SS is smaller compared to its upstream value, consistent with the axial evolution of the surface shear stress shown in Fig. 3 and with the fact that the surface undergoes a RS transition at x = 0, i.e. z 02 < z 01 . Furthermore, it is seen that the T SS varies linearly from its value at the bottom wall to the top of the IBL (marked by dashed gray lines). This indicates that several early analytical models (Panofsky and Townsend 1964;Plate and Hidy 1967;Taylor 1969) were based on an incorrect assumption of constant shear stress within the IBL, but supports the assumption made in recent models (Chamorro and Porté-Agel 2009;Ghaisas 2020).
The turbulence intensity (TI) is defined as the ratio of the standard deviation of the velocity to the mean velocity, Unlike the T SS, the T I downstream of the step is very much sensitive to the wall models as shown in Fig.6. The profile upstream of the roughness jump has larger values of T I close to the wall, once again consistent with the fact that the configuration being studied is a RS transition. The T I value at the wall obtained using the APA model is larger than that obtained using the BZ model, consistent with the surface shear stresses obtained using these two wall models (see Fig. 3). In each panel of Fig. 6, the influence of the changed roughness on the downstream profiles is seen to be prominent near the wall but it disappears after a certain height similar to the mean velocity profiles. This again indicates the presence of an IBL within which the effects of the changed surface roughness are confined.
To study the evolution of the internal boundary layer, the IBL heights extracted from the LES data based on two turbulent statistics, T SS and T I , are presented in Fig. 7. The dashed vertical line marks the jump in aerodynamic roughness from rough to smooth. For each profile, the IBL height is determined as the smallest distance from the bottom wall where the upstream and downstream profiles differ by less than 10%. The IBL profiles for both T SS and T I are insensitive to the wall models as seen from Fig. 7a, b. Also, it is evident from Fig. 7c that the IBL based on T SS profiles are similar to those based on T I profiles. Finally, the empirical relation proposed by Elliott (1958), Eq. (10), for the IBL height is seen to be quite accurate for m = 83.3.

Validation with Experiments of Li et al. (2021) and Sensitivity of APA wall model toĄ
s discussed in Sect. 2.2, the ratio of the equilibrium boundary layer height (δ e ) to the internal boundary layer height (δ i ), i.e. α = δ e /δ i , is an input parameter to the APA model.  To match these experimental conditions, we use z 01 = 6.726 × 10 −5 m and z 02 = 3.203 × 10 −6 m (m = 21). This section also serves to further validate our numerical framework since the LES results are being compared to an additional experimental data set. Figure 8 shows that the wall shear stresses obtained from the APA wall model are clearly sensitive to the parameter α. This is contrary to what was suggested, but not explicitly shown, in the study of Abkar and Porté-Agel (2012). For the simulations corresponding to the Li et al. (2021) data, LES results with α = 0.108 agree best with the experiments, while for the Chamorro and Porté-Agel (2009) cases, the agreement is best for α = 0.027. This is seen more quantitatively in Table 2, where the L 2 norms of the relative errors, expressed as a percentage, are shown. The tabulated values are calculated as: where N is the number of measurement points. The error norms are smallest (around 1-2% each) for α = 0.108 and 0.027 respectively for the two data sets. These results suggest that, for accurately predicting the wall shear stress evolution, the value α = 0.027 recommended by Abkar and Porté-Agel (2012) is not appropriate for all conditions but needs to be tuned as a function of z 01 and m. Figure 8 also shows that the surface shear stress values attained far downstream of the transition are independent of the value used for α. To understand the behaviour of the APA wall model in the intermediate region further, profiles of the blending function, λ(x, z/2), and the spanwise-averaged velocity at the first off-wall grid point,ũ 1 (x, z/2), are shown in Fig. 9 for the Chamorro and Porté-Agel (2009) cases. An increase in α leads to a significant decrease in λ for a given x, consistent with the model used for λ. In other words, increasing α leads to the λ profile approaching its far-downstream value of 0 faster. In conjunction with Eq. (9), this would at first suggest that the surface shear stress approaches its far-downstream value at a smaller x location for a larger value of α. Fig. 8b however shows the opposite is true, namely the surface shear stress approaches its far-downstream value faster for smaller α. This is explained by the small differences in the evolution of the streamwise velocity at the first off-wall grid point (Fig. 9b) due to different α values. Figure 10 compares the mean streamwise velocity profiles at several downstream locations with the measurements of Li et al. (2021). The LES results for all values of α are in very good agreement with the experimental data. The vertical profiles of turbulence intensity are compared in Fig. 11. For all values of α, the LES and experimental results are in good agreement close to the bottom wall and outside of the IBL. There is some under-prediction of the TI around the IBL height particularly closer to the location of the surface roughness jump. The under-prediction is about 12% at x/δ = 0.84 but reduces to 9% at x/δ = 3.49. Fig. 11 Vertical profiles of T I at different locations downstream of the abrupt surface roughness transition on grid G3 for m = 21 and z 01 = 6.726 × 10 −5 m using the APA wall model with different α compared to the experimental data of Li et al. (2021); 'Expt' stands for experimental data; 'upst' stands for upstream and 'dnst' stands for downstream Figure 12 shows that the height of the IBL is predicted accurately by the LES at all downstream distances. Figures 10, 11 and 12 also indicate that the mean velocity, turbulence intensity and the IBL height are not sensitive to the α value used in the APA model.
To summarize the results of this subsection, our numerical framework is capable of reproducing the results corresponding to the experiments of Li et al. (2021) in addition to those of Chamorro and Porté-Agel (2009). The wall shear stress profiles are sensitive to the α value used in the APA wall model, while other quantities are insensitive to the value of α. Similar agreement between the LES results and the experiments of Li et al. (2021) is obtained on using the BZ wall model. A similar level of sensitivity of the APA wall model results to α is found in simulations corresponding to the data of Chamorro and Porté-Agel (2009). Results for both these are not shown here for brevity.

Sensitivity to Roughness Ratio
Sensitivity to the roughness ratio is studied by analysing simulation results for m = 20, 83.3 and 125 and z 01 = 0.5 mm. The APA wall model with α = 0.027 is used for all the runs. Figure 13a shows that for a configuration with smaller m, the shear stress downstream of the surface roughness jump is larger. This is consistent with intuition since a smaller m = z 01 /z 02 for the same z 01 implies a downstream surface that exerts more resistance to the flow. The change in the downstream wall shear stress between m = 83.3 and 125 is smaller Due to smaller wall shear stress for higher m values, the flow accelerates faster after the roughness jump as seen in Fig. 14a. It is observed that for the largest two values of m studied here, the change inũ 1 is insignificant. The same trends are seen to hold for the T SS and T I profiles (Fig. 14b, c, respectively). A lower value of m indicates a rougher surface after the step jump, which increases the turbulent statistics in the flow. As a result, T SS and T I have larger magnitudes within the IBL for smaller values of m. Similar toũ 1 , the changes in T SS and T I are insignificant between m = 83.3 and 125. Finally, Fig. 13b shows that the IBL height evolution is very similar for all three surface roughness ratios.

Evaluation of IBL Height Models
Several analytical models (Panofsky and Townsend 1964;Chamorro and Porté-Agel 2009;Ghaisas 2020) for the mean velocity downstream of a surface roughness jump as well as the APA wall model require the internal boundary layer height as an input. A number of empirical and/or physics-based models have in turn been developed for the IBL height. Figure 15 compares the IBL height obtained from the LES using the APA wall model with predictions of five different IBL models. A quantitative assessment of the errors between the LES results and the IBL model predictions is shown in Table 4. Details of the models evaluated here are given in Table 3. Some of these models require the solution of a nonlinear equation which is achieved using the 'fsolve' function in Matlab. The prediction of the Elliott (1958) model serves as a good starting guess for the root-finding procedure at each x. The models by Wood (1982) and by Jegede and Foken (1999) agree with the LES results till about x/δ ≈ 2.5 but under-predict the LES results beyond this. These two models show significant errors, between 11 and 15%, for the three surface roughness ratios and depending upon whether the T I or the T SS profiles are used to calculate the IBL heights from the LES results. The model by Panofsky and Dutton (1984) shows similar under-prediction beyond x/δ ≈ 2.5 for m = 83.3 and 125 but shows good agreement throughout the domain for m = 20. This model is accurate for m = 20 (errors less than 6%) but not for the larger values of m (errors more than 10%). In contrast, the Elliott (1958) model shows overall good agreement for the larger two values of m (errors less than roughly 7%) but shows a significant over-prediction for m = 20 (error around 20%). The model by Savelyev and Taylor (2005) shows small over-predictions before, and small under-predictions after, x/δ ≈ 2.5 for m = 83.3 and 125, and small over-predictions throughout the domain for m = 20. The error norms are almost always less than 10%, so this model provides reasonable accuracy over a range of m values. Table 3 for different m values using the APA model on grid G3 for z 01 = 0.5 mm Table 3 List of IBL models evaluated using the LES data Author(s) IBL model Elliott (1958)

Model for Turbulence Intensity
Development of an analytical model for the turbulence intensity downstream of a step change in surface roughness is pursued in this section. To motivate the idea, Fig. 16 shows the T I profiles at different distances downstream of the surface roughness jump along with the profile averaged over the upstream portion. Further, the gray dashed line in Fig. 16 represents the T I profile obtained from a simulation of the flow over a homogeneous surface with roughness z 02 (i.e. with the same surface roughness as in the downstream portion of the heterogeneous case). This simulation is henceforth denoted as 'far-downstream', since it is expected that at sufficiently far downstream of the roughness jump, the flow would have adjusted fully to the new surface conditions (roughness z 02 ) and would have no imprint of the abrupt roughness transition. It is observed that as the downstream distance changes from x/δ = 1 to 3, the T I profile gradually departs from the upstream profile and approaches the far-downstream profile.  The observation in Fig. 16 that the T I profiles downstream of the roughness jump are bounded by the upstream and far-downstream T I profiles is utilized to develop a simple analytical model for the T I . The T I at a downstream location can be expressed as a weighted average of the upstream and far-downstream T I profiles, which can be arranged as: In the above equation, T I upstream and T I f ar−downstream are not evolving with the streamwise distance and are functions only of z since they come from simulations of flow over homogeneously rough surfaces with roughnesses z 01 and z 02 respectively. The empirical, reverse-logarithmic-law model (Stevens et al. 2018) or experimentally obtained profiles can be easily used for these two quantities in place of the simulation results. Figure 17a shows vertical profiles of φ(x, z) extracted from our LES results using Eq. 14 at representative downstream locations of x/δ = 1, 2 and 3 for m = 83.3. As expected, φ is bounded between 0 and 1. A phenomenological model is developed for the weighting function considering it to be dependent on the downstream distance x and the IBL height, Eq. 15 ensures that φ M O DE L goes to 0 at z = δ i and to √ C at z = δ e . For the current study a value of C = 0.8 is taken.
A further correction is required at the first vertical point from the wall, as this point is heavily influenced by the wall model in the LES simulations. Equation 15 is multiplied by 0.85 at the first computational point from the wall. To close this model, we use the Elliott (1958) relation for specifying the IBL height and set α = δ e /δ i = 0.027. Figure 17a shows that this model for φ is in fair agreement with the LES results. In particular, the variation of φ with z below the IBL height is captured well by the model for all downstream locations. Using this modelled profile for φ in Eq. 13 gives a model for T I downstream of a step change in surface roughness. Comparisons between the LES results and the model predictions are shown in Fig. 17b. It is clear that the proposed model predicts the T I very well at different downstream locations. The modelled T I profile shows a small Model for φ is Eq. 15 and for T I is Eq. 13. In each model, δ i is obtained using Elliot's relation kink near the top of the IBL, but the agreement with the LES results over the major portion of the domain is very good. The kink in the T I profiles is related to the fact that φ M O DE L goes to zero at z = δ i with a slope discontinuity, while φ L E S goes to zero smoothly a certain distance above the IBL height. A model for φ that decays smoothly to zero above the IBL height may be constructed in the future to reduce the kink in the T I predictions. The simple model proposed here is in good agreement with the LES results with the maximum relative error between the LES results and model predictions being 4% close to the top of the IBL.
The model is tested against LES data for other roughness ratios as presented in Fig. 18 at a downstream location of x/δ = 3.0. The profiles of φ and T I from the model are seen to be in good agreement with the LES results for these cases as well, indicating that the framework developed here is applicable for a range of surface roughness ratios. Sensitivity of the model for the turbulence intensity to the choice of IBL height model is shown in Fig. 19. Table 5 shows the maximum relative error between the T I predicted by model using different IBL models and the LES results. Using either the Panofsky and Dutton (1984) model or the Savelyev and Taylor (2005) model for the IBL height leads to fairly good predictions of the turbulence intensity profiles at different downstream locations for all three roughness ratios studied here, with the errors ranging from roughly 4-9 %. The kink close  to the IBL height is more pronounced when the IBL height is modelled using the relations proposed by Panofsky and Dutton (1984) and is smallest on using the Elliott (1958) relation.

Conclusion
The flow over a heterogeneously rough surface, with an abrupt change in the aerodynamic roughness, is studied here using large eddy simulations. Simulations are carried out using two wall models (BZ and APA), different upstream and downstream surface roughness values, three ratios of the upstream to downstream roughness (m = z 01 /z 02 ), different grid sizes and different values of the ratio α = δ e /δ i , which is an input to the APA wall model. The LES data are compared with appropriate results from the previously reported wind-tunnel experiments of Chamorro and Porté-Agel (2009) and Li et al. (2021). Different turbulent statistics of the ABL flow are found to be sensitive to the wall models to different extents. Specifically, the wall shear stress and turbulence intensity (T I ) profile show a large sensitivity to the wall model, with the APA model giving larger values for both, and being in better agreement with the experimental results. The mean velocity profile is affected by the wall model to a lesser extent while the profile of the total shear stress (T SS) is almost insensitive to the wall model except for very close to the bottom wall. The internal boundary layer height, defined as the height above the bottom surface above which the upstream and downstream profiles are the same, is largely insensitive to the wall model as well as to the quantity (either T SS or T I ) used to define it.
Our LES framework is validated against the wind-tunnel measurements of Chamorro and Porté-Agel (2009) as well as of Li et al. (2019) where δ/k is appropriately large. The wall shear stress on using the APA model is dependent on the ratio, α, of the equilibrium boundary layer height to the internal boundary layer (IBL) height. Our results show that the Chamorro and Porté-Agel (2009) and Li et al. (2021) measurements are predicted correctly on using α = 0.027 and 0.108, respectively. The other quantities, namely mean velocity, T I , T SS and the IBL height, are insensitive to the value of α and are in fairly good agreement with the experimental results.
The flow is found to be sensitive to the ratio m of the upstream and downstream roughness values. As m is increased, the downstream surface becomes more smooth and exerts lesser drag force on the flow. This leads to smaller surface shear stresses and turbulence intensities as well as to larger acceleration of the flow close to the wall. The IBL heights calculated based on T SS and T I profiles are found to be independent of the surface roughness ratio. Different analytical models for the IBL height evolution are evaluated. The widely-used Elliott (1958) empirical relation is found to be accurate for higher values of m, but is found to over-predict the IBL heights for smaller m. The model proposed by Panofsky and Dutton (1984) is found to be accurate for the smallest m but under-predicts the results for larger m. The Savelyev and Taylor (2005) model is found to be in reasonable agreement with the LES results for the IBL height over a broad range of m values.
An analytical model is proposed for the turbulence intensity downstream of the surface roughness jump. This model predicts the T I as a weighted average between the upstream T I profile and the T I profile far downstream of the surface roughness jump. The weighting function, φ(x, z), is determined by a simple relation and requires the IBL height as an input. Reasonably accurate predictions for the T I are obtained on using any of the three models mentioned above for the IBL height. Nominally, the IBL height given by the Elliott (1958) relation gives good prediction of the turbulence intensity at all downstream locations and over a broad range of surface roughness ratios. The simple TI model proposed here can serve as a building block in, e.g., modelling the TI in wind farms sited on heterogeneously rough surfaces.
The work presented in this paper can be extended along several directions. More experiments and/or wall-resolved LES at surface roughness ratios greater than m = 83.3 and for other z 01 values need to be carried out that will enable development of a methodology of specifying the input α to the APA wall model for a given combination of m and z 01 . The work here focused only on rough-to-smooth transition, and can be extended to smooth-to-rough transitions as well. Finally, studies of surface heterogeneities in the presence of other fea-tures, such as a hill, or one or more wind turbines, as well as other configurations of surface roughness heterogeneities can also be carried out.
Sensitivity of different statistics of the ABL flow to the grid resolution used in the LES are studied first. Besides using two wall models and three different grids, the results are also compared with experimental data of Chamorro and Porté-Agel (2009) wherever available. The surface roughness values are z 01 = 0.5 mm and z 02 = z 01 /83.3. The APA model is utilized with α = 0.027. Figure 20 shows the surface shear stress after the change in surface roughness for different grid sizes using the BZ and APA wall models for m = 83.3. Here, the shear stress at the bottom wall downstream of the surface roughness jump (τ ) is normalized by the surface shear stress upstream of the jump (τ 0 ). The LES data show appreciable change in magnitude when compared between the 128 × 32 × 32 and 192 × 64 × 64 grid cases. An additional simulation for grid size of 240 × 80 × 80 is also compared and it is observed that there is little change in magnitude when compared with the 192 × 64 × 64 grid.
The temporally and spanwise averaged streamwise velocities at two downstream locations (x/δ = 0.5, 1.0) after the roughness jump are presented in Fig. 21. These profiles are almost insensitive to the grid resolution for both models. Small differences are seen close to the top of the domain, where the velocity profiles are seen to agree better with the upstream logarithmic law, Eq. (1), with increasing grid resolution. Closer to the bottom boundary, the velocity accelerates due to the reduced surface roughness. This acceleration is the same for all grids for the BZ as well as APA wall models. The total shear stress (T SS) and the streamwise turbulence intensity (T I ) are shown at two downstream locations for the APA wall model in Fig. 22. The T I increases when going from the coarsest to the intermediate grid, but is unchanged over the two finest grids employed here. This indicates that a computational grid with 240 × 80 × 80 points is sufficient to obtain grid-independent results. Consequently, all simulations reported in Sect. 3 utilize these many grid points.