Simulating the Unsteady Stable Boundary Layer With a Stochastic Stability Equation

Turbulence in very stable boundary layers is typically unsteady and intermittent. The study implements a stochastic modeling approach to represent unsteady mixing possibly associated with intermittency of turbulence and with unresolved fluid motions such as dirty waves or drainage flows. The stochastic parameterization is introduced by randomizing the mixing lengthscale used in a Reynolds average Navier‐Stokes (RANS) model with turbulent kinetic energy closure, resulting in a stochastic unsteady RANS model. The randomization alters the turbulent momentum diffusion and accounts for sporadic events of possibly unknown origin that cause unsteady mixing. The paper shows how the proposed stochastic parameterization can be integrated into a RANS model used in weather‐forecasting and its impact is analyzed using neutrally and stably stratified idealized numerical case studies. The simulations show that the framework can successfully model intermittent mixing in stably stratified conditions, and does not alter the representation of neutrally stratified conditions. It could thus present a way forward for dealing with the complexities of unsteady flows in numerical weather prediction or climate models.


Introduction
The representation of the atmospheric boundary layer in stably stratified conditions is an intricate problem for numerical weather prediction (NWP) and climate models (Holtslag et al., 2013;Sandu et al., 2013).Stably stratified conditions can occur at nighttime when radiative cooling of the surface is predominant, or when warm air is advected over a cold surface, for example, over snow or ice.Such conditions favor model biases, a prominent example being systematic errors in the near-surface temperature (Davy & Esau, 2014;Esau et al., 2018;Køltzow et al., 2019).The different processes occurring at the interface between the surface and the lower atmosphere interact in complex ways, making the identification of the main source of error challenging.Observational studies highlight very large horizontal variability of atmospheric quantities in very stable conditions, in most cases occurring on scales that are often below typical model resolutions (Mahrt et al., 2021).Model errors have been related to shortcomings in the calculation of turbulent fluxes, radiative fluxes or ground heat fluxes, as well as to an overestimated heat capacity of a too deep boundary layer, preventing a sufficiently fast reaction of the near-surface temperature (Esau et al., 2018;Sandu et al., 2013;Tjernström et al., 2005).
Turbulence in the stable boundary layer (SBL) is generated by shear production, while its development is inhibited by buoyant forces.Due to this interplay, flow regimes with different physical and dynamical characteristics exist (Mahrt, 2014;van de Wiel & Moene, 2003).Fully turbulent SBL, also coined as weakly stable boundary layers, are rather well described by similarity theory, but the very SBL with intermittent turbulence is less well understood (Grachev et al., 2005;LeMone et al., 2018;Mahrt, 2014).At high stability, non-turbulent processes become more important, and the flow is characterized by strong non-stationarity (Mahrt & Bou-Zeid, 2020).For example, larger scale wave-like motions can interact in complex ways and contribute to intermittent turbulence (Cava et al., 2019).Non-turbulent flow features smaller than those traditionally classified as mesoscales, denoted as submesoscale motions, exist under all atmospheric stratifications for weak winds (Anfossi et al., 2005), but exert a critical influence under strong stratification.In these conditions characterized by a large Richardson number, turbulence production is closely related to local short-term accelerations associated with submeso motions (Boyko & Vercauteren, 2020;Lan et al., 2022;Mahrt, 2011).Approaches to parameterize non-turbulent motions are being developed, including the quasi-normal scale elimination (QNSE;Sukoriansky et al., 2005) that includes breaking gravity waves, or a quantification of wave drag due to small scale orography (Steeneveld et al., 2009).Another closure approach is based on the total turbulent energy that considers the potential energy due to density fluctuations of the fluid in addition to the traditional consideration of turbulent kinetic energy (Mauritsen et al., 2007;Zilitinkevich et al., 2007).A unified treatment of non-stationary turbulence in very stable conditions is however lacking (Edwards et al., 2020;LeMone et al., 2018).
With weak winds and clear-sky conditions, associated with strong stability, the atmosphere may become decoupled from the surface (Acevedo et al., 2016).This occurs when a layer near the surface becomes driven by radiation and soil thermal transport, while the surface turbulent heat flux is too weak to sustain the energy demand of the surface (Van de Wiel et al., 2012).In NWP, the decoupling can occur in very localized regions with a high spatial variability, and the positive feedback between weakening turbulence and radiative cooling can lead to further rapid cooling in decoupled regions (Kähnert et al., 2022).To avoid such decoupling and so-called runaway cooling to become unphysically important in models, operational parameterization schemes have implemented rather high levels of turbulent mixing (Cuxart et al., 2006;Derbyshire, 1999;Louis, 1979).This practice is often justified by the need to account for the numerous processes impacting mixing that are not resolved in NWP and climate models, such as unresolved surface heterogeneity or topography, and internal gravity waves.This enhancement of turbulent mixing is typically calibrated to reduce the activity of synoptic systems and improve model scores, with the negative consequence that NWP and climate models simulate too deep boundary layers, too weak low-level jets or wind veering with height (Sandu et al., 2013).Additionally, the dimensionless parameter used to characterize the stability of a rather large-scale numerical grid cell (e.g., the Richardson number) becomes less representative of the local-scale dimensionless stability when the flow is strongly stable, leading to model uncertainty.
In an effort to model the variability of mixing related to intermittency of turbulence, internal or related to submeso motions, Boyko and Vercauteren (2023) devised a stochastic extension to MONIN-OBUKHOV Similarity Theory (MOST) that is able to model intermittent turbulent bursts.The proposed approach keeps the physical basis of MOST untouched, assuming a gradient-diffusion model in which the diffusivity scales with an appropriate lengthscale incorporating the influence of dimensionless stability.It extends MOST by treating the stability correction and thus the mixing lengthscale as a time-continuous stochastic variable, thereby enabling the representation of unsteady mixing.There may be intrinsic limits in such a gradient-diffusion model structure, even when the diffusion coefficient is stochastic, since the validity of similarity theory is questionable in strongly stable cases (Grachev et al., 2013).However turbulence parameterization schemes used in operational NWP models were shown to reasonably capture the physics of the SBLs for a variety of forcing provided they do not apply excessive vertical mixing (Baas et al., 2018(Baas et al., , 2019;;Cuxart et al., 2006).Using tools from uncertainty quantification, Audouin et al. (2021) concluded that model deficiencies reflect a poor parameterization calibration rather than intrinsic limits of the parameterization formulation.These authors further suggested a framework combining single-column models and large eddy simulations to improve the calibration of SBL model parameters.In the observational study presented in Boyko and Vercauteren (2023), the calibration of a proposed time continuous Stochastic Stability Equation (SSE) is analyzed statistically using field observations and inverse modeling methods.The results highlight scaling of the stochastic model parameters with dimensionless atmospheric stability, providing a closed-form parametrization of turbulence that enables explicit treatment of the uncertainty of the fluxes to be modeled.Due to the time-continuous model structure, the proposed stochastic extension of MOST enables the representation of localized bursts of turbulence through a stochastic model.It could provide an alternative to the practice of imposing enhanced mixing, as one could instead keep the average mixing to lower levels justified by observations, but enhance it locally through intermittent bursts.Such a stochastic parameterization of turbulence can also provide much needed uncertainty estimations, and may also be needed to better represent the mean state and SBL regime transitions that can occur via inherent nonlinear processes (Berner et al., 2017;Van de Wiel et al., 2017).
In this study, the stochastic representation of the mixing length is implemented into a REYNOLDS-averaged NAVIER-STOKES (RANS) model.The momentum and heat diffusivity, as well as all the state variables are predicted from a stochastic mixing length according to the SSE introduced by Boyko and Vercauteren (2023).The fact that stochastic perturbations are introduced enables intermittency to be modeled.The study investigates the impact of the suggested stochastic scheme on a range of numerical case studies to evaluate the robustness of the proposed framework.The stochastic model, coined as Stochastically Unsteady REYNOLDS-averaged NAVIER-STOKES Equations (SURANS), is presented in Section 2 and its numerical implementation in introduced in Section 3. Section 4 presents the results of numerical case studies, which include a neutrally stratified boundary layer, a stably stratified boundary layer, followed by a case with variable geostrophic wind and radiative forcing.A summary and conclusions are given in Section 5.

Deterministic Model
For a flat surface in a dry atmosphere, assuming horizontal homogeneity, neglecting radiative flux divergence, applying the BOUSSINESQ approximation, and using a turbulence closure model based on eddy diffusivities (where w′u′ = K m ∂u ∂z , w′v′ = K m ∂v ∂z , and w′θ′ = K h ∂θ ∂z ), the idealized SBL can be represented by the following RANS model (Stull, 1988): where u, v are the mean (Reynolds averaged) horizontal wind components and θ is the mean potential temperature.The horizontal pressure gradient is prescribed through the geostrophic velocity above the SBL, whose wind components are (u g , v g ), and f c is the CORIOLIS parameter.The ground surface temperature θ g is the bottom boundary condition of Equation 3 and its evolution is modeled using a force-restore method (Acevedo et al., 2021;Garratt, 1994;Stull, 1988).The thermal capacity of the soil per unit area is denoted with C g .The soil heat transfer coefficient κ m = 1.18ω is related to the Earth's angular frequency ω.H 0 = ρc p w′θ′ 0 is the surface sensible heat flux, where ρ is the air density and c p is the specific heat of air at constant pressure, and R n is the net radiation.The temperature below the surface θ s at some finite depth is nearly constant and fluctuates on a seasonal scale.It is, therefore, deemed fixed for the simulation of individual nights.
Closing the model requires further specification of the eddy diffusivities K m and K h .Many operational NWP schemes use first-order schemes, in which the eddy diffusivity depends on the wind speed, a specified mixing lengthscale and a stability function (Cuxart et al., 2006).Higher-order schemes add more prognostic equations to the model to compute turbulent quantities.A common choice is that of a 1.5 order closure, in which a prognostic equation is used only for the evolution of the Turbulence Kinetic Energy (TKE), e.In this case, which will be further developed in the following model extension, the eddy diffusivity for momentum is expressed as follow: where l m stands for the momentum mixing length and α is a modeling constant (Cuxart et al., 2006;Rodrigo & Anderson, 2013).In the evolution equation for e, Equation 6, P E represents the production of TKE and ϵ its dissipation rate.Turbulent kinetic energy is produced through wind shear and buoyancy, hence where g = 9.81 m s 2 is the gravitational acceleration, Θ 0 = 300 K is a reference potential temperature.The dissipation rate ϵ is modeled using a dissipation length l ϵ , which is assumed equal to the mixing length in our study, that is, l ϵ = l m , leading to where α ɛ is a modeling constant set to α ɛ = 0.1 in this study (Cuxart et al., 2006;Rodrigo & Anderson, 2013).The turbulent PRANDTL number Pr t = K m K h can be used to obtain K h from K m and in the following, it is set to one for simplicity.A detailed presentation of several operational 1.5 order schemes can be found in Cuxart et al. (2006), where it can be seen that schemes differ in the values selected for the constants, in the parameterization used for the mixing lengths and in the stability functions used to scale the eddy diffusivities according to the static stability.

Stochastic Extension
The model extension implemented in this work, denoted as SURANS model, is developed as a set of prognostic equations for simulating unsteady intermittent turbulent mixing in the SBL.The main difference to the RANS model is a stochastic extension of MOST in the form of a SSE representing the evolution of a stability correction variable.The SSE derives from a data-driven modeling approach introduced by Boyko and Vercauteren (2023) with the goal of modeling the variability of turbulent fluxes due to the influence of unresolved submesoscale motions and more generally to turbulence intermittency.The SSE is limited at this stage of research to the nearsurface boundary layer where field observations were analyzed, and hence the following numerical implementation is meant to serve as a proof-of-concept where the effect of intermittent mixing is modeled to a certain maximum height above the surface.In reality, submesoscale motions tend to be more relevant at higher levels, rather than near the surface, hence extending the modeling approach to higher levels would be desirable.The height-limited implementation is chosen because the SSE was calibrated based on measurements up to 30 m at one field site (Boyko & Vercauteren, 2023), and we chose to leave extensions to higher levels for future research.This choice has impact on the type of nonstationarity that will be simulated by the model, and that aspect is discussed further below.The impact of the selected modeling strategy is analyzed based on selected numerical case studies.The set of equations forming the SURANS model complements the model (Equations 1-6) as follows: The SSE as Equation 14is the novel contribution to the classical RANS model and implements a time varying stochastic stability correction variable ϕ introduced in Boyko and Vercauteren (2023).This model and all variables included in the equation will be introduced in detail further below.Relaxation terms N u = (u u g )/τ r and N v = (v v g )/τ r are added to the momentum equations in Equations 9 and 10, where τ r is the relaxation time.
Those nudge the solution toward the geostrophic wind and are used to damp inertial oscillations that become too important when turbulent mixing is weak, which is likely unphysical.The value of τ r is set in a range of 3-6 hr, such that the solution is largely controlled by Equations 9-14 and only mildly nudged toward the geostrophic forcing (u g , v g ).The prognostic Equation 12describes the evolution of the TKE according to the model introduced in Section 2.1.Next, the gradient Ri number is used: The eddy diffusivities K m = K h are modeled according to Equation 5 with a parameterized turbulent mixing length.The chosen parameterization is similar to the analytical expression suggested by Blackadar (1962) for neutral ABLs, and extended by Delage (1974) to account for stability: where κ is the von Kármán constant and with the difference that φ(t, Ri), which will be properly defined in Equation 20, follows from the SSE and thus is a nondimensional stochastic process that replaces the use of the dimensionless shear in the original formulation (see e.g., Rodrigo & Anderson, 2013; Equation 17).Following Rodrigo and Anderson (2013), the value λ b , which restrains the size of the largest turbulent eddies in neutral stratification is parametrized as: where f c = 2 ω sin(φ), with φ = 40°N and ω = 7.27 × 10 5 s 1 .The stochastic variable φ(t, Ri) is constructed using a mixture of deterministic and stochastic formalism.Equation 14 determines the stability correction value ϕ from the surface up to some chosen height z s (set as z s = 50 m), above which a traditional scaling function ϕ f (Ri) is in operation, here taken as (Cuxart et al., 2006): Finally, the descriptions above and below z s are matched through the logistic sigmoid function: where z s is the sigmoid's midpoint, and k s = 0.1 is the steepness of the curve, which regulates the sharpness of transition from ϕ to ϕ f at the height z s .Then the linear-convex composite is defined: and is inserted into Equation 16.Due to the stochasticity of φ, the mixing length l m and hence the entire turbulence closure become stochastic.The stochastic process accounts for the variation of the mixing length hypothesized to be related to intermittency of turbulence and to submesoscale mixing events (Boyko & Vercauteren, 2023).
The stochastic process ϕ is expressed by the prognostic Equation 14 and represent the randomized stability correction.In this equation, the first term (multiplied by dt) is a deterministic tendency in the differential equation, and the second term (multiplied by dW t ) is the stochastic forcing of the differential equation.The process dW t corresponds to increments of a Wiener process, in other words to white noise.Finally, the data-driven scaling functions (Λ, V, Σ) obtained in Boyko and Vercauteren (2023) scale the model coefficients with the Ri number: V(Ri) = 10 (0.4 log 10 (Ri)+0.2) , ( 22) where σ s (see Equation 23) regulates the intensity of the stochastic component of Equation 14.The parameter σ s can be adjusted in the range [ 1,0].The value σ s = 1 equals to the considerably low intensity of the noise, such that the solution of Equation 14 becomes nearly deterministic.The value σ s = 0 corresponds to the level of the Fluxes Over Snow Surfaces Phase II (FLOSS2) data set and models relatively intense perturbations.All details related to the data-driven identification of the scaling function are given in Boyko and Vercauteren (2023) and are not repeated here.One important aspect is that the model gives an expected value equal to one when the number Ri approaches zero, to be matching the neutral stability case in this limiting condition.However, one issue remains in the limit of higher wind speeds in neutral conditions, and that is that the intensity of the fluctuations becomes too large for larger wind speeds.This limitation is discussed in Boyko and Vercauteren (2023) and further research is needed to find a better solution to ensure correct limiting behavior.As a practical solution, a threshold can be implemented to limit the intensity of the noise.An example realization of the SSE for different levels of σ s can be visualized in that paper, Figure 6.Finally, the data-driven identification of the parameters was done based on hourly time units.The constant τ h = 3,600 in Equation 14 transforms the units of the equation into seconds for the numerical implementation.Consider that due to E[(dW t ) 2 ] = dt, the process dW t has the units of ̅̅̅̅̅̅̅̅̅ time √ (Horsthemke, 1984), and hence the transformation of units for the noise (stochastic) term is different than in the drift (deterministic) term.

Discretization
Equations 9-14 are discretized and solved using the Finite Element Method (FEM) library FEniCS (Alnaes et al., 2015;Logg et al., 2012), which performs the discretization of the nonlinear system using the FEM.Dunbar et al. (2008) also applied the FEM to simulate the SBL and showed that an adaptive grid refinement approach significantly increases the accuracy of the solution.Nevertheless, the adaptive grid technique is not used here.Instead, a fine grid resolution is set and found to be affordable for the single-column proof-of-concept study done here.Equation 13 is discretized with the explicit EULER method in time.The stochastic Equation 14 is discretized with the MILSTEIN method in time (Lord et al., 2014).
Two different numerical grids are used in the discretization.For the variables u, v, θ, and e, a power-three transform on the z-axis is imposed to improve the resolution of the gradients in the vicinity of the surface.Such a non-equidistant grid cannot be used to solve the stochastic Equation 14 due to the sampling algorithm, which utilizes a FOURIER transform.The FOURIER transform is used because the sampling procedure of the noise process uses a correlation lengthscale, such that the random perturbations are correlated in space.The interested reader is referred to Boyko (2022) for full details on this implementation and on the definition of the correlation lengthscale.Furthermore, since the stochastic perturbations are included in the lower portion of the boundary layer (z < 50 m) the stochastic grid is confined to the lower portion of the computational domain.This saves computational resources and improves the vertical resolution of the stochastic perturbations.Figure 1 shows the description of the numerical grids along with the computation steps to obtain the hybrid stochastic mixing length correction φ defined by Equation 20.
The total domain is organized into three sub-layers, as indicated on the left in Figure 1.The stochastic layer reaches up to the height z s = 50 m.In this sub-domain, the dynamic of the stability correction variable is entirely determined by the Stochastic Differential Equation (Equation 14).From the height of z s up to the height 1.2z s < z p < 2.0z s , the stochastic fade-out layer is defined.The layer is responsible for the smooth transition from the stochastic to the deterministic value.The transition layer is also responsible for providing sufficient buffer length needed by the sampling algorithm to obtain random structures which do not re-enter the domain at the surface s 0 .Indeed, without a buffer layer, the stochastic structures would re-enter at the surface due to periodicity assumptions of the Fourier transform used to sample to stochastic process.A linear-convex combination is performed between the stochastic ϕ and the deterministic ϕ f variables (see Equation 20; also marked in Figure 1 with the step 5).The height z s = 50 m characterizes the smooth blending between stochastic ϕ and the deterministic ϕ f variables.Its value is set slightly larger than the measurement tower that was used to calibrate the stochastic part of the model.Hence, only the lowest 50 m of the simulations have a randomized stability correction in the application of MOST.

Initial and Boundary Conditions
Initial conditions are set following logarithmic profiles in neutral conditions, with: where u * ,init = 0.5 C f u 2 g ) 1/ 2 .Here C f ≈ 4 × 10 4 is a tuning parameter and is adjusted such that u g is obtained at the model top.The initial profile for e is estimated following Parente et al. (2011): The coefficients a 1 and a 2 are estimated using the following boundary values:  19).
e(z = z 0 ,t = 0) = u 2 * ,init (0.087) 1/2 , (27) where H is the domain height.The initial profile of the potential temperature is constant Θ 0 = 300 K up to a certain height H c = 200 m and then increases according to the dry adiabatic lapse rate Γ = 0.01 K m 1 as used by Sorbjan (2012): Regarding the boundary conditions, for the wind components no-slip conditions (DIRICHLET condition) are set at the surface, while at the top boundary, the vertical gradients are set to zero (NEUMANN condition).A lapse rate is imposed as upper boundary condition for the potential temperature.The values of parameters of the SURANS model used in the numerical cae studies are summarized in Table 1.

Numerical Case Studies
Idealized numerical case studies are used to test the SURANS model, validate the numerical stability of the proposed stochastic turbulence closure scheme and study the resulting differences to the classical RANS model with a 1.5 order closure.The impact of the stochastic perturbations that induce intermittency and unsteady mixing is analyzed by comparison to the unperturbed model in three numerical experiments differing in stability conditions.The neutral stratification is studied first.This study is a validation case where no stability correction is needed for the mixing length, hence the ensemble mean of the SURANS model should match the RANS model.
Next, the strongly SBL with intermittent mixing is analyzed.The SURANS model reproduces an intermittent TKE state.When analyzing this intermittent state, the ensemble mean is not a representative measure due to non-Gaussian statistics.A more appropriate measure is the central tendency (the value corresponding to the maximum of the Probability Density Function [PDF]), and its evolution is used to evaluate the performance of the models.
Those two studies are performed for a quasi-stationary case, where the geostrophic forcing and the soil properties are constant in time.The stochastic perturbations may also alter the solution under conditions with variable forcing, and this aspect is analyzed in a third numerical study.

Neutral Boundary Layer
As a first numerical experiment, a neutral boundary layer is simulated with the SURANS model.This experiment validates that the central tendency of the SURANS model, that is, the most probable value of an ensemble of realizations, is equivalent to the RANS solution.The initial conditions are set as neutral profiles as described in Section 3 and the simulation period is set to 15 hr.The solver specific settings for this experiment are given in Table 2 and the rest in Table 1.The forcing parameters are set to be constant.The stratification is controlled with two parameters of the surface energy balance implemented in Equation 13, namely the net radiation R n , and the restoring temperature θ s , which together control the degree of surface cooling.To simulate neutral conditions the net radiation is set to 0, hence forbidding radiative cooling.The restoring temperature θ s is set equal to the initial air temperature, ensuring strictly neutral stratification.The distributions of the TKE from the 100 SURANS simulations are close to being Gaussian, but more importantly, those are symmetrical.This symmetry indicates that the modeled type of turbulence is such that the perturbed solutions maintain their path around the central tendency, which itself is very close to the deterministic RANS solution.Hence in this neutral case, the stochastically added effect of unresolved random mixing events is small enough that the TKE remains in statistical equilibrium in the perturbed model.As shown in the next stably stratified experiments, the equilibrium becomes weaker and more sensitive to the perturbations at a larger Ri number, leading to turbulence intermittency.

Stably Stratified Boundary Layer
The next experiment considers a stably stratified boundary layer in the presence of random mixing events.Similar to the neutral case, the initial conditions are given in Section 3, and the simulation period is set to 15 hr.The solver-specific settings for this experiment are given in Table 3 and the rest in Table 1.The forcing parameters are set to be constant.The stratification is imposed with two mechanisms, the first being the difference between the restoring (soil) temperature of 290 K and the potential temperature of the air 300 K, and the second being a radiative cooling enhancing the stratification.The net radiation of 30 W m 2 is selected following Acevedo et al. (2021) and considered as the FLOSS2 data set average value.This setup may describe a typical cloud-free night in springtime.Note.The parameters marked with "-" are given individually in the following case studies.
Figure 3 illustrates the solution of the SURANS model.The TKE at the height z = 20 m is compared against the solution of the RANS model using several statistical metrics.In Figure 3a, a characteristic signature of intermittent TKE simulated with the SURANS model is highlighted in blue.Thin gray lines display other realizations of the stochastic model.Note the two different types of spikes found at t = 6 hr and t = 10 hr.Their magnitude is significantly larger than the ensemble mean (solid yellow) and the central tendency (solid red).The duration of these events is approximately 1 hr and falls within the characteristic range of sub-mesoscale motions (Mahrt, 2014;Vercauteren et al., 2016).A qualitative relationship between the duration of intermittent events with higher turbulent intensity and submesoscale wind velocity fluctuations was highlighted in a previous observational study (Boyko & Vercauteren, 2020).
The ensemble mean TKE of the simulations, shown in Figure 3, is slightly above the RANS prediction.However, the central tendency is significantly smaller and indicates that it is likely to observe an absence of turbulent mixing.The heavy tail in the ensemble distributions is significant and related to sporadic rare events.Some realizations of the model (not shown) predict a low TKE level throughout the entire simulation period.The wide variety of TKE signatures highlights the representative capabilities of the stochastic model.The central tendency is estimated based on the TKE distribution obtained through 100 model runs at t = 14 hr and shown in Figure 3b.The solid black line represents the prediction of the RANS model for comparison.The solid yellow line is the ensemble mean of the SURANS model, and the solid red line is the central tendency.The central tendency is estimated from the PDF, which is fitted to the histogram by applying the Kernel Density Estimation (KDE) method (Scott, 2015).The estimation is poor and violates the boundary condition on the left side.Nevertheless, the KDE is a time-efficient method to approximate the most probable value.The  Journal of Geophysical Research: Atmospheres 10.1029/2023JD039370 histogram indicates a smaller value of the central tendency than the estimated one.A better estimation can be achieved if a specific distribution type is assumed.However, one should keep in mind that the distribution type is influenced by the stratification.This dependence makes the fitting task less trivial and we refrain from using more complex estimation approaches for studying the distributions of the TKE.
Figure 4 shows a selected realization of the ensemble of simulations including clearly intermittent features.The largest intensity of each burst of TKE is found at the surface.The stochastic correction of the turbulent diffusion can in principle lead to intermittent patches detached from the ground (see Figure 6 in Boyko and Vercauteren (2023) for such an example), as is found to occur in observations (see e.g., Sun et al., 2002).Still, in the simulation we cannot find any turbulent patches that are clearly detached from the surface.The bursts are absent aloft because the turbulent diffusion is multiplied with the gradient of the mean wind, and hence the spatial distribution of the TKE is intrinsically constraint by the wind gradient.A slight inclination (as somebody brushed it from left to right) in the bursts is also present.Some events show that turbulence is still maintained away from the surface (see Figure 4a t = 3.5 hr and t = 5 hr), leading to TKE that is decoupled from the surface.Here the flow is forced with a steady mean wind.Changing the forcing changes the gradient away from the surface and could provide room for the stochastic perturbations to appear at higher levels due to localized shear accelerations.In fact, intermittent turbulent patches could either originate at the surface and move upwards, or could be triggered aloft and propagate toward the surface (Allouche et al., 2022;van der Linden et al., 2020).The chosen restriction of the stochastic forcing to the lower levels of the model likely biases the solution toward the surface-generated origin of patches, and extending the use of a similar stochastic modeling

Journal of Geophysical Research: Atmospheres
10.1029/2023JD039370 approach to higher levels above the surface would allow to study the propagation of events originating at different heights.This question is outside the scope of the current proof-of-concept study, but it should be clear that the restricted modeling choice impacts the type of intermittent events that can be simulated.
The impact of the randomized model on the temperature evolution is visualized in Figure 5.It is evident that in the case of stochastic perturbations, the mixing is performed faster.The mixing rate is higher, leading to a temperature profile which is better mixed, and the temperature inversion is also shifted up and is less abrupt.The temperature profile changes its shape in an unsteady way (compare Figures 5a and 5b), related to the activity of the intermittent burst periods (see Figure 4a).The RANS model on the other hand shows a sharp drop in temperature away from the ground, preventing heat exchanges with the higher levels.The temperature inversion is thus qualitatively different in both models.Quantitatively, the temperature obtained by the RANS simulation at the ground after 14 hr is lower by approximately 3 K compared to the SURANS.This shows that the stochastic perturbations prevent the runaway cooling occurring at the surface.
The profiles of the dominant wind velocity component u are visualized in Figure 6.A repeating pattern of the TKE bursts is visible in Figure 6a, comparable to the pattern seen in Figure 4a.The dominant stochastic turbulent diffusion dictates the boundary layer shape as a consequence of random mixing events.Figure 7 shows the evolution of the profiles for the Ri number.The SURANS model predicts a strongly unsteady local Ri number, but the bulk Ri number is computed for the layer between the z 0 level and z = 80 m.Deviations are found during random mixing events when the temperature profile is mixed sporadically, reducing the local bulk Ri number (compare with Figure 5).Journal of Geophysical Research: Atmospheres 10.1029/2023JD039370

Variable Geostrophic Wind and Net Radiation
In the last case study, a time varying forcing scenario is considered, thereby studying the impact of the stochastic perturbations during transient states.The initial conditions are given in Section 3, and the simulation period is set to 30 hr, which is longer than the average nighttime.In this experiment, the focus lies on computing the transitions between weakly and strongly SBL, as in Acevedo et al. (2021) and Maroneze et al. (2019).The novelty of this study is that random mixing events are included in the model, representing unresolved features of the flow.The nonstationary forcing is chosen such that the geostrophic wind increases gradually at some given time, while the radiative cooling increases once from 0, to go back to a 0 value later in the simulation.The simulation thereby covers four possible forcing combinations, alternating in time as shown in Figure 8a.The solver-specific settings for this experiment are given in Table 4 and the rest in Table 1.
The temporal evolution of the TKE at the height of z = 9 m is shown in Figure 8b, with additional exerts showing the profiles for the variables e, θ and U = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅ u 2 + v 2 √ at three different times (note the arrows in Figure 8b).The quantities visualized in Figure 8 are: • The 100 realizations of the SURANS model (gray, thin lines).
• The central tendency (solid red), estimated as the most probable value from the fitted distribution (see Figure 3).• The noise-free limit of the SURANS model.In this case, the stochastic equation is solved once with a sufficiently low value of the noise, such that the dynamical evolution can be considered deterministic (solid yellow line).The noise-free limit is introduced to eliminate the effect of the difference between the MOST stability function and the deterministic steady-state of the prognostic Equation 14(the expected value of the  random variable).One realization of the SURANS model is emphasized to highlight the rare events during the stable low-wind conditions (solid black line).• The prediction of the RANS model (solid blue line).To study the impact of the applied perturbations, we first compare a solution of SURANS in the noise-free limit with the central tendency estimated from the 100 realizations of the stochastic model (see yellow and red lines in Figure 8).There are no significant differences in the TKE (see Figure 8b and the corresponding profiles).However, there is a substantial impact of the applied perturbations on temperature and velocity profiles.With stochasticity, the temperature is mixed more effectively during the stably stratified period and the mixing extends above the average boundary layer height (see Figure 8 panel ( 1)).The central tendency of the velocity profile experiences a deceleration compared to the noise-free limit.For higher geostrophic winds in the second visualized period, the perturbation of the turbulent diffusion is propagated to the top of the boundary layer (200 m), although the actual perturbations are limited at 50 m.
As a next step, we compare the results of the SURANS and the RANS models.The RANS solution (blue line) predicts higher levels of TKE than the central tendency (red line) obtained by the SURANS model, throughout the entire simulation.Despite this lower level of TKE simulated by the SURANS model, transport of temperature and velocity is enhanced (see Figure 8 panels (1)).This nontrivial effect may result from non-equilibrium statistics in the stochastic formulation of the turbulent mixing length.The variability of results is visualized through the gray area in Figure 8, representing the 0.05-0.95quantile range of the 100 different model runs.For stable stratification (see Figure 8 panels (1)), the quantile range for the TKE is asymmetrical, showing the largest spread closest to the surface.In neutral conditions, the quantile range is symmetrical (see Figure 8 panels ( 3)).The model ensemble spread for the TKE profile is significantly different than the ensemble spread for the temperature and velocity profiles.The largest ensemble spread for temperature and velocity profiles is found in the middle of the boundary layer, with lower spread at the surface and the boundary layer top.
Observing the individual simulation paths (see Figure 8b thin gray lines), the impact of the random perturbations on the transition periods can be analyzed.The inset in Figure 8b highlights a transition from neutral to stable stratification induced by the onset of radiative cooling.The central tendency and the noise-free limit of the SURANS model overlap during the transition.However, multiple individual realizations (thin gray lines) show a pronounced tendency to delay the transition rather than induce early transition.In contrast, by transitioning from low wind to high wind (see Figure 8b from t = 15 to t = 20), the solution paths can show both early and delayed transitions.The individual simulation paths also show that during this period where u g increases, the variance in the TKE increases as well.When radiative cooling is interrupted (see t = 25 hr), the variance reduces to some lower value.The reason for this is the parametrization of the noise term in the stochastic equation, which only scales with the Ri number.This scaling was identified by Boyko and Vercauteren (2023), but possibly other dependencies could be investigated.The relative differences in space and time of solutions obtained through the SURANS and RANS models are shown in Figure 9, where panel a shows the differences in TKE.The transition from blue to white color (no difference) indicates approximately the boundary layer height.The boundary layer grows after t = 15 as the geostrophic wind is increased.For the time t > 25 hr the radiative cooling is interrupted, and the central tendency of the SURANS model becomes very similar to the RANS solution.For the time t > 6 hr the value of the TKE predicted by the SURANS model is 50% smaller than predicted by RANS on average, indicating a shallower boundary layer (as seen in the TKE profiles of Figure 8).Figure 9b shows the relative difference in the temperature.Within the boundary layer (where relative differences in TKE are found) the differences between SURANS and RANS are insignificant.At the boundary layer top, the SURANS model deviates from the RANS solution.For the stably stratified conditions (t ∈ (5, 15)), the central tendency of the SURANS solution predicts almost a 200% lower value of the temperature than the RANS model for a large area above the boundary layer (see the blue area in Figure 9b).At the same time, the differences at the surface are relatively small.This can be explained by the enhanced transport due to intermittent turbulence.By construction, the stochastic perturbations start to fade away above z > 50 m.At the same time, the boundary layer height is approximately 25 m, such that the stochastic perturbations determine the mixing of temperature.The red area at the top of the boundary layer in Figure 9b for t > 17 hr (the high-wind regime) means that the central tendency of the SURANS model is predicting an increased value of the temperature relative to the RANS model.Hence, the errors produced in the stable Reference potential Temperature (K) Θ 0 300 Net radiation (W m 2 ) R n Variable (see Figure 8a) Geostrophic wind (m/s) u g Variable (see Figure 8a) regime (5 < t < 15 hr) are propagated into the high-wind regime (t > 20 hr) at the boundary layer top.This findings suggest that the altered transport of temperature and possibly moisture (although not included in this model) may impact the creation of clouds in the early morning with increasing geostrophic winds.

Summary and Conclusions
A SSE, suggested by Boyko and Vercauteren (2023) to introduce a stochastic parameterization of unsteady turbulence, was implemented and tested in this study.The previous data-driven analyses showed that the stochastic model for turbulent mixing could in principle accommodate both the short-term intermittent behavior of turbulence and the long-term averaged mixing, as validated against field measurements.The stochastic model parameters in the SURANS model were found to scale with the local gradient Ri number (Boyko & Vercauteren, 2023).As a result, the intermittent statistical properties of the modeled TKE are changing continuously as a function of flow stability.In this paper, the stochastic parameterization was implemented in a SURANS singlecolumn model extended from a RANS model with 1.5 closure.The SSE can in principle also be used in a firstorder closure model.The impact of the randomized model was evaluated through selected idealized numerical case studies with varying stability conditions.In the current implementation, the stochastic equation is confined to the lower portion of the boundary layer and is blended with a deterministic model above.It is unknown at this stage if the proposed closure is locally valid in the outer boundary layer.
The proposed framework was found to be numerically stable.In the strongly stable condition it is advisable to use an adaptive time stepping in the time integration to avoid abrupt numerical instabilities.These instabilities come from the strong stratification in combination with the stochastic events.Due to the randomness of the stochastic events it can happen that negative TKE is induced.Any mechanism preventing the solver to run negative TKE values is necessary for strongly stable conditions.
In neutral conditions, the stochastic parameterization was found not to have a significant impact on the statistical properties of the modeled flow, simply introducing limited variability compared to the RANS reference model.Within the regime of strong stratification, the SURANS model adequately represents intermittent TKE patterns.The intermittent mixing events affect the boundary layer height.In conditions of weak stratification and large geostrophic wind speeds, the SURANS model appears to show unrealistically large variance, indicating that further model tuning may be necessary.For practical application it is advisable to limit the noise intensity in the SSE by some critical geostrophic wind, for example.In stably stratified conditions, the SURANS model shows enhanced mixing properties in comparison to a RANS with a linear stability correction function.The temperature profile is mixed faster and reaches over larger heights.In comparison to the RANS solution, the stochastic model predicts lower temperature value just above the shallow, stably stratified boundary layer.The effect of stochastic diffusion reaches beyond the limiting height of the perturbations.This results in qualitatively different profiles compared to the RANS solutions in the outer boundary layer.Furthermore, the boundary layer height becomes highly variable in strongly SBL and is determined by the random turbulent mixing events.
The presented SURANS model shows the potential to be used as an exploratory or even predictive tool.To investigate the use of the SSE for less idealized setups, future studies should validate the performance of the SURANS in controlled case studies using observational data.Future research is needed to investigate if the SURANS approach can be used as an alternative to the implementation of rather high levels of turbulent mixing to avoid the problem of runaway cooling in strongly stable conditions.

Figure 1 .
Figure 1.The computation of the hybrid stochastic mixed length correction φ using two different grids.The z-grid in red is non-equidistant and is used to solve the variables u, v, θ, e.The s-grid in blue is equidistant and is used to solve the stochastic variables ϕ.The circled numbers below mark the five steps to calculate the value of φ. (1) Calculate the Ri number on the grid z. (2) Interpolate the Ri number on the equidistant s-grid.(3) Evolve the stochastic variable ϕ to the next time step by solving the Stochastic Stability Equation.(4) Interpolate ϕ to the non-equidistant grid within the height z p .(5) Compute the linerconvex combination between the deterministic ϕ f and stochastic ϕ variables on the z grid using the sigmoid function sig(z) (see Equation19).

Figure 2
Figure2shows the comparison of the TKE at three different heights (z = 0.5, 70, 150 m) for simulations with and without the stochastic mixing induced by the SSE.A quasi steady-state solution is reached approximately after 6 hr with the RANS model.The central tendency of the SURANS model, which is estimated from averaging over 100 realizations, is nearly identical to the solution of the RANS model.The regularity of the sample paths (indicated with the thin colored lines) varies across the height.More rapid fluctuations are found closer to the surface (sample paths in gray), and smooth oscillations with smaller variances occur at z = 150 m (sample paths in green).The stochastic mixing length equation is only active up to the height z = 100 m.As indicated by the sample paths in green (z = 150 m), the variability induced at the surface is propagating into the upper levels of the boundary layer.Hence the stochastic MOST impacts the upper boundary layer.

Figure 2 .
Figure 2. Comparison of the predicted Turbulence Kinetic Energy (TKE) by the Stochastically Unsteady Reynolds-averaged Navier-Stokes Equations (SURANS) and Reynolds average Navier-Stokes (RANS) models in the condition of neutral stratification (Ri = 0) for three heights (z = 0.5, 70, 150 m).The evolution of TKE is shown in (a) and the corresponding color legend is given in (b).Panel (a) shows the RANS solution with a solid black line.The many lines in different colors indicate the 100 realizations of the SURANS model for their heights.The central tendency of the SURANS model is indicated by a dashed red line.The respective probability distribution of the TKE ensemble at t = 14 hr is given in panel (b).

Figure 3 .
Figure 3.Comparison of Stochastically Unsteady Reynolds-averaged Navier-Stokes Equations (SURANS) and Reynolds average Navier-Stokes (RANS) models predicted Turbulence Kinetic Energy (TKE) under the condition of strongly stable stratification (mass Ri ≈ 0.6) for height z = 20 m.For the visualization of the Ri number profiles, see Figure 7.The evolution of the TKE is shown in (a).The ensemble distribution of 100 sample paths of the SURANS model at t = 14 hr is shown in panel (b) along with the fitted probability density function (solid gray line) using a Kernel Density Estimation method.The thin gray lines show the 100 realizations of the SURANS model.

Figure 4 .
Figure 4. Temporal evolution of the profiles of Turbulence Kinetic Energy (TKE) for a realization of the Stochastically Unsteady Reynolds-averaged Navier-Stokes Equations (a) and Reynolds average Navier-Stokes (b) models.The color bar applies to both panels.

Figure 5 .
Figure 5. Temporal evolution of the temperature profiles for one realization of the Stochastically Unsteady Reynoldsaveraged Navier-Stokes Equations (a) and Reynolds average Navier-Stokes (b) models.The color bar is valid for both panels.

Figure 7 .
Figure 7. Temporal evolution of the Ri number profile for a realization of the Stochastically Unsteady Reynolds-averaged Navier-Stokes Equations (SURANS) (a) and Reynolds average Navier-Stokes (RANS) (b).The color bar applies to panels (a) and (b).Panel (c) shows the bulk Ri number calculated from the surface to z = 80 m.

Figure 6 .
Figure 6.Temporal evolution of the wind profiles (u component) for one realization of the Stochastically Unsteady Reynoldsaveraged Navier-Stokes Equations and Reynolds average Navier-Stokes models.The color bar is valid for both panels.

Figure 8 .
Figure 8. Solution of the Stochastically Unsteady Reynolds-averaged Navier-Stokes Equations (SURANS) model with variable forcing parameters R n (net radiation) and the geostrophic wind u g .The total simulation period is 30 hr.The nudging time scale is set to 5 hr.Panel (b) shows the evolution of Turbulence Kinetic Energy (TKE) at 9 m for 100 realizations, marked with gray lines.The zoom area highlights the transition to stable stratification in weak winds by increasing net radiation.The evolution of the forcing is shown in (a).The sub-images in (b) show profiles of the variables at three different times marked with black dots in (b).The SURANS profiles represent the central tendency, with the gray area showing the quantile range.The boundary layer height z bl used for normalization is 50 m (first period), 200 m (second period), 240 m (third period).

Figure 9 .
Figure 9.The relative difference in profiles between the Stochastically Unsteady Reynolds-averaged Navier-Stokes Equations (SURANS) (the central tendency of 100 realizations) and the Reynolds average Navier-Stokes (RANS) model according to the numerical study in Figure 8.The red color denotes the area where the variables of the SURANS model have larger magnitude than those of the RANS model.The white color denotes no differences.Panel (a) shows the Turbulence Kinetic Energy (TKE), (b) the temperature ( 2 [K] < (T 300) [K] < 0 [K]), and (c) the dominant u component of the wind.The forcing variables change with time and are shown in Figure 8a.Condition of stable stratification for t ∈ (5, 15) hr and condition of high wind pressure for t ∈ (5, 15) hr.The equation used to compute the differences is: (SURANS RANS)/ |RANS|.

Table 1
Summary of the Parameter Values of the Stochastically Unsteady Reynolds-Averaged Navier-Stokes Equations Solver

Table 2
Relevant Solver Settings for the Numerical Study of the Neutral Layer

Table 3
Relevant Solver Settings for the Numerical Study of the Stably Stratified Boundary Layer

Table 4
Relevant Solver Settings for the Numerical Study With Unsteady Forcing Variables