SUMMERTIME POST-COLD-FRONTAL MARINE STRATOCUMULUS TRANSITION PROCESSES OVER THE EASTERN NORTH ATLANTIC

OF THE DISSERTATION Summertime Post-Cold-Frontal Marine Stratocumulus Transition Processes over the Eastern North Atlantic by Melissa Kazemi Rad Dissertation Director: Mark A. Miller The Marine stratocumulus cloud system is a major component of the Earth’s energy budget. Mid-latitude stratocumulus are known to transition from a single, continuous cloud layer to a hybrid configuration that includes both stratocumulus and cumulus, and eventually to trade cumulus toward the tropics. Stratocumulus transitions are often observed in the wake of cold air outbreaks in the mid-latitude summertime marine boundary layer (MBL). Cloud morphology associated with two summertime cold fronts over the Eastern North Atlantic (ENA) is investigated using high resolution simulations from the Weather Research and Forecasting (WRF) model and observations from the Atmospheric Radiation Measurement

The Marine stratocumulus cloud system is a major component of the Earth's energy budget. Mid-latitude stratocumulus are known to transition from a single, continuous cloud layer to a hybrid configuration that includes both stratocumulus and cumulus, and eventually to trade cumulus toward the tropics. Stratocumulus   components of the planetary radiation budget that act to cool the ocean surface by significantly impeding incoming shortwave radiation due to their high albedo, while minimally affecting longwave radiation losses to space, since their cloud-top temperatures are only slightly cooler than the ocean surface beneath them. They are also potentially susceptible to global climate change. Due to their climatological importance, it is essential to understand underlying processes during various stages in their life cycle, and accurately parameterize them in regional and climate models. An extensive and variable MBL cloud system is present over the Eastern North Atlantic (ENA) that undergoes a morphological transition with latitude. This region can be qualitatively explained as the following: moderate or strong subsidence associated with a subtropical anticyclone produces adiabatic compression that leads to a sharp difference between the potential temperature of air in the lower troposphere and the potential temperature of the air immediately above the ocean surface with the former temperature being much greater than the latter. The very large heat capacity of the ocean leads to the formation of a stable layer adjacent the ocean surface that is capped by a temperature inversion. The primary factors that enhance this inversion and regulate the depth of the marine boundary layer include weak turbulence induced by wind shear near the ocean surface and cloud-top radiative cooling that produces negatively buoyant parcels, which provide the bulk of the turbulent kinetic energy required to maintain a mixed layer. Clouds over the northern reaches of the ENA tend to be layered and tend to be comprised of shallow turbulent moist stratus decks overlain by a strong temperature inversion and warm dry air, while those to the south tend to be more broken. Conceptual models of this cloud transition depict reductions in inversion strength and increasing sea-surface temperatures toward the Intertropical Convergence Zone accompanying the evolution to broken cloud structure. These large-scale changes in the upper and lower boundary conditions of the MBL are related to the latitudinal evolution of semi-permanent subtropical high-pressure systems that are part of the Hadley circulation. Given the importance of marine boundary layer clouds in the Earth's radiation budget, changes in the location or function of the ENA cloud transition region in response to anthropogenic warming would likely impart significant changes in the regional and planetary radiation budgets (Bretherton et al., 2013).
Stratocumulus, shallow cumulus, and complex combinations of these cloud types, often exhibiting mesoscale organization, are present in the transition MBL cloudiness over the ENA. Several studies have detailed the MBL cloud morphology, thermodynamic environment, and precipitation during the summertime , Miller and Albrecht (1995); Wood (2005); Dong et al. (2014a,b); Wood et al. (2015); and others), while additional studies have examined the basic life cycle of clouds through the transition process (Bretherton and Pincus (1995); de Roode and Duynkerke (1997) Changes in cloud morphology over the ENA are often related to a process known as decoupling (Nicholls, 1984), whereby the MBL is composed of two layers separated by an intermittent, weak inversion: a moist sub-cloud layer immediately above the ocean surface and a cloud layer above capped by the marine inversion . The Atlantic Stratocumulus Experiment (ASTEX) was an important study that was aimed at investigating processes responsible for the transition between solid stratocumulus and trade cumulus cloud regimes over the ENA . It was carried out during June 1992, utilizing island, aircraft, ship, and satellite measurements. The results from ASTEX highlighted the important role of decoupling in the atmospheric boundary layer in stratocumulus cloud transitions, as well as deepening of the boundary layer, solar absorption in clouds, and drizzle, all of which impact the decoupling process.
As the decoupling process takes its course, the solid deck of stratocumulus transitions into mesoscale, cellular structures that can be maintained by mesoscale convective systems that result from the accumulated moisture in the boundary layer.
A deepening-warming decoupling process first proposed by  is thought to be a key initiator of the transition from a sheet of stratocumulus to more broken cloud structures. According to this hypothesis, the cloud-base buoyancy flux increases in proportion to latent heating, while the MBL average buoyancy flux does not, leading to decoupling. At the onset of decoupling and the transition in cloud structure, the lifting condensation levels (LCLs) in updrafts and downdrafts gradually diverge and the cloud layer dries.
Some of the earliest regional simulations included a two-layer, 1-D model of the stratocumulus transition, which produced a regional gradient in the cloud structure (Wang et al., 1993), and a mesoscale model with 5-km horizontal resolution that was used to test cloudiness parameterizations (Mocko and Cotton, 1995).
Regional simulations designed to test the impacts of drizzle and microphysics on mesoscale organization showed that mesoscale structures were created in simulations using 2-km horizontal resolution when many physical processes were resolved, especially those involving drizzle production, but not when the resolution was increased to 18-km and parameterizations were required (Mechem and Kogan, 2003). A regional simulation using a mesoscale model with a resolution compatible with current climate models (60-km), which used MBL turbulence and shallow convection parameterizations, was able to successfully simulate the stratocumulus to cumulus transition (Bretherton et al. (2013); ).
A key finding of this study was that a shallow convection scheme was required to faithfully reproduce the stratocumulus to cumulus transition. A later study demonstrated that mesoscale models with coarser resolution are sensitive to the parameterization of shallow convection and to the parameterization of penetrative mixing at the top of the cumuli, which acts to deepen the trade inversion . Sensitivities in the representation of cloud-top entrainment were also reported in another regional study (Wang et al., 1993).
Cascading cold frontal cloud bands are integral components of the ENA transition region and are observed in all seasons (Reed (1960); Giangrande et al.
( 2019)). Cold fronts in this region have higher average maximum and minimum integrated water vapor contents than North Pacific fronts in the summertime, a difference that can be attributed to higher sea-surface temperatures and, consequently, higher evaporation and moisture flux in the region (McMurdie and Katsaros, 1991). A recent study by  documented the climatology of post-cold-frontal conditions over the ENA associated with 77 cold frontal passages, 15 of which occurred during the summertime. They reported that "postcold frontal periods (including all seasons) were more dynamically active, drier, and more unstable" than periods that did not involve frontal passages. Naud et al.
(2018) also found that post-cold-frontal MBL clouds are deeper, colder, and have higher cloud-tops, and possess higher liquid water paths than non-post-frontal clouds. Cloud depth was suggested in their study to be more dependent upon surface forcing than to inversion strength when winds were northerly. Another finding in their study was that the air-sea temperature gradients in post-coldfrontal situations likely enhance cloud liquid water path and that the weaker inversions and deeper MBLs that are observed in these environments tend to be more decoupled.
These studies motivate us to investigate the applicability of the

Scientific Objectives
The goal of this study is to improve understanding of the decoupling process and its impact on MBL cloud morphology in summertime post-cold-frontal air masses traversing the Eastern North Atlantic and beyond. Another goal is to determine if the deepening-warming hypothesis can be used to explain MBL cloud transitions in post-cold frontal air masses. One particular focus is the feedback between decoupling and the surface latent heat flux, which is thought to be a driver of the transition process in the deepening-warming schema. In addition, the simulations presented herein do not employ a shallow convective parameterization and are thus intended to reduce uncertainties associated with lower resolution models that require such a parameterization. Convection arises in a grid cell in our simulations primarily as a consequence of the interactions between the sub-grid scale turbulence and surface flux parameterizations, which are driven by the external conditions supplied to the cell. The cloud fraction reported at each height in each grid cell ranges from zero to one. If the cloud fraction is one at a given height, we view the convection as being "explicitly resolved" at that height, but if it is less than one, the sub-grid scale structure of the convection is determined by the turbulence scheme and is not explicitly resolved.  Lilly (1968). This early study focused on cloud-top radiative heat loss as the main contributor to the maintenance of the sub-cloud mixed layer that spans between 500 m to 1 km or more from the ocean surface, rather than the surface shear-generated turbulence. Therefore, the study primarily considered cloud-top instability in the event of entrainment of dry air above cloud top into the cloud layer as the driver of the transition to broken cloud structure. The results of the Lilly (1968) study are applicable to non-precipitating cases consisting of thin liquid water clouds. Such conditions are typically encountered in coastal ocean upwelling regions, for example immediately adjacent to the California or Portuguese coastlines or in the summertime high latitude regions.
When air from the free troposphere penetrates into the cloud layer, its temperature is reduced due to evaporation of cloud water droplets into the dry air.
This mixed air will be negatively buoyant and sink further in the cloud layer if it reaches saturation at a colder temperature than the surrounding air. If this process occurs spontaneously, it represents an instability that will continue until the cloud vanishes. The condition for the stability of the cloud-top layer with respect to penetration of the dry air aloft in Lilly (1968) was chosen such that the overlying air would have equal or higher wet-bulb potential temperature, θ w , than that at the cloud top. By definition, the wet-bulb potential temperature is the temperature an air parcel would have if it were cooled until reaching saturation and adiabatically brought to a reference pressure. Thus, the temperature of the entrained parcel after reaching saturation through the evaporation of surrounding cloud water would not be lower than the latter, prohibiting negative buoyancy that would result in unstable sinking of the parcel as a convective downdraft.
Thus, the entrained parcel would remain in the cloud and become part of it as a result of entrainment.
The results of the moist cloud model by Lilly (1968) showed that in steady state, the cloud layer is maintained by active convection, with the virtual heat flux moving upwards. This convective energy is supplied by the latent heat of evaporation from the sea surface which is later released as it is condensed in the cloud base. Whereas, the sub-cloud layer is driven by downward virtual heat flux, transporting sensible heat into the sea, since the sea surface is virtually colder than the surface air.
The Lilly (1968) 1-D model successfully predicted a positive vertical gradient in θ w at the inversion consistent with observations. One of the main findings of the study was that the general character of the steady-state solution was irrespective of the choice of an entrainment hypothesis within the prescribed bounds. It was conjectured that the existence of a low, radiatively effective cloud cover is key to maintaining such a strong inversion characteristic of the marine layer of coastal California and similar regions. Without this element, a 15 − 20 degree inversion at a height of 500 − 1000 meters seemed impossible to maintain by kinetic energy for any observed combination of subsidence, convection, shearing turbulence. Furthermore, Lilly (1968) speculated that this cloud layer might be a residue of the frontal disturbances that commonly pass through the Gulf of Alaska in summer, or simply generated by evaporative moistening of a previously clear mixed layer. In other words, the original genesis of the cloud layer described in the model was unknown. Randall (1980b) further investigated the criterion for cloud-top instability with respect to entrainment of the cloud-top air that could lead to buoyant convection in the cloud layer, the so-called "Cloud-Top Entrainment Instability" (CTEI).
The underlying motivation of this investigation was an attempt to uncover the processes that lead to the transition of a relatively steady-state stratocumulus deck into broken cloud elements. This CTEI study focused on the entrainment process itself, regardless of the cloud-top radiative cooling that is known to stimulate convection and turbulence in the cloud layer (Lilly (1968); Randall (1980a)).
The CTEI criterion was linked to whether entrainment leads to the destruction or generation of turbulence kinetic energy in a stratocumulus cloud. In a statically stable layer, the buoyancy force tends to attenuate or almost completely prevent the penetration of the dry air from aloft into the cloud layer. The work against this force is supplied by Turbulence Kinetic Energy (TKE). Randall (1980b) considered a uniform turbulent and saturated cloud layer underlying a quiet and unsaturated layer. Virtual dry static energy, s v , was employed in the study to account for the effect of vapor and liquid water in the cloud parcels on their relative buoyancy, defined as where T v is the virtual temperature, which is In this expression, q and l are water vapor and liquid water mixing ratios, respectively, c p is the specific heat at constant pressure, and δ = 0.608. As conservative variables, the total water mixing ratio (q + l) and the moist static energy h ≡ c p T + gz + Lq were used in this study, where L is the latent heat of vaporization.
Entrainment of cloud-top air into the cloud layer would promote convection, or equivalently increase turbulent fluxes of s v , F sv , at the highest level within the cloud under the condition that From Eq. 2.3 and Eq. 2.4, the criterion for the onset of this type of cloud-top instability, termed by Randall (1980b) as "Conditional Instability of the First Kind Upside Down", or CIFKU, can be written as According to Eq. 2.6, the requirement for the onset of this type of cloud-top instability is fulfilled when ∆s v is smaller than a critical value, (∆s v ) crit , which has an inverse relationship with the relative humidity of the air atop the cloud layer. This imposes a more stringent requirement than that of Lilly (1968), which considers decrease in the wet bulb potential temperature with height at cloud-top as the criterion for cloud-top instability, which is equivalent to ∆h < 0, while with CIFKU ∆h needs to be smaller than a negative value; ∆h < L∆(q + l)/β < 0.   Randall (1980a) In a separate study, Deardorff (1980) found similar results to Randall (1980b) for the cloud-top instability criterion with respect to entrainment of dry air into the cloud top.
The vertical structure of the turbulence in the marine boundary layer and its potential impacts was originally addressed by Turton and Nicholls (1987) and Nicholls and Turton (1986) in a series of studies. Nicholls and Turton (1986)  moistening it through detrainment, and either immediately resupply the stratocumulus beneath the inversion with liquid water and aerosols, or slowly attenuate entrainment in the desiccated layer until subsequent convective elements were able to resupply the layer of stratocumulus. The models of Nicholls and Turton (1986) were partially validated during ASTEX, and the progression of decoupling leading to moisture accumulation, CAPE accumulation, the LCL entering layer near the ocean surface, and the onset of cumulus convection were observed by Miller and Albrecht (1995) and by Miller et al. (1998).
This intermittent decoupling, sometimes referred to as a "cumulus coupled" marine boundary layer, often exhibits a level of mesoscale organization (Miller and Albrecht (1995); and others). Mesoscale clusters of cumulus clouds, some reaching into the overlaying stratocumulus and resupplying it, are often observed in the marine boundary layer in the ENA region and beyond. These mesoscale clusters, sometimes referred to as "mushrooms" or "shrooms" are characterized by an extensive stratocumulus shield being fed by cumulus clouds from below. They produce extensive drizzle storms and occasionally evaporation-driven penetrating downdrafts in their wake. They resemble a miniature version of marine tropical deep convection. The association between "shrooms" and marine fronts observed during ASTEX was never quite clear. At least one such structure was clearly associated with a marine cold front as observed during ASTEX (M.A. Miller, personal communication), but others appeared to exist independent of the frontal zone itself, but often in the post-cold-frontal air masses that passed when winds in the ENA were from the north.
Focus on cloud-top mixing as the primary driver of the marine stratocumulus to cumulus transition faded shortly after ASTEX as the role of cumulus convection in the transition process became more obvious. A series of 1-D modeling efforts led to the supposition that the transition process was driven from below by latent heat release associated with the cumulus convection. These 1-D simulations successfully simulated a transition from solid stratocumulus to cumulus and represented aspects of the intervening transition process in which cumulus and stratocumulus coexist in a "cumulus coupled" environment. The results of these 1-D models led to the introduction of the "deepening-warming" hypothesis, in which it is postulated that that the warming and moistening of the air mass adjacent to the ocean surface in conjunction with weakening subsidence leads to the observed latitudinal cloud transition. A hallmark of the hypothesis is the minimal role of cloud-top entrainment driven by LW cooling in the transition process.
The latitudinal deepening-warming has heretofore been viewed as a climatological process and the mesoscale clusters of cumulus rising into stratocumulus in a more or less local framework. There is no mathematical connection in the 1-D models between the mesoscale organization in the structures observed over the ENA and the deepening-warming hypothesis. Moreover, the association between the hypothesis, the clusters, and the frequently observed frontal structures over the ENA and over the Earth's mid-latitude ocean regions, is largely unknown.
In essence, the tools to make such connections, high-resolution models capable of resolving these clusters and their relationship to the deepening-warming hypothesis and frontal structures, have been heretofore non-existent. Recent advances in computing power now enable such simulations and they are used in this dissertation as a first attempt to address the interplay between these major components of the cloud transition region over the ENA.

Marine Frontal Clouds
A prominent feature of cyclonic storms in the mid-latitudes is the resulting cloud systems that are characteristic of these regions. These cloud structures greatly alter vertical fluxes of latent heat, momentum, and moisture, and comprise a major part of the precipitation and cloud radiative forcing in the mid-latitudes, all of which have substantial climatic implications (Stewart et al., 1998). These cyclones form as a result of baroclinic instabilities present in these regions. During their life cycle, strong oceanic convective systems evolve into stratiform structures (e.g. ; Hoskins (1982); ).
Fronts are transition zones between air masses with contrasting thermal properties. Frontal zones are often characterized by significant weather, including sudden changes in cloud structures and precipitation, as well as drastic gradients of temperature, moisture, wind direction, and horizontal and vertical wind shear (Carlson, 1991). The significant role of fronts in cyclogenesis has been been well established (Bjerknes (1919); Hoskins (1982); Bluestein (1986)). The clouds associated with mid-latitude cyclones are organized by the large-scale atmospheric dynamics, and cloud structures are different in each of the three (or more) air masses that are associated with the fronts. Fronts are generally synoptic-scale features, with mesoscale and cloud-scale motions taking place within them that act to modify moisture and temperature profiles, and consequently, are important features that require quantitative analysis (Stewart et al., 1998). Therefore, it is crucial to properly represent these heat and momentum transfers associated with the fronts and the corresponding air masses in climate models in order to adequately predict larger-scale features, and their impact on climate change.
However, a major source of uncertainty in climate model for predicting future climate can be traced to model representations of low-altitude clouds (Bony and Dufresne (2005); Medeiros et al. (2008); Williams and Webb (2009)) and marine stratocumulus clouds comprise a substantial portion of these cloud systems. Marine stratocumulus cloud systems are ubiquitous beneath anticyclones that occur under the subsiding branch of the Hadley circulation, as well as in post-coldfrontal regions of mid-latitude cyclones (Field and Wood (2007)). Trenberth and Fasullo (2010) discovered that most general circulation models suffer from negative cloud-cover bias that leads to higher rates of shortwave radiation absorption over the southern oceans. This bias is mostly associated with cold sectors of extratropical cyclone regions, and has been found to exist in the phase 5 of the Chapter 3

Overview of Data and Methodology
Two principal resources were utilized in this work: data collected using the suite of sensors deployed at the ENA site and the WRF model. These tools are described in the sections that follow.

Observation Data from ENA
Observations used in this study were collected at the Atmospheric Radiation When the DL pulse reaches the boundary of a cloud, its pulse is scattered such that it operates as a cloud-base height detector. Extremely thin clouds may allow the laser pulse to transmit through without being decimated thereby enabling additional cloud boundary detection above. Two minute averages of instantaneous vertical velocity measured by the DL that are centered on the half-hour output times of the WRF model are used in this study.
A K-a band, Zenith-pointing Radar (KAZR) operating at λ = 8.0 mm detects hydrometeors to a height of 18 km with an initial measurement height of 72 m.
A profile of cloud location and hydrometeor Doppler vertical velocity is collected every 2-seconds in all but heavy precipitation when its beam is attenuated by liquid water absorption. When a cloud is drizzling and the droplet spectrum is bi-modal, the Doppler velocity becomes ambiguous unless spectral processing techniques are applied. These techniques may, at times, introduce considerable uncertainty in the estimation of the Doppler velocity of cloud droplets, so incloud KAZR measurements of the in-cloud vertical velocity are not used in this study. The KAZR data are used to compute the average, height-dependent cloud fraction over 30-min averaging interval.

Parameterizations Used in WRF
The focus of the WRF simulations in this work was to better understand the life cycle of marine stratocumulus in the post-cold-frontal environment. To achieve the best possible representation of the observed conditions, various combinations of parameterizations and dynamic options were tested. Test simulations were conduced to the extent that the computational resources available at the time allowed. We found the combination of the following options that rendered the closest results to the satellite images and observation data from the site for the aforementioned periods. Table. 3.1 shows a summary of the parameterization schemes used in this study, as well as their corresponding WRF options. More details on the complete set of options can be found in Appendix A. The following sections discuss the choice of parameterizations outlined in Table. 3.1 in greater detail.
It must be noted that, since the duration of the simulations were not long enough (on the order of months or years), time-varying data such as monthly albedo values, sea-surface temperature (SST), and vegetation fraction were unnecessary and these fields remained constant and equal to their input values throughout the simulation period.

Cloud Microphysics Parameterization
An important factor determining the accuracy of atmospheric climate model simulations is representing clouds as correctly as possible. A critical aspect that controls these representations is cloud microphysics and precipitation processes.
However, accurate parameterization of these processes in climate models is still a work in progress, mainly because the cloud microphysical processes governing organization of hydrometeors are rather complex (Morrison and Milbrandt, 2015). There are two general approaches to representing cloud particles and their microphysical properties in climate models: the bulk method and the explicit bin-resolving method. Even though the latter approach leads to more accurate solutions to drop size distribution (DSD) than the former, their relatively high Land Surface NOAH sf surface physics 2 computational cost in climate models is a prohibitive factor. Therefore, the simpler bulk microphysical parameterizations that predict several DSD moments have been extensively implemented in mesoscale as well as some GCM models (Lim and Hong, 2010).
There are two main categories of bulk microphysics scheme: single-moment and multiple-moment. A single-moment scheme uses a distribution function, such as an exponential or gamma function, to represent hydrometeor size distributions for different categories of particle types, and then predicts the mixing ratios of these hydrometeors (Kessler (1969) In this study, test simulations were run with the following double-moment bulk microphysical parameterizations in WRF: Morrison double-moment scheme (Morrison et al., 2009), WRF Double-Moment 6-class (WDM6) scheme (Lim and Hong, 2010), and Thompson aerosol-aware scheme Thompson and Eidhammer (2014). The Morrison scheme predicts the mixing ratios and number concentrations of cloud droplets, cloud rain, graupel, ice, and snow. This scheme is an improvement on the parameterization proposed by Morrison and Pinto (2005), since prognostic variables for the mixing ratio and number concentration of graupel are also added to the model. WDM6 scheme supplemented double-moment predictions of warm-phase mixing ratios and number concentrations of cloud and rainwater to WRF single-moment 6-class (WSM6) scheme (Hong et al. (2004); Hong and Lim (2006)), while adopting the same ice-phase microphysics of Hong et al. (2004). An important feature of this scheme is that it incorporates a prognostic variable for CCN that helps to activate cloud water and ice particles. The activation process in WDM6 is governed by the interaction between the number of activated CCN and supersaturation (Twomey (1959); Khairoutdinov and Kogan (2000)). Thus, this scheme handles the process of the growth of cloud water by first computing the cloud drop nucleation process before proceeding to the condensation phase. However, WDM6 ignores all other CCN sink-source terms except the CCN activation and cloud droplet evaporation, and the CCN activation process adopted follows a simplified form of the Kohler equation which requires accurate estimation of the environmental supersaturation value, hence the CCN activation process is not adequately represented (Lim and Hong, 2010).
The importance of cloud-aerosol interaction in cloud microphysics and radiative properties modeling is demonstrated by Ramanathan et al. (2001), Wang (2005, and Khain et al. (2008). It is crucial to accurately represent the role of CCNs (known as aerosol effect) in cloud microphysics parameterizations, as it not only directly dictates the rate of generation of cloud water and ice particles, but has many indirect effects, as well. Among the most prominent studies on indirect aerosol effect are the first and second indirect effects proposed by Twomey et al. (1974) and Albrecht (1989), respectively. The first aerosol effect refers to the increased cloud-top albedo due to generation of larger numbers of smaller droplets for a given cloud liquid water content as a result of higher aerosol concentration.
The smaller sized cloud droplets, on the other hand, lead to reduced precipitation, causing the cloud to have a longer lifetime. This is known as the second indirect effect. However, it is a complex task to accurately quantify the aerosol effect on microphysical processes that, in turn, impact cloud radiation, precipitation, and dynamics, due to various interactions between aerosol and cloud ice particles, as well as many feedback processes involving other aspects of cloud dynamics (Levin and Cotton, 2008).
The improved Thompson aerosol-aware bulk microphysics scheme (Thompson and Eidhammer, 2014), that is based on the parameterization proposed by Thompson et al. (2008), incorporates the cloud and ice droplet nucleation processes by CCNs and ice nuclei (IN) aerosols, respectively, and explicitly predicts the number concentrations of cloud water droplets and the available CCN and IN. The Thompson et al. (2008) scheme includes treatments of microphysical processes for cloud water, cloud ice, rain, snow, and a hybrid graupel-hail category. Implementations of this scheme in WRF demonstrated good performance in terms of forecasting precipitation when compared to precipitation measurements (Liu et al., 2011). Thompson  (2012)). Thompson aerosol-aware scheme was also designed to be well-suited for use in high-resolution simulations. Thompson and Eidhammer (2014) conducted a series of experiments using high-resolution (4-km grid spacing) 72-hour WRF simulations of a winter storm event over the entire contiguous United States to evaluate the performance of this new aerosol-aware scheme and the effect of changing hygroscopic aerosol concentrations on the development processes of clouds and precipitation compared to the non-aerosol scheme of Thompson et al. (2008). The results showed increased number concentration of cloud droplets with significant drop in their mean size, heightened cloud albedo consistent with the first aerosol indirect effect, as well as reduced warm-rain rates that were expected on the basis of the second aerosol indirect effect, following increased aerosol concentration.
Using this last microphysical parameterization in our simulations yielded the best results, specifically in terms of maintenance of non-precipitating low marine stratocumulus decks. We think this is mostly due to the more comprehensive treatment of the aerosol effects on cloud water and ice nucleation processes, cloud droplet size and concentration distribution, and precipitation processes.
Therefore, the Thompson aerosol aware was the microphysical scheme that we ultimately opted to use in our simulations.

Planetary Boundary Layer (PBL) Parameterization
A crucial aspect of mesoscale models that significantly impacts forecast inaccuracy is the representation of lower tropospheric turbulent mixing that dictates the vertical thermodynamic and kinematic atmospheric profiles in this region.  2010)). The lower tropospheric air is in close contact with and directly affected by the earth's surface through heat, moisture and momentum exchanges with the surface. Turbulent eddies in this region lead to a mixing process that accommodates these exchanges. Since the spatio-temporal scales of these eddies are smaller than the grid scales and temporal resolution in most mesoscale models, they cannot be explicitly resolved. That is why mesoscale models need to employ parameterizations to represent the subgrid-scale lower tropospheric turbulent motions. Fig. 3.2 depicts the temporal evolution of the PBL. The surface layer is the closest to the earth's surface. It is in this region that the vertical shear is the strongest and molecular diffusion can be considered in the momentum equation (Holton, 1973). During the daytime, convection fueled by the heat fluxes from the surface lead to a convective mixed layer. The resulting convective plumes can generate cumulus cloud layers at the intersection with the entrainment zone which lies above the capping inversion. As daytime comes to a close and solar heating recedes, the surface starts to cool, and a stable nocturnal boundary layer forms in the lower PBL, with the remnants of the earlier convective layer residing atop this layer. In other words, the upper part of the PBL becomes decoupled from the surface.
PBL parameterization schemes estimate turbulence by decomposing equations of motion into mean and perturbation components, where the former represents the time-averaged variables that define the background atmospheric state, while the latter captures the turbulent fluctuations, or deviations, from the mean state.
The problem that is encountered with respect to the perturbations is that there are always more unknown than known terms, and the order of the unknown terms surpasses all the other terms. Therefore, in order to be able to solve the set of equations, n th -order turbulence closure is needed to provide the relationship between unknown terms of moment n+1 and known terms of lower moments. One major area where various PBL schemes differ is the order of turbulence closure.
The other significant distinction among them is whether they employ a local or non-local mixing approach. Local closure schemes only consider vertical levels adjacent to a given point to determine variables at that point, while this is not a necessary requirement in the case of non-local closure schemes.
Since the largest eddies in the PBL are responsible for the vertical mixing process throughout this region for the most part, and are mainly untouched by the local fluctuations in static stability, non-local schemes more accurately represent the mixing processes in the PBL (Stensrud, 2009). During daytime, heat from the surface is transported upward by these large eddies regardless of the localized stability maxima, and it is these eddies that can find their way to the top of the mixed layer, thereby entraining the air from the free atmosphere into the mixed layer, and deepening the PBL (Stull, 1988). However, using higher orders of turbulence closure (Mellor and Yamada (1982), Nakanishi and Niino (2009)) at a higher computational cost can improve local schemes.
Mellor (1973) followed Rotta's energy redistribution hypothesis (Rotta, 1951) and Kolmogoroff's hypothesis of local, small-scale isotropy (Kolmogorov, 1941) to approximate higher moment terms that appear in the set of prognostic equations. Rotta's model simplifies the energy redistribution term that appears in the Reynolds stress equation by suggesting that it may be proportional to Reynolds stress and mean wind shear. Following the assumption that the constitutive coefficients are isotropic tensors The same logic is applied to the equivalent term in the heat conduction moments equation where q ≡ (u 2 i ) 1 2 , and the lengths l 1 and l 2 , and the constant C are empirically determined. Mellor and Yamada (1974) also used Kolmogoroff's model to estimate the dissipation terms, as follow However, since an isotropic first-order tensor is non-existent, Other simplifications involve ignoring pressure diffusional terms pu i and pθ following Hanjalić and Launder (1972), and scaling 3rd-order moment turbulent velocity diffusion terms to 2nd-order gradients (Eqs. 3.9, 3.10, 3.11), where λ 1 , λ 2 and λ 3 are length parameters.
Mellor and Yamada (1974) offered a hierarchy of four PBL turbulence closure models based on the turbulent moment equations derived by Mellor (1973). The different levels of closure models are obtained by ordering terms as products of the degree of anisotropy, a, with a = 0 denoting isotropy, and q 3 /Λ. a is assumed to be small, and have the same order as the non-dimensional measure of the departure from isotropy, a ij The Level 4 model retains all the terms in the Reynolds stress model equations presented by Mellor (1973) after implementing the aforementioned simplifications to Eqs. 3.1, 3.2, and 3.3, except for those including the Coriolis parameter, since their values are comparatively small in the boundary layer. The only further assumption was that all length scales are proportional to the boundary layer length scale, l, as such: where A 1 , A 2 , B 1 , and B 2 were empirically determined, and l is the boundary layer length scale. Following Blackadar (1962) where k is the von Karman constant and z is the height from the surface. As z approaches to zero, l ∼ kz, while as z increases to infinity, i.e. exits the bounds of the PBL, l ∼ l 0 . Mellor and Yamada (1974) adopted the following expression for l 0 that incorporates features of the turbulence field, where α is a constant that is empirically determined, with a suggested value of 0.1 by Mellor and Yamada (1974). Finally, it was assumed that λ 3 = λ 2 = λ 1 , and equal to 0.23l following Mellor and Herring (1973).
Under the assumption that diffusion and advection terms are O(a), more precisely O(a q 3 Λ ), neglecting O(a 2 ) terms such as time-rate of change, advection and diffusion terms for the anisotropic components of turbulence moments yields the Level 3 model. This closure model retains prognostic equations only for TKE and θ 2 . Mellor and Yamada (1982) introduced a modification of Level 3 closure model, referred to as level 2.5, that serves as an intermediate step between Level 3 and Level 2 models. In this model, the material derivative and diffusion terms in the prognostic equation for θ 2 are ignored, as opposed to neglecting advection and diffusion terms from the set of prognostic equations (Level 2) (Mellor and Yamada, 1974). This model still benefits from retaining the isotropic components of transient and diffusive turbulent processes as Level 3 model, while reducing its computational cost by eliminating one differential equation.

Mellor-Yamada Nakanishi and Niino Level 3 (MYNN3)
Nakanishi and Niino (2004) (Sun and Ogura, 1980), as well as the decaying magnitude of TKE in a stably stratified nocturnal PBL (Turton and Brown, 1987), based on comparisons with large-eddy simulation (LES) data. This new model performed better than the improved version of the MY 2.5 Level model, although the latter is computationally less expensive.
Despite superior performance of the Level 3 Nakanishi and Niino (2004) model compared to the MY Level 3 model, it suffered from computational instability during the simulation. This problem was later addressed by Nakanishi and Niino (2006) by ensuring that their model could be realized under the full breadth of conditions encountered in the PBL. They proposed a scheme that imposed a set of constraints on the turbulent terms such as the temperature and velocity variances and the length scale in the case of turbulence. This model is heretofore referred to as MYNN. In the model proposed by Nakanishi and Niino (2004), S M and S H are stability functions for momentum, and heat and moisture, respectively, and components of the eddy diffusivity coefficients of the turbulent flux terms in the Level 3 model. In other words, they added a stability correction to regulate the turbulence so as to prevent runaway turbulence under conditions of rapid PBL growth. These stability correction terms are written as the sum of their expressions in the 2.5 Level model and the difference between Level 3 and Level 2.5 models, denoted by the subscript 2.5 and a prime, respectively (Eqs. 3.18).
The denominators of the stability functions, D and D are functions of the stability ratio G H , defined as where L is the master length scale, q 2 /2 TKE per unit mass, Θ 0 the reference potential temperature, Θ l the liquid water potential temperature, Q w the total water content, and β θ and β q are functions of water vapor fluctuations determined from the condensation process (Nakanishi and Niino, 2004).
In this scenario, D and D , may change signs with changes in G H due to variations in stability in the flow, leading to a singularity problem. For unstable stratified conditions, constraints proposed by Helfand and Labraga (1988) prevent the singularity problem. In stable stratified boundary layers, D can become nonpositive. To avoid this, Nakanishi and Niino (2006) found the following restriction on L/q, Since the vertical gradient of virtual potential temperature, ∂Θ v /∂z can be written as the sum of the terms in parentheses in Eq. 3.20, the condition on L prescribed by Eq. 3.20 reduces to L ≤ q/N , where N is the Brunt-Väisälä frequency. This condition constrains the vertical distance that a parcel of air can move upward against the buoyancy force in a stable layer.
Perhaps, the main improvement of MYNN over MY is the specification of L.
Unlike Mellor (1973) that disregarded effects of stability on L, MYNN proposed a new specification for L, such that it is constrained by the smallest of L S , L T , where L S is length scale in the surface layer that is proportional to kz, L T is similar to Eq. 3.16, except that α = 0.23, which depends on the length of the turbulent field in the PBL, and L B is the buoyancy length scale and proportional to q/N . Nakanishi and Niino (2009) son (1989) and Gayno (1994). They all employed moist-conserved variables to handle changes in boundary layer stratification due to saturation. However, GB01 succeeded at producing sharp inversions, in contrast to the diffused inversions above SCBLs rendered by other models, using an explicit entrainment closure model for determining entrainment diffusivity at the boundary of convective layers.
The moist-conserved variables used in UW scheme are total specific humidity where q v , q l , and q i are the specific humidity of water vapor, liquid water, and ice, respectively, and L and L f are the latent heats of vaporization and freezing. Moreover, since conditions are well-mixed within a turbulent layers, profiles of momentum and scalar variables remain rather unaltered by non-local scalar transport. Hence, diffusivities of these variables in this region was based on local TKE (e) estimation (Bretherton and Park, 2009). Turbulent fluxes of horizontal velocity (u and v) and moist-conserved variables s l and q t are represented as downgradient diffusion with eddy diffusivity and viscosity defined as where l is the turbulent master length scale characterizing TKE dissipation, where Ri CL is the bulk Richardson number for the convective turbulent layer.
Based on this expression, η is much larger in turbulent layers that are strongly convective, where Ri CL takes on negative values.
The main difference between UW parameterization scheme (Bretherton and Park, 2009), and MY and MYNN, is that in the UW scheme TKE is diagnostic rather than prognostic. UW neglects TKE storage, and models the production, transport and dissipation components of the horizontal mean TKE equation. The buoyancy and shear TKE production terms, B and P s , are prescribed as is the vertical shear squared, b buoyancy perturbation, and N 2 the buoyancy frequency, which is formulated in terms of vertical gradients of s l and q t , with their coefficients depending on cloud fraction. TKE dissipation rate, D, is parameterized by Bretherton and Park (2009) as where b 1 = 5.8. The TKE turbulent transport term, T e , the convergence of the sum of vertical TKE flux and work done by pressure, is only considered in convective layers in UW, and is modeled in terms of turbulent-layer-mean TKE < e > (Bretherton and Park, 2009): where a e is TKE relaxation rate. Bretherton and Park (2009) empirically determined a e = 1 for a convective layer from LES data.
An important feature of the UW model is that when buoyancy fluxes become negative in the turbulent layer, it decouples and the conditionally unstable convective layer with moist saturated updrafts and unsaturated dry downdrafts that ensue are handled by cumulus parameterization. In a stably stratified turbulent layer, the UW model ignores turbulent TKE transport and storage, and the closure model simplifies to a first-order one. UW follows the Nicholls and Turton (1986) Niino, 2006, 2009). These tests indicated that MYNN3 produced the closest qualitative representation of the clouds as compared to satellite observations and data collected at ACE-ENA and was the PBL scheme that was ultimately used. This PBL scheme is based on the Mellor-Yamada model (MY) Yamada (1974, 1982), both of which are prognostic turbulent kinetic energy (TKE) schemes. To review, TKE quantifies the magnitude of turbulence in the PBL, which is comprised of contributions from vertical wind shear, buoyancy, turbulent transport, and dampening due to molecular viscosity. An alternative PBL scheme that uses MY as its base and predicts TKE values is the Mellor-Yamada-Janjic (MYJ) Janjić (1994). It is a 1.5-order closure scheme that is an improvement upon MY Level 2.5 turbulent closure model. A 1.5-order closure parameterization scheme incorporates diagnostic second-order moments for variables such as variances of potential temperature and mixing ratios, and first-order moments for other variables to predict secondorder TKE (Cohen et al., 2015). As noted below, this alternative was dismissed on the basis of a past study. Huang et al. (2013) used WRF in single column mode to conduct three case studies over subtropical oceanic regions to assess the performance of various PBL schemes. Included in their study were MYNN, MYJ, and YSU (Yonsei University PBL ) parameterizations, as well as the total energy-mass flux (TEMF) (Angevine et al., 2010), and the results were evaluated according to their ability to represent stratocumulus and shallow cumulus clouds when compared to LES data from those regions. Their study concluded that in all three cases, MYNN better simulated vertical profiles of liquid water and potential temperature, and provided the least-biased thermodynamic structures than MYJ.
In our study, we also utilized the following optional settings in MYNN3: TKE advection, stochastic eddy diffusivity mass flux (StEM) Sušelj et al. (2013), momentum and TKE transport in mass flux scheme, a modified cloud-scale mixing length Ito et al. (2015), and cloud-top radiational cooling contributing to TKE production.  These structures are altogether absent in the UW plot in Fig. 3.3.
As discussed in the preceding subsections, the MYNN stability-adjusted specification of the turbulent master length scale, l, incorporates a buoyancy length scale, L B , while UW follows the approach taken by Mellor and Yamada (1974) to

Radiation Parameterization
It is crucial to correctly estimate the atmospheric radiative forcing in climate models in order to ensure their prediction accuracy. Among the various radiation transfer models available in WRF, the Rapid Radiative Transfer Model for GCMs (RRTMG) (Iacono et al., 2008), which is a modified version of RRTM (Mlawer et al., 1997) that can be readily applied to GCMs, was used in the longwave and shortwave radiation calculations in the model simulations in this study. RRTM and RRTMG are based on the line-by-line radiative transfer model (LBLRTM), which is an efficient and accurate model that applies to the full spectral range (Clough et al. (1992); Clough and Iacono (1995)). LBLRTM accounts for all significant molecular absorbers, and includes the continua of carbon dioxide, oxygen, nitrogen, ozone, and Rayleigh scattering extinction.
RRTM was developed for both the longwave (LW) and shortwave (SW) regions, and designed with the goal of increasing the accuracy of radiative fluxes and cooling rates (Mlawer et al., 1997). In order to accurately account for radiative transfer due to scattering in the atmosphere, a radiative transfer model  2003)), which is a statistical technique, has been adopted by RRTMG LW and RRTMG SW.
The advantage of using RRTMG as the radiation transfer scheme in our simulations, aside from its accuracy, was that we took advantage of the option provided by WRF to input climatological water and ice-friendly aerosols data (aer opt = 3 ), which is compatible with Thompson microphysics. This option is particularly useful because the quality of representing radiative processes associated with water vapor, ozone, and long-lived greenhouse gases, directly impacts the accuracy of cloud simulations in the model.

Land Surface and Surface Layer Parameterizations
The Fortunately, in the case of this study, the complexities of land surface variability and hydrology processes are substantially mitigated, due to the fact that the majority of the simulation domain is located over the ocean. The large heat capacity of the ocean prevents drastic diurnal temperature gradients as those characteristic of land surfaces. Also, the surface properties of the ocean is considerably more stable than those on land; multiple factors affecting land surface features, such as snow cover, soil texture, vegetation, orography, seasonal surface emissivity, and urban land-use, are absent in the case of the ocean.
The land surfaces on the islands within the simulation domains in this study were parameterized using the unified Noah land surface model (Tewari et al., 2004). This scheme has four soil layers (10, 30, 60, and 100 cm thick), and solves for prognostic land states, such as surface skin temperature, total and liquid soil moisture, and soil temperature at each layer, canopy water content, and physical snow depth. The initial conditions for these prognostic equations are provided by the preprocessed input data generated by WRF Preprocessing System (WPS).
The coupling between the surface skin layer and the first atmospheric model level is parameterized by a surface layer scheme. This scheme provides exchange coefficients for sensible and latent heat and moisture to the land surface scheme, and the friction stress and water/surface fluxes of heat and moisture directly to the PBL scheme. These fluxes are particularly sensitive to the treatment of the thermal and moisture roughness lengths, that are estimated by the surface layer scheme. Since MYNN3 was used as the PBL scheme, the MYNN surface layer was required for simulations in this study.

Cumulus Parameterization
Processes associated with cumulus clouds substantially impact the large-scale temperature and moisture fields. Cumulus-induced subsidence leads to largescale warming and drying of the air below, while cumulus detrainment results in large-scale cooling and moistening (Arakawa and Schubert, 1974). In the atmosphere, shallow and deep convective clouds usually coexist with each other.
Shallow convection operates on lower altitudes than those where precipitation processes are activated, resulting in non-precipitating, low-level cumulus clouds.
These convective updrafts strongly modify the depth of the boundary layer, temperature, moisture, wind, and cloud structures by taking the surface mixed-layer air to the lower free troposphere ; Kain (2004);Wang et al. (2004a); Wang et al. (2004b)).
The role of a cumulus parameterization scheme in climate models is to estimate the rate of convective precipitation, latent heat release, as well as the convective vertical redistribution of heat, moisture, and momentum, in the subgrid scale (Kain and Fritsch, 1990). However, because the resolution of the parent domain in our simulations was fine enough to explicitly resolve deep convection (< 10km), no cumulus parameterization scheme was required in our case. Thomas and Schultz (2019a) offered the following five components for defining and diagnosing fronts: a) a thermodynamic quantity or quantities, such as (equivalent) potential temperature, and wind, for identifying fronts; b) a function(s) acting on a quantity (e.g. frontogenesis and gradient functions), to create the field used to detecting fronts; c) a level(s) to operate the function(s) on the quantity, including surface layer, 850 hPa, or 1 km above sea level; d) a threshold field value(s) beyond which a frontal structure is recognized; e) identifying the category (e.g.

Identifying Fronts
warm or cold) and the spatial scales of the frontal zone.
In this study, fronts in the vicinity of the Azores and at the ENA site were In association with the frontal passage in the 2018 case, wind speed changed from 9 to about 1 ms −1 . Relative humidity from observed values dropped from about 90% to a little above 50%, while those from simulation decreased from 100% to 80%. Also, surface temperature decreased from 25 • C to 20 • C, and 23 • C to 16 • C for those obtained from measurements and WRF model, respectively.  Potential temperature and horizontal wind fields at 850 hPa during the passage of fronts over the Azores, which lie in the geometric center of the plots, for the two simulations are shown in Fig. 3.7. In both cases, sharp gradients in both variables are evident in the first two plots in the leftmost column. Similar to the time-series analysis of the five surface variables used to detect fronts discussed above, it can be seen in Fig. 3.7 that the frontal passage of July 2018 seems to be more vigorous than that of July 2017, judging from the larger spatial gradient in horizontal wind direction and the potential temperature, as well as southward advancement of the front in 2018. This analysis further bolsters those obtained from surface variables.
Overall, the WRF simulations for the two cases reasonably reproduce the surface observations of the frontal and post-frontal environments at the ENA site, as demonstrated in Fig. 3.6. Island effects, which have been reported in past studies from islands in this region (Miller and Albrecht, 1995), are evident in the observations and in the simulations. The daytime cold bias in the temperature simulations, which also leads to increases in the surface relative humidity, is likely attributable to sampling biases because a portion of the model's grid cell lies over the ocean. Given that WRF model simulations did a reasonable job of reproducing surface changes associated with the passage of these weak cold fronts, we are relatively confident that the simulated frontal and post-frontal environments over the entirety of the ENA domain bear some resemblance to reality.        Almost half a day after the onset of the breakup of the stratocumulus deck, the column transitioned into widely scattered, large mesoscale cloud clusters (4.3-4.4, rightmost column). These clusters contained relatively strong deep cumulus updrafts rising to nearly 1.5 km, and laterally detrained stratocumulus. These elements extend at least 300 m deeper than the broken clouds earlier along the trajectory (center column). A primary source of transporting moisture from the PBL to the upper troposphere, particularly in the Tropics and mid-latitudes in the summer, are these cumulus convective clouds (Xu and Randall, 1996). The anvil clouds produced by the lateral detrainment of cumulus clouds also comprise the major component of stratiform precipitation associated with mesoscale convective systems (Houze Jr, 1977).   Classifying the clouds in these regions as residing in a coupled or decoupled environment required an objective criteria. It must be noted, however, that we adhere to the notion that decoupling is a continuum and there is no rigid threshold for clouds to instantaneously change from a coupled state to a decoupled one.
What we seek is a measure that can separate the clouds that are either still coupled to the sub-cloud layer or in nascent stages of decoupling from those that are further into the decoupling process, as efficiently as possible. We label them as "coupled" and "uncoupled" clouds, respectively, for want of better terms that can convey our definitions of these two groups explained in the previous statement.
There are various decoupling diagnostic measures proposed in the literature. The decoupling metrics explored in this study are discussed in the following sections.

Minimal Decoupling Criteria
BW97 proposed a "minimal decoupling criterion" for the mixed-layer model of a cloud-topped boundary layer (CTBL) driven radiatively atop areas with fixed SSTs, while ignoring drizzle. This criterion for negative sub-cloud buoyancy fluxes is given by where ∆F R is the net radiative flux divergence across the stratocumulus cloudtopped boundary layer. In their model, BW97 assumed that ∆F R is the sole diabatic forcing in the CTBL and that its value is fixed. This flux divergence is given by F + R −F R (0), where F + R is the net radiative flux atop the cloud, and F R (0), the net upward radiative flux at surface. In Eq. 5.1, F L (0) is latent heat flux from surface, A, a non-dimensional entrainment efficiency, δz c , the cloud thickness, and z i , the inversion height. Lastly, η ≈ 0.9 is a thermodynamic coefficient that is weakly dependent on pressure and temperature.
This criterion was defined for typical values in subtropical stratocumulus and pertains to the moment when the sub-cloud buoyancy flux at the cloud base becomes negative. As acknowledged in BW97, this is the "broadest possible criterion for decoupling in that any negative buoyancy fluxes below cloud base indicate decoupling." According to Eq. 5.1, decoupling initiates when the ratio of radiative flux divergence in the CTBL to the surface latent heat flux (∆F R /F L (0)) becomes smaller than a threshold value. This threshold is directly proportional to δz c /z i , that reflects the fraction of the boundary layer occupied by the cloud layer, and the parameter representing the efficiency of entrainment, A. The threshold value of 0.45 in Eq. 5.1 is obtained for the typical case in which the cloud layer fills approximately 25% of the depth of the boundary layer, and A = 2.
Adopting this diagnostic proved particularly challenging in this study, due to the fact that radiative fluxes at cloud top are not provided by WRF, therefore the net cloud-top radiative flux, F + R , and consequently, the radiative flux divergence across the cloud layer cannot be readily computed from the WRF simulations.
The only approach for estimating the radiative flux divergence at cloud top would be to use the cloud-top radiative fluxes with those at TOA (that are available from WRF outputs) as a surrogate for this variable. To do this, it needs to be assumed that the longwave (LW) emission and shortwave (SW) albedo at TOA are very close to their magnitudes at cloud-top, and that there are no high-altitude clouds above those in the vicinity of the boundary layer that are the subjects of this study in the domain of interest. A closer inspection revealed that the approximated flux ratio would not respond strongly to the transition from coupled to decoupled clouds, which may have partly been due to the implementation of these coarse approximations. Jones et al. (2011) proposed three separate decoupling diagnostics that were investigated in this study. The first decoupling diagnostic proposed by Jones et al.

Jones Decoupling Criteria
(2011) takes the vertical structure of the boundary layer into account in the context of changes in two moist-conserved variables, θ and q T , between the surface and inversion, where θ is the liquid potential temperature and q T , the total water mixing ratio. These definitions are where θ is the potential temperature, L the latent heat of vaporization, c p the dry air specific heat at constant pressure, and q and q v are liquid water and water vapor mixing ratios, respectively. When the MBL is well-mixed, θ and q T are approximately constant beneath the capping inversion. Conversely, vertical gradients in θ and q T are observed when the MBL is decoupled. This decoupling diagnostic quantifies the vertical gradients of moisture and temperature, and is given by where variables with subscripts "bot" and "top" are calculated as the mean over the lower and upper 25% of the boundary layer below the inversion. Note that the subscript T has been omitted in Eq. 5.4 and hereafter in this section. Jones et al. (2011) found profiles where ∆q > 0.5 gkg −1 and ∆θ l > 0.5K to be decoupled, based on data collected during the VOCALS Regional Experiment (VOCALS-REx) in October-November 2008 in the MBL over the Southeast Pacific.
These metrics have two advantages in comparison to the one discussed earlier in this section; they capture the definition of decoupling in MBL clouds, hence the decoupling mechanism, in their formulation, and both can be easily and directly calculated from the output of WRF simulations. Therefore, they appeared promising metrics to adopt for the purposes of this study.
The second decoupling diagnostic suggested by Jones et al. (2011) is a subcloud decoupling measure prescribed by the difference between LCL and cloudbase height, z b , which is expressed a The final decoupling diagnostic proposed by Jones et al. (2011) is referred to as the "mixed layer cloud thickness" and is the difference between the inversion height (z i ) and LCL, or When the MBL is well-mixed, ∆z M is equivalent to the cloud thickness. Latent heat released due to condensation within moist updrafts generates a buoyancy flux within the cloud layer. Thus, a thicker cloud layer can have larger values of vertically integrated buoyancy flux along its depth, and consequently, more in-cloud turbulence. A more turbulent cloud layer generates more entrainment, which in turn facilitates mixing of dry and warm air above the inversion with the in-cloud air mass and helps the cloud layer become decoupled from the sub-cloud layer. Therefore, larger values of ∆z M should correspond to more decoupled cloud layers.
The procedure used in Jones et al. (2011) was adopted to determine which of the diagnostics was best suited for use with the WRF simulations, so we tested the second and third decoupling criteria against ∆q, which represented the first diagnostic. In order to calculate ∆q, the capping inversion height z i is required.
Because the capping inversion in the simulations is not sharp in most cases, two methods have been adopted to identify and cross-check z i . The first method is from Jones et al. (2011), where z i is determined as the vertical level with the minimum temperature and relative humidity of at least 45%. The second approach identifies the inversion height as the level where d dz θ l and d dz q T are maximum, while the relative humidity is greater than or equal to 45%. The results were almost the same in both cases for both 2017 and 2018 simulations, and wherever they did not coincide, the z i from the first approach was chosen.

U R N A L O F T H E A T M O S P H E R I C S C I E N C E S VOLUME 77
Downloaded from http://journals.ametsoc.org/jas/article-pdf/77/6/201 (b) Figure 5.4: Composite profiles for a) ∆q < 1.6(g/kg) and ∆q > 2.0(g/kg) on 1800 UTC 17 July 2017 (leftmost two plots), and b) ∆q < 1.7(g/kg) and ∆q > 2.1(g/kg) 1130 UTC 29 July 2018 (rightmost two plots). q T is indicated by the navy blue line, and θ by the black line. The green line shows 10 times the cloud fraction for better visualization. The inversion height z i and LCL are also shown in dashed black and red lines, respectively.
boundary layer become less well-mixed (i.e. the cloud element becomes less coupled, until the cloud-base and in-cloud regions form separate stable entities).
Therefore, we suspect that there is no exact threshold value for distinguishing between the decoupled and coupled clouds, and we adhere to this assumption in our attempts to classify clouds in the simulations according to their coupling state.
Hence, the plots designated as "coupled" and "decoupled" clouds in the following sections are not necessarily those that reside in absolutely well-mixed conditions, or are completely decoupled from the sub-cloud region, respectively, but rather encompass examples that fall somewhere in the decoupling range, leaning towards being more coupled or decoupled than otherwise, respectively.
With this explanation, for both simulations, we designated cloud layers with ∆q < 1.6 as coupled and those with ∆q > 2.1 as decoupled, and neglect those that fall somewhere in between these values. While these threshold values have an element of subjectivity, they are consistent with the decoupling continuum that exists during the decoupling process in both simulations.

Buoyant Production of TKE (QB) and Vertical Velocity (W)
Prior to exploring the decoupling diagnostics discussed above, the decoupling measure that was adopted for use in this study was the sub- This is surprisingly at odds with the general understanding of how overcast stratocumulus clouds are maintained, which is mainly by cloud-top processes (BW97).
As the breakup process initiates, QB values below the cloud base are reduced to the point that the positively buoyant layer periodically detaches completely from 2018 simulations. Vertical velocity at cloud base seemed to be a slightly preferable decoupling metric relative to QB because it exhibited a slightly linear relationship with ∆q values smaller than 1.5 g/kg, especially in 2018, but no robust correlation was found between these two variables. Based on the results in Fig. 5.7, points with W b greater than approximately 5 cm/s roughly correspond with ∆q values less than 2 g/kg. In the same interval, there are a considerable number of points with smaller W b values (mainly negative). Therefore, with W b as the decoupling criteria, many examples that would otherwise fall under the category of "coupled" clouds were classified as "decoupled".
The same conclusion can be reached about QB as a decoupling metric based on with ∆q across all ranges. Therefore, it can be concluded that QB would not be an adequate decoupling metric to be used in this study. We assume that a   The Lagrangian derivatives of CLD B and LH following the cloud base are laid out in Eqs. 6.1-6.2, respectively: The first term on the right-hand side of both equations above represents the local time derivative, and the rest of the terms characterize advection of each variable by the wind. The overbars in the right-hand side of Eq. 6.1 indicate that first cloud-base cloud fraction at each model vertical level that fell within the lowest 100 m of the cloud layer was used to estimate the local time derivative, as well as the advection terms for each of the zonal, meridional, and vertical components with the corresponding velocities at those levels (u, v, and w, respectively).

Correlations between lifting condensation level (LCL) and cloud-base cloud fraction
Latent heat flux modulates the LCL and past studies have also reported strong connections between the height of the LCL and the occurrence of coupled regions within a transition cloud field (Miller and Albrecht, 1995).  The results shown in Fig. 6.1 suggest that the changes in the surface LH following the cloud base may be partly responsible for the initiation of the decoupling process. There are at least two mechanisms that could lead to this link. One potential mechanism echoes the deepening-warming hypothesis; the LH release at cloud base becomes large enough that it fuels much larger updrafts within the cloud, which in turn increase mixing at cloud top and initiate downdrafts that penetrate the cloud base and dry the air below, locally decoupling the cloud from the surface. Another possibility is that warming dominates moistening, deepening the boundary layer and consequently raising the LCL. A higher LCL leads to local moisture accumulation, which eventually lowers the LCL and enables cumulus to from at the top of the decoupled layer and rise to the base of the inversion. This process requires no inherent contribution from cloud-top instabilities. In decoupled regions, the cloud is no longer substantially supplied with buoyant TKE and moisture from below and, as such, exhibits limited response to changes in the surface LH and moisture flux. While it is possible that both of the LH-related breakup scenarios suggested above operate depending upon conditions, cloud-top processes are also explored to better understand stratocumulus breakup and changes in the decoupled cloud structures beyond the breakup point.
6.3 Correlations between horizontal wind speed and relative humidity, and surface latent heat flux An increase in LH followed by the initiation of decoupling at cloud base should correspond with reduction in surface relative humidity (RH) due to cloud-top mixing, as described by the deepening-warming hypothesis proposed by BW97.

Concentrating on the Lagrangian correlations between RH and LH and UV10
and LH because of their consistency and integrative nature, we conclude that both variables are contributors to the observed behavior in the Lagrangian changes in LH. Further, we note that their respective relationships with LH are consistent with the expected behavior of these surface variables in the context of BW97.
We should also add that even though cloud-top mixing may perhaps play a role in reducing RH near the ocean surface under less well-mixed MBL conditions, we cannot definitely link these changes exclusively with cloud-top processes. The difference between correlation coefficients for coupled and decoupled clouds in the

CLOUD-TOP BUOYANCY INSTABILITY
Among the most established cloud-top processes in the literature responsible for stratocumulus breakup are Randall's "Conditional Instability of the First Kind Upside-Down (CIFKU)" (Randall, 1980b) and Deardorff's "Cloud Top Entrainment Instability (CTEI)" (Deardorff, 1980). The instability criteria for stratocumulus-topped boundary layer inversions in both of these studies state that the differences in values of virtual dry static energy or equivalent potential temperature, respectively, between the layer immediately atop the cloud above the inversion and the in-cloud layer below inversion must be smaller than a threshold value In this expression q s is the saturation mixing ratio, q is the mixing ratio, q w is the total water specific humidity, and B+ and B refer to above-cloud and cloud-top layers, respectively. where ∆ is the difference of the respective variables across the cloud-topped inversion, q t , c p , and L are the total water mixing ratio, specific heat of air at constant pressure and latent heat of water evaporation, respectively. There are different threshold values for κ proposed in different studies; the instability criteria found by Randall (1980b) and Deardorff (1980) is κ > 0.23, MacVean and Mason (1990) suggest that cloud-top buoyancy reversal is only observed for κ > 0.7, while a number of studies, e.g. Kuo and Schubert (1988) and Moeng (1986) have found no relationship between κ and stratocumlus breakup at cloud-top.
Each of these cloud-top buoyancy instability criteria requires a reference for calculating values corresponding to level B+, which is immediately above the cloud top beyond the inversion. The inversion height is determined using the two methods explained in Section 4. Thus, B is the last cloudy level below the inversion and B+ is the interpolated value of the desired variable 10 m above B. lation with this variable. It can be inferred that, as long as the cloud layer is still coupled to the sub-cloud layer, cloud-top entrainment due to buoyancy instabilities have little effect in diluting the cloud from above. As the decoupling process progresses, however, these instabilities become more effective at eroding the cloud layer from above, and as the instabilities strengthen, the clouds become even further diluted at cloud top.
A shortcoming of the CTEI criterion is that in climate models and non-LES regional models, accurately calculating the cloud-top jumps in θ e and q t requires assumptions and is not straightforward. The mixing processes that affect CTEI generally occur at scales on the order of a few meters, which are unresolved in our simulations and even in some LES models (e.g. Mellado (2017)). Hence, the cloud-top mixing processes associated with CTEI in our simulations arise from the MYNN turbulence parameterization, which is responsible for mixing free tropospheric air into stratocumulus. Our analysis of CTEI therefore represents the relationship between the parameterized cloud-top mixing and the simulated stratocumulus transition. Simulations from WRF in our study have an average vertical resolution of 80 m at marine stratocumulus cloud top. Therefore, ∆θ e and ∆q t need to be interpolated at cloud-top, which results in a rather rough estimation of κ that may affect the correlation coefficients with respect to CLD U .
It is unclear whether this shortcoming significantly impacts our result, but other studies have found discouraging correlations (Stevens (2010), Mellado (2010), Mellado et al. (2009)). Klein and Hartmann (1993) introduced LTS as a measure of lower-tropospheric instability, defined as the difference between potential temperature at 700 hPa and the surface, in an attempt to uncover the environmental variables that regulate low-level stratiform cloud amount. They found a reasonable correlation between LTS and regional mean cloud amount. LTS has been widely adopted by some models in boundary layer cloud parameterizatios (e.g. Rasch and Kristjánsson (1998); Köhler et al. (2011)). However, Wood and Bretherton (2006) pointed out that since LTS was designed and studied for subtropical regions, it might not be appropriate to use as a parameterization for low clouds in the mid-latitudes.

Estimated Inversion Strength
They proposed "estimated inversion strength" (EIS) as a new measure of the inversion strength for use in parameterizations of the low cloud cover fraction (CF), especially in subtropics and mid-latitudes. The EIS improves upon LTS by accounting for the moist-adiabatic potential temperature gradient at 850 hPa. Naud et al. (2016) reinforced the conclusion of Wood and Bretherton (2006) by demonstrating improved correlations between post-cold-frontal CF and EIS relative to correlations with LTS over mid-latitude oceans in both hemispheres. The EIS is an approximation of potential temperature, θ, at 850 hPa level which roughly corresponds to the typical inversion height for mid-latitude stratocumulus-topped boundary layers. Following Eqs.4 and 5 from Wood and Bretherton (2006), EIS is calculated as, where LT S = θ 700 − θ 0 , z 700 is the height of the 700 hPa isobar, Γ 850 m is the moist-adiabatic potential temperature gradient at T = (T 0 +T 700 )/2 and, pressure P = 850 hPa, as a measure of the moist adiabatic lapse rate at 850 hPa, which is given by In this expression, q s is the saturation mixing ratio, R a and R v are the gas constants for dry and water vapor, respectively, and L v is the latent heat of vaporization.
Correlations between EIS and cloud-top fraction (CLD U ) as a measure of effectiveness of the contribution of EIS changes to MSC breakup are shown in

A New Measure of Inversion Strength Designed for use in Large Scale Models
Because the existing measures of inversion strength studied earlier either fall short of explaining the breakup in the WRF simulations (EIS) or require interpolations to compute (the buoyancy reversal criteria, κ), we investigated an alternative measure which can be more easily calculated using data from climate and regional models with coarser vertical resolution than LES models.
Equivalent potential temperature (θ e ) is another moist-conserved variable that incorporates the heating effect of latent heat release on potential temperature.
As an entrained parcel ascends adiabatically across the inversion, a decrease in θ e indicates a reduction in the latent heating due to reduced moisture content.
Therefore, vertical gradients in θ e indicate moisture changes during adiabatic ascent and descent across the inversion. On the other hand, vertical gradients in the saturated equivalent potential temperature (θ es ) indicate temperature changes during adiabatic ascents and descents across the inversion. Typically, θ es above inversion is higher than below it because the air is warmer, while the opposite is true for θ e , since conditions above inversion are drier. The difference between θ es and θ e at the layers atop and below the inversion indicates how much drier and warmer the air above the inversion is relative to that just below the inversion.
Warmer and drier conditions indicated by a greater difference between θ es and θ e atop the inversion lead to reduced instability of the inversion layer to downward motions from above. Therefore, the inversion is strengthened, i.e. there is a greater resistance to entrainment of air above inversion into the layer below. The reverse is true for smaller values of this difference; the inversion becomes more susceptible to downward entrainment of air above that acts to dilute the cloud layer below. Therefore, the profile of these two variables provides a sense of the strength of the inversion in terms of the contrast between both the moisture and temperature of the two air masses. θ es sharply diverge beyond the inversion height, while having the same value in the in-cloud region. The cloud element at point II is in the broken cloud field downstream of the stratocumulus deck, and is being diluted from above as well as from below, judging by the cloud-fraction values. Unlike point I, θ e and θ es start diverging slowly in the cloud region for point II. This suggests that the inversion has become weak enough, allowing the dry air to penetrate through the cloud, thereby destroying it from above.
Since the model is prescribed with discrete vertical levels, the layer above inversion may be at a different altitude depending on the inversion height. Therefore, in order to offset the bias of the different distances between inversion and the layer above it, the vertical gradient of θ es -θ e at the inversion is used as a measure of the inversion strength. This variable is hereafter labelled "INV " INV is calculated for each example represented in Fig. 7.3, and as expected, it is greater in the case of a stronger inversion; 0.125 and 0.025 for the points I and II, corresponding to strong and weak inversions, respectively. decoupled ones. However, a positive trend can be seen between those with CLD U smaller than 0.5 and INV, which can be related to how strongly coupled they are to the surface. In other words, those clouds are likely in the beginning of the decoupling process, and are more susceptible to changes in inversion strength at cloud-top.
In conclusion, INV performs similar to the buoyancy reversal criteria (κ) in terms of quantifying the impact of inversion instabilities on the breakup of decoupled clouds, while being easier to calculate than the latter in non-LES models.
And because it has proved more successful than EIS, it can be viewed as an alternative metric for explaining the relationship between cloud-top instabilities and cloud-top erosion. Also, it can be interpreted that cloud-top instability and inversion strength measures may perhaps be byproducts of the breakup processes at cloud base. The main evidence is that they only seem to be correlated with breakup at cloud top when the cloud has been further decoupled after cloud-base processes such as local fluctuations in LH beneath the cloud have eroded it from below.
Chapter 8 Decoupling of the cloud layer from the surface supply of moisture is known to be an important process in this region, especially as it relates to cloud transitions.

Discussion and Conclusions
We tested three diagnostic decoupling parameters suggested by Jones et al. (2011), as well as sub-cloud buoyant production of TKE following , to classify the simulated cloud cover across the decoupling continuum and found that the diagnostic based upon the vertical moisture gradient yielded the best results. The Lagrangian cloud profiles as the MBL cloud layer evolved suggested changes at cloud base and cloud top as the clouds became more broken, so we examined processes at cloud base and cloud top separately in the context of more coupled and more decoupled cloud populations over a domain at the onset of the transition from solid deck of marine stratocumulus to broken cloud structures, for each of the 2017 and 2018 simulations.
A finding of particular importance was that cloud-base erosion in the simulations was mainly the result of local fluctuations in LH, which preceded any changes at cloud top as the clouds became more broken. Correlations between the absolute value of LH and changes in the cloud fraction at cloud base were substantially lower than those between the Lagrangian derivatives of LH and cloud fraction at cloud base for more coupled clouds. Similar correlations to those involving LH were found when analyzing the LCL in the context of cloud-base changes, which is expected given that the surface LH and surface LCL are related. The relationship between the Lagrangian derivatives of 10-m wind velocity and the surface RH and LH following the cloud base were strong and positive for both coupled and decoupled clouds, uniformly exceeding 0.7, the implication being that changes in the surface wind and humidity conditions were likely influencing the cloud-base erosion process through their effect on LH as long as the cloud element was still coupled to the sub-cloud region. Changes in the surface relative humidity and LH were virtually simultaneous making it difficult to determine if the LH affected the humidity as a consequence of cloud-top mixing following deepening-warming or if surface warming led to a lower RH, which subsequently leads to an increase in the LH. In the latter case, a potential feedback loop would exist.
We examined the relationship between CTEI and EIS in the context of decoupling on cloud-top cloud fraction and found that CTEI was moderately correlated with changes in cloud-top cloud fraction only for more decoupled clouds, while EIS exhibited weak correlations for both coupled and decoupled clouds. In addition, we tested a new cloud-top dilution diagnostic based upon estimating the inversion strength using the vertical gradients of the difference between the equivalent and saturation equivalent potential temperature in the vicinity of the capping inversion height and found that it preformed similar to CTEI as an indicator of cloud breakup at cloud top, while being readily computed from soundings. In summary, the degree of decoupling was found to significantly influence the extent to which both cloud-base and cloud-top processes affect the erosion of cloud elements at cloud base and cloud top, respectively. Our findings suggest that cloud-top processes are able to substantially alter the cloud structure only after significant decoupling has occurred.
Our results suggest that the deepening-warming process proposed by Bretherton and  may be operating in the post-cold-frontal environment in ENA, which may broaden its applicability in the marine atmosphere. A variation of the deepening-warming hypothesis in which the warming initially dominates the moistening leading to a deeper MBL and higher LCL could also explain some of our findings. In that scenario, a higher LCL resulting from warming leads to moisture accumulation in the sub-cloud layer, which subsequently drives the LCL lower and initiates convection. We view this possibility as a potential subset of the original hypothesis rather than a separate process.
The latent-heat driven, deepening-warming decoupling process is a key insti- Future studies will hopefully address these other transition modes.
An important observation stemming from the WRF simulation experiments was that production and maintenance of clouds, specifically those at the wake of the cold fronts in both cases, proved particularly sensitive to the choice of the planetary boundary layer parameterization. Many test runs were conducted with UW PBL scheme (Bretherton and Park, 2009), and the prognostic TKE PBL schemes Mellor-Yamada Nakanishi and Niino Level 2.5 and 3 (MYNN2.5 and MYNN3) Niino, 2006, 2009), with only the latter capable of producing cloud structures that resembled the satellite images available at around the same time for both simulations. A series of studies by Nakanishi and Niino (Nakanishi and Niino (2004); Nakanishi and Niino (2006); Nakanishi and Niino (2009)) used LES results to improve the second-order turbulence model of Mellor and Yamada (Mellor and Yamada, 1974), leading to a better representation of the planetary boundary layer (PBL) and marine fog in a mesoscale model. In the absence of a detailed comparison between the performance of these two PBL schemes, we can only assume that it may partly be due to the prescription of turbulent master length scale, l, by MYNN, that incorporates buoyancy length scale in the specification of l proposed by Mellor and Yamada (1974), while the UW scheme follows the approach adopted by Mellor and Yamada (1974).
A basic answer to the broader question of the association between the deepeningwarming hypothesis, mesoscale clusters of cumulus-coupled stratocumulus, and air-mass modifications is at least partially answered by the present analysis of the two simulations. The transition process, as depicted by these simulations, has at least three definable stages: the initiation of intrusions of unsaturated air (i.e. clear patches) at the beginning of the transition process, a period of chaotic cloudiness with disorganized and somewhat random patches of hybrid combinations of stratocumulus and cumulus, and a final configuration in which the chaotic cloudiness becomes organized into scattered patches of mesoscale cumulus-coupled stratocumulus that are organized. The exact relationship between the deepening and warming process and the development of mesoscale organization is not addressed by the current analysis, though answers may be discovered in future analyses of these simulations.
One element of this study that deserves special mention is the compatibility of the simulated cloud structure to the observations at the ENA site. Past comparisons of model results with detailed process-level observations have often suffered from limitations that prevented the regional scale cloud evolution to be studied in a manner consistent with observations that could be used to track the integrity of the model solutions. For example, LES simulations could not capture regional variability within an air mass and regional simulations that featured resolutions on the scale of tens of kilometers could not resolve convection and, accordingly, were unable to reproduce the observed mesoscale structure. The simulations presented herein are among the first to bridge this disconnect; cloud-scale solutions that can be evaluated at a specific site, but represent the regional variability.
To the extent that these two cases represent different cold front climatology and produce transitions at different latitudes, they suggest a sensitivity to changing post-cold-frontal air mass characteristics in a warming climate. Warming in Northern Canada and an increasingly ice-free Arctic are likely to produce less impressive cold air outbreaks during summer and, at least initially, potentially smaller gradients of latent heat flux across the ENA. A reduction in the gradient of the latent heat flux in this region suggests that the stratocumulus transition might occur at lower latitudes than at present until the warming of the ocean, which takes longer, potentially offsets the process. In any event, a near-term reduction in the gradient, which seems quite feasible given the current amplitude and speed of the Arctic warming, could produce (or is already producing) a decoupling-induced, negative cloud feedback during the summertime.
Holistically, WRF simulations combined with observations from the ENA surface site is a potent research combination. The model provides hints at the specific data and process combinations that might reveal a much more statistically defensible modus operandi for the transition region when applied to a longer term record from the ENA site. Future combined observational and modeling studies are recommended to further clarify the connection between the surface fluxes, decoupling, and cloud transitions.
Appendix A

WRF Settings
Here, the namelist file used to run the 2018 simulation is presented, along with detailed explanations of variables not addressed in the "WRF Parameterizations" section, with the aim of assisting replication of simulations used in this study. The 2017 WRF simulation had identical settings, except for the start and end dates of the simulation. Important namelist options in each section that required further explanations are discussed in the following sections.
A.1 Time Control interval seconds: Refers to the interval between input files produced by WPS.
Since the initial and boundary condition FNL data obtained had intervals of three hours, the preprocessed input data by WPS (metgrid files) are also 3-hourly.
history interval: Indicates the intervals at which WRF will write out the simulation data. For the parent and child domains, intervals of 60 and 30 minutes, were used in our simulations, respectively.
frames per outfile: We chose a high enough number for this variable, because it was more manageable in this case when all the output intervals are written in a single file. The output files for the second domain for both simulations are about 1 TB, which did not pose any memory restraints for us.
restart: We experienced difficulty using the "restart" option, which in essence allows one to segment the simulation into different runs, so that when the program encounters an error, previous simulations prior to that time would not be lost.
The problem with using this option in our case was that the model would crash after only a few time steps when using a restart file. It must be noted here that a lot of debugging went into making sure that the model would run without any error for the entirety of the simulation run, which lasted about 9 hours for a 72-hr simulation.
A.2 Domains eta levels: The eta levels were obtained using WRF Domain Wizard software, version 2.84 (NOAA, 2013) , which is a very useful tool to easily set up a domain.
The eta levels were estimated using hyperbolic tangential stretching option, that creates uniformly spaced layers in the PBL, allowing them to gradually stretch apart to maximum distance at mid-levels (around 700 hPa), and then compresses them further aloft. The vertical spacing between the lowest levels were set to start at 15 m.
feedback: This setup determine a two-way or one-way nested run. In the former scenario, the values of the variables in the coarse domain are overwritten by those in the nest at points located at the intersection of the two domains, when this function is activated. Otherwise, the nesting process will be referred to as "one-way". A two-way nested run will be more computationally expensive, but ultimately more accurate than a one-way one. In our simulations, we benefited from this setting.
wif input opt: This option was activated to enable processing of Water/icefriendly aerosol input from metgrid, used in the Thompson aerosol aware microphysics scheme (mp physics = 28). respectively. Usually, a larger value is assigned to nproc y than the nproc x.
Arriving at the optimal decomposition was one of the biggest challenges we faced when running the simulations. It is closely dependent on the total number of computational nodes, size of the domains, and the output intervals designated to each domain. After many trial and errors, we arrived at the decompositions specified in the namelist. We used a total number of 150 compute nodes, (150×36 cores), dedicating 75 × 60 to computation and the rest, 75 × 12 to the I/O tasks, explained in "Namelist Quilt" section.
A.3 Physics use aero icbc: This option is used with aerosol-aware Thompson scheme with water-and ice-friendly aerosol climatology, and when set to "True", allows the model to use climatological aerosol input from WPS.
icloud bl: Activating this option allows subgrid-scale clouds from the PBL scheme to be coupled to the radiation schemes. This options is specific to MYNN2.5 and MYNN3 PBL schemes.
bl mynn cloudpdf: In order to represent subgrid-scale clouds, three probability distribution functions (PDF) are provided. Rather than resorting to multiple schemes to classify and represent different cloud types, which is a complex approach and normally does not lead to smooth cloud transitions, statistical methods using PDFs based on turbulence closure schemes (Mellor (1977); Sommeria and Deardorff (1977)) can be used to arrive at a unified scheme to represent all types of clouds. Kuwano-Yoshida et al. (2010) proposed an improved cloud PDF scheme based on the assumption of joint-Gaussian probabilities for the liquid water potential temperature and total water content, with standard deviations estimated from the Mellor-Yamada level-2 turbulence scheme, and improved turbulence closure parameters and mixing length. Because it is a diagnostic secondorder turbulence scheme that relies on resolved gradients, rather than higher order moments, it is computationally efficient, and results in improved representation of marine boundary layer clouds (Kuwano-Yoshida et al., 2010).
bl mynn mixlength: An improved, scale-aware mixing length formulation for MYNN proposed by Ito et al. (2015) has been used that takes into account the cloud-specific mixing length, as well. This scheme is designed with the aim of resolving turbulent eddies that are on the order of the model grid resolution (∼ 1 km), and has demonstrated improved performance in comparison to the MYNN3 scheme (Ito et al., 2015).
icloud: This option controls whether the effect of clouds is considered on the optical depth calculations in shortwave and longwave radiation, as well as the representation of cloud fraction in the model. The method proposed by Xu and Randall (1996) has been adopted here, which takes subgrid-scale cloud amounts into account by using statistical distribution of cloud water and cloud ice on the subgrid-scale, rather than threshold approaches that give either 0 or 1 values for cloud fraction.
o3input: When set equal to "2", this option inputs climatological Ozone data adapted from the Community Atmosphere Model (CAM) radiation scheme, as well as aerosol data for RRTMG radiation scheme. In contrast to the default Ozone data used in the scheme that is only a function of height, this Ozone dataset has additional latitudinal and temporal (monthly) variations.
A.4 Dynamics w damping: Vertical velocity damping has been activated to ensure the model remains numerically stable throughout the run, and avoid CFL errors due to locally large vertical velocities.
dampcoef A Rayleigh damper with damping coefficient (inverse damping time scale) of 0.2 is used in the top 5 km to damp spurious waves in the stratosphere.
use theta m: Activating this option incorporates the effect of moisture on pressure in small intervals.
A.5 Namelist Quilt nio tasks per group & nio groups: Quilting enables dedicating compute cores specifically to the task of writing out the simulation data at output intervals. This options works with MPI configuration. The total number of I/O tasks is equal to nio tasks per group × nio groups, where nio tasks per group is the number of processors used for I/O quilting per I/O group which is specified by nio groups. According to Balle and Johnsen (2016), nio tasks per group cannot be greater than nproc y, and ideally nproc y should be an exact multiple of nio tasks per group.
This I/O configuration, together with the compute task configuration (nproc x and nproc y) gave the best results in terms of reducing the output overheat at each output interval, and balancing I/O and compute resources.
TURB OPTION: Four options are available for parameterizing the turbulent wind field in the PBL; option 0 completely ignores turbulence, option 1 uses the Hanna scheme Hanna et al. (1982) for the PBL turbulent mixing, and option 2, the prognostic turbulent kinetic energy (TKE) from WRF. Option 3 is similar to option 2, but partitions TKE obtained from WRF such that turblent energy production and dissipation are equal.  found that options 2 and 3 cannot keep the tracer well-mixed in the planetary boundary layer. They found out that when the Hanna scheme for turbulence parameterization is used in conjuction with skewed turbulence (rather than Gaussian) in the convective boundary layer (activated by the swtich CBL = 1), the tracer particles are almost homogenously well-mixed throughtout the depth of the convective boundary layer. They recommend this setting and it is what we used in all trajectories in this article. LU OPTION: In this study, the same land-use scheme used in WRF has been implemented in FLEXPART-WRF with option 1.