Realism of Lagrangian Large Eddy Simulations Driven by Reanalysis Meteorology: Tracking a Pocket of Open Cells Under a Biomass Burning Aerosol Layer.

An approach to drive Lagrangian large eddy simulation (LES) of boundary layer clouds with reanalysis data is presented and evaluated using satellite (Spinning Enhanced Visible and Infrared Imager, SEVIRI) and aircraft (Cloud-Aerosol-Radiation Interactions and Forcing, CLARIFY) measurements. The simulations follow trajectories of the boundary layer flow. They track the formation and evolution of a pocket of open cells (POC) underneath a biomass burning aerosol layer in the free troposphere. The simulations reproduce the evolution of observed stratocumulus cloud morphology, cloud optical depth, and cloud drop effective radius, and capture the timing of the cloud state transition from closed to open cells seen in the satellite imagery on the three considered trajectories. They reproduce a biomass burning aerosol layer identified by the in-situ aircraft measurements above the inversion of the POC. Entrainment of aerosol from the biomass burning layer into the POC is limited to the extent of having no impact on cloud- or boundary layer properties, in agreement with the CLARIFY observations. The two-moment bin microphysics scheme used in the simulations reproduces the in-situ cloud microphysical properties reasonably well. A two-moment bulk microphysics scheme reproduces the satellite observations in the non-precipitating closed-cell state, but overestimates liquid water path and cloud optical depth in the precipitating open-cell state due to insufficient surface precipitation. A boundary layer cold and dry bias occurring in LES can be counteracted by reducing the grid aspect ratio and by tightening the large scale wind speed nudging towards the surface.

Introduction

The climatically important cloud decks in the eastern subtropical oceans undergo an evolution during their passage to the equator. They begin their journey commonly as shallow marine stratus, grow into stratocumulus, and mature into trade cumuli. In this evolution, they can transition directly from an overcast state into trade cumuli, or through stages of less or more organized stratocumulus states, associated with different modes of boundary layer circulation (Wood, 2012).

The direct transition proceeds from a shallow, well‐mixed stratocumulus‐topped boundary layer to a deeper, decoupled boundary layer with cumulus rising into stratocumulus. This is followed by the dissipation of the overlying stratocumulus deck, which leaves behind a trade cumulus cloud field (Krueger et al., 1995a, 1995b). The underlying mechanism is the deepening of the boundary layer which is accompanied by a warming and decoupling, without precipitation (Bretherton & Wyant, 1997; Wyant et al., 1997). Free tropospheric aerosol can modulate the transition by radiative and microphysical mechanisms (Yamaguchi et al., 2015; Zhou et al., 2017). Evidence is growing that precipitation is also capable of driving the transition (Sarkar et al., 2020; Yamaguchi et al., 2017).

The staged transition can pass through the organized closed‐ and open cell stratocumulus states (Agee, 1984, 1987; Agee et al., 1973; Atkinson & Zhang, 1996). The closed‐cell state has a cloud fraction with a median of 0.9, while open cells exhibit a markedly lower cloud fraction with a median of about 0.5 (Wood & Hartmann, 2006), as well as a smaller cloud radiative effect (Goren & Rosenfeld, 2014). Observational (Bretherton et al., 2010; Comstock et al., 2005; Stevens et al., 2005; Wood et al., 2008, Wood, Bretherton, et al., 2011, Wood, Mechoso, et al., 2011) and modeling (Savic‐Jovcic & Stevens, 2008; Wang & Feingold, 2009a; Xue et al., 2008) studies show that precipitation is a necessary but not a sufficient condition for the transition from closed to open cells: precipitation needs to be sufficiently strong over a sufficiently large area, or have a spatial distribution that is conducive for the transition to occur (Yamaguchi & Feingold, 2015). Precipitation also maintains the open‐cell state and its spatial and temporal oscillations (Feingold et al., 2010).

Satellite imagery (Agee, 1987) indicates that the preferred cloud state evolution is from the closed‐to the open‐cell state in the cloud sheets of the subtropical eastern oceans. The reverse transition, from the open‐to the closed‐cell state, has been proposed by Rosenfeld et al. (2006) and identified in satellite observations by Goren and Rosenfeld (2012), at locations where aerosol particles from ship exhaust entered the cloud deck. This reverse transition occurs less readily in subtropical stratocumulus decks because it requires restoration of liquid water and cloud top cooling of the closed‐cell state by sustained suppression of precipitation with a substantial aerosol source (Feingold et al., 2015).

The onset and progress of cloud state transitions is tied to the state of the atmosphere and ocean, such as sea surface temperature, subsidence, lower tropospheric stability, free tropospheric humidity, and boundary layer depth (Agee, 1987; Bretherton & Wyant, 1997; Chung & Teixeira, 2012; Eastman & Wood, 2016, 2018; Eastman et al., 2017; Mauger & Norris, 2010; Mechem et al., 2012; Pincus et al., 1997; Sandu & Stevens, 2011; Sandu et al., 2010; van der Dussen et al., 2016; Wang et al., 2010; Wyant et al., 1997). When precipitation drives the transition, higher aerosol levels delay the onset, as found in simulations (Mechem et al., 2012; Wang et al., 2010; Yamaguchi & Feingold, 2015) and satellite observations (Gryspeerdt et al., 2014).

Atmospheric, oceanic, or aerosol conditions may hence shift the boundary between overcast and broken clouds up‐ or downstream, and reduce or increase the size of areas with high cloud fraction in the subtropical cloud sheets of the eastern oceans. Turbulence‐resolving simulations and analysis of emergent constraints using observations show a robust positive cloud feedback to climate change with a contribution from a faster transition from stratocumulus to cumulus as climate warms (Nuijens & Siebesma, 2019). Goren et al. (2019) showed, using satellite data and Lagrangian large eddy simulations (LES) driven by reanalysis meteorology, that the timing of the closed‐to open‐cell transition varies systematically with aerosol concentration, with higher aerosol concentrations delaying the transition, even in polluted conditions. Christensen et al. (2020) analyzed satellite data along Lagrangian trajectories spanning several days along stratus‐to‐cumulus transition. They found that clouds forming on relatively polluted trajectories tend to have higher cloud albedo and cloud fraction compared with unpolluted trajectories. The response of cloud state transitions to environmental conditions therefore connects anthropogenic climate change, aerosol emissions, and Earth's radiation balance.

Low clouds represent a challenge to the fidelity of climate models (Bony & Dufresne, 2005; Lin et al., 2014; Vial et al., 2013; Williams & Webb, 2009), and cloud state transitions contribute to the challenge. Teixeira et al. (2011) evaluated an array of models along a Pacific Ocean cross section, from the stratocumulus regions off the coast of California, across the shallow convection‐dominated trade winds, to the deep convection regions of the intertropical convergence zone. They found that the stratocumulus‐to‐cumulus transition occurred too early along the trade wind Lagrangian trajectory. The transition also occurred either too abruptly or too smoothly, depending on model, with observations in‐between the extremes.

LES is the numerical tool of choice for the study of boundary layer clouds. In the Eulerian framework, in which the simulation domain is stationary relative to the surface, it performs well against surface‐based remote sensing and aircraft in‐situ observations (Berner et al., 2011; Kazil et al., 2011; Yamaguchi et al., 2013). Atlas et al. (2020) studied how well Eulerian LES driven by reanalysis meteorology, represent observed boundary layer structures and low clouds over the summertime Southern Ocean. LES also capture well the observed boundary layer and cloud state when following the course of a ship, provided horizontal winds relative to the simulation domain and thus large‐scale advective forcings are not too strong (McGibbon & Bretherton, 2017).

Eulerian LES contains advective terms which represent horizontal transport of heat and moisture. With these terms, Eulerian LES is able to represent the evolution of an advecting cloud field, provided the upstream state of the boundary layer at an earlier time is not too different from the state in the Eulerian simulation at that time. This needs not hold in reality: for example, the local and the upstream mesoscale flow can differ substantially, and the advective tendencies of heat and moisture may be unable to force the Eulerian simulation towards the advected state.

Lagrangian LES, in contrast, tracks the cloud state evolution in a specific air mass, a more realistic approach compared to Eulerian LES. Lagrangian LES is, however, subject to uncertainty arising from a number of open questions that do not affect Eulerian LES. One such question pertains to the optimal height above the surface at which the velocity of the Lagrangian simulation domain is calculated from the large scale wind field. Another question is the role and representation of horizontal advection of heat and moisture due to vertical shear in the boundary layer, which may alter boundary layer properties even at zero net horizontal advection in the boundary layer.

We note that strictly speaking, the LES in this and other works are quasi‐Lagrangian. The qualification ”quasi” is often omitted for convenience. Truly Lagrangian simulations have no advective interaction between the simulation domain and its exterior. This condition is not met, for example, when applying tendencies or nudging to represent horizontal advection in the free troposphere, which moves, in general, in a different direction than the boundary layer.

Lagrangian LES have been used extensively to study boundary layer cloud state transitions (de Roode et al., 2016; Krueger et al., 1995a, 1995b; Sandu & Stevens, 2011; Wyant et al., 1997; Yamaguchi et al., 2017; Yamaguchi & Feingold, 2015). These Lagrangian LES employed idealized initial and boundary conditions and forcings. Sandu and Stevens (2011) showed, using forcing composited from a set of trajectories in a reanalysis meteorology, that Lagrangian LES can represent well the major features of the observed stratocumulus‐to‐cumulus transition. van der Dussen et al. (2013) conducted an intercomparison of LES of a Lagrangian transition from a vertically well‐mixed stratocumulus‐topped boundary layer to a state in which shallow cumuli penetrate an overlying layer of thin and broken stratocumulus. They constructed the initial and forcing data from reanalysis meteorology and aircraft measurements, and showed that the model results agree remarkably well with aircraft observations. Neggers et al. (2019) studied Arctic cloudy mixed layers using Lagrangian LES driven with forcings and boundary conditions estimated from analysis and forecast products of the European Centre for Medium‐Range Weather Forecasts (ECMWF), and with calibrated initial conditions to reproduce ship‐based observations.

We have developed an approach to drive Lagrangian LES with reanalysis data to study the evolution of boundary layer clouds and their state transitions. It differentiates itself from the approaches of Sandu and Stevens (2011), van der Dussen et al. (2013), and Neggers et al. (2019) by the absence of special treatment of the initial and forcing data, such as compositing or the inclusion of observations. Furthermore, the LES has greater freedom to shape boundary layer and cloud properties in that hydrometeor numbers are not fixed, but calculated interactively (e.g., Yamaguchi et al., 2015). This enables the interaction of cloud microphysics with the aerosol population, required for microphysically driven cloud state transitions. The approach is hence more general and allows for a more stringent standard in evaluating Lagrangian LES of boundary layer clouds against observations compared to previous works. The fidelity of the reanalysis meteorology used to drive the simulations is a necessary but not a sufficient pre‐requisite for the LES to reproduce the observations. The approach was used by Goren et al. (2019) to study the evolution and response to anthropogenic aerosol of a mid‐latitude cloud deck in continental outflow, and by Narenpitak et al. (2021) to study the transition from the ”sugar” to the ”flower” trade cumulus cloud state.

The purpose of the current work is to document and evaluate the approach, based on two‐day simulations of a sub‐tropical cloud state transition, using satellite observations covering the simulation period, and aircraft profiles at its end. We show that the approach realistically simulates the observed clouds and their evolution, and determine key elements in the model formulation and simulation setup that are essential for its fidelity. We examine uncertainties and biases and identify research and development needs for Lagrangian LES driven by reanalysis meteorology.

The paper is organized as follows: Section 2 introduces the methods and data. Section 3 presents the simulation results and their evaluation, and explores the role of model and simulation setup. Section 4 discusses uncertainties and biases, and research and development needs. A summary and conclusions are given in Section 5.

Methods and Data

Observed Cloud State Evolution and Trajectories

We study a pocket of open cells (POC) sampled during flight C052 of the Cloud‐Aerosol‐Radiation Interactions and Forcing (CLARIFY) campaign (Abel et al., 2020; Haywood et al., 2021). The cloud state evolution is documented with imagery from the Spinning Enhanced Visible and Infrared Imager (SEVIRI) onboard the Meteosat Second Generation (MSG) satellite in Figure 1, and the Movie S1. Simulations in this work follow three distinct boundary layer air mass trajectories (green, blue, and red) that begin on September 3, 2017, 14:45:00 UTC and end on September 5, 2017, 17:00:00 UTC. We determined the trajectories from the wind field of the fifth generation of the ECMWF atmospheric reanalysis (ERA5, Hersbach et al., 2020) at a resolution of 0.3°, using the Hybrid Single Particle Lagrangian Integrated Trajectory Model (HYSPLIT, Stein et al., 2015). The trajectories are located at a constant height of 500 m above sea level. This height was chosen to distance the trajectory from shear effects near the surface and the inversion. This choice of trajectory height is supported by very small or absent lateral motion of the cloud field relative to the simulation domain in the cases considered in this work, indicating that the bulk of the boundary layer moves with the velocity at the chosen height. A different trajectory height may be appropriate for different cases.

Figure 1

Meteosat Second Generation Spinning Enhanced Visible and Infrared Imager imagery, with simulation domains on the green, blue, and red trajectory, to scale.

The transition from the closed‐to the open‐cell stratocumulus state occurs at different times on each trajectory. On September 4, 2017, 06:00:00 UTC (Figure 1a), the POC has begun to form on the red trajectory. The stratocumulus deck is still in the overcast, closed cell state at this time on the blue and green trajectories. By September 4, 2017, 14:30:00 UTC (Figure 1b), an open cell state has formed on the red trajectory. The transition from the closed‐to the open cell state is in progress on the blue trajectory at this time, while a closed‐cell state is still present on the green trajectory. The next day, on September 5, 2017, 10:45:00 UTC (Figure 1c), an open cell state is present on each trajectory, with differences in morphology: open cells are distinctly smaller on the red trajectory compared to the blue and green trajectories. On September 5, 2017, 14:00:00 UTC (Figure 1d), the POC is beginning to dissipate.

Model and Simulation Setup

We use the System for Atmospheric Modeling (SAM; Khairoutdinov & Randall, 2003), version 6.10.10, with periodic lateral boundary conditions.

Large Scale Meteorology

We use meteorology and sea surface temperature from ERA5 at 0.3° resolution to drive the simulations. ERA5 assimilates radiosonde profiles and satellite radiances (Hersbach et al., 2020), which helps to capture the effect of phenomena that are not represented by the underlying model, such as heating due to absorption of radiation by aerosol. Adebiyi et al. (2015) showed that the ERA‐Interim reanalysis, the predecessor of ERA5, captures thermodynamic profiles measured by radiosondes in the South‐East Atlantic better than other reanalysis products when compared to radiosonde measurements, under the caveat that the evaluated reanalyses, to different degrees, assimilate radiosonde data.

We use ERA5 temperature and moisture profiles to initialize the simulations, and nudge mean temperature and water vapor in the free troposphere towards ERA5 with Newtonian relaxation. The nudging tendencies are calculated from the difference of the horizontal mean in the simulation and the ERA5 value, hence variance of temperature and water vapor is preserved. Nudging begins 100 m above the inversion in the simulation or in ERA5, whichever is higher. From this nudging base level, the nudging tendencies increase smoothly over a height interval of 500 m from a value of zero to a value corresponding to the nudging time scale of 1,800 s. The inversion is diagnosed at the height of the maximum vertical gradient of liquid water static energy in the simulations, and at the height of the maximum vertical gradient of liquid water potential temperature in ERA5.

We nudge mean horizontal wind speed towards ERA5 at all levels with Newtonian relaxation. When model levels are located below the lowest level of ERA5, we extrapolate the ERA5 wind speed towards the surface assuming a logarithmic wind profile. Nudging tendencies are calculated from the difference of the horizontal mean and the ERA5 value, which preserves the variance of horizontal wind speed. The nudging time scale is 10 s between the surface and 500 m, and 1,800 s above 1,000 m, with a smooth interpolation in‐between. The short nudging time scale near the surface counters deceleration by surface drag and maintains the mean wind speed close to ERA5 values. The more relaxed nudging above 500 m allows the simulations to establish their own wind speed structure around an inversion height of their choice, rather than conforming to the wind speed structure at the inversion height of ERA5. The impact and relevance of the short time scale for mean wind speed nudging is discussed in Section 3.7.2. We apply ERA5 profiles of vertical velocity (subsidence) to temperature, water vapor, and aerosol.

Cloud Microphysics

We represent cloud microphysical processes with a bin or a bulk scheme. The bin microphysics scheme is the Tel Aviv University two‐moment bin microphysics model (Feingold et al., 1996; Tzivion et al., 1987). The hydrometeor size distribution is divided into 33 bins with mass doubling from one bin to the next. The smallest droplet radius is 1.56 μm. Cloud and rain hydrometeors are distinguished for diagnostic purposes by a threshold radius of 25 μm. Supersaturation is calculated based on the balance of dynamical and microphysical source and sink terms over the course of a time step (Clark, 1973). Activation of aerosol is based on the predicted supersaturation. Condensation and evaporation are computed via vapor diffusion to/from drops using the method of Stevens et al. (1996). Collection processes are based on Tzivion et al. (1987) and breakup processes on Feingold et al. (1988). The collection kernels are based on collision efficiencies after Hall (1980) as well as coalescence efficiencies for drizzle (Ochs et al., 1986) and raindrops (Low & List, 1982). In the rain drop regime where drops are unstable enough to be able to breakup as a result of binary collisions, the breakup efficiency is assumed to be 1 minus the coalescence efficiency. Drop sedimentation is computed with a first‐order upwind scheme. The bin microphysics scheme as implemented in SAM is described in further detail by Yamaguchi et al. (2019).

The bulk microphysics is a two‐moment bin‐emulating method (Feingold et al., 1998; Wang & Feingold, 2009a, 2009b) that calculates mass and number of hydrometeors. Cloud and rain water modes are represented using lognormal functions with fixed geometric standard deviation of 1.2. The threshold between the two modes is a radius of 25 μm. Supersaturation and aerosol activation are calculated as in the bin microphysics scheme. Condensation and evaporation are calculated analytically. Sedimentation of mass and number mixing ratios is calculated from mass‐ and number‐weighted average sedimentation velocities, respectively, and for each hydrometeor mode. Hydrometeor breakup is not implemented. The bulk microphysics as implemented in SAM is described in further detail by Yamaguchi et al. (2017).

In both microphysics schemes, advection is applied to the total mass mixing ratio (sum of vapor and condensate) and total number concentration (sum of aerosol and hydrometeors). Water vapor mixing ratio and aerosol number concentration are diagnostic variables. This implementation implicitly maintains the budget of both mass mixing ratio and number concentration through cloud microphysical processes. Further details of the implementation are given in Yamaguchi et al. (2019).

In both microphysics schemes, aerosol particles are activated in supersaturated conditions and removed from the aerosol population, increasing the hydrometeor number by the same amount. Collision‐coalescence reduces the hydrometeor number, thereby allowing for cloud processing of the aerosol. Upon evaporation, hydrometeors release one aerosol particle for each evaporated drop (Mitra et al., 1992). Surface precipitation removes hydrometeors and the corresponding number of aerosol particles from the atmosphere.

Aerosol

We use a simplified representation of the aerosol size distribution with a lognormal mode with a geometric‐mean diameter D  = 200 nm and a geometric standard deviation σ = 1.5. These parameters are consistent with the aerosol accumulation mode size distribution measured by the Passive Cavity Aerosol Spectrometer Probe (PCASP, Rosenberg et al., 2012) during CLARIFY flight C052 in both the overcast stratiform region surrounding the POC (D  = 186 nm, σ = 1.51), and within the free‐tropospheric biomass burning aerosol layer above the stratiform and POC cloud regimes (D  = 206 nm, σ = 1.53). The surface flux of ocean‐emitted aerosol is calculated with the parameterization of sea salt aerosol production of Clarke et al. (2006). The whitecap fraction is parameterized with the expression of Monahan et al. (1986) as a function of wind speed at 10 m above the ocean surface.

Radiation

Radiation is computed every 10 s from the distribution of temperature, gas phase constituents, and liquid water mass mixing ratio and cloud drop effective radius, with the Rapid Radiative Transfer Model (RRTMG, Iacono et al., 2008; Mlawer et al., 1997). Between the top of the model domain and the top of the atmosphere radiation is calculated with profiles of temperature, water vapor, and ozone from ERA5. The ocean surface albedo is set to 0.06, emissivity to 0.95. CO2 is set to the September 2017 value of 403 ppm (McGee, 2020).

In the case considered in this work, free tropospheric biomass burning aerosol remains above the POC inversion (Abel et al., 2020). It will be shown that this is also the case in the simulations (Sec. 3.4). Aerosol in the boundary layer is assumed to be predominantly sea spray, with negligible interaction with radiation. The interaction between aerosol and radiation is therefore not treated explicitly in the simulations, and the representation of effects from heating due to the absorption of radiation by biomass aerosol in the free troposphere is delegated to ERA5 and its assimilation of radiosonde and satellite data (Hersbach et al., 2020).

Numerics

SAM solves the anelastic system of equations using the finite difference approximation formulated on the Arakawa C grid with a height coordinate. Velocity components are predicted using the third‐order Adams‐Bashforth scheme (Durran, 1991) and the second order center advection scheme. Scalars are integrated with the forward in time, monotonic fifth‐order advection scheme of Yamaguchi et al. (2011). Diffusion is explicitly computed with eddy coefficients obtained from the 1.5‐order subgrid scale turbulence kinetic energy (TKE) scheme based on Deardorff (1980). Surface fluxes of sensible and latent heat and of momentum are computed based on Monin‐Obukhov similarity (Monin & Obukhov, 1954). The time step of 1 s is dynamically shortened by SAM to meet the Courant‐Friedrichs‐Levy condition.

Advection due to subsidence is solved in SAM with the advective form of the transport equation. This approach is known to preserve shape but not mass, resulting in spurious sink or source terms in the presence of velocity and tracer gradients. Strong tracer gradients exist at the upper and lower boundary of biomass burning layers. We maintain vertically integrated tracer mass in the layer by re‐normalizing its vertically integrated tracer mass in the free troposphere after the model applies advection due to subsidence. The correction is not applied outside the biomass burning layer or in the boundary layer, where mixing due to turbulence quickly dissipates strong gradients.

Domain Size, Grid, and Sampling

Sufficiently large horizontal domains sizes are required to capture mesoscale organization and the associated properties of the open‐cell (Feingold et al., 2010) and closed‐cell (Kazil et al., 2017) stratocumulus cloud state. We use domain sizes of 76.8 × 76.8 km2 and 38.4 × 38.4 km2 depending on the open cell size seen in the SEVIRI imagery on a given trajectory. The simulations employ a horizontal grid spacing of dx = dy = 200 m, and a vertical grid with variable spacing. At the surface, the thickness of the first three (mass) levels is dz1 = 35 m, dz2 = 22.5 m, and dz3 = 12.5 m. dz is 10 m to an altitude of 1,965 m, and 20 m to an altitude of 4,025 m, and coarsens thereabove by 10% per level to the domain top at 7,000 m. A grid with a constant dz = 10 m from the surface to 1,965 m and otherwise identical grid structure is also tested. 3D fields are saved every hour, 2D fields and domain mean profiles every minute. The results are sampled as a function of fractional day of year d, with d = 0 corresponding to January 1, 00h00m00s.

Simulations

The simulations are run from September 3, 2017, 14:45:00 UTC to September 5, 2017, 17:00:00 UTC (d = 245.61458 to 247.70833). We analyze results starting on September 4, 2017, 06:00:00 UTC (d = 246.25), allowing 15.25 hr for spin‐up. The simulations and their setup are listed in Table 1. The simulations G , B , and R run on the green, blue, and red trajectory, respectively; i is the simulation number. The simulations G and B use a 76.8 km domain to capture the larger open cell size on their trajectories (Figure 1). Smaller open cells are present along the red trajectory, and the simulations R use a 38.4 km domain.

Table 1

Simulations and Setup

Simulation	Trajectory	Initial aerosol (# mg−1)			Biomass Burning layer (m)	Cloud microphysics	Grid structure near surface			τuv0−500 (s)	Domain size (km2)
							dz1	dz2	dz3
		BL	FT	BB layer			(m)
G0	Green	145	37.5	–	–	Bin	35	22.5	12.5	10	76.8 × 76.8
G1	Green	145	37.5	4,700	3,100–3,700	Bin	35	22.5	12.5	10	76.8 × 76.8
G2	Green	145	37.5	4,700	3,100–3,700	Bulk	35	22.5	12.5	10	76.8 × 76.8
G3	Green	145	37.5	4,700	3,100–3,700	Bulk	10	10	10	10	76.8 × 76.8
G4	Green	145	37.5	4,700	3,100–3,700	Bulk	10	10	10	1,800	76.8 × 76.8
B1	Blue	145	37.5	4,700	3,100–3,700	Bin	35	22.5	12.5	10	76.8 × 76.8
B2	Blue	145	37.5	4,700	3,100–3,700	Bulk	35	22.5	12.5	10	76.8 × 76.8
R1	Red	115	37.5	4,700	3,100–3,700	Bin	35	22.5	12.5	10	38.4 × 38.4
R2	Red	115	37.5	4,700	3,100–3,700	Bulk	35	22.5	12.5	10	38.4 × 38.4
R3	Red	145	37.5	4,700	3,100–3,700	Bin	35	22.5	12.5	10	38.4 × 38.4

Note. dz is the thickness of the model level i, counted from the surface. is the nudging time constant for domain mean horizontal wind speed between the surface and 500 m. BL = boundary layer, FT = free troposphere.

Space‐borne lidar (Cloud‐Aerosol Lidar with Orthogonal Polarization, CALIOP) measurements on September 3, 2017, 01:35:00 UTC, 13 hr before the start of the simulations, show an aerosol layer approximately between 3–4.5 km, upstream of the POC locations sampled by CLARIFY flight C052 (Abel et al., 2020, Figure 4d). Based on this observation, we initialize the aerosol profiles with a biomass burning aerosol layer between 3,100–3,700 m (Figure S1 in Supporting Information S1). The mean aerosol number mixing ratio in the layer is set to 4,700 mg−1, corresponding to a mean number concentration of 3,850 cm−3. Outside the biomass burning aerosol layer, the initial aerosol number mixing ratio is 37.5 mg−1 in the free troposphere. In the boundary layer, we set the aerosol number mixing ratio to 145 mg−1 on the green and blue trajectories, and to 115 mg−1 on the red trajectory. The lower value is motivated by the earlier onset of open cell formation on the red trajectory (Figure 1), which indicates the presence of stronger precipitation and aerosol reduction by wet scavenging, and hence a more depleted aerosol population compared to the blue and green trajectories. The aerosol size distribution is described in Sec. 2.2.3.

Figure 4

Cloud drop effective radius (r eff) and optical depth (τ) in the simulations along the green (a), (b), blue (c), (d), and red (e), (f) trajectory. SEVIRI retrievals at the trajectory locations are shown in magenta (ORAC, Thomas et al., 2009) and black (Peers et al., 2019; Peers et al., 2021), with squares representing the mean and whiskers the lower and upper standard deviation. Simulations are listed in Table 1. Gray shading indicates nighttime.

Satellite Data

We use SEVIRI measurements processed with the algorithm developed by Peers et al. (2019), and with the Optimal Retrieval for Aerosol and Cloud (ORAC, Thomas et al., 2009) algorithm. The cloud retrievals are aggregated at hourly intervals over a 1° × 1° region that moves along the Lagrangian trajectories, and the mean and standard deviation of the aggregated data are used to evaluate the simulations. For more information on the data extraction process, see Christensen et al. (2020).

The Peers et al. (2019) algorithm accounts for absorbing aerosols located above clouds. The presence of absorbing aerosol above clouds has a small effect on retrieved cloud drop effective radius r eff, but cloud optical depth τ is underestimated by 35% when ignoring above‐cloud aerosol (Peers et al., 2021). The retrieved cloud properties are only weakly sensitive to assumptions on the properties of the absorbing aerosol, with biases lower than 6% in τ and 3% in r eff. The retrieved cloud properties match well MODIS retrievals and in‐situ measurements from the CLARIFY field campaign (Peers et al., 2021).

The ORAC algorithm uses an optimal estimation technique applied to two visible (0.64 and 0.84 μm), two near infrared (1.6 and 3.9 μm) and seven infrared channels (6.2, 7.3, 8.7, 9.7, 10.8, 12.0 and 13 μm) to retrieve r eff and τ at the native resolution of the SEVIRI instrument (3.5 km at nadir). The r eff retrievals operate on the 1.6 μm band, the τ retrievals use the visible channels. ORAC provides top and bottom of atmosphere broadband radiative fluxes that were recently used in aerosol‐cloud interaction studies (Christensen et al., 2017; Neubauer et al., 2017), and is described in detail in Sus et al. (2018) and McGarragh et al. (2018). ORAC has been evaluated with ground‐based measurements (Stengel et al., 2020) and the top of atmosphere fluxes agree to within 3%. The uncertainty under ideal conditions, for example, unbroken closed‐cell stratocumulus cloud decks, in droplet r eff and τ is approximately 30%. Uncertainties are considerably larger in broken cloudy conditions due to issues involving three‐dimensional radiative transfer and photon leakage out of the sides of clouds (Coakley et al., 2005). The main difference in the retrieved cloud properties between the ORAC applied to SEVIRI and MODerate Resolution Imaging Spectroradiometer collection six products (MODIS, Platnick et al., 2017) from satellites Terra and Aqua is the broader range in solar and satellite zenith angles, as well as the broader range covered by the lookup tables used for cloud drop effective radius (5 for MODIS; 1 ORAC) and cloud optical depth (τ < 100 MODIS; τ < 250 ORAC). The two products broadly agree, particularly for homogenous low‐level stratocumulus cloud layers. ORAC does not account for the effect of absorbing aerosol located above clouds. More information on the ORAC cloud retrieval algorithm is given in Sus et al. (2018) and McGarragh et al. (2018).

In‐Situ Data

We use in‐situ data collected during the CLARIFY flight C052 on its profiles P1–P7 (Abel et al., 2020). P1 to P7 sampled the open cell region within the POC (Figure S2 in Supporting Information S1). P1 was a descent from 7,150 m altitude to 35 m above the sea‐surface, enabling both the free‐tropospheric biomass burning aerosol layer and the boundary layer to be characterized. Profiles P2 to P7 then measured the boundary layer vertical profile on a track through the POC, sampling from altitudes close to the surface to the lower free‐troposphere above the trade‐wind inversion.

Simulations on the green trajectory enclose the CLARIFY flight C052 profile P7, approximately on September 5, 2017 16:30:00 UTC. We evaluate simulations on the green trajectory with boundary layer temperature, water vapor, aerosol concentration, and hydrometeor properties at its intersect with profile P7. We also evaluate the simulations with hydrometeor properties aggregated over the profiles P1 to P7. Measurements taken along the profile P1 are used to evaluate biomass burning aerosol concentrations in the free troposphere.

Results

Cloud State Transition From Closed‐to Open Cells

Figure 2 shows snapshots of the cloud state in the simulations G1, B1, and R1, at the time and locations of the satellite imagery in Figure 1. The simulated cloud state evolution is also shown in Movie S1. The cloud deck starts out overcast on September 4, 2017, 06:00:00 UTC (Figures 2a, 2e, and 2i) in all three simulations, and transitions into a broken, open‐cell state. The transition takes place at a different time on each trajectory: it occurs the latest in G1, earlier in B1, and the soonest in R1. The cloud deck is homogeneous on September 4, 2017, 06:00:00 UTC in G1 and B1 (Figures 2a and 2e), while in R1, individual locations with elevated cloud optical depth dot the cloud field (Figure 2i), indicating localized cloud thickening, drizzle, and the onset of the transition. By September 4, 2017, 14:30:00 UTC, cloud breakup has set in (Figures 2b, 2f, and 2j). Open cells are present in all three simulations on September 5, 2017 (Figures 2c, 2d, 2g, 2h, 2k, and 2l). Visual comparison of these snapshots (Figure 2) with the cloud deck at the corresponding locations in the satellite imagery (Figure 1) shows that the simulations match the observed cloud state evolution, including the timing of the transition on the three the trajectories. This is most evident in the satellite image of September 4, 2017, 14:30:00 UTC, when the cloud deck is mostly overcast on the green trajectory, partly broken on the blue trajectory, and fully broken on the red trajectory (Figure 1b), as in the simulations (Figures 2b, 2f, and 2j).

Figure 2

Optical depth of liquid water in simulation G1, green trajectory (a–d), B1, blue trajectory (e–h), and R1, red trajectory (i–l), at the times and locations shown in Figure 1.

Figure 3 shows the time series in the simulations G1, B1, and R1 from September 4, 2017, 06:00:00 UTC (d = 246.25). The different timing of the transition from closed‐to open cells between the trajectories is evident in cloud fraction (Figure 3a), rain water path (Figure 3c), and surface precipitation (Figure 3d). Cloud fraction drops and rain water path and surface precipitation rise the earliest in simulation R1, and the latest in G1. The early onset of the transition in simulation R1 is caused by its lower initial boundary layer aerosol concentration (Section 2.3). A lowered aerosol concentration at the outset of the simulation may arise for meteorological reasons farther upstream, such as a moister boundary layer with a higher liquid water path and enhanced wet scavenging. It may also be caused by variability in aerosol itself, without a contribution from meteorology. Simulation G1 and B1 have identical initial aerosol concentrations in the boundary layer (Table 1), yet the transition is delayed in G1 relative to B1, consistent with the satellite imagery (Figure 1). It is hence meteorology that determines the timing of the transition in G1 and B1, a hint that ERA5 may capture spatial variability in meteorology that drives the formation of this POC. The different timing of the transition is also apparent in aerosol (Figure 3e) and cloud and rain drop (Figures 3f and 3g) number concentrations: On the green trajectory aerosol removal by cloud scavenging is slowest, resulting in higher aerosol and cloud drop concentrations throughout the simulation (Figure 3e). Faster cloud scavenging on the blue trajectory results in lower aerosol concentrations, and the lowest aerosol concentrations are present on the red trajectory.

Figure 3

Time series in the simulations G1 (green), B1 (blue), and R1 (red). Gray shading indicates nighttime.

Evaluation With Satellite Data

Figure 4 evaluates the simulated r eff and τ with SEVIRI Peers and ORAC retrievals along the trajectories. The r eff in the simulations was calculated as the ratio of the third and second moments of the hydrometeor size distribution, after the moments were averaged over one optical depth at cloud top at locations where τ ≥ 1. τ in the simulations was calculated from the hydrometeor size distribution at locations where cloud τ ≥ 1. This section evaluates the simulations G1, B1, and R1, which use bin microphysics. The other simulations shown in Figure 4 are discussed in Sections 3.4, 3.5, and 3.6.

The simulations G1, B1, and R1 capture the evolution of r eff (Figures 4a, 4c, and 4e) and τ (Figures 4b, 4d, and 4f) retrieved by the satellite instrument over the two day simulation period. r eff evolves from smaller values on the first day to larger values on the second day. This increase of hydrometeor size reflects the transition from a non‐precipitating closed‐cell state on the first day to a broken, precipitating open‐cell state on the second day. The simulations show a daytime dip in τ on the first day. This daytime dip is driven by a combination of insolation warming and precipitation. On the second day, when the cloud deck is in the open‐cell state, the observed τ assumes very low values.

The simulations G1, B1, and R1 are in overall good agreement with the SEVIRI Peers r eff and τ, but exhibit biases and mismatches due to model and retrieval uncertainties. On the first day, when the cloud deck is in the closed‐cell stratocumulus cloud state, the simulated r eff and τ are consistent with the SEVIRI Peers retrievals, but r eff is biased low in R1 (Figure 4e), and τ in G1 (Figure 4b) and B1 (Figure 4d). The likely cause is the finite vertical resolution and associated numerical diffusion in the simulations, which causes spurious entrainment drying across the strong inversion of the closed‐cell stratocumulus cloud state, thereby reducing hydrometeor size and mass. On the second day, when the cloud deck is in the open‐cell stratocumulus cloud state, the simulated r eff and τ are in very good agreement with the SEVIRI Peers retrievals, except in R1 (Figure 4e), when the SEVIRI Peers retrieval gives very high r eff values, up to 60 μm. These high values may be an artifact of the data filter used by the algorithm, which rejects pixels identified as partly cloudy and/or associated with cloud edges, and heterogeneous clouds in the SEVIRI data aggregated at 0.1°×0.1° resolution (Peers et al., 2021). Such pixels would be associated with smaller r eff compared to fully cloudy pixels, and their rejection would result in an overestimation of r eff. The ORAC algorithm uses a data filter that captures more broken clouds, and produces smaller r eff values compared to the very high Peers values on the second day (Figure 4e). The simulation R1 is in line with the smaller r eff values produced by ORAC.

The ORAC algorithm does not account for absorbing aerosol above clouds which were present in the free troposphere during CLARIFY flight C052. The ORAC retrieval gives very similar r eff values as the Peers retrieval on the first day of the simulations, when the cloud deck is in the closed‐cell stratocumulus clouds state (Figures 4a, 4c, and 4e). The Peers retrieval gives generally higher values on the second day, when the cloud deck is in the open‐cell stratocumulus cloud state. However, above‐cloud absorbing aerosol has only a small effect on retrieved r eff values (Peers et al., 2021). The higher r eff values of the Peers retrieval on the second day may hence in general be caused by its data filter, discussed above. The ORAC retrieval gives systematically lower τ values than the Peers retrieval (Figures 4b, 4d, and 4f). This low bias is small on the first day and larger on the second day, when it assumes values that are by and large consistent with an underestimation of 35% caused by ignoring above‐cloud aerosol (Peers et al., 2021).

A comparison of simulated r eff and τ with MODIS collection 6 products (Platnick et al., 2017) is shown in Figure S3 in Supporting Information S1 for completeness. r eff and τ in G1, B1, and R1 agree overall well with the MODIS retrieval, with similar biases as seen relative to the SEVIRI data.

Evaluation With In Situ Data

Figures 5, 6, 7 evaluate the simulation G1 at the location where its trajectory crosses the path of CLARIFY flight C052 (Abel et al., 2020). The location of the simulation domain and of the CLARIFY flight C052 profiles P1–P7 that provide measurements are shown in Figure S2 in Supporting Information S1.

Figure 5

Domain mean profiles in simulation G1 (green) on September 5, 2017, 16h45m00s UTC, ERA5 profiles (black), interpolated to the vertical grid of the simulations, on September 5, 2017, 16h32m30s UTC, and in‐situ measurements from CLARIFY flight C052. CLARIFY C052 profile P7 (5 September 2017, 16:27:49–16:36:33 UTC, magenta) is located at the intercept of the simulation trajectory and the CLARIFY C052 flight path. CLARIFY flight C052 profile P1 (September 5, 2017, 15h25m18s–15h50m53s, gray) is located upstream along the flight path. The location of the simulation domain, the CLARIFY flight C052 path and its profile P1 and P7 are shown in Figure S2 in Supporting Information S1.

Figure 6

Cloud properties (a), (b), (c) and frequency of occurrence of cloud and rain water (d). Simulation G1 (green) is shown on September 5, 2017, 16h45m00s UTC, at the intercept of the simulation trajectory with the path of CLARIFY flight C052. Measurements at the intercept (magenta) were taken during the C052 flight profile P7 (September 5, 2017, 16:27:49–16:36:33 UTC). Measurements from C052 flight profiles P1–P7 (median and interquartile range, yellow/black dots with whiskers) represent a longer flight segment, extending upstream of the intercept (5 September 2017, 15h44m10s–16h39m41s UTC). CLARIFY flight C052 profiles P1–P7 mean cloud (gray) and rain (black) fractions are shown in panel (d) The location of the simulation domain at the intercept with CLARIFY flight C052, and locations of the profiles P1–P7 are shown in Figure S2 in Supporting Information S1.

Figure 7

Same as Figure 6 but rain properties.

Aerosol, Temperature, and Water Vapor

Figure 5 compares aerosol concentrations, temperature, and water vapor in simulation G1 with CLARIFY C052 data. Simulation G1 accurately reproduces the observed aerosol profile in the boundary layer (below about 1,800 m), including its slight negative gradient with altitude which arises from sea surface emissions and depletion in the cloud layer. The depletion of boundary layer aerosol by cloud processes is evident in the low observed aerosol concentration (Figure 5a, magenta) at the location of the simulation domain (CLARIFY profile P7), where open cells have existed for a longer time, relative to the less depleted aerosol concentration (Figure 5a, gray) farther upstream (CLARIFY profile P1), where open cells have formed more recently. The simulation matches the vertical distribution of the biomass burning layer above the inversion, with the exception of a “bite‐out” in the lower free troposphere around 2,200 m, and peak aerosol concentrations around 3,200 m (Figure 5a).

The simulated temperature (Figure 5b) and water vapor (Figure 5c) profiles reproduce well the qualitative features of the observations, with quantitative biases. The simulated inversion misses the observed inversion by only 100 m, but the boundary layer has a cold and dry bias. Under the assumption that simulated and observed variability in these quantities is comparable, this is unlikely due to a sampling bias at levels where the observed profiles are outside the 1st–99th percentile range of the simulated values (Figures 5b and 5c). The assumption may not apply, since the observations, which sample a limited volume of the boundary layer, show variability that is comparable or greater than variability over the entire simulation domain. Hence variability in the simulations may be too small and the observed profiles not statistically representative.

The slightly low inversion, and the cold and dry bias in the boundary layer may be caused by a horizontal grid spacing that is too coarse. A finer horizontal grid spacing would reduce numerical diffusion of vertical momentum and hence strengthen turbulence. This would enhance mixing and reduce the dry bias in the upper boundary layer by transporting moisture from the surface to higher levels. The stronger turbulence would also drive entrainment, lift the inversion, and warm the boundary layer. This response to a refinement of grid spacing would in part be offset by adjustments in surface fluxes, cloud water content, and radiative heating and cooling. Simulation grid effects are discussed in more detail in Section 4.1.

ERA5 places the inversion at the observed height, but has a very strong cold (Figure 5b) and dry (Figure 5c) bias in the upper boundary layer. Since it reproduces both temperature and water vapor well near the surface, the cold and dry bias farther aloft may arise from insufficient boundary layer turbulence and mixing. ERA5 has a cold and moist bias in the lowermost free troposphere, at about 2,000 m, relative to the observations (Figures 5b and 5c). The warmer temperature in the observed lowermost free troposphere may be caused by heating from absorption of radiation by the biomass burning aerosol layer, which may not not be fully captured in ERA5 via data assimilation. The ERA5 biases in the lowermost free troposphere may propagate into our simulations by affecting inversion stability, and entrainment of heat and moisture into the boundary layer.

Cloud and Rain Properties

Figures 6 and 7 compare cloud and rain properties in simulation G1 with measurements from CLARIFY flight C052. The simulation results are evaluated at the intercept of the simulation trajectory with the path of flight C052 with measurements from that location (flight profile P7), and with measurements from a longer flight segment that extends upstream of the intercept (flight profiles P1–P7). The simulation results were processed to emulate the sampling by the aircraft instruments using thresholds given in Abel et al. (2020): Cloud properties represent hydrometeors up to 25 μm radius, sampled from locations where their liquid water content is 0.01 gm−3, and such locations contribute to the calculation of cloud fraction. Rain properties represent hydrometeors of 30 μm in radius or larger, sampled from locations where their concentration exceeds 1 L−1, and such locations contribute to the calculation of rain fraction.

Simulated cloud properties are consistent with the observations at the intercept of the simulation trajectory with the path of CLARIFY flight C052 (flight profile P7, Figures 6a–6c). Large scatter in the measurements arises from the profile P7 extending almost across the simulation domain, sampling different locations and cloud elements in the cloud field (Figure S2 in Supporting Information S1). Despite the scatter in the observed cloud water (Figure 6a) and cloud drop number (Figure 6b), the model and observations show good agreement in the vertical structure of cloud drop mean volume radius, which increases with height (Figure 6c).

Simulated rain water and rain drop number are in reasonable agreement with the observations at the intercept, given the significant scatter in the measurements (Figures 7a and 7b). Observed rain drop mean volume radii are often around 35 μm, smaller compared to the simulated values across the domain (Figure 7c). Simulated rain rates are consistent with measured values (Figure 7d), with the caveat of large scatter in the observations that likely arises from the sampling of different locations in the cloud field.

Since the aircraft profile P7 represents only a very small sample volume relative to the simulation domain, we compare the simulation with observations aggregated over the flight profiles P1–P7 (Figures 6 and 7). Simulated cloud water and cloud drop number are in good agreement below 1,250 m (Figures 6a and 6b). Above 1,250 m, observed values of cloud water and cloud drop number are systematically higher than the simulation results, but the distributions of simulated and observed data exhibit a large overlap. The higher observed values can have several causes. The leading potential cause is that the clouds sampled on the profiles P1–P6 are systematically richer in liquid water and cloud drops compared to the clouds sampled on the profile P7 at the simulated location. This is supported the SEVIRI 1.6 μm channel image at time of the profiles P1–P7 (Figure S2 in Supporting Information S1), which shows that the flight segment upstream (east‐southeast) of profile P7 crossed brighter cloud elements, indicating higher liquid water than is present on flight profile P7, or within the simulation domain. In agreement, measurements upstream (farther east) of profile P7 show higher values of cloud water and cloud drop number (Figure S4 in Supporting Information S1). Differences between the simulated cloud water and cloud drop number and the measurements along the profiles P1–P7 hence likely arise from different conditions and cloud state along the flight segment upstream of the simulated location. Despite these differences, simulated cloud drop mean volume radii are in good agreement with the observations along the profiles P1–P7 (Figure 6c). The agreement holds up to the cloud top region, which is consistent with the good agreement between the simulated cloud drop effective radii in G1 and SEVIRI measurements at the corresponding time, d = 247.69 (Figure 4a).

Simulated rain water and rain drop number are in good agreement with the observations aggregated along the CLARIFY flight C052 profiles P1–P7, with observed values generally higher than the simulations at altitudes above 1,500 m (Figures 7a and 7b). The simulation captures well the vertical structure in the observed rain drop number, which shows lower values near the surface and higher values near cloud top (Figure 7b). The simulation closely reproduces observed rain drop mean volume radii (Figure 7c) and rain rates (Figure 7d), both in terms of absolute values and vertical structure.

Figure 6d shows the cloud and rain fractions from simulation G1 and along the CLARIFY flight C052 profiles P1–P7. Observed values are larger than the simulated values for both cloud and rain fraction. Visual inspection of the satellite imagery shows that the aircraft frequently crossed cloudy areas upstream (east‐southeast) of the simulation domain (Figure S2 in Supporting Information S1). Frequent cloud encounters by the aircraft upstream (east‐southeast) of the simulation domain are documented in measurements of cloud water and cloud drop number (Figure S4 in Supporting Information S1). We hence attribute the higher observed cloud and rain fractions to different conditions and cloud state along the flight segment upstream of the simulated location. Despite these differences, simulation and observations agree on a higher rain fraction compared to the cloud fraction.

In summary, the simulation is in reasonable agreement with the in‐situ observations of cloud and rain water mass and number, with the main limitation of the evaluation being the sparseness of the observations relative to the variability in the cloud field on the scale of the simulation domain, and the associated scatter. To reduce the uncertainty in the evaluation of the simulation from this scatter, we compared the simulation with observations aggregated over a longer flight segment. The simulation matches this larger sample better, although in the upper boundary layer, it systematically underestimates cloud water mass and number, and to some extent rain water mass and number. These biases are likely caused by different conditions and cloud state along the longer flight segment compared to the simulated location. The simulation does well in reproducing profiles of cloud and rain drop mean volume radii, rain water and rain drop number, and rain rate. This, together with the good agreement with the satellite cloud optical depth and cloud drop effective radii (Sec. 3.2), indicates that the model performs well, and that disagreements with the in‐situ data are largely due to in‐situ undersampling and different conditions at the simulated and observed locations.

Impact of the Free Tropospheric Biomass Burning Layer Observed During CLARIFY Flight C052

Based on measurements during CLARIFY flight C052, Abel et al. (2020) found weaker entrainment of free‐tropospheric biomass burning aerosol into the boundary layer of the underlying POC compared to the surrounding closed‐cell cloud deck. We hence tested the impact of entrainment of aerosol from the free‐tropospheric biomass burning layer on cloud‐ and boundary layer properties underneath.

Simulation G0 and G1 have identical setups, but simulation G1 is initialized with a biomass burning layer in the free troposphere and simulation G0 without it (Table 1 and Figure S5a in Supporting Information S1). The two simulations produce nearly identical evolutions of r eff and τ (Figures 4a and 4b). The time series of other cloud‐ and boundary layer properties are also nearly identical except for variability in the rain water path and surface precipitation on short time scales (Figure S6 in Supporting Information S1) At the time of the intercept of the simulation trajectory with the path of flight C052 with measurements from that location (flight profile P7), the simulations produce identical profiles of aerosol, temperature, and water vapor, except in the free troposphere, where simulation G0 exhibits the free tropospheric background aerosol concentration, while simulation G1 matches the aerosol concentrations observed in biomass burning layer above the inversion (Figure S5b–S5c in Supporting Information S1). We conclude that entrainment of aerosol from the biomass burning layer overlying the POC sampled during CLARIFY flight C052 is limited to the extent of having no impact on cloud‐ or boundary layer properties. This is in agreement with the observations of Abel et al. (2020).

Bin‐ and Bulk Microphysics

Time series of r eff and τ calculated with the bin (G1, B1, R1) and the bulk (G2, B2, R2) cloud microphysics scheme (Section 2.2.2) are shown in Figure 4. The r eff time series shows no systematic difference between the two microphysics schemes (Figures 4a, 4c, and 4e). Systematic differences are present in the τ time series (Figures 4b, 4d, and 4f): the bulk scheme produces high values starting at night and over the course of the second day. The elevated τ values are due to a higher rain water path in the bulk scheme, which dominates liquid water path, and is caused by lower surface precipitation during the night and the second day in the simulations (Figure S8, S9, and S10 in Supporting Information S1).

Figure 8 shows the hydrometeor mass distribution averaged over the last 12 hr of simulation G1 (bin microphysics) and simulation G2 (bulk microphysics), in the center of the boundary layer and at the surface. The cloud deck is in the open‐cell state at this stage of the simulations. The bin microphysics produces a rain mode with hydrometeors that are approximately one order of magnitude larger than those in the rain mode of the bulk microphysics. Between the center of the boundary layer and the surface, the rain mode moves to larger sizes in the bin microphysics and to smaller sizes in the bulk microphysics. The larger hydrometeors in the rain mode of the bin microphysics and their faster fall speeds are responsible for a larger precipitation flux near the surface compared to the bulk microphysics (Figure 9a). In turn, less liquid water is retained in simulation G1 and more in simulation G2 (Figure 9b).

Figure 8

Hydrometeor mass distribution, averaged over the last 12 hr (d = 247.2 to 247.7), in simulation G1 (green) using bin microphysics, and simulation G2 (dark green) using bulk microphysics, in the center of the boundary layer (solid) and at the surface (dotted).

Figure 9

Vertical profiles, averaged over the last 12 hr (d = 247.2 to 247.7) of simulation G1 (green) using bin microphysics, and simulation G2 (dark green) using bulk microphysics.

Hence when precipitation is present, the bulk microphysics overestimates liquid water compared to the bin microphysics, and produces a high bias in τ relative to the SEVIRI observations, because of an underestimation of surface precipitation due to insufficient formation of large rain drops. In the bulk microphysics, the largest hydrometeors in the rain mode evaporate as they travel towards the surface, in contrast to the bulk microphysics (Figure 8). The evaporation of rain with the bulk microphysics results in a moist (Figure 9c) and cold (Figure 9d) bias in the lower regions of the boundary layer relative to the bin microphysics.

Initial Aerosol Concentration

On the red trajectory, a lower boundary layer aerosol concentration is used compared to the green and blue trajectories to initialize simulations. The rationale for the lower value on the red trajectory is given in Section 2.3. Here we expand on this rationale by comparing simulation R1, which uses the lower initial aerosol concentration of 115 mg−1, with simulation R3, which uses the higher value of 145 mg−1, in the context of the SEVIRI satellite observations. Time series of r eff and τ in simulations along the red trajectory are shown in Figures 4e and 4f. Simulation R3 underestimates the observed r eff on the first day, while simulation R1 produces larger r eff values that match the observations better (Figure 4e). No distinguishing difference exists between the simulations in terms of r eff on the second day. τ is consistent in both simulations with the observations on both days (Figure 4f). The better agreement of R1 compared to R3 with the observed r eff on the first day supports the use of a lower initial boundary layer aerosol concentration for simulations on the red trajectory.

Role of Simulation Setup

We determine key elements of the simulation setup for the ability of the simulations to reproduce the observations. These elements have in common that they act on the surface flux and vertical transport of water vapor, and thereby alter the thermodynamic properties of the boundary layer. The analysis uses simulations with the bulk microphysics scheme.

Vertical Grid Spacing and Ventilation of the Surface Layer

The simulations in this work employ a vertical grid that coarsens towards the surface in the lowermost three levels. To illustrate its effect, we compare simulation G2, which uses the grid coarsening towards the surface, with simulation G3, which uses a constant grid down to the surface. The constant grid has a finer spacing near the surface (Table 1 and Section 2.2.6).

The surface latent heat flux is lower in G3 compared to G2, and G3 has a drier boundary layer (Figures 10a and 10b). Yet, G3 has a moister surface layer compared to G2 (Figures S10a, S10b, S11a, and S11b in Supporting Information S1), indicating suppressed surface ventilation. The surface sensible heat flux, in contrast, is nearly identical in G2 and G3 during the first daytime period of the simulations (Figure 10c), and the surface layer in G3 is warmer by only fractions of a degree during this time (Figures S10 c and S11 c in Supporting Information S1). We hypothesize that surface layer warming due to a weaker ventilation of the surface in G3 is offset by longwave radiative cooling.

Figure 10

Time series in simulation G2 (dark green), G3 (aquamarine), and G4 (beige). Gray shading indicates nighttime.

Following the first daytime period, the surface sensible heat flux (Figure 10c) rises sooner in G2 (at nightfall) compared to G3 (towards dawn), because in the moister boundary layer in G2, rain and surface precipitation form sooner (Figures 10f and 10g). The associated evaporation results in a greater temperature contrast near the surface between G2 and G3 (Figures S10d and S11d in Supporting Information S1). The boundary layer as a whole, however, is warmer in G2 compared to G3 (Figure 10d), owing to its higher surface precipitation.

Hence a constant, refined vertical grid near the surface hampers ventilation of the surface layer, causing it to moisten. This in turn suppresses the surface latent heat flux, which alters the thermodynamic properties of the boundary layer: The boundary layer is drier and surface precipitation delayed and suppressed. This results in a warmer surface layer but a generally cooler boundary layer. Overall, the refined vertical grid near the surface causes a dry and cold bias in the boundary layer. Coarsening the vertical grid towards the surface reduces this effect.

Wind Speed Nudging Near the Surface

Mean horizontal wind speed is maintained in the simulations by nudging towards ERA5 wind speed profiles. To counter deceleration by surface drag away from the ERA5 wind speed, the nudging time constant tightens towards the surface (Section 2.2.1). To illustrate the impact, we compare results obtained with the tighter nudging towards the surface (simulation G3) against results obtained with more relaxed nudging at all levels (simulation G4). With the relaxed nudging, wind speed near the surface is slower in G4 compared to G3 (Figure S12 in Supporting Information S1). The surface latent heat flux falls in response, which renders the boundary layer drier (Figures 10a and 10b). The surface sensible heat flux is largely insensitive to the surface wind speed reduction, but the boundary layer in G4 is cooler compared to G3 (Figures 10c and 10d). This response arises from the complex interactions connecting dynamics, surface fluxes, cloud state, radiative cooling, and entrainment warming.

Figure S10b, and S10d in Supporting Information S1 shows that simulation G4 has a larger dry and cold bias, and a stronger low inversion bias compared to simulation G3, relative to the in‐situ observations. This supports the use of a short time scale for mean horizontal wind speed near the surface. In summary, relaxed nudging of the mean horizontal wind speed towards ERA5 reduces surface wind speed and causes a dry and cold bias in the boundary layer. This can be counteracted by tighter nudging towards the surface.

Discussion

The cold and dry bias in the simulated boundary layer relative to the in‐situ observations found in this work indicates remaining model uncertainties and potential for improvement. Areas that contribute to these uncertainties, and where improvements are possible, are discussed in the following, with comments on future high‐resolution global models.

Grid Anisotropy

Non‐isotropic grids with large grid aspect ratios enable large simulation domains that capture the stratocumulus mesoscale structure and its effect on cloud properties (Kazil et al., 2017). In this work we used an aspect ratio of 20 across the boundary layer, except towards the surface, where it falls to 5.7. This reduction of the grid aspect ratio towards the surface, implemented by a coarsening of the vertical grid spacing, improves surface ventilation and reduces a boundary layer cold and dry bias (Section 3.7.1). A cold and dry bias was also found by Vogel et al. (2020) in LES of shallow cumuli with an aspect ratio of 7.8 near the surface. This raises the question whether in general, large grid aspect ratios near the surface should be avoided.

Nishizawa et al. (2015) investigated the role of LES grid aspect ratio at fixed surface heat flux for turbulence in the dry boundary layer. In their simulations, which reduced the grid aspect ratio from 20 to 10 and from 6 to 2, implemented by coarsening the vertical grid spacing at all levels, the vertical component of grid‐resolved TKE increased at all scales in the surface layer. This translates to better surface ventilation. Higher up in the boundary layer, the reduction in aspect ratio had only a small effect on the vertical component of grid‐resolved TKE. These findings provide more general support for improving surface ventilation by using smaller grid aspect ratios towards the surface.

The dependence of LES results on the grid aspect ratio represents an uncertainty that will, over time, diminish as increasing computing power enables smaller grid aspect ratios on large domains. Concurrently, the issue will arise in global models as increasing computing power enables finer grids spacings. Once grid spacings are too fine for boundary layer parameterizations to be applicable, and to compensate a suppression of surface ventilation by a large grid aspect ratio, coarsening the vertical grid spacing towards the surface, as done in this work, may offset artifacts.

The subgrid scale TKE scheme (Section 2.2.5) used in our simulations does not account for different grid spacings in the horizontal and vertical dimension. Sub‐grid scale turbulence parameterizations that account for grid anisotropy may constitute a better approach. Nishizawa et al. (2015) demonstrated the importance of parameterizing the LES mixing length as a function of grid aspect ratio, and of using an appropriate LES filter length to isolate sub‐grid scales in order to obtain the theoretical scaling of TKE with wavenumber. Aspect ratio‐aware sub‐grid scale parameterizations may help reduce or eliminate the dependence of results on the grid aspect ratio in LES, and later serve in high‐resolution global simulations with non‐isotropic grids.

Aerosol and Cloud Microphysics

Simulated cloud properties are sensitive to the representation of the aerosol size distribution and of activation and condensational growth (e.g., Feingold & Kreidenweis, 2002). The simplified representation of the aerosol size distribution used in this work (Section 2.2.3) may affect its response to activation and cloud processing, with potential consequences for subsequent activation and cloud microphysical processes. Representing the aerosol size distribution with a bin scheme, for example, could enable a more detailed and potentially more accurate response of the aerosol size distribution to activation and cloud processing. However, we expect the uncertainty due to the representation of the aerosol size distribution to be relatively small owing to the buffering of aerosol activation by supersaturation, i.e., overactivation suppresses supersaturation, which self‐corrects the strength of activation.

The two‐moment bin microphysics scheme used in our simulations performs well relative to the observations (Sections 3.2 and 3.3), but the remaining biases and deviations relative to the observations could, potentially, arise from its limitations. One such limitation is the artificial broadening of hydrometeor size distributions. The use of a two‐moment bin scheme reduces the broadening considerably but is not immune to numerical diffusion artifacts (Witte et al., 2019). The broadening arises from numerical diffusion caused by the remapping of the hydrometeor size distribution after growth and collisions (see, for example, Khain et al., 2008, and references therein), and due to numerical diffusion associated with advection (Morrison et al., 2018). Morrison et al. concluded that Eulerian dynamical models, such as most LES using bin microphysics, may be unable to investigate the physical mechanisms for size distribution broadening, even though they may reasonably simulate overall size distribution characteristics. More advanced representations of the hydrometeor size distribution and processes could identify and reduce or eliminate potential artifacts. Lagrangian cloud microphysics schemes (e.g., Grabowski et al., 2018), in combination with a linear eddy model to represent unresolved turbulent mixing at the sub‐grid scale of LES (Hoffmann et al., 2019) can eliminate issues affecting other microphysics schemes and their calculation of droplet activation and growth. However, such highly accurate solutions will remain computationally too expensive for some time for typical LES applications, and even more so for climate models. This represents an opportunity for further research and development of microphysics schemes that aim for reduced artifacts at manageable computational expense, such as three‐moment bulk schemes (e.g., Paukert et al., 2019).

Large Scale Meteorology

Improved understanding of the role of assumptions and methods used in the implementation of large scale meteorology in Lagrangian LES driven by reanalysis meteorology may improve the approach. In this work, the mean LES temperature and moisture profiles in the free troposphere are nudged towards the reanalysis. An alternative is the application of tendencies of temperature and moisture due to horizontal advection from the reanalysis to the LES temperature and moisture profiles. In the boundary layer, the tendencies will vanish in good approximation as the LES domain moves with the boundary layer air mass. Around the inversion and above, these tendencies will be different from zero and could be used instead of nudging. Using horizontal advective tendencies instead of nudging would allow the LES radiation scheme to act on temperature in the free troposphere. This would, for example, enable the study of the effect of absorbing aerosol layers in the free troposphere. This approach could, however, overestimate heating by the absorbing aerosol, as its effect could already be partially present in the horizontal advective temperature tendency from the reanalysis, as a result of data assimilation. A further potential downside of using tendencies instead of nudging is that differences between the radiation schemes in the LES and reanalysis model may lead to inconsistent free tropospheric temperature profiles between the LES and reanalysis, with potential consequences for the LES results.

The mean horizontal wind speed in the simulations in this work is nudged towards the reanalysis wind speed with a short nudging time scale near the surface, to offset spurious slowing by surface drag and to drive appropriate surface fluxes. Higher up, a longer nudging time scale is used to allow the LES to establish its own wind speed profile around the inversion, as opposed to being forced by the wind speed profile around the inversion in the reanalysis. Still, shear in the mean horizontal wind speed around the inversion in the reanalysis may affect the mean wind speed profile in the LES and possibly lead to artifacts, in particular if the inversion in the reanalysis is located at a different height than the inversion in the LES. An alternative approach is to use the horizontal pressure gradient from the reanalysis to let the LES generate its own mean horizontal wind field. This may reduce or eliminate artifacts that potentially arise from nudging towards the mean horizontal wind speed in the reanalysis.

A need for development is present in the treatment of tracer advection due to subsidence when tracers have strong vertical gradients, such as aerosol layers. The numerical treatment of advection by subsidence by the model used in this work preserves shape, but not mass. We conserve tracer mass by re‐normalizing its vertically integrated value in the free troposphere after subsidence is applied (Section 2.2.5). A better solution is needed in the form of an advection scheme that maintains both the shape and mass of free‐tropospheric tracer distributions against advection by subsidence.

The reanalysis meteorology that drives the Lagrangian LES is itself a source of uncertainty. ERA5 performs better relative to its predecessor ERA‐Interim, and ERA5 deviations from observations just prior to their assimilation are decreasing over the reanalysis period. Yet, for example, the 30‐day mean of the ERA5 standard deviation from observed 2 m relative humidity just prior to its assimilation is 9–10% in 2017, the year of our simulations (Hersbach et al., 2020). Larger uncertainty should be expected at locations where observations are not assimilated, on shorter time scales, and in quantities that are not constrained by data assimilation. Subsidence, for example, has been found to exhibit large variability among different reanalyses as well as biases relative to observations (Uma et al., 2021). However, the overall good agreement of our simulation results with the observations indicates that ERA5 characterizes large scale meteorology well in the considered case.

Summary and Conclusions

In this work we presented and evaluated an approach to drive Lagrangian LES with reanalysis meteorology to simulate boundary layer clouds. The simulation domains follow trajectories of the boundary layer flow. The simulated case is a sub‐tropical transition from a closed‐to an open‐cell stratocumulus cloud state over a period of two days, which occurred during the formation and evolution a pocket of open cells (POC) underneath a free‐tropospheric biomass burning aerosol layer.

The simulations were evaluated with retrievals of cloud optical depth τ and cloud drop effective radius r eff from the SEVIRI instrument on board the MSG satellite (Christensen et al., 2020; Peers et al., 2019, 2021), and with aircraft in‐situ measurements from the CLARIFY field campaign (Abel et al., 2020). The simulations reproduce the observed cloud morphology, τ, and r eff observed by the satellite in the overcast, closed‐cell stratocumulus cloud state on the first day of the simulations and in the broken, open cell state on the second day. They capture the timing of the cloud state transition from the closed to the open cell state seen in the satellite imagery on the three considered trajectories. The simulated inversion height of the open‐cell state matches the the aircraft data, but the boundary layer has a cold and dry bias relative to the in‐situ measurements.

We found two key elements in the simulation setup that contribute to the cold and dry bias of the open cell state: first, large grid aspect ratios, needed to cover large domains, suppress ventilation of the surface layer. Reducing the grid aspect ratio towards the surface by coarsening the vertical grid spacing improves surface ventilation and reduces this cold and dry bias. Second, the use of a short time scale for the nudging of mean horizontal wind speed towards the reanalysis near the surface maintains mean wind speed close to the reanalysis values. This maintains the surface fluxes of sensible and latent heat and warms and moistens the boundary layer. The remaining cold and dry bias in the simulated boundary layer likely includes contributions from the still anisotropic grid, from the treatment of cloud microphysics, and from uncertainty in the reanalysis meteorology used to drive the simulations.

The simulations closely reproduce a biomass burning aerosol layer identified by the in‐situ aircraft measurements just above the inversion of the POC, as well as the aerosol concentration in the boundary layer. Simulations with and without the biomass burning layer produce nearly identical results. We conclude that entrainment of aerosol from the biomass burning layer overlying the POC is limited to the extent of having no impact on cloud‐ or boundary layer properties. This is in agreement with observations from the CLARIFY field campaign, which found only limited entrainment of biomass burning aerosol into the boundary layer (Abel et al., 2020).

Simulated mass and number of cloud and rain are consistent with the in‐situ aircraft measurements. Simulated cloud and rain drop sizes, as well as rain rates are in good agreement with the observations. Based on analysis of the satellite imagery and the in‐situ data, we conclude that aggregation of measurements along the aircraft flight path, which reduces noise but merges data from conditions with different cloud properties, is a key contribution to differences between simulated and observed hydrometeor properties.

Simulations using a numerically efficient two‐moment bulk microphysics scheme, instead of the two‐moment bin microphysics scheme, reproduce the satellite r eff and τ in the non‐precipitating closed‐cell state of the simulations well. However, they overestimate τ in the precipitating, open‐cell state. The cause is an insufficient formation of large rain drops, which results in an underestimation of surface precipitation and overestimation of liquid water path.

In summary, we find that Lagrangian LES, driven by reanalysis meteorology, are capable of realistically simulating boundary layer clouds. Owing to its ability to reproduce real‐world cases, the approach is suited to investigate and explain observed phenomena, such as in the context of field campaigns. This potential for realism, together with a spatially and temporally highly resolved output, also renders the approach suitable as a framework for the development of process representations, such as cloud microphysics schemes, and of single column models and retrieval algorithms for remote sensing instruments.

The challenges facing Lagrangian LES driven by reanalysis meteorology, such as the dependence of the results on the grid aspect ratio, will diminish over time, possibly due to the use of improved sub‐grid scale turbulence parameterizations that account for grid anisotropy, and certainly as increasing computing power will enable smaller grid aspect ratios on large domains. These challenges will, however, with increasing computing power and finer grids, eventually arise in global models. They will be compounded by boundary layer parameterizations being applicable only as long as the grid spacings are not too fine. The development and use of Lagrangian LES driven by reanalysis meteorology can thus pave the way for the development of future global models.

Supporting Information S1

Click here for additional data file.

Movie S1

Click here for additional data file.