An idealized description for the diurnal cycle of the dry atmospheric boundary layer

We present a conceptual model for the diurnal cycle of the dry atmospheric boundary layer (ABL). It may serveasaframeworkforfuturenumericalstudiesonthetransitionaldynamicsthatcharacterizetheABLoverland.Theconceptualmodelenablesustodeﬁneexpressionsforrelevantphysicalscalesasafunctionofthe most prominent forcing parameters and the low degree of complexity facilitates a dimensionless description. This is useful to help generalize boundary layer dynamics that occur on a diurnal time scale. Further, the model’sapplicationfornumericalstudiesisillustratedhereinwithtwoexamples:asingle-column-modelstudythatassessestheeffectofwindforcingonthemaincharacteristicsofthediurnalcycle,andalarge-eddy-simulation study on the daily evolution of turbulence under weak-wind-forcing conditions. The results from these studies sketch the general evolution of the present set of diurnal-cycle systems in more detail. We discuss how the setups are able to reproduce well-known dynamical features of the ABL and also highlight lim- itations, where the simple conceptual system is unable to describe realistic ABL behavior. We conclude that the present conceptual model has an interesting balance between model-system complexity and physical realism, such that it is useful for future idealized studies on the diurnal cycle of the ABL.


Introduction
Apart from the Antarctic summer and wintertime, the atmospheric boundary layer (ABL) over land is characterized by a diurnal cycle. As such, a thorough understanding of the boundary layer dynamics on a 24-h time scale is highly desirable. A classic result on the diurnal evolution of the ABL structure is presented in the well-known book of Stull (1988). The seminal depiction of the ABL's diurnal evolution is powerful for didactic purposes and provides a context for indepth research on the various ABL archetypes. For example, the daytime convective boundary layer (CBL) and the nighttime stable boundary layer (SBL) are often studied separately, due to their dynamical differences. As such, conceptual and quantitative models for the ABL typically assume (quasi) stationary dynamics. This approach has proven successful and has yielded detailed knowledge on the complex dynamics that are encountered Denotes content that is immediately available upon publication as open access. in the ABL. Further, studies focusing on the earlymorning transition (e.g., Angevine et al. 2001;Beare 2008) and late-afternoon transition (e.g., Grimsdell and Angevine 2002;van Heerwaarden and Mellado 2016) aim to bridge the gap in our understanding that lays between convective and stable boundary layer conditions. However, proper understanding of the interactions between even the most prominent processes that govern the dynamical evolution of the ABL during a diurnal cycle, including the transitional dynamics, remains elusive (Lothon et al. 2014;LeMone et al. 2019). For this purpose, we build upon the idealized diurnal cycle as presented in Stull (1988) and introduce a simplified description of the governing processes. The additional value is that the concept yields quantitative expressions for relevant physical scales and a dimensionless framework that may be useful for future investigations. Following the conceptual development, we emphasize numerical modeling as a tool to study the dynamical behavior of the simplified framework in more detail. We conduct exemplifying numerical studies that illustrate the potential for more in-depth research via this route.
When the aim is to accurately describe local weather patterns, all relevant physical processes need to be taken into account. Alternatively, in order to gain understanding of the dynamics, idealized model cases have been studied in virtually countless works regarding the ABL. The decrease in the model-system complexity allows to better distinguish the interactions between the processes that remain represented in the model. However, interpretation of the results requires special attention when placing them in the context of reality and usually comes with disclaimers regarding the validity of the assumptions that were made. As such, both conceptual and numerical models have to strike a balance between physical realism and model-system complexity and they typically need to be tailored for a specific research purpose. Herein we introduce such a simplified model, aiming to capture the key dynamical aspects of a diurnal cycle of the dry ABL. We introduce the set of parameters that govern the dynamics, and formulate expressions for relevant characteristic physical scales such as boundary layer height, wind speed and temperature. Furthermore, the low degree of complexity enables to describe the system with few dimensionless parameters.
Apart from the a priori analysis, it is tempting to study the dynamics of the present conceptual model for the diurnal cycle in more detail. Herein we aim to do so using numerical methods, where the computer model of the ABL is setup following the present conceptual/simplified descriptions of the relevant forcing mechanisms. More specifically, we conduct a single-column-model (SCM) study on the effect of the pressure-gradient forcing and study the evolution of the atmospheric turbulent structures using a turbulence-resolving large-eddy simulation (LES). These preliminary numerical model studies illustrate how the present concepts can form a useful framework for future model studies on the diurnal cycle. The results exemplify how relevant and wellknown features of the ABL are dynamically retrieved by the model setup. Further, we also highlight some of the limitations of the present concepts with respect to their ability to describe reality.
From the numerical-modeling perspective, we argue that there exist a bifurcation regarding the description of the surface boundary in model scenarios. On one hand, there are studies at the idealized side of the spectrum, that artificially prescribe the boundary conditions of the ABL regarding the thermodynamic variable at the underlying surface via the evolution of the surface sensible heat flux or the evolution of the surface temperature. These studies typically employ turbulence resolving methods such as LES or direct numerical simulation (DNS) (e.g., Sullivan et al. 2016;Haghshenas and Mellado 2019) techniques. Of course, in reality, both surface flux and temperature are internal variables, that will actively respond on atmospheric dynamics itself. Rather than prescribing the surface conditions beforehand, there are other numerical studies that assess the surface temperature and the heat flux as parameters that are integrally part of the soil-atmosphere continuum. This requires the evaluation of the various terms in the surface energy budget (SEB) equation, which in turn increases the model's complexity by either resolving more physical processes or introducing more extensive closures for their parameterization. The advantage being that the ABL can then develop according to the external forcings of the system. Such models typically rely on Reynolds-averaged Navier-Stokes (RANS) modeling (Reynolds 1895), either in a full three-dimensional (3D) numerical model (e.g., Greve et al. 2013) or using a singlecolumn model (e.g., Baas et al. 2018). The more realistic setups, which aim to include all the relevant processes for weather prediction, have also become the realm of LES (e.g., Bertoldi et al. 2007;Liu and Shao 2013;Heinze et al. 2017;Maronga and Bosveld 2017). In this work we introduce a simple SEB that takes the midway between both approaches. The setup is simple and does not require the model to resolve or parameterize much physical processes other than turbulent transport. This is advantageous as it means that the setup can be readily run with existing numerical-solver codes and has limited degrees of freedom. On the other hand, it does not require a prescribed formulation of the heat flux or temperature at the surface, which is interesting from a physical point of view. This paper delineates between the description of the conceptual model, including the a priori analysis, and the two numerical studies. The description of the system and the accompanying idealizations are described in section 2a. This scenario is a priori analyzed, resulting in a conceptual model, the identification of relevant scales and the dimensionless groups that describe the system in section 2b. Next, numerical-method-specific details and choices for the two aforementioned approaches are discussed in section 3. The numerically obtained results are then presented and discussed in section 4. Finally, the conclusions of the work herein are presented in section 5. This is supplemented with the appendix that details on the performed numerical simulations, and lists links to online locations that present a precise description of the methods that were used.

A simplified model for the diurnal cycle and its nondimensional representation
In this section, we discuss the idealized description of the physical forcings and the initial conditions that govern the present conceptual model for the diurnal cycle. Next, the scenario is a priori analyzed and casted into a dimensionless form. Thereafter, the choices that are specific to the numerical methodologies are presented for a singlecolumn (RANS) model and the LES study. Furthermore, some considerations for a DNS are presented. Although, this work does not present any results obtained from DNS, those considerations are relevant as they illustrate how the chosen numerical method can introduce methodspecific dimensionless groups that govern the dynamics.

a. Forcing and initial conditions
The forcing of air mass in the ABL is governed by large-scale ½O (1000) km horizontal circulations and their dynamics. Herein, we only aim to model the ABL over a comparatively small horizontal extent and hence we need to parameterize and idealize the mean wind forcing. Following countless numbers of works, the most prominent wind-forcing mechanism is described as a constant horizontal pressure gradient (2= h P). Furthermore, we include the effect of background rotation with respect to the Coriolis parameter f. With a reference air density r and for nonvanishing f, we may conveniently write the pressure-gradient forcing as a velocity vector U geo that is known as the geostrophic wind, wherek is the unit vector in the direction of the zenith. The velocity components in the two lateral ({x, y}) and vertical (z) directions are labeled {u, y, w}, respectively, and are initialized in geostrophic balance: fu, y, wg 5 fU geo , 0, 0g, where U geo 5 kU geo k, meaning that the coordinate system is chosen such that the pressure gradient points in the negative y direction for positive values of f. The persistence of a constant-with-height and constantin-time pressure gradient and a geostrophic balance in the free atmosphere is part of the model idealization and does not represent the general case of momentum forcing in the ABL. Furthermore, we state explicitly that the velocities at the surface vanish due to the small, yet finite molecular viscosity of the air in the atmosphere. However, the implementation of the surface boundary is argued to be a method-specific feature and the associated treatment is discussed for the different numerical approaches separately in section 3. The forcing of the thermodynamic variable in the ABL is considered to occur at the underlying surface and is described via the sensible heat flux H. In reality, the partitioning of the available energy among the various terms of the SEB remains a topic of discussion, and is a major challenge for modeling. Therefore, in view of the desired simplicity of our description, we adopt a simplified version of the SEB were we only consider the net radiation Q n , the soil heat flux F soil and the sensible heat flux, where Q n is defined positive toward the surface, and F soil and H are defined positive when directed away from the surface, into the soil and atmosphere, respectively. For generality, we use buoyancy b as our thermodynamic variable, and it is related to the potential temperature u according to where g is the acceleration due to gravity and u ref a constant reference potential temperature such that u 21 ref is a sufficient approximation of the thermal expansion coefficient of the (dry) air in the atmosphere. We reexpress the SEB Eq. (3) for the corresponding buoyancy fluxes. We write Q*, G, and B as the buoyancy-flux equivalents to replace the heat fluxes Q n , F soil , and H, respectively: These buoyancy-flux terms differ by a factor of u ref rC p g 21 from their heat-flux counterparts, where rC p is the heat capacity of air at constant pressure, with (constant) density r. In this work, we prescribe the diurnal evolution of Q*(t) and calculate G at the surface based on a simple closure so that B can be evaluated as a residual. For the soil buoyancy flux G, we adopt a ''lumped-parameter view'' (van de Wiel et al. 2017). The concept is that there exists a negative feedback of heat based on the surface temperature T s . The main assumption is that this term can be described using an effective feedbacktemperature scale T d and coupling strength l such that the feedback heat flux F can be expressed according to Here l represents the effective combined coupling strength of the surface temperature with T d due to conductive and radiative processes in the soil and atmosphere, respectively, and hence the name ''lumped parameter.'' In this work we only use the concept to model the soil heat flux (F soil 5 F) such that l and T d are model parameters that only concern the soil. This corresponds to a zero-layer model of the soil. The temperature scale T d can be interpreted as the temperature within the soil at the depth d where the diurnal temperature variations are damped out, and l represents the effective soil and vegetation conductivity. Following Eq. (6), G is expressed as where L is the coupling parameter for buoyancy corresponding to l, b surf and b d are the buoyancy-equivalentto-temperature of T s and T d , respectively; see Eq. (4). In practice, heat storage effects (i.e., history), the different properties of the soil layers, the vegetation properties, the phase changes of water, the water-table depth and the actual temperature variations within the soil are important aspects that determine the response of the soil heat flux as a function of T s over time. To stress our simplification, Fig. 1 shows the diagnosed l from data obtained at Cabauw, the Netherlands (van Ulden and Wieringa 1996). Parameter l is calculated for each 10-min period based on Eq. (6), where we have used the soil temperature at a depth of 30 cm below the surface, the surface (skin) temperature and the heat flux at 5 cm below the surface as T d , T s , and F soil , respectively. The data show that the concept of a constant-in-time lumped parameter value is not consistent with the observations. Not only does there appear to be a wide spread in the diagnosed l values, there appears to be a distinctive diurnal pattern (which may be influenced by phase-lag/history effects). The authors stress that our simple evaluation of G is not an alternative to a more elaborated description of the soil physics in models (e.g., Chen et al. 1997;van Tiggelen 2018). We argue that it is a specific shade of gray in the balance between modelsystem complexity and physical realism, and the results presented in section 4 show that it yields model results with some degree of realism and also highlight some prominent limitations.
The temporal variation in the net radiative flux at the surface Q n , which is herein considered to be horizontally homogeneous, is the main driver of the diurnal cycle. Flux Q n is often expressed in models as the sum of its longwave and shortwave, downwelling and upwelling components. These four terms can be modeled individually using parameterized descriptions for radiation. However, in this work we model a dry, cloud-free day and do not consider radiative divergence within the atmosphere. Due to the smooth and well known characteristics of the evolution of the solar zenith angle, and the limited temperature range within a day, we assume that Q n has some generic features that may be exploited to arrive at a simple prescription for its evolution in time t. Therefore, we plot the diurnal evolution of the net radiation for a clear-sky day at Cabauw, the Netherlands, as was observed on 17 August 2016, in Fig. 2. The well-known diurnal pattern is characterized by positive values that reach a maximum during daytime, and during the nighttime, Q n is more constant in time with negative values. Based on this prototypical evolution we choose to prescribe Q* 5 Q*(t), according to the following functional form: where B 0 and B 1 are buoyancy-flux scales for the daytime and nighttime, respectively (with B 0 . 0 . B 1 ), max [a, b] is an operator that selects the maximum value between the dummy variables a and b, and T is the time scale associated with the periodicity of the diurnal cycle (i.e., T 5 24 h on Earth). This formulation implies that the daytime section of the cycle (i.e., Q* . 0) is exactly equal to the nighttime part of the day (i.e., Q* , 0). This is typically only accurate for a specific time of the year depending on, among other variables, the location's latitude. For the initialization of the buoyancy field, we consider a stably stratified atmosphere associated with the Brunt-Väisälä frequency N (cf. Garcia and Mellado 2014;van Heerwaarden and Mellado 2016;Mellado et al. 2017 where z is the height above the surface. Note that the set of buoyancy scales {b d , b surf,t50 } as used in Eqs. (7) and (9) are important system parameters. However, in order to limit the degrees of freedom of the setups, we only consider the case where that is, equal to the buoyancy associated to the reference potential temperature u ref ; see Eq. (4). Furthermore, the finite molecular viscosity n and scalar diffusivity k of the air govern the diffusive processes at the smallest length scales of the atmosphere. However, the treatment of their exact values is considered to be a method specific consideration and it is therefore discussed in section 3.

b. Physical scales and dimensionless groups
In our a priori analysis of the system, we distinguish between relevant scales for the daytime (0 , t , T/2) CBL and the nighttime (T/2 , t , T) SBL. We assume that B 0 . 0 . B 1 and kB 0 k ) kB 1 k, as is hinted in for the surface buoyancy for the CBL and SBL period, respectively, Noting that in general, the assumption of vanishing sensible heat flux is only accurate for the very stable boundary layer regime (Howell and Sun 1999). The assumption is particularly unsuitable for the clear-sky convection-driven ABL at mid day (i.e., when Q* 5 B 0 ). Therefore, we will motivate a more relevant alternative for the daytime buoyancy scale: it is well known that the convective boundary layer is characterized by a wellmixed region with a characteristic buoyancy scale. When we consider the convective limit, with no coupling to the underlying soil (L / 0), then this buoyancy scale can be evaluated using the total available buoyancy provided by Q* to uniformly ''heat'' the ABL's initial buoyancy profile. If we assume a well-mixed layer in the ABL, and the persistence of the initial stratification aloft, then a CBL with height L c can be associated with a convective buoyancy scale b c,Q* , The total time-integrated buoyancy flux during a single day is limited and therefore an upper bound for the values of L c and b c,Q* exist. We equate the red-shaded areas in Figs. 2 and 3 (encroachment principle); Combining the equations results in an expression for the CBL height scale L c and the corresponding buoyancy scale b c,Q* , For the SBL, the typical boundary layer height is a result from the internal dynamics: for example, a very SBL (VSBL) is typically shallower compared to a weakly SBL (WSBL). Hence, we cannot a priori formulate a relevant expression for the SBL height. However, if we follow a few assumptions, we can find a reference length scale for the depth of the SBL, denoted L s . d The SBL starts to cool from the surface under a wellmixed layer with a buoyancy that is O ( b c,Q* ) d The surface buoyancy of the SBL will decrease over the course of the night to be O (b s,L ).
with b diurnal a diurnal-buoyancy-range scale based on the first two assumption in this list, Following the same reasoning as with Eq. (13), we assume the blue-shaded areas of Figs. 2 and 3 to be equal, and therefore we write Apart from the geostrophic velocity scale U geo , a relevant velocity-fluctuation scale due to convection can be expressed according to Deardorff (1970): Similarly, we express a reference velocity scale for the SBL (cf. van de Wiel et al. 2012;van Hooijdonk et al. 2015), The present system is governed by the set of seven parameters {B 0 , B 1 , T, N, f, L, U geo }. With two base units (i.e., for length and time), we can identify sets of five independent dimensionless groups that govern the system dynamics (Buckingham 1914). An example of such a set is Here we have chosen to nondimensionalize B 1 , N, f, and L with the daytime forcing parameters B 0 and T, and U geo with the Deardorff velocity scale U c . They are formulated such that the values of these groups are typically larger than unity (see section 3). The reference potential temperature u ref , the acceleration due to gravity g and properties of the air (r and C p ) are not additional variables that determine the system dynamics. This is due to the fact that the heat fluxes are expressed in terms of buoyancy fluxes, the case setup and corresponding assumptions. Note that the set of groups in Eqs. (21) is based on an arbitrary choice out of infinitely many (equivalent) possibilities and any set of five orthogonal groups can be re-expressed in terms of those presented in Eqs. (21) (Vaschy 1892;Buckingham 1914).

Numerical case studies
This section presents the setups for two example studies on the diurnal cycle using both an SCM and an LES approach. Additionally, relevant dimensionless groups for a DNS study are discussed.
a. An example single-column-model study For this exemplifying study, we are inspired by the work of van der Linden et al. (2017), who developed a multiyear climatology of the nocturnal boundary layer at Cabauw, from the perspective of large-scale pressuregradient forcing. As discussed in depth, our model setup does not aim to capture the conditions to accurately model the ABL at Cabauw. Yet, it is interesting to see what wellknown features of the diurnal cycle are realistically reproduced with by the model. As such, we run a suite of model cases where we keep the values for P 1 , P 2 , P 3 , and P 4 fixed [see Eqs. (21)], and vary the parameter that controls the pressure-gradient forcing P 5 .
Associated with the pressure-gradient forcing is the geostrophic wind U geo ; see Eq. (1). It is well known that with increasing U geo , the wind shear also increases, which affects the mixing within the ABL. The shear owes its existence to the combination of a mean wind and the fact that velocities vanish at the surface. Since our single-column model does not resolve the turbulent flow near the surface, we rely on a closure that parameterizes the effect of a rough surface as a drag force on the flow. This warrants the introduction of a so-called roughness length for momentum z 0,m to our setup. We therefore add a sixth dimensionless group that governs the model system, Furthermore, to compute the buoyancy flux B, the soil flux G needs to be evaluated using the surface buoyancy b surf . In our model b surf is found by assuming a logarithmic variation of the buoyancy with height on the subgrid scale within the lowest two grid cells. A profile may be fitted and the buoyancy can be evaluated at the surface (see the appendix for details). This introduces a roughness length for scalars z 0,h , and a corresponding dimensionless group, It is commonly accepted that P 7 # 1, and there is considerable debate in the literature what it value should be. For simplicity, we choose a value of P 7 5 1 in this work.
Noting that z 0,h is related to b surf and that this buoyancy value should be interpreted as an effective surface buoyancy with respect to Eq. (7). We choose values for the six groups of Eqs. (21) and (22) according to the physical parameters listed in Table 1: and the following values for the dimensionless groups as presented in Eqs. (21) and (22): A SCM solves an evolution equation for the relevant atmospheric vertical profiles according to parameterized descriptions for turbulent transport (Betts 1986). We use the single-column model that is described in van Hooft et al. (2018a), based on the Basilisk toolbox that is available at www.basilisk.fr (Popinet 2015). The descriptions for the turbulent transport are based on simple first-order, local, stability-dependent K-diffusion closures (Louis et al. 1982;Troen and Mahrt 1986;Holtslag and Boville 1993). The mesh-element sizes are varied adaptively, based on the evolution of the numerical solution. For clarity of presentation, this section does not address the details of our model. Therefore, a more elaborated description is presented in the aforementioned work of van Hooft et al. (2018a) and the appendix. The implementation are documented and freely available. Links to their corresponding online locations are given in the appendix.
b. An example large-eddy simulation study The results of the SCM study are supplemented with a 3D turbulence-resolving study that concerns a similar setup as the SCM study, but using a small value for the wind-forcing parameter, P 5 5 0.5. This value is chosen small in order to limit the effects of shear on the ABL dynamics, yet large enough to ensure that the closures for the surface transport (i.e., a logarithmic ''law of the wall'') remains valid (Sullivan and Patton 2011;Ansorge 2019). The values of the other dimensionless groups are taken from the SCM study [see Eqs. (25)]. With an LES approach, one aims to resolve the most dominant (turbulent) flow structures in the atmosphere explicitly, assuming that 1) these govern the overall turbulent dynamics of the ABL, 2) these large turbulent structures decay into smaller structures via a predominantly downward cascade, and 3) these smaller-scale structures are universal and can be parameterized with sufficient accuracy. Such an approach may be preferred over RANS modeling as it yields 3D solution data that explicitly includes the internal variability of the ABL system and reduces the reliance on closures for turbulent transport. However, not the full range of turbulent motions are resolved and this requires an associated model for the subgrid-scale (SGS) turbulent transport. In this work, we parameterize the effects of the SGS turbulent motions to be locally diffusive with an effective eddy viscosity that is calculated based on the formulation of Vreman (2004). Furthermore, at the approach of the surface, the size of the dominant turbulent eddies decreases to such small (viscous and roughness) scales, that the effect of the underlying surface also needs to be parameterized similar to the SCM approach. The LES domain is a cube of size L 3 0 , with L 0 5 3L c and the lateral boundaries prescribe periodicity of the solution in the horizontal directions. The resulting convective boundary layer depth will grow to be '0.8L (see section 4). The aspect ratio of the horizontal periodicity and the vertical length scale (approximately 3.75:1) is in the ''danger zone,'' as there is the risk of large plumes interacting with themselves across the periodic boundaries (Schmidt and Schumann 1989;van Heerwaarden and Mellado 2016).  van Hooft et al. 2018b). The mesh resolution is an important parameter when assessing the fidelity of an LES result (Sullivan and Patton 2011;van Stratum and Stevens 2015). As such, the details for the LES setup, including the resulting mesh-size distribution, are presented in the appendix.

c. Toward a direct numerical simulation study
For the 3D LES study we have adopted a parameterized description for the drag at the surface and the mixing by SGS turbulence. It is well known that the assumptions for an LES are not generally valid and that SGS models may fail to describe the turbulent mixing with sufficient accuracy for all atmospheric conditions (Mellado et al. 2018;Ansorge 2019). Alternatively, one may aim to solve the set of equations for fluid motion directly, without adopting a closure for the turbulent transport (Moin and Mahesh 1998). This entails that the flow needs to be resolved down to the length scales of the turbulent cascade where gradients are sufficiently large such that the molecular diffusion becomes an effective dissipation mechanism for the second-order moments of momentum and buoyancy. This warrants the introduction of the fluid's viscosity n and thermal diffusivity k, which introduces a dimensionless group that is known as the molecular Prandtl number: From the seminal work of Kolmogorov (1941), we learned that the smallest relevant flow-structure scales for the inertial subrange of turbulence are dictated by the dissipation rate of turbulent kinetic energy « and the fluid's molecular viscosity n. It is known as the Kolmogorov, viscous, or dissipation length scale h and may be written based on dimensional analysis (Buckingham 1914), To assess the feasibility of performing a DNS of the present diurnal cycle setup, in the limiting case where there is a total absence of wind forcing (i.e., U geo 5 0), we follow the works of Garcia and Mellado (2014) and van Heerwaarden and Mellado (2016) and approximate the maximum value of « as a fraction of order one of B 0 . We can then define a dimensionless group as the four-thirds power of the scale separation that is known as the Reynolds number (Re), Typical atmospheric conditions are characterized by U c 5 O (1) m s 21 , L c 5 O (1) km, and furthermore, n ' 1.5 3 10 25 m 2 s 21 meaning that Re 5 O (10 8 ) and the minimum resolution ½D min 5 O (h) would then need to be Â D min 5 O (10 26 L c ) Ã on the order of millimeters. This is not feasible for 3D turbulence resolving studies and may remain elusive for a considerable period in the future (Bou-Zeid 2015). Therefore, in virtually all studies that employ DNS for ABL flows, an offer is made with respect to the Reynolds number, assuming that for sufficiently large values of Re, certain relevant (scaled) statistics become insensitive to its exact value, a concept known as ''Reynolds-number similarity.'' Similar to the LES approach, it is assumed that the full depth of the inertial subrange of isotropic turbulence does not need to be resolved to obtain accuracy for the lower-order solution statistics. However, with the available resources for this preliminary study, we cannot perform a DNS with a reasonable value for Re of the full diurnal cycle and hence, this work does not present any results obtained with the DNS technique. Further, heat fluxes in such a simplified representation of the ABL do not scale inviscidly (cf. van Heerwaarden and Mellado 2016;van Hooijdonk et al. 2018), that is, they have to be normalized with account for changes in Re. This requires a thorough analysis of the viscid scaling behavior of P 1 2 P 5 .

a. Single-column model results
Since the SCM approach employs semiempirical descriptions of the turbulent transport, the obtained evolution of the vertical profiles are subject to the particular choices regarding the used closures. As such, the results presented in this section are only intended to provide a qualitative overview of the modeled dynamics. Figure 4 presents the results from the simulation run with P 5 5 3.5 and provides a prototypical example of the evolution of the vertical profiles with a time interval Dt 5 T/12 (corresponding to a 2-hourly interval). We observe that the for this case, the depth of the well-mixed layer grows during the daytime to the order of L c with a characteristic buoyancy that is a fraction of order one of b c,Q* . During the night, an SBL forms that has a depth on the order of L s [cf. Eq. (18)]. As the SBL is shear driven, the shapes of the nighttime profiles are more sensitively responding to the value of P 5 compared to their daytime counterparts (not shown). Figure 5 shows the partitioning of the available buoyancy Q* between the buoyancy flux B and the soil flux G for two runs with different pressure-gradient forcing (i.e., P 5 5 {2, 5}). It can be observed that for the daytime results, there are hardly any differences, as thermodynamics are dominated by convection. At the end of the day time (t ' T/2), the buoyancy flux changes from positive to negative values before Q* does so (cf. van der Linden et al. 2017). At night however (i.e., B , 0), the weak wind speeds cause the buoyancy flux to be close to zero and the radiative buoyancy/energy loss is compensated by the soil flux G, corresponding to a VSBL. Alternatively, for the strong-wind case (i.e., P 5 5 5), the soil flux remains small in the SBL and the radiative loss of energy is mostly compensated by the buoyancy flux B, corresponding to a WSBL [see, e.g., observational data presented in van de Wiel et al. (2003)]. Remarkably, the simple setup appears sufficient to retrieve the dynamical difference between the VSBL and the WSBL due to variations in the pressure-gradient forcing. Inspired by this result, we explore the transition between the VSBL and WSBL as a function of P 5 and plot the difference of buoyancy at the surface and at z 5 L c /20, denoted as Db. This inversion is averaged over hour 19 of simulation (i.e., in the middle of the night) and the results are plotted in Fig. 6a. It appears that for the low pressure-gradient FIG. 4. The solution profiles obtained for the single-column model run with P 5 5 3.5. The plots display the profiles of (a) the daytime the buoyancy b and (b) the horizontal wind magnitude U. (c),(d) As in (a) and (b), but at nighttime. Notice the different scaling of the axis for daytime and nighttime data. For reference; all panels show the corresponding profile at t 5 T/2 (red line).
forcing the inversion is of the order one of the anticipated diurnal buoyancy range b diurnal . In literature that concern observational analysis, the large-scale pressuregradient forcing is commonly not available and the inversion strength is often plotted against a wind speed within the boundary layer. As such, we average the wind speed at z 5 L c /20 [ 50 m over the same period and the results are plotted in Fig. 6b. Here we see that the transition between the WSBL and VSBL appears much sharper. This is due to the nonlinear response of the wind speeds within the SBL to U geo . For a more in-depth analysis of this aspect, see Baas et al. (2019).
Finally, To point out a limitation of the model scenario with respect to its ability to describe reality, we focus on the evolution of the buoyancy near the surface [b(z, t) 5 b(L c /50, t)] in Fig. 7 for the run with a value of P 5 5 4. For a comparison, we plot the evolution of the absolute temperature at z 5 20 m for five consecutive days at the CEASAR observational site in Cabauw, The Netherlands. Data are plotted for 19-23 August 2018. These days are associated with a heat wave and a drought in the Netherlands. These data therefore provide a reference for diurnal cycles with a limited influence of the moist dynamics. Neglecting the actual diurnal temperature range, we see that the overall shape of the observations is captured to some degree in the model. However, during the transition from the CBL to the SBL, the cooling rate predicted by the model is (comparatively) too large. This may be due to the fact that the used closures fail to represent the mixing of decaying turbulence and/or the fact that the soil heat flux description in our model is particularly unsuitable for this period of the day (cf. Fig. 1). To assess this discrepancy, we have run the model with a full (i.e., layered) soil model such that heat-storage effects are explicitly taken into account. An alternative model where we have used the so-called ''enhanced mixing'' formulation under stable conditions is also run (Holtslag and Boville 1993). The cooling rate at t ' T/2 is smaller than with the present L-based description of the soil flux G [cf. Eq. (7)]. This is due to FIG. 5. The diagnosed buoyancy flux B, the soil heat flux G, and the prescribed net radiation Q* obtained from two runs, each with different P 5 value (P 5 5 {2, 5}). The differences due to the wind forcing are most notable in the stable boundary layer regime (i.e., t . T/2). After the occurrence of a buoyancy flux minimum, the sensible heat flux vanishes for the case with the weaker wind forcing, whereas it remains finite for the case with the stronger forcing case.
FIG. 6. The buoyancy inversion Db between the surface and at z 5 L c /20 averaged over simulation hour 19 obtained from runs with different pressure gradient forcings (i.e., P 5 ). The value of Db is plotted against (a) the pressure forcing parameter P 5 and (b) the hour-19 average of the horizontal wind speed at z 5 L c /20 [U z(Lc/20) ]. The dimensionless numbers on the lower and left-hand-side x and y axis are converted to their meteorological equivalents using the values in Table 1 and are displayed on the corresponding upper and right-hand-side axes (gray).

DECEMBER 2019 V A N H O O F T E T A L .
significant effects of soil-heat storage: the temperature of the upper soil is higher than that of, for example, the deep soil T d in Eq. (6) (see Fig. 7c). Therefore, heat is released into the atmosphere. It can also be observed that the model using the enhanced mixing formulation does indeed predict a slower cooling rate. However, even with this mixing, the rapid cooling during the transition remains, which stresses the importance of an accurate representation of ground heat storage. As such, the present lumped parameter model may be a poor representative of reality during the transition period itself (cf. Fig. 1).

b. Large-eddy simulation results
In this section we give an overview of the evolution of the resolved turbulent structures in the ABL during the diurnal cycle with the present framework, containing coupled surface boundary conditions. Due to our low value for P 5 , turbulence is weak during the night and the main actors in the SEB are Q* and G (see Fig. 5). Correspondingly, there are no appreciable turbulence levels in the simulated SBL and the authors argue that in the absence of large eddies, the assumptions that underlie the LES approach may not hold for this period. Therefore, we focus our analysis on the growth and decay of the CBL and we only briefly consider the day-night transition. The SBL itself is left for a more in-depth future study. As we have not prescribed the evolution of buoyancy at the surface, the internal variability of the atmospheric dynamics give rise to a heterogeneous surface buoyancy structure. We plot horizontal slices of the surface buoyancy for the part of the day where the CBL is growing in Fig. 8. It appears that initially, at t/T 5 1/48, the surface buoyancy structures are organized in elongated ''streets,'' whereas at later stages, the footprint of horizontally isotropic convection cells are observed. This can be explained by the fact that in the early morning, convection is still weak compared to the shear in the boundary layer. It is known that shear dominated convection tends to organize in elongated surface ''streaks'' (Adrian 2007) or convective ''rolls'' (Coleman et al. 1994;De Roode et al. 2004;Conzemius and Fedorovich 2008;Gibbs and Fedorovich 2014). As the buoyancy flux and the boundary layer height increase, this balance shifts to a convection dominated regime. Furthermore, with an increase of the CBL depth over time, the size of the convection cells also grows and this is reflected by the structures in the surface buoyancy. Apart from the largest-cell structures, smaller filaments can be observed. Within the CBL, thermal plumes rise and this process may be visualized by plotting 3D renderings of an isosurface of the buoyancy colored with the local vertical velocity. Figure 9 presents such rendering at two stages for the growth of the CBL. The isosurface of the buoyancy is chosen such that it is representative of cores of the thermal plumes and a layer with an identical buoyancy value is found aloft in the inversion layer. It can be observed that at t/T 5 1/24 there exist many convective plumes within the domain that cause a wave pattern in the overlaying inversion layer. At a later stage (i.e., t/T 5 6/24) the thermal plumes have grown in size. The snapshot in Fig. 8b) also depicts the process of entrainment (arrow) as we see the isosurface of the buoyancy in the inversion layer being entrained downward to be within the CBL. Finally, 1700 consecutive snapshots of the turbulent structures in the solution field are combined into a movie. It is made available via Vimeo (https://vimeo.com/292329175).  7. (a) The evolution of the temperature at z 5 20 m for five consecutive clear-sky and relatively dry days with a pronounced diurnal pattern, at Cabauw, the Netherlands and (b) the modeled buoyancy at z 5 L c /50 [ 20 m obtained from three different models (see text section 3). It appears that at the onset of the stable boundary layer (t ' T/2), only the model that includes heat-storage effects (yellow line) is able to predict a realistic cooling rate. (c) The corresponding soil-buoyancy profiles are displayed for various times during the day-to-night transition. These profiles reveal that after the transition, the soil is releasing heat into the atmosphere.
An important characteristic of convective boundary layers is that it may display self-similar behavior over a range of CBL heights L c and convective velocity scales U c (Jonker et al. 2013). For the well-mixed-layer height growing into a linear stratification L c , we follow the definition of van Heerwaarden and Mellado (2016), Based on this length scale, we define the instantaneous Deardorff convective velocity scale, The temporal evolution of these scales is plotted in Fig. 10a and are used to rescale the vertical profiles of the vertical velocity variance hw 2 i in Fig. 10b. Here we see that for the largest part of the daytime a self-similar profile is diagnosed that is described with appreciable accuracy by the well-known functional shape as listed in Troen and Mahrt (1986). However, the profiles from the early morning and the later afternoon show clear deviations form the self-similar profile. The shallow early-morning CBL may suffer from the limited grid resolution and may also be affected by the relative importance of shear (Conzemius and Fedorovich 2008). Furthermore, during the late-afternoon period, the heat flux vanishes and the convective turbulent structures decay. Hence, it is not expected that the CBL remains self similar nor that Deardorff scale U c is relevant for that period. To obtain an estimate for the growth and decay of the vertical motions, the height-integrated vertical velocity variance e w 5 Ð ztop z 0,h hw 2 i dz , averaged for a period Dt 5 T/48 is presented in Fig. 10c. In the simulation, the decay of e w during the day-to-night transition appears to be linear with time. This does not follow the conclusions from any of the more idealized studies on the decay of convective turbulence (Nieuwstadt and Brost 1986;Sorbjan 1997;van Heerwaarden and Mellado 2016), where several negative scaling exponents with time are proposed to describe the decay of convection. This aspect is not further discussed herein.
Finally, we briefly discuss the surface structure for the buoyancy and surface heat flux of the present setup. FIG. 9. Isosurface renderings of the buoyancy field for (a) t/T 5 1/24 and (b) t/T 5 6/24, with b iso /b c,Q * 5 f0:125, 0:64g, respectively. The values for the iso surfaces are (subjectively) chosen such that they depict the upward convective motions and the surfaces are colored with the local vertical velocity from red (upward) to green (w 5 0) to blue (downward). The black arrow in (b) annotates the process of boundary layer entrainment.  (30) and (b) the vertical z profiles of the velocity variance hw 2 i, plotted for different times and normalized with the external scales L c and U c , respectively. The dashed black line indicates a fit (by eye) of a functional relation for the self-similar profile based on Troen and Mahrt (1986). (c) The growth and the subsequent decay of the resolved turbulent motions are quantified by the total height-integrated vertical velocity variance e w .
We note that there is quite some freedom on how to implement the surface boundary conditions when using a SEB equation. In this work, B is computed as a residual in the SEB, B 5 Q* 2 G, were G is evaluated via the buoyancy at the surface b surf . In turn, the surface buoyancy is computed from fitting a log-linear profile to the local buoyancy field in the lowest two grid points and evaluating the corresponding buoyancy at a height of z 5 z 0,h (see the appendix for details). Even though assuming a logarithmic profile seems to be a robust choice, the main disadvantage of this approach is that the exact value of z 0,h becomes an important parameter to obtain b surf . Alternatively, a method taken from the TESSEL/HTESSEL land surface model (Balsamo et al. 2009;Maronga and Gehrke-Scharf 2018) has successfully been employed in the LES study by Maronga and Bosveld (2017). Here the surface temperature is computed as the temperature that closes a (more extensive) SEB equation, where the individual terms are a function of T surf . Remarkable, an inconsistency arises between the two methods: with the present model formulation, the correlation between the instantaneous fluctuations in the surface buoyancy field and the surface buoyancy flux (i.e., b 0 surf and B 0 ) are linked via L. This is because the regions with a warmer surface correspond to a heat flux toward the soil G, and vice versa. Hence, Alternatively, the TESSEL/HTESSEL formulation computes the sensible heat flux H according to a resistance law (i.e., bulk transfer model), where u surf and u 1 are the potential temperature at the surface and the first model level, respectively, and r is the aerodynamic resistance which is computed based on a Monin-Obukhov-style closure. In practice, Eq. (32) causes a positive correlation between the fluctuations in H and in u surf , as warmer surfaces transport heat more effectively toward the atmosphere (van Tiggelen 2018). Although this seems reasonable, we argue that the formulation of Eq. (32) was designed for weather forecasting models where the resolved spatial scales are much larger than individual convection cells. Whereas the LES approach applies its closures for the variations within a convection cell, which are governed by different dynamics. Figure 11 shows an illustrative representation for the convergence of warm air near the surface toward the base of a convection cell. At the base, the surface temperature is at a local maximum and at some height above the surface, this is where the heat flux is indeed large. However, it is not obvious how the fluctuations, albeit small, in temperature and sensible heat flux at the surface relate to each other. With the advent of convection permitting weather models (Prein et al. 2015) and the usage of LES for weather prediction (e.g., Heinze et al. 2017), a proper description is needed. This may be based on the results from a dedicated observational campaign and/or direct numerical simulations of the near-surface heat transport that includes the heat exchange with the soil.

Conclusions
In this work, we have presented a conceptual model for the diurnal cycle of the dry atmospheric boundary layer (ABL). The corresponding forcing parameters can be described with only five dimensionless groups and this provides a convenient framework to study cause-and-effect relations for the effects of the various forcing mechanisms on the ABL dynamics. Relevant scales for the boundary layer height and atmospheric buoyancy values are derived from the model. Furthermore, the model served as a setup for numerical simulations. This aspect is illustrated for the atmospheric single-column model and a large-eddy simulation techniques.
The results from the preliminary study that used a single-column model showed that the difference between the weakly stable and very stable boundary layer are dynamically retrieved for variations in the wind forcing. This is a key result of the present work and shows that interesting physics can be present in a model that only has a low degree of complexity and low number of input parameters. A large-eddy simulation for a diurnal cycle with weak-wind forcing showed how, for this case, convective structures grow and retain their statistically self-similar structure throughout the day. Furthermore, the decay of residual turbulence in the late-afternoon-to-evening transition was quantified. The integrated vertical velocity variance appeared to decay linearly with time.
Limitations of the present description were also identified. Due to the nonphysical nature of the lumped parameter L to describe the response of the soil to varying surface temperatures, so-called history effects are not modeled. These effects appear to play an important role in the day-night transition. Furthermore, for the three-dimensional simulations, the model predicts a heterogeneous surface buoyancy and surface buoyancy flux. However, by design, deviations from the means of these two quantities are perfectly anti correlated. This may not be accurate, but we argue it remains unclear how to accurately model the soil response in an LES.
We conclude that the present conceptual model provides additional insight to the classical (Stullian) view of the diurnal cycle (Stull 1988), as it enables to define expressions for relevant physical scales that characterize the diurnal cycle. The dimensionless description can be useful to generalize the dynamical behavior of the ABL. Finally, we foresee that more in-depth numerical studies, based on the present idealized forcings for the ABL, can help to better understand the interactions between the most prominent processes that govern the diurnal cycle. provided by the European Research Council via the consolidator grant (648666). We thank Fred Bosveld for the data from KNMI's CEASAR database. Furthermore, we acknowledge that we are indebted to Stéphane Popinet for the Basilisk code (www.basilisk.fr) and to those who contributed to GNU/Linux and/or the GNU-compiler collection (www.gnu.org).

Numerical Model Setups
In this section we discuss the details regarding the single-column model (SCM) and large-eddy simulation (LES) setups. The exact implementation of the described approaches are documented and made available online: The SCM solves an evolution equation for the vertical profiles of the horizontal velocity components u and y and the buoyancy b according to a parameterized description for turbulent mixing: where U geo is the geostrophic wind, f the Coriolis parameter [cf. Eq. (1)] and K the eddy diffusivity. The latter is parameterized according to with l the Blackadar length scale and V* a velocity scale due to shear and convection. For l we write with k 5 0.4 the von Kármán constant. For V* we depart from the work of van Hooft et al. (2018a) as we do not only include the effects of shear, but also model mixing due to convective motions. We follow the work of Troen and Mahrt (1986), where w c is the local convective velocity scale that is parameterized using where L c and U c are computed via Eqs. (29) and (30). The shear magnitude S is evaluated as and f(Ri) is the stability-dependent mixing function, that is, based on the gradient Richardson number (Ri), For the enhanced mixing experiment, we use the so-called ''long tail'' formulation (Louis et al. 1982), The drag at the bottom surface is based on a near-neutral logarithmic law of the wall, and we express the friction velocity u t as function of the bottom-gridcell-averaged value of u (labeled u 1 ): which follows directly from integrating the law of the wall over the lowest cell with size D ) z 0,m . A logarithmic law of the wall is also applied for the buoyancy profile in the lowest two grid cells. This results in an expression for b surf , with b 1 and b 2 the cell averaged values of b in the lowest two cells. The value of c results from an integration exercise, where we have assumed that the lowest two cells have a vertical size of D and D ) z 0,h , the roughness length for heat. The surface buoyancy flux B 5 Q*(t) 2 G can then be readily evaluated. The domain has a height of z top 5 3L c and the mesh-element sizes are adaptively varied down to a maximum resolution of D min 5 z top /512 [ 6 m, based on the estimated discretization error in the representation of the solution fields for u, y, and b (Popinet 2015; van Hooft et al. 2018a). The corresponding refinement-criteria values are z u,y 5 U geo /20 and z b 5 b diurnal /50. These values resulted from a convergence study and aim to strike an arbitrary balance between accuracy and the number of grid cells (van Hooft et al. 2018a).

THE LAYERED SOIL MODEL
For the SCM run with a layered soil model (see Fig. 7), the flux (G 5 G soil ) is described via a conduction layer on top of a soil: The layers have a buoyancy-equivalent-to-heatcapacity (rc y ) soil 5 1000 3 (rc p ) air , a diffusivity k soil 5 5 3 10 25 m 2 s 21 and the conductivity A 5 2 3 10 23 m s 21 .
The eddy diffusivity K is evaluated according to the subgrid-scale model of Vreman (2004) and at the bottom  (2017) and Popinet (2018). The spatial discretization follows the finite-volume methodology, employing cubic cells to mesh the cubic (3L c 3 3L c 3 3L c ) domain using an octree grid structure. A damping layer in active in the top half of the domain. The octree grid allows gridcell sizes to vary by factors of 2 (Popinet 2015). To resolve the largest eddies, those that dominate the dynamics, the tree grid is locally refined and coarsened adaptive to the prevailing situation. Each time step, a multiresolution analysis is performed to estimate the local deviation from a spatially linearly varying test function. This deviation is computed for all velocity components x u,y,w and the buoyancy field x b . Grid refinement is carried out when the x value exceeds the corresponding refinement criterion (x . z), with a maximum resolution corresponding to a 256 3 -cell equidistant grid. Cells may be coarsened when x , (2/3)z, with the additional requirement that the coarse cell ''fits'' the hierarchical grid structure (van Hooft et al. 2018b).
To properly take into account the time-varying nature of the daytime convection, the refinement criteria are dynamically evaluated as fractions of the Deardorff velocity and buoyancy fluctuation scale (for B . 0), z u,y,w 5 max U c c u , z min,u ! , Here c u and c b are dimensionless constants that may be tuned to balance accuracy (via grid resolution) versus computational speed. We have chosen c u 5 2 and c b 5 1 and furthermore, z min,u 5 U c /10 and z min,b 5 (B 2 0 /L c ) 1/3 /5.
The authors of this work realize that hinting at the accuracy of a simulation via z is unconventional compared to listing a fixed mesh-element size (see van Hooft et al. 2018a). Therefore, Fig. A1 displays the horizontally averaged effective resolution, presented via the corresponding number of equidistant grid cells. Furthermore, Fig. A2 presents a horizontal and vertical slice of the grid at midday (i.e., t/T 5 6/24). It can be observed that the mesh elements are dynamically focused on the surface layer, the entrainment and inversion layer and toward the decay of convective turbulence. Furthermore, At the late stage of the simulation, the internal waves in the free troposphere are the most dynamically active process. The grid slices in Fig. A2 reveal that within the CBL, a high resolution is employed to resolve the cores of the thermal plumes.