The Earth Model Column Collaboratory (EMC²) v1.1: an open-source ground-based lidar and radar instrument simulator and subcolumn generator for large-scale models

Climate models are essential for our comprehensive understanding of Earth's atmosphere and can provide critical insights on future changes decades ahead. Because of these critical roles, today's climate models are continuously being developed and evaluated using constraining observations and measurements obtained by satellites, airborne, and ground-based instruments. Instrument simulators can provide a bridge between the measured or retrieved quantities and their sampling in models and field observations while considering instrument sensitivity limitations. Here we present the Earth Model Column Collaboratory (EMC²), an open-source ground-based lidar and radar instrument simulator and subcolumn generator, specifically designed for large-scale models, in particular climate models, but also applicable to high-resolution model output. EMC² provides a flexible framework enabling direct comparison of model output with ground-based observations, including generation of subcolumns that may statistically represent finer model spatial resolutions. In addition, EMC² emulates ground-based (and air- or space-borne) measurements while remaining faithful to large-scale models' physical assumptions implemented in their cloud or radiation schemes. The simulator uses either single particle or bulk particle size distribution lookup tables, depending on the selected scheme approach, to perform the forward calculations. To facilitate model evaluation, EMC² also includes three hydrometeor classification methods, namely, radar- and sounding-based cloud and precipitation detection and classification, lidar-based phase classification, and a Cloud Feedback Model Intercomparison Project Observational Simulator Package (COSP) lidar simulator emulator. The software is written in Python, is easy to use, and can be straightforwardly customized for different models, radars, and lidars.

Following the description of the logic, functionality, features, and software structure of EMC², we present a case study of highly supercooled mixed-phase cloud based on measurements from the U.S. Department of Energy Atmospheric Radiation Measurement (ARM) West Antarctic Radiation Experiment (AWARE). We compare observations with the application of EMC² to outputs from four configurations of the NASA Goddard Institute for Space Studies (GISS) climate model (ModelE3) in single-column model (SCM) mode and from a large-eddy simulation (LES) model. We show that two of the four ModelE3 configurations can form and maintain highly supercooled precipitating cloud for several hours, consistent with observations and LES. While our focus is on one of these ModelE3 configurations, which performed slightly better in this case study, both of these configurations and the LES results post-processed with EMC² generally provide reasonable agreement with observed lidar and radar variables. As briefly demonstrated here, EMC² can provide a lightweight and flexible framework for comparing the results of both large-scale and high-resolution models directly with observations, with relatively little overhead and multiple options for achieving consistency with model microphysical or radiation scheme physics.

1 Introduction

The representation of cloud processes in large-scale models is continuously advancing, conceptually, and in the level of details and complexity implemented in the micro- and macro-physical schemes (e.g., Lin et al., 2019; Cesana et al., 2019). These improvements are reflected in the accuracy of the resulting model output (e.g., Klein et al., 2013; Lin et al., 2019; Myers et al., 2021; Wang et al., 2019), yet results still show large inter-model variability (e.g., Zelinka et al., 2020). This variability results from, among other sources, model weaknesses concerning their ability to predict atmospheric state variables and processes such as cloud geometrical and optical thicknesses (e.g., Cesana and Waliser, 2016; Klein et al., 2013), formation and transition of marine stratocumulus clouds (e.g., Rémillard and Tselioudis, 2015; Lin et al., 2014; Cesana et al., 2019), and the formation and maintenance of supercooled water (e.g., Cesana et al., 2012; Tan and Storelvmo, 2016; Silber et al., 2019b).

Meaningful model evaluation benefits from a direct (“apples-to-apples”) comparison with observations (e.g., Bodas-Salcedo et al., 2014; Suzuki et al., 2015). For the evaluation of cloud representation, model output is often compared with active remote-sensing measurements from instrumentation such as lidars and radars, which provide information on the spatial structure of clouds and some direct indications about active microphysical processes. However, model evaluation is challenging because of observational detectability constraints (e.g., signal extinction) and lack of retrievals or large uncertainties in some microphysical and atmospheric state quantities by these instruments, for example, hydrometeor number concentration or water content. In addition, spatial resolution differences between a model's output and an observing instrument's measurement resolution present an additional difficulty.

To bridge the gap between large-scale models such as weather or climate models and active remote-sensing observations, instrument simulators with different purposes have been developed to estimate observed parameters using model output. For example, the Cloud-resolving model Radar SIMulator (CR-SIM; Oue et al., 2020) was developed to emulate zenith-pointing and scanning radar and lidar variables using high-resolution model output, with considerations of hydrometeor shape and the resulting scattering calculations (see also Mech et al., 2020). The Cloud Feedback Model Intercomparison Project Observational Simulator Package (COSP; Bodas-Salcedo et al., 2011; Swales et al., 2018), on the other hand, was developed to operate over large-scale model output targeting satellite data as observational constraints, although expansions for the emulation of ground-based radars and lidars have been developed (e.g., Zhang et al., 2018; Kuma et al., 2020). Among a variety of reasons such as the demanding computation associated with the emulation of satellite measurements and the ability to use and output detailed data to and from simulators, COSP is typically implemented on-line into models' code to facilitate output.

To account for spatial resolution discrepancies, which are typically accentuated in the case of large-scale models due to their coarser resolution, some model evaluation studies and forward simulators emulate a higher spatial resolution by generating subcolumns (e.g., Bodas-Salcedo et al., 2008, 2011; Chepfer et al., 2008; Klein and Jakob, 1999; Lamer, 2019; Stephens et al., 2010; Swales et al., 2018; Webb et al., 2001). Statistics calculated using multiple generated subcolumns, which are nominally faithful to the processed model's physics, can be directly compared with the associated observations, thereby mitigating spatial resolution biases and errors.

Here we present the Earth Model Column Collaboratory (EMC²), an open-source ground-based lidar and radar simulator and subcolumn generator, which is designed to operate over large-scale model output while being faithful to the physics implemented in models' microphysics or radiation schemes but can also be applied to high-resolution model output. EMC² enables detailed evaluation of atmospheric thermodynamic profile and cloud properties extracted from local, regional, and global simulation outputs against long-term ground-based, air- or space-borne datasets. The software is written in Python, allowing quick installation and providing customizable operation (for scattering calculations, etc.). In Sect. 2, we briefly illustrate the Python code and software followed by a detailed description of the subcolumn generator, the forward calculations, and classification routines currently implemented in EMC². In Sect. 3 we demonstrate the use of EMC² by comparing observations of a case study of highly supercooled drizzling cloud over West Antarctica (see Silber et al., 2019a) with application of EMC² to outputs from a large-eddy simulation (LES) and the NASA Goddard Institute for Space Studies (GISS) ModelE3 climate model (see Cesana et al., 2019) in single-column model (SCM) mode. Finally, Sect. 4 provides a summary of the code features and case study demonstration.

2 EMC² description

2.1 Software description

EMC² depicts a workflow for comparing forward calculated radar or lidar variables generated from large-scale model output with radar or lidar measurements. Figure 1 shows a flowchart example of this workflow for using EMC² to compare ModelE3's output with high spectral resolution lidar (HSRL) measurements. The workflow starts with the Model class incorporated within the emc2.core module. The Model class contains model output field namelists and default hydrometeor parameters (Table 1). Using Python's class inheritance, EMC² allows the creation of a custom class specifying a given model's namelists and parameters (a ModelE3 class in this example), which ensures that the model output can be standardized and used by the other modules in EMC². Once loaded through the Model class internal methods, model output data are stored within the Model object using the xarray dataset module (Hoyer and Hamman, 2017).

Figure 1 Standard workflow of EMC² utilization for direct comparison between observations and model output. The flowchart exemplifies the use of EMC² to compare the NASA GISS ModelE3 climate model output with corresponding high spectral resolution lidar (HSRL) measurements. N_s designates the specified number of subcolumns.

The Model object is then input to the subcolumn generator (Sect. 2.2). Note that subcolumn generator processing can be practically skipped by setting the number of subcolumns (N_s) to 1. The results of the subcolumn generator are stored in the xarray dataset contained within the Model object. Here, we introduce the Instrument class. Similar to the Model class, the Instrument class contains relevant information about the instrument being simulated (some of which is listed in Table 2) as well as the single-particle and bulk scattering calculation lookup tables (LUTs) (see Sect. 2.3 and 2.4). Currently, zenith-pointing instrument properties and scattering calculation LUTs are available for various lidars and radars operated by the Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) climate research facility: the elastic micropulse lidar (MPL) and HSRL, both operating at a wavelength of 532 nm (e.g., Flynn et al., 2007; Eloranta, 2005), the 910 nm CL31 ceilometer (Morris, 2016), the 1064 nm HSRL elastic channel (see Razenkov and Eloranta, 2018), the elastic channel of ARM's Raman lidar operating at 355 nm (e.g., Newsom, 2009), the C-band scanning ARM precipitation radar (C-SAPR; starting in EMC² v. 1.2; Widener and Bharadwaj, 2012), the X-band scanning ARM cloud radar (XSACR; Widener et al., 2012b), the Ka-band ARM zenith radar (KAZR; Widener et al., 2012a), and the W-band ARM cloud radar (WACR and M-WACR; Widener and Johnson, 2006). Similar to the Model class, the Instrument class allows EMC² to be tailored to other radars and lidars deployed at different sites, which does not confine the analysis of measurements and model output to specific ARM instruments or sites, as long as the required parameters and suitable scattering LUTs are provided. Thus, the various parameters mentioned in the next subsections (e.g., all parameters shown in Tables 1 and 2) as well as the LUTs can be easily specified and set to match configurations and assumptions implemented in different large-scale models, as well as complex scattering models more commonly implemented in cloud-resolving and LES models.

Thorsen and Fu (2015, Table 3)Nott and Duck (2011)Silber et al. (2019a)Morrison and Gettelman (2008, Table 2)

Table 1 Hydrometeor class parameter values implemented in EMC².

Following the subcolumn generator process, the Instrument and Model objects are then input to the lidar (Sect. 2.3) or radar (Sect. 2.4) simulators (lidar simulator in the case of Fig. 1). The forward calculation results are stored in the same xarray dataset in the Model object. Simulated hydrometeor classification (Sect. 2.5) can be performed following the completion of the forward calculations and stored in the xarray dataset. For comparison and visualization of these results, EMC² uses the Atmospheric Community Toolkit (ACT; Theisen et al., 2020). Thus, a SubcolumnDisplay object, inherited from ACT's Display object, contains the necessary methods for quick visualization of the simulated instrument variables. In addition, the SubcolumnDisplay object also contains several internal methods for generating curtain and profile plots of observational and simulated data stored in the Instrument or Model objects, allowing masking of simulated signals below instrument detectability, for example. The figures presented in the next section (Sect. 3) show examples of EMC²'s visualization capabilities. Finally, since the data are in the xarray dataset format, EMC² also contains all of xarray's analysis and visualization capabilities for these simulated datasets.

Widener and Johnson (2006); Widener et al. (2012a, b)Widener and Johnson (2006); Widener et al. (2012a, b)Widener and Mead (2004)Falconi et al. (2018)Burns et al. (2016)Falconi et al. (2018)Burns et al. (2016)

Table 2 Radar instruments and some of their characteristics currently implemented in EMC². The ARM SGP, ENA, NSA, AWR, and MOS site abbreviations denote the Southern Great Plains, Eastern North Atlantic, North Slope of Alaska, the AWARE campaign (at McMurdo Station, Antarctica), and the Multidisciplinary Drifting Observatory for the Study of Arctic Climate (MOSAiC) Expedition, respectively. The calculations of the minimum detectable signal ($Z_{{\mathrm{e}}_{min}}$) for KAZR are based on the analysis in Silber et al. (2018a). All $Z_{{\mathrm{e}}_{min}}$ values correspond to 2 s integration time except for the AWR and MOS WACR, the values of which correspond to 0.2 s integration time.

n/a – not applicable

EMC² incorporates a suite of unit tests for each function using the pytest testing tool (https://pytext.readthedocs.io/, last access: 1 December 2021) to inspect the integrity and functionality of the code. These unit tests are combined with continuous code integration using TravisCI integration service (https://travis-ci.com/, last access: 1 December 2021), which runs the unit tests every time a developer submits a pull request on GitHub. If the unit test passes with the developer's changes to the code, then the changes are approved to be a part of EMC². Documentation is also automatically generated from the metadata strings in each subroutine to ensure that each part of the code is well documented.

2.2 Allocation of hydrometeors to subcolumns

The following simulator description assumes large-scale model convective and/or stratiform cloud scheme outputs containing four hydrometeor classes: cloud water (cl), cloud ice (ci), rain (precipitating liquid; pl), and snow (precipitating ice; pi). While these four hydrometeor classes are widely used in various microphysics schemes, we note that EMC² can be generally adapted to cases in which fewer or additional classes are used. We note that all the acronyms, abbreviations, and symbols used throughout this study are listed in Appendix A.

The subcolumn allocation and forward calculations in EMC² can be performed using two main approaches; the first follows the assumptions and general logic implemented in large-scale model radiation schemes while the second follows the assumptions and logic in model microphysics schemes. In short, the radiation scheme approach largely treats hydrometeor fractions in a generalized manner and utilizes bulk (hydrometeor population) scattering lookup tables (LUTs) generated using specific size distribution assumptions. The microphysics scheme approach, on the other hand, includes full integration of single-particle scattering LUTs using hydrometeor particle size distributions prognosed by models with consideration of sub-grid hydrometeor class fraction variability assumptions.

We note that the detailed description below of the radiation and microphysics approaches is congruent with the current implementation of these approaches in the GISS ModelE3 climate model. However, the core of these approaches is similar in other climate and Earth system models (ESMs), and EMC² can be easily adapted to fit specific variations in a model assumptions (see Sect. 2.1). For example, the microphysics approach currently operates only on stratiform microphysics scheme output using a two-moment bulk scheme (Gettelman and Morrison, 2015, hereafter MG2) that has been implemented in climate models such as the Community Earth System Model version 2 (CESM2) Community Atmosphere Model version 6 (CAM6) (Danabasoglu et al., 2020), the Energy Exascale Earth System Model (E3SM; Golaz et al., 2019), and ModelE3.

Prior to the radar and lidar forward calculations, EMC² generates subcolumns for each model output column. These subcolumns emulate a higher model spatial resolution, which partially reconciles the locality of ground-based measurements and allows a more robust statistical model evaluation. Subcolumns are generated and populated with hydrometeors from the top-down using the maximum-random overlap approach (Tian and Curry, 1989; Fan et al., 2011; Hillman et al., 2018) similar to the description by Lamer (2019, chap. 6). EMC² translates hydrometeor fractions in the model grid to a binary set of hydrometeor-containing and hydrometeor-free subcolumn bins. That is, given a specified number of subcolumns (N_s; determined by the user), the total number of hydrometeor-filled subcolumn bins at model level h and time step t is equal to the rounded value of $N_{\mathrm{s}}\times f_{\mathrm{hyd}}(h,t)$, where f_hyd is the volume fraction of a hydrometeor class (e.g., f_cl, f_pi) or a generalized hydrometeor fraction (f_gen) used in the model radiation scheme, at the same model level and time step. Here and henceforth, we assume for simplicity an SCM output (no horizontal coordinate dimensions), even though EMC² can generally operate not only on SCM simulation output but also on global simulation output. Note that EMC² does not calculate f_gen, but can use it if it is a model output field. The f_gen parameter is a generalization of cloud and precipitation fractions denoting air volume in which scattering and absorption by hydrometeors influence radiative calculations for a model. (In the case of ModelE3, for example, f_gen corresponds to cloud fraction where the layer-average cloud opacity exceeds that of precipitation, or vice versa.)

The following steps are applied in order to populate subcolumns (number $\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{s}}$) with hydrometeors:

1. Convective cloud hydrometeors (cl and ci) are allocated to the first subcolumns (lowest index; if $f_{\mathrm{hyd}}(h,t)>\mathrm{0}$), thereby generating cloud-containing subcolumns with maximum convective geometrical cloud depths.

2. Stratiform cloud hydrometeors are allocated to subcolumn bins unoccupied by convective cloud hydrometeors following the maximum-random overlap approach. Thus, stratiform clouds at model level h are first randomly allocated to subcolumn bins with overlying stratiform clouds (at level h+1), followed, if necessary, by random allocation to subcolumns with cloud-free bins directly above at level h+1. This order of processing where clouds are preferentially extended vertically conforms with assumptions that are often implemented in large-scale model radiative transfer calculations. (Note that our objective is to be faithful to model physics, even though in the case of the maximum-random overlap approach, it could produce, in certain cases, larger model biases compared to observations (e.g., Hillman et al., 2018).) During this hydrometeor allocation step, maximum overlap between liquid and ice cloud hydrometeors is maintained in grid cells containing both cloud phases. For example, if $f_{\mathrm{ci}}(h,t)>f_{\mathrm{cl}}(h,t)>\mathrm{0}$, all cl-containing subcolumn bins belonging to the same grid cell will necessarily contain ci as well (mixed-phase), but there will be additional subcolumn bins in that grid cell containing only ci hydrometeors. In the case where overlying stratiform hydrometeors exist, the overlying layer phase is not a factor of consideration, such that a subcolumn bin containing ci may be located right above a subcolumn bin containing cl or both cl and ci (mixed-phase), and vice versa.

3. Convective and stratiform precipitating hydrometeors (pl and pi) are allocated to subcolumns without convective-stratiform restrictions, such that convective and stratiform precipitation may co-exist in the same subcolumn bin. Similar to the stratiform cloud allocation, precipitation is allocated with maximum-random overlap; i.e., precipitation is first extended vertically. Thus, either stratiform or convective cloud hydrometeors may exist in the same subcolumn bins as convective and/or stratiform precipitating hydrometeors (for example, describing a scenario of snow generated by a stratiform cloud seeding an underlying drizzling stratocumulus cloud field). If after the vertical extension of overlying precipitation some subcolumn bins are still to be populated with precipitating hydrometeors (e.g., in a case where $f_{\mathrm{pl}}(h,t)=f_{\mathrm{pl}}(h+\mathrm{1},t)=\mathrm{0}$ and $f_{\mathrm{pi}}(h,t)>f_{\mathrm{pi}}(h+\mathrm{1},t)>\mathrm{0}$), these hydrometeors are randomly allocated to cloudy grid cells of the same type (convective or stratiform), followed by random allocation of hydrometeors to cloud-free subcolumn bins. Similar to the allocation of stratiform cloud hydrometeors, maximum overlap is maintained between liquid and ice precipitation, and overlying precipitation phase is not a factor of consideration during the step of vertical extension of precipitation.

Once the subcolumns are populated with hydrometeors, per hydrometeor class except for stratiform cl in the case of the microphysics approach, hydrometeor mixing ratio is set in every hydrometeor-containing subcolumn bin by $q_{\mathrm{hyd}}(s,h,t)={\overline{q}}_{\mathrm{hyd}}(h,t)/f_{\mathrm{hyd}}(h,t)$, where q_hyd and ${\overline{q}}_{\mathrm{hyd}}$ designate the mixing ratio of a hydrometeor class (e.g., q_cl, q_ci) in subcolumn bin s and a corresponding model grid cell mean, respectively. In the case of cl when using the microphysics approach, at every model level h and time step t, $q_{\mathrm{cl}}(s,h,t)$ is randomly set in cl-containing subcolumn bins such that it would comply with the sub-grid variability gamma distribution described by Morrison and Gettelman (2008, Eq. 8) while adjusting the values in the last cl-containing subcolumn bin such that hydrometeor mass is conserved (as in the case of other hydrometeor classes), i.e.,

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{N_{\mathrm{s}}}{q}_{\mathrm{hyd}}(i,h,t)}{N_{\mathrm{s}}}}={\displaystyle \frac{{\overline{q}}_{\mathrm{hyd}}(h,t)}{f_{\mathrm{hyd}}(h,t)}}.\end{array}$$

We note that sub-grid scale variability of cloud water in ModelE3 is tied to the sub-grid scale variability of moisture rather than set at a fixed value as in Morrison and Gettelman (2008).

In stratiform hydrometeor-containing subcolumn bins, hydrometeor number concentration is set for every hydrometeor class by $N_{\mathrm{hyd}}(s,h,t)={\overline{N}}_{\mathrm{hyd}}(h,t)/f_{\mathrm{hyd}}(h,t)$, where N_hyd and ${\overline{N}}_{\mathrm{hyd}}$ designate the number concentration of a hydrometeor class in subcolumn bin s and a corresponding model grid cell, respectively. EMC² assumes that convective schemes do not diagnose ${\overline{N}}_{\mathrm{hyd}}$, and hence, this information is currently not produced by the simulator.

2.3 Forward calculation of lidar variables

2.3.1 Microphysics approach

In the microphysics approach, applicable only to stratiform hydrometeors, per hydrometeor diameter D, EMC² calculates the hydrometeor size distribution, ${\mathit{\varphi }}_{\mathrm{hyd}}(D,s,h,t)$, defined by q_hyd and N_hyd, fully consistent with the MG2 scheme (see Morrison et al., 2009, Eqs. 1–3). Using LUTs containing full Mie calculations for spheres (following Bohren and Huffman, 1983, Appendix A) of single particle extinction and backscatter efficiencies at lidar operating wavelength λ_l ($Q_{{\mathrm{e}}_{\mathrm{hyd}}}(D,{\mathit{\lambda }}_{\mathrm{l}})$ and $Q_{{\mathrm{bs}}_{\mathrm{hyd}}}(D,{\mathit{\lambda }}_{\mathrm{l}})$, respectively), the lidar particulate extinction cross-section (${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}(s,h,t)$) and backscatter cross-section (${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{hyd}}}(s,h,t)$) are calculated in every hydrometeor-bearing subcolumn bin by

$$\begin{array}{}\text{(2a)}& \begin{array}{rl}{\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}(s,h,t)& ={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\underset{D_{min}}{\overset{D_{max}}{\int }}{\mathit{\varphi }}_{\mathrm{hyd}}(D,s,h,t)Q_{{\mathrm{e}}_{\mathrm{hyd}}}(D,{\mathit{\lambda }}_{\mathrm{l}})D^{\mathrm{2}}\mathrm{d}D\\ & \approx {\displaystyle \frac{\mathit{\pi }}{\mathrm{8}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{D}}-\mathrm{1}}\left({\mathit{\varphi }}_{\mathrm{hyd}}\right(D_i,s,h,t\left)Q_{{\mathrm{e}}_{\mathrm{hyd}}}\right(D_i,{\mathit{\lambda }}_{\mathrm{l}})D_i^{\mathrm{2}}\\ & +{\mathit{\varphi }}_{\mathrm{hyd}}(D_{i+\mathrm{1}},s,h,t)Q_{{\mathrm{e}}_{\mathrm{hyd}}}\\ & (D_{i+\mathrm{1}},{\mathit{\lambda }}_{\mathrm{l}})D_{i+\mathrm{1}}^{\mathrm{2}})\mathrm{\Delta }{D}_{i,i+\mathrm{1}}\end{array}\end{array}$$

$$\begin{array}{}\text{(2b)}& \begin{array}{rl}{\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{hyd}}}(s,h,t)& ={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\underset{D_{min}}{\overset{D_{max}}{\int }}{\mathit{\varphi }}_{\mathrm{hyd}}(D,s,h,t)Q_{{\mathrm{bs}}_{\mathrm{hyd}}}(D,{\mathit{\lambda }}_{\mathrm{l}})D^{\mathrm{2}}\mathrm{d}D\\ & \approx {\displaystyle \frac{\mathit{\pi }}{\mathrm{8}}}\sum _{i=\mathrm{1}}^{N_{\mathrm{D}}-\mathrm{1}}\left({\mathit{\varphi }}_{\mathrm{hyd}}\right(D_i,s,h,t\left)Q_{{\mathrm{bs}}_{\mathrm{hyd}}}\right(D_i,{\mathit{\lambda }}_{\mathrm{l}})D_i^{\mathrm{2}}\\ & +{\mathit{\varphi }}_{\mathrm{hyd}}(D_{i+\mathrm{1}},s,h,t)Q_{{\mathrm{bs}}_{\mathrm{hyd}}}\\ & (D_{i+\mathrm{1}},{\mathit{\lambda }}_{\mathrm{l}})D_{i+\mathrm{1}}^{\mathrm{2}})\mathrm{\Delta }{D}_{i,i+\mathrm{1}},\end{array}\end{array}$$

where we use the trapezoidal rule for discrete integration over a series of D values, which can be unevenly spaced by $\mathrm{\Delta }{D}_{i,i+\mathrm{1}}=D_{i+\mathrm{1}}-D_i$, while noting that $D_{\mathrm{1}}=D_{min}$ and $D_{N_{\mathrm{D}}}=D_{max}$, where N_D is the number of diameters for which $Q_{{\mathrm{e}}_{\mathrm{hyd}}}$, $Q_{{\mathrm{bs}}_{\mathrm{hyd}}}$, and ϕ_hyd are calculated. In the case of the Mie calculations currently available in EMC², D₁=0.1 µm, $D_{N_{\mathrm{D}}}=\mathrm{1}$ cm, and $\mathrm{\Delta }{D}_{i,i+\mathrm{1}}$ is constant and equals 0.1 µm. The complex refractive indices (m_hyd(λ_l)) used for liquid hydrometeors in the Mie calculations can be taken from Segelstein (1981, Table 1) or Rowe et al. (2020, for a temperature of −10 ^∘C). Refractive indices for ice hydrometeors are taken from Warren and Brandt (2008). The Maxwell Garnett equation (Bohren and Battan, 1980, Eq. 1) for a mixture of ice and air is used to calculate the effective m_hyd for ci and pi based on the ice densities implemented in EMC² for ModelE3 (Table 1) relative to bulk ice density of 917 kg m⁻³.

The total α_p and β_p (${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{tot}}}(s,h,t)$ and ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}(s,h,t)$, respectively) are calculated as the sum of each of these variables for cl, ci, pl, and pi. The lidar linear depolarization ratio (LDR) is estimated by weighting fixed LDR values (per hydrometeor class; see Table 1) with the relative contribution of ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ to ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$.

The cumulative optical thickness (from the surface upward) at the base of a given subcolumn bin, τ_hyd, is calculated by

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{\mathrm{hyd}}(s,h,t)=\sum _{i=\mathrm{2}}^h{\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}(s,i-\mathrm{1},t)\mathrm{\Delta }z(i-\mathrm{1},t),\end{array}$$

where Δz(h,t) denotes the geometrical thickness of model level h at time step t and ${\mathit{\tau }}_{\mathrm{hyd}}(s,h=\mathrm{1},t)=\mathrm{0}$. The total integrated optical thickness, ${\mathit{\tau }}_{\mathrm{tot}}(s,h,t)$, being the sum of τ_hyd for cl, ci, pl, and pi, is used to estimate the level of full lidar signal attenuation (received signal not detectable by the simulated instrument), the value of which can be used to constrain comparisons between model output forward calculations and observations. Lidar signal extinction at visible wavelengths typically occurs at an optical thickness of 3–5 (e.g., Sokolowsky et al., 2020, Fig. 4), and hence, EMC² assumes by default that the lidar signal is extinct at a level where τ_tot=4 to produce a lidar signal extinction mask. We note that EMC² allows calculating τ_hyd from the top-down (i.e., at the top of a given model layer), thereby enabling simulation of airborne and spaceborne lidar measurements.

In some cases, the observational dataset only consists of measurements made by elastic lidars that do not allow direct retrieval of α_p and β_p. In these cases, observations can be compared with a diagnostic attenuated ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$, ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot},\mathrm{att}}}$, which is often referred to as the normalized relative backscatter (NRB; Campbell et al., 2002) after several lidar signal corrections are applied. ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot},\mathrm{att}}}$ is calculated here by

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot},\mathrm{att}}}(s,h,t)& =\left[{\mathit{\beta }}_{\mathrm{p}}\right(s,h,t)+{\mathit{\beta }}_{\mathrm{m}}(s,h,t\left)\right]\\ & T_{\mathrm{m}}^{\mathrm{2}}(s,h,t){\exp }^{-\mathrm{2}\mathit{\eta }{\mathit{\tau }}_{\mathrm{tot}}(s,h,t)},\end{array}\end{array}$$

where β_m and $T_{\mathrm{m}}^{\mathrm{2}}$ are the molecular backscatter cross-section and the two-way transmittance calculated following Penndorf (1957), and η is the multiple scattering coefficient. The value of η is set to 1 by default, effectively implying no multiple scattering by hydrometeors or perfect multiple scattering treatment. Such an assumption is not realistic in the vast majority of cases even though effective η values closer to 1 are more likely given the common narrow field-of-view of ground-based (and airborne) lidars. While a value of 1 is most likely unrealistic, it precludes the introduction of additional confounding factors stemming from the dependence of multiple scattering effects on the hydrometeor particle size distributions (e.g., Eloranta, 1998), which could significantly impact model–observation comparisons of ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot},\mathrm{att}}}$ (not shown). Moreover, this dependence of multiple scattering on particle size distributions suggests that a fixed η value smaller than 1 would impact the faithfulness of the lidar simulator to the model physics. That said, we note that η can be manually set to other fixed values based on the physical assumptions made or certain empirical results (e.g., Winker, 2003), and that the determination of the extinction level based on τ_tot is independent of η.

2.3.2 Radiation approach

With the radiation approach, applicable to both stratiform and convective cloud scheme output, forward calculations rely on bulk scattering LUTs. Therefore, this approach is more than 2 orders of magnitude faster than the microphysics approach, thereby rendering EMC² more suitable for the analysis of large model output datasets. Using this approach, a geometric cross-sectional area for each hydrometeor-bearing subcolumn bin is first calculated assuming geometric scatterers:

$$\begin{array}{}\text{(5)}& A_{\mathrm{hyd}}(s,h,t)={\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}{\displaystyle \frac{q_{\mathrm{hyd}}(s,h,t){\mathit{\rho }}_a(h,t)}{{\mathit{\rho }}_{\mathrm{b}}{r}_{{\mathrm{e}}_{\mathrm{hyd}}}(h,t)}},\end{array}$$

where ρ_a is the density of air, ρ_b is the bulk density of the hydrometeor class phase (1000 and 917 kg m⁻³ is the case of liquid and solid water, respectively), and $r_{{\mathrm{e}}_{\mathrm{hyd}}}$ is the effective radius of a hydrometeor class in the model grid cell, provided as a model output field. The ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ and ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ are then calculated by

$$\begin{array}{}\text{(6a)}& {\displaystyle }{\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}(s,h,t)=Q_{\mathrm{e},{\mathrm{vol}}_{\mathrm{hyd}}}\left(r_{{\mathrm{e}}_{\mathrm{hyd}}}^{*}\right(h,t),{\mathit{\lambda }}_{\mathrm{l}})A_{\mathrm{hyd}}(s,h,t)\text{(6b)}& {\displaystyle }{\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{hyd}}}(s,h,t)=Q_{\mathrm{bs},{\mathrm{vol}}_{\mathrm{hyd}}}\left(r_{{\mathrm{e}}_{\mathrm{hyd}}}^{*}\right(h,t),{\mathit{\lambda }}_{\mathrm{l}})A_{\mathrm{hyd}}(s,h,t),\end{array}$$

where $Q_{\mathrm{e},{\mathrm{vol}}_{\mathrm{hyd}}}\left({\mathit{\lambda }}_{\mathrm{l}}\right)$ and $Q_{\mathrm{bs},{\mathrm{vol}}_{\mathrm{hyd}}}\left({\mathit{\lambda }}_{\mathrm{l}}\right)$ represent in this case bulk efficiencies per unit volume taken from LUTs, in which they are provided as function of $r_{{\mathrm{e}}_{\mathrm{hyd}}}$. In Eqs. (6a) and (6b), however, $Q_{\mathrm{e},{\mathrm{vol}}_{\mathrm{hyd}}}$ and $Q_{\mathrm{bs},{\mathrm{vol}}_{\mathrm{hyd}}}$ are functions of the adjusted effective radius, $r_{{\mathrm{e}}_{\mathrm{hyd}}}^{*}$, which equals $r_{{\mathrm{e}}_{\mathrm{hyd}}}$ in all hydrometeor classes and cloud types except for stratiform cloud ice and snow. In these two cases, $r_{{\mathrm{e}}_{\mathrm{hyd}}}^{*}=r_{{\mathrm{e}}_{\mathrm{hyd}}}{\mathrm{\Phi }}_{\mathrm{hyd}}^{*}$ and ${\mathrm{\Phi }}_{\mathrm{hyd}}^{*}={\mathrm{\Phi }}_{\mathrm{hyd}}\frac{{\mathit{\rho }}_{\mathrm{hyd}}}{{\mathit{\rho }}_{\mathrm{b}}}+(\mathrm{1}-{\mathrm{\Phi }}_{\mathrm{hyd}})(\frac{{\mathit{\rho }}_{\mathrm{hyd}}}{{\mathit{\rho }}_{\mathrm{b}}}{)}^{\frac{\mathrm{1}}{\mathrm{3}}}$, where ρ_hyd is the hydrometeor class density (see Table 1), and Φ_hyd is a constant fluffiness factor. Φ_hyd is used such that a value of 0 represents an equivalent mass bulk sphere while a value of 1 represents a fluffy sphere with an equivalent maximum dimension as in Gettelman and Morrison (2015). In the case of ModelE3, for example, Φ_hyd is set by default to an intermediate value of 0.5 and generally serves as one of many tuning parameters

The default $Q_{\mathrm{e},{\mathrm{vol}}_{\mathrm{hyd}}}$ and $Q_{\mathrm{bs},{\mathrm{vol}}_{\mathrm{hyd}}}$ in Eqs. (6a) and (6b) (respectively) implemented in EMC² were calculated using single-particle full Mie calculations in the case of liquid hydrometeors and single-particle scattering LUTs for a severely roughened eight-column ice aggregate (Yang et al., 2013) in the case of solid hydrometeors. These ice aggregate scattering calculations have been shown by Holz et al. (2016) to enable a closure between infrared Moderate-Resolution Imaging Spectroradiometer (MODIS; Platnick et al., 2003) and visible Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP; Winker et al., 2009) satellite ice optical thickness retrievals and were included in the MODIS collection 6 (C6) cloud product (Platnick et al., 2017). In order to calculate the $Q_{\mathrm{e},{\mathrm{vol}}_{\mathrm{hyd}}}$ and $Q_{\mathrm{bs},{\mathrm{vol}}_{\mathrm{hyd}}}$ LUTs, we assumed the same gamma distribution parameters as those implemented in the C6 dataset (see Hansen, 1971, Eq. 1), consistent with the bulk LUTs utilized by ModelE3's radiation scheme.

Following the calculations of ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ and ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{hyd}}}$, the total variables ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{tot}}}$ and ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$ as well as τ_hyd, τ_tot, and ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot},\mathrm{att}}}$ are calculated as in the microphysics approach (Sect. 2.3.1).

2.4 Forward calculation of radar variables

2.4.1 Microphysics approach

When the microphysics approach is selected in EMC², the first three radar moments are calculated for each hydrometeor class in every hydrometeor-bearing subcolumn bin, that is, the equivalent reflectivity factor ($Z_{{\mathrm{e}}_{\mathrm{hyd}}}$), the mean Doppler velocity ($V_{D_{\mathrm{hyd}}}$) and the Doppler spectral width (${\mathit{\sigma }}_{D_{\mathrm{hyd}}}$) as well as total radar moment variables ($Z_{{\mathrm{e}}_{\mathrm{tot}}}$, $V_{D_{\mathrm{tot}}}$, and ${\mathit{\sigma }}_{D_{\mathrm{tot}}}$). Full Mie calculation LUTs for the emulated radars (Table 2) are first used to calculate ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ at the radar operating wavelength λ_r following Eq. (2b). The m_hyd(λ_r) used for liquid in the Mie calculations can be taken from Segelstein (1981, Table 1) or Turner et al. (2016, using a temperature of −10 ^∘C), while m_hyd(λ_r) values for ice are taken from Mätzler (2006, chap. 5.3; using temperature of −10 ^∘C). $Z_{{\mathrm{e}}_{\mathrm{hyd}}}$ in every hydrometeor-containing subcolumn bin is then calculated (in linear units) using (Doviak and Zrnić, 1993, Eq. 4.33)

$$\begin{array}{}\text{(7)}& Z_{{\mathrm{e}}_{\mathrm{hyd}}}(s,h,t)={\displaystyle \frac{{\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{hyd}}}(s,h,t){\mathit{\lambda }}_{\mathrm{r}}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}|K_{\mathrm{w}}{|}^{\mathrm{2}}}},\end{array}$$

where $|K_{\mathrm{w}}{|}^{\mathrm{2}}$ is the dielectric factor for water used in the raw radar observational processing (see Table 2). Using the resultant $Z_{{\mathrm{e}}_{\mathrm{hyd}}}$, $V_{D_{\mathrm{hyd}}}$ is then calculated by implementing the hydrometeor class terminal velocities parametrization used in the MG2 scheme (cf. Morrison and Gettelman, 2008, Table 2). In the calculation of $V_{D_{\mathrm{hyd}}}$ we neglect the model grid cell vertical wind, w, which predominantly has little impact on the $V_{D_{\mathrm{hyd}}}$ value, especially at coarser spatial and temporal resolutions typical to large-scale models. Finally, ${\mathit{\sigma }}_{D_{\mathrm{hyd}}}$ is calculated, while we note that beamwidth and turbulent broadening (e.g., Chen et al., 2018) are omitted from this calculation but will be added in future work.

To allow a valid comparison between the forward calculations and observations, signal attenuation is considered in the calculation of the attenuated $Z_{{\mathrm{e}}_{\mathrm{tot}}}$, $Z_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}$:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{Z}_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}(s,h,t)& =\mathrm{10}{\log }_{\mathrm{10}}\left(Z_{{\mathrm{e}}_{\mathrm{tot}}}\right(s,h,t\left)\right)\\ & -\mathrm{2}\left(Y_{{\mathrm{hyd}}_{\mathrm{tot}}}\right(s,h,t)+Y_{\mathrm{gas}}(s,h,t\left)\right),\end{array}\end{array}$$

where Y_gas and $Y_{{\mathrm{hyd}}_{\mathrm{tot}}}$ are the one-way integrated attenuation at the base of the subcolumn bin (in dB) due to atmospheric gases (O₂ and H₂O; see Ulaby et al., 1981, Sect. 5.3–5.5) and all hydrometeors, respectively. $Y_{{\mathrm{hyd}}_{\mathrm{tot}}}$ is calculated using

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{Y}_{{\mathrm{hyd}}_{\mathrm{tot}}}(s,h,t)& =\mathrm{10}{\log }_{\mathrm{10}}\left(e\right)\sum _{i=\mathrm{2}}^h{\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{tot}}}(s,i-\mathrm{1},t)\\ & \mathrm{\Delta }z(i-\mathrm{1},t),\end{array}\end{array}$$

where ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{tot}}}$ is determined by summing ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ based on Eq. (2a) calculated at λ_r over all hydrometeor classes while setting $Y_{{\mathrm{hyd}}_{\mathrm{tot}}}(s,h=\mathrm{1},t)=\mathrm{0}$.

The vertical profile of the minimum detectable equivalent reflectivity factor, $Z_{{\mathrm{e}}_{min}}$, is calculated by

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{Z}_{{\mathrm{e}}_{min}}(h,t)& =Z_{{\mathrm{e}}_{min}}(z=\mathrm{1000}\mathrm{m})\\ & +\mathrm{20}{\log }_{\mathrm{10}}\left(z\right(h,t)/\mathrm{1000}),\end{array}\end{array}$$

where z(h,t) is the height at the base of model level h (in meters) at time step t and $Z_{{\mathrm{e}}_{min}}(z=\mathrm{1000}\mathrm{m})$ is the minimum detectable signal at 1 km (using the highest sensitivity mode; Table 2). When compared with observations, subcolumn bins where $Z_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}(s,h,t)<Z_{{\mathrm{e}}_{min}}(h,t)$ can be treated as returned signal below the radar noise floor and, hence, are effectively considered hydrometeor-free.

2.4.2 Radiation approach

When the radiation approach is selected in EMC², forward radar calculations using bulk LUTs are limited to the zeroth radar moment (Z_e) due to a set of limitations:

1. Large-scale model radiation schemes are not informed with hydrometeor fall velocities. Moreover, fall velocity parametrizations in microphysics schemes do not necessarily fully overlap with the hydrometeor size and shape assumptions implemented in radiation schemes.

2. Noting that EMC² operates off-line, hydrometeor class fall velocities are typically reported in model outputs as weighted means. Because not all cloud schemes enable back-tracing of hydrometeor class fall velocities as a function of particle diameter using analytical expressions and weighted output fields (e.g., the convective cloud scheme in ModelE3; see Elsaesser et al., 2017), hydrometeor class fall velocities per subcolumn bin cannot be straightforwardly reproduced.

3. The total radar moments include combinations of the different hydrometeor class mixing ratios and, hence, cannot be determined using a single set of bulk LUTs per hydrometeor class.

Thus, EMC² calculates only the $Z_{{\mathrm{e}}_{\mathrm{hyd}}}$, $Z_{{\mathrm{e}}_{\mathrm{tot}}}$, and $Z_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}$. Equations (6b) and (7) are used to calculate $Z_{{\mathrm{e}}_{\mathrm{hyd}}}$ with bulk scattering LUTs for λ_r based on full single-particle Mie scattering calculations in the case of liquid hydrometeors and single-particle scattering calculations for the Yang et al. (2013) eight-column aggregate at a temperature of 270 K (see Ding et al., 2017) in the case of solid hydrometeors. Similar to the implementation of the radiation approach in the forward lidar calculations (Sect. 2.3.2), $r_{{\mathrm{e}}_{\mathrm{hyd}}}$ considers the fluffiness factor in the case of solid (ice) hydrometeors. Finally, $Z_{{\mathrm{e}}_{\mathrm{tot}}}$ and $Z_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}$ are calculated similar to the microphysics approach using Eqs. (8) and (9) with ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ (and ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{tot}}}$) calculated using Eq. (6a).

2.5 Hydrometeor classifications

Once the total lidar and/or radar variables are calculated, EMC² can be used to classify the subcolumn simulator output. Classification masks can serve as tools for direct comparisons between the simulator output and observational data utilizing similar classification methodologies, some of which can be used to calculate water phase ratios. Currently, EMC² incorporates three hydrometeor classification methods: the radar-sounding cloud and precipitation detection and classification, the modified fixed lidar variable threshold phase classification, and the COSP lidar simulator emulator (henceforth referred to as the COSP emulator; “emulator” is used here in its generic sense rather than a machine-learning context, and there is no training involved).

The radar-sounding cloud and precipitation detection method (Silber et al., 2021a; see also Vassel et al., 2019) emulates the combined use of relative humidity with respect to water (RH) sounding measurements for the detection of liquid-bearing clouds and radar echoes (received signal above the radar noise floor) for the detection of precipitating hydrometeors. In the case of large-scale model output, RH below 100 % in a model grid cell does not necessarily indicate a lack of cloud water, because of implemented assumptions concerning the sub-grid distribution of cloud water content (e.g., Smith, 1990). Therefore, the approach most consistent with the observational method is to simply use the cloud water-bearing subcolumn bins (see Sect. 2.2) to classify the subcolumn bin as “cloud”. “Precipitation” are those subcolumns bins in which $Z_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}\ge Z_{{\mathrm{e}}_{min}}$. Subcolumn bins that can be classified as both “cloud” and “precipitation” are set as “mixed”. We note that at temperatures below 0 ^∘C, the “mixed” classification type becomes more likely to represent a mixed-phase cloud with decreasing temperatures but, in general, may represent bins containing only liquid hydrometeors.

The modified fixed lidar variable threshold phase classification method is similar to previous studies that incorporated fixed LDR and ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$ thresholds to classify hydrometeor-bearing air volumes to “liquid” and “ice-only” using lidar measurements (e.g., Shupe, 2007; Thorsen and Fu, 2015). By default, however, EMC² includes two additional “undefined” classes that cover intermediate regions in the LDR-${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$, such that a subcolumn bin classified as “undefined1” has a higher probability that it includes some amount of liquid water while “undefined2” is more likely to only contain ice hydrometers (see Sect. 3.3 for discussion and illustration of the default thresholds). The notion behind the addition of these two “undefined” classes is the fixed-threshold method limitations that could originate in the following:

1. Drizzle- or rain-bearing air volumes may produce moderate ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$ and LDR on the order of 0.1, especially when ice hydrometeors are present in the same air volumes (e.g., Derr et al., 1976; Sassen, 2003; Silber et al., 2019a).

2. Cases can occur in which an entire relatively tenuous liquid-bearing cloud layer or other cases with liquid-bearing air volumes just above cloud base occur concurrently with ice precipitation sufficiently concentrated and intense to generate a combination of high ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$ and LDR that can reach values of 0.10–0.20 (e.g., Derr et al., 1976; Silber et al., 2018b, Fig. S1). This influence of ice hydrometeors also applies to rain-bearing air volumes.

3. Horizontally oriented ice hydrometeors may produce low (moderate) LDR (${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$) via specular reflection even in cases when the lidar is titled up to several degrees off-zenith (e.g., Noel et al., 2002; Silber et al., 2018b, Appendix A). Note that this limitation only applies to observations, since specular reflection and ice particle canting angles are not represented in large-scale models.

Adaptive fixed thresholds that vary with site, instrument, and period of study (e.g., Silber et al., 2018a; Zaremba et al., 2020) or lidar ratio constraints (e.g., Thorsen and Fu, 2015) can compensate for some of these limitations. However, these approaches cannot be objectively translated to the model output domain to enable a direct comparison between the observations and the simulator output. Therefore, the modified fixed threshold routine, which largely agrees with existing measurements yet acknowledges both model and observational uncertainties, may allow better direct comparisons to be made.

The emulator of the COSP lidar simulator follows the same equations and logic of the on-line lidar simulator (Cesana and Chepfer, 2013) implemented in numerous climate models. In short, the attenuated total backscatter (ATB) calculated in the COSP emulator routine while assuming η=0.7 is used to calculate the lidar scattering ratio (the ratio of total to molecular attenuated backscatter) for the detection of hydrometeors in subcolumns by selecting scattering ratio values larger than 5. Calculated cross-polar ATB as a function of the total ATB is then used to classify the detected hydrometeors into liquid or ice, based on an empirical phase discrimination line. As the last step of this classification method, hydrometeors below (top-down lidar view) or above (bottom-up lidar view) a subcolumn bin with scattering ratio larger than 30 are classified as “undefined”. Note that, unlike the on-line COSP lidar simulator, this emulator operates using the model vertical levels and does not interpolate the model output onto an evenly spaced vertical grid. As with all other EMC² forward calculations and classification routines, this emulator can operate using a top-down viewing approach thus providing a bridge between the COSP simulator and EMC².

Finally, EMC² also includes internal functions to calculate mass or frequency phase ratios using each of these hydrometeor classification methods, providing metrics to compare model output with observations or with outputs from other models.

3 Case study example: highly supercooled Antarctic cloud

To demonstrate the application of EMC² and its output using the different forward calculation approaches, here we describe and analyze a Lagrangian LES case study (Silber et al., 2019a) adjusted for running and testing the ModelE3 climate model (as well as other climate models) in SCM mode.

3.1 Case description

As described by Silber et al. (2019a), the stratiform cloud event that we compare with model simulations was observed over McMurdo Station, Antarctica, as part of the ARM West Antarctic Radiation Experiment (AWARE; Lubin et al., 2020), on 16 August 2016. During the event, cyclone-driven wind confluence with southwesterly katabatic flow resulted in relatively warm and moist marine air convergence along the Ross Ice Shelf coast, part of which was advected towards the McMurdo Station measurement site (Fig. 2a). This air convergence induced a widespread (>1000 km) cloud field evident by Clouds and the Earth’s Radiant Energy System (CERES) measurements indicating an extensive region with enhanced top of atmosphere (TOA) upwelling longwave radiation (dashed green shape in Fig. 2b; note that a surface-based temperature inversion is common during the austral winter resulting in smaller TOA radiation fluxes).

Figure 2 (Top) Sea level pressure (contours; land areas are masked), total condensate water path (logarithmic color scale), and surface winds (quivers) resolved by ERA5 (Hersbach et al., 2020) for 16 August 2016, at 10:00–11:00 UTC. The magenta-filled marker denotes McMurdo Station. The dashed green shape highlights a widespread cloud field along the Ross Ice Shelf (yellow shape) coast associated with the ground-based observations. The inset panel shows a topographic map of Antarctica (the red box highlights the region depicted in the main panel). (Bottom) Top of atmosphere upwelling longwave radiation measured by CERES Aqua on 16 August 2016, at 10:44–10:50 UTC. The 1 arcmin topographic data were developed by the National Geophysical Data Center (NGDC Amante, 2009) and are freely available at https://www.ngdc.noaa.gov/mgg/global/global.html (last access: 1 December 2021). Reproduced from Silber et al. (2019a).

Over McMurdo Station, a decoupled persistent mixed-phase cloud with temperatures as low as −29 ^∘C was observed for ∼ 39 h. The observed cloud was nearly continuously precipitating ice particles and was also drizzling for more than 7 h, concluded from a comprehensive analysis of sounding, HSRL, and KAZR measurements (see Silber et al., 2019a).

3.2 ModelE3 SCM configuration

Based on Lagrangian simulations constrained by the remote-sensing observations, Silber et al. (2019a) postulated that the activated ice-nucleating particle (INP) and cloud condensation nuclei (CCN) concentrations during the event were on the order of 0.2 L⁻¹ and 20 cm⁻³, respectively, to enable drizzle to be produced and precipitate along with ice precipitation below the highly supercooled cloud base. Their simulations were 9 h in duration, starting on 16 August 2016, at 01:00 UTC and ending at 10:00 UTC. By imposing large-scale vertical wind extracted from back-trajectory calculations, they emulated the transport of the cloud layer, initially forming in a stable layer, towards McMurdo Station. The end of the simulation at 10:00 UTC designated the time at which the cloud field reached the fixed observational site at McMurdo, and hence, statistics of that hour of observations (10:00–11:00 UTC) were compared with the model output.

Here, we slightly adjusted the case study initialization files to enable running this case in a climate model's SCM mode while using the same single hour of radar and lidar observations and the Distributed Hydrodynamic Aerosol and Radiative Modeling Application (DHARMA) model (Stevens et al., 2002) baseline LES output (see Silber et al., 2019a) as benchmarks. Namely, we simplified the profiles of vertical motion (as in Silber et al., 2020) and then converted the height coordinates of the initial sounding (see Silber et al., 2019a) and vertical wind time series to pressure coordinates (these converted sounding and forcing files are available at https://doi.org/10.17632/gz4gdn3jvz.1, Silber, 2021). The utilization of these files to run the highly supercooled cloud case study enables testing of an SCM and, hence, the cloud schemes implemented in a climate model.

The LES is initialized here with activated INP concentration of 0.1 L⁻¹ and with both cloud ice and snow hydrometeor classes (whereas Silber et al. (2019a) used only a single ice class), and the simplified vertical motion profiles are then imposed. These few LES adjustments make the model configuration more consistent with climate model microphysics and SCM initialization while effectively resulting in the same hydrometeor content and cloud evolution as the baseline simulation presented in Silber et al. (2019a) (not shown). Thus, the SCM is run equivalently to the LES, that is, using the same initial sounding and forcing files, setting the coordinates to −77.85^∘ S, 166.72^∘ E, and prescribing a monomodal log-normal aerosol particle concentration of 20 cm⁻³ with a mean radius of 0.1 µm, geometrical standard deviation of 2, and a hygroscopicity parameter of 0.4. We note that activated INP concentrations are not prescribed in the SCM simulations because (a) ModelE3's final configurations are defined by specific values of certain model tuning parameters (among others) associated with INP parameterization, and (b) diagnostically prescribing the activated INP concentration is not faithful to the temperature-dependent approach implemented in the model and, hence, would not necessarily be informative of true climate model weaknesses.

3.3 Comparison between observations and ModelE3 SCM using EMC²

The following figures show the EMC² forward calculation results using the DHARMA LES simulation and ModelE3's configuration Tun3, one of four configurations of ModelE3 derived in part via a machine learning parameter tuning approach that will be described in a manuscript in preparation, to be included in the Coupled Model Intercomparison Project Phase 6 (CMIP6). The SCM using this configuration was able to generate a cloud-top inversion and turbulent layer via cloud-top radiative cooling and produced the best agreement with the observations and the LES relative to the other three configurations (see Appendix B). Whereas here we examine application of EMC² to a single ModelE3 configuration in SCM mode in a case where we can also compare with LES, we note that EMC² is designed to enable detailed evaluation of atmospheric thermodynamic profile and cloud properties extracted from global simulations of ModelE3 configurations and other climate models against long-term datasets at fixed sites in future dedicated work.

The left panels in Fig. 3 show the mixing ratios of the four hydrometeor classes evolving through the simulated SCM case study. These depicted mixing ratios are the output of the SCM's stratiform cloud scheme, but because this simulation does not generate any convective hydrometeors (as expected), they also represent the total mixing ratios. Cloud water mass (top panel) dominates over the other hydrometeor classes in hydrometeor-bearing model grid cells through much of the simulation, even at lower levels in which the cloud water fraction is rather low (Fig. 3, right panels). Rain (effectively drizzle) is produced by the model as well but has a relatively smaller mass compared with that of snow generated. This reduced amount of rainwater in the SCM simulation relative to rainwater dominance in the simulations of Silber et al. (2019a) is largely the result of the different autoconversion parameterization schemes implemented in ModelE3 (Seifert and Beheng, 2001) relative to that implemented by default in the DHARMA LES (Khairoutdinov and Kogan, 2000), which produces significantly smaller rain mass mixing ratios in this case (not shown), contrary to some previous studies (e.g., Heiblum et al., 2016; Xiao et al., 2021). Understanding the source of this differing autoconversion parameterization behavior requires a dedicated study that is a beyond the scope of this article.

Figure 3 ModelE3 configuration Tun3 single-column model (SCM) output for the Antarctic cloud case study (16 August 2016) showing mixing ratio (left) and hydrometeor fraction (right) time–height curtain plots of (from top to bottom) cloud water, rain, cloud ice, and snow. The generalized hydrometeor fraction used in the radiation approach (see Sect. 2) is depicted in the bottom right panel.

Figures 4 and 5 depict the HSRL and KAZR variables observed during a single hour over McMurdo Station and simulated with EMC² using the DHARMA LES three-dimensional output at the end of that simulation (without using the subcolumn generator), and by applying each of the three approaches on ModelE3 configuration Tun3 SCM output for 05:00 UTC. Because our goal in this section is to demonstrate that EMC² can reasonably match cloud observations given comparable input, we present the EMC²-processed SCM output 4 h into the simulation when cloud top heights are similar to observed instead of the end of the simulation (the SCM develops the supercooled cloud faster than the baseline LES; see Appendix B). When evaluating the processed model output against the observations, we essentially exchange temporal resolution with spatial resolution (three-dimensional model domain in the case of the LES) or an emulated spatial resolution (in the case of the SCM). We set the number of SCM subcolumns (N_s) to 100, which effectively represents hydrometeor class fractions up to the second decimal point and enables drawing more robust statistics by emulating sub-grid variability of all hydrometeor classes combined.

Figure 4 Observed and simulated HSRL particulate backscatter cross-section (${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$; top row), particulate extinction cross-section (${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{tot}}}$; middle row), and linear depolarization ratio (LDR; bottom row) using the DHARMA LES output and ModelE3 configuration Tun3 SCM output. The columns show (from left to right) the observations from McMurdo Station between 10:00–11:00 UTC, DHARMA LES three-dimensional output (horizontal dimensions stacked onto a single domain column dimension) processed using EMC² microphysics with the subcolumn generator turned off (N_s=1), ModelE3 EMC² output using the microphysics approach, and EMC² output using the radiation approach. The DHARMA LES output corresponds to 10:00 UTC (arrival of cloud field at McMurdo Station; see Sect. 3.2), while ModelE3's EMC² panels depict the SCM output for 05:00 UTC (see text) processed with EMC² using 100 subcolumns. A mask denoting full lidar signal attenuation generated using a total accumulated optical thickness (τ_tot) condition of τ_tot>4 is plotted over the simulated data (hatched areas).

Figure 5 As in Fig. 4, but with the KAZR attenuated equivalent reflectivity factor ($Z_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}$; top row), mean Doppler velocity ($V_{D_{\mathrm{tot}}}$; middle row), and spectrum width (${\mathit{\sigma }}_{D_{\mathrm{tot}}}$; bottom row). Note that only $Z_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}$ is calculated in EMC²'s radiation approach (see Sect. 2.4.2). A radar signal-to-noise ratio (SNR) mask of $Z_{{\mathrm{e}}_{min}}>Z_{{\mathrm{e}}_{\mathrm{tot}}}$ ($Z_{{\mathrm{e}}_{min}}$ is the minimum detectable signal on each model level) is applied to the plotted data (hatched areas).

The processed LES output exhibits good apparent agreement with the observations, evident by the comparable lidar and radar variable values and their horizontal variability (Figs. 4 and 5), the vertical cloud structure and boundaries, as well as the full lidar signal attenuation near cloud top (Fig. 4). A multi-layer cloud water structure developed by the LES is suggested by the intermittent breaks in the large ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$ and ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{tot}}}$ values, supported by a separate analysis of ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ and ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ (not shown). This multi-layer structure is also indicated by the lidar observations (Fig. 4) and was comprehensively discussed by Silber et al. (2019a).

Using the microphysics approach, the SCM sub-grid variability is more pronounced relative to the radiation approach owing to the implementation of cloud water sub-grid variability (as defined in the MG2 microphysics scheme), which is not considered in the radiation approach. Evaluation of heights with full cloud cover, indicated by the large ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$ and ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{tot}}}$ values (see also Fig. 3, right), suggests that the SCM has an apparent reasonable agreement with the observations there. ModelE3's macrophysics scheme leads to notably ragged liquid cloud base heights that are evident in the SCM fields compared to the more uniform liquid cloud base in both observations and LES, as discussed further below. Relatively enhanced ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{tot}}}$ values are produced in some SCM subcolumns, which leads to full lidar signal attenuation, as also indicated by the HSRL measurements (Fig. 4). Full lidar signal attenuation does not occur in this case using the radiation approach because of the uniform distribution of cloud water between subcolumns commensurate with τ_tot values just below the full attenuation threshold of 4 (not shown). The subcolumns with full lidar signal attenuation in the microphysics approach call for the use of radar measurements for cloud-top detection (e.g., Fig. 5 discussed below).

Except for the overestimated LDR values below the supercooled cloud layer relative to the observations, the values of the other two lidar variables and the general scenario structure appear to agree with the observations and LES in both the microphysics and radiation approaches, which exhibit an encouraging consistency with each other (Fig. 4). The comparison of the subcolumn-averaged lidar-variable profiles illustrated in Fig. 6 allows a more quantitative comparison. This figure better indicates that the DHARMA LES ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{tot}}}$ within the main liquid-bearing cloud layer is within the range of the observations and their uncertainty, consistent with the conclusions of Silber et al. (2019a) that activated CCN number concentrations were low during this highly supercooled drizzle event.

Figure 6 Mean profiles of observed (time-averaged) and simulated (subcolumn-averaged) HSRL and KAZR variables calculated using EMC²'s microphysics (top) and radiation (bottom) approaches as shown in Figs. 4 and 5 (full lidar signal attenuation and radar SNR masks are applied prior to averaging). Shaded regions denote the mean ±1 standard deviation.

The EMC² lidar output also highlights a few model weaknesses. For example, using the microphysics approach, cloud base (at a height of ∼ 1.0 km) is highly variable and may extend nearly down to the surface, which is in contradiction to the observations, where cloud base variability is on the order of a few hundred meters throughout the depicted period (Fig. 4). Using the radiation approach, on the other hand, breaks in hydrometeor cover are seen below the fully overcast layer, as a result of the generalized hydrometeor fraction used in the radiative transfer calculations, which can be lower or higher than the associated hydrometeor class fraction (e.g., Fig. 3, right). Lesser (greater) generalized hydrometeor fractions relative to the associated hydrometeor class fraction therefore imply greater (smaller) subcolumn bin mixing ratios (see Eq. 1) and, hence, enhanced (diminished) associated ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ and ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ values. Thus, while the occurrence of low-level cloud water-bearing bins produces full attenuation of the simulated lidar signal in some subcolumns using the microphysics approach (Fig. 4), the smaller ${\mathit{\alpha }}_{{\mathrm{p}}_{\mathrm{hyd}}}$ ensuing from the larger low-level generalized cloud fraction relative to the cloud water fraction (see Fig. 3, right) causes no low-level (or any level in this case) signal extinction when the radiation approach is used (see Fig. 4). Note that in both the microphysics and radiation approaches the subcolumn representation of hydrometeor mass remains consistent with the model output variables; i.e., Eq. (1) holds for each hydrometeor class.

The microphysics and radiation approaches exhibit $Z_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}$ values that are too large, especially at higher levels (Figs. 5 and 6). EMC²'s radar processing using the microphysics approach provides the $V_{D_{\mathrm{tot}}}$ and ${\mathit{\sigma }}_{D_{\mathrm{tot}}}$ variables in addition to the $Z_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}$ variable. We note that because spectral broadening terms other than the microphysical broadening are currently not considered in EMC², the calculated ${\mathit{\sigma }}_{D_{\mathrm{tot}}}$ generally serves as a lower bound for comparison purposes; i.e., the simulated ${\mathit{\sigma }}_{D_{\mathrm{tot}}}$ values need to be smaller to some extent than the observed values. As indicated from Fig. 5, both $V_{D_{\mathrm{tot}}}$ and ${\mathit{\sigma }}_{D_{\mathrm{tot}}}$ show grossly reasonable correspondence between the observations and the SCM, except at heights above ∼ 0.8 km, where $V_{D_{\mathrm{tot}}}$ values show large deviations (Fig. 6). A separate analysis (not shown) suggests that, in this case, these deviations are mainly the result of relatively fast fall velocities (see Table 1) and the dominance of large snow hydrometeors over $Z_{{\mathrm{e}}_{\mathrm{tot},\mathrm{att}}}$ at these levels and not the product of vertical air motion being convolved into the radar moments only in the observations.

Figure 7 delineates the three hydrometeor classification methods currently implemented in EMC² applied over the processed model output, with the COSP emulator “observing” the model domain from the top down, similar to the on-line simulator. Congruent with the description above of the calculated radar and lidar variables, the radar-sounding and modified fixed lidar threshold methods show the domination of liquid hydrometeor classes above ∼ 1.0 km, the prevalence of precipitating hydrometeors at lower levels (with lower occurrence using the radiation scheme approach), and the occasional full attenuation of the simulated lidar signal in the case of the lidar-based classification method. The COSP emulator detects a clear liquid-bearing subcolumn bin signal at cloud tops. However, because the cloud-top layer is highly reflective, generating large lidar scattering ratio values (the ratio of total to molecular attenuated backscatter), most of the underlying layers are either classified as “undefined” or generate signals too weak to be detected.

Figure 7 Hydrometeor classification of ModelE3 configuration Tun3 SCM output for 05:00 UTC processed using EMC²'s microphysics (left) and radiation (middle) approaches while considering signal attenuation and detectability: (from top to bottom) radar-sounding method, modified fixed lidar variable threshold method, and the COSP lidar simulator emulator (top-down view). (Right) Brief summary of each classification method.

Using each of the three classification methods, phase ratio statistics can be generated with EMC² offering a method for analyzing the SCM simulation. Figure 8 portrays the temporal evolution of the SCM simulation from the view of the simulated instruments and classification methods using the microphysics or radiation approaches. When radar-sounding and fixed lidar threshold methods are applied while using the microphysics approach, the evolution of the simulated cloud appears self-consistent between the two methods and generally follows the prototypical appearance of nearly continuously precipitating liquid-bearing cloud layers with weakly varying cloud base height (e.g., de Boer et al., 2011; Fridlind and Ackerman, 2018; Silber et al., 2021a). Here, the radar-sounding skill is associated with the diminished cloud water fraction relative to the large fraction of the other hydrometeor classes (Fig. 3, right) and the method's ability to correctly detect cloud water layers. The capabilities of this method should nonetheless be considered carefully when directly compared with observations because realistic sounding profiles typically lack the fine temporal resolution emulated by EMC² here, and while in situ observations can provide a robust characterization of liquid-bearing cloud layers, they can also produce sporadic false positive or negative cloud detections (e.g., Silber et al., 2020, Fig. S1; Vassel et al., 2019). In the case of the modified fixed lidar threshold method, the low-level phase ratio skill originates from the hydrometeor-bearing subcolumn bins being largely classified as “ice”. This classification decision is the result of the low ${\mathit{\beta }}_{{\mathrm{p}}_{\mathrm{tot}}}$ and moderate-to-high LDR (e.g., Fig. 6) produced by the prevalence of ice hydrometeors relative to cloud water (Fig. 3, right) and the greater mass (and likely volume due to the spherical representation) of these hydrometeors relative to rain (Fig. 3, left).

Figure 8 Time–height curtain plots showing the liquid-bearing frequency phase ratio (calculated over all hydrometeor-bearing subcolumn bins) of ModelE3 configuration Tun3 SCM simulation using the microphysics (left) and radiation (right) approaches. The classification methods used to calculate the phase ratios are (from top to bottom) the radar-sounding method, the modified fixed lidar variable threshold method, the COSP lidar simulator emulator (top-down view), and the on-line COSP lidar simulator implemented in ModelE3, which is only processed using ModelE3's radiation approach. In this figure, the “mixed” class of the radar-sounding method is counted as liquid, while the “undefined” classes in the other two methods are treated as “non-liquid” even though in some cases they are more likely to be liquid-bearing (e.g., the “undefined1” class in the modified fixed lidar variable threshold method). The lower-left panel shows the mass phase ratio calculated using the ModelE3 output without EMC² processing.

The COSP emulator using the microphysics approach with a top-down view is consistent with the example in Fig. 7, in which hydrometeor detection is limited by the optically thick and highly reflective cloud-top region, resembling observational retrievals (e.g., Cesana and Chepfer, 2013; Cesana et al., 2016).

Using the radiation approach, phase ratios more frequently show sharper transitions between the extreme values (all liquid or all non-liquid) stemming from the utilization of the generalized hydrometeor fraction. In the case of the radar-sounding method, for example, there is a distinct dominance of liquid-bearing bins, which only require any amount of cloud water to be classified as such. This dominance originates in the common occurrence of some cloud water mass mixing ratios in model grid cells throughout the SCM simulation (Fig. 3, left) combined with the implementation of the generalized hydrometeor fraction, which necessarily increases the overlap between cloud water and other hydrometeor classes. The limited number of model grid cells in which the subcolumn bins exhibit a more mixed behavior are the result of the randomized component of the subcolumn generator, which does not necessarily require overlap between cloud and precipitating hydrometeors (see Sect. 2.2). Based on these classification results, we suggest that the radar-sounding method could lead to hydrometeor classification favoring liquid-bearing classes when the radiation approach or similar model output with a generalized hydrometeor fraction is used.

The modified fixed lidar threshold method is relatively consistent with the microphysics approach (Fig. 8), even though some times (mainly around 06:00 and 10:00 UTC) are characterized by greater (smaller) relative liquid occurrence at lower (higher) model levels (phase ratio values closer to 0.5). These phase ratio differences relative to the microphysics approach are the result of the convolution of the generalized hydrometeor fraction and its deviation relative to the cloud water fraction (Fig. 3, right).

Using a top-down view, the COSP emulator agrees with both its application using the microphysics approach as well as the COSP output from the on-line simulator implemented in ModelE3, which utilizes the radiation approach (Fig. 8. Unlike the on-line simulator, the emulator detects some “undefined” hydrometeors at low levels (down to the surface), which can be explained by the lack of model interpolation onto a uniform vertical grid and/or small differences in the compiled simulator code related to signal attenuation (e.g., the accumulation of optical thickness).

We note that the treatment of the COSP emulator's “undefined” subcolumn bins as “ice” to produce phase ratio statistics leaves the impression that only ice hydrometeors exist below cloud top. However, a rather different impression of mostly liquid water dominance, though not as stark as in the radar-sounding method using the radiation approach, is perceived when the mass phase ratio calculated using the raw SCM output is examined (Fig. 8, lower left). Contrary to the COSP emulator, treating “undefined” bins as “ice” in the modified fixed lidar threshold method increases its apparent frequency phase ratio agreement with the mass phase ratio in multiple time–height bins. Phase classification depends on instrument measurement characteristics and limitations and hydrometeor properties such as their class, relative mixture with other hydrometeor classes, as well as their size distributions. Therefore, such an apparent agreement between different variables and phase occurrence metrics, as well as between the same variables and metrics based on different instruments and/or methodologies, should be taken with a grain of salt (cf. Cesana et al., 2021; see also Silber et al., 2021c).

Figure 8 demonstrates the sensitivity of phase ratio statistics to the classification method, the viewing direction of the examined instrument, and the method by which “liquid” and “non-liquid” or “ice” classes are being counted. It shows that the use of forward simulators alone is not a guarantee for an “apples-to-apples” comparison, which requires matching processing steps to ensure its robustness.

4 Summary

EMC² provides an easy to use and flexible framework for the analysis of large-scale model output and its direct comparison with ground-based observations via the generation of subcolumns intended to explicitly represent a sub-grid scale variability, and the simulation of ground-based (and air- or space-borne) radars and lidars. EMC²'s framework is already tailored to the MG2 two-moment microphysics while using single-particle scattering LUTs and has the proper infrastructure required for it to be customized to other similar schemes, as well as high-resolution model output. EMC²'s option for using radiation scheme logic in the subcolumn generator and simulator enables direct comparison with other on-line active instrument simulators (e.g., the COSP lidar simulator) with a bottom-up or top-down view option that can bridge between different methodologies by evaluating differences between the outputs resulting from their implementation.

Because it is generally designed to emulate ground-based systems, EMC² is suitable for the evaluation of column output extracted from global simulations against long-term ground-based datasets. The general adaptability of the software code to other climate models and instruments via the model and instrument Python classes renders EMC² as a flexible framework to enable consistent and reproducible post-processing methods and evaluation across multiple models.

An AWARE case study was used to illustrate the application of EMC² to LES and SCM simulations of a highly supercooled Antarctic cloud, including the utility of the program for hydrometeor classification using radar-sounding, lidar variable thresholds, and COSP emulator methods. The ModelE3 SCM using configuration Tun3 showed general agreement with the observations at the examined simulation time as well as with the baseline DHARMA LES used to develop this case study (see Silber et al., 2019a). The LES output can be processed with EMC² after a few adaptations made only to an inherited Model class (see Sect. 2.1). Thus, although it was developed for large-scale models, EMC² can also be used to compare cloud resolving or LES models with observations (as shown for DHARMA). EMC² also allows the implementation of advanced scattering model calculations in the forward calculations via customized LUTs that could be matched to some scattering assumptions made by models, for example, the implementation of the MODIS C6 calculations in both ModelE3's radiation scheme and EMC².

The AWARE case study presented here is suitable for simulation by any global model in SCM mode (see input file repository specified under code and data availability). Case study observations, as well as ModelE3 SCM and DHARMA LES inputs and outputs from EMC² used to produce all examples above, are also provided for step-by-step illustration (see code and data availability). We plan that additional case study examples will similarly be provided to illustrate results under differing cloud regimes.

Planned future additions to EMC² include an extension to ground-based scanning radars, a Mie scattering calculator, spectral broadening estimates for the radar simulator, and a multiple-scattering model for the lidar simulator, all of which will be configured for consistency with model physics and output fields. We invite the community to take advantage of the framework provided by EMC² and to contribute to its further development and applications.

Appendix A: Lists of acronyms, abbreviations, and symbols

Gettelman and Morrison (2015)

Table A1 Acronyms and subscript abbreviations used in this article.

Table A2 Symbols of variables and parameters used in this article and their units (unless explicitly stated otherwise in the text).

Appendix B: ModelE3 SCM output using the four different model configurations

Figure B1 Time series showing (a) liquid water path (LWP) and (b) ice water path (IWP) in each case study simulation using the DHARMA LES (baseline simulation in Silber et al., 2019a), and the four ModelE3 configurations (Tun1, Tun2, Tun3, and Phys). The Eulerian retrievals shown for 10:00 UTC (error bars denote uncertainty) correspond to the time at which the Lagrangian simulated domain approaches McMurdo Station (see Sect. 3.2). (c, d, e) Temperature, turbulent kinetic energy, and cloud water mixing ratio profiles at the end of the SCM and DHARMA LES simulations.

Figure B1 illustrates time series of liquid water path (LWP) and ice water path (IWP) from the ModelE3 SCM case study output using configurations Tun1, Tun2, Tun3, and Phys (see Sect. 3.3), as well as the temperature, turbulent kinetic energy (TKE), and cloud water mixing ratio profiles at the end of the simulation. Out of these four model configurations, only configurations Tun3 and Phys maintain substantial amounts of LWP in this highly supercooled cloud case study and agree with both the depicted DHARMA LES output and the observed LWP retrievals (see Silber et al., 2019a). Moreover, the SCM using each of these two configurations is able to develop a cloud-top inversion and TKE, both of which are driven by radiative cooling of cloud water, consistent with the DHARMA LES (see Fig. B1, bottom panels) and various polar cloud observations (e.g., Morrison et al., 2012; Silber et al., 2019a). All of ModelE3's configurations except for Tun1, which generates the largest amount of ice, are within the range of IWP estimated retrieval uncertainty. Out of those three model configurations, Tun3 is closest to the retrieved IWP and best matches the cloud formation and evolution in the LES. Because the supercooled cloud is developed faster in the SCM than in the LES model (see Fig. B1, top and middle), we choose to demonstrate EMC² using the SCM output corresponding to 05:00 UTC. Figure B2 shows the mean lidar and radar variable profiles at that model time step from EMC² for each of the four model configurations together with the time-averaged observed profiles. Using either the microphysics or radiation approach, configurations Tun3 and Phys best match the observations. Most of the radar and lidar variables calculated using these two configurations are largely consistent with each other, but overall, the Tun3 SCM output shows the best agreement with the mean observed profiles.

To summarize, two of ModelE3's final four configurations, that is, configurations Tun3 and Phys, show reasonable agreement with both the observed quantities as well as the LES output variables. Because configuration Tun3 performs slightly better, we focus on this model configuration for detailed comparison with LES and observations.

Figure B2 Mean profiles of observed (time-averaged) and simulated (subcolumn-averaged) HSRL and KAZR variables calculated using EMC² microphysics (top) and radiation (bottom) approaches using the four ModelE3 configurations (see legend), without applying a lidar and radar extinction or SNR masks, respectively.