# Generative convective parametrization of dry atmospheric boundary layer

Florian Heyder  Institut für Thermo- und Fluiddynamik, Technische Universität
Ilmenau, Postfach 100565, D-98684 Ilmenau, Germany  Juan Pedro Mellado
Meteorologisches Institut, Universität Hamburg, Bundesstraße 55, D-20146
Hamburg, Germany  Jörg Schumacher  Institut für Thermo- und Fluiddynamik,
Technische Universität Ilmenau, Postfach 100565, D-98684 Ilmenau, Germany
Tandon School of Engineering, New York University, New York City, NY 11201,
USA

###### Abstract

Turbulence parametrizations will remain a necessary building block in
kilometer-scale Earth system models. In convective boundary layers, where the
mean vertical gradients of conserved properties such as potential temperature
and moisture are approximately zero, the standard ansatz which relates
turbulent fluxes to mean vertical gradients via an eddy diffusivity has to be
extended by mass flux parametrizations for the typically asymmetric up- and
downdrafts in the atmospheric boundary layer. In this work, we present a
parametrization for a dry convective boundary layer based on a generative
adversarial network. The model incorporates the physics of self-similar layer
growth following from the classical mixed layer theory by Deardorff. This
enhances the training data base of the generative machine learning algorithm
and thus significantly improves the predicted statistics of the synthetically
generated turbulence fields at different heights inside the boundary layer.
The algorithm training is based on fully three-dimensional direct numerical
simulation data. Differently to stochastic parametrizations, our model is able
to predict the highly non-Gaussian transient statistics of buoyancy
fluctuations, vertical velocity, and buoyancy flux at different heights thus
also capturing the fastest thermals penetrating into the stabilized top
region. The results of our generative algorithm agree with standard two-
equation or multi-plume stochastic mass-flux schemes. The present
parametrization provides additionally the granule-type horizontal organization
of the turbulent convection which cannot be obtained in any of the other model
closures. Our work paves the way to efficient data-driven convective
parametrizations in other natural flows, such as moist convection, upper ocean
mixing, or convection in stellar interiors.

##  I Introduction

Convective phenomena occur at various spatial and temporal scales in the
atmosphere, ranging from meter-scale thermals at the surface Lenschow and
Stephens (1980) via individual clouds or cloud layers Bodenschatz et al.
(2010); Mellado (2017) and mesoscale convection Atkinson and Wu Zhang (1996);
Bony et al. (2020) to global-scale patterns, such as the monthly timescale
Madden-Julian oscillation in the tropics Madden and Julian (1972); Miura et
al. (2007). Capturing this spatial heterogeneity, the transition between
different convective regimes, and the resulting fluxes remains a challenge and
requires turbulence parametrizations of the dynamics at the scales which
remain unresolved in the numerical models. Mass-flux parametrizations are one
specific way of representing unresolved convective sub-grid scale processes in
atmospheric turbulence Ertel (1942); Deardorff (1972); Arakawa (1974), see
also LeMone et al. (2019); they belong to the class of counter-gradient models
which become necessary in the absence of a mean gradient of conserved
properties Tiedtke (1989); Siebesma and Cuijpers (1995); Želi et al. (2021).
These parametrizations of the dynamics in atmospheric boundary layers (ABL)
capture the effective vertical transport of heat, moisture, and momentum in
the lower atmosphere in a set of reduced equations that are solved together
with the transport equations for energy and total water content Wyngaard
(2010). By incorporating mass-flux schemes, kilometer-scale global circulation
models can couple the complex interactions between large-scale dynamics and
smaller-scale convective processes at mesoscales, essential for improved
predictions of shallow cloud patterns Klingebiel et al. (2021); Vogel et al.
(2022) and resulting rainfall distributions Ahlgrimm and Forbes (2014); Becker
et al. (2021) or short-term climate change prognoses Hohenegger et al. (2023).
While these physical parametrization models have proven to be valuable tools,
they have some limitations resulting from a set of various free parameters
that have to be adjusted to control the physical processes Siebesma et al.
(2007). Furthermore, they do not provide any information on the horizontal
organisation of the updrafts which determines the cloud bases that form on the
top of the dry layer.

In the past years, machine learning (ML) methods have shown promise as an
alternative and data-driven approach to model convective processes even though
their generalization capabilities remain often limited Brenowitz and
Bretherton (2018); Reichstein et al. (2019); Kashinath et al. (2021). These
algorithms have been applied to classify convection patterns Fonda et al.
(2019); Rasp (2020); they served as subgrid-scale models in large-eddy
simulations Beck and Kurz (2021); Yuval et al. (2021); Cheng et al. (2022);
Guan et al. (2023), learned the nudging tendencies for the correction of
physical parametrizations Bretherton et al. (2022), or were used as reduced-
order models in the form of recurrent neural networks to reproduce low-order
turbulence statistics Heyder and Schumacher (2021); Heyder et al. (2022). Deep
generative ML algorithms, e.g., in the form of Generative Adversarial Networks
(GANs), create synthetic data which correspond to convection phenomena, such
as cloud cover Dai and Hemri (2021), vertical cloud structure Leinonen et al.
(2019), or precipitation distributions Price and Rasp (2022). Shamekh and
Gentine used a variational autoencoder, a further generative machine learning
model, to approximate the mean vertical profiles of the buoyancy flux due to
wind shear and convection Shamekh and Gentine (2023). This specific class of
machine learning methods opens thus doors to access the horizontal
organisation of the turbulent convection and sets the stage and motivation for
the present investigation.

In this work, we present a data-driven mass-flux parametrization of the
buoyancy flux in the form of a GAN Goodfellow et al. (2014, 2016). We
therefore consider a paradigm for the dry atmospheric boundary layer, a shear-
free flux-driven convective boundary layer. In a high-fidelity three-
dimensional direct numerical simulation (DNS), we fully resolve the transient
growth of the ABL Mellado and Ansorge (2012); Garcia and Mellado (2014). The
DNS snaphots provide training data for our semi-supervised ML method Brunton
et al. (2020). The GAN architecture consists of two competing deep neural
networks; (1) the generator is here a U-shaped convolutional neural network,
in short U-Net Ronneberger et al. (2015), (2) the discriminator is a deep
neural network. The generator models the transient growth of the ABL by
inclusion of the physical laws of self-similar growth of the layer Deardorff
(1972). In this way, we augment the training data base and thus enhance the
prediction capability of the ML method significantly. The GAN is set up for
different heights in the ABL, it provides the complete non-Gaussian statistics
of the vertical velocity, buoyancy fluctuations, and buoyancy flux at
reference heights in the sub-grid scale range. Our model is also capable to
generate synthetic data with the right area fraction of updraft regions. It is
found to agree with a mass-flux parametrization in the form of a set of
reduced equations with adjustable parameters. The latter is for example
applied for the atmospheric boundary layer in the Integrated Forecasting
System, a global circulation model of European Centre for Medium-Range Weather
Forecasts Siebesma et al. (2007); Witte et al. (2022).

Figure 1: Self-similar temporal growth of the atmospheric boundary layer.
(A-C) Contours of the vertical velocity fluctuations $w^{\prime}$ viewed from
the top onto the whole horizontal domain for three different times. (D-F)
Corresponding turbulent kinetic energy, ${\rm
TKE}=(u^{\prime\,2}+v^{\prime\,2}+w^{\prime\,2})/2$, in a vertical contour
plot. The horizontal extension of the cut is indicated in (A-C). (G) Plane
averaged vertical profiles of the buoyancy versus height at three time
instants from panels (A-C).

##  II Convective boundary layer

We solve the three-dimensional Boussinesq equations which couple the two
turbulent fields, the velocity field ${\bf v}=(u,v,w)$ and the buoyancy field
$b$, in DNS that resolve all physically relevant scales down to the
dissipation scales. In an atmospheric boundary layer without phase changes,
the buoyancy is directly proportional to the virtual potential temperature
field $\theta_{\mathrm{v}}$ and given by

$b({\bf x},t)=g\frac{\theta_{\mathrm{v}}({\bf
x},t)-\theta_{\mathrm{v,0}}}{\theta_{\mathrm{v,0}}}\,,$ (1)

with the acceleration due to gravity $g$ and a reference temperature
$\theta_{\mathrm{v,0}}$. Convection proceeds in a transiently and self-
similarly growing convective boundary layer (CBL) with an encroachment height
$h(t)$ that is to a good approximation given by Deardorff (1972); Garcia and
Mellado (2014)

$h(t)=L_{0}\left[2(1+Re_{0}^{-1})N_{0}(t-t_{0})\right]^{1/2}\,.$ (2)

This height contains the Brunt-Väisälä frequency $N_{0}$, the reference time
$t_{0}$, the Reynolds number $Re_{0}$, and the reference Ozimidov-type scale
$L_{0}=(B_{0}/N_{0}^{3})^{1/2}$. It is driven by prescribed buoyancy fluxes at
the bottom and top of the layer. The prescribed constant flux at the bottom is
$B_{\mathrm{v,0}}$. Consequently, the convective Rayleigh number
$Ra_{c}(t)=B_{0}h(t)^{4}/(\nu\kappa^{2})$ grows with time and reaches at an
intermediate time a value of $3\times 10^{8}$ in the present case. Here, $\nu$
and $\kappa$ are the kinematic viscosity and buoyancy diffusivity,
respectively. The Prandtl number is $Pr=\nu/\kappa=1$. Further details on the
model equations, the essential dimensionless parameters, and the numerical
simulation technique are provided in the  SI Appendix.

The buoyancy flux is the central transport measure. It couples buoyancy and
vertical velocity fluctuations, i.e., the deviations from the respective
vertical mean profiles of both fields. In an eddy-diffusivity mass-flux
approach, it is modeled as follows,

$\rho_{0}\langle
w^{\prime}b^{\prime}(z)\rangle_{x,y}=-\rho_{0}\kappa_{t}\frac{\partial\langle
b\rangle_{x,y}}{\partial z}+M_{u}(\langle b(z)\rangle_{u}-\langle
b(z)\rangle_{x,y})\,,$ (3)

where $\rho_{0}$ is the constant reference density in the convective boundary
layer and $M_{u}$ the mass flux parameter. Here, $\langle\cdot\rangle_{x,y}$
are horizontal plane averages. The first term in (3) is the standard
Boussinesq term for a turbulent stress which contains the turbulent
diffusivity $\kappa_{t}$. The second term on the right hand side stands for
the mass-flux parametrization which has to be included when the mean buoyancy
gradient vanishes for the well-mixed layer, as being the case here.

##  III Up -and downdraft regions

Figure 1 illustrates the temporal growth of the convective boundary layer.
While panels (A–C) provide a view from the top for the vertical velocity
fluctuations with an aggregation of the convection cells, panels (D–F) display
the turbulent kinetic energy in a vertical cut illustrating the transient
growth process. The remaining panel (G) of Fig. 1 demonstrates that this
growth follows a self-similar process and a rescaling of lengths with $h(t)$
and buoyancy with $N_{0}^{2}h(t)$ collapses the vertical mean buoyancy
profiles at different times, $\langle b(z,t)\rangle_{x,y}$. It is exactly this
self-similar behaviour given by (2) for the growth of the mesoscale patterns
in vertical and horizontal directions which will be used to augment the
training data base for the generative network.

Figure 2: Statistics of the up- and down-welling fluid motion in the
convective boundary layer (top $95\%$ of $w$-values). (A) Probability density
functions (PDFs) of the area $A_{u}$ of the up-welling fluid motion normalized
by the total horizontal area $A_{\rm total}$ as a function of $h/L_{0}$ (see
color legend in (B)). The two insets display a cutout of the updrafts in a
fixed horizontal plane $z=0.5$ at different times. (B) PDF of the area $A_{d}$
of the down-welling fluid motion normalized by the total horizontal area
$A_{\rm total}$ as a function of time. (C) Horizontal cutout which explains
the decomposition of $A_{\rm total}$ into regions of strongest updrafts
$A_{u}$, strongest downdrafts $A_{d}$, and the intermediate motions $A_{i}$.
Note that the union of the regions of the strongest downdrafts and
intermediate motions yields the environment $A_{e}$ for the mass-flux
formulation.

The prefactor $M_{u}$ in the second term of (3) incorporates the specifics of
the upward transport. As detailed in SI Appendix, the unknown mass-flux
parameter $M_{u}$ is estimated as a product of the updraft area fraction,
$a_{u}$, and the excess of the mean vertical velocity in the updraft regions,
$\langle w\rangle_{u}$, over the mean upward motion, i.e.,

$M_{u}\approx\rho_{0}a_{u}(\langle w\rangle_{u}-\langle w\rangle_{x,y})\,.$
(4)

The horizontal cross section of the layer $A_{\rm total}$ is decomposed into
disjoint strongest updraft (u) and downdraft (d) regions, as well as the
remaining intermediate motions region (i), i.e., $A_{\rm total}=A_{u}\cup
A_{d}\cup A_{i}$. We have thus further refined the environment region, which
is typically taken as $A_{e}=A_{i}\cup A_{d}$, e.g., when the mass-flux models
are discussed. This step is taken to analyse the asymmetry between up-and
down-welling motion of same strength. While $\langle\cdot\rangle_{x,y}$ are
horizontal averages with respect to the whole cross section $A_{\rm total}$,
$\langle\cdot\rangle_{u}$ denotes an average over the updraft and regions, see
also (4).

Figure 2 underlines our observation from Fig. 1 more quantitatively. The
buoyancy flux is arranged in narrower updraft and somewhat broader downdraft
regions in the convective boundary layer. Here, we show the probability
density functions (PDFs) for both cases at different times of the transient
growth. A direct comparison of the top and bottom panels shows that the tails
of the downdraft distribution have somewhat broader tails. Consequently,
downdrafts occupy slightly broader areas which underlines the top-down
asymmetry of the convective motion. Moreover, it can be seen that the PDFs
collapse increasingly better on each other in the core and parts of the tails
as $h/L_{0}$ grows. The insets in the Fig. 1(A) underline that the cell width
of the updraft patterns grows with time. This indicates the regime, where the
boundary layer growth rate becomes small compared to the turbulent mixing rate
by the convective motion.

Figure 3: Sketch of the architecture of the applied generative adversarial
network. The generator is composed of a U-shaped deep network, a U-Net.
Details, training procedure, and Wasserstein loss metrics are detailed in the
SI text.

##  IV Generative machine learning model

###  IV.1 Architecture

We introduce in brief the generative ML model as well the augmentation of the
training database. Figure 3 illustrates the concept of the generative
adversarial network approach: two convolutional neural networks compete in an
adversarial process and reach a final equilibrium state, similar to the Nash
equilibrium in a two-player non-cooperative game Nash (1950). Here, the
generator network has the architecture of a U-Net Ronneberger et al. (2015),
which applies additional residual connections between the contraction and
expansion paths of the network and is used for complex segmentation tasks with
few training examples Ronneberger et al. (2015); Fonda et al. (2019). It is
seeded by random latent variables and outputs a horizontal slice of
fluctuations of vertical velocity ${w^{\prime}}^{\rm GAN}(x,y)$ and buoyancy
${b^{\prime}}^{\rm GAN}(x,y)$. The second network, the discriminator or
critic, receives both the synthetic snapshot and a true DNS snapshot
${w^{\prime}}(x,y)$, ${b^{\prime}}(x,y)$. It then emits a score indicating
whether the current input snapshot stems from the DNS, labeled as real, or the
generator, labeled as fake. Finally, both networks are trained via stochastic
gradient descent. While the generator is trained to fool the critic network,
the latter aims to classify both snapshots correctly.

Figure 4: Comparison of the generated turbulent fields with the ground truth.
(A–C) Horizontal cross sections at $z=0.5h=0.584$ at $h/L_{0}=48.41$ of the
vertical velocity fluctuations $w^{\prime}$ (A), the buoyancy fluctuations
$b^{\prime}$ (B), and the local convective heat flux $w^{\prime}b^{\prime}$
(C) from DNS which is the ground truth (GT). (D–F) Corresponding cross
sections which have been generated by the generative adversarial network
(GAN). (G-I) Comparison of the probability density functions of the fields.

Here, we use the Wasserstein–GAN formulation introduced by Arjovsky et al.
Arjovsky et al. (2017) which applies a specific loss function for the
optimization of the networks. The Wasserstein loss helps to overcome an
initial training phase faster, where the generator network outputs sub-optimal
snapshots, which are easily identified by the critic. Furthermore, it
increases the overall stability of the network training process. Details on
the Wasserstein–GAN implementation and the training can be found in the  SI
Appendix. After the successful training, one uses the optimized generator
network only to produce synthetic snapshots of $w^{\prime}$ and $b^{\prime}$
at a given height and for different times without knowing the underlying
Boussinesq equations. In this way, we directly obtain the local buoyancy flux
field, $w^{\prime}b^{\prime}(x,y,z_{0},t)$, at a given height $z_{0}$ inside
the CBL and hence its horizontal average $\langle
w^{\prime}b^{\prime}\rangle_{x,y}$ without the need of a physical
parametrization model as described in (3). More importantly, this data-driven
convective parametrization will be generalizable, i.e., incorporating the
transient dynamics of the CBL since the network is trained with data that are
rescaled in correspondence with the self-similar mixing layer theory which we
detail in the following.

###  IV.2 Physics-informed data augmentation

We augment the training data in a re-normalization procedure that removes the
major first-order effects of the CBL growth. For this, we use snapshots in the
statistical quasi-steady phase of the CBL, where the spatial patterns of the
up- and downdrafts are not affected by the very initial fast CBL growth, see
again Fig. 2. In this regime, a new set of variables can be introduced that
embeds, to leading order, the transient temporal growth of the boundary layer
height. To this end, we rescale the turbulence fields with the similarity
variables which correspond to the mixing layer theory by Deardorff Deardorff
(1970). This gives the following transformation rules,

$\tilde{w}(\tilde{\bf x},\tilde{h}(t))=\frac{w({\bf
x},t)}{(B_{0}h(t))^{1/3}}\,,\;\;\;\;\tilde{b}(\tilde{\bf
x},\tilde{h}(t))=\frac{b({\bf x},t)}{(B_{0}^{2}/h(t))^{1/3}}\,,$ (5)

with rescaled time $\tilde{h}(t)=h(t)/L_{0}$ and the normalized coordinates
$\tilde{\bf x}={\bf x}/h(t)$. Note that this transformation leads to a
decreasing horizontal extent as the boundary layer grows. The training data of
the GAN is then chosen to be a horizontal slice at a plane $z/h(t)$ with
constant horizontal extent in normalized coordinates for all snapshots. The
latter step then involves a cropping procedure for most snapshots, which in
turn results in similar-sized spatial patterns of $\tilde{w}$ and $\tilde{b}$
and thus of $\tilde{w}^{\prime}\tilde{b}^{\prime}$. In summary, this results
in an augmented training database that encodes the effects of a growing
boundary layer to the first order. The similar spatial features are then
easily picked up by the generator network. Finally, we reverse the re-
normalization procedure to obtain synthetic snapshots at any given rescaled
time $h/L_{0}$. Further details on the augmentation are found in SI Appendix.

Figure 5: Prediction of the time evolution of the turbulence statistics by the
generative adversarial network (GAN). The probability density functions of the
vertical velocity fluctuations (A), the buoyancy fluctuations (B), and the
buoyancy flux (C) are shown for different times $h/L_{0}$. All quantities are
normalized by their corresponding root mean square values. The ground truth
(GT) is given in blue, the GAN result in red. The bottom row panels display
the time evolution of the moments of the buoyancy flux. These are the
fluctuations (D), the skewness (E), and the flatness (F). The blue line is for
the GT, the orange line for the GAN. Data are for $z/h=0.5$.

##  V Generated subgrid-scale fields

###  V.1 Structure and statistics

Figure 4(A–F) demonstrates that the GAN can produce strikingly real horizontal
slices of $w^{\prime}$, $b^{\prime}$, and $w^{\prime}b^{\prime}$ that
reproduce all spatial features of all (ground truth) fields, here for
$\tilde{z}=0.5$ at $h/L_{0}=48.41$. Moreover, the generated data reproduce the
highly non-Gaussian, skewed PDFs of all three fields almost perfectly, see
panels (G–I). Even far-tail events are captured fairly well in all three
cases. This has direct implications for a machine learning-based convective
parameterization, as information on the spatial organization of individual
plumes becomes available in contrast to the physical models based on means and
Gaussian statistics. Additionally, the derived PDFs allow to access the high-
order statistics of the fields, such as the strongest updrafts that shoot into
the stably stratified layer above the CBL. Similar results can be observed for
the distinct patterns and PDFs at $\tilde{z}=0.2$, slightly above the surface
layer, and for $\tilde{z}=1.0$, in the entrainment zone, as well as at
different times of the boundary layer growth, see  SI Appendix. For these
heights, we train two additional GANs with data from the respective horizontal
planes. In this way, the ML approach enables a fast generation of spatial and
statistical information at three important and representative heights, the
surface and mixed layer and the entrainment zone.

###  V.2 Transient time evolution

Figure 5 reports the turbulence statistics with respect to time. In panels
(A–C) of the figure, we show the PDFs of $w^{\prime}/\sigma_{w^{\prime}}$,
$b^{\prime}/\sigma_{b^{\prime}}$, and
$w^{\prime}b^{\prime}/\sigma_{w^{\prime}b^{\prime}}$ where the corresponding
$\sigma^{2}_{p}(t)=\langle p^{2}\rangle_{x,y}$ for
$p=\\{w^{\prime},b^{\prime},w^{\prime}b^{\prime}\\}$. It is seen that the GAN
generates the correct time evolution of the statistics and the corresponding
moments. In panels (D–F), we show the time evolution of the fluctuations
$\sigma_{w^{\prime}b^{\prime}}^{2}=\langle(w^{\prime}b^{\prime})^{2}\rangle_{x,y}$
and normalized third- and fourth order statistical moments of the buoyancy
flux versus time. The latter two are the skewness $S_{w^{\prime}b^{\prime}}$
and the flatness $F_{w^{\prime}b^{\prime}}$ which are given by

$S_{w^{\prime}b^{\prime}}(t)=\frac{\langle(w^{\prime}b^{\prime})^{3}\rangle_{x,y}}{\sigma_{w^{\prime}b^{\prime}}^{3}}\quad\mbox{and}\quad
F_{w^{\prime}b^{\prime}}(t)=\frac{\langle(w^{\prime}b^{\prime})^{4}\rangle_{x,y}}{\sigma_{w^{\prime}b^{\prime}}^{4}}\,.$
(6)

The figure shows that both normalized moments deviate strongly from the
Gaussian values which would follow to $S_{w^{\prime}b^{\prime}}=0$ and
$F_{w^{\prime}b^{\prime}}=3$. Again, the GAN reproduces the slow growth of the
moments fairly well. Note also that the fluctuations of the time series grows
with time since the sampling area for taking statistics shrinks with respect
to time. Flatness valus of the order of 50 indicate a highly intermittent
transport across the growing CBL.

###  V.3 Comparison with eddy diffusivity mass-flux

Finally, the GAN networks enable a computation of the mean turbulent buoyancy
flux $\langle w^{\prime}b^{\prime}\rangle$ at reference heights inside the
CBL. We now compare their performance with the ground truth, as well as the
estimations of two frequently used closures, see (3). Es mentioned before,
they are termed eddy diffusivity mass-flux (EDMF) and contain the mass-flux
parametrization as a second term next to the Boussinesq term which always
connects subgrid-scale flux (or stress) with a mean gradient via an eddy (or
turbulent) viscosity.

To this end, we applied an EDMF scheme which follows from a steady plume model
Siebesma et al. (2007). As detailed in the SI Appendix, the model comprises
(3) and the following two equations

$\displaystyle\frac{\partial\langle b\rangle_{u}}{\partial z}$
$\displaystyle=\epsilon(\langle b\rangle_{x,y}-\langle b\rangle_{u})\,,$ (7)
$\displaystyle\frac{1}{2}(1-2\mu)\frac{\partial\langle
w\rangle_{u}^{2}}{\partial z}$ $\displaystyle+C_{1}\epsilon\langle
w\rangle_{u}^{2}=\langle b\rangle_{u}-\langle b\rangle_{x,y}\,$ (8)

These equations still contain the adjustable model parameters $\mu$, $C_{1}$,
and $\epsilon$. Typically $\mu\approx 0.15$ and $C_{1}\approx 0.5$ are chosen.
The yet unknown entrainment rate $\epsilon$ is also height-dependent and
modeled by the following expression Siebesma et al. (2007)

$\displaystyle\epsilon(z)\simeq 0.4\left(\frac{1}{z}+\frac{1}{h-z}\right)\,.$
(9)

Equations (7)–(9) for the variables $\langle b(z)\rangle_{u}$ and $\langle
w(z)\rangle_{u}$ together with (3) are then solved numerically. Closure EDMF
takes finally $M_{u}\simeq a_{u}w_{u}(z)$ Soares et al. (2004) where
$a_{u}\approx 0.05$. See again SI Appendix where we also provide the model
formula for the turbulent diffusivity $\kappa_{t}$. The comparison of all
models is shown in Fig. 6 for three different times $h/L_{0}$. We find
excellent agreement between the three flux values predicted by the GAN and the
DNS profiles. Our ML results are also close to the profiles of the EDMF
parameterization.

##  VI Final discussion and outlook

We have presented a parametrization for a dry convective boundary layer on
sub-kilometer scales. Our convective parametrization is based on a generative
adversarial network with Wasserstein loss metrics which provides the
horizontal spatial organization and the statistics of vertical velocity
component and the buoyancy, and thus the one of the local buoyancy flux at
different heights inside the CBL. Our machine learning algorithm is able to
predict the statistical distributions of the buoyancy flux in the course of
the diurnal cycle which causes a growth of the layer height $h(t)$ and an
increasing entrainment into the stable atmosphere on top of the CBL. The
highly non-Gaussian and intermittent statistics of this process is
successfully reproduced by our generative convective parametrization.

We have used the self-similar growth of the CBL in correspondance with the
classical mixed-layer theory by Deardorff to augment the training data base
for the algorithm. In this way, we incorporated he leading-order physical
effects of the boundary layer growth into the ML algorithm. This enabled us to
generate synthetic snapshots of the growing CBL at three different heights
which the highly asymmetric turbulence when comparing updrafts in the form of
a narrow-ridge network and downdrafts in the form of broader downwellings in
between. The U-net architecture of the generator is able to generate highly
segmented fields with the correct spatial organisation of the convective
turbulence on the basis of a relatively small number of training data. This is
a significant improvement over conventional parameterization techniques like
the eddy diffusivity mass-flux, which cannot provide information on the
horizontal organization of strong individual plumes that overshoot in the
stable layer at the top and determine the entrainment rate. Moreover, as
opposed to stochastic parameterization schemes which are built on Gaussian
statistics, our ML approach reproduces the non-Gaussian distributions of
$w^{\prime}$, $b^{\prime}$, and $w^{\prime}b^{\prime}$. The trained generator
network is generalizable to different surface buoyancy fluxes and
stratification strengths of the free atmosphere, since the training is done in
self-similar units of the mixed layer theory.

Several directions for future research are possible from this point. The first
is to include wind shear into the present convective boundary layer model and
to extend the GAN. The cloud layer which typically start at the top of the CBL
can be included by a switch from the single buoyancy field as given by (1) to
a buoyancy, which is given by the liquid water potential temperature
$\theta_{l}$, in combination with the total water content. A further direction
is related to other applications of convection that include compressibility
effects, e.g. for astrophysical flows where data-driven mass-flux
parametrizations might substitute classical mixing length closures Kupka and
Muthsam (2017). These studies are underway and will be reported elsewhere.

Figure 6: Vertical profile of the buoyancy flux $\langle
w^{\prime}b^{\prime}\rangle/B_{0}$ at three values of $h/L_{0}$. The eddy
diffusivity mass-flux scheme uses a mass flux that is proportional to the
updraft velocity $M_{u}\sim w_{u}$, which is obtained from a plume model
(green curve), see also the SI Appendix. The generative network results at
three heights $z/h=0.2$, $0.5$, and $1.0$ are marked as red symbols.

## Acknowledgments

F.H. is supported by the project No. P2018-02-001 ”DeepTurb – Deep Learning in
and of Turbulence” of the Carl Zeiss Foundation. Partial support for J.P.M.
was provided by grant PID2019-105162RB-I00 funded by
MCIN/AEI/10.13039/501100011033. The authors gratefully acknowledge the Gauss
Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project
by providing computing time on the GCS Supercomputer SuperMUC-NG at Leibniz
Supercomputing Centre (www.lrz.de). We thank Christopher Müller for his
initial help with the machine learning algorithm.

##  Appendix A Convective boundary layer model

In this section, we collect the equations of motion which model the convective
boundary layer together with the most important dimensionless parameters and
the direct numerical simulation (DNS) method.

###  A.1 Equations and dimensionless parameters

We consider the Boussinesq limit of the Navier–Stokes equations which is given
by

$\displaystyle\frac{\partial{\bf v}}{\partial t}+{\bf\nabla}\cdot({\bf
v}\otimes{\bf v})$ $\displaystyle=-{\bf\nabla}\pi+\nu{\bf\nabla}^{2}{\bf
v}+b{\bf e}_{z}\,,$ (10) $\displaystyle{\bf\nabla}\cdot{\bf v}$
$\displaystyle=0\,,$ (11) $\displaystyle\frac{\partial b}{\partial
t}+{\bf\nabla}\cdot({\bf v}b)$ $\displaystyle=\kappa{\bf\nabla}^{2}b\,,$ (12)

where $\pi$ is a modified dynamic pressure or kinematic pressure, i.e., a
deviation from the pressure in the hydrostatic balance divided by a reference
density. Quantity $\nu$ is the kinematic viscosity, $\kappa$ is the buoyancy
diffusivity.

The velocity field ${\bf v}({\bf x},t)=(u,v,w)$ satisfies no-slip boundary
conditions at the bottom and free-slip boundary conditions at the top Garcia
and Mellado (2014). The buoyancy field $b({\bf x},t)$ satisfies Neumann
boundary conditions at the bottom and the top. The surface buoyancy flux is
given by $B_{0}$. Periodic boundary conditions are imposed in the horizontal
directions for all fields. The initial condition is ${\bf v}=0$ and, for the
buoyancy,

$b_{\rm bg}(z)=N_{0}^{2}z+f(z)\,,$ (13)

where $f(z)$ is the variation near the bottom surface to satisfy
$-\kappa\partial b/\partial z=B_{0}$. The function $f(z)$ is only non-zero
very near the surface and we consider long enough times for the effects of
$f(z)$ on the results to be negligibly small. The quantity $N_{0}$ is the
Brunt-Väisälä frequency in (13). The Prandtl number $Pr$ and the buoyancy
Reynolds number are the two central dimensionless parameters in this setup.
They are given by

$Pr=\frac{\nu}{\kappa}\,,\quad\quad Re_{0}=\frac{B_{0}}{\nu N_{0}^{2}}\,.$
(14)

Times are expressed in units of $N_{0}^{-1}$ and lengths in units of the
Ozmidov-type scale which is given by

$L_{0}=\left(\frac{B_{0}}{N_{0}^{3}}\right)^{1/2}\,.$ (15)

The Prandtl number is $Pr=1$ in the present DNS.

###  A.2 Convective boundary layer and encroachment height

In the convective boundary layer, there is no upper lid that confines the
motion; a turbulent boundary layer forms next to the bottom surface. The
difference to standard Rayleigh-Bénard convection is that here the boundary
layer is only caused by turbulent free convection, forced by a constant heat
flux at the bottom surface. The convective boundary layer (CBL) warms and
thickens with respect to time. The top boundary of the computational domain is
placed far enough to not influence the boundary-layer properties, and
otherwise the height of the computational domain becomes unimportant for the
analysis.

Except for the CBL, the density in the region above the boundary layer is
stably stratified, which hinders the growth of the boundary layer. This region
above the boundary layer is referred to as free troposphere or free
atmosphere; it is equivalent to the free-stream in a shear-driven boundary
layer. One effect of this stratification is that the transition region from
the turbulent boundary layer to the non-turbulent free stream is thinner than
in the case without stratification. This is because the stratification of the
free troposphere hinders the vertical displacement of fluid particles. This
transition region is referred to as the entrainment zone. Another effect is
that this stratification hinders the growth of the CBL thickness, to the point
that the characteristic time for changes of mean properties is small compared
to the large eddy turnover time, and a statistically quasi-steady state is
observed.

This stratification of the free-stream region is imposed by setting the
initial condition of the buoyancy field as $b_{\mathrm{bg}}=N_{0}^{2}z$, see
(13). The transition region between the turbulent boundary layer and the free
tropopause (i.e., the entrainment zone) occurs where the mean profile of
buoyancy, $\langle b(z)\rangle$, changes from the vertically constant value
inside the turbulent CBL to the linear profile $b_{\mathrm{bg}}=N_{0}^{2}z$ in
the non-turbulent free troposphere, see also Fig. 1(G) in the main text. We
refer to the vertical location of this thin transition region as CBL top.
There are various ways to define the CBL top and thus the thickness of the
whole convective boundary layer. To a first approximation, they all are
proportional to the encroachment height which is defined as Carson and Smith
(1975); Garcia and Mellado (2014)

$h(t)=\left\\{2N_{0}^{-2}\int_{0}^{z_{\infty}}(\langle
b(z)\rangle_{x,y}-b_{\mathrm{bg}}(z))\mathrm{d}z\right\\}^{1/2}\,.$ (16)

From the point of view of dimensional analysis, we have substituted the cell
height in a Rayleigh-Bénard setup by $L_{0}$ which was defined in (15). It
characterizes the thickness of the transition region at the top of the CBL.
The reference length $L_{0}$ is constant, but $h$ increases with time $t$.
Consider the mean buoyancy equation with statistical homogeneity in the
horizontal directions. Then follows from (12)

$\frac{\partial}{\partial t}[\langle b\rangle_{x,y}-b_{\rm
bg}(z)]=-\frac{\partial}{\partial z}\left[\langle
w^{\prime}b^{\prime}\rangle_{x,y}-\kappa\frac{\partial\langle
b\rangle_{x.y}}{\partial z}\right]\,,$ (17)

where we have added the time-independent profile $b_{\rm bg}(z)$ on the left
hand side. Integrate (17) from $z=0$ to $z=z_{\infty}$ and use the fact that
far away from the surface $\langle b\rangle_{x,y}\to b_{\rm
bg}(z)=N_{0}^{2}z$. This results to

$\frac{\partial}{\partial t}\int_{0}^{z_{\infty}}[\langle
b\rangle_{x,y}-b_{\rm bg}(z)]dz=\kappa
N_{0}^{2}+B_{0}=B_{0}\left(1+\frac{1}{PrRe_{0}}\right)\,.$ (18)

Thus an additional integration with respect to time, together with (16), gives

$\int_{0}^{z_{\infty}}[\langle b\rangle_{x,y}-b_{\rm
bg}(z)]dz=B_{0}\left(1+\frac{1}{PrRe_{0}}\right)(t-t_{0})=\frac{h^{2}N_{0}^{2}}{2}\,.$
(19)

Encroachment height and Ozmidov-type scale, see (15), are then related
according to

$\frac{h(t)}{L_{0}}=\left[2N_{0}(t-t_{0})\left(1+\frac{1}{PrRe_{0}}\right)\right]^{1/2}\;,$
(20)

where $t_{0}$ is an integration constant. See also (12) in the main text for
$Pr=1$. The growth rate of $h(t)$ is referred to as mean entrainment velocity

$w_{\mathrm{e}}\equiv\frac{\mathrm{dh}}{\mathrm{d}t}\,.$ (21)

The convective Rayleigh number, which is given by

$Ra_{c}(t)=\frac{B_{0}h(t)^{4}}{\nu\kappa^{2}}\,,$ (22)

increases as $h(t)$ increases with time $t$.

###  A.3 Direct numerical simulations

Discretization in space is performed using sixth-order spectral-like compact
finite differences on a structured Cartesian grid Mellado and Ansorge (2012).
A low-storage fourth-order accurate Runge–Kutta scheme is used for time
stepping. The discrete solenoidal constraint is satisfied to machine accuracy
using a Fourier decomposition along the periodic horizontal planes in $x$-$y$
directions and a factorization of the resulting set of equations along the
vertical coordinate $z$. The local ratio between the vertical grid spacing,
$\Delta z$, and the Kolmogorov scale, $\eta_{K}$, is $\Delta
z/\eta_{K}\lesssim 1.2$.

##  Appendix B Mass-flux parametrization

In this section, we describe in brief the basics of the mass-flux
parametrization which is used in many global numerical models of the Earth
system, such as the Integrated Forecasting System of the European Centre for
Medium-Range Weather Forecasts (www.ecmwf.int/en/publications/ifs-
documentation). The mass-flux parametrization is one part of the eddy
diffusivity mass-flux which models the vertical profiles of the subgrid-scale
fluxes of liquid-water potential temperature $\theta_{l}$ (convective heat
flux) or total water content $q_{t}$ (moisture flux) across the atmosphere.

###  B.1 Buoyancy flux closure by eddy diffusivity mass-flux

In the present proof-of-concept study we omit phase changes and cloud
formation. The liquid-water potential temperature $\theta_{l}$ reduces then to
the virtual potential temperature $\theta$ and the buoyancy $b$ becomes

$b({\bf x},t)=g\frac{\theta_{\rm v}({\bf x},t)-\theta_{{\rm
v},0}}{\theta_{{\rm v},0}}\,,$ (23)

with the acceleration due to gravity $g$ and the virtual potential temperature
at the surface $\theta_{{\rm v},0}$. The eddy diffusivity mass-flux closure
for the sub-grid scale buoyancy flux is given by Siebesma and Holtslag (1996);
Siebesma et al. (2007); Witte et al. (2022),

$\rho_{0}\langle
w^{\prime}b^{\prime}(z)\rangle_{x,y}=-\rho_{0}\kappa_{t}\frac{\partial\langle
b\rangle_{x,y}}{\partial z}+M_{u}(\langle b(z)\rangle_{u}-\langle
b(z)\rangle_{x,y})\,.$ (24)

Here, $\rho_{0}$ is the constant reference density in the lower atmospheric
boundary layer and $M_{u}$ is the unknown mass flux. As already said in the
main text, the horizontal cross section $A$ can be decomposed into disjoint
updraft regions (u) and the remaining environment (e), i.e., $A_{\rm
total}=A_{u}\cup A_{e}$. For simplicity we will denote $A_{\rm total}$ simply
by $A$ to ease notation. The first term in (24) is the standard Boussinesq
term for a turbulent stress which contains the eddy diffusivity $\kappa_{t}$
and the mean buoyancy gradient. The additional second term comprises the mass-
flux parametrization. In detail,

$\displaystyle\langle w^{\prime}b^{\prime}(z)\rangle_{x,y}$
$\displaystyle=\frac{1}{A}\int_{A}(w-\langle w\rangle_{x,y})(b-\langle
b\rangle_{x,y})\,dA\,,$ (25) $\displaystyle a_{u}\langle
w^{\prime}b^{\prime}(z)\rangle_{u}$
$\displaystyle=\frac{1}{A}\int_{A}(w-\langle w\rangle_{u})(b-\langle
b\rangle_{u})\,dA\,,$ (26) $\displaystyle(1-a_{u})\langle
w^{\prime}b^{\prime}(z)\rangle_{e}$
$\displaystyle=\frac{1}{A}\int_{A}(w-\langle w\rangle_{e})(b-\langle
b\rangle_{e})\,dA\,,$ (27)

with the area fractions $a_{u}=A_{u}/A$ and $1-a_{u}=A_{e}/A$. We neglect the
first gradient term on the right hand side of (24) for a while. The second
term of (24) can be rewritten to

$\langle w^{\prime}b^{\prime}\rangle_{x,y}=a_{u}\langle
w^{\prime}b^{\prime}\rangle_{u}+(1-a_{u})\langle
w^{\prime}b^{\prime}\rangle_{e}+a_{u}(1-a_{u})[\langle w\rangle_{u}-\langle
w\rangle_{e}][\langle b\rangle_{u}-\langle b\rangle_{e}]\,.$ (28)

The first term quantified the buoyancy flux connected with $A_{u}$, the second
the one with $A_{e}$, and the third is an exchange term due to organized
turbulence which connects both. The following approximations can be made: (i)
$a_{u}\ll 1$, in practice mostly about 5 %; (ii) $\langle
w\rangle_{e}\simeq\langle w\rangle_{x,y}$, and (iii) $\langle
b\rangle_{e}\simeq\langle b\rangle_{x,y}$. This results to the mass-flux
parametrization

$\langle
w^{\prime}b^{\prime}(z)\rangle_{x,y}\simeq\frac{M_{u}}{\rho_{0}}(\langle
b(z)\rangle_{u}-\langle b(z)\rangle_{x,y})\quad\mbox{with}\quad
M_{u}=\rho_{0}a_{u}[\langle w(z)\rangle_{u}-\langle w(z)\rangle_{x,y}]\,,$
(29)

with the undetermined mass flux $M_{u}$. Recall also that the mass flux
carries the physical dimension of kg m-2s-1. A plume model is required to
determine this unknown quantity which is discussed in the following.

###  B.2 Steady plume model

(12) of the present model is the counterpart of the set of prognostic
equations in a global circulation model. We apply a Reynolds decomposition and
obtain the following unclosed coarse-grid equation Siebesma and Cuijpers
(1995)

$\frac{\partial\langle b(z)\rangle_{x,y}}{\partial
t}=-\frac{1}{\rho_{0}}\frac{\partial}{\partial z}[\rho_{0}\langle
w^{\prime}b^{\prime}(z)\rangle_{x,y}]+\frac{\partial\langle
b(z)\rangle_{x,y}}{\partial t}\Bigg{|}_{\rm forcing}\,.$ (30)

All remaining terms due to horizontal mean advection and diffusion are
summarized in the last term of (30) as a forcing term. Note that we get in
this way an effective vertical transport model without any information on the
spatial horizontal mesoscale organisation of the convective turbulence which
the DNS in the main text display. The area decomposition into $A_{u}$ and
$A_{e}$ gives

$\langle b(z)\rangle_{x,y}=a_{u}\langle b(z)\rangle_{u}+(1-a_{u})\langle
b(z)\rangle_{e}\,.$ (31)

Equations (28) and (31) are now plugged into (30). We use that
$a_{u}(1-a_{u})\approx a_{u}$ since $a_{u}\ll 1$ and apply additionally the
approximations (ii) and (iii) from the last subsection. This leads to the
following updraft budget equation (omitting the $z$ argument),

$\displaystyle\rho_{0}\frac{\partial}{\partial t}[a_{u}\langle
b\rangle_{u}]=-\frac{\partial}{\partial z}[\rho_{0}a_{u}\langle
w^{\prime}b^{\prime}\rangle_{u}]-\frac{\partial}{\partial
z}[\underbrace{\rho_{0}a_{u}(\langle w\rangle_{u}-\langle
w\rangle_{x,y})}_{=M_{u}}\langle b\rangle_{u}]+E\langle b\rangle_{e}-D\langle
b\rangle_{u}+\rho_{0}a_{u}\frac{\partial\langle b\rangle_{x,y}}{\partial
t}\Bigg{|}_{\rm forcing}\,,$ (32)

and the corresponding downdraft budget equation

$\displaystyle\rho_{0}\frac{\partial}{\partial t}[(1-a_{u})\langle
b\rangle_{e}]=-\frac{\partial}{\partial z}[\rho_{0}(1-a_{u})\langle
w^{\prime}b^{\prime}\rangle_{e}]+\frac{\partial}{\partial z}[M_{u}\langle
b\rangle_{e}]-E\langle b\rangle_{e}+D\langle
b\rangle_{u}+\rho_{0}(1-a_{u})\frac{\partial\langle b\rangle_{x,y}}{\partial
t}\Bigg{|}_{\rm forcing}\,.$ (33)

We have added two terms in both equations to express the mass exchange due to
entrainment and detrainment with rates $E$ and $D$, respectively. Adding both
equations brings us back to (30). Furthermore, a continuity equation is
required which can be obtained from one of the budget equations by
substituting $\langle b\rangle$ with 1. This gives

$\rho_{0}\frac{\partial a_{u}}{\partial t}=-\frac{\partial M_{u}}{\partial
z}+E-D\,.$ (34)

Next a (statistically) steady state is assumed such that (32) and (34)
simplify to, see also ref. Siebesma et al. (2007),

$\displaystyle\frac{\partial}{\partial z}[M_{u}\langle b\rangle_{u}]$
$\displaystyle=-\frac{\partial}{\partial z}[\rho_{0}a_{u}\langle
w^{\prime}b^{\prime}\rangle_{u}]+M_{u}\epsilon\langle
b\rangle_{e}-M_{u}\delta\langle b\rangle_{u}\,,$ (35)
$\displaystyle\frac{\partial M_{u}}{\partial z}$
$\displaystyle=M_{u}(\epsilon-\delta)\,,$ (36)

where fractional entrainment and detrainment rates follow from
$E=M_{u}\epsilon$ and $D=M_{u}\delta$, respectively. When the first term on
the right hand side of (35) is set to zero, both equations can be combined to
give the following simple updraft equation

$\frac{\partial\langle b\rangle_{u}}{\partial z}=\epsilon(\langle
b\rangle_{x,y}-\langle b\rangle_{u})\,.$ (37)

Furthermore, we need an equation for $\langle w\rangle_{u}$ such that we can
close the buoyancy flux model Schumann and Moeng (1991); Siebesma et al.
(2007). This is obtained from the momentum balance for the updrafts, assuming
again a steady regime. This gives

$-\frac{1}{2}\frac{\partial\langle w\rangle_{u}^{2}}{\partial
z}-C_{1}\epsilon\langle
w\rangle_{u}^{2}-\frac{1}{\rho_{0}}\frac{\partial\langle
p\rangle_{x,y}}{\partial z}+\langle b\rangle_{u}-\langle b\rangle_{x,y}=0\,.$
(38)

This equation went again through a number of simplification steps: (i) the
first term results from the nonlinear advection of $\langle w\rangle_{u}$,
(ii) the second term encloses all remaining advection terms, (iii) the
effective entrainment rate which is connected with the turbulent mixing in
(ii) is assumed to be directly proportional to the entrainment for $\langle
b\rangle_{u}$ in (37). (iv) The pressure gradient term is approximated finally
to be $\partial(\mu\langle w\rangle_{u})/\partial z$. This is the last step
and the mass-flux parametrization comprises the following coupled system of
equations which would have to be solved together with the mean prognostic
equation for $\langle b\rangle_{x,y}$,

$\displaystyle\frac{\partial\langle b\rangle_{u}}{\partial z}$
$\displaystyle=\epsilon(\langle b\rangle_{x,y}-\langle b\rangle_{u})\,,$ (39)
$\displaystyle\frac{1}{2}(1-2\mu)\frac{\partial\langle
w\rangle_{u}^{2}}{\partial z}+C_{1}\epsilon\langle w\rangle_{u}^{2}$
$\displaystyle=\langle b\rangle_{u}-\langle b\rangle_{x,y}\,.$ (40)

If we know $\langle w\rangle_{u}$ and $\langle b\rangle_{u}$ then we can
calculate the second term on the right hand side of the EDMF model in (24) for
a given $a_{u}$. The model parameters $\epsilon$, $C_{1}$, and $\mu$ are
however unknown. Large eddy simulations studies by Siebesma et al. Siebesma et
al. (2007) suggest $C_{1}\simeq 0.5$ and $\mu\simeq 0.15$. The fractional
entrainment rate $\epsilon$, which will be height dependent, is more recently
included as a stochastic multi-plume model, i.e., for each plume $i$ the EDMF
model equations are solved. The entrainment rates $\epsilon_{i}$ are chosen
from a Poisson distribution Witte et al. (2022).

###  B.3 Numerical Implementation

The eddy-diffusivity $\kappa_{t}$ is implemented by a simple profile method
Holtslag and Duynkerke (1998) and results to

$\displaystyle\kappa_{t}(z)=k\left(39k\frac{z}{h}\right)^{1/3}\frac{z}{h}\left(1-\frac{z}{h}\right)^{2}w_{\ast}h\,,$
(41)

where $k\approx 0.4$ is the von Kármán constant. The mass-flux contributions
were obtained by integrating equations (39) and (40) using a fourth-order
Runge-Kutta scheme. We choose the starting height $z_{0}>0$ together with the
boundary conditions

$\displaystyle\langle b(z_{0})\rangle_{u}$ $\displaystyle=\langle
b(z_{0})\rangle+\alpha B_{0}/\sigma_{w},$ (42) $\displaystyle\langle
w(z_{0})\rangle_{u}^{2}$ $\displaystyle=\sigma_{w}^{2}(z_{0}),$ (43)

where we make use of an empirical standard deviation $\sigma_{w}$ as in ref.
Holtslag and Moeng (1991) which is given by

$\displaystyle\sigma_{w}(z)\simeq
1.3\left(0.6\frac{z}{h}\right)^{1/3}\left(1-\frac{z}{h}\right)^{1/2}w_{\ast},$
(44)

which fits well for the present direct numerical simulations as seen in Fig.
7.

Figure 7: Standard deviation of vertical velocity $w$ obtained from the direct
numerical simulation (blue) in comparison to the expression (44) (orange). The
GAN values are shown in red.

Furthermore, we use an entrainment rate profile as suggested in van Ulden
(1997); Siebesma et al. (2007) which is given by

$\displaystyle\epsilon\simeq 0.4\left(\frac{1}{z}+\frac{1}{h-z}\right).$ (45)

These input lead to the profiles in Fig. 5 of the main text which we used to
compare this closure with our machine learning results.

##  Appendix C Generative Adversarial Network

In this section, we describe the machine learning algorithm which is used to
generate synthetic data. We explain the architecture of the network, details
on the Wasserstein loss metrics, and our augmentation of the training data
base.

###  C.1 Generator and discriminator architecture

A Generative Adversarial Network (GAN) Goodfellow et al. (2014) is a neural
network comprising two interconnected and adversarial components: a generator
and a discriminator (or critic). The generator generates synthetic data, in
our case, horizontal slices of an atmospheric boundary layer simulation. The
discriminator’s role is to distinguish between the synthetic data produced by
the generator and real data from the DNS. The two components are trained
simultaneously in a competitive manner, with the generator striving to
generate data that can fool the critic and the critic continuously improving
its ability to classify the data correctly. This adversarial relationship
drives the GAN to produce increasingly realistic, high-quality synthetic data.
The training of the network is completed once the Nash equilibrium is reached
Nash (1950); the discriminator is unable to distinguish synthetically
generated and ground truth data. Then the generator can be decoupled from the
discriminator and operate autonomously.

Specifically, we use the Wasserstein–GAN (WGAN) formulation from Arjovsky et
al. Arjovsky et al. (2017), which solves the vanishing gradient problem of the
generator of the original GAN model. This is done using the Wasserstein
distance (or Earth-movers distance) between the synthetic and real data in the
WGAN loss function. In this way, the generator network will overcome the
initial phase, where the discriminator will quickly identify the generated
data samples and the GAN training ends in favor of the critic network.
Furthermore, the WGAN shows greater stability and convergence during the
training.

Here, we use a generator with a U-shaped deep neural network structure, in
short U-Net Ronneberger et al. (2015); Fonda et al. (2019). It has seven
contraction and seven subsequent expansion layers. After each contraction
layer, the spatial dimensions are halved via a convolution operation, while
the channel dimensions are doubled. The opposite is done in the expansion path
via transposed convolutions. Both paths are connected via residual
connections, allowing for a direct information flow between two layers in the
contraction and expansion paths of the U-shape. The activation function was
chosen as the Parametric Rectified Linear Unit He et al. (2015) (PReLU)

$\displaystyle{\rm PReLU}(x)=\left\\{\begin{array}[]{lr}x,&\text{if }x\geq
0\\\ ax,&{\rm otherwise}\end{array}\right.$ (48)

where the slope $a$ is a learnable parameter. The discriminator network
consists of eight convolution layers, which reduce the input data to a scalar
score. Here, the PReLU activation was used as well. Finally, to increase the
variance of the network, we use Dropout layers, which set the layer output to
zero with a probability of $0.3$.

###  C.2 Loss metrics

In the WGAN formulation, the output of the critic corresponds to a score that
grades the origin of the input. Therefore, the critic is trained to output a
higher score for the DNS input $\mathbf{x}$ than for its synthetic counterpart
$\mathbf{x}_{G}$. The generator on the other hand, is trained to maximize the
critic’s output for $\mathbf{x}_{G}$. The Wasserstein distance is incorporated
into the discriminator loss function using the Kantorovich-Rubinstein duality
Villani (2008). However, this restricts the critic $D(\cdot)$ to be a
1-Lipschitz function (i.e., the Lipschitz constant is one), which can be
enforced by an additional gradient penalty term Gulrajani et al. (2017). The
WGAN discriminator loss function then reads

$\displaystyle\mathcal{L}_{\rm D}=\underset{\mathbf{x}_{G}\sim P_{\rm
G}}{\mathbb{E}}\left[D(\mathbf{x}_{G})\right]-\underset{\mathbf{x}\sim P_{\rm
r}}{\mathbb{E}}\left[D(\mathbf{x})\right]+\lambda\underset{\hat{\mathbf{x}}\sim
P_{\hat{\mathbf{x}}}}{\mathbb{E}}\left[\left(\|\nabla_{\hat{\mathbf{x}}}D\left(\hat{\mathbf{x}}\right)\|_{2}-1\right)^{2}\right],$
(49)

where $\mathbf{x}$ is sampled from the real, i.e., DNS data distribution
$P_{\rm r}$ and $\mathbf{x}_{G}$ from the generator model distribution $P_{\rm
G}$. The latter is given by $G\left(\mathbf{\xi}\right)$ with the random
latent variable $\mathbf{\xi}\sim\mathcal{N}(0,1)$. The last term in (49) is
the gradient penalty term, which enforces the 1-Lipschitz condition on the
discriminator. The random samples $\hat{\mathbf{x}}$ are sampled from straight
lines $\chi\mathbf{x}+(1-\chi)\mathbf{x}_{G}$ between pairs of $\mathbf{x}$
and $\mathbf{x}_{G}$, as proposed in by Gulrajani et al. Gulrajani et al.
(2017). Here, $\chi\sim U[0,1]$, i.e., uniformly distributed between 0 and 1.
The gradient penalty weight $\lambda$ is a hyperparameter. Moreover, the
generator loss is given by

$\displaystyle\mathcal{L}_{\rm G}=-\underset{\mathbf{x}_{G}\sim P_{\rm
G}}{\mathbb{E}}[D(\mathbf{x}_{G})].$ (50)

Finally, the Wasserstein loss function, which is given by

$\displaystyle\mathcal{L}_{\rm WGAN}=\underset{\mathbf{x}\sim P_{\rm
r}}{\mathbb{E}}[D(\mathbf{x})]-\underset{\mathbf{x}_{G}\sim P_{\rm
G}}{\mathbb{E}}[D(\mathbf{x}_{G})]\,,$ (51)

reflects the WGAN performance. We summarize the WGAN training procedure in
Algorithm 1.

###  C.3 Physics-informed augmentation of training data base

As discussed in the main text, we renormalize the DNS data to leverage the
self-similar behavior of the CBL. For this, we renormalize the fluctuations of
vertical velocity and buoyancy fields by their corresponding convective scales
Deardorff (1970), see eqns. (4) in the main text, where the CBL height $h(t)$
is given by definition (1) in the main text. Further, the coordinates are
measured in units of $h(t)$, i.e., $\tilde{x}=x/h$, $\tilde{y}=y/h$,
$\tilde{z}=z/h$. By doing so, we fix the CBL height. This has implications for
the horizontal extent of the normalized snapshots as the aspect ratio of the
domain decreases over time. Hence, to compare two snapshots at two different
$h/L_{0}$ with the same horizontal extend, one has to crop the snapshot with
the smaller value of $h(t)/L_{0}$. In this way, we remove the first-order
effects of the transient CBL growth in a fixed plane $\tilde{z}={\rm const}$.
This process is illustrated in Fig. 8.

Figure 8: Renormalization and cropping procedure illustrated for the vertical
velocity fluctuations $\tilde{w}^{\prime}$ in plane $\tilde{z}=0.5$. (A–D)
Original simulation snapshots of $w^{\prime}(x,y)$. (E–H) Rescaled
fluctuations $\tilde{w}^{\prime}(\tilde{x},\tilde{y})$ in the new coordinate
system $(\tilde{x},\tilde{y},\tilde{z})$. (I–L)
$\tilde{w}^{\prime}(\tilde{x},\tilde{y})$ cropped to the domain
$\tilde{x},\tilde{y}\in[-\tilde{L},\tilde{L}]$ with $\tilde{L}\approx 8.5$ at
time $h/L_{0}=48.41$ of panel (L). The data preparation for the buoyancy
fluctuations proceeds along the same lines.

In panels (I–L), it can be seen that for sufficiently large values of
$h/L_{0}$, this procedure produces horizontal slices with comparable flow
patterns. The same does hold for the buoyancy fluctuations $b^{\prime}$ (not
shown). It is these flow patterns that the GAN will be trained with to
generate synthetic data.

The cutoff value of $h/L_{0}$, after which all re-normalized snapshots can be
considered for the training routine, is determined by the following procedure.
First, the individual up- and downdraft areas $A_{u}$ and $A_{d}$, given by
values $w>0.95\max(w)$ and $w<0.95\min(w)$ respectively, are retrieved using
the algorithm introduced in Suzuki and be (1985) via the OpenCV Itseez (2015)
library implementation. Finally, we computed the Kullback-Leibler divergence
(KLD) between the PDF $p$ of up- or downdrafts at $h/L_{0}$ with the final PDF
$p_{\rm final}$ at $h/L_{0}=48.41$. For the updrafts, it is then given by

$\displaystyle D_{\mathrm{KL}}(p\|p_{\rm
final})=\sum_{A_{u}}p(A_{u})\log\left(\frac{p(A_{u})}{p_{\rm
final}(A_{u})}\right)\,.$ (52)

Figure 9 shows the KLD for the three planes over $h/L_{0}$. The slope of the
KLD becomes less steep near $h/L_{0}=35.87$ (start of green shaded area).
After this value, we consider all snapshots for the GAN training.

Figure 9: Kullback-Leibler divergence between the probability density
functions of normalized updraft (downdraft) areas $A_{u}$ ($A_{d}$) at
$h/L_{0}$ and the corresponding final probability density function. The green
shaded area marks the snapshots that were considered for the training of the
GAN.

###  C.4 Network training

The normalized and cropped horizontal slices at the three heights are
subsequently used to train three independent generative networks. For each
height, $550$ snapshots are available for training. Moreover, we use data
augmentation to increase the number of training samples to $32\cdot
550=17600$. For this, we consider a quarter of the renormalized DNS data as
discriminator input. Hence, also the generator will be trained to generate a
quarter of the original horizontal extent. Furthermore, we add rotations by
angles of $\pi/2,\pi,3\pi/2$, reflections about the $\tilde{x}$\- and
$\tilde{y}$-axes, as well as combinations of $\pi/2$-rotations and
$\tilde{x}/\tilde{y}$-reflections. This significantly increased the number of
training samples, which resulted in better generator performances.

The GANs were trained for $2000$ epochs using the Adaptive moments (Adam)
optimizer Kingma and Ba (2017). An early stopping routine ceased training
after $300$ epochs of no improvement of the Wasserstein metric
$\mathcal{L}_{\rm WGAN}$. Each run was performed on two NVIDIA A100 Tensor-
Core GPUs and took up to 24 h to finish. The choice of GAN hyperparameters for
the three planes is listed in Tab. 1. Figure 10 shows additionally the losses
for both, the generator and the discriminator, together with the Wasserstein
loss for all three GANs.

plane $\tilde{z}$ | learning rate | batch size | $N_{\rm D}/N_{\rm G}$ | $\lambda$
---|---|---|---|---
$0.2$ | $2\cdot 10^{-5}$ | $64$ | $12$ | $11$
$0.5$ | $2\cdot 10^{-5}$ | $64$ | $11$ | $10$
$1.0$ | $5\cdot 10^{-5}$ | $64$ | $90$ | $10$
Table 1: Neural network hyperparameters. The learning rate regulates the speed
of the gradient descent update step. The batch size refers to the number of
training examples processed together in a single forward and backward pass.
The parameter $N_{\rm D}/N_{\rm G}$ signifies the degree to which the critic
network undergoes additional training compared to the generator network. The
gradient penalty parameter $\lambda$ controls the strength of the 1-Lipschitz
condition in (49). Figure 10: Training loss of the three neural networks at
height $z/h=0.2$ (A), $z/h=0.5$ (B) and $z/h=1.0$ (C). The discriminator
(blue) and generator (orange) loss $\mathcal{L}_{\rm D}$ and $\mathcal{L}_{\rm
G}$ converge towards individual values. The Wasserstein loss $\mathcal{L}_{\rm
WGAN}$ (green) measures the quality of the GAN results. Hence the parameters
of the GAN are saved at the minimum Wasserstein loss value.

##  Appendix D Additional results

In this section, we supply further results that supplement the discussion of
the main text. Supporting Figs. 11, 12, and 13 compare output of the GAN with
the ground truth from the direct numerical simulations at different vertical
positions $z$.

Figure 11: Horizontal slices at $h/L_{0}=34.41$ in plane $z=0.5h=0.43$.
Compare to Fig. 4 in the main text. Figure 12: Horizontal slices at
$h/L_{0}=48.41$ in plane $z=0.2h=0.233$. Compare to Fig. 4 in the main text.
Figure 13: Horizontal slices at $h/L_{0}=48.41$ in plane $z=1.0h=1.167$.
Compare to Fig. 4 in the main text.

## References

  * Lenschow and Stephens (1980) D. H. Lenschow and P. L. Stephens, Boundary Layer Meteorol. 19, 509 (1980).
  * Bodenschatz et al. (2010) E. Bodenschatz, S. P. Malinowski, R. A. Shaw, and F. Stratmann, Science 327, 970 (2010).
  * Mellado (2017) J. P. Mellado, Annu. Rev. Fluid Mech. 49, 145 (2017).
  * Atkinson and Wu Zhang (1996) B. W. Atkinson and J. Wu Zhang, Rev. Geophys. 34, 403 (1996).
  * Bony et al. (2020) S. Bony, H. Schulz, J. Vial, and B. Stevens, Geophys. Res. Lett. 47, e2019GL085988 (2020).
  * Madden and Julian (1972) R. A. Madden and P. R. Julian, J. Atmos. Sci. 29, 1109 (1972).
  * Miura et al. (2007) H. Miura, M. Satoh, T. Nasuno, A. T. Noda, and K. Oouchi, Science 318, 1763 (2007).
  * Ertel (1942) H. Ertel, Meteor. Z. 59, 1690 (1942).
  * Deardorff (1972) J. W. Deardorff, J. Geophys. Res. 77, 5900 (1972).
  * Arakawa (1974) W. H. Arakawa, A.and Schubert, J. Atmos. Sci. 31, 674 (1974).
  * LeMone et al. (2019) M. A. LeMone, W. M. Angevine, C. S. Bretherton, F. Chen, J. Dudhia, E. Fedorovich, K. B. Katsaros, D. H. Lenschow, L. Mahrt, E. G. Patton, et al., Meteorol. Monogr. 59, 1 (2019).
  * Tiedtke (1989) M. Tiedtke, Mon. Weather Rev. 117, 1779 (1989).
  * Siebesma and Cuijpers (1995) A. P. Siebesma and J. W. M. Cuijpers, J. Atmos. Sci. 52, 650 (1995).
  * Želi et al. (2021) V. Želi, G. Brethouwer, S. Wallin, and A. V. Johansson, Boundary Layer Meteorol. 178, 487 (2021).
  * Wyngaard (2010) J. C. Wyngaard, _Turbulence in the Atmosphere_ (Cambridge University Press, Cambridge, UK, 2010).
  * Klingebiel et al. (2021) M. Klingebiel, H. Konow, and B. Stevens, J. Atmos. Sci. 78, 3205 (2021).
  * Vogel et al. (2022) R. Vogel, A. L. Albright, J. Vial, G. George, B. Stevens, and S. Bony, Nature 612, 696 (2022).
  * Ahlgrimm and Forbes (2014) M. Ahlgrimm and R. Forbes, Mon. Weather Rev. 142, 668 (2014).
  * Becker et al. (2021) T. Becker, P. Bechtold, and I. Sandu, Q. J. R. Meteorol. Soc. 4388–4407 (2021).
  * Hohenegger et al. (2023) C. Hohenegger, P. Korn, L. Linardakis, R. Redler, R. Schnur, P. Adamidis, J. Bao, S. Bastin, M. Behravesh, M. Bergemann, et al., Geosci. Model Dev. 16, 779 (2023).
  * Siebesma et al. (2007) A. P. Siebesma, P. M. M. Soares, and J. Teixeira, J. Atmos. Sci. 64, 1230 (2007).
  * Brenowitz and Bretherton (2018) N. D. Brenowitz and C. S. Bretherton, Geophys. Res. Lett. 45, 6289 (2018).
  * Reichstein et al. (2019) M. Reichstein, G. Camps-Valls, B. Stevens, M. Jung, J. Denzler, N. Carvalhais, and Prabhat, Nature 566, 195 (2019).
  * Kashinath et al. (2021) K. Kashinath, M. Mustafa, A. Albert, J.-L. Wu, C. Jiang, S. Esmaeilzadeh, K. Azizzadenesheli, R. Wang, A. Chattopadhyay, A. Singh, et al., Phil. Trans. R. Soc. A 379, 20200093 (2021).
  * Fonda et al. (2019) E. Fonda, A. Pandey, J. Schumacher, and K. R. Sreenivasan, Proc. Natl. Acad. Sci. USA 116, 8667 (2019).
  * Rasp (2020) S. Rasp, Bull. Am. Meteorol. Soc. 101, E1980 (2020).
  * Beck and Kurz (2021) A. Beck and M. Kurz, GAMM Mitteilungen 44, e202100002 (2021).
  * Yuval et al. (2021) J. Yuval, P. A. O’Gorman, and C. N. Hill, Geophys. Res. Lett. 48, e2020GL091363 (2021).
  * Cheng et al. (2022) Y. Cheng, M. G. Giometto, P. Kauffmann, L. Lin, C. Cao, C. Zupnick, H. Li, Q. Li, Y. Huang, R. Abernathey, et al., J. Adv. Model. Earth Syst. 14, e2021MS002847 (2022).
  * Guan et al. (2023) Y. Guan, A. Subel, A. Chattopadhyay, and P. Hassanzadeh, Physica D 443, 133568 (2023).
  * Bretherton et al. (2022) C. S. Bretherton, B. Henn, A. Kwa, N. D. Brenowitz, O. Watt-Meyer, J. McGibbon, W. A. Perkins, S. K. Clark, and L. Harris, J. Adv. Model. Earth Syst. 14, e2021MS002794 (2022).
  * Heyder and Schumacher (2021) F. Heyder and J. Schumacher, Phys. Rev. E 103, 053107 (2021).
  * Heyder et al. (2022) F. Heyder, J. P. Mellado, and J. Schumacher, Phys. Rev. E 106, 055303 (2022).
  * Dai and Hemri (2021) Y. Dai and S. Hemri, Mon. Weather Rev. 149, 3923 (2021).
  * Leinonen et al. (2019) J. Leinonen, A. Guillaume, and T. Yuan, Geophys. Res. Lett. 46, 7035 (2019).
  * Price and Rasp (2022) I. Price and S. Rasp, PMLR 151, 10555 (2022).
  * Shamekh and Gentine (2023) S. Shamekh and P. Gentine, J. Adv. Model. Earth Syst. pp. 1–33, submitted (2023).
  * Goodfellow et al. (2014) I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, arXiv:1406.2661 (2014).
  * Goodfellow et al. (2016) I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_ (MIT Press, Cambridge, USA, 2016).
  * Mellado and Ansorge (2012) J. P. Mellado and C. Ansorge, Z. Angew. Math. Mech. 92, 1 (2012).
  * Garcia and Mellado (2014) J. R. Garcia and J. P. Mellado, J. Atmos. Sci. 71, 1935 (2014).
  * Brunton et al. (2020) S. L. Brunton, B. R. Noack, and P. Koumoutsakos, Annu. Rev. Fluid Mech. 52, 477 (2020).
  * Ronneberger et al. (2015) O. Ronneberger, P. Fischer, and T. Brox, Lecture Notes in Computational Science 9351, 234 (2015).
  * Witte et al. (2022) A. Witte, A. Herrington, J. Teixeira, M. J. Kurowski, M. J. Chinita, R. L. Storer, K. Suselj, G. Matheou, and J. Bacmeister, Mon. Weath. Rev. 150, 2339 (2022).
  * Nash (1950) J. F. Nash, Proc. Natl. Acad. Sci. USA 36, 48 (1950).
  * Arjovsky et al. (2017) M. Arjovsky, S. Chintala, and L. Bottou, arXiv:1701.07875 (2017).
  * Deardorff (1970) J. W. Deardorff, J. Atmos. Sci. 27, 1211 (1970).
  * Soares et al. (2004) P. M. M. Soares, P. M. A. Miranda, A. P. Siebesma, and J. Teixeira, Q. J. R. Soc. Meteorol. 130, 3365 (2004).
  * Kupka and Muthsam (2017) F. Kupka and H. J. Muthsam, Living Rev. Comput. Astrophys. 3, 1 (2017).
  * Carson and Smith (1975) D. J. Carson and F. B. Smith, in _Advances in Geophysics_ (Elsevier, 1975), vol. 18, pp. 111–124.
  * Siebesma and Holtslag (1996) A. P. Siebesma and A. A. M. Holtslag, J. Atmos. Sci. 53, 2354 (1996).
  * Schumann and Moeng (1991) U. Schumann and C.-H. Moeng, J. Atmos. Sci. 48, 1758 (1991).
  * Holtslag and Duynkerke (1998) A. Holtslag and P. G. Duynkerke (Royal Netherlands Academy of Arts and Sciences, 1998).
  * Holtslag and Moeng (1991) A. A. M. Holtslag and C.-H. Moeng, Journal of the Atmospheric Sciences 48, 1690 (1991), URL https://doi.org/10.1175/1520-0469(1991)048<1690:edacti>2.0.co;2.
  * van Ulden (1997) A. P. van Ulden, A. P.and Siebesma, Amer. Meteor. Soc p. 257–259 (1997).
  * He et al. (2015) K. He, X. Zhang, S. Ren, and J. Sun, arXiv:1502.01852 (2015).
  * Villani (2008) C. Villani, _Optimal Transport: Old and New_ , Grundlehren der mathematischen Wissenschaften (Springer, Berlin, 2008), ISBN 9783540710509.
  * Gulrajani et al. (2017) I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, and A. C. Courville, in _Advances in Neural Information Processing Systems_ , edited by I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Curran Associates, Inc., 2017), vol. 30.
  * Suzuki and be (1985) S. Suzuki and K. be, Comput. Graph. Image Process. 30, 32 (1985), URL https://doi.org/10.1016/0734-189x(85)90016-7.
  * Itseez (2015) Itseez, _Open source computer vision library_ , https://github.com/itseez/opencv (2015).
  * Kingma and Ba (2017) D. P. Kingma and J. Ba, arXiv:1412.6980 (2017).