# Large-eddy simulations of marine boundary-layer clouds associated with cold
air outbreaks during the ACTIVATE campaign– part 1: Case setup and
sensitivities to large-scale forcings 22footnotemark: 2

###### Abstract

Large-eddy simulation (LES) is able to capture key boundary-layer (BL)
turbulence and cloud processes. Yet, large-scale forcing and surface turbulent
fluxes of sensible and latent heat are often poorly constrained for LES
simulations. We derive these quantities from measurements and reanalysis
obtained for two cold air outbreak (CAO) events during Phase I of the Aerosol
Cloud meTeorology Interactions oVer the western ATlantic Experiment (ACTIVATE)
in February-March 2020. We study the two contrasting CAO cases by performing
LES simulations and test the sensitivity of BL structure and clouds to large-
scale forcings and turbulent heat fluxes. Profiles of atmospheric state and
large-scale divergence and surface turbulent heat fluxes obtained from the
reanalysis data ERA5 agree reasonably well with those derived from ACTIVATE
field measurements for both cases at the sampling time and location.
Therefore, we adopt the time evolving heat fluxes, wind and advective
tendencies profiles from ERA5 reanalysis data to drive the LES simulations. We
find that large-scale thermodynamic advective tendencies and wind relaxations
are important for the LES to capture the evolving observed BL meteorological
states characterized by the hourly ERA5 reanalysis data and validated by the
observations. We show that the divergence (or vertical advection) is important
in regulating the BL growth driven by surface heat fluxes in LES simulations.
The evolution of liquid water path is largely affected by the evolution of
surface heat fluxes. The liquid water path simulated in LES agrees reasonably
well with the ACTIVATE measurements. This study paves the path to investigate
aerosol-cloud-meteorology interactions using LES informed and evaluated by
ACTIVATE field measurements.

##  1 Introduction

When viewed from space, about $70\%$ of Earth’s surface is covered by clouds
(Schneider et al. 2017). Clouds, the regulator of the radiative heating and
cooling of the planet (Ramanathan et al. 1989), represent a major complication
in the current modeling of the climate system (Schneider et al. 2017; Stevens
and Bony 2013; Bony et al. 2017). One of the most challenging problems of
cloud–climate interactions is to understand how cloud microphysical processes
affect atmospheric water and radiation budgets, such as how precipitation
efficiency affects radiative properties of stratocumulous clouds (Boucher et
al. 2013).

The Western North Atlantic Ocean (WNAO) region has attracted decades of
atmospheric research due to the complex atmospheric system (Painemal et al.
2021), pollution outflow from North America (Corral et al. 2021), and
accessibility by aircraft and ships. However, the subject of Aerosol-Cloud-
Interaction (ACI) is the least investigated among all the field campaign
measurements over the WNAO (Sorooshian et al. 2020) partly because of the
complicated chemical, physical, and dynamical processes in this region. ACI
involves processes from the formation of nm-sized aerosols to the life cycle
of kilometer-sized clouds, which covers a scale range of about $10^{12}$. Such
a scale separation coupled with turbulence poses great challenge for both
measurements and numerical modeling. The spatial distribution of aerosols and
the ambient humidity fields determine the formation of cloud droplets and ice
crystals and their size distribution (Shaw 2003). Precipitation and radiative
properties of clouds are altered by the size distribution of particles. The
Aerosol Cloud meTeorology Interactions oVer the western ATlantic Experiment
(ACTIVATE) field campaign aims to tackle ACI by performing comprehensive
measurements of cloud macro/micro properties and atmospheric states using two
aircraft simultaneously, which can be used to evaluate and constrain
atmospheric models (Sorooshian et al. 2019, 2020).

Large eddies on the order of 100 m in the planetary boundary layer are
important for turbulent mixing, heat/moisture transport, and cloud formation.
Large-eddy simulation (LES) has been widely used to model marine boundary-
layer clouds and ACI (Bretherton et al. 1999; Stevens et al. 2002; Brown et
al. 2002; Ackerman et al. 2009; Wang et al. 2009; Wang and Feingold 2009; de
Roode et al. 2019). LES resolves the intermediate and large
($\mathcal{O}(100\,\rm{m})$) turbulent eddies and parameterizes smaller
scales. One of the most challenging problems of using idealized LES with
doubly periodic boundary conditions to represent evolving clouds is
realistically configuring large-scale forcings (e.g., horizontal advection
tendencies, divergence (D) of flow, and surface heat fluxes) that determine
the spatiotemporal variation of large-scale ambient conditions for the cloud
system. The “large-scale” here refers to scales of $50-500\,\rm{km}$ (Bony and
Stevens 2019). There are different ways to construct large-scale forcings for
typical LES domains. For example, $D$ profiles, horizontal advection
tendencies and heat fluxes can be obtained from analysis/reanalysis products
and other numerical models, as well as measurements of variables used for the
calculation. To validate $D$ and heat fluxes obtained from numerical models,
observations such as sounding profiles of atmospheric state and surface
temperature measurements are needed. Endo et al. (2015) constructed continuous
large-scale and surface forcings from reanalysis data to successfully simulate
continental boundary layer clouds during the Routine Atmospheric Radiation
Measurement Aerial Facility Clouds with Low Optical Water Depths Optical
Radiative Observations (RACORO) campaign using an idealized LES model. The
application of such large-scale forcing schemes for marine stratocumulus
clouds in the WNAO region, however, has not been reported, which is a focus of
this study.

Marine stratocumulus clouds associated with cold air outbreaks (CAOs) with
mesoscale (scales larger than a few kilometers) fluctuations are challenging
to represent in climate models. CAO occurs when cold air mass moves over a
warm sea surface, creating strong convection analogous to Rayleigh-Bénard
convection (Agee 1987). CAO events are characterized by stronger surface
latent heat fluxes (up to $1200\,\rm{Wm}^{-2}$) and subsidence of up to an
order of magnitude, compared to non-CAO cases due to the large temperature
difference (Agee 1987). This can makes it difficult to simulate convection and
clouds associated with CAOs. The ratio of buoyancy force to shear (Richardson
number), precipitation, and entrainment contribute to the topological cloud
structure of CAOs (e.g., cloud streets) (de Roode et al. 2019). The cloud roll
structure occurring during CAOs was first simulated by Liu et al. (2004) using
a cloud-resolving model able to capture the transition of clouds from two-
dimensional roll structure to three-dimensional closed cells. Gryschka and
Raasch (2005) performed the first LES to simulate CAO cloud roll structures.
They found that roll structures are sensitive to the spatial resolution of
LES. Tomassini et al. (2017) investigated how well a CAO event over the North
Atlantic Ocean is represented in global models as compared to LES. They found
that the global models employed in their study underestimate the amount of
cloud liquid water compared to LES results. More recently, de Roode et al.
(2019) performed an LES intercomparison study of a CAO case observed during
the CONSTRAIN campaign. They found that the evolution of the stratocumulus
cloud deck and the timing of its breakup differ significantly among seven LES
models. Tornow et al. (2021) investigated a marine CAO case in the
northwestern Atlantic and found that frozen hydrometeors accelerate the
transition of cloud decks into broken cloud streets.

In this first part of two companion studies, we first introduce two CAO cases
sampled during the 2020 winter deployment of ACTIVATE and describe the
numerical experiment setup for idealized LES to model the two cases. Then we
use divergence profiles and surface heat fluxes derived from ACTIVATE
dropsondes and sea surface temperature (SST) measurements to first evaluate
these quantities from ERA5 reanalysis data. We further examine the
sensitivities of LES results surface heat fluxes and large-scale thermodynamic
advective tendencies from ERA5. We adopt the same LES model and large-scale
forcing scheme as in Endo et al. (2015). This part focuses on contrasting
meteorological and large-scale ambient conditions between the two CAO cases
and their impact on boundary-layer (BL) clouds. The second part of the two
companion studies will tackle the aerosol-clouds-interactions.

##  2 Observations, reanalysis data, and LES numerical experiment setup

###  2.1 ACTIVATE campaign

The ACTIVATE field campaign aims to collect sufficient measurements to
understand interactions of marine boundary-layer clouds with meteorological
conditions and aerosol particles, which eventually leads to improved physical
understanding of cloud micro/macro processes and reduced uncertainty in their
representation in global climate models. A total of 150 coordinated flights
with two air-borne platforms is planned for three years (2020-2022) over the
western North Atlantic Ocean ($25^{\circ}-50^{\circ}$N,
$60^{\circ}-85^{\circ}$W) to characterize aerosol-cloud-meteorology
interactions in a systematic and simultaneous manner (Sorooshian et al. 2019).
This is being achieved by flying two aircraft simultaneously at different
altitudes. The low-flying HU-25 Falcon measures in-situ trace gases, aerosol,
clouds, precipitation, and meteorological properties below, in, and above
clouds. The higher-flying King Air above clouds simultaneously acquires remote
retrievals of aerosols and clouds while launching dropsondes.

During the two CAO process-study cases over the WNAO region on 28 February
2020 and 1 March 2020, corresponding to Research Flight $\\#10$ and $\\#13$
(hereafter RF10 and RF13), respectively, 11 dropsondes (model Vaisala NRD41)
were released from the King Air. Each of them provided vertical profiles of
air pressure $p$, temperature $T$, relative humidity RH, and horizontal
velocities $u$ and $v$ with a vertical resolution of 5-10 m and an accuracy
(with associated uncertainty) of 0.01 hPa ($\pm 0.5$ hPa), $0.01^{\circ}C$
($\pm 0.2^{\circ}C$), $0.01\%$ ($\pm 0.3\%$), and $0.01\,{\rm m\,s}^{-1}$
($\pm 0.5\,{\rm m\,s}^{-1}$) (National Center for Atmospheric Research 2021),
respectively. The King Air flew in a circular pattern with a diameter of about
$152\,{\rm km}$ to cover the largest enclosed area for dropsonde measurements
and to avoid sharp turns. Such a flight pattern for the dropsonde measurements
was first proposed by Lenschow et al. (1999). Dropsondes were released at a
height of about 8 km. We interpolate the measured data evenly with a vertical
spacing of 10 m for further analysis. Two contrasting CAO cases were observed
over the WNAO region on February 28 (RF10, dropsonde-circle center at
$33.66^{\circ}$N, $286.69^{\circ}$E) and March 1 (RF13, dropsonde-circle
center at $38.01^{\circ}$N, $288.36^{\circ}$E) as shown in Figure 1. Table 1
summarizes the start/end time of dropsonde measurements, the location, $10$-m
wind speed $U_{10\rm m}$, $q_{v,\rm{10m}}$, $T_{\rm{10m}}$, and ERA5 SST at
the center of dropsonde circle for the February 28 and March 1 cases,
respectively.

Table 1: This table lists the start-end time of dropsonde measurements, the location, $10$-m wind speed $U_{10\rm m}$, $q_{v,\rm{10m}}$, and $T_{\rm{10m}}$ at the center of dropsonde circle for the 28 February 2020 and 1 March 2020 cases. The corresponding SST from ERA5 reanalysis is also documented.  Case | Start-end times (UTC) | Lat (∘) | Lon (∘) | $U_{10\rm{m}}\,(\rm{m\,s}^{-1})$ | $q_{v,10\rm{m}}$ ($\rm{g\,kg}^{-1}$) | $T_{10\rm{m}}\,(K)$ | SST ($K$)
---|---|---|---|---|---|---|---
MW-IR | ERA5
0228 | 1540-1645 | 33.66 | 286.69 | 5.94 | 4.64 | 286.75 | 293.25 | 293.60
0301 | 1451-1547 | 38.01 | 288.36 | 11.74 | 2.79 | 275.71 | 286.84 | 286.88

Dropsonde measurements are used to characterize the meteorological conditions
and derive large-scale divergence and surface heat fluxes for both cases.
Cloud droplets and ice crystals were observed for both cases. The mean number
concentration of cloud droplets obtained from Fast Cloud Droplet Probe (FCDP)
measurement (Taylor et al. 2019; Knop et al. 2021) is about $\langle
N_{c}\rangle=650\,{\rm cm}^{-3}$ for the February 28 case and $\langle
N_{c}\rangle=450\,{\rm cm}^{-3}$ for the March 1 case. These values are
acquired by averaging in-cloud FCDP measurement with a lower cutoff of liquid
water path of $0.02\,{\rm g\,kg}^{-1}$ and effective diameter of $3.5\,\mu{\rm
m}$. There were also detailed measurements of aerosol particles including mass
and number concentration, composition, size distribution, hygroscopicity, and
optical properties. Given the focus of this study, we only use the mean cloud
drop number in our LES sensitivity simulations on meteorological conditions
and large-scale forcings. Liquid water path is retrieved from Research
Scanning Polarimeter (RSP) (Alexandrov et al. 2012, 2018). Given the
instantaneous field of view of 14 mrad, typical cloud tops (about 2 km), and a
flight altitude of the King Air during ACTIVATE (8 9 km), the nadir pixel size
of the RSP is approximately 100 m. To compare with the LES simulation with a
300 m horizontal grid spacing, we average the RSP sampling every 3 s, given
that the moving speed of King Air is about $100\,\rm{m\,s}^{-1}$.

###  2.2 ERA5 and MERRA-2 reanalysis data

The ERA5 reanalysis data are generated using the fifth generation of European
Centre for Medium-Range Weather Forecasts’s Integrated Forecast System
(Hersbach et al. 2020). We use the ERA5 hourly data at a horizontal resolution
of 31 km. For three-dimensional fields, there are 137 model levels up to a
height of 80 km. Since ERA5 only provide $\bar{D}$ in datasets with specified
pressure levels, we use pressure-level data for the comparison of $\bar{D}$
and the corresponding large-scale vertical velocity $w$. The ERA5 large-scale
forcings for the LES are obtained at the model levels instead of the pressure
levels. The ambient meteorological conditions for a given location during a
CAO usually evolve quickly due to strong winds and large surface heat fluxes
under winter mid-latitude weather disturbances. Since the measurement time
window for the two CAO events is about one hour during the ACTIVATE field
campaign, we are not able to use the measurements directly to drive the LES
simulation for many hours. We validate the ERA5 reanalysis data against the
limited field measurements and then use the evolving forcing conditions from
ERA5 reanalysis data to drive the LES simulations.

The Modern-Era Retrospective analysis for Research and Applications version 2
(MERRA-2) (Global Modeling and Assimilation Office (2015) GMAO) is also used
to compare with the ERA5 reanalysis and dropsonde measurements. The MERRA-2
reanalysis data are generated using the Goddard Earth Observation System
version 5 (GEOS-5) with its Data Assimilation System version 5.12.4 (Gelaro et
al. 2017). MERRA-2 has a horizontal resolution of $0.5^{\circ}\times
0.625^{\circ}$ with 72 model levels and provides 3-hourly datasets at 42
pressure levels and 1-hourly two-dimensional datasets.

###  2.3 Satellite measurements

We use daily sea surface temperature (SST) retrieved from microwave and
infrared based satellite measurements (MW-IR SST) produced by Remote Sensing
Systems (National Centers for Environmental Information 2008). The SST product
has a horizontal grid spacing of 9 km. This resolution is 3 times higher than
the SST from ERA5 reanalysis data.

###  2.4 LES numerical experiment setup

We use the Weather Research and Forecasting (WRF) model (Skamarock et al.
2019) in the idealized LES mode (WRF-LES) (Wang and Feingold 2009) to simulate
the two CAO cases and test the sensitivities of the marine BL and clouds to
large-scale forcing and heat fluxes. Doubly periodic boundary conditions are
employed in horizontal directions. The horizontal resolution is set to
$dx=dy=300\,\rm{m}$ with 200 lateral grid cells, which results in a horizontal
domain size of $L_{x}=L_{y}=60\,\rm{km}$. The domain height is $z_{\rm
top}=7\,\rm{km}$ with 153 vertical $\eta$-layers. The horizontal resolution of
300 m is quite coarse for LES but it has proven to be able to simulate the
formation and evolution of cloud cellular structures in marine stratocumulus
(Wang and Feingold 2009). The time step is set to $\Delta t=3\,{\rm s}$ in all
simulations. Simulations are initiated at 06:00 UTC to allow sufficient model
spin-up time before the WRF-LES results are evaluated against measurements
taken during 16:00-17:00 UTC on February 28 and 15:00-16:00 UTC on March 1.

The two-moment Morrison cloud microphysics scheme (Morrison et al. 2009) is
used. In this part of the study, a constant number concentration of cloud
droplets derived from in-situ measurements during the ACTIVATE campaign is
prescribed in the Morrison scheme to stay focused on cloud-meteorology
interactions. Both shortwave and longwave radiative schemes are originally
from the NCAR Community Atmosphere Model (CAM 3.0), which were used in
previous WRF-LES studies, such as Wang et al. (2009) and Wang and Feingold
(2009). Surface heat fluxes and SST are all prescribed in the model as the
boundary conditions at the sea surface.

LES with horizontally uniform initial conditions and periodic boundary
conditions cannot predict changes in atmospheric state at scales larger than
its domain size. Thus, we apply relaxation to horizontal wind components $u$
and $v$ and advective tendencies to potential temperature $\theta$ and water
vapor mixing ratio $q_{v}$ as forcing terms in the prognostic equations. We
adopt the same large-scale forcing and relaxation schemes as in Endo et al.
(2015). To derive the large-scale forcings, we simplify the governing equation
of $\theta$ and $q_{v}$ by removing sink and source terms as,

$\frac{\partial\theta}{\partial
t}=-{\bm{u}}\cdot{\bm{\nabla}}\theta=-u\frac{\partial\theta}{\partial
x}-v\frac{\partial\theta}{\partial y}-w\frac{\partial\theta}{\partial z},$ (1)
$\frac{\partial q_{v}}{\partial
t}=-{\bm{u}}\cdot{\bm{\nabla}}q_{v}=-u\frac{\partial q_{v}}{\partial
x}-v\frac{\partial q_{v}}{\partial y}-w\frac{\partial q_{v}}{\partial z}.$ (2)

Applying Reynolds decomposition to Equations (1) and (2), we obtain temporal
tendencies at large scales,

$\frac{\partial\bar{\theta}}{\partial
t}=-\bar{u}\frac{\partial\bar{\theta}}{\partial
x}-\bar{v}\frac{\partial\bar{\theta}}{\partial
y}-\bar{w}\frac{\partial\bar{\theta}}{\partial z},$ (3)

and

$\frac{\partial\bar{q}_{v}}{\partial
t}=-\bar{u}\frac{\partial\bar{q}_{v}}{\partial
x}-\bar{v}\frac{\partial\bar{q}_{v}}{\partial
y}-\bar{w}\frac{\partial\bar{q}_{v}}{\partial z}.$ (4)

The large-scale horizontal advective tendencies of $\bar{\theta}$ and
$\bar{q}_{v}$ are given by the first two terms of r.h.s of Equation (3) and
Equation (4), respectively. The third term in the r.h.s of Equation (3) and
Equation (4) represents the large-scale vertical advective tendencies. These
large-scale advective tendencies are obtained from the hourly ERA5 reanalysis
data. The horizontal wind components $u$ and $v$ are applied with a relaxation
strategy (i.e., nudging LES domain-average winds to a reference state), as
also used in previous LES studies (Wang and McFarquhar 2008; Endo et al.
2015), defined by

$\left.\frac{\partial u}{\partial t}\right|_{\rm R}=\frac{u_{\rm
ERA5}-\bar{u}}{\tau},$ (5) $\left.\frac{\partial v}{\partial t}\right|_{\rm
R}=\frac{v_{\rm ERA5}-\bar{v}}{\tau},$ (6)

where $\tau$ is the relaxation time scale, which is set to be 1 hour in this
study. The subscript “R” denotes the relaxation adjustment to the horizontal
wind components.

We acknowledge that applying relaxation of wind to WRF-LES lacks physical
judgment as also addressed in Endo et al. (2015). However, LES simulations of
horizontal winds with relaxation adjustments are found to be comparable with
the reanalysis and observational data. This is not new and has been used in
the single column model (Randall and Cripe 1999) and many LES works (Neggers
et al. 2012; Heinze et al. 2017) in the meteorology community.

We also test the sensitivities of WRF-LES results to prescribed surface heat
fluxes obtained from ERA5 reanalysis data. Table 2 lists parameters examined
in the sensitivity tests.

Table 2: List of WRF-LES simulations with different forcings.  Simulation |  $\frac{\partial\bar{\theta}}{\partial t}$ $\&$ $\frac{\partial\bar{q}_{v}}{\partial t}$ |  $u\&v$ relaxation | SHF ($\rm{W\,m}^{-2}$) | LHF ($\rm{W\,m}^{-2}$) | $\bar{D}$ |  $N_{c}$ [$\rm{cm}^{-3}$]
---|---|---|---|---|---|---
0228A | Yes | Yes | 79.91 | 305.02 | Yes | 650
0301A | Yes | Yes | 231.76 | 382.18 | Yes | 450
0228B | No | Yes | 79.91 | 305.02 | Yes | 650
0228C | Yes | No | 79.91 | 305.02 | Yes | 650
0228D | No | No | 79.91 | 305.02 | Yes | 650
0301B | Yes | Yes | 231.76 | LHF(t) | Yes | 450
0301C | Yes | Yes | SHF(t) | 382.18 | Yes | 450
0301D | Yes | Yes | SHF(t) | LHF(t) | Yes | 450
0301E | Yes | Yes | 231.76 | 382.18 | No | 450

##  3 Meteorological conditions and forcings for the two cases

###  3.1 Dropsonde measurements and derived divergence

Figure 1 shows the location of individual dropsondes and the center of
dropsonde circle on an ERA5 SST map for both cases. The nearest ERA5 grid
points to the dropsondes are also shown in gray open symbols which are used to
obtain the SST for the corresponding dropsondes. Clearly, the SST is much
warmer over the circle on February 28 than March 1.

Figure 1: Location of dropsondes for the (a): February 28 case and (b): March
1 case. Black open symbols are the location of dropsondes at the surface and
the gray ones are mapped locations from the ERA5 reanalysis data. Solid dots
represent the center of the dropsonde circle. The contour map shows ERA5 SST
in the measurement region. Black stars represent the location of Hampton, VA
on this map. The yellow curve shows the flight path.

Figure 2 shows the vertical profiles of RH, $q_{v}$, $\theta$, $u$, and $v$
from dropsonde measurements for the two cases. The February 28 case (RF10) is
characterized by a deeper boundary layer with a depth of about $2.8\,{\rm km}$
and a drier free troposphere compared to the March 1 case (RF13). Individual
RH and $q_{v}$ profiles show more fluctuations from the mean in the free
troposphere on March 1 than the February 28 case. The boundary layer for the
March 1 case is shallower. The magnitude of $u$ and $v$ increases rapidly with
height above the boundary layer, which is more profound on February 28,
showing a strong wind shear. The meteorological states evolve substantially
during the one-hour sampling time period of both cases, as indicated by the
contrast between the first dropsonde (blue curve) and the last one (red curve)
that were released roughly at the same location. The boundary layer became
deeper (shallower) with time on February 28 (March 1).

Figure 2: Profiles from dropsonde measurements for the February 28 case (upper
panel) and March 1 case (lower panel). The gray lines represent vertical
profiles measured from 11 dropsondes and the thick black lines represent the
corresponding mean profile. The blue and red curves represent the first and
last dropsonde, respectively, released at about the same location but one hour
apart.

The vertical velocity of airflow regulates the atmospheric water distribution
but is difficult to measure (Bony and Stevens 2019). The continuity equation
of nearly incompressible airflow with velocity $\bm{u}=\bm{u}(u,v,w)$ is given
by

${\bm{\nabla}}\cdot{\bm{u}}=\frac{\partial u}{\partial x}+\frac{\partial
v}{\partial y}+\frac{\partial w}{\partial z}=0.$ (7)

According to Equation (7), the vertical velocity $w$ can be expressed as

$w(z)=-\int_{0}^{z}\left(\frac{\partial u}{\partial x}+\frac{\partial
v}{\partial y}\right){\rm d}z.$ (8)

Divergence is defined as

$D=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}.$ (9)

Thus, Equation (9) can be written as

$w(z)=-\int_{0}^{z}D{\rm d}z.$ (10)

Therefore, $w$ can be indirectly obtained from the measured horizontal wind
components. In the atmospheric boundary layer, motion of airflow is
conventionally decomposed to large and small scales. By applying Reynolds
decomposition to Equation (7), we obtain the large-scale divergence,

$\bar{D}=\frac{\partial\bar{u}}{\partial x}+\frac{\partial\bar{v}}{\partial
y}.$ (11)

We follow the procedure described by Lenschow et al. (2007) to calculate the
divergence from dropsonde measurements, details of which are given in Appendix
A. Since we use the linear regression method to estimate $\bar{D}$, the
standard error $\sigma_{\epsilon}$ can be estimated. The large-scale
divergence $\bar{D}$ enters the third term in the r.h.s of Equation (3) and
Equation (4) through Equation (10). Thus the effect of $\bar{D}$ is taken into
account via vertical motion, $w$, in the vertical component of the advective
tendency.

We compare $\bar{D}$ estimated from the dropsonde measurements with the one
obtained from ERA5 reanalysis data for both cases as shown in Figure 3.
ERA5-$\bar{D}$ profile (red curve) averaged between 16:00 UTC and 17:00 UTC
(two vertical profiles) at the dropsonde center is able to capture the sign of
$\bar{D}$ vertical-structure measured by dropsonde (black curve) for the
February 28 case (Figure 3(a)) within 1 km above the surface. However, it
differs from the dropsonde measurements in both signs and magnitude above 1
km. More interestingly, ERA5-$\bar{D}$ profile evolves from convergence to
divergence within the boundary layer from 16:00 UTC (blue curve) to 17:00 UTC
(cyan curve). We also examine the ERA5 $\bar{D}$ averaged between 16:00 UTC
and 17:00 UTC at the location of each dropsonde as shown by coral-colored
dashed lines in Figure 3, which exhibits large spatial fluctuations. This
demonstrates strong spatial inhomogeneity in large-scale vertical motions over
this area on February 28, as also indicated by the sensitivity of $\bar{D}$ to
different dropsonde subsets as shown in Figure A2. The amplitude of the mean
and largest value of dropsonde-$\bar{D}$ for the February 28 case is
$\langle|\bar{D}|\rangle=1.38\times 10^{-5}\,{\rm s}^{-1}$ and $|\bar{D}|_{\rm
max}=4.06\times 10^{-5}\,{\rm s}^{-1}$, respectively. For the March 1, they
are $\langle|\bar{D}|\rangle=2.99\times 10^{-5}\,{\rm s}^{-1}$ and
$|\bar{D}|_{\rm max}=7.09\times 10^{-5}\,{\rm s}^{-1}$. These values are about
one order of magnitude larger than the ones from non-CAO marine BL cloud
regimes. The mean value of $\bar{D}$ inferred from the ensemble of radiosondes
during the Atlantic Trade-Wind Experiment (Augstein et al. 1973) and the
Barbados Oceanographic and Meteorological Experiment (Holland and Rasmusson
1973) is about $10^{-6}\,{\rm s}^{-1}$. Similar values were used in the case
studies during the VAMOS Ocean-Cloud-Atmosphere-Land Study (Rahn and Garreaud
2010; Wang et al. 2010) and Second Dynamics and Chemistry of Marine
Stratocumulus field study (DYCOMS II) field campaign (Wang and Feingold 2009).

Figure 3: Comparison of $\bar{D}$ profiles estimated from dropsonde
measurements with the one from ERA5 reanalysis data for the (a): February 28
case and (b): March 1 case. Black dotted lines represent $\bar{D}$ estimated
from 10 dropsondes with $\pm\sigma_{\epsilon}$ uncertainty (gray shaded area).
Blue lines represent ERA5-$\bar{D}$ profiles at 16:00 UTC for the February 28
case and 17:00 UTC for the March 1 case at the center of dropsondes. Cyan
lines represent ERA5-$\bar{D}$ profiles at 17:00 UTC for the February 28 case
and 16:00 UTC for the March 1 case at the center of dropsondes. Red solid-
dotted lines represent $\bar{D}$ from ERA5 reanalysis data averaged during the
measurement time (between blue lines and cyan lines). The dashed coral-colored
lines represent ERA5 reanalysis data at the location of individual dropsondes
averaged during the measurement time for each case.

Figure 3(b) shows the same comparison but with the dropsonde measurements
conducted on March 1. In this case, $\bar{D}$ obtained from ERA5 reanalysis
data (red solid-dotted line) at the center of dropsonde circle is able to
capture the general vertical structure of $\bar{D}$ estimated from dropsonde
measurements. Similar to the February 28 case, there is a strong
spatiotemporal variation in the ERA5 $\bar{D}$.

We further compare the large-scale vertical velocity $w$ (subsidence) with
MERRA-2 reanalysis data as shown in Figure 4(e) for both cases. For the
February 28 case, the $w$ profile from MERRA-2 is averaged between 15:00 UTC
and 18:00 UTC and the one from ERA5 is averaged between 16:00 UTC and 17:00
UTC to better match the dropsonde sampling time. Both the ERA5 and MERRA-2 can
reasonably capture the vertical profile of $w$ when compared with dropsonde
measurements for this case. For the March 1 case, the $w$ profile from MERRA-2
reanalysis data at 15:00 UTC is used to compare with dropsonde measurements
while the one from ERA5 reanalysis data is averaged between 15:00 UTC and
16:00 UTC. The ERA5 reanalysis data agree with the dropsonde measurements in
the sign but underestimate the magnitude. The MERRA-2 does not capture the
structure and magnitude of the vertical profile of $w$ well. Comparison of
$\theta$, $q_{v}$, $u$, and $v$ profiles is also shown in Figure 4(a)–(d).
MERRA-2 shows a slightly warmer boundary layer for the February 28 case while
the ERA5 shows a colder one. Both MERRA-2 and ERA5 data capture the $\theta$
profile well for the March 1 case. ERA5 yields a drier ($q_{v}$ profiles)
boundary layer while MERRA-2 capture the $q_{v}$ well compared to the
dropsonde measurements for both cases. The $u$ and $v$ profiles within the
boundary layer are represented well by MERRA-2 and ERA5 data for both cases,
given the large spread among the individual dropsondes for the circled area
(see Figure 2). The ERA5 captures those profiles above the boundary layer
better than the MERRA-2. Overall, comparing to MERRA-2, ERA5 profiles are more
consistent with the dropsonde measurements, as also shown in Chellappan et al.
(2021) for the broader WNAO region.

Since we aim to use the divergence as part of the large-scale forcings to
drive WRF-LES, the agreement of $\bar{D}$ (and the corresponding $w$) from
ERA5 with the estimates from dropsondes for the March 1 case affords
confidence to use the hourly ERA5 divergence to test the sensitivity of WRF-
LES to time-varying large-scale forcings.

Figure 4: Comparison among dropsonde measurements (black curves), ERA5 (red
curves), and MERRA-2 (blue curves) at the dropsonde center for both cases.
Solid symbols represent the February 28 case (ERA5 data averaged between 16:00
and 17:00 UTC and MERRA-2 averaged between 15:00 and 18:00 UTC) and the open
circles mark the March 1 case (ERA5 averaged between 15:00 and 16:00 UTC and
MERRA-2 at 15:00 UTC). Large-scale vertical velocity $w$ from ERA5 corresponds
to ERA5-$\bar{D}$ in Figure 3.

###  3.2 Surface heat fluxes

Turbulent sensible and latent heat fluxes at the surface are important
forcings to drive LES of boundary layer clouds. They are responsible for the
heat and moisture exchange between the ocean and atmosphere. Surface heat
fluxes are challenging to measure and estimate due to the nonlinear processes
involved. Therefore the so-called bulk aerodynamics parameterization has been
used to estimate surface heat fluxes. Bulk aerodynamics algorithms
parameterize the turbulence instability as well as the roughness length of the
wind speed, temperature and the water vapor mixing ratio (Zeng et al. 1998).
The surface sensible heat flux (SHF) and latent heat flux (LHF) are given by
(Smith 1988)

${\rm SHF}=C_{T}\rho c_{\rm p}U(T_{\rm s}-\theta),$ (12) ${\rm
LHF}=C_{E}L_{v}\rho U(q_{\rm s}-q),$ (13)

where $\rho$ is the air mass density, $c_{\rm p}$ is the specific heat of air,
$U$ is the mean wind speed at a particular reference height, $T_{\rm s}$ is
the SST, $\theta$ is the potential temperature at the reference height,
$q_{\rm s}$ is the saturated water vapor mixing ratio, $q$ is the water vapor
mixing ratio, and $L_{v}$ is the latent heat of evaporation. The sensible heat
flux coefficient $C_{T}$ and evaporation coefficient $C_{E}$ are determined by
the empirical Monin–Obukhov similarity theory. In this study, we adopt the
bulk aerodynamics algorithms developed by Zeng et al. (1998) to calculate heat
fluxes for comparison between ERA5 and dropsonde measurements, which we refer
to as “Z98” hereafter. The input parameters of the Z98 algorithm are SST, 10-m
wind speed $U_{10\rm m}$, temperature $T_{10\rm m}$, and the water vapor
mixing ratio $q_{v,10\rm m}$ as described by Equation (12) and Equation (13).

To derive the surface heat fluxes using dropsonde measurements and evaluate
ERA5 reanalysis data, we first compare the ERA5 SST to satellite measurement
for the two cases as a quality check. As shown in Figure 5, for both 28
February 2020 (black symbols) and 1 March 2020 cases (red symbols), SST from
ERA5 and the satellite measurement matches well at the center of dropsonde
circle. At the location of individual dropsondes the agreement is reasonable
on February 28 while several points on March 1 are quite off, which is likely
because of the location mismatch due to the resolution difference and sampling
area being near strong SST gradients (shown in Figure 1). The Normalized Root-
Mean-Square Error (NRMSE) is $0.1\%$ for the February 28 case and is $0.6\%$
for the March 1 case. This comparison suggests that the SST from ERA5 can also
be used as the initial input for our WRF-LES. ERA5 has assimilated the
Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) system for
hourly SST starting from 2007 (Hirahara et al. 2016). The OSTIA assimilated
the MW-IR measurements. Thus, the agreement between ERA5 and the satellite
retrievals is expected.

Figure 5:  Comparison of SST from satellite retrievals and ERA5 reanalysis
data for both February 28 (black) and March 1 (red) cases . Open symbols
represent SST at the position of each dropsonde. The two solid dots represent
SST at the center of dropsonde circles.

Next, we compare surface heat fluxes directly obtained from ERA5 reanalysis
data and the ones estimated from ACTIVATE measurements. Since the SST obtained
from ERA5 agrees with the satellite measurement, we try to examine if the ERA5
heat-fluxes can be reproduced from ERA5 SST and dropsonde measurements. First,
we use the Z98 algorithm to calculate heat fluxes based on SST, $U_{\rm
10m}^{\rm ERA5}$, $q_{v,{\rm 10m}}^{\rm ERA5}$ from ERA5 and $T_{\rm 10m}^{\rm
drop}$ from dropsonde. Figure 6(a) shows the comparison between the estimated
SHF and the ERA5 reanalysis data (red and black circles). The corresponding
NRMSE is $5.1\%$ and is $36.4\%$ for the March 1 and February 28 case,
respectively. The comparison of LHF is shown in Figure 6(b) with NRMSE of
$8.5\%$ and $30.0\%$ for the March 1 (red circles) and February 28 case (black
circles), respectively. These demonstrate a good agreement between the
estimated heat fluxes and the ERA5 reanalysis data on March 1 case, given the
large spread within the dropsonde circle. The agreement is particularly good
at the circle center (solid symbols). This suggests that the Z98 algorithm can
be used to calculate heat fluxes and the use of $T_{\rm 10m}^{\rm drop}$ is
justified in this study. It is evident that SHF and LHF calculated using Z98
are underestimated compared to ERA5 for the February 28 case. This is because
ERA5-$\theta$ is smaller than the dropsonde-$\theta$ within the boundary
layers as shown in Figure 4(a). We then use $U_{\rm 10m}^{\rm drop}$ and
$q_{v,{\rm 10m}}^{\rm drop}$ obtained from dropsonde measuremens to estimate
the fluxes (stars in Figure 6), which yields a NRMSE value of $11.5\%$ for SHF
and $19.3\%$ for LHF for the March 1 case, further indicating that ERA5 gives
a good estimate of turbulent heat fluxes for the March 1 case. For the
February 28 case, NRMSE of SHF and LHF calculated using $U_{\rm 10m}^{\rm
drop}$ and $q_{v,{\rm 10m}}^{\rm drop}$ are $33.2\%$ and $23.1\%$,
respectively. This underestimation is because that $q_{v}$ from ERA5 is
smaller than the one from dropsonde measurements within the boundary layer as
shown in Figure 4(b). We also compare heat fluxes between the ERA5 and MERRA-2
reanalysis data. Both SHF and LHF agree well between MERRA-2 and ERA5 for the
February 28 case. However, MERRA-2 underestimates SHF and LHF compared to ERA5
for the March 1 case (see appendix B).

Overall, by adopting the ERA5 SST and $U_{\rm 10m}$, $q_{\rm v,10m}$, and
$T_{\rm 10m}^{\rm drop}$ from dropsonde measurements to estimate heat fluxes
using Z98 algorithm, we are able to evaluate the heat fluxes from ERA5. The
time-varying ERA5 heat fluxes are then used in the WRF-LES sensitivity tests.

Figure 6: Comparison between the estimated heat fluxes and the one from ERA5
reanalysis data at 16:00 UTC for the February 28 case (black symbols) and
15:00 UTC for the March 1 case (red symbols). SST from ERA5 is adopted for the
calculation. Open symbols represent SST at the position of individual
dropsondes and solid symbols represent the center of the dropsonde circle. Two
sets of input variables are adopted to calculate the surface heat fluxes:
$U_{10\rm m}^{\rm ERA5}$ and $q_{v,10\rm m}^{\rm ERA5}$ (circles) and
$U_{10\rm m}^{\rm drop}$ and $q_{v,10\rm m}^{\rm drop}$ (stars).

##  4 WRF-LES sensitivities to large-scale forcings and contrast between the
two CAO cases

###  4.1 Sensitivities to large-scale advective tendencies and relaxation

In this section, we investigate how to better represent time-varying
meteorological states in idealized WRF-LES simulation applying either
advective tendencies to $\theta$ and $q_{v}$, relaxation to $u$ and $v$, or
both. Simulations are driven by constant surface fluxes ${\rm SHF}(t_{0})$ and
${\rm LHF}(t_{0})$. Here $t_{0}$ denotes the starting time of simulations.
Since we have shown in the previous section that ERA5 reanalysis data agree
well with the dropsonde measurements during the sampling time periods of the
two CAO cases, we adopt hourly $\theta$, $q_{v}$, $u$, and $v$ vertical
profiles from ERA5 reanalysis data and derive the corresponding vertical
profiles of advective tendencies and relaxation adjustments. The hourly
meteorological states simulated in WRF-LES are then validated against ERA5
reanalysis data and dropsonde measurements.

Figure 7 shows the hourly (rainbow-colored lines) input meteorological forcing
being obtained from ERA5 reanalysis data for the WRF-LES simulations. The
evolution of vertical profiles are averaged over a $2^{\circ}\times 2^{\circ}$
area centered at the middle of the dropsonde circle of each case. This
selected area sufficiently covers the dropsonde circle. Vertical profiles of
$\theta$, $q_{v}$, $u$, $v$, and $w$ obtained from ERA5 reanalysis data
averaged during the measurement time (black solid lines) agree reasonably well
with the dropsonde measurements (gray dashed lines) for the February 28 case
(upper row) and for the March 1 case (lower row). Vertical profiles of
advective tendencies of $\theta$ and $q_{v}$ (i.e.,
$\partial\bar{\theta}/\partial t$ and $\partial\bar{q}_{v}/\partial t$) are
calculated from $\theta_{\rm ERA5}$ and $q_{v,\rm ERA5}$. Vertical profiles of
$\theta_{\rm ERA5}$, $q_{v,\rm ERA5}$, $u_{\rm ERA5}$, and $v_{\rm ERA5}$ at
06:00 UTC are taken as the input sounding for WRF-LES simulations when the
simulation starts. We note that $w_{\rm ERA5}$ averaged over the
$2^{\circ}\times 2^{\circ}$ area at 15:00 UTC differs slightly from the one at
the dropsonde center shown by the red curve in Figure 4 as expected. This is
because of the strong spatial variation of $\bar{D}$ as shown in Figure 3. We
have tested the relaxation time scale $\tau$ of $u$ and $v$ for the February
28 case and found that WRF-LES simulations with $\tau$=30 min, 1 h, and 3 h
reveal almost identical vertical profiles and liquid water path (LWP).
Therefore we adopt $\tau$=1 h for all the simulations as the ERA5 reanalysis
data has a time resolution of one hour. The relaxation is applied to all
vertical layers of the LES domain.

Figure 7: Hourly meteorological state and forcing profiles for the February 28
case (simulation 0228A, upper panel) and March 1 case (simulation 0301A, lower
panel) from ERA5 reanalysis data averaged over a $2^{\circ}\times 2^{\circ}$
domain. The rainbow color scheme represents the time evolution (06:00-21:00
UTC): from purple to red. The averaged ERA5 reanalysis data over the
measurement time period are marked by black lines, which are compared with the
dropsonde measurements (dashed gray lines).

We first perform a simulation without applying advective tendencies of
$\theta$ and $q_{v}$ and relaxation of $u$ and $v$ (simulation 0228D) for the
February 28 case. It is shown by the blue curves in Figure 8 that such a
configuration yields vertical profiles that have a large deviation from the
ERA5 reanalysis data (cyan curves) and dropsonde measurements (grey curves).
The $\theta$ profile from WRF-LES differs considerably from ERA5 above the
boundary layer and the $q_{v}$ profile shows a more humid boundary layer than
the ERA5. The $u$ and $v$ profiles from WRF-LES deviate from the ERA5 and
dropsonde measurements. When $\partial\bar{\theta}/\partial t$ and
$\partial\bar{q}_{v}/\partial t$ are applied (simulation 0228C), $\theta$ and
$q_{v}$ profiles from WRF-LES agree well with the dropsonde measurements as
shown by the red curves in Figure 8. However, $u$ and $v$ profiles still
deviate from the ERA5 reanalysis data. We then only apply the $u$ and $v$
relaxation to the WRF-LES (simulation 0228B). As shown by the green curves of
Figure 8, $u$ and $v$ profiles from WRF-LES are in good agreement with
dropsonde measurements even though $\theta$ and $q_{v}$ profiles differ from
the measurements. This naturally leads to the configuration of applying
advective tendencies of $\theta$ and $q_{v}$ together with $u$ and $v$
relaxation to $u_{\rm ERA5}$ and $v_{\rm ERA5}$ (simulation 0228A). Such a
configuration leads to vertical profiles ($\theta$, $q_{v}$, $u$, and $v$)
that are comparable to ERA5 reanalysis data and dropsonde measurements as
shown by the black curves of Figure 8. Therefore, this combined forcing and
relaxation scheme is justified to simulate the two CAO cases. Evolution of the
vertical profiles for simulations 0228A, 0228B, 0228C, and 0228D are shown in
section C. As also shown in Figure 8, the magnitude and time evolution of
cloud water simulated by WRF-LES are sensitive to the boundary layer
meteorological conditions.

Figure 8: Domain-averaged vertical profiles for simulation 0228A (black
curve), 0228B (green curve), 0228C (red curve), and 0228D (blue curve) with
the corresponding input forcings shown in Figure 7 for the February 28 case
during the measurement time shown in Table 1. The cyan curve represent
profiles from the ERA5 reanalysis data averaged during the measurement time
and the grey curve represent the ones from dropsonde measurements. Heat fluxes
are from ERA5 reanalysis data: ${\rm SHF}(t_{0})=79.91\,\rm{W\,m}^{-2}$ and
${\rm LHF}(t_{0})=305.02\,\rm{W\,m}^{-2}$.

Applying $\partial\bar{\theta}/\partial t$ and $\partial\bar{q}_{v}/\partial
t$ leads to a colder and less humid boundary layer (compare simulation 0228C
and 0228D) but allows the boundary layer to grow higher. This results in a
deeper cloud layer with a reduced amount of the liquid water content $q_{c}$,
which is enhanced by a factor of two by applying $u$ and $v$ relaxation. We
then examine the time evolution of LWP. LWP peaks around 16:00 UTC and then
starts to decrease as shown in Figure 9, which is very likely due to the solar
heating. To validate the simulated LWP, we compare it with the RSP retrievals
during the ACTIVATE field campaign. As shown in Figure 10(a), the WRF-LES
simulation (averaged over the measurement time) agrees reasonably well with
the RSP measurement. This further illustrates that the WRF-LES is able to
capture the cloud formation and evolution in this case study. We also tested
the tendencies and relaxation forcing for the March 1 case, which yields the
same conclusion as for the February 28 case. Figure 10(b) shows that the WRF-
LES underestimates the frequency of lower LWP values (less than $100\,{\rm
g\,m}^{-2}$) but overestimates the frequency between $200-400\,{\rm
g\,m}^{-2}$ for the March 1 case.

Figure 9: Time series of domain-averaged LWP in the simulations with different
forcing options as shown in Figure 8.  Figure 10: Comparison of PDF of LWP
between RSP measurements and WRF-LES for simulation (a): 0228A and (b): 0301A.
The PDF for WRF-LES LWP is calculated from 3 snapshots during the measurement
time. The minimum value of LWP from the RSP measurement is taken as a lower
cutoff for the simulated LWP.

###  4.2 Sensitivities to large-scale divergence $\bar{D}$

We have shown that applying $\partial\bar{\theta}/\partial t$,
$\partial\bar{q}_{v}/\partial t$, $u$ relaxation to $u_{\rm ERA5}$, and $v$
relaxation to $v_{\rm ERA5}$ is essential to reproduce the time evolution of
meteorological states for the two CAO cases considered in this study. This
configuration is adopted to further test the sensitivities of WRF-LES results
to $\bar{D}$. We focus on the March 1 since it is more challenging to simulate
due to the large surface heat fluxes.

We perform two WRF-LES simulations with or without the large-scale vertical
velocity as a forcing (third term on the r.h.s of Equations (3) and (4)) for
the March 1 case. The forcing configuration for the baseline simulation 0301A
is the same as simulation 0228A. To examine the impact of large-scale
divergence separately, we conduct a simulation (0301E in Table 2) that
excludes the forcing term related to $\bar{D}$ (vertical component of the
advective tendencies) but keeps all other forcings the same as the baseline
(0301A). As shown in Figure 11, the initial model spin-up time (06:00-08:00
UTC) is characterized by a sharp increase of LWP from 0 to about
$200\,\rm{g\,m}^{-2}$ for the baseline simulation (blue). Without the time
varying large-scale divergence (cyan), LWP and IWP experience a larger
increase compared to the baseline (blue) during the initial time steps as
expected. LWP from simulation 0301E decreases and becomes smaller compared to
the one from simulation 0301A. This is because IWP from 0301E are about 10
times more than the ones from 0301A. In the absence of $\bar{D}$, the updraft
cooling is more profound. Thus, the ice formation is enhanced. The cyan curve
in Figure 12 shows the deviation of vertical profiles of the simulation 0301E
to 0301A averaged over the measurement time. As expected from a positive $w$,
a much deeper boundary layer is developed in the simulation without
subsidence. The moist air is mixed to over 3 km. This can be explained by the
fact that the lack of divergence breaks the balance between the BL growth
driven by surface fluxes and suppression due to subsidence. Figure 13 shows
the contribution of horizontal and vertical advective tendencies to
$\partial\theta/\partial t$ (upper row) and $\partial q_{v}/\partial t$ (lower
row) for simulation 0301A and 0301E, respectively. When $\bar{D}=0$
(simulation 0301E), only horizontal advective tendencies contribute to
$\partial\theta/\partial t$ by comparing the solid cyan curves in Figure
13(a)–(c). When $\bar{D}\neq 0$, it is evident that the vertical advective
term ($\bar{D}$) dominates the temperature and humidity changes due to large-
scale tendencies, especially near the inversion layer. The contribution from
horizontal terms are small when comparing simulations 0301A and 0301E. The
vertical profiles of $\partial\theta/\partial t$ and $\partial q_{v}/\partial
t$ only evolve slightly from 15:00 to 16:00 UTC because $\bar{D}$ does not
vary much for the March 1 case. Recall that for the February 28 case,
$\bar{D}$ evolves from convergence to divergence from 16:00 to 17:00 UTC
(Figure 3(a)). Consequently, the vertical profiles of $\partial\theta/\partial
t$ and $\partial q_{v}/\partial t$ are nealry zero as shown in Figure 13. When
comparing vertical profiles of $\partial\theta/\partial t$ and $\partial
q_{v}/\partial t$ for the March 1 (blue lines) and February 28 (black lines)
cases, we see that $\bar{D}$ has a more profound impact for the March 1 case.

Figure 11:  Time series of domain-averaged (a): LWP (in-cloud liquid water and
rain) and (b): IWP (ice, graupel, and snow) from WRF-LES simulations with
different heat fluxes and large-scale divergence as indicated by the legends
for the March 1 case. Heat fluxes are from ERA5 reanalysis data, the values of
which are $\rm{SHF}(t_{0})=231.76\,\rm{W\,m}^{-2}$ and
$\rm{LHF}(t_{0})=382.18\,\rm{W\,m}^{-2}$. Here $t_{0}$ denotes the starting
time of the WRF-LES simulation, which is 06:00 UTC.  Figure 12: Deviation of
vertical profiles of simulation 0301B, 0301C, 0301D, and 0301E from the
baseline simulation 0301A averaged over the measurement time (3 snapshots over
15:00 to 16:00 UTC).  Figure 13: Contributions of large-scale forcings to
$\partial\theta/\partial t$ (upper row) and $\partial q_{v}/\partial t$ (lower
row) for simulation 0301A (blue curves) and 0301E (cyan curves) at 15:00 UTC
(solid lines) and 16:00 UTC (dashed lines). Black lines represent the ones for
simulation 0228A at 16:00 UTC (solid line) and 17:00 UTC (dashed line). LS,
LSHOR, and LSVER denote large-scale forcing due to total, horizontal, and
vertical advective tendencies, respectively.

To conclude, $\bar{D}$ likely has a strong control on the spin-up and time
evolution of the boundary layer in WRF-LES simulations of the two CAO events
we explored here, especially for the March 1 case. Therefore, including a
time-varying $\bar{D}$ profile to drive WRF-LES is necessary for simulating
the fast-evolving CAO events over WNAO.

###  4.3 Sensitivities to surface heat fluxes

To test the sensitivities to surface fluxes, we perform another three WRF-LES
simulations using the same forcing configuration as the baseline simulation
0301A but with temporally varying and spatially uniform surface heat fluxes
for the March 1 case. The time series of such surface heat fluxes obtained
from ERA5 reanalysis data at the center of dropsonde circle is shown in Figure
14. When the WRF-LES simulation is forced by ${\rm SHF}(t_{0})$ and ${\rm
LHF}(t)$ (red curve), LWP evolves in the same pattern as the baseline but with
larger values between 10:00-16:00 UTC as shown in Figure 11. This is because
${\rm LHF}(t)$ is larger than ${\rm LHF}(t_{0})$ until 15:00 UTC, as shown in
Figure 14. Overall, simulations driven by ${\rm LHF}(t)$ result in more LWP
compared with the one by ${\rm LHF}(t_{0})$. Simulations forced by ${\rm
SHF}(t)$ and ${\rm LHF}(t_{0})$ (green curve) exhibit the same trend as the
one by ${\rm SHF}(t_{0})$ and ${\rm LHF}(t)$. When the time-varying ${\rm
SHF}(t)$ and ${\rm LHF}(t)$ are both applied to the WRF-LES simulation (black
curve), the initial increase in ${\rm SHF}(t)$, as compared to ${\rm
SHF}(t_{0})$ (red), does not have an impact on the LWP. Overall, the LWP
simulated by WRF-LES evolves monotonically with the evolution of ${\rm
SHF}(t)$ and ${\rm LHF}(t)$ in the parameter range we have explored. Since the
forcing ${\rm SHF}(t)$ and ${\rm LHF}(t)$ only vary slightly, the mean LWP
values do not show a significant difference when comparing the four WRF-LES
simulations. We also compare the Ice Water Path (IWP) as shown in Figure
11(b). The evolution of these quantities follow the same trend as LWP. Figure
12 shows the corresponding deviations of vertical profiles of simulation
0301B, 0301C, 0301D, and 0301E from the baseline simulation 0301A. These
profiles are averaged over the measurement time (3 snapshots over 15:00 to
16:00 UTC). Differences at the inversion layer (about 2 km) are the most
pronounced. The green curves (${\rm SHF}(t)$, ${\rm LHF}(t_{0})$) deviate the
least from the blue curves (baseline simulation) while the red (${\rm
SHF}(t_{0})$, ${\rm LHF}(t)$) and black (${\rm SHF}(t)$, ${\rm LHF}(t)$)
curves diverge the most within the boundary layer . The red and black curves
are almost identical except for the slight difference in $q_{c}$. This
indicates that the effect of SHF variations on LWP is not obvious.

Figure 14: Surface heat fluxes from ERA5 reanalysis data at the center of
dropsonde circle on February 28 (black) and March 1 (red). Solid dots and
squares represent heat fluxes at the center of dropsonde circles calculated
based on the Z98 algorithm.

##  5 Discussion and conclusions

We have reported two contrasting cold air outbreak (CAO) cases observed during
the ACTIVATE field campaign and the corresponding WRF-LES modeling of them.
The February 28 case is characterized by weaker turbulent surface heat fluxes
(${\rm SHF=79.91\,W\,m}^{-2}$ and ${\rm LHF=305.02\,W\,m}^{-2}$) than those of
the March 1 case (${\rm SHF=231.76\,W\,m}^{-2}$ and ${\rm
LHF=382.18\,W\,m}^{-2}$). The divergence is on the order of $10^{-5}\,{\rm
s}^{-1}$ for both cases, which is about 10 times larger than common marine
cases (e.g., $\bar{D}\approx 10^{-6}\,{\rm s}^{-1}$ in the DYCOMS-II case
simulated by Wang and Feingold (2009)) and about two times larger than the CAO
case in de Roode et al. (2019). A deeper, warmer, and more humid boundary
layer was observed for the February 28 event than the one on March 1.

To examine and validate different prescribed forcing options to drive WRF-LES,
we first evaluate divergence obtained from the ERA5 reanalysis data against
the one derived from dropsonde measurements for the two CAO cases. The
divergence profile and the corresponding vertical velocity obtained from ERA5
reanalysis data at the center of dropsonde circle are able to capture the
structure of the ones estimated from dropsonde measurements for the March 1
case. This gives us the confidence to adopt the time-varying divergence
profiles from ERA5 to drive our WRF-LES simulations.

Since the surface turbulent heat fluxes are partly determined by SST, we
compare SST from ERA5 to the one from satellite retrievals. They agree very
well for both February 28 and March 1 cases. Therefore, SST from ERA5 together
with the 10 m temperature, water vapor mixing ratio, and wind speed from
dropsonde measurements are used to calculate heat fluxes for the March 1 case
and those from ERA5 for the February 28 case using the bulk aerodynamic
algorithms from Zeng et al. (1998). The estimated sensible and latent heat
fluxes agree well with the ones directly obtained from ERA5 reanalysis data
for the March 1 case case. They are underestimated by about $30\%$ compared to
the ERA5 heat fluxes for the February 28 case.

By applying the surface heat fluxes, large-scale temperature and moisture
advective tendencies, and wind relaxation adjustments from ERA5 to the WRF-LES
simulations, the simulated meteorological states for both CAO cases match the
ERA5 reanalysis data and the ACTIVATE field campaign measurements. We also
conduct WRF-LES sensitivity simulations on the surface fluxes and divergence
and find that the divergence is important in suppressing the evolution of the
boundary layer and achieves the observed states of the boundary layer for this
case, while surface heat fluxes are more influential for the simulated LWP.
The PDF of LWP produced from our WRF-LES simulations agrees reasonably well
with the measured ones from the ACTIVATE campaign for both the February 28
case. Since the large-scale tendencies profiles vary with time for the two CAO
cases, it is important to apply time-varying tendencies to the WRF-LES instead
of constant ones.

In summary, with initial conditions, large-scale forcings, and turbulent
surface heat fluxes obtained from ERA5 and validated by ACTIVATE airborne
measurements, WRF-LES is able to reproduce the observed boundary-layer
meteorological states and LWP for two contrasting CAO cases. This manifests
the meteorological impact on marine boundary layer and clouds associated with
CAO over WNAO. This study (Part 1) paves the path to further investigation of
aerosol effects on cloud microphysics during the CAO events to be reported in
the forthcoming companion paper (Part 2).

###### Acknowledgements.

This work was supported through the ACTIVATE Earth Venture Suborbital-3
(EVS-3) investigation, which is funded by NASA’s Earth Science Division and
managed through the Earth System Science Pathfinder Program Office. The
Pacific Northwest National Laboratory (PNNL) is operated for the U.S.
Department of Energy by Battelle Memorial Institute under contract DE-
AC05-76RLO1830. We wish to thank the pilots and aircraft maintenance personnel
of NASA Langley Research Services Directorate for their work in conducting the
ACTIVATE flights. We thank Andrew S Ackerman for discussions. The source code
used for the simulations of this study, the Weather Research and Forecasting
(WRF) model, is freely available on https://github.com/wrf-model/WRF. The
simulations were performed using resources available through Research
Computing at PNNL. ACTIVATE data are publicly available at: https://www-
air.larc.nasa.gov/cgi-bin/ArcView/activate.2019 MW-IR SST are produced by
Remote Sensing Systems and sponsored by NASA. Data are available at
www.remss.com

##  Appendix A Dropsonde measurements

This appendix is to review the method being adopted to calculate divergence
$D$ from the dropsonde measurements and to test statistical convergence of $D$
to the number of dropsondes used in the calculation.

The integral form of Equation (7) is

$\int{\bm{\nabla}}\cdot{\bm{u}}{\rm d}S=0.$ (14)

Thus,

$\int\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial
y}\right){\rm d}S=-\int\frac{\partial w}{\partial z}{\rm d}S=-A\frac{\partial
w}{\partial z}.$ (15)

According to the Stokes theorem, Equation (15) can be written as

$\frac{\partial w}{\partial z}=-D=-\frac{1}{A}\oint_{l}v_{n}{\rm d}l.$ (16)

Based on Equation (16), Lenschow et al. (2007) noted that the most efficient
flight track is a circle-like shape since the circle has the largest enclosed
area of any closed curve and the turning rate of flight is slow. Such a method
to calculate $D$ is called the “linear integral method”. This method requires
a closed and circular flight track and a linearly evolving wind speed, none of
which can be satisfied. Therefore, Lenschow et al. (2007) developed the
“regression method”, which can alleviate the requirements of the linear
integral method. The first order Taylor expansion of horizontal wind velocity
${\bm{v}}={\bm{v}}(u,v)$ at the dropsondes center is

${\bm{v}}={\bm{v}_{0}}+\frac{\partial{\bm{v}}}{\partial x}\Delta
x+\frac{\partial{\bm{v}}}{\partial y}\Delta y+\frac{\partial{\bm{v}}}{\partial
t}\Delta t,$ (17)

where $\Delta x$ and $\Delta y$ are the eastward and northwest displacements
from the center of dropsondes. Assuming a stationary state, the term
$\frac{\partial{\bm{v}}}{\partial t}\Delta t$ can be neglected. This
assumption suggests that all the dropsondes should be released simultaneously,
which is not feasible experimentally. The sampling lag in space and time
between different dropsondes may cause error of calculating $\bar{D}$.
However, Bony and Stevens (2019) demonstrated that the stationarity assumption
is not bad.

We then test the sensitivity of $\bar{D}$ to the number of dropsondes used in
the calculation. Figure A1 shows the circular distribution of 10 dropsondes
(black solid dots) for the February 28 case. Here, only 10 dropsondes are used
to calculate the divergence because two (out of the 11 dropsondes) were
released at the same location. We select two subsets as shown in Figure A1(a).
The corresponding $\bar{D}$ is shown in Figure A2(a). $\bar{D}$ derived from
the 5-dropsonde circle agrees with that from the 10-dropsonde circle. We also
test four subsets of 4-dropsonde circles, as shown in Figure A1(b), and the
corresponding $\bar{D}$ is shown in Figure A2(b). Same as the 5-dropsonde
circles, 4-dropsonde subset agree with the 10-dropsonde circle even though the
individual subsets exhibit differences.

We also apply the same analysis to dropsonde measurements being carried out on
the March 1 case as shown in Figures A3 and A4. $\bar{D}$ derived from the
5-dropsonde circle agrees with that from the 10-dropsonde circle above the
inversion layer but differs within the boundary layer. The difference is even
larger between the two different sets of 5-dropsonde circles (red and blue
curves in Figure A4(a)) in the boundary layer. Interestingly, $\bar{D}$ from
the 4-dropsonde subset (Figure A4(b)) is closer to the one derived from the
10-dropsonde. Therefore, reducing the total number of dropsondes in a circle
results in statistical uncertainties. Bony and Stevens (2019) suggested that
at least 12 dropsondes are needed to estimate $\bar{D}$. More dropsondes can
indeed improve the accuracy of the estimation in the tropics as shown in
Figure 5 of Bony and Stevens (2019). However, the additional two dropsondes
are not expected to make a big difference.

A1Latitude and longitude coordinates of dropsondes released on February 28
case. Black dots represent the 10 dropsondes. (a): Red circles and blue
triangles represent two subsets of 5 dropsondes, respectively. (b): Red
circles, blue triangles, cyan triangles, and red crosses represent four
subsets of 4 dropsondes, respectively.

A2Corresponding $\bar{D}$ profiles derived from dropsonde measurements shown
in Figure A.

A3Same as Figure A1 but for the March 1 case.

A4Corresponding $\bar{D}$ profiles derived from groups of dropsondes shown in
Figure A.

##  Appendix B Surface heat fluxes: ERA5 versus MERRA-2

Figure A5 shows the comparison of heat fluxes between ERA5 and MERRA-2
reanalysis data. MERRA-2 underestimates the heat fluxes compared to ERA5.

A5Comparison of heat fluxes between ERA5 and MERRA-2 reanalysis data at 15:00
UTC for the February 28 case (black squares) anad 16:00 UTC for the March 1
case (red squares). Heat fluxes from the MERRA-2 reanalysis data are averaged
between 16:30 and 17:30 for the February 28 case and between 14:30 and 15:30
for the March 1 case.

##  Appendix C Vertical profiles for LES with different forcing

Figure A6 shows the evolution of vertical profiles of the meteorological state
for simulations shown in Figure 8.

A6Evolution of domain-averaged vertical profiles from the WRF-LES simulation
with the corresponding input forcings shown in Figure 7 for the February 28
case. The rainbow color scheme represents the time evolution (06:00-21:00
UTC): from purple to red. The solid black line represents the ERA5 reanalysis
data and the dashed one represent the WRF-LES simulation averaged during the
measurement time. The grey curve represents the dropsonde measurement. From
the top to the bottom, rows represent the simulation with no forcing, only
advective tendencies ($\partial\bar{\theta}/\partial t$ and
$\partial\bar{q}_{v}/\partial t$), only relaxation of $u$ and $v$, and
advective tendencies ($\partial\bar{\theta}/\partial t$
$\partial\bar{q}_{v}/\partial t$) plus $\bar{u}$ and $\bar{v}$ relaxation to
$u_{\rm ERA5}$ and $v_{\rm ERA5}$, respectively.

## References

  * Ackerman et al. (2009) Ackerman, A. S., and Coauthors, 2009: Large-eddy simulations of a drizzling, stratocumulus-topped marine boundary layer. Monthly Weather Review, 137 (3), 1083–1110.
  * Agee (1987) Agee, E. M., 1987: Mesoscale cellular convection over the oceans. Dynamics of atmospheres and oceans, 10 (4), 317–341.
  * Alexandrov et al. (2012) Alexandrov, M. D., B. Cairns, C. Emde, A. S. Ackerman, and B. van Diedenhoven, 2012: Accuracy assessments of cloud droplet size retrievals from polarized reflectance measurements by the research scanning polarimeter. Remote sensing of environment, 125, 92–111.
  * Alexandrov et al. (2018) Alexandrov, M. D., and Coauthors, 2018: Retrievals of cloud droplet size from the research scanning polarimeter data: Validation using in situ measurements. Remote Sensing of Environment, 210, 76–95.
  * Augstein et al. (1973) Augstein, E., H. Riehl, F. Ostapoff, and V. Wagner, 1973: Mass and energy transports in an undisturbed atlantic trade-wind flow. Monthly Weather Review, 101 (2), 101–111.
  * Bony and Stevens (2019) Bony, S., and B. Stevens, 2019: Measuring area-averaged vertical motions with dropsondes. Journal of the Atmospheric Sciences, 76 (3), 767–783.
  * Bony et al. (2017) Bony, S., B. Stevens, and D. Carlson, 2017: WMO Bulletin, 66, 8–11.
  * Boucher et al. (2013) Boucher, O., and Coauthors, 2013: Clouds and aerosols. Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, 571–657.
  * Bretherton et al. (1999) Bretherton, C. S., S. K. Krueger, M. C. Wyant, P. Bechtold, E. Van Meijgaard, B. Stevens, and J. Teixeira, 1999: A gcss boundary-layer cloud model intercomparison study of the first astex lagrangian experiment. Boundary-Layer Meteorology, 93 (3), 341–380.
  * Brown et al. (2002) Brown, A., and Coauthors, 2002: Large-eddy simulation of the diurnal cycle of shallow cumulus convection over land. Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, applied meteorology and physical oceanography, 128 (582), 1075–1093.
  * Chellappan et al. (2021) Chellappan, S., and Coauthors, 2021: On assessing era5 and merra2 representations of cold-air outbreaks across the gulf stream. Geophysical Research Letters.
  * Corral et al. (2021) Corral, A. F., and Coauthors, 2021: An overview of atmospheric features over the western north atlantic ocean and north american east coast–part 1: Analysis of aerosols, gases, and wet deposition chemistry. Journal of Geophysical Research: Atmospheres, e2020JD032592.
  * de Roode et al. (2019) de Roode, S. R., and Coauthors, 2019: Turbulent transport in the gray zone: A large eddy model intercomparison study of the constrain cold air outbreak case. Journal of Advances in Modeling Earth Systems, 11 (3), 597–623.
  * Endo et al. (2015) Endo, S., and Coauthors, 2015: RACORO continental boundary layer cloud investigations: 2. large-eddy simulations of cumulus clouds and evaluation with in situ and ground-based observations. Journal of Geophysical Research: Atmospheres, 120 (12), 5993–6014.
  * Gelaro et al. (2017) Gelaro, R., and Coauthors, 2017: The modern-era retrospective analysis for research and applications, version 2 (merra-2). Journal of climate, 30 (14), 5419–5454.
  * Global Modeling and Assimilation Office (2015) (GMAO) Global Modeling and Assimilation Office (GMAO), 2015: Merra-2 inst3$\\_$3d$\\_$asm$\\_$nv: 3d, 3-hourly, instantaneous, model-level, assimilation, assimilated meteorological fields v5.12.4. Greenbelt, MD, USA, Goddard Earth Sciences Data and Information Services Center (GES DISC), URL https://disc.gsfc.nasa.gov/datasets/M2I3NVASM˙5.12.4/summary, accessed March 10, 2021, 10.5067/WWQSXQ8IVFW8.
  * Gryschka and Raasch (2005) Gryschka, M., and S. Raasch, 2005: Roll convection during a cold air outbreak: A large eddy simulation with stationary model domain. Geophysical research letters, 32 (14).
  * Heinze et al. (2017) Heinze, R., C. Moseley, L. N. Böske, S. K. Muppa, V. Maurer, S. Raasch, and B. Stevens, 2017: Evaluation of large-eddy simulations forced with mesoscale model output for a multi-week period during a measurement campaign. Atmospheric Chemistry and Physics, 17 (11), 7083–7109.
  * Hersbach et al. (2020) Hersbach, H., and Coauthors, 2020: The era5 global reanalysis. Quarterly Journal of the Royal Meteorological Society, 146 (730), 1999–2049, https://doi.org/10.1002/qj.3803, URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3803.
  * Hirahara et al. (2016) Hirahara, S., M. A. Balmaseda, E. de Boisseson, and H. Hersbach, 2016: 26 sea surface temperature and sea ice concentration for era5.
  * Holland and Rasmusson (1973) Holland, J. Z., and E. M. Rasmusson, 1973: Measurements of the atmospheric mass, energy, and momentum budgets over a 500-kilometer square of tropical ocean. Monthly Weather Review, 101 (1), 44–55.
  * Knop et al. (2021) Knop, I., S. E. Bansmer, V. Hahn, and C. Voigt, 2021: Comparison of different droplet measurement techniques in the braunschweig icing wind tunnel. Atmospheric Measurement Techniques, 14 (2), 1761–1781.
  * Lenschow et al. (1999) Lenschow, D. H., P. B. Krummel, and S. T. Siems, 1999: Measuring entrainment, divergence, and vorticity on the mesoscale from aircraft. Journal of atmospheric and oceanic technology, 16 (10), 1384–1400.
  * Lenschow et al. (2007) Lenschow, D. H., V. Savic-Jovcic, and B. Stevens, 2007: Divergence and vorticity from aircraft air motion measurements. Journal of atmospheric and oceanic technology, 24 (12), 2062–2072.
  * Liu et al. (2004) Liu, A., G. Moore, K. Tsuboki, and I. Renfrew, 2004: A high-resolution simulation of convective roll clouds during a cold-air outbreak. Geophysical Research Letters, 31 (3).
  * Morrison et al. (2009) Morrison, H., G. Thompson, and V. Tatarskii, 2009: Impact of cloud microphysics on the development of trailing stratiform precipitation in a simulated squall line: Comparison of one-and two-moment schemes. Monthly weather review, 137 (3), 991–1007.
  * National Center for Atmospheric Research (2021) National Center for Atmospheric Research, 2021: Avaps dropsondes. URL https://www.eol.ucar.edu/content/avaps-dropsondes, accessed March 16, 2021.
  * National Centers for Environmental Information (2008) National Centers for Environmental Information, 2008: Remote sensing systems (remss). ghrsst level 4 mw$\\_$ir$\\_$oi global foundation sea surface temperature and analysis version 5.0 from remss (gds versions 1 and 2). URL https://www.ncei.noaa.gov/archive/accession/GHRSST-MW˙IR˙OI-REMSS-L4-GLOB, accessed March 16, 2021.
  * Neggers et al. (2012) Neggers, R. A., A. Siebesma, and T. Heus, 2012: Continuous single-column model evaluation at a permanent meteorological supersite. Bulletin of the American Meteorological Society, 93 (9), 1389–1400.
  * Painemal et al. (2021) Painemal, D., and Coauthors, 2021: An overview of atmospheric features over the western north atlantic ocean and north american east coast–part 2: Circulation, boundary layer, and clouds. Journal of Geophysical Research: Atmospheres, e2020JD033423.
  * Rahn and Garreaud (2010) Rahn, D. A., and R. Garreaud, 2010: Marine boundary layer over the subtropical southeast pacific during vocals-rex–part 1: Mean structure and diurnal cycle. Atmospheric Chemistry and Physics, 10 (10), 4491–4506.
  * Ramanathan et al. (1989) Ramanathan, V., R. D. Cess, E. F. Harrison, P. Minnis, B. R. Barkstrom, E. Ahmad, and D. Hartmann, 1989: Science, 243 (4887), 57–63.
  * Randall and Cripe (1999) Randall, D. A., and D. G. Cripe, 1999: Alternative methods for specification of observed forcing in single-column models and cloud system models. Journal of Geophysical Research: Atmospheres, 104 (D20), 24 527–24 545.
  * Schneider et al. (2017) Schneider, T., J. Teixeira, C. S. Bretherton, F. Brient, K. G. Pressel, C. Schär, and A. Siebesma, 2017: Nature Climate Change, 7 (1), 3–5.
  * Shaw (2003) Shaw, R. A., 2003: Particle-turbulence interactions in atmospheric clouds. Annu. Rev. Fluid Mech., 35 (1), 183–227.
  * Skamarock et al. (2019) Skamarock, W. C., and Coauthors, 2019: A description of the advanced research wrf model version 4. National Center for Atmospheric Research: Boulder, CO, USA, 145.
  * Smith (1988) Smith, S. D., 1988: Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature. Journal of Geophysical Research: Oceans, 93 (C12), 15 467–15 472.
  * Sorooshian et al. (2019) Sorooshian, A., and Coauthors, 2019: Aerosol–cloud–meteorology interaction airborne field investigations: Using lessons learned from the us west coast in the design of activate off the us east coast. Bulletin of the American Meteorological Society, 100 (8), 1511–1528.
  * Sorooshian et al. (2020) Sorooshian, A., and Coauthors, 2020: Atmospheric research over the western north atlantic ocean region and north american east coast: a review of past work and challenges ahead. Journal of Geophysical Research: Atmospheres, 125 (6), e2019JD031 626.
  * Stevens and Bony (2013) Stevens, B., and S. Bony, 2013: Science, 340 (6136), 1053–1054.
  * Stevens et al. (2002) Stevens, D. E., A. S. Ackerman, and C. S. Bretherton, 2002: Effects of domain size and numerical resolution on the simulation of shallow cumulus convection. Journal of the atmospheric sciences, 59 (23), 3285–3301.
  * Taylor et al. (2019) Taylor, J. W., and Coauthors, 2019: Aerosol influences on low-level clouds in the west african monsoon. Atmospheric Chemistry and Physics, 19 (13), 8503–8522.
  * Tomassini et al. (2017) Tomassini, L., P. R. Field, R. Honnert, S. Malardel, R. McTaggart-Cowan, K. Saitou, A. T. Noda, and A. Seifert, 2017: The “grey zone” cold air outbreak global model intercomparison: A cross evaluation using large-eddy simulations. Journal of Advances in Modeling Earth Systems, 9 (1), 39–64.
  * Tornow et al. (2021) Tornow, F., A. S. Ackerman, and A. M. Fridlind, 2021: Preconditioning of overcast-to-broken cloud transitions by riming in marine cold air outbreaks. Atmospheric Chemistry and Physics Discussions, 1–25.
  * Wang and Feingold (2009) Wang, H., and G. Feingold, 2009: Modeling mesoscale cellular structures and drizzle in marine stratocumulus. part i: Impact of drizzle on the formation and evolution of open cells. Journal of the Atmospheric Sciences, 66 (11), 3237–3256.
  * Wang et al. (2010) Wang, H., G. Feingold, R. Wood, and J. Kazil, 2010: Modelling microphysical and meteorological controls on precipitation and cloud cellular structures in southeast pacific stratocumulus. Atmospheric Chemistry and Physics, 10 (13), 6347–6362.
  * Wang and McFarquhar (2008) Wang, H., and G. M. McFarquhar, 2008: Modeling aerosol effects on shallow cumulus convection under various meteorological conditions observed over the indian ocean and implications for development of mass-flux parameterizations for climate models. Journal of Geophysical Research: Atmospheres, 113 (D20).
  * Wang et al. (2009) Wang, H., W. C. Skamarock, and G. Feingold, 2009: Evaluation of scalar advection schemes in the advanced research wrf model using large-eddy simulations of aerosol–cloud interactions. Monthly Weather Review, 137 (8), 2547–2558.
  * Zeng et al. (1998) Zeng, X., M. Zhao, and R. E. Dickinson, 1998: Intercomparison of bulk aerodynamic algorithms for the computation of sea surface fluxes using toga coare and tao data. Journal of Climate, 11 (10), 2628–2644.


$\displaystyle N_{2}(\psi)=\sum_{a,b,c=1,2}N^{abc}_{2}(\psi)\,,$ (D.5)

where the $N_{2}^{abc}(\psi)$, for $a,b,c=1,2$, denote the operators

$\displaystyle N^{abc}_{2}(\psi)=$ $\displaystyle\oint R_{a}\phi R_{b}\psi
R_{c}\,.$ (D.6)

We estimate the terms individually. Below, we use the estimate (see (4))

$\displaystyle\|(1-\Delta)^{\alpha}R_{i}(z)\|\lesssim\ $ $\displaystyle
d^{\alpha-1}\lesssim 1,\ i=1,2,$ (D.7)

for $\alpha\in[0,1]$ and $z\in\Gamma$, where

$\displaystyle d\equiv d(z):=\mathrm{dist}(z,\sigma(h_{\rm
per}))\geq\frac{1}{4}.$ (D.8)

Case 1: $(121)$ and $(212)$. We estimate the case for $(121)$, the other case
is done similarly. Since $P_{1}P_{2}=0$, we write

$\displaystyle N^{(121)}_{2}(\psi)=$ $\displaystyle\oint R_{1}\psi
R_{2}P_{2}\psi R_{1}$ (D.9) $\displaystyle=$ $\displaystyle\oint
R_{1}[P_{1},\psi]P_{2}R_{2}P_{2}[\psi,P_{1}]R_{1}\,.$ (D.10)

Applying Lemma 2.1 and equation (D.7) to the r.h.s. and using that the
operator norm is bounded by the $I^{2}$ norm, we find

$\displaystyle\|\operatorname{den}[N^{(121)}_{2}(\psi)]\|_{L^{2}}$
$\displaystyle\lesssim\|(1-\Delta)^{3/4+\epsilon}N^{(121)}_{2}(\phi)\|_{S^{2}}$
(D.11)
$\displaystyle\lesssim\left|\oint\right|d^{-3}\|[P_{1},\psi]P_{2}\|_{S^{2}}^{2}.$
(D.12)

where, recall,

$\displaystyle\left|\oint\right|:=$
$\displaystyle\frac{1}{2\pi}\int_{\Gamma}dz|f_{T}(z-\mu)|.$ (D.13)

A key observation allowing us to obtain an improved estimate is that the
commutators lead to gradient estimates:

######  Lemma D.4.

Let Assumption [A1] hold, we have the estimate

$\displaystyle\|[P_{i},\psi]\|_{S^{2}}\lesssim\|\nabla\psi\|_{L^{2}},\ i=1,2.$
(D.14)

###### Proof of Lemma D.4.

Since the identity commutes with any operator and $P_{2}=\mathbf{1}-P_{1}$
(see (D.4)), we prove the lemma for $P_{1}$ only. Since $h_{\rm per}$ (see
(3.8)) has a gap at $\mu$, the Cauchy integral formula implies

$\displaystyle P_{1}=\frac{1}{2\pi i}\int_{\Gamma_{1}}(z-h_{\rm
per})^{-1}=\frac{1}{2\pi i}\int_{\Gamma_{1}}r_{\rm per}(z)$ (D.15)

where $\Gamma_{1}$ is the contour $\\{t+i;-c\leq t<\mu\\}\cup\\{t-i;-c\leq
t<\mu\\}\cup\\{-c-it+(1-t)i:t\in[0,1]\\}\cup\\{\mu-it+(1-t)i:t\in[0,1]\\}$,
where $c>0$ is any constant such that $h_{\rm per}>-c+1$, and the contour is
traversed counter-clockwise. We see that

$\displaystyle[P_{1},\psi]=\frac{1}{2\pi i}\int_{\Gamma_{1}}[r_{\rm
per}(z),\psi]$ (D.16) $\displaystyle=$ $\displaystyle\frac{1}{2\pi
i}\int_{\Gamma_{1}}r_{\rm per}(z)[\nabla\cdot,\nabla\psi]r_{\rm per}(z)$
$\displaystyle+\frac{1}{2\pi i}\int_{\Gamma_{1}}r_{\rm
per}(z)(2\nabla\psi\cdot\nabla)r_{\rm per}(z).\,$ (D.17)

Lemma D.4 is now proved by an application of the Kato-Seiler-Simon inequality
((2.14)) to (D.17) and noting that $\Gamma_{1}$ is compact and has length
$O(1)$. ∎

Using Lemma D.4 and estimates (3.6) and (D.8) in (D.12) yields that

$\displaystyle\|\operatorname{den}[N^{(121)}_{2}(\psi)]\|_{L^{2}}\lesssim\|\nabla\psi\|_{L^{2}}^{2}\,.$
(D.18)

Case 2: $(112),(211),(122),$ $(221)$. We estimate the case for $(112)$, the
other cases are done similarly. Again, since $P_{1}P_{2}=0$, we write

$\displaystyle N^{(112)}_{2}(\psi)=$ $\displaystyle\oint R_{1}\psi R_{1}\psi
R_{2}$ (D.19) $\displaystyle=$ $\displaystyle\oint R_{1}\psi
R_{1}[\psi,P_{1}]R_{2}\,.$ (D.20)

Using Lemma 2.1 as in with $N^{(121)}_{2}(\psi)$ in (D.11), we estimate (D.20)
as

$\displaystyle\|\operatorname{den}[N^{(112)}_{2}(\psi)]\|_{L^{2}}$
$\displaystyle\lesssim\left|\oint\right|d^{-1}\|\psi
R_{1}\|\|[\psi,P_{1}]R_{2}\|_{S^{2}}.$ (D.21)

where $\left|\oint\right|$ is defined in (D.13). By the inequality
$\|A\|\leq\|A\|_{I^{p}}$, for any $p<\infty$ for any operator $A$ on
$L^{2}(\mathbb{R}^{3})$, and the Kato-Seiler-Simon inequality (2.14), we find
$\|\psi R_{1}\|\leq\|\psi R_{1}\|_{S^{6}}\lesssim\|\psi\|_{L^{6}}$. Using
this, together with Lemma D.4, in (D.21), we obtain

$\displaystyle\|\operatorname{den}[N^{(112)}_{2}(\psi)]\|_{L^{2}}\lesssim$
$\displaystyle\left|\oint\right|d^{-2}\|\psi\|_{L^{6}}\|\nabla\psi\|_{L^{2}}.$
(D.22)

Combining this with (3.6), (D.8) and Hardy-Littlewood’s inequality (2.18)
gives

$\displaystyle\|\operatorname{den}[N^{(112)}_{2}(\psi)]\|_{L^{2}}\lesssim$
$\displaystyle\|\nabla\psi\|_{L^{2}}^{2}.$ (D.23)

Case 3: $(111)$ and $(222)$. We use the $L^{2}$-$L^{2}$ duality to estimate
the $L^{2}$ norm of $\operatorname{den}[N^{(qqq)}_{2}(\psi)],q=1,2$. We have,
by (2.5) and definition (5.12),

$\displaystyle\|\operatorname{den}[N^{(qqq)}_{2}(\psi)]\|_{L^{2}}$
$\displaystyle=\sup_{\|f\|_{L^{2}}=1}|\int
f\operatorname{den}[N^{(qqq)}_{2}(\psi)]|$
$\displaystyle=\sup_{\|f\|_{L^{2}}=1}|\mathrm{Tr}[fN^{(qqq)}_{2}(\psi)]|$
$\displaystyle=\sup_{\|f\|_{L^{2}}=1}|\oint\mathrm{Tr}(fR_{q}\psi R_{q}\psi
R_{q})|.$ (D.24)

(In the last two lines, $f$ is considered as a multiplication operator.) To
show that the integral on the r.h.s. converges absolutely, we follow the
arguments in (5.17)-(5.20) to prove, for $q=0,1$,

$\displaystyle\big{|}\mathrm{Tr}(fR_{q}(\psi R_{q})^{2})\big{|}\lesssim$
$\displaystyle
d^{-4/3}\|f\|_{L^{2}}\|\nabla\psi\|_{L^{2}}^{4/3}\|\psi\|_{L^{2}}^{2/3}.$
(D.25)

Due to definition (D.8) of $d\equiv d(z)$, this shows that the integral on the
r.h.s. of (D.24) converges absolutely.

######  Lemma D.5.

For $q=1,2$, we have

$\displaystyle\oint$ $\displaystyle\mathrm{Tr}[fR_{q}gR_{q}hR_{q}]$
$\displaystyle=\frac{1}{2\pi
i}\int_{\Gamma_{q}}(f_{T}(z)-1)\mathrm{Tr}[fR_{q}gR_{q}hR_{q}].$ (D.26)

###### Proof.

Note that the contour $\Gamma$ in Figure 2 is the union of two disjoint
contours, $\Gamma=\Gamma_{1}\cup\Gamma_{2}$, with $\Gamma_{1}$ being the
closed contour and $\Gamma_{2}$ unbounded one (i.e. the parts of $\Gamma$ with
$\mathrm{Re}\,z<\mu$ and $\mathrm{Re}\,z>\mu$). Hence, since $P_{1}/P_{2}$
projects to the spectrum of $h_{\rm per}$ on the left/right of $\mu$, the
Cauchy theorem implies that

$\displaystyle\int_{\Gamma}dz\,$
$\displaystyle\mathrm{Tr}[fR_{q}gR_{q}hR_{q}]=\int_{\Gamma_{q}}dz\,\mathrm{Tr}[fR_{q}gR_{q}hR_{q}].$
(D.27)

Next, we note that, by Bloch’s theory, we have

$\displaystyle\int_{\Gamma_{q}}dz\,\mathrm{Tr}[fR_{q}gR_{q}hR_{q}]$
$\displaystyle=\int_{\Gamma_{q}}dz\,\int_{(\Omega^{*})^{3}}d\hat{k}d\hat{k}_{1}d\hat{k}_{2}$
(D.28) $\displaystyle\times\mathrm{Tr}_{L^{2}_{\rm
per}}f_{k-k_{1}}(R_{q})_{k_{1}}g_{k_{1}-k_{2}}(R_{q})_{k_{2}}h_{k_{2}-k}(R_{q})_{k}.$
(D.29)

Computing the trace in the complete orthonormal basis of eigenvectors
$\varphi_{m,k}$ of $(R_{q})_{k}$ (with eigevalues $\lambda_{m,k}$) and
inserting the complete orthonormal bases of eigenvectors $\varphi_{n,k_{1}}$
and $\varphi_{r,k_{2}}$ of $(R_{q})_{k_{1}}$ and $(R_{q})_{k_{2}}$ (with
eigevalues $\lambda_{n,k_{1}}$ and $\lambda_{r,k_{2}}$) into (D.29), we see
that

$\displaystyle\int_{\Gamma_{q}}\mathrm{Tr}[fR_{q}g$ $\displaystyle
R_{q}hR_{q}]=\sum_{m,n,r}\int_{(\Omega^{*})^{3}}d\hat{k}d\hat{k}_{1}d\hat{k}_{2}$
(D.30)
$\displaystyle\times\langle\varphi_{m,k},f_{k-k_{1}}\varphi_{n,k_{1}}\rangle\langle\varphi_{n,k_{1}},g_{k_{1}-k_{2}}\varphi_{r,k_{2}}\rangle\langle\varphi_{r,k_{2}},h_{k_{2}-k}\varphi_{m,k}\rangle$
(D.31)
$\displaystyle\times\int_{\Gamma_{q}}dz\,\frac{1}{(z-\lambda_{m,k})(z-\lambda_{n,k_{1}})(z-\lambda_{r,k_{2}})}.$
(D.32)

For $q=1$, since $P_{1}$ projects to the spectrum of $h_{\rm per}$ on the left
of $\mu$, we see that $\lambda_{m,k},\lambda_{n,k_{1}},\lambda_{r,k_{2}}<\mu$.
In particular, these eigenvalues are in the left closed contour in Figure 2.
Consequently, Cauchy’s integral formula shows that the integal in (D.32) is
identically zero. Similar argument applies to $q=2$. This shows that

$\displaystyle\frac{1}{2\pi
i}\int_{\Gamma_{q}}\mathrm{Tr}[fR_{q}gR_{q}hR_{q}]=0.$ (D.33)

Thus (D.26) follows. ∎

Using the explicit form of the Fermi-Diract distribution $f_{T}$ in (1.2), we
see that

$\displaystyle|f_{T}(z-\mu)-1|=\frac{e^{\beta(\mathrm{Re}\,z-\mu)}}{|1+e^{\beta(z-\mu)}|}\,.$
(D.34)

By condition (1.26) and by the choice of the contour, $\Gamma$, in Figure 2,
we see that the if $z\in\Gamma$, then $\mathrm{Re}\,z$ is at least at the
distance $\geq 1$ from $\mu$. Hence, for $z$ in a contour $\Gamma$, (D.34)
implies that

$\displaystyle|f_{T}(z-\mu)-1|\lesssim e^{-\beta}.$ (D.35)

Applying estimates (D.35) and (D.25) to the r.h.s. of (D.26) and recalling the
definition (D.8) of $d\equiv d(z)\geq\frac{1}{4}$, we arrive at the inequality

$\displaystyle|\oint$ $\displaystyle\mathrm{Tr}[fR_{q}gR_{q}hR_{q}]|$
$\displaystyle\lesssim
e^{-\beta}\|f\|_{L^{2}}\|\nabla\psi\|_{L^{2}}^{4/3}\|\psi\|_{L^{2}}^{2/3}.$
(D.36)

This inequality, together with the relation (D.24), gives

$\displaystyle\|\operatorname{den}[N^{(qqq)}_{2}(\phi)]\|_{L^{2}}\lesssim$
$\displaystyle e^{-\beta}\|\nabla\phi\|_{L^{2}}^{4/3}\|\phi\|_{L^{2}}^{2/3}.$
(D.37)

Inequalities (D.18), (D.23) and (D.37) imply estimate (D.3) for $k=2$.

Now we estimate $N_{k}$ for $k>2$. By (3.6) and (3.10), it suffices to
estimate $\operatorname{den}[R(\phi R)^{k}]$ where $R=r_{\rm per}(z)$ is given
in (3.8). Using Lemma 2.1, we see that

$\displaystyle\|\operatorname{den}[R(\phi R)^{k}]\|_{L^{2}}\lesssim$
$\displaystyle\|(1-\Delta)^{4/3+\epsilon}R(\phi R)^{k}\|_{S^{2}}$ (D.38)
$\displaystyle\lesssim$ $\displaystyle\|(\phi R)^{k}\|_{S^{2}}.$ (D.39)

Using Hölder’s inequality with
$\frac{1}{2}=\frac{1}{6}+\frac{1}{6}+\frac{1}{6}+$ another $k$ terms of
$\frac{1}{\infty}$, (D.39) becomes

$\displaystyle\|\operatorname{den}[R(\phi R)^{k}]\|_{L^{2}}\leq$
$\displaystyle\|\phi R\|_{S^{6}}^{3}\|\phi R\|^{k-3}$ (D.40)
$\displaystyle\leq$ $\displaystyle\|\phi R\|_{S^{6}}^{k},$ (D.41)

where the last line follows since $\|\cdot\|\leq\|\cdot\|_{S^{p}}$ for
$p<\infty$. Combining with Kato-Seiler-Simon’s inequality (2.14) and Hardy-
Littlewood’s inequality (2.18), (D.41) implies (D.3) for $k\geq 3$. ∎

The rest of the proof of Proposition D.1 proceeds as the proof of in
Proposition 5.3. ∎

## References

  * [1] A. Anantharaman, E. Cancès, Existence of minimizers for Kohn–Sham models in quantum chemistry, Ann. I. H. Poincaré – AN 26, 2425–2455 (2009).
  * [2] V. Bach, S. Breteaux, Th. Chen, J. M. Fröhlich and I.M. Sigal: The time-dependent Hartree-Fock-Bogoliubov equations for bosons, J. Evol Eqn, 2020 (to appear), arXiv:1602.05171v2.
  * [3] V. Bach, E. H. Lieb, and J. P. Solovej, Generalized Hartree-Fock theory and the Hubbard model, J. Stat. Phys., 76, 3–89 (1994).
  * [4] C. Brislawn, Kernels of trace class operators, Proceedings AMS 104, No. 4 (1988).
  * [5] E. Cancès, A. Deleurence, M. Lewin, A new approach to the modeling of local defects in crystals: the reduced Hartree-Fock case, Comm. Math. Phys., Vol. 281, Issue 1, 129–177 (2008).
  * [6] E. Cancès, A. Deleurence, M. Lewin, Non-perturbative embedding of local defects in crystalline materials. J. Phys.: Condens. Matter, 20:294213 (2008).
  * [7] E. Cancès, M. Lewin, The dielectric permittivity of crystals in the reduced Hartree-Fock approximation, Arch. Ration. Mech. Anal. 197, Issue 1, 139–177 (2010).
  * [8] E. Cancès, M. Lewin and G. Stoltz, The microscopic origin of the macroscopic dielectric permittivity of crystals, in: Lecture Notes in Computational Science and Engineering 82, Springer 2011.
  * [9] E. Cancès, G. Stoltz, A mathematical formulation of the random phase approximation for crystals, Ann. I. H. Poincaré – AN 29 , Issue 6 , 887 - 925 (2012).
  * [10] I. Catto, C. Le Bris, P.-L. Lions, On the thermodynamic limit for Hartree-Fock type models. Ann. I. H. Poincaré – AN 18, Issue 6, 687–760 (2001).
  * [11] I. Catto, C. Le Bris, P. -L. Lions, On some periodic Hartree type models, Ann. I. H. Poincaré – AN 19, Issue 2, 143 - 190 (2002).
  * [12] I. Chenn and I. M. Sigal, On the Bogolubov-de Gennes equations, [arXiv:1701.06080v2] (2019).
  * [13] I. Chenn and I. M. Sigal, On effective PDEs of quantum physics, in New Tools for Nonlinear PDEs and Application, M. D’Abbicco et al. (eds.), Birkhäuser series ”Trends in Mathematics”, https://doi.org/.
  * [14] H. Cycon, R. Froese, W. Kirsch, B. Simon, Schrödinger Operators (with Applications to Quantum Mechanics and Global Geometry). Springer: Berlin (1987).
  * [15] W. E, J. Lu, Electronic structure of smoothly deformed crystals: Cauchy-Born Rule for the nonlinear tight-binding model, Commun. Pure Applied Math. 63(11), 1432–1468 (2010).
  * [16] W. E, J. Lu, The electronic structure of smoothly deformed crystals: Wannier functions and the Cauchy-Born rule, Arch. Ration. Mech. Anal. 199 no. 2, pp. 407–433 (2011).
  * [17] W. E, J. Lu, The Kohn-Sham equation for deformed crystals, Memoirs AMS (2013).
  * [18] W. E, J. Lu, X. Yang, Effective Maxwell equations from time-dependent density functional theory. Acta. Math. Sin.-English Ser. 27: 339–368 (2011).
  * [19] F. Fogolari, A. Brigo and H. Molinari, The Poisson-Boltzmann equation for biomolecular electrostatics: a tool for structural biology, J. Mol. Recognit. 15: 377–392 (2002).
  * [20] S. J. Gustafson, I. M. Sigal, Mathematical Concepts of Quantum Mechanics Second Edition. Universitext, Springer-Verlag 2011.
  * [21] Ch. Hainzl, M. Lewin, E. Seré, Existence of atoms and molecules in the mean-field approximation of no-photon quantum electrodynamics. Arch. Ration. Mech. Anal 192, 3, 453-499 (2009).
  * [22] C. Le Bris and P.-L. Lions, From atoms to crystals: a mathematical journey, Journal: Bull. AMS 42, 291–363 (2005).
  * [23] A. Levitt, Screening in the finite-temperature reduced Hartree-Fock model, arXiv:1810.03342v1 (2018).
  * [24] E. Lieb, M. Loss. Analysis, 2nd Edition. AMS Press, 2001.
  * [25] E. H. Lieb, The stability of matter: From atoms to stars, Bull. AMS 22, 1-49 (1990).
  * [26] E.H. Lieb, B. Simon. The Hartree-Fock theory for Coulomb systems, Commun. Math. Phys. 53, no.3, 185–194 (1977).
  * [27] G. Lindblad, Expectations and entropy inequalities for finite quantum systems. Commun. Math. Phys. 39, 111–119 (1974).
  * [28] P.L. Lions, Hartree-Fock and related equations. Nonlinear Partial Differential Equations and Their Applications. Collège de France Seminar, Vol. IX. Pitman Res. Notes Math. Ser. 181, 304–333 (1988).
  * [29] P.L. Lions, Solutions of Hartree-Fock equations for Coulomb systems, Comm. Math. Phys. 109, 33-97 (1987).
  * [30] P. A. Markowich, G. Rein and G. Wolansky, Existence and nonlinear stability of stationary states of the Schrödinger-Poisson system, J. Stat. Phys. 106, Issue 5–6, pp 1221–1239 (2002).
  * [31] N. D. Mermin, Thermal properties of the inhomogeneous electron gas, Phys. Rev. 137, A: 1441 (1965).
  * [32] F. Nier, A variational formulation of Schrödinger-Poisson systems in dimension $d\leq 3$. Commun. PDEs 18(7 and 8): 1125–1147 (1993).
  * [33] E. Prodan and P. Nordlander. On the Kohn-Sham equations with periodic background potential. J. Stat. Phys. 111, issue 3-4, 967-992 (2003).
  * [34] M. Reed, B. Simon, Method of Modern Mathematical Physics IV: Analysis of Operators. Academic Press, Inc. (Londond) Ltd. (1978).
  * [35] C. Sevik, C. Bulutay. Theoretical study of the insulating oxides and nitrate: $SiO_{2}$, $GeO_{2}$, $Al_{2}O_{3}$, $Si_{3}N_{4}$, and $Ge_{3}N_{4}$, arXiv:cond-mat/0610176v2 (2008).
  * [36] B. Simon. Trace Ideals and Their Applications. 2nd Edition. AMS Press, 2005.