Neutrally- and stably-stratified boundary layers adjustments to a step change in surface roughness

Abstract

In this work, we study the development of the internal boundary layer (IBL) induced by a surface roughness discontinuity, where the downstream surface has a roughness length greater than that upstream. The work is carried out in the EnFlo meteorological wind tunnel, at the University of Surrey, in both thermally neutral and stable cases with varying degrees of stability. For the neutrally-stratified boundary layer, the IBL development in the log-law region shows good agreement with the diffusion model proposed by Panofsky and Dutton (Atmospheric turbulence, Wiley, New York, 1984) provided that a modified origin condition is introduced and its growth rate is dictated by a constant diffusion term. However, the model over-predicts the growth of the IBL in the outer layer, where the IBL depth grows slowly with fetch following a power function with exponent n being 0.61 (whereas the original model prescribes $n\approx 0.8$). For the stably-stratified boundary layers, n is found to further reduce as the bulk Richardson number, $\textrm{Ri}_\textrm{b}$, increases. The analysis of the top region of the IBL shows that the slow growth rate is due to a combination of the decay of the diffusion term and a significantly negative mean wall-normal velocity, which transports fluid elements towards the wall. Considering these two effects, a modified diffusion model is proposed which well captures the growth of the IBL for both neutrally and stably-stratified boundary layers.

Graphical abstract

1 Introduction

Heterogeneity of the surface topology, such as that caused by a change in surface roughness conditions along the wind direction, will induce a boundary layer adjustment region below an interface known as the internal boundary layer (IBL) (Garratt 1994). Such a layer arises ubiquitously in nature when the wind blows from forests to prairie, when the offshore atmospheric boundary layer encounters the coast line, etc. The adjustment of the atmospheric boundary layer has repercussions on many industrial applications, e.g. the local wind pattern ventilation rates, the wind loading on buildings, or the energy contained in the wind available to wind turbines. Therefore, understanding the relation between the IBL and the surface geometrical properties, as well as modelling the interface development, is advantageous in many applications where interpreting and forecasting the complicated interaction between the boundary layer and the surface morphology is of importance.

Much of the work, both experimental and numerical, has focused on understanding the development of the IBL due to an abrupt change in surface roughness. The first comprehensive observation of the development of the IBL in the laboratory was conducted by Antonia and Luxton (1971a), Antonia and Luxton (1971b), Antonia and Luxton (1972). In their first experiment (Antonia and Luxton 1971a), where the IBL was induced by a surface change from a smooth wall to a rough one with recessed roughness elements ($S \rightarrow R$), the height of the IBL $\delta _i$ was found to grow with fetch x in a power function fashion, that is, ${\delta _i} {\propto } x^n$ with the value of n being 0.79. Whereas in another experiment, similarly of $S\rightarrow R$(Antonia and Luxton 1971b), they observed $n=0.5$ when the roughness elements were protruding into the incoming flow. They suggested the difference in the growth rate of the IBL between upstanding and recessed roughnesses may be ascribed to the perturbation induced by the first row of roughness elements in the former. In their third experiment of $R \rightarrow S$ transition, with the upstream rough surface covered with the same recessed elements as those in their first experiment, $n=0.43$ was reported (Antonia and Luxton 1972).

Consistently, $n\approx 0.8$ was also reported for cases of $S \rightarrow R$ or $R_1 \rightarrow R_2$ (where $R_1$ stands for a rough upstream surface smoother than the rough downstream surface $R_2$) by Andreopoulos and Wood (1982), Efros and Krogstad (2011) and Gul and Ganapathisubramani (2022) in the laboratory where the downstream surfaces were constructed using sandpapers or recessed roughness elements. Much smaller values of the exponent ($0.21<n<{0.58}$) have, however, been observed in other $S \rightarrow R$ or $R_1 \rightarrow R_2$ scenarios when the roughness elements were raised. Cheng and Castro (2002) found $n=0.33$ for a step change from a 3D rough surface to a rougher one covered with upstanding ribs. Numerical work includes those considering the change from a smooth surface to a rough one with 2D upstanding rods by Lee (2015) and a sinusoidal roughness by Rouhi et al. (2019) using direct numerical simulation (DNS), those over arrays of cuboids by Tomas et al. (2016) and Sessa et al. (2018), and those from a smooth surface to a forest pattern by Lopes et al. (2015) using large-eddy simulation (LES). Sessa et al. (2020) also studied the influence of stable stratification on the IBL’s growth rate developed over cuboids elements using LES and found that the thickness of IBL was significantly reduced by stable thermal stratification. A large variation of the value of n was also reported in cases of $R_2 \rightarrow R_1$ or $R \rightarrow S$. Rouhi et al. (2019) reported $0.34<n<0.51$ when the flow went from a sinusoidal roughness to a smooth surface using DNS. Ismail et al. (2018) obtained $n=0.41$ for a canonical channel flow from 2D ribs to a smooth surface. Li et al. (2021) found $n\approx 0.8$ when the rough surface was created by sandpaper roughness in a wind tunnel. Bou-Zeid et al. (2004) using LES, found that the IBL development from $R_2 \rightarrow R_1$ and $R_1 \rightarrow R_2$ collapsed onto each other with $n=0.79$ for long enough roughness fetch.

To interpret the development of the IBL, several models have been proposed, which mainly stem from two ideas. The first one is originated by Elliott (1958) studying the lower atmospheric layer, where the mean streamwise velocity is constructed as a piecewise function, with the layer above the IBL having a upwind-logarithmic profile and the layer below the IBL assumed to be a logarithmic function with an artificial friction velocity. The model was further developed by Panofsky and Townsend (1964), Antonia and Luxton (1972), Ismail et al. (2018), Rouhi et al. (2019). Li et al. (2022) extended the model to boundary layers with finite thickness through introducing a wake function to the outer streamwise velocity profile, and prescribing the decay of the shear stress.

The second approach, based on the diffusion analogue of plumes, was introduced by Miyake (1965). This theory prescribes the growth rate of the IBL, $\frac{d\delta _i}{dt}$, to be determined by the vertical diffusion,

$$\begin{aligned} \frac{\text {d}\delta _i}{\text {d}t}=U(\delta _i)\frac{\text {d}\delta _i}{\text {d}x}=A\sigma _w(\delta _i), \end{aligned}$$

(1)

where $U(\delta _i)$ and $\sigma _w(\delta _i)$ represent, respectively, the mean streamwise velocity and the standard deviation of vertical velocity along the top of the IBL, and A is a constant. Considering a case of $R_1\rightarrow R_2$ where the friction velocity and roughness length of the upstream surface are $u_{*1}$, $z_{01}$ and those of downstream one are $u_{*2}$ and $z_{02}$, respectively, the model has an analytical solution when the IBL is within the logarithmic layer developing over the local surface condition. Assuming that the diffusion process is controlled by the downstream surface conditions and that $\sigma _w = c u_{*2}$, substituting $U(\delta _i)=u_{*2}/\kappa \ln (\delta _i/z_{02})$ into Eq. (1) with the initial condition $x=0, {\delta _i=\delta _i^{0}}$ yields

$$\begin{aligned} \delta _i \;{\ln }\;{\frac{\delta _i}{z_{02}}}-\delta _i-\delta _i ^0\; {\ln }\;{\frac{\delta _i^0}{z_{02}}}+\delta _i ^0=Ac\kappa x, \end{aligned}$$

(2)

where the parameter $\kappa$ represents the von Kármán coefficient (having a value of around 0.41), the constants c and A are usually assumed to be around 1.25 and 1, respectively. While we acknowledge the body of work on the exact value of the von Kármán constant and the discussion in the field regarding its constancy, this is not of relevance in the current work, hence it is not discussed here. The above solution yields ${\delta }_i {\sim } x^{n}$ with $n\approx 0.8$ in the limit of a long fetch. Panofsky and Dutton (1984) suggested that the initial condition for $x=0$ should be ${\delta _i^0=z_{02}}$. Jackson (1976) proposed that $\delta _i^0$ should be $h=h_2+2(d_2-d_1)$ with $h_2$ being the crest height of the rougher surface, and $d_2$($d_1$) the zero-displacement downstream (upstream) of the surface discontinuity. Savelyev and Taylor (2001) found the relation between A (in Eq. (1)) and the roughness length ratio, $M=\ln ({z_{02}/z_{01}})$, to be $A=1+0.1M$ through fitting of experimental data in both $S\rightarrow R$($R_1 \rightarrow R_2$) and $R\rightarrow S$ ($R_2 \rightarrow R_1$) scenarios. In their later work, the term 0.1M was proposed to be related to the streamline displacement, which involves the mean vertical velocity due to the flow acceleration or deceleration after the roughness change (Savelyev and Taylor 2005). Thus, they proposed a modified model,

$$\begin{aligned} \frac{U_1(\delta _i)d{\delta _i}}{dx}=A_1{{\sigma }_{w1}({\delta _i})}+A_2W({\delta }_i) \end{aligned}$$

(3)

where the subscript 1 refers to quantities in the upwind surface and $A_1$, $A_2$ are constant parameters. In this model the value of the mean vertical velocity W at the top of the IBL is assumed to be positive and it is modelled to be proportional to $\delta _i/x$ or $\text {d}\delta _i/\text {d}x$ using the continuity equation. They further extended the foregoing model to the system with a sharp change in surface heat flux—this is however outside the scope of the current work.

To resolve the inconsistency in the value of n, Rouhi et al. (2019) showed that this exponent is sensitive to the definition of the IBL (i.e. the procedure to determine its depth) and may decrease with decreasing the Reynolds number (Re). Except for these two possible alternatives, comparing previous work, the value of n seems to be highly dependent on the type of roughness elements used. One clue is in the two experiments conducted by Antonia and Luxton (1971a) and Antonia and Luxton (1971b) where the same procedure is employed to calculate the IBL depth and at a the same Re. The former work employs recessed roughness elements, meanwhile protruding elements are used in the latter. The exponent, n, was found to be 0.79 and 0.43, respectively. In addition, Sessa et al. (2018) reported a small value of n ($\approx 0.2$) in a system with protruding roughness elements independently from the definition/methodology to calculate the IBL depth.

Regardless of the many efforts devoted to modelling the development of IBLs, their seemingly slow growth (i.e. small values of n) when the roughness elements downstream to the surface discontinuity are protruding has yet to be understood, and shows a large disagreement with the prediction of the diffusion models which prescribe $n{\approx }0.8$ for long development fetch. (The latest modified diffusion model proposed by Savelyev and Taylor Savelyev and Taylor (2005) gives $n\ge 0.8$ for $S\rightarrow R$ ($R_1 \rightarrow R_2$) in a long fetch.)

However, upstanding roughness elements appear ubiquitous in natural systems, e.g. in the flow transition from rural to urban areas, or from farmlands to woodlands Understanding the physical cause of the slow growth of the IBL is fundamental to modelling local flowfield (and therefore climate) in these scenarios.

This work is, therefore, focused on studying the growth rate of the IBL induced by an abrupt change in surface roughness from $R_1\rightarrow R_2$ where the upstream surface is typical of atmospheric offshore site, and the rougher downstream surface conditions are typical of onshore urban roughness (with protruding roughness elements in both cases). Of interest are the cases of thermally neutral or stably-stratified boundary layers.

2 Experimental setup

2.1 Experimental facility

Fig. 1

a Schematic of the experimental setup in the wind tunnel. The left subplot shows the inlet temperature profile used in the stably-stratified boundary layers. The right subplot shows the surface roughness conditions. The sonic anemometer mounted at $z=1$ m, $y=1$ m and 5 ms from the inlet section provides a reference velocity of the incoming flow, $U_{ref}$. b The setup of the LDA probe, the thermistor and the cold wire. c Measurement locations and the coordinate system in the wind tunnel. The origin is set at the location of the roughness step change. The short (long) lines denote the locations of the upstream (downstream) roughness elements. The red dots represent the measurement locations for vertical profiles, some of which are labelled

Experiments were conducted in the EnFlo meteorological wind tunnel at the University of Surrey, which is a suck-down open-return type. The uniqueness of this facility is that it allows for the generation of boundary layers with non-neutral thermal stratification, via a combination of heating and cooling of the floor of the tunnel and the imposition of a desired temperature inlet profile. Stable flows are obtained by cooling the floor of the tunnel with imposed inlet temperature profiles, while convective cases can be achieved by heating the floor instead of cooling. 15 heating elements (drawing 405 kW) are uniformly spaced vertically and positioned at the inlet of the test section, able to generate the desired temperature profiles. The working section of the wind tunnel is 20 m long, 3.5 m wide and 1.5 m high. The wind speed ranges from 0.3 to 2.5 $\mathrm {m/s}$. A sonic anemometer mounted above the developing boundary layer at 5 m from the inlet of the tunnel provides a reference wind speed $U_\textrm{ref}$ of the incoming flow throughout the duration of the experiments. This is set to 1 or 1.5 $\mathrm {m/s}$ as further discussed in Sect. 2.4.

2.2 Rough surfaces and boundary layers

In order to study the internal boundary layer development from an abrupt change in surface conditions, two kinds of roughness elements are used as shown in Fig. 1a. The upstream rough surface stretches for 11 m from the inlet section, and comprises of roughness elements of 50 $\times$ 16 $\times$ 5 mm, standing on the cross area of 50 $\times$ 5 mm in 50% staggered pattern with spacing of 510 mm in the spanwise direction and 360 mm in the streamwise (both centre-to-centre). The downstream rough surface covers the next 7 m downstream from the step change in surface condition ($x=0$) and employs roughness elements of 80 $\times$ 20 $\times$ 2 mm, standing on the area of 80 $\times$ 2 mm to form of 50% staggered pattern with spacing of 240 mm (centre-to-centre) both in the spanwise and streamwise directions. Both surfaces are fairly open arrangements, as detailed by the values of the respective roughness lengths in Sect. 3. These rough walls, in isolation, were fully characterised, respectively, by Hancock and Hayden (2018) and Marucci et al. (2018). To be noted that of all roughness change problems of practical interest, here we have considered a relatively simple one, but the choice of upstanding elements is not an unusual condition in nature. Furthermore, the surface discontinuity considered here is a pure roughness step change, and it ignores variations in other flow properties which often characterise planetary boundary layers. These choices are in line with the aims of the manuscript, and their repercussions on the modelling are further discussed in Sect. 4.

At the inlet of the working section, 13 Irwin-spires with height of 600 mm are placed uniformly in the spanwise direction to generate an artificially-thickened boundary layer, whose depth is around 550 mm. The thermally stable stratified boundary layers were obtained by cooling at the floor the heated flow coming from the inlet of the wind tunnel, where 15 electrical heaters imposed a vertical temperature profile as shown, for example, in the subplot in Fig. 1b. The floor of the wind tunnel was kept un-cooled over the first 4 m from the inlet of the tunnel, while the rest was chilled to a prescribed surface temperature $\Theta _0$, which varied across cases to maintain the prescribed temperature difference $\Delta \Theta =16 K$. This is in accordance with recommendation from previous work by Hancock and Hayden (2018) that found this un-cooled length to be crucial for the thermal and mechanical properties of the boundary layer of the approach flow reaching equilibrium in a short fetch (around 9 m) on the offshore surface, varying smoothly with the wall-normal direction, and comparing well to meteorological data and surface scaling, hence the 11 m of offshore surface upstream to the roughness change is employed here.

2.3 Experimental methods

The coordinate system used herein is shown Fig. 1a and c; its origin is set in front of the first row of the downstream roughness along the central line of the wind tunnel floor. The streamwise velocity $u=U+u'$ and vertical velocity $w=W+w'$ along the streamwise direction x and the vertical direction z are measured by two-component laser-Doppler anemometer (LDA). U and W denotes the mean (time-averaged) values, $u'$ and $w'$ represent the velocity fluctuations. A Dantec 160 mm focal length Fibre Flow probe connected to the LDA system was mounted onto the traverse system of the tunnel, as shown in Fig. 1b. The LDA system was set to operate at a minimum sampling rate of 100 Hz, which can be fulfilled through controlling the seeding concentration automatically by an in-house LabView software. The beam separation on the probe measuring the streamwise velocity was calculated by an independent 1D standard LDA prior to the measurement.

A cold-wire anemometer and a thermistor were mounted on the same traverse system as the LDA probe. In order to measure the temperature as close to the measurement volume of the LDA, while ensuring that the temperature and velocity correlation was free of significant blockage effects, the location of the cold wire was carefully adjusted to be exactly at 4 mm downstream of the velocity measurement volume. An acceptable value of this distance should be comparable to the length calculated from the sampling frequency of the LDA and the minimal local mean wind speed. The thermistor is also placed 4 mm downstream of the measurement volume of the LDA but shifted 15 mm in y away from the cold-wire. The same setup of cold-wire and the thermistor has been amply applied in previous work (Hancock and Hayden 2018; Heist 1998). The cold wire is calibrated by a pre-calibrated thermistor. The mean potential temperature $\Theta$ was measured by the thermistor and the temperature fluctuations around this mean, $\theta$, were obtained through the cold wire. The sampling rate of the cold wire is set to be 1000 Hz, which is much larger than that of the LDA to allow for time interpolation and calculation of heat fluxes. The bias error in the mean velocity is within $\pm 1\%$, and $\pm 0.1\%$ for $\Theta$. For the second-order moments, e.g. any correlation between $u'$, $w'$ and $\theta$, the error is within $\pm 7\%$.

2.4 Case studies

Two neutral cases were studied in this work with $U_\mathrm{{ref}}=1.5$ m/s and $U_\mathrm{{ref}}=1$ m/s, yielding $\mathrm{{Re}}_{\delta }=4.5{\times }10^4$ and $3.0{\times }10^4$, respectively, with $\mathrm{{Re}}_{\delta }=\delta _0 U_\mathrm{{ref}}/{\nu }$ and $\delta _0$ indicating the boundary layer thickness of the approach flow. To study the impact of stable stratification on the IBL development, two cases with weak and moderate stable stratification were studied. The weak stable case is obtained by introducing the inlet temperature profile in Fig. 1a for $\mathrm{{Re}}_{\delta }=4.5{\times }10^4$, and bulk Richardson number $\mathrm{{Ri}}_\textrm{b}=\frac{g(\Theta _{\delta }-\Theta _0)\delta _0}{\Theta _0 U_\delta ^2}=0.13$, where g is the gravitational acceleration, $\Theta _{\delta }$ ($U_{\delta }$) the mean temperature (streamwise velocity) at the top of the boundary layer, and $\Theta _0$ the mean temperature on the floor. The moderate stable case was obtained by reducing the freestream wind speed to $U_\textrm{ref}=1$ m/s, while keeping the same inlet temperature profile as the weaker stable case, which yielded $\mathrm{{Ri}}_\textrm{b}=0.27$. The effects of this reduction in the Reynolds number are further discussed in Sect. 3. The details of experimental parameters and some of the quantities of interest for all four cases examined herein are listed in Table 1.

Figure 1c shows the 14 locations along the centre line of the tunnel where the measurements of vertical profiles were conducted. These included 40 height levels with a logarithmic-uniform spacing, starting at $z=50$ mm and ending at $z=800$ mm. To inform on the homogeneity of the IBL along the spanwise direction, vertical profiles at further 8 spanwise locations were measured both in neutral and stable cases at $x=720$ mm and $x=5880$ mm. The former streamwise location is just in front of a roughness element and the latter one is in the middle of two rows of roughness elements, as shown in Fig. 1c. These measurements are not included here as they confirmed the trends highlighted by the detailed measurements on the tunnel centre-line, which are discussed at length in the following.

3 Results and discussion

3.1 Changes in surface properties

Figure 2a shows profiles of $U/U_{\infty }$ as a function of $z/\delta$ for the neutral stratification at $\textrm{Re}_\delta =4.5\times 10^4$. In this study the thickness of the boundary layer $\delta$ for each vertical profile is defined as the height where U reaches 99$\%$ of the free stream velocity $U_{\infty }$, which is defined as the value of U at $z=800$ mm. Based on the observation that the vertical profile of mean streamwise velocity at $x=-0.38$ m collapses onto that at $x=-0.76$ m, the approach flows can be assumed to have reached equilibrium with the underlying surface and it can be considered fully developed. It is visible in Fig. 2a that in front of the first row of the downstream roughness elements (just before $x=0$ m), the vertical profile collapses with that at $x=-0.76$ m indicating no adjustment of the mean streamwise velocity within measurement region. With the increasing fetch in x, $U/U_{\infty }$ in the lower layer ($0.15<z/\delta <0.5$) starts to deviate from the curve at $x=-0.76$ m. For $x=3$ m, the curve of $U/U_{\infty }$ at small z merges with the profile at $x=5.88$ m implying that the flow in the lower layer has reached a new equilibrium which is now controlled by the downstream rough surface. To evaluate how the logarithmic region within the boundary layer is modified by the step change in surface properties, the friction velocity $u_*$ is estimated through linear extrapolation of the shear stress $\overline{u'w'}$ to the ground (Marucci et al. 2018), that is, $u_*=\sqrt{-\overline{u'w'}_0}$. Here, the overbar denotes time averaged quantities, and the subscript 0 indicates the value extrapolated to the ground ($z=0$). The effective roughness length is obtained by curve-fitting the data in the logarithmic region ($0.1<z/\delta <0.2$) using $U=u_*/\kappa \ln \frac{z-d}{z_0}$, where d is the displacement height, as used in Jackson (1981). The use of the formula is justified by the fully rough character of the wall surfaces at the studied Reynolds numbers. The solid lines in Fig. 2a show these fitting functions to the data.

Fig. 2

The vertical profiles of $U/U_{\infty }$ for a $\textrm{Re}_{\delta }=4.5\times 10^4$ and $\textrm{Ri}_\textrm{b}=0$ (Neutral—Case 1), and b $\textrm{Re}_{\delta }=4.5\times 10^4$ and $\textrm{Ri}_\textrm{b}=0.13$ (Stable—Case 2). The solid lines in a denote logarithmic fits at $x=-0.76$ m and $x=5.88$ m. The solid curves in b denote the fits into the data at $x=-0.76$ m and $x=5.88$ m using Eq. (5). See legend in a for symbols. Dashed lines denote the edge of the boundary layer

When the stable stratification comes into effect, the Monin–Obukhov similarity (Garratt 1984) dictates the scaling of the lower part of the boundary layer, which gives

$$\begin{aligned} 1+{\beta }_m{\zeta } =\frac{\kappa z}{u_*}\frac{\partial U}{\partial z}, \end{aligned}$$

(4)

where $\beta _m =8$ for the stably-stratified boundary layer (Hancock and Hayden 2018), and ${\zeta }=z/L_0$. The length scale $L_0$ is the surface Obukhov length defined as $L_0=u_*^2 {\Theta _0}/({g\kappa \theta _*})$, with the friction temperature ${\theta _*}$ determined from the relation ${\theta _*}=-\frac{\overline{w'\theta }_0}{u_*}$. The value of $\overline{w'\theta }_0$ is estimated through extrapolating $\overline{w'\theta }$ to the ground using a linear function. Integrating the above equation gives

$$\begin{aligned} U(z)=u_*/\kappa \left[ \ln {\frac{z-d}{z_0}}+{\beta _m}\frac{z-d-z_0}{L_0}\right] . \end{aligned}$$

(5)

During the fitting process, the displacement thickness, d, for both neutral and stable cases is set to zero, as suggested in Hancock and Hayden (2018), Marucci et al. (2018) for the same rough surfaces. This is justified by the sparseness of the roughness elements which induce isolated flow regime (Grimmond and Oke 1999).

Figure 2b shows an example of the vertical profiles of streamwise velocity in stable stratification (for $\textrm{Ri}_\textrm{b}=0.13$). The solid curves are fitting functions of the profiles at $x=-0.76$ m and $x=5.88$ m using Eq. (5) with ${\beta }_m=8$. Figure 3a shows the variation of friction velocity $u_*$ with fetch scaled by the upstream friction velocity $u_{*1}$, where $u_{*1}$ is defined as the $u_*$ at $x/{\delta _0}=-1.8$. After the step change in surface roughness, $u_*/u_{*1}$ starts to increase slowly, followed by a sharp increase as the flow adjust to the new boundary condition at the wall, and eventually it levels off around $x=6{\delta _0}$ for neutral flows. For stable flows, $u_*/u_{*1}$ overshoots after the initial increase and oscillates around its equilibrium value, where the flow has finally adjusted to the downstream roughness, hence scaling with $u_{*2}$. Especially for the case at $\textrm{Ri}_\textrm{b}=0.27$, the oscillation is found to be large in amplitude. Figure 3b shows the variation of the ratio of the roughness lengths along fetch in the form $\ln \;(z_0/z_{01})$, where $z_{01}=z_0(x=-1.8\delta _0)$. The picture is similar to that drawn for the friction velocity ratio $u_*/u_{*1}$, where oscillations in $\ln (z_0/z_{01})$ also take place for stable cases. Due to these oscillations over the long fetch, the equilibrium values of $u_*/u_{*1}$ ($z_{0}/z_{01}$) are here defined as the spatial-averaged values, that is, $\langle u_*(x)/u_{*1}\rangle$ ($\langle z_{0}(x)/z_{01}\rangle$), where $\langle .\rangle$ denotes streamwise spatial average of the quantity of interest within $6\delta _0<x<12\delta _0$. Regardless of the oscillations, it is clearly observable that both the values of $u_*/u_{*1}$ and that of $z_0/z_{01}$ rise sharply at the step change and then reach a plateau—this implies that a new equilibrium has been reached in the logarithmic layer formed over the downstream surface after the change in surface roughness. The friction velocity ${u_{*2}}$ of downstream surface is defined as the value of $\langle u_*(x)\rangle$, and similarly for the roughness length of the downstream surface $z_{02}$. For the stable cases, the surface Obukhov length of the downstream surface $L_{02}$ is also defined as $\langle L_{0}\rangle$, and for the upstream surface $L_{01}=L_0(x=-1.8\delta _0)$. The parameters of surface properties are listed in Table 1, where the value of $u_{*1}$ and $z_{01}$ for $\textrm{Re}_\delta =4.5 \times 10^4$ and $\textrm{Ri}_\textrm{b}=0.13$ are consistent with the results in Hancock and Hayden (2018) using the same roughness setup in the approaching flow. Typical fitting procedure uncertainty for rough-wall boundary layer parameters (e.g. $u_*$) is of the order of $\pm 10\%$ (Schultz and Flack 2005).

Comparing the surface properties of neutrally-stratified flows with those of stably-stratified flows for the same Reynolds number, the friction velocity and roughness length are both found to be reduced by the thermal stability. However, the ratios $u_{*2}/u_{*1}$ and $M=\ln (z_{02}/z_{01})$ are slightly enhanced by thermal stability.

Fig. 3

Scaled a friction velocity $u_*/u_{*1}$, and b roughness length $\ln {(z_{0}/z_{01})}$ as a function of fetch $x/\delta _0$. Here $u_{*1}=u_*({x/\delta _0=-1.8})$ and $z_{01}=z_0({x/\delta _0=-1.8})$. Symbols: ${\bigcirc }$ $\textrm{Re}=4.5{\times }10^4$ and $\textrm{Ri}_b=0$, ${\Box }$ $\textrm{Re}=3.0{\times }10^4$ and $\textrm{Ri}_\textrm{b}=0$, ${\triangledown }$ $\textrm{Re} = 4.5{\times }10^4$ $\textrm{Ri}_\textrm{b} = 0.13$ and ${\Diamond }$ $\textrm{Re}=2.9{\times }10^4$ and $\textrm{Ri}_\textrm{b}=0.27$

Table 1 Summary of the experimental parameters and quantities of interest

3.2 IBL growth

There are several approaches to determine the depth of the internal boundary layer. Based on dimensional analysis, Antonia and Luxton (1971a) plotted U as a function of $y^{0.5}$ and found a ’knee’ point where a sharp change in the slope of $U(y^{0.5})$ takes place. They determined the location of the ’knee’ point through fitting straight lines to the two regions near this point and defined the intersection as the thickness of the IBL. Efros and Krogstad (2011) defined the IBL thickness through fitting linear functions to ${\sigma }_u^2(z)$ in the internal layer and outer layer, where ${\sigma }_u^2(z)$ indicates the variance of streamwise velocity. The intersection of these two lines is defined as the thickness of the internal boundary layer. Sessa et al. (2018) extended the definition of IBL to the variance of vertical velocity ${\sigma }_w^2(z)$ and verified that the power exponent n obtained by the above three methods are close to each other. In the current work, we apply the latter two methods to determine the height of IBLs.

Figure 4 shows the development of ${\sigma }_u^2(z), {\sigma }_w^2(z)$ over the downstream roughness compared to those at $x=-0.76$ m. Along the fetch, ${\sigma }_u^2(z)$ and ${\sigma }_w^2(z)$ in the upper layer ($z>0.18\delta$ at $x=0.36$ m) fall onto their counterparts at $x=-0.76$ m. While in the lower layer ($z<0.18\delta$ at $x=0.36$ m) their values are enlarged due to the effect of the changing underlying roughness conditions and they start to deviate from the upstream profile at the vertical location $z={\delta }_i(x)$, which is determined by the intersection of two linear fitting functions.

Fig. 4

The vertical profiles of a ${\sigma }_u^2/u_{*1}^2$ and b ${\sigma }_w^2/u_{*1}^2$ at various x locations. The solid lines are linear fits to data to define the thickness of IBLs. The error bars represent the standard deviation of $\sigma _u^2$ and $\sigma _w^2$ in a and b, respectively. Neutral—Case 1 ($\textrm{Re}_{\delta }=4.5{\times }10^4$ and $\textrm{Ri}_\textrm{b}=0$)

Fig. 5

Contour plots of a U, b W, c $\sigma _u^2$, and d $\sigma _w^2$. The dots on the solid curves represent the upper limit of the IBL, which is obtained from the vertical profiles of $\sigma _u^2$. Neutral—Case 1 ($\textrm{Re}_{\delta }=4.5{\times }10^4$ and $\textrm{Ri}_\textrm{b}=0$.)

Fig. 6

Contour plots of a U, b W, c $\Theta$, d $\sigma _u^2$, e $\sigma _w^2$ and f $\overline{-w' \theta }$. The dots on the solid curves represent the upper limit of the IBL, which is obtained from the vertical profiles of $\sigma _u^2$. Stable—Case 2 ($\textrm{Re}_{\delta }=4.5{\times }10^4$ and $\textrm{Ri}_\textrm{b}=0.13$)

The development of IBLs across the step change in surface conditions is illustrated by black dots in Fig. 5 for $\textrm{Re}_\delta =4.5\times 10^4$ and $\textrm{Ri}_\textrm{b}=0$ and in Fig. 6 for $\textrm{Ri}_\textrm{b}=0.13$. For both neutral and stable cases, the mean streamwise velocity U is reduced near the wall due to the blockage imposed by the new surface condition, and it is accelerated around the top of the boundary layer due to its lower part adjusting to a rougher surface. Meanwhile, the disruption induced by the first row of the downstream roughness elements causes a local positive mean vertical velocity around $x=0$ m in proximity of the wall. The vertical extent of positive W reaches a height of 10$h_2$ (200 mm) as shown in Fig. 5b, and it is found independent of thermal instability when comparing it with Fig. 6b. This is an indication that the wall-normal velocity disturbance is more significantly linked to the geometry of the roughness than it is to the state of the boundary layer–this is perhaps not surprising. A second structure characterised by a positive W is also visible in the contour plot in Figs. 5 and 6, and appears to be induced by the second row of roughness elements at $x=0.24$ m. Its vertical extent is reduced to around 150 mm for both neutral and stable stratification. Such positive W structure in the contour plots disappears at $x=0.72$ m, which is located at the front of the forth row of the downstream roughness elements (see Fig. 1c). Further downstream, the value of W becomes negative. The positive and negative signs of W are related to the deceleration of U near the ground and the acceleration aloft as prescribed by continuity.

With respect to the variance of streamwise velocity and vertical velocity in Fig. 5c, d, their values on the downstream roughness surfaces are around 2 times larger than those over the approaching flow (i.e. upstream roughness). For the stable stratification, the negative heat flux, $\overline{w' \theta }$, becomes stronger in magnitude after the abrupt change in surface roughness, leading to a smaller temperature near the floor when compared to that in the upstream roughness.

To study the growth curves of the IBL, its depth, $\delta _i$, is plotted as a function of x for four cases in Fig. 7. Since the power-exponent of $\delta _i(x)$ is sensitive to the choice of the origin location, Jackson (1976) suggested correcting the original height of the IBL by the value of $h_2+2(d_2-d_1)$, where $h_2$ denotes the crest height of the downstream roughness elements, $d_2$($d_1$) are the zero-displacements of the downstream(upstream) surfaces, respectively. The growth curves of IBLs after the origin correction of two neutral cases show a good concurrence, implying the development of IBLs is independent of $\textrm{Re}_{\delta }$ in our experiment. Since the reduction in the Reynolds number has no significant effect on the IBL development, the moderately stable case ($\textrm{Ri}_\textrm{b}=0.27$) is obtained by reducing $\textrm{Re}_{\delta }$ while keeping the temperature difference $\Delta {\Theta }$ fixed to that of the $\textrm{Ri}_\textrm{b}=0.13$ case.

Fig. 7

The thickness of IBLs after the origin correction $(\delta _i-h_2)/h_2$ vs. $x/h_2$, where $h_2$ denotes the height of the downstream roughness elements. Symbols represent the same conditions as in Fig. 3. Filled symbols are determined from $\sigma _w^2$ and open ones from $\sigma _u^2$. The dashed line represents the prediction from Eq. (6). The solid lines are fits to data using power functions. The dashed-dotted line denotes the spatially-averaged upper edge of logarithmic layer for the approaching neutrally-stratified boundary layers

The growth curve within the logarithmic region ($\delta _i<0.2{\delta _0}$) is found to be well predicted by the diffusion model proposed by Panofsky and Dutton (1984) provided that the IBL grows from $z=h_2$ at $x=0$ m. This can be expressed by the following formula,

$$\begin{aligned} \delta _i \ln {\frac{\delta _i}{z_{02}}}-\delta _i-h_2\ln {\frac{h_2}{z_{02}}}+h_2=0.5x, \end{aligned}$$

(6)

which is shown by the dashed line in Fig. 7. However, it is obvious from the same figure, that when reaching the outer layer (see data above the dashed-dotted line), the measured IBL depth is shallower than the prediction of Eq. (6), where all data points depart and lay below the dashed line. The growth rate of the IBL, however, is still well represented by a power function for all cases as discussed in the following. Through fitting data of two neutral cases data using a power function, the power-exponent n is found to reduce to around 0.61 in the outer layer. The slow growth rate of the IBL leads to the value of $\delta _i$ at $x=5.88$ m that is $41\%$ lower than that predicted by Eq. (6). For stable stratification, the IBL become even shallower; for $\textrm{Ri}_\textrm{b}=0.13$ the value of $\delta _i$ is reduced by $30\%$ compared to that in the neutral cases and by $54\%$ for $\textrm{Ri}_\textrm{b}=0.27$ (at $x=5.88$m). The power-exponent n is reduced to 0.46 for $\textrm{Ri}_\textrm{b}=0.13$ and 0.39 for $\textrm{Ri}_\textrm{b}=0.27$.

Finally, a note on the influence of the adopted identification procedure for the IBL estimation and the overestimation of the model in Panofsky and Dutton (1984) when compared to the experimental data. Bou-Zeid et al. (2004) studied IBLs (from $R_1 \rightarrow R_2$ and $R_2 \rightarrow R_1$) and found that although the growth rate of the IBL identified by multiple methods varied to a large extent, all of them could still be well described by the model of Panofsky and Dutton (1984). This is in contrast with the findings of this work. Here, the difference between the values of $\delta _i$ determined from $\sigma _w^2$ and $\sigma _u^2$ is within $5\%$. Taking Case 1 for instance, the value n in the outer layer determined from $\sigma _u^2$ is 0.60 and becomes 0.62 from $\sigma _w^2$, yielding a difference around $3\%$. Based on these findings, we use the $\delta _i$ and n determined from $\sigma _u^2$ hereafter. In addition, since any rough surface is inherently three-dimensional, we performed checks on the spanwise homogeneity of IBL depth at $x=0.72$ m and $x=5.88$ m. The difference in $\delta _i$ at various y locations is around $12.4\%$ for the neutral case (Case 3) and $5.3\%$ for the stable case (Case 4) at $x=0.72$ m. This difference becomes $7.6\%$ for Case 3 and $3.6\%$ for Case 4 at $x=5.88$ m. Due to these small differences in the IBL thickness in the y direction, we mainly focus on the IBL development along the centre line of the wind tunnel. We also tested employing the mean streamwise velocity profiles to identify the IBL as in Antonia and Luxton (1971b), with negligible effects on the growth curves of $\delta _i$, which largely collapsed onto those derived from $\sigma _u^2$. Therefore, we conclude that the slow growth of the IBL in this work is of a physical nature and does not depend on the adopted methodology. Sect. 3.3 introduces and discusses a needed modified model to capture the underlying physics.

3.3 A modified diffusion model

To interpret the measured slow growth rate of the IBL when compared to the model prediction, we will have a closer look to the diffusion term. Figure 8 shows the plots of ${\sigma }_w/U$ as a function of $U/U_\textrm{ref}$ along the top part of the IBLs. The dashed line denotes the assumption used in the customary diffusion model, for which ${\sigma }_w/U {\propto } U_\textrm{ref}/U$. For neutral cases, the data within log-law region (see square and round symbols data for $U/U_\textrm{ref}<0.75$) fall onto the predicted relation so that the growth curve of the IBL is well predicted by the original diffusion model (Eq. (6)). However, in the outer layer ${\sigma }_w/U$ decays much faster than $U^{-1}$ leading to weaker diffusion effects. It is found that, for all cases examined herein, the values of $\sigma _w/U$ at the top of the IBL can be well described by a power function of U in the following form:

$$\begin{aligned} \sigma _w/U=\sigma _0({U/U_\textrm{ref}})^{-\beta }. \end{aligned}$$

(7)

The parameters $\sigma _0$ and $\beta$ are found to be a function of the thermal state of the boundary layer. For neutral cases, these parameters assume the values of $\sigma _0=0.039$ and $\beta =1.9$. For stably-stratified boundary layers, the curves of ${\sigma }_w/U$ are much lower than those of the neutral stratification, implying that the diffusion effects are suppressed by the stable stratification. The value of $\beta$ is found to be 2.7 for both cases. The smaller values of $\sigma _w/U$ for $\textrm{Ri}_\textrm{b}=0.27$ are achieved by a smaller value of $\sigma _0$ when compared to those for $\textrm{Ri}_\textrm{b}=0.13$. These fitting parameters are reported in Table 2.

Considering the decay law of vertical diffusion and its departure from the customary diffusion model just discussed, we propose to replace the term ${\sigma }_w/U$ in Eq. (1) with Eq. (7). The comparison between data and the model curve after the proposed correction of diffusion term is shown by grey solid curves in Fig. 9. It is notable that the predicted IBL depth becomes shallower after the correction of the diffusion term, especially for stable stratification cases. However, the grey curves still do not satisfactorily fill the gap between the measured data and the predictions of the modified model.

By analysing the mean flow fields, one can note that there exists vertical advection effects downstream of the step change in surface roughness. These are shown as negative mean vertical velocity in Fig.  5b. Figure 8b presents the ratio W/U at the top of the IBL and how this ratio varies with the spatial gradient of the IBL, $\textrm{d}\delta _i/\textrm{d}x$. With an increase in fetch (i.e. smaller $\textrm{d}\delta _i/{\rm d}x$), the vertical advection effects described by the term W/U becomes stronger. For instance, the magnitude of W/U reaches $40\%$ of the value of $\sigma _w/U$ at $x=5.88$ m for neutral flows and this percentage increases to around $95\%$ for the most stable case ($\textrm{Ri}_\textrm{b}=0.27$). In first approximation, we can assume this effect to be linear and have indicated this in Fig. 8b by the black solid line, where the vertical advection term is modelled as:

$$\begin{aligned} W/U=w_0+\gamma (\textrm{d}\delta _i/\textrm{d}x), \end{aligned}$$

(8)

with constant parameters $w_0=-0.02$ and $\gamma =0.14$. In the above relation, the increase of W/U with $\textrm{d}\delta _i/\textrm{d}x$ is expected (Savelyev and Taylor 2005), whereas $w_0$ plays a key role here. This offset might be dictated by the type of roughness used, for instance, one might expect the value of $w_0\rightarrow 0$ for skimming flow over recessed roughnesses. We note that there is an uncertainty in determining the absolute value of $w_0$ in experiments due to any small misalignment of the LDA beams and the free-stream wind direction. The consequence can be seen in contour plots in Figs. 5 and 6. Careful analysis of the upstream flow field, in particular the vertical profiles of W/U, enabled compensation. To understand the significance of the two different terms (diffusion and advection), we can also combine the vertical advection effects into the modified diffusion model in Fig. 8c. Here, the spatial gradient of the IBL depth, $\textrm{d}\delta _i/\textrm{d}x$, is found to be adequately captured by the solid line, which represents a fitting curve in the form:

$$\begin{aligned} \textrm{d}{\delta }_i/\textrm{d}x=C({\sigma }_w+W)/U, \end{aligned}$$

(9)

with the constant parameter C being 1.58. The decent concurrence of all cases considered herein onto this fitting line implies that this further modified model considering the mean vertical velocity W, should now be able to interpret—and adequately capture—the growth of the IBL (in the outer layer) for both neutrally- and stably-stratified flows. Combining the interpretation and findings discussed in Fig. 8 (via fitting Eqs. (7) and (8)) as well as the upstream profile of mean streamwise velocity into Eq. (9), a final modified model to predict the growth of the IBL, can be expressed as follows:

$$\begin{aligned} {\int _{\delta _i^0}^{\delta _i}}\frac{1-C\gamma }{{C\sigma _0 \left[ \frac{u_*}{\kappa }\left(\ln {\frac{z}{z_0}}+\beta _m\frac{z-z_0}{L_0}\right) \right] }^{-\beta }+Cw_0}\textrm{d}{z}=x-x_0, \end{aligned}$$

(10)

where $\beta _m =0$ for neutral cases and $\beta _m=8$ for stable cases as in Hancock and Hayden (2018). The initial condition $(x_0, \delta _i^0)$ is defined as the quantities at the location where the IBL reaches the upper boundary of the logarithmic region of the approach flow (downstream roughness). Figure 9 confirms how the curves generated from Eq. (10), here shown in blue, can finally well capture the measured growth rates of IBLs, which was over-predicted with the original model (Eq. (6)). Given the relatively crude approximation of the fitting of W/U in Fig. 8b, the exact value of $w_0$ for each case in Fig. 9 was slightly modified to optimally fit the data. The parameters used in the new proposed model are listed in Table 2.

Fig. 8

a $\sigma _w/U$ as a function of $U/U_\textrm{ref}$ along the top of the IBLs. The dashed line is a power function with power exponent of $-1$. The solid lines are data fits using a power function in the form $\sigma _w/U={\sigma _0(U/U_{ref})^{-\beta }}$. The fitting coefficients $({\sigma _0},\beta )$ are included in  2. b W/U along the top of the IBLs as a function of $\textrm{d}\delta _i/\textrm{d}x$. The solid line is a linear fit to the data in the form $W/U=-0.02+0.14(\textrm{d}\delta _i/\textrm{d}x)$. c The plot of $\textrm{d}\delta _i/\textrm{d}x$ as a function of $(\sigma _w+W)/U$. The solid line is a fit to the data using Eq. (9) with $C=$1.58. Symbols represent the same conditions as in Fig. 3

Fig. 9

Comparison of growth curves of IBLs with models. Symbols represent the experimental data with the same conditions as in Fig. 3. The dashed line denotes the predicted curve generated from Eq. (6). The grey curves are generated from the model considering the decay of the diffusion term in Eq. (7). The blue solid curves are generated from Eq.  (10), which considers the decay of the diffusion effect as well as the vertical advection effect. In each set of model curves, from top to bottom are for $\textrm{Ri}_\textrm{b}=$0, 0.13 and 0.27. Grey curves are discussed in Sect. 3.3

Table 2 Fitting parameters for the new modified model

4 Conclusions

This paper investigated the IBL development from an abrupt change in surface roughness experimentally. Both thermally-neutral and two cases of increasingly stably-stratified boundary layers ($\textrm{Ri}_\textrm{b}=0.13,0.27$) were considered in our meteorological wind tunnel. Through comparing two neutral boundary layers at different Reynolds numbers, the IBL growth curve was found to be nearly independent of $\textrm{Re}_{\delta }$ within the range investigated ($2.9\times 10^4<\textrm{Re}_{\delta }<4.5\times 10^4$). For neutral flows, both the depth and the growth rate of the IBL are well predicted by the diffusion model of Panofsky and Dutton (1984) provided that the IBL is embedded within the logarithmic layer of the boundary layer developed over the upstream roughness, and that an appropriate modified origin is considered. However, when the IBL grows into the outer layer of the approach flow, its depth and growth rate are over-predicted by the original diffusion model. The growth is found to be represented by a power function with exponent n being 0.61, which is slower than the model prediction ($n=0.8$). For thermally stable flows, the power exponent in the outer layer was found to decrease with $\textrm{Ri}_\textrm{b}$ number, where n was reduced to 0.46 for $\textrm{Ri}_\textrm{b}=0.13$ and further decreased to 0.39 when $\textrm{Ri}_\textrm{b}=0.27$.

The analysis on the spatial IBL growth, $\textrm{d}\delta _i$/dx, along the top of the IBL show that two causes lead to its slow growth for neutral flows compared to the model prediction. On the one hand, the diffusion term, $\sigma _w/U$, decays with $U/U_\textrm{ref}$ faster than the model prediction ($({U/U_\textrm{ref}})^{-1}$). On the other hand, the vertical advection effects described by $W/U<0$ become significant and comparable to the magnitude of the diffusion term in the outer layer, particularly in the far field from the step change in surface roughness. Within experimental uncertainty, this vertical advection term is found to be nearly independent of the thermal stability. Thus, compared to the neutral flows, the reduced growth rates of IBLs in stable cases are mainly due to the suppression of diffusion effect (smaller $\sigma _w/U$). The suppression of turbulent fluctuations is a well-known effect of thermal stability (Marucci et al. 2018), so this is perhaps not surprising. Based on these observations, we proposed a modified diffusion model which includes an additional term, W/U, which adequately captures the measured growth of the IBLs for both neutrally and stably-stratified flows. This newly proposed modified diffusion model is semi-empirical, and relies on measurements of the mean vertical velocity and its standard deviation, to accurately predict the depth of the IBL developing over step change in surface roughness. It must be stressed that, although the model has been developed based only on the surface/thermal conditions used herein, and therefore the numerical parameters in Table 2 are tailored to these conditions, the methodology introduced in this work is generic and does not rely on particular surface conditions. As such, it should be equally applicable to $S \rightarrow R$, $R \rightarrow S$ and $R_1 \rightarrow R_2$ step changes in roughness. Furthermore, there are strong indications (see Sect. 3) that this methodology could be extended to more stable flows (as the fitting parameters in Eq. (7) are a function of the thermal state of the boundary layer), and that the fitting parameters for neutral conditions are somewhat universal (see values for Case 1 and Case 3 in Table 2) regardless of the roughness conditions.

Further work is needed to develop a truly predictive model, which does not necessitate of detailed velocity measurements. For this to be successful, one would hope to be able to predict the growth of the IBL given appropriate initial conditions as well as the aerodynamic properties of, at least, the upstream surface (e.g. roughness length $z_0$ and friction velocity $u_*$). This, it is hypothesised, should be possible by using a framework similar to that in Li et al. (2022).