Cloud-Top Entrainment in Mixed-Phase Stratocumulus and Its Process-Level Representation in Large-Eddy Simulation

: Cloud-top entrainment is a crucial process for the evolution of stratocumulus and is driven by interactions of radiation, microphysics, and turbulence on scales reaching down to less than one meter. Regardless of this fact, most large-eddy simulation studies still apply a horizontal resolution of tens of meters, not resolving these interactions sufﬁciently. Here, based on an extensive observational campaign, we deﬁne a weak-shear benchmark scenario for large-eddy simulation over Arctic ice and for the ﬁrst time perform large-eddy simulation of mixed-phase stratocumulus with horizontal resolutions of 35, 10, and 3.5 m. Thereby, we investigate the processes contributing to cloud-top entrainment and their role for the evolution of stratocumulus with a particular focus on resolution sensitivity. First, we ﬁnd that a horizontal grid spacing larger than 10 m insufﬁciently represents the effects of small-scale microphysical cooling and turbulent engulfment on cloud-top entrainment. Indeed, the small size of energy-containing eddies—a consequence of the intense stratiﬁcation in the vicinity of the cloud-top region—violates the underlying assumptions of subgridscale models by buoyant suppression of eddies at the large-eddy simulation ﬁlter scale. Second, the decrease in cloud-top entrainment due to these insufﬁciently represented processes results in 15% less cloud water after 6 h of simulation and a corresponding optical thinning of the cloud. Third, we show that the applied nonequilibrium microphysics cause microphysical heating beneath the cloud top, which partly counteracts the evaporative cooling.


Introduction
Stratocumulus (Sc) is the most common cloud type in the Arctic atmospheric boundary layer (ABL) (Eastman and Warren 2010). It may cover tens of thousands of square kilometers horizontally but is vertically confined to few hundreds of meters, often only tens of meters. Indeed, processes that drive Sc such as cloud-top entrainment and turbulence act on even smaller scales of meters to millimeters (Wood 2012;Mellado 2017). Consequently, mesoscale and climate models must parameterize Sc, thus enhancing corresponding uncertainties in such models (Tjernström et al. 2008;Pithan et al. 2014;Schneider et al. 2017). This delicacy is not only a consequence of Sc parameterization, but Sc also intensely interacts with other components of the Arctic climate system and sensitively responds to changes within them. The evolution and optical properties of Sc strongly depend on anthropogenic and natural aerosols (Possner et al. 2017;Solomon et al. 2018), surface latent and sensible heat fluxes (Gultepe et al. 2003;Solomon et al. 2014), and atmospheric state variables such as temperature, humidity, and wind (Wood 2012). In turn, its coverage and horizontal extent make Sc a key agent in Arctic climate change for its direct and indirect contribution to the surface energy and water budgets. These interactions between Sc and other components of the Arctic climate system result in feedback mechanisms (Kay et al. 2016;Goosse et al. 2018) such as the cloud-optical-depth feedback (Zelinka et al. 2012;Ceppi et al. 2017) and the cloud-sea ice feedback (Kay and Gettelman 2009;Morrison et al. 2018). While some feedbacks are well understood, others remain puzzling. Indeed, several interactions that are understood and quantified for strongly simplified scenarios become puzzling as soon as more complex setups are investigated.
Arctic Sc is matter of recent simulation-based studies (Klein et al. 2009;Pithan et al. 2014) as well as field campaigns (Verlinde et al. 2007;Tjernström et al. 2014) that affirm its importance for an understanding of the changing Arctic climate. These studies also reveal that there is a particular need for better understanding of processes driving Sc. Field campaigns are one option to improve this process-level understanding. In situ observations provide locally highly resolved data but often lack adequate spatial coverage. On the contrary, remote sensing offers sufficient spatial coverage but often lacks local resolution. Sometimes precision becomes an issue in both approaches. A tailor-made combination of them allows for deep insight to the processes from both a large-scale and local perspective. Recently, such a combination of different measurement techniques was shown to improve process-level understanding of Arctic Sc through the campaigns Arctic Cloud Observations Using Airborne Measurements during Polar Day (ACLOUD) and Physical Feedbacks of Arctic Boundary Layer, Sea Ice, Cloud and Aerosol (PASCAL) (Wendisch et al. 2017. Comparable field campaigns in the Arctic are scarce for observation-related challenges involved: low temperature, high relative humidity, and lack of daylight. In lower latitudes, such extensively equipped field campaigns are more common (Albrecht et al. 1988;Stevens et al. 2003;Zhou et al. 2004;Sorooshian et al. 2018). In the framework of these campaigns, large-eddy simulation (LES) is already a well-established research tool for process-level investigation of Sc (Stevens et al. 2005). Indeed, observation and LES fertilize and complement each other. LES applies microphysical parameterizations which are often calibrated using observational data. In turn, it may support field campaigns by forecasting conditions or suitable flight tracks, for instance, to reduce instrument icing or to ascertain adequate statistical coverage. The reason LES became such an essential tool for the investigation of Sc is that it allows one to modify in a straightforward and controlled manner the atmospheric state, cloud properties, and surface conditions. Such precise control of the initial and boundary conditions motivates the use of LES as a virtual laboratory where effects of particular changes can be studied and mechanistically attributed. While undoubtedly present in the field, such effects may be hardly, if at all, observable. In any case, the isolation of effects and cause-effect attribution is highly problematic in the field for the manifold interactions and broadscale nature of meteorological processes. Additionally, LES can unify the advantages of a high temporal and local resolution as well as sufficient spatial coverage. With recent advances in cost-and energy-efficient high-performance computing, the cost for such simulation is relatively low when compared to comprehensively equipped field campaigns, such as ACLOUD or PASCAL.
In practice, there is a dichotomy in numerical simulation between spatial coverage and high resolution that is overcome here by extensive use of computational resources. Indeed, one to two decades ago representation of cloud-top entrainment was problematic in many studies and known to be underresolved, even in carefully designed dedicated LES intercomparisons such as Stevens et al. (2005). Since an investigation of the idealized so-called smoke cloud , vertical resolution on the order of few meters and higher is applied. While an adequate process-level representation of cloud-top entrainment evidently requires a small vertical grid spacing, results concerning sensitivity toward horizontal resolution were contradictory at that time (Stevens and Bretherton 1999;Lewellen and Lewellen 1998). This ambiguity is shown to be a consequence of-consciously or unconsciously-tuning cloud-top entrainment by choice of grid resolution and aspect ratio to achieve a sufficient representation of the ABL and its cloudiness (Pedersen et al. 2016;Mellado et al. 2018). As discussed in Pedersen et al. (2016), higher horizontal resolution allows more small-scale turbulence which enhances cloud-top entrainment. Contrary, increased vertical resolution strengthens the capping inversion which hinders entrainment. This, however, depends on the definition of the capping inversion itself which can be either sharp or smooth. Pedersen et al. (2016) further show that domain size has only a minor effect on the results as long as it is large enough to encompass the largest updrafts. These results indicate that the use of strongly anisotropic small-sized grids substantially reduces computational demand but at the same time neglects physical processes that act on small scales but drive the system.
More recently, LES-based studies (Matheou and Chung 2014;Matheou 2018) and direct-numerical simulation studies (de Lozar and Mellado 2017;Mellado et al. 2018) investigated Sc with approximately isotropic grids up to a resolution of 1 m (Mellado et al. 2018). These studies provide deep process-level insight to cloud-top entrainment and demonstrate that specific configurations of LES may only work for particular purposes. For example, the lower-resolved tuned setups mentioned above are sufficient in a forecasting context where many of the relevant properties are governed by the boundary conditions (Schemann and Ebell 2020). But, their utility for process-level analysis is strongly limited. Although these numerical and also recent observational studies deepen our understanding of Sc cloud-top entrainment (Mellado 2017), most of them focus on simplified scenarios. They are either idealized or substantially simplified by their reliance on equilibrium microphysics, neglect of solar radiation, being at near-steady state, or confinement to a homogeneous surface or the open ocean. Comparable studies of more complex transitional scenarios are lacking to the authors' knowledge, but they are essential to eventually predict why Sc forms, breaks up, or vanishes.
LES studies on Arctic Sc are rare if compared to their lowerlatitude counterpart; to date, they rely on strongly anisotropic grids with low horizontal resolution. Most of them focus on the complex microphysical processes, i.e., the interaction of aerosols with distributions of liquid and solid hydrometeor species (Klein et al. 2009;Ovchinnikov et al. 2014;Solomon et al. 2014;Kaul et al. 2015;Young et al. 2017). Those sensitivity studies show that LES ensembles exhibit substantial spread for mixedphase scenarios while they tend to converge more for liquidonly scenarios. The complex interactions in multiphase microphysical processes have been identified as one cause for this intermodel spread. However, as mentioned before, studies at lower latitudes indicate that higher horizontal resolution is required to gain process-level representation of the clouddriving mechanisms. On top, the relatively small size of eddies in the Arctic ABL-as compared to its lower-latitude counterpart-also calls for higher resolution, in particular under neutral to stable conditions which frequently occur (Beare et al. 2006). Stable conditions are of particular relevance near the top of Sc where entrainment takes place. Thus, higher horizontal resolution ensures a realistic turbulence representation within the whole Arctic ABL and at the same time helps to physically understand cloud-driving processes such as cloud-top entrainment.
Existing studies on Arctic Sc call for extensive process-level research and identify particular challenges to traditional cloudresolving LES (such as Stevens et al. 2001Stevens et al. , 2005Neggers et al. 2011). First, as a consequence of stable stratification and shallow ABL height, characteristic eddies are smaller than in lower latitudes, which requires higher resolution. Second, Arctic Sc is predominantly of mixed phase, and therefore microphysical processes are more complex than in liquid-only Sc. In this study, we transfer the idea of applying highly resolved isotropic grids from the lower latitudes to a semi-idealized Arctic mixed-phase Sc scenario. While the chosen scenario is inspired by observations from ACLOUD and PASCAL, a best-possible reconstruction of these observations is not the primary goal of this study. Instead, we strive for a realistic setup that can be used as a benchmark scenario-and that can in principle be validated by observation-wherein to study clouddriving processes. This allows us to analyze the processes determining cloud-top entrainment in mixed-phase Sc and to investigate them through a resolution-sensitivity study. We hence focus on time-slab-averaged profiles of radiative and microphysical cooling rates and turbulence as well as their spatial variability and representativity.
In pursuit of this endeavor, we modify the LES mode of the Weather Research and Forecasting Model (WRF, version 4.0.3) to run at meter-scale resolution, to account for large-scale subsidence, and to appropriately initiate incloud turbulence (section 2). The modified model is then applied to the benchmark scenario defined for mixed-phase Sc over ice, inspired by observational data from the recent field campaigns ACLOUD and PASCAL (section 3). We analyze three simulations of this scenario at a horizontal resolution of 35, 10, and 3.5 m where the vertical resolution is about 3 m (section 4), and we discuss implications of our results (section 5). Full three-dimensional fields and statistics of each of the simulations are available in an open-data repository (Rauterkus and Ansorge 2019).

Model description and modification
WRF 4.0.3 solves the compressible, nonhydrostatic prognostic flux-form equations for dry air surface pressure; moist potential temperature u m ; geopotential f; Cartesian velocities u, y, and w [v 5 (u, y, w) T ]; and various scalars such as tracer or moisture variables (Skamarock et al. 2019). We use a terrain-following coordinate system and apply a fifth-order centered-difference advection scheme except for vertical scalar advection where we apply third-order centered differences. Following Pressel et al. (2017), we use monotonic advection for all scalar variables to eliminate unphysical drying above the cloud layer. Time integration is carried out by a third-order Runge-Kutta scheme separating low-frequency meteorological from high-frequency acoustic modes with the help of a time-split integration approach described in Wicker and Skamarock (2002). Although compressibility and some missing LES features count against the usage of WRF, it offers several advantages outweighing these drawbacks. In particular, WRF is open source and used by a large community. It runs and scales well across several computing architectures and provides the capability to simulate at the mesoscale. Thus, a transfer of results and parameterizations between LES and mesoscale simulation is straightforward. WRF also contains various wellestablished and documented physical parameterizations, especially for mixed-phase microphysics.
We use here the two-moment microphysical approach described in Morrison et al. (2009). It represents microphysical processes of four hydrometeor classes: cloud water, rain, cloud ice, and snow (we do not activate the graupel option). The cloud-droplet concentration N c remains as a fixed external parameter that can be used to compensate for problematic behavior of the microphysical parameterization. All other hydrometeor classes are described by prognostic mixing ratio and particle number concentration. Morrison and Pinto (2006), Morrison et al. (2008), and Hines and Bromwich (2017) point out the strengths of this parameterization and its sensitivity toward cloud droplet concentration and the amount of ice nuclei (IN) in the case of mixed-phase Sc. Based on observations, we set N c 5 75 cm 23 and IN 5 3.2 L 21 (Dupuy et al. 2019). As in Morrison et al. (2008) and following the default settings of the microphysical parameterization for Arctic Sc, we use IN about one order of magnitude higher than observed.
For radiation, we apply the RRTMG scheme (Iacono et al. 2008). We keep the surface properties constant throughout the simulation. Thus, we do not use a land model here and parameterize the surface fluxes via the revised version of the MM5 similarity scheme (Jiménez et al. 2012). A Smagorinsky model (Smagorinsky 1963) which is constrained by the Smagorinsky coefficient C S 5 0.25 and the Brunt-Väisälä frequency is applied for the subgridscale (SGS) turbulence closure. Based on this closure, we derive an estimator of ensemble-averaged SGS turbulent kinetic energy hTKE s i to assess the relative importance of unresolved turbulent mixing (appendix A), where S is the filtered rate of strain and D 5 ffiffiffiffiffiffiffiffiffiffiffiffiffi , with D x,y the horizontal resolution and D z the vertical resolution. Furthermore, we modify WRFs LES mode to allow for simulation of a mixedphase cloud at meter scale: 1) a second turbulence-seeding layer within the cloud, 2) a distinct microphysical time step, and 3) the inclusion of large-scale subsidence.
1) The default turbulence-seeding layer of WRF-LES adds to the four lowest model levels a pseudorandomized temperature perturbation field with values in the range of 60.05 K. This perturbation is constant within each column. In our case, it takes about 20-30 min for the shear-driven turbulence initiated at the surface to notably affect the whole ABL. After 1 h we consider the ABL to be well mixed. Both periods do not show notable sensitivity toward horizontal resolution in our scenario. Without additional turbulence seeding in the cloud, the cloud top might-due to radiative cooling-experience substantial artificial cooling and would just thicken until sheardriven turbulence initiated at the surface reaches the cloud. To counteract such artificial behavior, we implement a second turbulence-seeding layer within the cloud where a pseudorandomized temperature perturbation field with values in the range of 60.1 K is added to a level close below the expected cloud top; the same perturbation with opposite sign is added to the level below to approximately conserve the mean temperature thus keeping the amount of condensed water approximately constant. This additional turbulence seeding causes the development of convectively driven turbulence in the cloud within a few minutes. 2) We implement a distinct microphysical time step that acts similar to the radiative time step, which already exists in WRF. Thus, it allows for less frequent calls to the microphysical parameterization. This is not done here to speed up the calculation but to keep the simulation in a numerically stable regime. Without this distinct time step, the frequent calls to the microphysical parameterization would cause density-pressure waves which intensify until they violate the Courant-Friedrichs-Lewy condition at our fixed time step (Table 2). These waves originate from the frequent changes in hydrometeor species and temperature and occur only for time steps much smaller than 1 s where indeed assumptions underlying the microphysical scheme become questionable. 3) We represent large-scale subsidence by replacing w 1 w 1w in every prognostic equation in Skamarock et al. [2019, their Eqs. (2.8)-(2.14)], wherew represents the large-scale subsidence and F X the physical source term of a variable X (appendix B): with g the gravitational constant, m d being proportional to column mass of dry air, and capitalized letters being masscoupled variables, for example U 5 m d u. We apply largescale subsidence to all prognostic variables (including the microphysical ones), not just u, y, potential temperature, and scalars as in Yamaguchi and Feingold (2012). For largescale subsidence, we assume incompressibility and a constant horizontal divergence across the domain. Thus, we can link large-scale horizontal divergence D and large-scale subsidence: Within WRF, we implement the modifications in Eqs. (2a)-(2f) as a new physical parameterization, though it is a dynamic one, and discretize by Euler-forward differences in space. The temporal integration is carried out by the Runge-Kutta scheme mentioned above since the physical parameterizations only provide tendencies. We evaluated the modifications by applying the modified model to the well-investigated liquid-only Dynamics and Chemistry of Marine Stratocumulus (DYCOMS-II) RF01 scenario (Stevens et al. 2003). Results were compared to other LES studies of this scenario (Stevens et al. 2005;Yamaguchi and Feingold 2012;Matheou and Chung 2014), and with the presented setup we reproduce previous work. A corresponding resolution-sensitivity study produces results comparable to the ones we present in this paper.

Campaign and scenario description
The campaigns ACLOUD and PASCAL were carried out near Svalbard between 23 May and 26 June 2017 and focused on the Arctic ABL, in particular on clouds and their interaction with the surface. Variety and topicality of these observations set an ideal basis to define a benchmark scenario for investigation of Arctic mixed-phase Sc. We use here data from two research aircraft (Polar 5, Polar 6), one research vessel (RV Polarstern), and ground-based observations. Wendisch et al. (2019) give a summary of the first results and more detailed information on both campaigns and the instruments used. The radio sounding and aircraft data used in this work are described by Ehrlich et al. (2019); aircraft data on mean and turbulent meteorological quantities are made available by Hartmann et al. (2019), and Knudsen et al. (2018) complement this campaign overview by a summary of the observed synoptic development.
We define the benchmark scenario based on 5 June 2017 (Wendisch et al. 2019, their Fig. 18). It covers the period between two radio soundings (Schmithüsen 2017a,b), launched around 1100 and 1700 UTC from RV Polarstern. The scenario lies within the second synoptic key period as identified by Knudsen et al. (2018). Polar 5 and Polar 6 performed the research flight ACLOUD RF13 and flew multiple parallel triangular patterns close to RV Polarstern within, beneath, and above the cloud around 1200 UTC (Ehrlich et al. 2018a,b). Figure 1 sets the observations in relation to the position of RV Polarstern and gives more detailed information about location and time of the measurements we refer to. Further, it shows the low-level cloud deck in the ERA5 (European Centre for Medium-Range Weather Forecasts 2017), which observations confirm and identify as mostly liquid mixed-phase Sc between 100 and 500 m altitude.
As mentioned above, a best-possible convergence to the observations is not the objective of this study; rather, we focus on the representation and sensitivity of Sc-driving processes under realistic conditions. Thus, the initialization profiles are an idealized approximation of the 1100 UTC sounding, taking also other observations into account: where z i 5 470 m is the estimate for ABL height. Figure 2 shows observed potential temperature u and water vapor mixing ratio r y as compared to those calculated from Eqs. (4a) and (4b) assuming a moist-adiabatic cloud and well-mixed ABL. The radio sounding indicates weak wind in the lowest kilometer of the atmosphere which we approximate as u 5 1.5 m s 21 and y 5 20.6 m s 21 . The surface properties of the simulated cases are given in Table 1 along with a reference or data source, respectively. Contrary to the other parameters in this table, the latent and sensible heat fluxes are not an input parameter but instead calculated via the surface parameterization mentioned in section 2. They vary locally and temporarily throughout the simulations and also show slight differences between each simulation. These differences are, however, less than 1 W m 22 and can be estimated from the data provided in the accompanying data publication (Rauterkus and Ansorge 2019). Observational data of atmospheric heat fluxes is described by Egerer et al. (2019, their Fig. 11) and supports our weak positive surface fluxes around 11.5 W m 22 . We note that measurements are ambiguous regarding the existence of a moisture inversion above or near the top of the ABL. Even data of single methods do not allow for a definite conclusion regarding the existence or absence of such an inversion (red and blue dot in Fig. 2a). It thus remains open whether or not the moisture inversion (around z ' 400 m in airplane data and around z ' 500-600 m in the radio sounding) is an observational artifact (Dirksen et al. 2014), a local-, or a large-scale phenomenon as it may occur in the Arctic (Naakka et al. 2018). Recent LES studies indicate that a moisture inversion may evolve when Sc decouples from the surface. In this decoupled state of Sc, precipitation may lead to a drying of the ABL and therefore to a formation of a moisture inversion above (Neggers et al. 2019). For all observations available as well as throughout the simulation, the Sc is, however, always coupled to the surface. We thus decide to proceed without a specified moisture inversion in the initial condition which is consistent with a number of Polar 6 samples and the dropsonde.
We do not expect this interpretation of observational data to affect our results on process-level representativity. While we do expect different moisture fluxes, LES of Arctic Sc seems to differ at the quasi-equilibrium state only if no moisture source is available at all (Solomon et al. 2014). In our benchmark scenario, however, the surface latent heat flux provides moisture, and the Sc is located much closer to the surface than the one studied in Solomon et al. (2014). Additionally, Mellado et al. (2017) showed analytically for a cloud-free convective ABL that the moisture flux through the entrainment layer needs to be of similar magnitude and opposite sign to the one at the surface to have a drying effect on the ABL; this is not the case here, as in our benchmark scenario the magnitude of the moisture flux from the inversion would be much smaller due to the relatively small moisture jump.
To investigate the benchmark scenario, we simulate three cases that differ in horizontal resolution. The time steps-fixed in each simulation-are changed accordingly to satisfy the Courant-Friedrichs-Lewy condition, and the number of grid points is adapted such that the domain size remains constant ( Table 2). All three cases apply doubly periodic lateral boundary conditions and have an identical vertical resolution of D z , 3 m in and right above the ABL. In the nonturbulent region above the ABL, we gradually increase vertical resolution to about D z ' 20m at the domain top at about z top ' 850 m.

Results
We cover with the three cases one order of magnitude in horizontal resolution and simulate for the first time an Arctic Sc using an isotropic grid at the meter scale. We therefore focus on aspects that benefit from such high resolution (section 1): representation of turbulence (section 4b), representation of observed cloud and atmospheric state properties (section 4c), and scale dependence of individual processes of cloud-top entrainment (section 4d). Before we commence with the process-level analysis based on this high-resolution dataset, we assess our data from the perspective of representativity in the context of observations that were used to define the setup (section 4a). The simulation data consist of instantaneous three-dimensional snapshots every 5 min. Time series of the corresponding slab-averaged first three statistic moments and full three-dimensional snapshots at the initialization and end are made accessible by Rauterkus and Ansorge (2019).

a. Representativity of the results
If not otherwise mentioned, we use here time-slab-averaged profiles of the simulation data. The observations considered are, on the contrary, point measurements representative of a limited period and a finite volume. This benchmark scenario is idealized (section 3). In particular the geostrophic forcing and FIG. 1. Study area and flight path of Polar 6 during ACLOUD RF13. The shading represents low-level cloud cover (white: .95%; dark gray: ,50%) and the dashed line contours 50% sea ice coverage (where coverage increases northward), both derived from ERA5 at 1100 UTC. The observations shown in Fig. 2 are also symbolized: the green lines represent the flight track of Polar 6, the pink star marks the coordinates of the dropsonde (started from Polar 5), and the black star marks the coordinates of RV Polarstern where the radiosoundings were started and our simulation is placed.
large-scale divergence are constant in time, initialization profiles are smooth, and double periodic lateral boundary conditions prevent external lateral forcing of the ABL. These aspects require careful analysis and explanation of how the simulation results are related to the observations from 1) radio soundings and 2) Cloudnet data.
1) Two radio soundings (1100 and 1700 UTC) are used to initialize the simulation (section 3) and to compare the evolution of humidity, temperature, and wind between observation and simulation. They provide temporally and locally highly resolved data. However, they are instantaneous point measurements and thus affected by turbulent motion and local conditions. Both are filtered out when using time-slab averages as commonly done in the analysis of LES data. For comparison between model and radio sounding, we average over the last hour of each case, which corresponds to the period from 1600 to 1700 UTC. We do so because we force the model with a large-scale divergence averaged over the whole simulation period. Therefore, we expect the ABL height and related properties to match the observations best at the end of the simulation. 2) Cloudnet is an algorithm that provides various cloud properties based on observational data, such as data from Doppler cloud radar, microwave radiometer, and lidar (Illingworth et al. 2007). We use the derived liquid-and ice-water content, LWC and IWC. The Cloudnet output consists of time series and is provided by the Leibniz Institute for Tropospheric Research for the PASCAL campaign. Provided are point measurements in time and space, averaged over 30 s. However, Cloudnet data still exhibit a large temporal variability (Fig. 3a), manifest for instance in gaps between clouds (e.g., around 1330 UTC) and periods with lower liquid water path Observations and initial profiles of (a) water vapor mixing ratio r y and (b) potential temperature u. The green error bars mark the minimum and maximum of individual Polar 6 flight legs (lowest leg below cloud, highest leg above cloud, other legs within cloud) and the opaque red and blue dots encircle ascents or descents where the moisture inversion was observed or not observed by Polar 6. The measurements from Polar 6 are not single-valued due to multiple ascents and descents in the flight pattern.
LWP (e.g., from 1630 UTC onward). Time-slab averages from our LES do not capture this variability. While, for the large amount of data, we do not enable similar temporally resolved output for the LES, we can use Taylor's frozenturbulence hypothesis (Taylor 1938;Moin 2009) and convert a streamwise-sampled line from a horizontal slice to a time series using the height-local mean wind as advection velocity (Fig. 3b). We thus demonstrate that our simulation and the Cloudnet data share the same statistical behavior regarding their variability. Since simulation and Cloudnet data have comparable averaged LWP for the period from 1430 to 1530 UTC, we perform our comparison between both data series for that period.

b. Representation of turbulence
The general idea of LES is to resolve eddies up to the inertial subrange of turbulent motion where theories exist on how turbulence behaves and eventually dissipates. Thus, LES requires to explicitly resolve the eddies that carry most of the turbulent kinetic energy (TKE). The smallest eddies are filtered by low-pass filtering the governing equations which drastically reduces computational cost. Thereby, TKE can be separated in an unfiltered or resolved part, TKE r , and a filtered or SGS part, TKE s . The larger the contribution of TKE s is, the more important become the choices regarding the SGS model. The Smagorinsky model applied here (section 2, appendix A) simulates the SGS dissipation based on the stratification, resolved strain tensor, and equilibrium assumptions on the cascade of energy through the inertial subrange (Smagorinsky 1963;Pope 2000). That is, it requires a finite range of scales below the largest scales in the inertial subrange to be explicitly resolved throughout the simulated domain.
In most altitudes, TKE r contributes to about 95% to the total TKE, TKE t (Fig. 4c). This does, however, not hold for the cloud-top region where the amount of TKE r accounts for only about 19% of the total TKE in case #35.0m and reaches a maximum of about 80% in case #03.5m (Table 2). These values illustrate the substantial contribution of the SGS model to the simulated atmospheric state. Therefore, we assess whether or not the assumption of a resolved part of the inertial subrange is met for the three cases by analyzing their spectra of kinetic energy, given in Figs. 4a and 4b.
The wavelength beyond which numerical dissipation notably affects turbulent motion is-in accordance with Skamarock (2004)-about 7 times the grid resolution for both the total and the vertical energy spectrum (Figs. 4a,b). The share of the vertical fluctuation within the boundary layer varies about 40.0% (37.8%-47.5%) for the three cases and two heights (60 and 220 m above ground) investigated-illustrating the convective nature of our benchmark scenario, despite the  2. Simulation settings and cloud and turbulence properties averaged over the last hour of each simulation. For the simulation settings D x,y is the horizontal resolution, N x,y 3 N z are the amount of grid points, Dt ACC is the time step for the acoustical mode, Dt is the time step for fast physics, Dt MP for microphysics, and Dt R is the time step for radiation. For averaged properties, max(TKE s /TKE t ) is the vertical maximum of horizontally averaged SGS TKE ratio, LWP and IWP denote the horizontally averaged liquid and ice water paths, d y /d z is the horizontally averaged void fraction of the cloud, and d z is the horizontally averaged cloud depth. for case #10.0m, and up to 29.6% for case #03.5m. Although the suppression of turbulent motion is larger at small wavenumbers (large eddies) when compared to horizontal motion, it affects the whole spectrum and is caused by the strong stratification at cloud top which suppresses turbulence, especially in the vertical and for large scales. This is nicely illustrated by the wavenumber dependence of anisotropy a w (k) 5 E w (k)/E uyw (k) (Fig. 4c). We find a peak in the relative role of vertical motion in the spectrum around k r ' 0.01-0.03 m 21 , i.e., for scales on the order of 30-100 m that we can only simulate for the cases using a resolution of 10 and 3.5 m. This presence of vertical mixing is crucial for cloud-top entrainment and the maximum in this range of scales underlines the role of small-scale turbulence in supporting vertical penetration through the strongly stratified cloud-top interface: It shows that the resolved scales seem to properly feature the entrainment behavior only for the case with 3.5 m resolution. We, however, note that this seemingly crucial scale range is explicitly resolved for case #10.0m, too, which may already help because it suppresses inappropriate (diffusive) parameterization at this scale. Figures 4a and 4b show further that case #35.0m does not resolve the inertial subrange while case #03.5m resolves parts of it at wavenumbers k ' 3 3 10 22 m 21 , close before numerical dissipation gets notable. Case #10.0m resolves the inertial subrange, if at all, only barely with the inertial 25/3 power-law scaling being a tangent at wavenumbers around k ' 2 3 10 22 m 21 . In this range of wavenumbers, numerical dissipation starts to act and thus we cannot finally decide whether the beginning of the inertial 25/3 power-law scaling arises for physical or numerical reasons.

c. Bulk evolution of the simulations
Despite the indications that the SGS model cannot be fully trusted for the two lower-resolved cases, their profiles of atmospheric state variables evolve alike the ones of case #03.5m. Especially the profiles of potential temperature u and water vapor mixing ratio r y show no notable spread among each other (Fig. 5). Over time, the ABL moistens slightly and shrinks in height as it is also indicated by the observational trend. The cloud layer remains coupled to the underneath ice, and a moisture inversion does not form over the simulation period. Thus, our choice to not specify a moisture inversion is consistent with the relations used for initialization [Eqs. (4a) and (4b)] and also with the evolution of profiles over the course of simulation.
While u and r y agree well with the observations, there is a discrepancy between simulated and observed wind for two reasons. First, we compare instantaneous and local observations with time-slab-averaged profiles from LES (section 4a). Second, our simulation is at the beginning of a warm-air advection period. This large-scale phenomenon is not represented in our setup and the main reason for the different wind patterns between simulation and observations. During the simulation period, the warm-air advection sets in and the mean wind turns from northwesterly to southerly and increases from around 2 m s 21 to around 3 m s 21 within the ABL. Since the effect of warm-air advection is cumulative and it merely sets in, the differences in thermodynamic variables between observation and simulation remain small (Knudsen et al. 2018, their Fig. 4).
The primary ice particle formation takes place within the first hour and consumes water from the initialized liquid Sc (Fig. 6). This process is a consequence of our initial condition that does only specify total water mixing ratio in the form of liquid water instead of partitioning water to individual hydrometeor classes. It results in a simulated mixed-phase Sc that has a liquid mass fraction of about 95%. This fraction is a consequence of constraints from the atmospheric state, and it is also indicated by Cloudnet data. Over the course of the simulation LWP and the ice water path, IWP, decrease slowly in all three cases. This decrease is, however, not signifying a dissolving cloud here; it rather results from decreasing ABL and cloud height as consequence of the large-scale subsidence. This strength of the decrease directly results from our choice for an intense large-scale subsidence and is in accordance with the 1700 UTC radio sounding. All three cases show a similar evolution of the cloud quantities, and the simulated amplitudes of LWC and IWC match Cloudnet data in general FIG. 5. Radiosounding of (a) water vapor mixing ratio r y and (b) potential temperature u compared to profiles of the three simulations #35.0m, #10.0m, and #03.5m. Dashed horizontal lines symbolize cloud boundaries and the solid line the cloud maximum for case #03.5m, according to their definitions in Eq. (5).
FIG. 6. Temporal evolution of horizontally averaged simulated LWP (solid) and IWP (dashed) for the three simulated cases #35.0m, #10.0m, and #03.5m. DECEMBER 2020 R A U T E R K U S A N D A N S O R G E (Fig. 7). They, however, differ substantially below 150 m where Cloudnet data indicate near-zero IWC while all three cases predict ice particles reaching the surface. Reasons for this discrepancy can evolve from the microphysical parameterization through a too-high fall speed or a too-slow melt process of ice particles. The simulated precipitation rate (liquid and solid) is, however, smaller than 0.003 mm h 21 and thus very likely too low for detection by the Cloudnet procedures. Therefore, we conclude that the low IWC is a consequence of the correctly represented slow ice-particle sedimentation-a common feature of Arctic Sc. The vertical positioning of the Sc shows a small offset in LWC and IWC of a few tens of meters with respect to Cloudnet data. This offset is only one to two grid points given the vertical resolution of Cloudnet and thus mostly results from gridresolution uncertainties. It also is an averaging artifact as the cloud-top parameter derived from Cloudnet is much closer to the simulated cloud top (dashed horizontal lines in Fig. 7). Another potential cause for a slight vertical offset is our idealized setup (section 4a).
Comparing the different cases among each other, the evolution of the cloud quantities diverges between case #35.0m and the higher-resolved cases #10.0m and #03.5m. This difference increases the more the longer the simulation takes (Fig. 6). It illustrates the delicacy of process-level representation in cloud-top entrainment: while the profiles of temperature and humidity are to zero-order and first-order constrained by the ABL energetics, a physically accurate cloud representation requires a realistic coverage of the scales where the cloud-top entrainment occurs. This delicacy is an inherent problem when looking at ABL clouds as they only exist in the range of saturation from an averaged ABL perspective. After 6 h of the simulation, the difference results in a decrease of about 10%-15% in LWP and IWP. This decrease in the cloud water and ice is caused by a slightly drier and warmer ABL in the higher resolved cases. While differences in the profiles of temperature and humidity are minimal and barely visible in Fig. 5, the nonlinear relation to cloud quantities enhances them in the profiles of LWC and IWC in Fig. 7. Below a horizontal resolution of D x,y 5 10 m, we find that the cloud profiles stay nearly constant, which indicates their convergence in zero-and first-order statistics. Other profiles, such as temperature and humidity, already converge at an even lower resolution.

d. Cloud-top entrainment
While tiny given the precision of even the most recent observational tools, the difference in the cloud profiles between case #35.0m and the higher-resolved ones constitutes a substantial difference in cloud dynamics. This difference may cause dissipation of the cloud or formation of drizzle where this should not happen. Such misrepresentation would have an enormous impact on ABL dynamics in a prognostic setting which warrants further investigation. As stated in section 4c, the reduction in condensed water is related to a slightly warmer and drier ABL. We show in this section that this warming and drying of the ABL is due to enhanced cloud-top entrainment of air from above the ABL. This results in a slightly higher ABL, notable in all variables.
The main processes contributing to cloud-top entrainment are wind shear, radiative, and microphysical cooling (Mellado 2017). In this section, we analyze which of them causes the difference in cloud-top entrainment due to resolution and why this is so. Schulz and Mellado (2018) analyze in which scenarios wind shear is a driver of cloud-top entrainment. Since our scenario is convectively driven and has only a weak jump in wind velocity of about 0.2 m s 21 at the ABL top, we focus here on the radiative and microphysical cooling component of cloud-top entrainment, in form of their corresponding cooling rates R and E. To analyze the cloud structure, we follow Matheou (2018) and define cloud depth d z and void depth d y within a cloud: Here, z t and z b are the column cloud-top and cloud-bottom altitudes, where we use the condition r l 1 r s $ 0.01 g kg 21 for the existence of a cloud with r l and r s the liquid and solid water mixing ratios, respectively. Cloud depth is zero in the absence of a cloud and void depth is one. Their ratio d y /d z is a measure of how much dry air is engulfed within a column of the cloud as a consequence of cloud-top entrainment.
Profiles of R and E show that microphysical cooling is confined to the immediate vicinity of the cloud top while radiative cooling and microphysical heating occur below, namely, within the uppermost cloud region (Fig. 8a). Besides the slight microphysical heating in the cloud layer, this is in excellent agreement with theoretical considerations by Mellado (2017) and even more pronounced when considering the vertical coordinate z 2 z «,t , i.e., when using height relative to the columnlocal cloud top as vertical coordinate (Fig. 8c) with EPS the machine epsilon for a single-precision floatingpoint number. Columns with z «,t 5 0 m are not considered in the statistic moments and account for far less than 1% of the horizontal area across all three cases. Using the relative coordinate z 2 z «,t , we avoid averaging values at the cloud top with values from the middle of or outside the cloud. In consequence, the profiles averaged this way better reflect the local processes and peak more intensely. Notably, the profile of the microphysical cooling increases by a factor of about 5 and sharpens while the radiative cooling profile increases only by about 25% and keeps its shape (Figs. 8a,c). This different behavior among both processes results from underlying-different-physical mechanisms. Microphysical cooling and heating occurs where hydrometeor phase changes take place: condensation and evaporation are primarily confined to the cloud boundary. Hydrometeor phase changes within the cloud tend to balance each other due to up-and downdrafts (Fig. 10). Therefore, averaging values from different positions relative to the cloud top drastically weakens the resulting microphysical cooling at the cloud interface. Radiative cooling, on the contrary, is linked to the optical thickness of the cloud and thus not an interfacial process. It starts to act at the cloud top but strengthens as the cloud optical thickness increases since each cloud layer transmits, reflects, and emits radiation. This leads to a much smoother local vertical profile of radiative cooling in comparison with that of microphysical cooling. Except for a local maximum directly above the Sc, related to mostly frozen cloud fractions above the main part of the cloud, our radiative cooling profile follows the well-known shape of this quantity (Stevens et al. 2005;Mellado 2017). However, compared to lower-latitude Sc, its peak is less pronounced and the radiative cooling occurs over a deeper layer of the Sc. We attribute both phenomena to the LWP which is lower than for low-latitude Sc where cloud-optical thickness is higher. This illustrates a particular peculiarity of Arctic Sc, manifesting in strong dynamic differences when compared to low-latitude Sc, and it highlights the importance for dedicated studies investigating the dynamics of these clouds.
Increasing horizontal resolution from cases #35.0m to #10.0m results in intensified microphysical cooling (higher resolved gradients) and weakened radiative cooling (less cloud thickness). That means that higher resolution tends to intensify the cooling at the cloud top-both absolutely and in relation to other relevant processes-and therefore enhances the cloud-top entrainment due to buoyancy reversal (Mellado et al. 2009). Thereby, the corresponding differences between the cases #10.0m and #03.5m are smaller than the ones between cases #35.0m and #10.0m and mostly of opposite sign. We therefore doubt that all processes are sufficiently resolved in case #10.0m, although its time-slab-averaged zero-and first-order statistics do not differ from the ones of case #03.5m. A nonlinear interaction between turbulence and microphysics on scales around D x,y 5 10 m that is not parameterized must cause the opposite sign of the changes between the individual cases' radiative and microphysical cooling rates. For further investigation of this feature, we turn to investigation of spatial variability and microphysics.
Instantaneous values of microphysical cooling at the cloud top locally reach values of up to 2300 K h 21 but also microphysical heating at a rate of up to 100 K h 21 (Fig. 9a). These values are FIG. 9. Probability distributions in absolute numbers of (a) radiative and (b) microphysical cooling rates (R and E, respectively) from case #03.5m, accumulated over the last hour of simulation. The bin sizes for R and E are 0.00625 and 0.667 K h 21 , respectively. The dashed horizontal lines symbolize cloud bottom and top and the solid line marks the cloud maximum for case #03.5m, according to their definitions in Eq. (5). strongly variable in time and space and occur in all three cases, i.e., irrespective of horizontal resolution. They result from different microphysical processes in up-and downdrafts where vertical velocities reach values of about 61 m s 21 in combination with the time ''lag'' of the microphysical parameterization. In updrafts, the lag intensifies supersaturation of an air parcel being lifted until it enters the temperature inversion above the cloud (Fig. 10). Thus, in the uppermost cloud layer, supersaturation is largest and causes microphysical heating around 5-10 m below the cloud top. In the inversion, all condensed water evaporates and causes intense cooling. Air masses being engulfed from above into the cloud are warm and also enhance evaporation locally. Thus, they also contribute to the microphysical cooling maximum at the cloud top. At lower altitude, the microphysical time-slabaveraged tendencies are about one order of magnitude smaller and the diabatic cooling in downdrafts tends to approximately balance the diabatic heating in updrafts, resulting in a consistent, albeit small, microphysical heating of around 0.5-1 K h 21 across most of the cloud layer (Fig. 10). The radiative cooling is almost constant within the middle and lower cloud layers and exhibits only small spatial variation as compared to microphysical cooling within the upper cloud layers (Fig. 9).
In section 4b we discuss the requirement of high resolution to resolve the energy spectrum sufficiently. The same holds for the spatial variability of cloud-top cooling (Figs. 8b,d and 9), which in turn triggers turbulent motion. The higher the resolution is, the larger are horizontal gradients and in-plane variances. This effect is of particular importance at the cloud top where a strong stratification suppresses turbulence and thus requires intense eddies to break through it and induce cloudtop entrainment (Garcia and Mellado 2014). The same process also takes place at lower height. However, there the ABL is neutrally stratified and well mixed. In this regime, higher variability in cooling may induce or enhance eddies, it is but of no benefit for the cloud-top entrainment since the turbulence in the well-mixed layer is of larger scale and already more intense.
Considering the energy spectra (Fig. 4) and void-depth ratio (Fig. 11b) altogether reveals another reason for increased cloud-top entrainment with higher resolution: Assuming an eddy to require about seven grid points (section 4b; Skamarock 2004), case #35.0m can only resolve eddies which are larger than '250 m while case #03.5m resolves eddies down to '25 m. The encroachment of air above the ABL into the Sc occurs to a large part in very localized and small regions at the cloud edges which often span only few tens of meters (Fig. 11). Since cloudtop entrainment is the only process that can trigger this encroachment, we suggest that resolving these regions is of particular importance to resolve cloud-top entrainment.
To analyze if and how strong encroachment depends on the horizontal resolution, we consider the production terms of › t u 02 , i.e., of temperature variance: Here, X 0 denotes the horizontal deviation of the variable X from its mean, and R, E, A, D, S, and N refer to the cooling rates induced by radiation (R), microphysics (E), advection (A), large-scale subsidence (D), the SGS model (S), and metric and acoustic terms (N), respectively. We give the corresponding profiles in Fig. 12, not showing the profiles related to metric terms and acoustics (N). Vertical advection induced by large-scale subsidence is the strongest contributor to › t u and also contributes notably to › t u 02 , but it is insensitive to horizontal resolution. The term N, estimated as a residuum, shows sensitivity to horizontal resolution but it is very small in magnitude compared to the other terms and thus not considered henceforth. As discussed before, microphysical and radiative cooling produce temperature variance in heights around 400-420 and 370-410 m, respectively. This production is for both processes close to cloud top but it is only the microphysical process which shows a strong sensitivity to horizontal resolution. Figure 12 shows that both advection and the SGS consume temperature variance by mixing air masses thus reducing vertical gradients at the cloud top. The mixing of air masses into the ABL is called encroachment and is dominated by vertical motion in our benchmark scenario. Further, advective and SGS elimination of temperature variance show a high sensitivity toward horizontal resolution. Since the sensitivity is of opposite sign and comparable magnitude to the one of microphysical temperature variance production, we conclude that it is the microphysical cooling which causes enhanced entrainment by encroachment on both the resolved scales and the parameterized SGS. In a quantitative sense, we find an approximately threefold increase of averaged void-depth ratio d y /d z from case #35.0m to case #10.0m and a further approximately 20% increase from case #10.0m to case #03.5m (Table 2) we interpret the much smaller difference from case #10.0m to case #03.5m as an indication of convergence. All three cases show a cloud cover exceeding 99% and thus the increase in cloud-void ratio illustrates that the process is only marginally represented in the low-resolved case #35.0m.
In conclusion, our results show that, although the time-slabaveraged profiles of many atmospheric quantities converge for a horizontal resolution of 10 m and below, the simulation does not converge in terms of process representativity. Indeed, processes have to be taken into account at such high resolutions which are not of importance at lower resolution. Thereby we observe a beginning process-level convergence around our finest resolved case #03.5m since the changes from case #35.0m to case #10.0m are much more intense than the ones from case #10.0m to case #03.5m.

Discussion and implications
We defined and simulated a benchmark scenario for Arctic mixed-phase Sc based on the recent field campaigns ACLOUD and PASCAL. Therefore, we modified WRF-LES 4.0.3 to simulate mixed-phase Sc at meter-scale resolution and taking into account large-scale subsidence in all prognostic variables. A detailed analysis of the microphysical cooling profiles unveils a characteristic ''S'' shape with slight warming through microphysics just below the cloud top. This shape differs from the prototypical profiles that emerged from equilibrium microphysical considerations over the past years (Mellado 2017). We suggest that this is due to (i) accumulated supersaturation in updrafts being released close to the cloud-top region and (ii) different microphysical processes acting in updrafts and downdrafts. Both aspects illustrate the importance of nonequilibrium microphysical effects if such subtleties of the cloud-top region are of interest.
Our resolution-sensitivity study covers one magnitude in horizontal resolution and includes an approximately isotropic grid at D x,y 5 3.5 m and D z , 3 m. We show that increasing the horizontal resolution from 35 to 3.5 m results in more spatial variability which improves three aspects of the processlevel representativity of an LES with regard to cloud-top entrainment.
d First, it allows us to resolve parts of the inertial cascade of turbulent motion without notable numerical dissipation. This holds for the whole domain but it is particularly crucial for the cloud top where the model heavily relies on the SGS parameterization. The wavenumber dependence of anisotropy shown in section 4b suggests that resolution of motion at scales O(10) m is crucial to correctly simulate entrainment; parameterization of this behavior based on resolved-scale properties does not seem possible using standard turbulence parameterization. At the same time, the convergence of bulk and entrainment statistics beginning for case #10.0m DECEMBER 2020 R A U T E R K U S A N D A N S O R G E (sections 4c and 4d) also suggests that it is sufficient to resolve part of the dynamics at this scale-on the order of few meters in the case of the benchmark scenario introduced here-for a correct representation of entrainment.
d Second, the extended range of resolved turbulence allows smaller turbulent eddies which in turn produce more engulfment of free-tropospheric air into the cloud layer.
d Third, a higher resolution increases the resolved spatial variability in the atmospheric variables. These variations increase by nonlinear relationships to an even higher spatial variability in microphysical cooling which, in turn, causes buoyant instability and therefore triggers turbulent eddies.
All three aspects are closely connected through production of and being driven by turbulence which outlines the importance of a valid SGS parameterization. Further, we find that the dominant role of cloud-top interfacial microphysical processes found here is a peculiarity of the Arctic environment where low optical thickness of the cloud renders radiation as clouddriving process overall relatively less important in comparison with low-latitude Sc. We show-in accordance with Mellado et al. (2018)-that LES can converge in zero-order and first-order statistics of timeslab-averaged variables even though underlying assumptions are violated or crucial processes nonresolved. In our scenario, this applies for a horizontal resolution of 10 m and below. Coarser resolutions, still frequently applied, cause a slightly cooler and more moist atmosphere and therefore an increase in liquidand ice-water paths LWP and IWP of about 15% after 6 h of simulation. Those differences are rather small but can lead to a collapse or breakup of Sc if located close to its transition regime. Considering observational uncertainty and the fact that other variables vary even less due to resolution, we agree that lower resolution still can get the mean profiles right. However, these profiles are mainly constrained by boundary conditions and external forcing. We thus suggest that ''correct'' slab-averaged profiles are not necessarily proof for the physical correctness of process-level representation in LES.
Acknowledgments. RR gratefully acknowledges the funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)-project ID 268020496-TRR 172-within the Transregional Collaborative Research Center ''ArctiC Amplification: Climate Relevant Atmospheric and SurfaCe Processes, and Feedback Mechanisms (AC) 3 .'' CA was partly funded by the German Research Council (DFG) through the Transregional Research Collaborative 32 ''Patterns in Land-Surface-Atmosphere Systems (TR32),'' and supported by the European Research Council (ERC; Grant Agreement 853147) and through a UoC PostDoc Grant by the University of Cologne as part of its institutional strategy for the support of early stage researchers within the excellence initiative. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time through the John von Neumann Institute for Computing (NIC) project HKU24 on the GCS Supercomputers JUQUEEN and JUWELS at Jülich Supercomputing Centre (JSC). We further thank the Alfred Wegener Institute (AWI), the PS106/1 crew, and the ACLOUD science teams for making the field campaign happen. We especially thank Ulrich Löhnert and Yaping Shao for their input to our work and many informative discussions, Christof Lüpkes for providing and discussing noseboom data from Polar 6, Regis Dupuy and Stephan Mertes for providing information on aerosols, and Johannes Stapf for providing information on the surface properties.
Data availability statement. The datasets generated and analyzed during the current study are available from the accompanying data publication (Rauterkus and Ansorge 2019).

Calculation of Subgridscale Turbulent Kinetic Energy
The Smagorinsky closure expresses SGS diffusivity based on the strain tensor and thus does not carry any explicit information regarding the SGS TKE. We, however, believe that this is an important quantity to assess whether or not the LES hypotheses are valid and give in this appendix an estimate of the SGS TKE that is consistent with the Smagorinsky closure by applying some integral relations introduced in Pope (2000) and considering them in the local limit.
We estimate the ensemble-averaged SGS TKE hTKE s i by integration of the TKE spectrum E(k) from the cutoff wavenumber k c 5 1/D to infinity, i.e., where C 5 1.5 and D 5 ffiffiffiffiffiffiffiffiffiffiffiffiffi D 2 x,y D z 3 q . Assuming-in accordance with the underlying assumptions of the Smagorinsky closure-local balance of SGS turbulence, we can replace dissipation « by the production hPi and express it by the SGS production term from the Smagorinsky model [Pope 2000, their Eq. (13.129)] to yield where S is the filtered rate of strain and C S the Smagorinsky constant. Note that this replacement of dissipation at SGS by resolved-scale production involves the assumption that all dissipation happens at the SGS. We note that Eq. (A2) is only valid for the ensemble average, but the right-hand side can be evaluated point-wise to yield local estimates of TKE s . These local estimates provide an unbiased estimate, i.e., the slabaveraged TKE s may be approximated to arbitrary precision when the sample size is increased.

Inclusion of Large-Scale Subsidence
For the inclusion of large-scale subsidencew, we add it to the vertical Cartesian velocity w in every prognostic equation: with m d being proportional to column mass of dry air and v 5 (u, y, w) T where u and y are the horizontal Cartesian velocities, maps to since › z m d 5 0 for the terrain-following coordinate. Accordingly, their Eq. (2.11) for moist potential temperature u m with F X the physical source term of a variable X and the capitalized letter being the corresponding mass-coupled variable X 5 m d x, maps to and substituting the modified continuity equation results in The last summand in this expression describes the effect of the introduced large-scale subsidence on the tendency of moist potential temperature. Since the inclusion of large-scale subsidence directly affects the vertical velocity, this modification would require a modification of the dynamic core of WRF-LES. In agreement with work done before Yamaguchi and Feingold (2012), we included the effect of large-scale subsidence as a physical source term, thus For the other prognostic variables similar expressions can be found and Eqs. (2.8)-(2.14) in Skamarock et al. (2019) read as F W 1 F W 2 m d (w› z w 1w› zw 1 w› zw ) , with g the gravitational constant, f the geopotential, and q m the different mixing ratios, here for cloud ice, snow, cloud water, and rain.