1. Introduction
When an atmospheric boundary layer flows over a surface with different properties, be it the roughness, temperature, moisture or all of the aforementioned factors, an internally developing boundary layer will be produced. Meteorological examples include wind flowing over different types of vegetation/media, surface morphologies (rural to urban or vice versa), etc. An overview of internal boundary layers in the context of field measurements, laboratory experiments, numerical simulations and analytical models was provided by the seminal paper of Garratt (Reference Garratt1990). Studying the statistical properties of the internally developing flow has implications for applications such as wind-energy utilisation and meteorological predictions in situations where there is a sudden change in surface properties.
To investigate the effect of surface aerodynamic properties’ variation on the boundary layer development, a commonly used model is a surface constructed with a streamwise step change in roughness. Developed from the discontinuity in surface roughness, the flow inside a layer, i.e. the internal boundary layer (IBL), will adjust its characteristics to the underlying new surface, while the flow aloft retains the upwind characteristics. The produced IBL usually includes two layers: the one just beneath its upper edge is in a non-equilibrium state, and the one below contains flows fully adapted to the underlying surface condition.
A large amount of previous work about the IBL focuses on its growth rate and on models to predict its interface with the outer flow. The thickness of the IBL, 
$\delta _i$, is found to grow with fetch 
$x$ as a power function with an exponent 
$n$, namely, 
$\delta _i \approx x^n$. The power exponent 
$n$ varies within the range of 
$[0.2,0.8]$ in the existing literature, and it is dependent on the roughness arrangement and sensitive to the estimation procedure. The approaches to identify the IBL can be categorised into two branches according to the variables of interest, which are the mean streamwise velocity (Elliott Reference Elliott1958; Antonia & Luxton Reference Antonia and Luxton1971a; Cheng & Castro Reference Cheng and Castro2002; Gul & Ganapathisubramani Reference Gul and Ganapathisubramani2022) and normal stress (Efros & Krogstad Reference Efros and Krogstad2011; Li et al. Reference Li, de Silva, Chung, Pullin, Marusic and Hutchins2021). Thus, recent work applied distinctive identification procedures to compare the growth curves of IBLs (see Rouhi, Chung & Hutchins (Reference Rouhi, Chung and Hutchins2019), Bou-Zeid, Meneveau & Parlange (Reference Bou-Zeid, Meneveau and Parlange2004), Sessa, Xie & Herring (Reference Sessa, Xie and Herring2018) and Li et al. (Reference Li, de Silva, Chung, Pullin, Marusic and Hutchins2021) amongst others). The growth curves of IBLs for a long fetch with 
$n$ close to 0.8 were anticipated by the original model proposed by Elliott (Reference Elliott1958) and by the diffusion analogue in Miyake (Reference Miyake1965), which were further developed by Savelyev & Taylor (Reference Savelyev and Taylor2005). However, slower growth of the IBLs (
$n<0.8$) is generally observed when there is a change from a rough to a smooth surface (Antonia & Luxton Reference Antonia and Luxton1972) as well as from a smooth (or smoother than the downstream surface) to protruding-roughness surface (Antonia & Luxton Reference Antonia and Luxton1971a; Cheng & Castro Reference Cheng and Castro2002; Lee Reference Lee2015; Sessa et al. Reference Sessa, Xie and Herring2018; Ding et al. Reference Ding, Placidi, Carpentieri and Robins2023), which is troublesome to explain/reproduce with these early analytical models. Progress in interpreting the appearance of small 
$n$ values includes the study of a flow from a rough to a smooth surface by Li et al. (Reference Li, de Silva, Chung, Pullin, Marusic and Hutchins2022), where they incorporated the finite IBL thickness and thus included the shear stress decay into the model of Elliott (Reference Elliott1958). In addition, Ding et al. (Reference Ding, Placidi, Carpentieri and Robins2023) incorporated the decay laws of the diffusion and the vertical advection terms into the original diffusion model to interpret the small exponent 
$n$ that arises for both neutrally and stably stratified boundary layers developing from a rough to a much rougher surface.
In contrast with their growth rate, the dynamics of IBLs has been studied less. Quadrant analysis, describing the contribution to the shear stress from four quadrants of the streamwise and wall-normal velocity fluctuations (
$u$ and 
$w$, respectively) serves as a powerful tool to study the dynamics of turbulent boundary layers (Lu & Willmarth Reference Lu and Willmarth1973). Through studying the recovery of a developing flow from the two-dimensional (2-D) square-cylinder ribbed surface to a smooth surface, Ismail, Zaki & Durbin (Reference Ismail, Zaki and Durbin2018) found that, among four quadrant events, ejections 
$(u<0,w>0)$ recover at the fastest rate followed by sweeps 
$(u>0,w<0)$. Gul & Ganapathisubramani (Reference Gul and Ganapathisubramani2022) found that the contribution from the second and fourth quadrants can be affected significantly by the step change in surface roughness. In their cases, from a sandpaper roughness to a smoother surface, ejection events were observed to diminish and sweep events enhance. Near the wall, ejection (sweep) events first overshoot (undershoot) and then undershoot (overshoot) before reaching their equilibrium state when the streamwise locations of the two surfaces were inverted.
Hanson & Ganapathisubramani (Reference Hanson and Ganapathisubramani2016) and Li et al. (Reference Li, de Silva, Chung, Pullin, Marusic and Hutchins2021) studied the power spectrum density (PSD) of streamwise velocity fluctuations during the recovery of a rough flow across a smooth surface. Li et al. (Reference Li, de Silva, Chung, Pullin, Marusic and Hutchins2021) found that the difference in the PSD between the developing flow and the fully developed flow on the smooth wall showed a footprint of energy excess. Such energy excess arose in restricted time scales ranging from the Kolmogorov time scale to the time scale associated with the large-scale motion aloft (Ismail et al. Reference Ismail, Zaki and Durbin2018).
Although considerable progress has been made in understanding the recovery of the equilibrium layer over smooth surfaces, the dynamics related to the development of the non-equilibrium layer over a rough surface remains less explored. Other than the aforementioned quadrant analysis by Gul & Ganapathisubramani (Reference Gul and Ganapathisubramani2022), another paramount contribution is from Antonia & Luxton (Reference Antonia and Luxton1971a), who performed turbulent energy balance analysis in the IBL. One of their key findings is that the turbulent diffusion term reaches a peak just beneath the IBL, playing a key role in the growth of the IBL and reflected in the large skewness of wall-normal velocity observed in this layer.
Thermal stability, an inevitable ingredient in atmospheric fluid motion in nature, modifies the near-wall structure, and hence the turbulent properties of the boundary layer. Due to the suppression of the vertical mixing by thermal stability (Williams et al. Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017; Marucci, Carpentieri & Hayden Reference Marucci, Carpentieri and Hayden2018), it was found that the IBL becomes much shallower when compared with neutrally stratified boundary layers (Sessa et al. Reference Sessa, Xie and Herring2018; Ding et al. Reference Ding, Placidi, Carpentieri and Robins2023), hence the slow growth rates of stable IBLs. However, the statistical properties during the development of stably stratified IBLs remain largely unexplored and their dynamics is poorly understood.
This study aims to understand the dynamics of IBLs induced by an abrupt change in surface roughness, and especially the impact of thermal stability (see the experimental set-up in § 2). Firstly, the skin-friction coefficient for stably stratified boundary layers is studied in § 3.1, followed by the outer-layer similarity in § 3.2. Then, by investigating high-order moments of the velocity in § 3.3, a non-adjusted region characterised by strong skewness of wall-normal velocity is discussed. Quadrant analysis is performed in § 3.4. The results associated with thermal properties are given in § 3.5. Conclusions are presented in § 4.
2. Experimental set-up and parameters
2.1. Facility and boundary layer generation
The experiments were conducted in the EnFlo wind tunnel, which is a suck-down open-circuit wind tunnel with a working section size of 
$1.5\ {\rm m} \times 3.5\ {\rm m} \times 20\ {\rm m}$ in the wall-normal, spanwise and streamwise directions, respectively. The wind tunnel is capable of simulating non-neutrally vertically stratified boundary layers. These are achieved through 15 vertical levels of independent heaters (total 450 kW) at the inlet section. In addition, the floor of the wind tunnel is made of 20 panels connected to a recirculating chilled water system and can be cooled/heated independently.
The tunnel floor was covered by two types of three-dimensional (3-D) roughness elements, both of which were arranged in staggered patterns, as shown in figure 1. The first 11 metre fetch from the inlet section was covered by the roughness element ‘
$A$’ of height 
$h_1=16$ mm. The following 9 metre rougher surface was made by using roughness element ‘
$B$’, with a height 
$h_2=20$ mm. Further details on their dimension and spacing are discussed in Ding et al. (Reference Ding, Placidi, Carpentieri and Robins2023), where this set-up was first used. The resulting roughness lengths, evaluated with the methodology described in § 2.3, are 
$z_{01}=0.1\ {\rm mm}$ for the upstream and 
$z_{02}=1\ {\rm mm}$ for the downstream surface in neutral conditions, respectively. The coordinate system is defined in figure 1, where 
$x$, and 
$z$ represent the streamwise and wall-normal directions, respectively. The origin of this system is located at the intersection of the centre line of the wind tunnel and the first row of the downstream roughness, on the tunnel floor (i.e. at the location of the discontinuity in surface roughness). The homogeneity of flows in the spanwise direction is rigorously verified in Appendix B.
Figure 1. Schematic of the experimental set-up. Symbols referring to the measurement locations of the vertical profiles are used consistently throughout the paper.
Thick boundary layers were generated by 13 Irwin spires of height of 600 mm spaced uniformly at the working section inlet. This is a widely used technique to scale the atmospheric boundary layer in the laboratory. The reference wind speed 
$U_{ref}$ at 
$z=1$ m and 5 m downstream from the inlet section was monitored by a sonic anemometer. To generate stably stratified boundary layers, the incoming flow was heated at a prescribed temperature profile at the inlet and developed over the cooled floor which is chilled to a constant value of surface potential temperature, 
$\varTheta _0$, except for the first 5 m. The uncooled fetch is necessary for the incoming flow to reach equilibrium at short fetch and to obtain smoothly varying boundary layer quantities (within 9 m, see Hancock & Hayden (Reference Hancock and Hayden2018), who utilised that same set-up).
The measurement of streamwise and wall-normal velocities (
$\tilde {u},\tilde {w}$) was performed using a two-component laser Doppler anemometer (LDA). Two probes of diameter 27 mm and a focal length of 160 mm were used. Mean and fluctuating absolute temperatures were recorded by a thermistor and a cold wire, respectively, which were mounted onto the same traverse as the LDA probe; then corresponding potential temperatures, 
$\varTheta$ and 
$\theta$ were calculated. The cold wire was placed 4 mm downstream of the LDA measurement volume (estimated to be 1.57 mm long in the spanwise direction and with a diameter of 0.074 mm) and the thermistor was placed adjacent to it (15 mm) to avoid any interference with the LDA, as described by Marucci et al. (Reference Marucci, Carpentieri and Hayden2018). The bias error in the mean velocity is within 
$\pm 1\,\%$, and 
$\pm 0.1\,\%$ for 
$\varTheta$. For any cross- or self-correlation between 
$u$, 
$w$ and 
$\theta$, the error is within 
$\pm 7\,\%$. The measured locations of vertical profiles are presented as vertical dashed lines in figure 1, while the symbols on the top of these profiles are used throughout the paper to denote their streamwise positions.
2.2. Cases studied
Four different cases are considered in this work. Two neutral boundary layers with Reynolds number 
$Re_{\delta }=4.5\times 10^4$ and 
$Re_{\delta }=3.0\times 10^4$, achieved by varying the reference wind speed, were studied. Here, 
$Re_{\delta }=U_\infty \delta _0/\nu$, with 
$U_\infty$ denoting the free-stream velocity at a height of 
$800\ {\rm mm}$ from the floor, 
$\delta _0$ being the depth where the mean streamwise velocity 
$U$ reaches 
$0.99U_{\infty }$ and 
$\nu$ being the kinematic viscosity. Two stably stratified boundary layers with different degrees of thermal stability were investigated, as measured by the bulk Richardson number
\begin{equation} Ri_b=\frac { g (\varTheta_{\delta}-\varTheta_0)\delta_0}{\varTheta_0 U^2_{\infty}}, \end{equation}
where 
$g$ denotes the gravitational acceleration and 
$\varTheta _0$ (
$\varTheta _{\delta }$) the mean potential temperature of the floor (the top of the boundary layer). For the weakly stable case, 
$Ri_b=0.13$ and 
$Re_{\delta }=4.5\times 10^4$, whereas, for the moderately stable case, 
$Ri_b=0.27$ was achieved by reducing the wind speed, giving 
$Re_{\delta }=3.0\times 10^4$, providing 
$Re_{\delta }$-matched conditions in stable and neutral boundary layers. Further boundary layer properties are reported in table 1.
Table 1. Summary of experimental parameters. The equivalent sand grain roughness is determined by 
$k_s^+=({z_0u_*}/{\nu })\textrm {e}^{\kappa A_{AF}}$, with 
$A_{AF}=8.5$ (Jiménez Reference Jiménez2004).
2.3. Roughness length and friction velocity determination
Through varying the wind speed and thermal stability, the aerodynamic properties of the upstream and downstream surfaces are altered accordingly. The friction velocities 
$u_*$ of the upstream and downstream surfaces are determined by linearly extrapolating 
$\overline {uw}$ within the logarithmic region 
$0.1< z/\delta <0.2$ to the ground for both neutral and stable conditions. Here, the overbar denotes time averaging. The incoming flow is in a state of equilibrium, as discussed in Ding et al. (Reference Ding, Placidi, Carpentieri and Robins2023), and its friction velocity, 
$u_{*1}$, is determined at 
$x=-0.76\ {\rm m}$. With regards to the downstream surface, the friction velocity 
$u_{*2}$ is determined by the mean value of the local friction velocity in the range 
$2.28\ {\rm m}< x< 5.88\ {\rm m}$ based on the observation that the local friction velocity becomes invariant with fetch for 
$x>2.28\ {\rm m}$. According to the study of Placidi & Ganapathisubramani (Reference Placidi and Ganapathisubramani2015), the discrepancy of using the least-square fit procedure to determine 
$u_*$ when compared with that obtained using direct measurements is around 
$10\,\%$. In neutrally stratified boundary layers, the roughness length 
$z_0$ was determined by least-squared fitting the mean streamwise velocity 
$U$ in the logarithmic region (
$0.1< z/\delta <0.2$) using the expression 
$U=u_*/\kappa \ln {(z-d)/z_0}$, with the von Kármán constant 
$\kappa =0.41$, which is applicable in rough-wall boundary layers (Snyder & Castro Reference Snyder and Castro2002). The zero-plane displacement 
$d=0$ arises from the sparseness of the roughness elements regardless of thermal stability. The rough surfaces in figure 1 resemble that of Snyder & Castro (Reference Snyder and Castro2002), where the roughness Reynolds number, 
$Re_* =z_0u_*/\nu$, for the upstream roughness is 
$0.4< Re_*<0.6$ and 
$Re_*>1$ for the downstream roughness (neutral cases 1 and 3). Thus, the approach flow is close to an aerodynamically smooth boundary layer and the downstream flow is fully rough.
For the stable case, the similarity in the surface layer is dictated by the Monin–Obukhov similarity theory (MOST) (Monin & Obukhov Reference Monin and Obukhov1954), which prescribes the vertical profile of the mean streamwise velocity in the surface layer
\begin{equation} U(z) = u_*/\kappa \left[ \ln(z/z_0)+\varPsi_m(z) \right]\!. \end{equation}
According to the MOST, 
$\varPsi _m(z)=\beta _m(z-z_0)/L_0$, which measures the modification to the law of the wall. Here, 
$L_0$ is the surface Obukhov length, that is, 
$L_0=u_*^2\varTheta _0/(g\kappa \theta _*)$. In the expression, 
$\theta _*$, the friction temperature, is determined by the ratio of the wall heat flux to the friction velocity, i.e. 
$\overline {w\theta }_0/u_*$. The parameter 
$\beta _m$ is a constant value, which was found to be 8 in the case studied (Hancock & Hayden Reference Hancock and Hayden2018). The roughness length 
$z_0$ in stable cases was determined by fitting the above vertical profiles of 
$U$ in the surface layer instead. The parameters describing the surface properties are listed in table 1.
Two different methodologies were employed to determine the friction velocity as it responds to the change in roughness. The first methodology, based on the model in Elliott (Reference Elliott1958), has been widely used to determine the response of the wall shear stress to a roughness change (Gul & Ganapathisubramani Reference Gul and Ganapathisubramani2022; Li et al. Reference Li, de Silva, Chung, Pullin, Marusic and Hutchins2022). This model assumes that the flow within the IBL reaches equilibrium promptly, hence there is a constant roughness length of the downstream surface 
$z_{02}$, while the mean streamwise velocity assumes the form of a piecewise function
\begin{equation} U(z)= \begin{cases} \dfrac{u_{*1}}{\kappa}\ln(z/z_{01}) , & z \ge \delta_i, \\[10pt] \dfrac{u_{*}}{\kappa}\ln(z/z_{02}) , & z \le \delta_i. \end{cases} \end{equation}
For the continuity of 
$U$ at 
$z=\delta _i$, the relation between 
$u_*$ and 
$\delta _i$ needs to be
\begin{equation} u_*=u_{*1}\ln(\delta_i/z_{01})/\ln(\delta_i/z_{02}). \end{equation}Here, we extend this methodology to stably stratified boundary layers, under the two following assumptions. Firstly, the shallow IBL is assumed to be located within the surface layer. Secondly, MOST is assumed for velocity and temperature, in the forms of piecewise functions
\begin{gather} U(z)= \begin{cases} \dfrac{u_{*1}}{\kappa}\left[\ln(z/z_{01})+\dfrac{\beta_m(z-z_{01})}{L_{01}}\right], & z \ge \delta_i, \\[10pt] \dfrac{u_{*}}{\kappa}\left[\ln(z/z_{02})+\dfrac{\beta_m(z-z_{02})}{L_{02}}\right], & z \le \delta_i, \end{cases} \end{gather}
\begin{gather} \varTheta(z)= \begin{cases} \dfrac{\theta_{*1}}{\kappa}\left[\ln(z/z_{0h1})+\dfrac{\beta_h(z-z_{0h1})}{L_{01}}\right], & z \ge \delta_i, \\[10pt] \dfrac{\theta_{*}}{\kappa}\left[\ln(z/z_{0h2})+\dfrac{\beta_h(z-z_{0h2})}{L_{02}}\right], & z \le \delta_i. \end{cases} \end{gather}
The kinematic (thermal) roughness length, 
$z_{02}~(z_{0h2})$ downstream of the step is determined at the furthest measurement location and kept invariant during the fitting procedure, so was the downstream Obukhov length 
$L_{02}$. Here, 
$\beta _m=8$ and 
$\beta _h=16$ (Hancock & Hayden Reference Hancock and Hayden2018). For each 
$x$, the vertical fitting regions are divided into two sections by the thickness of IBL, which is identified as the height of local variance of the streamwise velocity beginning to deviate from its counterpart upstream of the roughness change (Ding et al. Reference Ding, Placidi, Carpentieri and Robins2023). The second common methodology (Cheng & Castro Reference Cheng and Castro2002) deduces the friction velocity from 
$\overline {uw}$, as mentioned in the context of equilibrium boundary layers. The comparison of the obtained friction velocity from two methodologies is conducted in § 3.1 and some complementary results are present in Appendix A.
3. Results and discussion
3.1. Skin-friction coefficient
For a boundary layer in equilibrium, the skin-friction coefficient 
$C_f=2\rho _w u_*^2/(\rho _\infty U_{\infty }^2)$, where 
$\rho _w, \rho _{\infty }$ denote the density at the wall and in the free stream, is found to decrease with an increasing degree of thermal stability and Reynolds number (Williams et al. Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017). This is also the case in the current work. Throughout the paper, 
$C_f=2u_*^2/U_{\infty }^2$ is used for the neutrally stratified boundary layers. Herein, we calculated 
$\rho _w$ using the surface temperature, which was monitored by thermistors attached to the tunnel floor. The floor temperature is 
$15\,^{\circ }{\rm C}$ for case 2 and 
$14\,^\circ {\rm C}$ for case 4, leading to 
$\rho _w=1.225\ {\rm kg}\ {\rm m}^{-3}$ for case 2 and 
$\rho _w=1.221\ {\rm kg}\ {\rm m}^{-3}$ for case 4.
3.1.1. Skin-friction coefficient development
Figure 2(a) shows the behaviour of the skin-friction coefficient ratio 
$C_f/C_{f1}$ across the roughness change, where 
$C_{f1}$ represents the reference value at 
$x=-0.76$ m. An overshoot of 
$C_f/C_{f1}$ to the roughness change takes place for all cases before the skin-friction coefficient adjusts to the new underlying surface and its degree is enhanced slightly by thermal stratification. This overshoot is more pronounced in the results obtained by Elliott's methodology (solid symbols) when compared with those determined from 
$-\overline {uw}$ (empty symbols), especially very near the roughness change. Here, we discuss the accuracy of the two methodologies in determining 
$u_*$ following the roughness change. Li et al. (Reference Li, de Silva, Chung, Pullin, Marusic and Hutchins2022) found the predicted behaviour of 
$u_*$ across the transition in roughness with Elliott's methodology to be 90 % accurate when compared with direct measurements provided that (i) the thickness of the IBL was less than 0.6
$\delta$ and (ii) that the point of interest was far enough from the roughness change (
$x^+>5000$). The measurement region in our study is from 1.32 m (1.8 m) onward for case 1 (3), which meets the requirement of 
$x^+>5000$. Therefore, we expect the estimation of the friction velocity for neutral cases at 
$x>1.32$ m to be acceptable when using Elliott's model. On the other hand, using a linear extrapolation to determine 
$u_*$ after the roughness change is highly dependent on the existence of an equilibrium logarithmic layer. Rouhi et al. (Reference Rouhi, Chung and Hutchins2019) reported that the recovery region of the boundary layer ‘reaches the beginning of the logarithmic layer after a fetch of 2.5
$h$’ for a transition from smooth to rough surfaces, implying that the recovery of the logarithmic layer takes place after 2.5 times the boundary layer thickness. This adjustment region is – no doubt – roughness-property dependent. In figure 2 of Ding et al. (Reference Ding, Placidi, Carpentieri and Robins2023), we observed that the mean streamwise velocity profiles in the lower region (in the range 
$0.1< z/\delta <0.2$) collapse onto each other after 
$x=2.28$ m; this fetch for log-law recovery is consistent with that in Rouhi et al. (Reference Rouhi, Chung and Hutchins2019). Furthermore, the difference between the estimation of the skin-friction coefficient from these two methodologies nearly vanishes around 
$3\delta _0$, as shown in the inset of figure 2(a). This suggests that the estimation of 
$u_*$ by the extrapolation method enables a reasonable evaluation of the friction velocity in this region. It must be noted, however, that, near the roughness change, the actual value of 
$u_*$ is difficult to determine accurately – as highlighted by the discrepancy between the methodologies adopted herein – but it is expected to be somewhere between the two sets of estimated values. We fully acknowledge the uncertainty in determining this quantity in the near-step-change region, but this is not the focus of the work here.
Figure 2. Variation of the skin-friction coefficient (a) with fetch, and (b) with the displacement thickness. The empty symbols in (a) are from the extrapolation method and solid symbols are from Elliott's method. The inset in (a) shows the difference of 
$C_f$ between two methodologies. In (b), the darkness of colour decreases with fetch, and the black solid and dotted curves are generated from (3.6) with 
$\varPi =0.55$ and 
$\varPi =0.7$, respectively. The blue and light blue curves are from (3.6) with 
$K=3.88$ for 
$Ri_b=0.13$ and 
$K=7.85$ for 
$Ri_b=0.27$, respectively.
3.1.2. Skin-friction coefficient as a function of displacement thickness
The relation between 
$C_f$ and the displacement thickness, 
$\delta ^*$, or the momentum thickness, 
$\theta$, customarily defined as 
$\delta ^*=\int _0^{\infty }(1-U/U_{\infty }) \,{\rm {d}}z$ and 
$\theta =\int _0^{\infty }U/U_{\infty }(1-U/U_{\infty })\,{\rm d}z$, respectively, reflects properties of developing boundary layers. Clauser (Reference Clauser2002) and Rotta (Reference Rotta1962) derived an analytical relationship between 
$C_f$ and 
$\theta$ with neutral-thermal stratification, which was revisited by Castro (Reference Castro2007) for its application on fully rough boundary layers. Here, we extend the methodology in Castro (Reference Castro2007) to stably stratified conditions. Two assumptions associated with the similarity in the surface and outer layers are applied.
Firstly, given that the MOST is applicable within the surface layer, using the error function to restrict the effective region of MOST modifies the gradient of 
$\varPsi _m(z)$ into
\begin{equation} \frac{{\rm d}\varPsi_m}{{\rm d}z}=\frac{\beta_m}{L_0}\left[0.5-0.5\,{\rm{erf}}\left(\frac{z-z_m}{s}\right)\right], \end{equation}
where the parameter 
$z_m$ denotes the height of the surface layer and 
$s$ measures the height of the transition region from the surface layer to the outer layer. In our cases, 
$z_m/L_0\approx 0.2$. Secondly, in the outer region of the weakly stable case, the modification is so small that its wake function has nearly the same form as that of the neutrally stratified boundary layer. Using the general wake function 
$\mathscr {W}$ (derived empirically by Coles (Reference Coles1956) in table 1 of his paper), the complete profile of the mean streamwise velocity reads
\begin{equation} U(z) = \frac{u_*}{\kappa} \ln{\frac{z}{z_0}}+ \frac{u_*\beta_m}{\kappa L_0} F(z,z_m,s)+\frac{u_*}{\kappa}\varPi \mathscr{W}(z/\delta), \end{equation}
where 
$F(z,z_m,s)=0.5\int _{z_0}^z(1-\textrm {erf}({z-z_m}/{s})) \,\textrm {d}z$, with the parameter 
$\varPi$ representing the strength of the wake function. The above formulation can return to that of neutral conditions when the second term on the right-hand side vanishes owing to the Obukhov length 
$L_0$ approaching infinity.
When 
$z=\delta$ in (3.2), the velocity at the top of the boundary layer (BL) reads
\begin{equation} U_e = \frac{u_*}{\kappa} \ln{\frac{\delta}{z_0}}+ \frac{u_*\beta_m}{\kappa L_0} F_{\delta}+\frac{u_*\varPi}{\kappa} \mathscr{W}(1), \end{equation}
where 
$F_{\delta }$ denotes the value of the function 
$F(z,z_m,s)$ at the top of the BL, that is, 
$F_{\delta }\equiv F(\delta,z_m,s)$. Thus, the velocity profile in the defect form is given by
\begin{equation} \frac{U_e-U}{u_*} =-\frac{1}{\kappa} \ln{\frac{z}{\delta}}+ \frac{\beta_m}{\kappa L_0} \left[F_{\delta}-F(z)\right]+\frac{\varPi}{\kappa}\left[\mathscr{W}(1)-\mathscr{W}(z/\delta)\right], \end{equation}
where the wake function 
$\mathscr {W}(z/\delta )$ has the same boundary and normalising conditions as those in the neutral case, namely, 
$\mathscr {W}(1)=2$, 
$\int _0^{1}\mathscr {W}d(y/\delta )=1$ and 
$\mathscr {W}(0)=0$. Figures 3(e) and 3(g) show the theoretical curves of velocity defect for the stably stratified BLs. The integral of the aforementioned equation establishes the relationship between the displacement thickness and the BL thickness, namely,
\begin{equation} \frac{U_e\delta^*}{u_*\delta} = \frac{1}{\kappa}+\frac{2\varPi}{\kappa}+ \frac{\beta_m}{\kappa L_0}\left[F_{\delta}-\int_0^1F(z)\,\textrm{d}(z/\delta)\right]={\rm const}. \end{equation}
Replacing 
$\delta$ in (3.3) with its function of 
$U_e/u_*$ and 
$\delta ^*$ in (3.5), after rearrangement, the relation of 
$C_f$ and 
$\delta ^*$, now for any thermal stratification, is prescribed by
\begin{equation} \frac{\delta^*}{z_0}=\sqrt{\frac{\rho_\infty}{\rho_w} \frac{C_f}{2}}\exp\left({\kappa\left(\sqrt{\frac{2}{C_f}\frac{\rho_w}{\rho_\infty}}-K\right)}\right)\!, \end{equation}
with the constant value 
$K={2\varPi }/{\kappa }-({1}/{\kappa })\ln {({1+\varPi +\tilde {F}\beta _m/L_0})/{\kappa }}+{\beta _m F_{\delta }}/{\kappa L_0}$ and 
$\tilde {F}=F_{\delta }-\int _0^1F(z,z_c,s)\,\textrm {d}(z)$.
Figure 3. Velocity defect 
$(U_\infty -U)^+$ as a function of (a,c,e,g) 
$z/\delta$ and (b,d,f,h) 
$z/\varDelta$ at several streamwise locations. The inset plots in (a,c,e,g) only show the data outside the IBL for 
$x>0$ and the complete profile at 
$x=-0.76\ {\rm m}$. The inset plots in (b,d,f,h) only show the data within the IBL for 
$x>0$ and the complete profile at 
$x=-0.76\ {\rm m}$. (a,b) Case 1; (c,d) case 3; (e,f) case 2; (g,h) case 4. The solid curves represent the theoretical expression given by (3.4) with the wake strength 
$\varPi =0.70$ for (a–d), 1.0 for (e,f) and 1.2 for (g,h). Colours/symbols are defined in figure 1. The local 
$u_*$ is determined from Elliott's model.
Note that the above equation is also applicable in neutrally stratified BLs with 
$K={2\varPi }/{\kappa }-({1}/{\kappa })\ln {(({1+\varPi })/{\kappa })}$ and density ratio 
${\rho _\infty }/{\rho _w}=1$ as 
$L _0 \to \infty$. Under neutral conditions, the formulation given in (3.6) recovers to that in (1.9) in Castro (Reference Castro2007) via (1.6) in the same paper.
Let us now apply this newly developed formulation to our data. In the neutrally stratified boundary layers (figure 2b), 
$C_f$ in the incoming flow faithfully follows (3.6) with 
$\varPi =0.55$. A sharp upward deviation takes place after the step change in surface roughness, which is expected when only parts of BLs have adjusted to the roughness change. For the stably stratified BLs, the curves of 
$C_f(\delta ^*)$ fall below the neutral cases due to the modulation of the shear stress by thermal stability; this discrepancy is enlarged by the thermal stability. This thermal stability effect on 
$C_f(\delta ^*)$ is well predicted by (3.6).
3.2. Outer-layer similarity
3.2.1. Velocity defect profiles
Velocity defect profiles have been widely studied to assess outer-layer similarity (Castro Reference Castro2007; Efros & Krogstad Reference Efros and Krogstad2011). Figure 3 shows the velocity defect scaled by the local friction velocity 
$(U_\infty -U)^+ \equiv (U_\infty -U)/u_*$ as a function of 
$z/\delta$. In the developing neutrally stratified flow in figures 3(a) and 3(c), the velocity defect profiles fall below that of the upstream flow (
$x=-0.76$ m). With an increase of 
$x$, the lower part of the profiles (within the IBL) lifts up gradually and approaches the counterpart upstream of the roughness change. Meanwhile, data beyond the IBL depth remain divergent from those of the upstream flow. The local friction velocity is not an adequate scaling parameter for the entire layer, the upstream friction velocity 
$u_{*1}$ should be employed for depicting the outer-layer similarity. As shown in insets of figures 3(a) and 3(c), 
$(U_\infty -U)/u_{*1}$ as a function of 
$z/\delta$ outside the IBL indeed collapses well onto the upstream reference curve. This is also the case for the stably stratified flows in figures 3(e) and 3(g).
Regarding the velocity defect within the IBL, the vertical length scale 
$\varDelta =\delta^* U_{\infty }/u_*$ (first proposed by Castro (Reference Castro2007) to describe the universal defect law for different types of rough-wall BLs) improves the collapse of data within IBLs for all cases, as shown in figures 3(b), 3(d), 3(f) and 3(h). The variation of 
$U$ within the IBL is proportional to the variation of the wall shear stress (Antonia & Luxton Reference Antonia and Luxton1971b); for this reason, 
$\varDelta$ incorporating the local friction velocity counteracts the change in 
$U$. Thus, 
$\varDelta$ acts as a better scale than 
$\delta$ for describing the universal defect law within the IBL. The quality of the data collapse in figure 21 shows that the scaling is largely insensitive to the methodology adopted for calculating 
$u_*$ for data 
$x>1.32$ m as the difference of 
$u_*$ for two methodologies vanishes at 
$x=2.28$ m (see figure 2). However, for stable cases, the improvement of collapse by using 
$\varDelta$ is not as significant when compared with neutral cases since 
$\varDelta$ does not incorporate the Obukhov length, which becomes one of the dominant length scales for stably stratified flows. A more appropriate length scale to describe the defect law in stable cases needs to be further studied.
3.2.2. Diagnostic plots
Figure 4 shows the response of diagnostic plots (
$\sigma _u/U$ as a function of 
$U/U_\infty$, where 
$\sigma _u$ represents the standard deviation of the streamwise velocity) to the step change in surface roughness. In the neutral approach flow (
$x=-0.76\ {\rm m}$), 
$\sigma _u/U$ in the outer layer (
$0.6< U/U_\infty <0.98$) follows the smooth-wall asymptote of Alfredsson, Segalini & Örlü (Reference Alfredsson, Segalini and Örlü2011) described by 
$\sigma _u/U=0.286\unicode{x2013}0.255U/U_\infty$. Given the sparse and small roughness characterising the upstream approach flow, it is within expectations that the diagnostic plots are close to those of a smooth surface.
Figure 4. Diagnostic curves for (a) case 1, (b) case 3, (c) case 2 and (d) case 4 at several streamwise locations. The dashed and the solid lines correspond to the smooth (Alfredsson et al. Reference Alfredsson, Segalini and Örlü2011) and fully rough (Castro et al. Reference Castro, Segalini and Alfredsson2013) asymptotes, respectively. The dotted lines in (c,d) are fits to the data in the forms of 
$\sigma _u/U=0.322\unicode{x2013}0.335U/U_\infty$ and 
$\sigma _u/U=0.314\unicode{x2013}0.342U/U_\infty$, respectively. Colours/symbols are defined in figure 1.
Adjusted to the rougher surface, the slope of the diagnostic curve in the lower region of the outer layer becomes 
$-0.38$, meanwhile this region gradually expands in height with fetch. The modification of the slope ceases at the height of the IBL. While the adjusted slope in the lower region aligns with that in the fully rough flows (solid line), its magnitude is lower compared with a fully rough asymptote (Castro, Segalini & Alfredsson Reference Castro, Segalini and Alfredsson2013). This discrepancy may be attributed to the relatively small Reynolds number. Gul & Ganapathisubramani (Reference Gul and Ganapathisubramani2021) studied the effect of the Reynolds number on the diagnostic curves by varying 
$Re_\delta$ from 
$3\times 10^4$ to 
$1.1\times 10^5$ on sandpaper-rough surfaces and found that all of them remain lower than the curve for fully rough flow found in Castro et al. (Reference Castro, Segalini and Alfredsson2013), even though the diagnostic curves become higher with increasing 
$Re_\delta$. In our study, 
$Re_\delta$ is around 
$10^4$, which is comparable.
The linear relation in the diagnostic plots modified to account for the thermal stability of the BL is shown in figure 4(c,d). In the regime of 
$U/U_\infty <0.75$ (
$z/\delta _0<0.2$), the diagnostic curves have a slope of 
$-0.255$, the same as that in neutral conditions, but the whole curve shifts downwards by an amount that increases with the degree of thermal stability. Moreover, such a linear relation terminates at 
$U/U_\infty =0.75$ and is seemingly replaced with a new linear relation with a smaller slope of 
$-0.171$ (
$-0.117$) in the upper region of the outer layer (
$0.75< U/U_\infty <0.96$) for 
$Ri_b=0.13$ (0.27). In contrast to the gradual change with fetch in neutral flow, the diagnostic curves for stable cases overshoot the asymptote following the roughness change, and then approach equilibrium states for longer fetch. Within the surface layer (
$U/U_\infty <0.75$), the linear relations at 
$x=5.88\ {\rm m}$ for both stable cases have a slope of 
$-0.34$, which is slightly smaller than that for the neutrally stratified flows. This linear relation intercepts with 
$U/U_\infty =0$ at around 0.322 (0.314) for 
$Ri_b=0.13$ (
$0.27$), which is much smaller than that in the neutral cases (0.4 for both neutral cases).
By comparing the diagnostic curves in the first linear region of the two stable cases with those of the neutral cases, it is clear that thermal stability remarkably shrinks the linear region as well as the intercept values, but does not alter the slope of the diagnostic curve.
3.3. Skewness and kurtosis
Skewness, 
$S_q={\overline {q^3}}/{\sigma _q^3}$, and kurtosis, 
$K_q={\overline {q^4}}/{\sigma _q^4}$, where 
$q$ represents the streamwise (wall-normal) velocity fluctuations, 
$u$ (
$w$), are associated with the intermittency of turbulent BLs, and their dependence on thermal stability is of practical interest in atmospheric turbulence (Maurizi Reference Maurizi2006).
In the neutral case, the skewness (kurtosis) of the streamwise velocity, i.e. 
$S_u$ (
$K_u$), varies insignificantly in response to the roughness change when comparing its values within the IBL with that at the same height but prior to the roughness change, as shown in figure 5(a) (figure 6a). In contrast, the skewness and kurtosis of the wall-normal component, namely 
$S_w$ and 
$K_w$, are modified in a remarkable way. Underneath the IBL, noticeable peaks arise in 
$S_w$ (e.g. see 
$z/\delta _0\approx 0.15$ for 
$x=0.72\ {\rm m}$ in figure 5b). Peaks appear at the same heights in 
$K_w$ with their values exceeding 4, implying a non-Gaussian probability density function of 
$w$. For stable cases in figures 5(d) to 5(f), the peaks of 
$S_w$ and 
$K_w$ become more noticeable within the IBL, and their magnitude gets amplified by the thermal stability. Furthermore, 
$S_u$ is distinguished from that in the neutral case and has a negative wide peak within the IBL, for instance, 
$z/\delta _0=0.42$ at 
$x=5.88\ {\rm m}$ in figure 5(c); peaks arise at a similar height in 
$K_u(z)$ (figure 6).
Figure 5. Vertical profiles of skewness of streamwise velocity 
$S_u$ (a,c,e) and wall-normal velocity 
$S_w$ (b,d,f) at various 
$x$ locations. Data of (a,b) are for case 1 (
$Ri_b=0, Re_{\delta }=4.5\times 10^4$), data of (c,d) are for case 2 (
$Ri_b=0.13, Re_{\delta }=4.5\times 10^4$) and data of (e,f) are for case 4 (
$Ri_b=0.27, Re_{\delta }=2.9\times 10^4$). Vertical dashed lines represent upper edges of IBLs determined from 
$\sigma _u$ for the data with the same colour shading. Colours/symbols are defined in figure 1.
Figure 6. Vertical profiles of kurtosis of streamwise velocity 
$K_u$ (a,c,e), wall-normal velocity 
$K_w$ (b,d,f) at various 
$x$ locations. Data of (a,b) are for case 1 (
$Ri_b=0, Re_{\delta }=4.5\times 10^4$). Data of (c,d) are for case 2 (
$Ri_b=0.13, Re_{\delta }=4.5\times 10^4$). Data of (e,f) are for case 4 (
$Ri_b=0.27, Re_{\delta }=2.9\times 10^4$). Vertical dashed lines from left to right present the upper edges of IBLs determined from 
$\sigma ^2_u$. Colours/symbols are defined in figure 1.
To estimate the quantitative thermal impact on high-order moments, the locations of 
$K_w$ peaks, 
$\delta _k$, their width, 
$l_k$, and magnitude, 
$K_w^{max}$, were studied. As shown in figure 7(a), 
$\delta _k$ is around 
$0.1\delta _0$ lower than the upper edge of the IBL and grows with the IBL depth in all cases. The regions of 
$K_w>K_w^c$ with a thickness of 
$l_k$ are illustrated by vertical bars in figure 7(a). Here, 
$K_w^c=3.3$ denotes the mean value of 
$K_w$ upstream of the roughness change. These regions are named non-adjusted regions hereafter as the notable peak in the skewness (kurtosis) of 
$w$ arises. Along the fetch, these regions expand in height from a thickness of 
$0.1\delta _0$ to 
$0.2 \sim 0.3\delta _0$. It is worth noting that their upper limits exceed the top edges of the IBL, which are identified from the merging point of 
$\sigma _u^2$ with its counterpart upstream to the roughness change in Ding et al. (Reference Ding, Placidi, Carpentieri and Robins2023).
Figure 7. (a) Height of the IBL, 
$\delta _i$ (closed symbols), and the locations of the maximum 
$K_w$, 
$\delta _k$ (open symbols), for several cases. The bar on each data shows the height range with 
$K_w$ larger than 
$K_w^c$. The curves for 
$Ri_b=0.13$ (0.27) are shifted downwards by 0.3 (0.6). (b) Variation of 
$\Delta K_w$ with 
$\delta _i/\delta _0$.
Figure 7(b) shows the interesting tendency of 
$\Delta K_w$ with 
$\delta _i/\delta$ increasing, where 
$\Delta K_w=K_w^{max}-K_c$. For the neutral case, 
$\Delta K_w$ shows non-monotonic variation with height such that it decreases with height when 
$\delta _i/\delta <0.4$, then reaches a minimum around 
$\delta _i/\delta \approx 0.4$ and finally increases with height. The decreasing tendency of 
$\Delta K_w$ with height is also observed in stable cases, where 
$\Delta K_w$ has larger magnitudes for the more stable case. The increasing tendency in the outer region is absent in stable cases due to the shallow nature of the IBL.
In studying the effect of thermal stability on the high-order moments, we note that the larger degree of thermal stability of case 4 (
$Ri_b=0.13, Re_{\delta }=2.9\times 10^5$) compared with case 2 (
$Ri_b=0.13, Re_{\delta }=4.5\times 10^5$) was achieved by reducing 
$Re_{\delta }$, which might affect the results. To examine this possibility the results (including 
$\delta _k$ and 
$\Delta K_w$) for case 1 (neutral, 
$Re_{\delta }=4.5\times 10^5$) and case 3 (neutral, 
$Re_{\delta }=2.9\times 10^5$) were compared, and they were found to vary little despite the differences in 
$Re_{\delta }$. Therefore, we believe that the differences in skewness and kurtosis between the two stable cases can be attributed to the variation in thermal stability and are independent of the Reynolds number.
3.4. Quadrant analysis
3.4.1. The 
$Q$ events upstream of the roughness change
Prior to the investigation into the response of the dynamics of BLs to a step change in roughness, we study the effect of thermal stability and of Reynolds number on the dynamical events characterising the surface upstream to the roughness change using quadrant analysis, as in Lu & Willmarth (Reference Lu and Willmarth1973). The definition of quadrant events is briefly introduced here
\begin{equation} \overline{uw}_i=\frac{1}{N}\varSigma_{j=1}^N h_i uw, \end{equation}
where 
$N$ is the sample size. The parameter 
$h_i=1$ if the pair 
$(u,w)$ locates in the 
$i$th quadrant and satisfies the condition 
$uw>H\sigma _{u}\sigma _{w}$, otherwise, 
$h_i=0$. The threshold 
$H$ denotes the size of the hyperbolic hole for the conditional statistics. The second quadrant 
$\overline {uw}_2$ with 
$(u<0,~w>0)$ is associated with low-momentum lifting (ejection 
$Q_2$) events, and the fourth quadrant 
$\overline {uw}_4$ with 
$(u>0,~w<0)$ is related to high-momentum injection (sweep 
$Q_4$) events. The 
$Q_4$ events are distributed most probably on the outer side of hairpin vortex legs while 
$Q_2$ events arise between two legs of hairpin vortices (Adrian, Meinhart & Tomkins Reference Adrian, Meinhart and Tomkins2000). These two events are major contributors to the Reynolds shear stress and the turbulent kinetic energy production. The 
$Q_1$ and 
$Q_3$ events are associated with (
$u>0,w>0$) and (
$u<0,w<0$), respectively.
Figure 8(a) shows 
$\overline {uw}_i^+(H=0)$ at 
$x=-0.72\ {\rm m}$, where the dominant contributions to the shear stress come from 
$Q_2$ and 
$Q_4$ (Lu & Willmarth Reference Lu and Willmarth1973). Comparing the data with the same 
$Re_\delta$ number, e.g. case 1 vs case 2 or case 3 vs case 4, it is apparent that the thermal stability alters the balance between 
$Q_2$ and 
$Q_4$ events. The 
$Q_2$ events are reduced to a greater extent than 
$Q_4$, while 
$Q_1$ and 
$Q_3$ are rarely altered by the thermal stability. These results are consistent with the findings in Williams et al. (Reference Williams, Hohman, Van Buren, Bou-Zeid and Smits2017), where they argued that the greater damping of 
$Q_2$ when compared with that of 
$Q_4$ is due to the higher temperature gradient characterising the near-wall region where lifting events (ejections) are generated. Therefore, the ratio of 
$\overline {uw}_2$ to 
$\overline {uw}_4$ should be reduced by the thermal stability, as shown in figure 8(b).
Figure 8. (a) Inner-scaled shear stress of four events as functions of height. Closed symbols are for 
$Re_{\delta }=4.5\times 10^4$, open symbols are for 
$Re_{\delta }=3\times 10^4$. (b) Ratio 
$\overline {uw}_2/\overline {uw}_4$ as a function of height. In (a), case 1 (black solid symbols); case 2 (blue solid symbols); case 3 (black empty symbols); case 4 (light blue empty symbols).
As for 
$Re_\delta$ effects, magnitudes of 
$\overline {uw}_i^+$ are reduced by more than half in case 3 when compared with those at the same height in case 1, as shown in figure 8(a). Considering that the BL with a higher 
$Re_{\delta }$ is able to yield more turbulent hairpin vortices, all events in case 1 (
$Re_\delta =4.5\times 10^4$) are thus expected to take place more frequently than those in case 3 (
$Re_\delta =2.9\times 10^4$) (Priyadarshana & Klewicki Reference Priyadarshana and Klewicki2004). Regardless of 
$Re_\delta$ effects, the ratio 
$\overline {uw}_2/\overline {uw}_4$ is invariant, as shown in figure 8(b), suggesting that the dynamics of hairpin vortices, hence the balance between 
$Q_2$ and 
$Q_4$ events, is independent of 
$Re_\delta$.
3.4.2. Response of $Q$ events to the roughness change
To investigate the response of each type of quadrant event to the surface roughness transition, we examine the relative contribution from each quadrant event to the total shear stress by looking at the ratio 
$r_{i,H}=|\overline {uw}_{i,H}/\overline {uw}|$. The results for case 4 with 
$H=0$ are presented in figure 9. Downstream of the roughness change, 
$r_{1,0}$ and 
$r_{3,0}$ within the IBL are reduced when compared with those at the same height ahead of the step, as shown by the darker shading region under the solid curve in figures 9(a) and 9(c). The contour plots of 
$r_{2,0}$ and 
$r_{4,0}$ within the IBL have complicated structures. For 
$r_{2,0}$, a layer of amplified magnitude and thickness 
$0.1\delta _0$ arises beneath the upper edge of the IBL, which collapses well with the location of the non-adjusted region (as discussed in § 3.3). The amplification of 
$r_{2,0}$ in such a layer suggests more frequent (or more intense) ejection events. In contrast, the relative contribution 
$r_{4,0}$ from 
$Q_4$ events in this layer is reduced, suggesting that sweep events occur less frequently or subside. Below the non-adjusted region, 
$r_{4,0}$ (
$r_{2,0}$) increases (decreases) in magnitude but still remains smaller (larger) than the counterpart upstream of the roughness change.
Figure 9. Contour plots of contribution from quadrant events to the total shear stress for case 4 
$(Ri_b=0.27)$ with 
$H=0$. The colour shading in (a–d) represents 
$r_i (i=1,2,3,4)$ respectively. The solid curves represent the upper edges of the IBL determined from 
$\sigma ^2_u$.
Given that quadrant events are dependent on the threshold value of 
$H$, we conducted the analysis by varying 
$H$ from 0 to 4. Figure 10 shows the comparison between the vertical profiles of 
$r_{i,H}$ at 
$x=3.72\ {\rm m}$ and their counterparts upstream of the roughness change 
$(x=-0.76\ {\rm m})$. Increasing 
$H$, the contributions from all events diminish. The difference between the two profiles (red and black lines) of 
$r_{i,H}$ diminishes with increasing 
$H$ for 
$Q_1$, 
$Q_3$ and 
$Q_4$ events, however, the difference of 
$r_{2,H}$ around 
$z=\delta _k=0.38\delta _0$ is maintained and even becomes greater. This result suggests that the amplification of 
$r_{2,H}$ in the non-adjusted region is ascribed to the appearance of extremely strong ejection events. This is supported by the local minimum (maximum) of 
$S_u$ (
$S_w$), which is observed in the same region. In neutrally stratified BLs, the modification of 
$Q_2$ and 
$Q_4$ events caused by the roughness change is minor or hardly discernible (at some 
$x$ locations for 
$H<2$). Increasing 
$H$, the enhancement in the contribution from 
$Q_2$ events becomes more noticeable.
Figure 10. Comparison of 
$\overline {uw}_{i,H}/\overline {uw}$ at 
$x=3.72\ {\rm m}$ (red curves) with that at 
$x=-0.76\ {\rm m}$ (black curves) with several 
$H$ values denoted in (a). The vertical dashed line denotes the upper edge of the IBL. Data for case 4.
3.4.3. Strong ejection events
A layer embedded underneath the IBL with strong ejection events is found to be universally present, both in neutral and stable cases. Figure 11 demonstrates the significant contrast of 
$Q_2$ events before and after the step change in roughness for the four cases, as evaluated by 
$\Delta r_{2,4}=r_{2,4}(x,z)-r_{2,4}(x=-0.76\ {\rm m}, z)$. In the case of 
$Ri_b=0.27$, as shown in figure 11(d), a noticeable compact layer with 
$\Delta r_{2,4}>0.15$ is located just below the upper edge of the IBL. The behaviour of this layer, including the variation of its width and strength of 
$\Delta r_{2,4}$ with height, shows consistency with that of 
$\Delta K_w$ in figure 7(b).
Figure 11. Contour plots of 
$\Delta r_{2,4}$ for (a) case 1, (b) case 3, (c) case 2 and (d) case 4. The solid curves denote the upper edges of the IBLs.
To explore the properties of these strong 
$Q_2$ events, we analyse the time series of 
$uw(t)$. In figure 12(b), the time trace of 
$-uw(t)$ at 
$z=\delta _k, x=3.72\ {\rm m}$ has several extremely sharp troughs (strong 
$Q_2$ events) when compared with the time trace of 
$uw(t)$ at the same height in the incoming flow (figure 12a). To investigate the time scales of these events, the complete structure of strong 
$Q_2$ events is extracted based on the following two-step procedure:
(i) The data with
$uw>4\sigma _u\sigma _w$ are identified and separated into isolated events (marked as red dots in figure 12a,b);(ii) for each event, its starting (
$t_s$) and ending ($t_e$) time stamps are determined as the time points when the event intersects with a threshold $\beta$ (see the green dashed line in figure 12c).
The value 
$H=4$ is chosen to ensure that the properties of strong 
$Q_2$ events are not statistically weakened by weaker events and 
$\beta =0.2$ is chosen to reduce the influence of the background noise in extracting the complete signal structure of 
$Q_2$ events.
Figure 12. Time series of 
$uw/\sigma _u\sigma _w$ at (a) 
$z=\delta _k$ at 
$x=-0.76\ {\rm m}$, and (b) 
$x=3.72\ {\rm m}$. (c) Partial zoomed-in view of (b). The red dots represent the data point identified within strong 
$Q_2$ events, 
$\tau$ denotes the time duration of a strong 
$Q_2$ event, 
$T$ denotes the time separation of two successive strong events and the horizontal dashed line has a value of 
$\beta =-0.2$. Data for case 4.
The time duration, 
$\tau$, of a single strong 
$Q_2$ event is defined as 
$\tau =t_e-t_s$ in figure 12(c), and 
$T$ represents the time separation between two successive maxima of 
$-uw/(\sigma _u\sigma _w)$. Figure 13 shows the event-averaged time duration 
$\langle \tau \rangle$ and the bursting period 
$\langle T\rangle$ of strong 
$Q_2$ events at several 
$x$ locations for 
$Ri_b=0.27$. In figure 13(a), it is noticeable that 
$\langle T\rangle$ is much smaller than that ahead of the step change with its minimum located at 
$z=\delta _k$ (the location of the maximum 
$K_w$, marked by arrows). Meanwhile, 
$\langle \tau \rangle$ reaches its minimum around 
$0.08\ \textrm {s}$ at 
$z\approx \delta _k$. Results suggest that strong 
$Q_2$ events are curtailed and take place more frequently in the non-adjusted region.
Figure 13. (a) Mean duration time 
$\langle \tau \rangle$, and (b) bursting period 
$\langle T \rangle$ of identified strong 
$Q_2$ events as a function of height at various streamwise locations. The arrows from left to right illustrate the locations of 
$\delta _k/\delta _0$ at 
$x=0.72$, 3.00 and 5.88 m, respectively. Colours/symbols are defined in figure 1.
The variation of the time scales associated with strong 
$Q_2$ events (at 
$z=\delta _k$) with fetch and the impact of thermal stability are summarised in figure 14. For all cases, 
$\langle \tau \rangle$ and 
$\langle T\rangle$ vary insignificantly with 
$x$ or with the thermal stability, having mean values of around 0.1 s and 2.5 s, respectively. To highlight the link between the time scale of strong 
$Q_2$ events to flow scales, the mixed time scale 
$\tau _m$ is estimated. In atmospheric BLs, 
$\tau _m$ is typically associated with the ‘mesolayer’, the length scale of which is the geometric mean of inner and outer scales and can be estimated by 
$\tau _m=[(\nu /u_*^2)(\delta /U_{ref})]^{1/2}$ (Alfredsson & Johansson Reference Alfredsson and Johansson1984). For the approach flow, 
$\tau _m\approx 0.05s$, which is of the same order as 
$\langle \tau \rangle$.
Figure 14. (a) Mean time duration 
$\langle \tau \rangle$, and (b) bursting period 
$\langle T \rangle$ of strong 
$Q_2$ events at 
$z=\delta _k$ as functions of 
$x$.
3.5. Thermal properties
3.5.1. Heat flux
Figure 15 shows the response of the vertical heat flux to the roughness change. The heat flux in the surface layer peaks just after the roughness change, before adjusting back down to the recovered state (
$x=5.88$ m). The heat flux recovery rate depends on the case under examination. The vertical profiles of heat flux downstream of the roughness change merge into those upstream at 
$\delta _\theta /\delta _0$, where 
$\delta _\theta$ is defined as thermal internal BL thickness. As shown in the inset plots in figure 15, the curve of 
$\delta _\theta$ (solid symbols) follows very closely that of 
$\delta _i$ (empty symbols) with a slight upward shift. This discrepancy is sensitive to the methodology used to evaluate both 
$\delta _\theta$ and 
$\delta _i$, but their growths are comparable.
Figure 15. Vertical profiles of heat flux at several 
$x$ for (a) case 2 and (b) case 4. Colours/symbols are defined in figure 1. Insets: the thickness of thermal IBLs 
$\delta _\theta$ (solid symbols) compared with 
$\delta _i$ (empty symbols).
To examine the influence of different events on the heat flux, we conduct quadrant analysis including the temperature, which divides the sample space (
$u,w,\theta$) into 8 quadrants, known as octant analysis after Suzuki, Suzuki & Sato (Reference Suzuki, Suzuki and Sato1988). Relevant events to the following discussion are: 
$i=2, u<0,w>0,\theta >0$ (warm ejections), 
$i=4, u>0,w<0,\theta >0$ (warm sweeps), 
$i=6, u<0,w>0,\theta <0$ (cold ejections) and 
$i=8, u>0,w<0,\theta <0$ (cold sweeps). The conditionally averaged heat flux is defined as 
$\overline {w\theta }_i=(({1}/{N})\sum _{j=1}^N h_i w\theta )$, where 
$N$ denotes the sample size. The quadrant parameter 
$h_i=1$ if 
$(u,w,\theta )$ is located in the 
$i$th quadrant and satisfies 
$uw>H\sigma _u\sigma _w$, and the contribution of each event is defined as 
$\overline {w\theta }_i/\bar {w\theta }$. Figure 16 presents the contributions from all events/interactions with their summation being 100 %. Taking 
$H=0$ for the incoming flow of case 4, cold ejections contribute approximately 70 % of the heat flux, which is slightly higher than the second-dominant warm sweeps. This result is consistent with the observation of the dominance of gradient transport in other studies (Wallace Reference Wallaces2016). The difference between these two types of events varies little with the threshold 
$H$. Just beneath the upper edge of the IBL at 
$x=3.72$ m, as shown in figure 16(d), the contribution from cold ejections increases up to 90 % and that from warm sweeps reduces to 40 %, while the contributions from other events remain nearly invariant for 
$H=0$. The difference between cold ejections and warm sweeps is maintained by increasing 
$H$. When 
$H=4$ (figure 16f), the contribution from warm sweeps and other events vanishes, while that from cold ejections remains substantial. Similar to the discussion in § 3.4 around the shear stress, ejections and sweeps are also unbalanced from the point of view of the heat flux and strong ejections are associated with the cold fluid parcels.
Figure 16. Vertical profiles of conditionally averaged heat flux based on octant analysis for case 4. (a–c) for 
$x=-0.76$ m and (d–f) for 
$x=3.72$ m. Panels show (a,d) 
$H=0$, (b,e) 
$H=2$, (c,f) 
$H=4$. The notation of the symbols: 
$i=1,u>0,w>0,\theta >0$ (solid square); 
$i=2,u<0,w>0,\theta >0$ (solid circle, warm ejections); 
$i=3,u<0,w<0,\theta >0$ (solid diamond); 
$i=4,u>0,w<0,\theta >0$ (solid triangle, warm sweeps); 
$i=5,u>0,w>0,\theta <0$ (empty square); 
$i=6,u<0,w>0,\theta <0$ (empty circle, cold ejections); 
$i=7,u<0,w<0,\theta <0$ (empty diamond); 
$i=8,u>0,w<0,\theta <0$ (empty triangle, cold sweeps). The vertical dashed lines in (d–f) represent the upper edge of the IBLs.
3.5.2. Conditional average of high-order moments
To shed light on the impact of thermal stability on turbulence statistics in the non-adjusted layer, we performed a conditional analysis of high-order moments conditioned with thermal properties. The conditionally averaged skewness and kurtosis of the vertical component are defined as 
$S_{w,i}=(({1}/{N})\sum _{j=1}^N h_i w^3)/\sigma _w^3$ and 
$K_{w,i}=(({1}/{N})\sum _{j=1}^N h_i w^4)/\sigma _w^4$. Their complements are also examined given the definition 
$\tilde {S}_{w,i}=(({1}/{N})\sum _{j=1}^N \tilde {h}_i w^3)/\sigma ^3_{w}$ and 
$\tilde {K}_{w,i}=(({1}/{N})\sum _{j=1}^N \tilde {h}_i w^4)/\sigma _w^4$, where 
$\tilde {h}_i=1$ for 
$(u,w)$ in the 
$i$th quadrant satisfying the condition of 
$uw \le H\sigma _u\sigma _w$. These 
$\tilde {}$ quantities represent contributions from events with small magnitudes of 
$uw$. e.g. the relation between 
$S_{w,i}$ and its complement is prescribed by 
$S_w=\sum _{i=1}^8 (S_{w,i}+ \tilde {S}_{w,i})$.
Figure 17(a) shows that 
$S_{w,6}$ (positive) and 
$S_{w,4}$ (negative) dominate over other components, suggesting that the major contributions to the high-order moments come from cold ejections and warm sweeps. Given that the summation of all components (solid curves) returns the vertical profile of 
$S_w$ observed in figure 5(f), the positive peak of 
$S_w$ is a result of the magnitude of 
$S_{w,6}$ exceeding 
$S_{w,4}$.
Figure 17. Conditionally averaged (a,c) skewness of vertical velocity and (b,d) kurtosis of vertical velocity for case 4 at 
$x=3.72$ m. Panels show (a,b) 
$H=0$; (c,d) 
$H=4$. The solid lines in (a,c) denote the summation of all quadrant components with 
$i\in [1,8]$. The vertical dashed line denotes the upper edge of IBL. Symbols are defined in figure 16. Insets of (c,d) show complementary components of skewness and kurtosis, i.e. 
$\tilde {S}_{w,i}$ and 
$\tilde {K}_{w,i}$.
Here, we also examine the result of 
$H=4$. Figure 17(c) shows that 
$S_{w,4}$ becomes zero while 
$S_{w,6}$ remains substantial. The inset plot shows that the summation of 
$\tilde {S}_{w,i}$ of all components becomes zero, suggesting a null contribution to 
$S_w$ from events with 
$H<4$. These results provide evidence of a firm link between strong cold ejections and the peak in 
$S_{w}$, as discussed in § 3.3. The dominant contribution from strong cold ejections on 
$K_w$ is also evident in figure 17(d) as the peak of 
$K_{w,6}$ is located at 
$\delta _k$. Based on the above analysis, we can conclude that the enhanced shear stress as well as the local maxima of skewness and kurtosis of velocity underneath the edge of the IBL can be ascribed to the strong ejections.
3.5.3. Local gradient Richardson number
In the final part of this section, we want to provide some insight into the strong ejection events discussed so far. Regarding neutral BLs, Antonia & Luxton (Reference Antonia and Luxton1971a) found that the large skewness of wall-normal velocity observed in this region can be ascribed to peaks in the turbulent diffusion term. Our analysis on conditionally averaged high-order moments provides a firm link between these strong ejections and high kurtosis/skewness, which suggests that strong ejections correspond with energy gained by diffusion.
In stable flows, the enhanced contribution of ejections can appear counter-intuitive as the thermal stratification is expected to dampen the fluctuations and weaken lifting events (due to the potential energy gradient). These effects can be quantified by the relative strength between the thermal stratification and the wall-normal shear, which is evaluated quantitatively by the gradient Richardson number 
$Ri_g=(g/\varTheta )(\partial \varTheta /\partial z)/(\partial U/\partial z)^2$. Miles (Reference Miles1961) and Howard (Reference Howard1961) proposed a necessary condition for instability in an inviscid, incompressible stratified flow, i.e. 
$Ri_g< Ri_c$, where the critical Richardson number 
$Ri_c=0.25$. Turbulent bursting corresponding to 
$Ri_g<0.25$ was observed in stable atmospheric BLs due to the strong shear imposed by lower-level jets (Ohya, Nakamura & Uchida Reference Ohya, Nakamura and Uchida2008). Given that the mean streamwise velocity in the lower region of our studied cases is vastly reduced due to the rougher downstream surface, it is interesting to examine the local shear induced by the roughness change in stable layers. Figure 18 presents contours of the local gradient Richardson number; here, regions of 
$Ri_g<0.25$ are highlighted. Upstream of the roughness change the region of 
$Ri_g<0.25$ is restricted to the surface layer, while after the roughness change, this region expands in height along the fetch, closely following the growth of the IBL (this is particularly obvious in figure 18b). This suggests that the deceleration of mean streamwise velocity within the IBL, induced by the roughness change, can trigger a strong shear, resulting in 
$Ri_g$ dropping below the critical value. This can lead to the modification of the turbulence dynamics within the IBL which, we infer, might be the reason for the enhanced contribution of ejections observed in the two stable cases. The critical Richardson number (
$Ri_c$) can vary between 0.2 and 1 depending on the scenarios studied and its determination in turbulent environments is particularly challenging (Galperin, Sukoriansky & Anderson Reference Galperin, Sukoriansky and Anderson2007), however, this is outside the scope of this manuscript. Similarly, a further investigation into the spatial organisation of flow structures, such as hairpin vortices, would be needed for a more thorough interpretation of the intense ejections found just below the IBL depth.
Figure 18. Contour plots of the local gradient Richardson number 
$Ri_g$ for case 2 (a) and case 4 (b). The solid curves represent the upper edges of IBLs.
4. Conclusions
In this paper, we studied the statistical properties of neutrally and stably stratified boundary layers in response to a step change in surface roughness from a rough to a rougher surface. By comparing the results of four cases with various 
$Re_\delta$ and 
$Ri_b$ values, the effects of thermal stability on the statistical properties of developing flows are studied. Our primary findings are listed below.
(i) The skin-friction relationship,
$C_f(\delta ^*)$, for the stably stratified BLs in equilibrium has been derived for the first time, extending the approach of Castro (Reference Castro2007) to stable flows. The new relation asymptotically approaches that in the neutrally stratified BL of Castro (Reference Castro2007) when the Obukhov length approaches infinity. The overshoot behaviour of $C_f(x)$ during the transition to a rougher downstream surface is slightly enhanced by the thermal stability.(ii) The normalised length scale
$\varDelta =\delta ^*u_*/(U_{\infty })$ is found to be able to collapse the velocity defect profiles for neutrally stratified IBLs onto those upstream of the roughness change. Influenced by the Obukhov length scale, $\varDelta$ becomes less suitable for collapsing the velocity defect profiles in stably stratified BLs. Further investigation into a more appropriate scaling is needed.(iii) The diagnostic plots are significantly altered by thermal stability. Upstream of the roughness change, the relation between
$\sigma _u/U$ and $U/U_\infty$ breaks into two linear regions. The first linear relation has a slope of $-0.255$ (the same as that in the neutral cases), while the magnitude of the slope in the second linear region becomes smaller. Downstream of the roughness change, the diagnostic curves in the lower region of the outer layer initially overshoot and then approach the fully developed state, in contrast to the more gradual approach in neutral cases.(iv) A region with large magnitudes of skewness and kurtosis of streamwise velocity and wall-normal velocity is identified just beneath the interface of the IBL with the outer flow for stably stratified flows. Quadrant analysis reveals that this feature/region is characterised by an increased occurrence of strong ejections (
$Q_2$) and a reduced occurrence of sweeps ($Q_4$) when compared with the flow upstream of the roughness change. This non-adjusted region with strong ejection events is present in both neutrally and stably stratified IBLs and becomes more noticeable with increased thermal stability.(v) In stably stratified BLs, these strong ejections are associated with cold fluid parcels and make major contributions to the heat flux, which is around twice that of the warm sweeps. These strong cold ejections account for the peaks of skewness and kurtosis observed in the non-adjusted region.
In conclusion, this work demonstrated the significant impact of thermal stability on the statistical properties and dynamics within the IBL.
Acknowledgements
S.-S.D. wants to thank Dr X. Jiang for insightful discussion.
Funding
This work is supported by NERC under the agreement ASSURE: Across-Scale processeS in URban Environments (NE/W002825/1).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The wind tunnel data utilised in this study can be accessed through the following link: https://doi.org/10.17605/OSF.IO/KH58X.
Author contributions
S.-S.D. conducted data analysis and wrote the first draft of the manuscript. All authors have commented/contributed to each version of this manuscript.
Appendix A. Mean velocity profiles – surface scaling
Figure 19(a) shows the inner-normalised mean velocity profiles using the 
$u_*$ value determined from 
$\overline {uw}$. A logarithmic layer is still identifiable in the lower part of the developing neutrally stratified BL where the data collapse well onto the logarithmic curve. With 
$x$ increasing, 
$U^+$ becomes progressively smaller due to the rougher downstream surface. Figure 19(b) shows scaling of the mean profiles in stable conditions. The log-corrected mean streamwise velocity near the wall follows 
$\beta _m(z-z_0)/z_0$, evidence of the existence of the surface layer (consider how, for the data at 
$x=5.88$ m, the region 
$0.05<(z-z_0)/L_0 <0.15$ aligns with the solid line). We note that this inner layer has a height of around 
$0.2L_0$ (
$0.22\delta _0$) in the incoming stable flow (see solid grey symbols). After the roughness change, this region grows in height (to around 
$0.46\delta _0$ at 
$x = 5.88$ m), suggesting a thicker surface layer for a rougher wall.
Figure 19. (a) Inner-scaled profiles of mean streamwise velocity, which is scaled by the local friction velocity determined by the extrapolation of 
$\overline {uw}$. The solid line represents the logarithmic function of 
$U^+=({1}/{\kappa })\ln (z/z_0)$. (b) Inner-scaled mean streamwise velocity after correction of law of the wall as a function of 
$(z-z_0)/L_0$. The solid line represents a linear function with a slope of 8 and intercept on the 
$y$-axis at 0. (a) Is for the neutral case (
$Ri_b=0,Re_{\delta }=4.5\times 10^4$) and (b) is for the stable case (
$Ri_b=0.13, Re_{\delta }=4.5\times 10^4$). Colours/symbols are defined in figure 1.
Figure 20(a) shows the velocity scaled by the local 
$u_*$ calculated from Elliott's model. The lower part of the developing neutrally stratified boundary layer collapses well onto the logarithmic curve when 
$z$ is scaled by the local roughness length 
$z_0$. With an increase in 
$x$ the wake region progressively lifts due to the rougher downstream surface.
Figure 20. (a) Inner-scaled profiles of mean streamwise velocity, which is scaled by the local friction velocity calculated from Elliott's model. The solid line represents the logarithmic function of 
$U^+=({1}/{\kappa })\ln (z/z_0)$. (b) Inner-scaled mean streamwise velocity after correction of law of the wall as a function of 
$(z-z_0)/L_0$. The solid line represents a linear function with a slope of 8 and intercept on the 
$y$-axis at 0. Panels show (a) 
$Ri_b=0,Re_{\delta }=4.5\times 10^4$ and (b) 
$Ri_b=0.13, Re_{\delta }=4.5\times 10^4$. Colours/symbols are defined in figure 1.
Figure 21. Velocity defect 
$(U_\infty -U)^+$ as a function of (a,c,e,g) 
$z/\delta$ and (b,d,f,h) 
$z/\varDelta$ at several streamwise locations. The inset plots in (a,c,e,g) only show the data outside the IBL for 
$x>0$ and the complete profile at 
$x=-0.76\ {\rm m}$. The inset plots in (b,d,f,h) only show the data within the IBL for 
$x>0$ and the complete profile at 
$x=-0.76\ {\rm m}$. (a,b) Case 1; (c,d) case 3; (e,f) case 2; (g,h) case 4. The solid curves represent the theoretical expression given by (3.4) with the wake strength 
$\varPi =0.70$ for (a–d), 1.0 for (e,f) and 1.2 for (g,h). Local 
$u_*$ is determined from the linear extrapolation method. Colours/symbols are defined in figure 1.
Finally, the mean velocity in defect form using 
$u_*$ from the linear extrapolation methodology is presented in figure 21. These data support that the collapse of the velocity defect discussed in § 3.2.1 is insensitive to the methodology adopted to calculate the friction velocity.
Appendix B. Two-dimensionality of the mean flow
Given that most of the analysis presented in the manuscript is performed along the centre line of the wind tunnel and given the 3-D nature of the roughness, it is important to assess the homogeneity of the IBL along the spanwise direction. Vertical profiles at eight different spanwise locations, spanning one period of the roughness elements arrangement, were measured both in neutral and stable cases at 
$x = 0.72$ mm and 
$x = 5.88$ mm for this purpose. The former streamwise location is just in front of a roughness element while the latter is in the middle of two rows of roughness elements, as shown in figure 1. Figure 22 shows an example of contour plots of the mean velocity, mean temperature, as well as the covariances in the 
$y\unicode{x2013}z$ plane for case 4 at 
$x=0.72$ m. Above the surface layer (
$z>0.1$ m), the homogeneity in the spanwise direction is evident, while the mean flow below this height is affected by the roughness arrangement. We note that, at 
$x=0.72$ m, the top of the IBL reaches the top of the surface layer. The difference in 
$\delta _i$ at various 
$y$ locations is around 12 % for the neutral case (case 3) and 
$\approx$5 % for the stable case (case 4) at 
$x = 0.72$ m. This difference decreases to 
$\approx$8 % for case 3 and 
$\approx$4 % for case 4 at 
$x = 5.88$ m. Regarding the surface properties, the difference between the friction velocity calculated along the wind tunnel centre line (at 
$y=0$) and its spanwise mean value for both cases at two streamwise locations is below 6 %. Given the intrinsic accuracy of the methodology used herein, the surface properties measured along the centre line of the wind tunnel can be considered to be representative, more generally, of the surface as a whole.
Figure 22. Contour plots of mean streamwise velocity 
$U$ (a), mean temperature 
$\varTheta$ (b), the shear stress 
$\overline {uw}$ (c) and the heat flux 
$\overline {w\theta }$ (d) in the 
$y\unicode{x2013}z$ plane. The black line and symbol represent the top of the IBL. Example for case 4 at 
$x=0.72$ m.
Open access
