The Cloud-resolving model Radar SIMulator (CR-SIM) Version 3.3: description and applications of a virtual observatory

Ground-based observatories use multisensor observations to characterize cloud and precipitation properties. One of the challenges is how to design strategies to best use these observations to understand these properties and evaluate weather and climate models. This paper introduces the Cloud-resolving model Radar SIMulator (CR-SIM), which uses output from high-resolution cloud-resolving models (CRMs) to emulate multiwavelength, zenith-pointing, and scanning radar observables and multisensor (radar and lidar) products. CR-SIM allows for direct comparison between an atmospheric model simulation and remote-sensing products using a forward-modeling framework consistent with the microphysical assumptions used in the atmospheric model. CR-SIM has the flexibility to easily incorporate additional microphysical modules, such as microphysical schemes and scattering calculations, and expand the applications to simulate multisensor retrieval products. In this paper, we present several applications of CR-SIM for evaluating the representativeness of cloud microphysics and dynamics in a CRM, quantifying uncertainties in radar–lidar integrated cloud products and multi-Doppler wind retrievals, and optimizing radar sampling strategy using observing system simulation experiments. These applications demonstrate CR-SIM as a virtual observatory operator on high-resolution model output for a consistent comparison between model results and observations to aid interpretation of the differences and improve understanding of the representativeness errors due to the sampling limitations of the ground-based measurements. CR-SIM is licensed under the GNU GPL package and both the software and the user guide are publicly available to the scientific community.

1 Introduction

Ground-based observatories offer an integrated view of cloud and precipitation systems complementary to that available from satellites with excellent vertical resolution, especially in the boundary layer, and an accompanying description of the large-scale forcing. Currently, a number of observatories are continuously operated in different climate regimes (Illingworth et al., 2007; Löhnert et al., 2015; Stevens et al., 2016; Mather et al., 2016) with evolving measurement capabilities. In the beginning, zenith-pointing cloud radars, lidars, and radiometers provided the primary cloud and precipitation measurements. Recently, the need to characterize the mesoscale organization of clouds and precipitation over a larger domain has heightened the sophistication and complexity of these observatories to go beyond single, one-dimensional profiling measurements. For example, the US Department of Energy (DOE) Atmospheric Radiation Measurement (ARM) observatories offer observations from distributed networks of profiling and scanning radars, lidars, and radiometers (Turner and Ellingson, 2016; North et al., 2017).

Multiparametric information from profiling and scanning radars, lidars, and radiometers has been used to retrieve cloud microphysical and kinematic properties, such as ice particle properties (e.g., Zhang et al., 2014; Kneifel et al., 2015, 2016; Matrosov et al., 2017; Von Lerber et al., 2017), and vertical velocities (North et al., 2017; Giangrande et al., 2016). However, the comparison between the retrieved observables (e.g., ice water content from radar reflectivity) and model-produced parameters often involves large uncertainties. Several factors, not limited to the nature of ground-based observations, challenge model evaluation using these observations. In many cases, the retrieval algorithms are based on statistical estimation of ill-posed inverse problems, and the results may not capture well the observed variability of natural data because of limitations from assumptions embedded in the retrieval algorithms (e.g., Szyrmer et al., 2012; Szyrmer and Zawadzki, 2014). Furthermore, determining critical parameters for model evaluation such as the cloud fraction profile requires complementary, synergistic observations from both radar and lidar. One such example is the Active Remote Sensing of Clouds (ARSCL; Clothiaux et al., 2001) product that combines radar and lidar data to determine the hydrometeor height distributions. Other examples of critical parameters that require a multisensor approach include cloud and precipitation classification schemes (Illingworth et al., 2007) and hydrometeor phase classification (e.g., Shupe, 2007; Luke et al., 2010; Lamer et al., 2018). So, how do we best compare such products developed using multiple sensors with different capabilities (i.e., sensitivity) with numerical model output? Further, challenges may arise from the sampling strategy used to obtain the observations. For example, a recent study by Oue et al. (2016) has shown that zenith-pointing observations from one location are inadequate to provide reliable cloud fraction profile estimates in a cumulus field. A similar investigation on 3D wind retrievals in deep convection using multi-Doppler radar techniques highlights similar deficiencies of our current observing systems (Oue et al., 2019a). How do we best quantify the measurement uncertainty introduced by the observational strategy?

In this paper, we introduce the Cloud-resolving model Radar SIMulator (CR-SIM), which has been continuously developed over the last 5 years to facilitate model–observation comparisons. CR-SIM applies forward simulators to atmospheric model output to simulate ground-based measurements. These simulations may be used (1) to compare directly to the measurements, which provides an apples-to-apples comparison of the observed variables, or (2) as input to retrieval algorithms to assess the retrieval methodology or sampling strategy using the original atmospheric model output as “truth”. In this study, the CR-SIM architecture and capabilities are presented along with a series of forward simulations that emphasize its capabilities. In particular, we highlight the applications of CR-SIM in investigations of observational uncertainties. Although accurate estimation of uncertainties in the retrieval products (e.g., ice water content, liquid water content, vertical velocity) is challenging, forward simulators allow us to emulate the observational retrieval products accounting for known error sources to understand the exact impacts of those error sources on the products by comparisons with the truth, which is usually the input model data. Observing system simulation experiments (OSSEs) take advantage of forward modeling to produce simulated measurements. The understanding from OSSEs help us to (i) evaluate the model simulations using the observations while accounting for the observation limitations, (ii) estimate uncertainties in retrieval techniques used, (iii) propose new retrieval techniques accounting for the uncertainties, and (iv) optimize new observation system strategies. This study demonstrates the application of the CR-SIM forward simulator in several OSSEs in which ARM multisensor products, such as cloud locations and vertical velocity, are evaluated by considering limitations inherently imposed by the nature of the observations. (A list of acronyms is provided in Appendix B for easy reference.)

2 Forward simulators

Forward simulators have been widely used to design observing systems and to provide an alternative path to model–observation comparisons by transforming the model geophysical quantities into remote-sensing observables. There are several sophisticated radar simulators which have been developed for specialized applications of interest. For example, Snyder et al. (2017a, b) simulate polarimetric radar characteristics of a supercell using radar forward simulators to understand the contribution of microphysical characteristics to the polarimetric properties and their wavelength dependency. They account for the water fraction of solid ice particles to realistically simulate differential reflectivity (Z_DR) columns, specific differential phase (K_DP) columns, and co-polar correlation coefficient (ρ_hv) rings in supercells. A cloud radar simulator developed by Zhang et al. (2018) is designed to simulate vertically pointing cloud radar reflectivity (e.g., Ka- and W-band radars) using global climate model (GCM) output, which is beneficial for comparison of datasets at different scales (cloud-scale observational data versus global-scale data). Matsui et al. (2019) simulate polarimetric precipitation radar-based hydrometeor classification, vertical velocity, and rain rate from CRM output to examine uncertainties in the retrieval algorithms and model microphysical parameterizations using the POLArimetric Radar Retrieval and Instrument Simulator (POLARRIS). The uncertainties are attributed to assumptions of hydrometeor particle size distribution, density, axis ratio, and/or canting angle. Lamer et al. (2018) developed a GCM-oriented ground-observation forward simulator ((GO)²-SIM), which is a comprehensive radar–lidar simulator for the GCMs that emulates radar Doppler spectra moments, lidar backscatter, and depolarization, and it provides synthetic estimates of mixed-phase cloud occurrence in the GCMs that are comparable to those estimated from observations using the same methodology.

CR-SIM can simulate the ideal multiwavelength radar and lidar observables and multisensor integrated products. The zenith-pointing and scanning radar observables include radar reflectivity, Doppler velocity, spectrum width, and polarimetric fields. Zenith-pointing lidar observables include lidar backscatter and extinction coefficient. The idea behind CR-SIM is to have a forward-model operator that provides idealized radar and lidar observables (i.e., actual observations after perfect quality control and correction for the total (two-way) attenuation) on the same grid as in the CRMs or large-eddy simulations (LESs) to facilitate model–observation comparisons. Further, the design is flexible enough to be coupled with different microphysical schemes and different scattering methods (e.g., T-matrix method by Mishchenko, 2000; discrete dipole approximation by Yurkin and Hoekstra, 2011).

The CR-SIM forward simulator is tailored to compute radar and lidar observables by integrating scattering properties over the discrete particle size distributions (PSDs) using a constant bin size for each hydrometeor (Table 1), based on the microphysical scheme used in the CRM or LES. The environmental variables are obtained from a mandatory set of model output variables consisting of pressure, temperature, dry air density, and height above sea level. The single-scattering properties are calculated using the T-matrix method and packaged as look-up tables (LUTs) in CR-SIM. The simulated idealized radar and lidar variables are provided at each model grid box and can be easily compared with real observations.

Table 1 The minimum and maximum sizes and bin spacing of simulated particles for each hydrometeor category from the bulk microphysics schemes except the P3 scheme. “Particle size” here refers to the maximum particle dimension.

2.1 Scattering properties

The LUTs consist of the complex scattering amplitudes S_ij of the 2×2 scattering matrix for single, nonspherical particles at fixed orientations with equally spaced particle sizes computed using the T-matrix code of Mishchenko and Travis (1998) and Mishchenko (2000). Following Ryzhkov et al. (2011), we assume that the scattering characteristics of arbitrarily oriented particles can be expressed by the scattering amplitudes f_a and f_b corresponding to the principal axes of spheroid when the electric vector of the illuminating electromagnetic wave is directed along the spheroid axes a and b, respectively. Each hydrometeor species is characterized by its equivalent volume diameter, $D=(a{b}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{3}}$, where a is the symmetry axis of the spheroid, and b is its transverse axis. In this convection, an oblate hydrometeor has a<b. The scattering amplitudes f_a and f_b correspond to S_vv and −S_hh computed by the Mishchenko (2000) T-matrix code, respectively (f_a≡S_vv; $f_b\equiv -S_{\mathrm{hh}})$. The minus sign in the expression for f_b is to account for switching from the forward-scattering convention used for the amplitude and phase matrix components in the T-matrix code to the backscatter-alignment convention used in the definitions of polarimetric radar measurements except for forward-scattering parameters (attenuation and phase shift).

The LUTs of scattering properties for single particles are constructed for each hydrometeor class corresponding to the CRM or LES simulation output (e.g., rain drop, snowflakes, cloud droplet, ice crystal, graupel) as a function of particle phase, bulk density, aspect ratio, size, and temperature. We used fixed values of these parameters depending on the “scattering type”, which must be assigned to the corresponding hydrometeor class. The bulk density used is as parameterized in the selected microphysics scheme, assuming spheroid particle shapes. For each hydrometeor class, the complex scattering amplitudes are calculated for 91 elevation angles from 0 to 90^∘ with a spacing of 1^∘; five radar frequencies at 3 GHz (S band), 5.5 GHz (C band), 9.5 GHz (X band), 35 GHz (Ka band), and 94 GHz (W band); different temperature ranges for the liquid hydrometeors; different particle densities for solid hydrometeors; and a few different values of particle aspect ratio. The scattering types built in the current LUTs and their parameter settings are presented in Table 2. For lidar scattering properties, the single-particle extinction σ_α and backscattering cross section σ_β for spherical cloud droplets and cloud ice are calculated using the BHMIE Mie code (Bohren and Huffman, 1998). CR-SIM operates for observables from the ceilometer (wavelength of 905 nm) and micro pulse lidar (MPL, wavelengths of 353 and 532 nm).

Table 2 Scattering types built in the current LUTs and settings for the scattering calculations. Diameter ranges are the same as shown in Table 1.

Although most of the parameters related to hydrometeor particles (e.g., particle bulk density, size) required in the scattering calculations can be obtained from the prognostic and diagnostic variables from the CRM or LES, aspect ratios and canting angles must be assumed in the simulator and as such are prescribed by the users. All liquid and ice hydrometeors are assumed to be oblate spheroids, except cloud droplets. Raindrops are represented as oblate spheroids with a size-dependent aspect ratio following an empirical equation as a function of particle diameter, based either on Brandes et al. (2002) or Andsager et al. (1999). A fixed aspect ratio is used for each solid hydrometeor category, but for graupel and hail the empirical expression proposed by Ryzhkov et al. (2011) is also available. For all model hydrometeors, the radar polarimetric variables (which depend on particle orientation) are calculated using complex scattering amplitudes from the pre-built LUTs, assuming a mean particle canting angle of 0^∘ (Ryzhkov, 2001) with a choice of the particle orientation distribution. The possible choices for the distribution of particle orientations are fully (three-dimensional) random orientation, random orientation in the horizontal plane, and a two-dimensional axisymmetric Gaussian distribution of orientations. In this paper, for all simulations, we used aspect ratios proposed by Brandes et al. (2002) for rain drops, 0.2 for cloud ice, 0.6 for snow, those by Ryzhkov et al. (2011) for graupel and hail, and the two-dimensional axisymmetric Gaussian distribution for all hydrometeor species.

A hydrometer in CR-SIM is either pure liquid or a mixture of ice and air. The composition of particles within a volume is represented in the scattering computations by an appropriate selection of the dielectric constant for different hydrometeor types. The dielectric constant of liquid particles is frequency dependent and temperature dependent (Ray, 1972). Ice phase hydrometeors are assumed to be composed of ice and air in an ice matrix, and their effective dielectric constant is computed using the Maxwell Garnet mixing formula (Maxwell Garnet, 1904). The output radar reflectivity (Z_hh) for all hydrometeor species is the equivalent radar reflectivity, for which the computations adopt a dielectric factor of 0.92 for all hydrometeor species. This choice of the dielectric factor ensures a convention that the definition of radar reflectivity reduces to the following form: $Z=\int N\left(D\right)D^{\mathrm{6}}\mathrm{d}D$ for (spherical) liquid particles, where D is the droplet diameter and N(D) is the droplet size distribution function.

The LUTs of scattering properties currently incorporated into CR-SIM were created using the T-matrix method where solid phase hydrometeors are represented as dielectrically dry oblate spheroids. Although these assumptions are rather simple compared to some other radar simulators that take into account complex electromagnetic scattering by mixed-phase hydrometeors or ice hydrometeors with possibly irregular shapes (e.g., Snyder et al., 2017a, b; Jiang et al., 2019), such complex electromagnetic scatters can be easily incorporated by adding LUTs of their scattering properties from different scattering calculation methods (e.g., Kneifel et al., 2017; Leinonen and Moisseev, 2015; Leinonen and Szyrmer, 2015; Lu et al., 2016) without any structural change to the code.

2.2 Calculations of radar and lidar observables

The PSD for each hydrometeor species is produced based on the model microphysics scheme. The incorporated microphysics schemes and corresponding CRM currently available in CR-SIM are listed in Table 3. CRMs coupled with bulk moment microphysics (i.e., single and multi-moments) prognose mixing ratio and/or the number concentration of each hydrometeor species. Using these parameters, the PSD is determined in combination with size distribution assumptions used in the microphysics scheme. Bin microphysics schemes explicitly calculate the evolution of the PSDs. Radar moment observables are computed by integrating scattering properties (see Appendix A) from the LUTs over the discrete PSD for each hydrometeor type following Ryzhkov et al. (2011) while accounting for an orientation distribution as described in Sect. 2.1. They are then integrated over all simulated hydrometeor species to produce a unique value for each observable at each grid box (see Tatarevic et al., 2019, for detailed descriptions). Computed radar variables are listed in Table 4.

Table 3 Incorporated microphysics schemes and corresponding CRMs.

Table 4 Computed radar variables. Units of the variables stored in the output files are provided in parenthesis.

Particle fall velocity, which is used for Doppler velocity and spectrum width computations, is parameterized as a function of particle diameter in the same manner as in the selected microphysics scheme. The fall velocity size relationship (V_f(D)) for each hydrometeor species is specified in the following form:

$$\begin{array}{}\text{(1)}& V_{\mathrm{f}}\left(D\right)=f_{\mathrm{c}}{a}_{\mathrm{v}}{D}^{b_{\text{v}}},\end{array}$$

where $f_{\mathrm{c}}={\left({\mathit{\rho }}_{\text{surf}}/\mathit{\rho }\right)}^k$ is the correction factor for air density with exponent k. The values for the prefactor a_v, the exponent b_v, and the exponent k vary according to the microphysics scheme and do not depend on particle orientation. The air density at sea level, ρ_surf, is computed for the first model level. The values of coefficients and detailed descriptions concerning the microphysics schemes are found in Tatarevic et al. (2019).

The versatility of the CR-SIM computational capabilities is briefly demonstrated by the following examples.

• Figure 1 shows an example of S-band (3 GHz) radar observables from CR-SIM for a mesoscale convective system (MCS) observed on 20 May 2011, during the Midlatitude Continental Convective Clouds Experiment (MC3E; Jensen et al., 2016). The convective system was simulated using the Weather Research Forecasting (WRF) model (Skamarock et al., 2008) with the Morrison 2-moment microphysics scheme, a horizontal resolution of 0.5 km, and the vertical resolution of approximately 0.25 km.

• CR-SIM includes a computation of the Doppler power spectra by introducing the method used in Kollias et al. (2014). Figure 2 shows examples of the Doppler spectra and its moments at S band for a WRF-simulated cumulus field. In the figure, a pulse repetition frequency (PRF) of 600 Hz is used, the noise power at 1 km is −40 dB, and the number of Doppler velocity bins is 256.

• CR-SIM also includes forward simulators for the ceilometer (wavelength of 905 nm) and ground-based micro pulse lidar (wavelengths of 532 and 353 nm). The lidar observables are computed for cloud ice and cloud droplet species (see Table 5 and Appendix A). Figure 3 shows an example of profiles calculated for lidar observables for a cumulus case from the LES ARM Symbiotic Simulation and Observation project (LASSO; Gustafson et al., 2020) using WRF coupled with the Morrison 2-moment microphysics scheme. In this simulation, typical profiles are presented for aerosol backscatter (β_aero), extinction coefficient (α_{ext_aero}), and molecular backscatter (β_mol) based on Spinhirne (1993). As expected, the lidar backscatter is significantly attenuated by cloud droplets, but the very high lidar backscatter at the interface between air and cloud can be used to detect cloud base height.

Figure 1 Radar observables produced by CR-SIM for a mesoscale convective system. The system was simulated using WRF with the Morrison 2-moment microphysics scheme at 1.8 km altitude. Shown are horizontal cross sections of (a) total hydrometeor content and (b) vertical air velocity from the WRF simulation. CR-SIM produces the following parameters for a scanning S-band (3 GHz) radar located at the center of the domain: (c) Z_hh, (d) Z_DR, (e) K_DP, (f) radar antenna elevation angle, (g) Doppler velocity, and (h) spectrum width at 12:18:00 UTC.

Figure 2 CR-SIM examples of radar observables for a shallow convection LASSO case from a WRF simulation coupled with the Thompson microphysics scheme. Shown are (a) simulated radar Doppler spectra, (b) model vertical velocity (w, solid line) and simulated mean Doppler velocity (V_dop, dashed line), (c) simulated spectrum width (SW, solid line) and simulated reflectivity-weighted velocity (V_fall, dashed line), and (d) simulated total reflectivity (Z_hh) at S band (3 GHz). In (a) and (b), a positive sign indicates upward motion, and in (c) a positive sign indicates downward motion (fall speed).

Table 5 Computed lidar variables.

Figure 3 Lidar observables from CR-SIM for a cumulus case from LASSO using WRF with the Morrison 2-moment microphysics scheme. Example of simulated vertical profiles are shown for β_total, β_{total_atten}, β_mol, and β_aero at a wavelength of 532 nm.

2.3 Instrument model

An instrument model is used to postprocess the CR-SIM results to account for the total (two-way) attenuation in a pathway and the effects of technological specifications on the observables, such as sampling volume and detector sensitivity. The standard output of CR-SIM consists of synthetic profiling radar and lidar observations at each grid box of the model domain and synthetic scanning radar observations for a radar positioned at a user-specified location inside the model domain. The output synthetic fields are artifact free, with no propagation (e.g., total attenuation in a pathway) or instrument sampling effects (e.g., antenna beamwidth, range-gate spacing). This approach is based on the notion that the real observations used for comparison against the synthetic simulated observables will have undergone rigorous postprocessing that mitigate attenuation effects to the extent possible, e.g., velocity folding. However, the user can emulate the true observations of a scanning radar while selecting the placement of a radar or a network of radars within the model domain; thus, imposing a volume coverage pattern (VCP) scan strategy. In this case, the idealized, standard CR-SIM output at the model grid resolution can be further used as input into a radar instrument model that is written specifically for the postprocessing of the CR-SIM radar simulations. The radar instrument model accounts for radar distance to the target, elevation as provided by the VCP, pulse length, range resolution, antenna beamwidth, and receiver noise, and it outputs the radar observables in radar polar coordinates. For the calculation of elevation angles, the Earth's surface is assumed to be flat, which is an acceptable assumption for general radar observation ranges ($<\sim \mathrm{90}\mathrm{km}$ for the vertical model grid spacing of >0.5 km). Gaussian functions are used as the antenna-weighting the range-weighting functions to estimate the contribution of the model grid observables to the radar polar coordinate system observables. Depending on the azimuthal resolution and the antenna beamwidth, this instrument model also accounts for the radar sampling resolution.

Radar reflectivity attenuated for hydrometeors can be computed using the integrated specific attenuation for hydrometeors along a radar beam path. The total hydrometeor attenuation (A_tot) at each grid box is then equal to twice the sum of the specific attenuation for all simulated hydrometeors (A_h) along a radar beam path from the location of the radar to a distance of the target at r in kilometers:

$$\begin{array}{}\text{(2)}& A_{\text{tot}}\left(r\right)=\mathrm{2}\underset{\mathrm{0}}{\overset{r}{\int }}{A}_{\mathrm{h}}\left(r\right)\mathrm{d}r,\end{array}$$

where A_tot is in decibel (dB) and A_h in decibel per kilometer (dB km⁻¹). Here, gaseous attenuation was not included. The observed reflectivity $Z_{\mathrm{hh}}^{\text{obs}}$ (logarithmic scale) is computed by subtracting A_tot from Z_hh on a logarithmic scale:

$$\begin{array}{}\text{(3)}& Z_{\mathrm{hh}}^{\text{obs}}\left(r\right)=Z_{\mathrm{hh}}\left(r\right)-A_{\text{tot}}\left(r\right).\end{array}$$

Like Z_hh, the attenuated differential reflectivity $Z_{\text{DR}}^{\text{obs}}$ on a logarithmic scale is calculated as

$$\begin{array}{}\text{(4)}& Z_{\text{DR}}^{\text{obs}}\left(r\right)=Z_{\text{DR}}\left(r\right)-\mathrm{2}\underset{\mathrm{0}}{\overset{r}{\int }}{A}_{\text{DP}}\left(r\right)\mathrm{d}r,\end{array}$$

where A_DP represents specific differential attenuation in decibel per kilometer (dB km⁻¹). The minimum detectable reflectivity Z_MIN (logarithmic) is applied with a constant C:

$$\begin{array}{}\text{(5)}& Z_{\text{MIN}}=C+\mathrm{20}{\log }_{\mathrm{10}}\left(r\right),\end{array}$$

where r is the radial distance in kilometers, and the constant C represents the minimum detectable signal at r=1 km for the pulse length selected by the user.

Figure 4 shows simulated range–height indicator (RHI) measurements at C and X bands (5.5 and 9 GHz, respectively) accounting for Z_MIN, total hydrometeor attenuation, and the radar range-gate sampling volume for convective cells associated with an MCS observed on 20 May 2011 during MC3E. The input convective system simulation data are the same as Fig. 1. The instrument specifications used for the RHI simulations are for the X-band radar, a beamwidth of 1^∘, range-gate spacing of 50 m, and a constant C of −50 dBZ for the Z_MIN calculation. The C-band radar specifications are a beamwidth of 1^∘, range-gate spacing of 120 m, and a constant C of −35 dBZ for the Z_MIN calculation. These specifications follow the X-band scanning ARM precipitation radar (X-SAPR) and C-band scanning ARM precipitation radar (C-SAPR) configurations at the ARM Southern Great Plains (SGP) site during MC3E. The results are reasonable, showing strong attenuation in Z_hh and Z_DR by rain at X band and relatively less at C band. The simulated K_DP at X band is approximately 1.6 times larger than that at C band because of the wavelength dependency.

Figure 4 Examples of C- and X-band (5.5 and 9 GHz, respectively) RHI scans with a beamwidth of 1^∘ produced using CR-SIM for a convective cell in a mesoscale convective system (MCS). The simulation uses WRF with the Morrison 2-moment microphysics scheme for an MCS on 20 May 2011 during MC3E. Shown are variables at X- and C-band frequencies 15–35 km from the radar as a function of height at 12:18:00 UTC: (top raw) Z_hh, (middle raw) Z_DR, and (bottom raw) K_DP. The figure shows (a) C-band variables without attenuation, (b) X-band variables with attenuation from a radar 24 km from the convective cell, (c) C-band variables with attenuation from a radar 24 km away, and (d) C-band variables from a radar located 59 km away with attenuation.

2.4 Code features

CR-SIM is written in Fortran 95 standard including all GNU extensions and parallelized with OpenMP. The input to CR-SIM is the output from a CRM or LES in NetCDF format. The output of CR-SIM is in NetCDF format and includes simulated observables for each hydrometeor species, and one total for all the hydrometeors. These features allow users to understand the contributions of each hydrometeor species to the radar observables for a sophisticated evaluation of microphysics schemes. The code includes various microphysics schemes as shown in Table 3. The code structure supports different CRMs or LESs; flexible microphysics package extensions; and diverse assumptions such as particle shape, density, and PSD for different hydrometeor categories in the models as well as inclusion of scattering properties computed by different methods.

The CR-SIM runtime depends on the computing power, number of threads used, simulation domain size, and the number of cloudy grid boxes. The simulations presented in this article were run on a computer having 500 GB memory and 24 processors each with 12 cores. For the MCS case in Fig. 1, having a simulation domain size of $\mathrm{667}\times \mathrm{667}\times \mathrm{12}$ (5.3 M grid points), the runtime is approximately 270 s using 16 threads.

The code has been released under GNU General Public License, and both the software and a detailed user guide are publicly available online (Tatarevic et al., 2019).

3 Sample applications of CR-SIM

In this section, several applications of CR-SIM are presented that highlight its capabilities. These applications are (i) a comparison of observed and modeled cloud fraction profiles (CFPs), (ii) a quantification of uncertainty in the estimate of domain-averaged CFP, (iii) an evaluation of a novel retrieval technique for the estimation of cloud fraction, (iv) an investigation of the impacts of limitations imposed by the nature of observations themselves on multi-Doppler wind retrieval techniques, and (v) an optimization of a new radar observation strategy for multi-Doppler wind retrievals.

Figure 5 shows a flow diagram of our application processes. First, the forward simulator produces idealized observables at each model grid box (the “Output 1” box in Fig. 5). In the second step, an instrument model is used to account for the instrument features (as described in Sect. 2.3). Third, the output from the instrument model (“Output 2”) is then used to retrieve the CFP (the retrieval model and “Output 3”) for a direct comparison and, most importantly, for a quantification of the uncertainties in the cloud fraction estimates and evaluation of the new retrieval technique (applications i–iii). On the other hand, the output from the instrument model is also used as an input for multi-Doppler wind retrieval model to investigate the uncertainty of the retrieval method and to optimize the new radar observation strategy (applications iv and v, with “Output 3”). The final step consists of a comparison of the retrieved quantities using a multi-Doppler wind retrieval against the input CRM or LES data and a quantitative estimation of uncertainties attributed to the observation limitations and the retrieval algorithms. In the following subsections, we briefly describe and summarize the findings of the studies using CR-SIM.

Figure 5 Diagram for CR-SIM and its applications. The diagram indicates the CR-SIM input and the different levels of output for the forward model, instrument model, and retrieval model.

3.1 Comparison of observed and modeled cloud fraction profiles

Measurements of the CFP are important to quantify the impact of shallow cumulus clouds on the grid-scale meteorological state, because the fractional cloudiness of a grid box has an impact on the radiative transfer (e.g., Albrecht, 1981; Larson et al., 2001) and the vertical cumulus mass flux (e.g., de Roode and Bretheton, 2003; van Stratum et al., 2014). Zenith-profiling cloud radar and lidar measurements traditionally have been used to provide CFP estimates (e.g., Hogan et al., 2001; Kollias et al., 2009; Remillard et al., 2013; Angevine et al., 2018). Typically, the profiling radar and lidar observations are combined synergistically to provide a hydrometeor mask such as those described in ARSCL (Clothiaux et al., 2000) and the CloudNet target classification (Illingworth et al., 2007). This approach takes advantage of the radar and lidar capabilities and maximizes our ability to detect thin cloud layers. However, the performance of the combined radar and lidar algorithm degrades at heights where the lidar observations are unavailable due to complete signal attenuation. These attenuation effects are naturally not represented in model output and thus may lead to large disagreements between observations and models.

Our focus is to use CR-SIM to generate a synthetic ARSCL product that is directly comparable to the ARSCL generated using measurements from the Ka-band ARM Zenith-pointing Radar (KAZR), ceilometer, and MPL. This analysis uses a shallow cumulus cloud field over ARM SGP site simulated by the LASSO project. The simulation is for the 27 June 2015 case and uses WRF-LES coupled with the Morrison 2-moment microphysics scheme (Morrison et al., 2005). The horizontal and vertical resolutions are 100 and 20 m, respectively, and the horizontal dimension of the simulation domain is 14.4 km.

First, the KAZR, ceilometer, and MPL measurements from the ARM SGP Central Facility are simulated using the CR-SIM forward simulator. The simulation output corresponds to the box “Output 2” in Fig. 5. Simulated KAZR reflectivity accounts for hydrometeor attenuation ($Z_{\mathrm{hh}}^{\text{obs}}$) and radar sensitivity (Z_MIN) as described by Eqs. (2, 3, and 5). The attenuated MPL hydrometeor backscatter (β_{hydro_atten}) is obtained by subtracting β_{aero_atten} and β_{mol_atten} from β_{total_atten}, since the MPL total backscatter includes aerosol backscatter and molecular backscatter (see Table 5). The value obtained for β_{hydro_atten} is below noise level if less than the unattenuated background scatter (β_aero+β_mol), which is used in this simulation as the minimum detectable MPL backscatter value. The ceilometer-detected first cloud base is estimated at each grid column following O'Connor et al. (2004). Using the simulated observables, we estimate cloud locations as provided by ARSCL (“Output 3” in Fig. 5). A grid box where either KAZR $Z_{\mathrm{hh}}^{\text{obs}}$ or MPL β_{hydro_atten} has a detectable value is indexed as a “cloudy” grid box, and grid boxes below the simulated ceilometer first-cloud base are indexed as “clear”.

An example of an ARSCL simulation is shown in Fig. 6 that uses the LASSO LES data as an input. The WRF simulation shows cumulus clouds below 5 km and cirrus clouds covering the entire domain at 12–14 km. In Fig. 6b–d, the limitation of each instrument is represented in the forward simulations. The simulated KAZR measurements can detect cumulus cloud layers but cannot detect cirrus clouds, due to their low reflectivity value (lower than Z_MIN). Instead, the cirrus clouds can be detected by the simulated MPL measurements. However, the cirrus clouds may be missed by both radar and lidar measurements when cumulus clouds are present, because the MPL signal becomes fully attenuated by the low-level clouds. Figure 6f and g show the domain-averaged CFPs from the LES hydrometeor mixing ratio and from the simulated ARSCL, which assumes the ARM instruments are located at every grid column (as shown in Fig. 6e). Comparison between the two CFP plots suggests that the ARSCL for this LASSO case underestimates cirrus CFPs by 20 %, likely due to lidar beam attenuation by lower-level cumulus clouds that have a horizontal fraction of 20 %.

Figure 6 Simulated vertically pointing radar and lidar measurements and the ARSCL product for a shallow convection case on 27 June 2015. (a–e) Vertical cross sections of (a) water content from the WRF model, (b) Ka-band (35 GHz) radar reflectivity accounting for radar sensitivity and attenuation, (c) MPL attenuated backscatter, (d) ceilometer backscatter (color bar) and first cloud base (gray dots), and (e) the ARSCL cloud mask. (f, g) Height versus time cross sections of domain-averaged cloud fraction from (f) WRF water content $>\mathrm{0.001}\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$ and (g) the simulated ARSCL product.

3.2 Uncertainty quantification of domain-averaged cloud fraction profile estimates

The ARSCL product is usually integrated for 1–3 h to provide an average CFP estimate for that time period. These CFP estimates are often compared with the model domain-averaged CFPs. However, the spatially heterogeneous distribution of the shallow cumulus clouds (Wood and Field, 2011) raises questions regarding the ability of short-term (1–3 h) zenith-profiling observations to provide adequate sampling of the cloud field. Uncertainties in the profiling measurement of cloud fractions are introduced by the limited sampling of a highly heterogeneous cloud field. We investigate these uncertainties as a function of the number of profiling sites and integration time using the CR-SIM virtual observations and the WRF LES simulation presented in the previous application. The WRF LES output is saved every 10 min, and CR-SIM is run for each output file. In this analysis, we assume that no cloud evolution occurs within a 10 min interval.

Figure 7a and b show the domain-averaged CFP from the simulated ARSCL and directly from the WRF LES using a cloud water content threshold of 0.01 g kg⁻¹. The colors indicate different integration periods. Note that the WRF dataset in this analysis is for a shallow convective case at SGP on 11 June 2016, different from the one used in Fig. 6, which has higher cumulus cloud tops. The simulated ARSCL CFP is in good agreement with the WRF CFP for each integrated period (compare Fig. 7a and b), indicating that uncertainties attributed to observation limitations (e.g., sensitivity and attenuation) are small. Thus, the limited spatial and/or temporal sampling is the major error source to consider when comparing the profiling measurement derived CFP with the domain-averaged WRF CFP.

Figure 7 Investigation of errors of cloud fraction profiles (CFPs) from profiling measurements. CFPs from single sites are estimated by integrating over time and then averaged. Shown are domain-averaged cloud fraction profiles (CFPs) from (a) WRF-simulated cloud water mixing ratio and (b) the simulated ARSCL product for a shallow convection case on 11 June 2016. Colors in (a) and (b) represent different integration time periods centered at 21:00:00 UTC. The minimum threshold for the WRF cloud water mixing ratio is 0.01 g kg⁻¹. (c–h) CFPs from the simulated ARSCL with different number of observation sites N and different integration periods T. The black line in (c–h) represents the domain-averaged CFP from the WRF-simulated cloud water mixing ratio, blue lines represent CFPs from individual observation sites, and the red line represents the mean CFP from averaging over the individual sites. (i, j) Root mean square error (i) and mean absolute percentage error (j) of the simulated ARSCL CFPs as a function of the number of observation sites and integration period.

To emulate vertical profiling measurements, we sampled data as follows. First, observation sites are randomly selected within the horizontal domain. Second, for each snapshot of the simulation, clouds over the observation sites are sampled as if the clouds are frozen in time and advected by the mean environmental wind. Thus, the columns are sampled along the direction of the horizontal wind over the advected distance (i.e., horizontal wind speed ×10 min), where the environmental horizontal wind at each snapshot is the mean horizontal wind across the simulation domain within the cumulus cloud layer (i.e., the layer between the mean cumulus cloud base and the maximum cumulus cloud top). Last, the CFP is estimated by varying both the number of observation sites and the integration period.

Figure 7c–h show the comparisons of the WRF CFPs and the simulated ARSCL CFPs for different number of observation sites (top row, c–e) and integration periods (bottom row, f–h). The center of the integration period is 21:00 UTC. Blue lines in each panel represent the simulated ARSCL CFPs integrated over time from each selected observation site for the period indicated, the red line represents the mean ARSCL CFP averaged over the sites, and the black line represents the domain-averaged WRF CFP integrated over the indicated period. Each panel shows that the CFPs at a single site (blue lines) have large uncertainties even though they are integrated over long periods, ranging from 30 to 180 min. Those uncertainties are reduced when averaging the CFP profiles across the different sites; consequently, the mean CFP (red line) becomes closer to the domain-averaged WRF CFP (black lines). However, it also becomes evident that a small number of observation sites (Fig. 7c) may not be adequate to estimate the true CFP.

Figure 7i and j show the root mean square error (RMSE) and mean absolute percentage error (MAPE) of the simulated ARSCL CFPs as a function of the number of observation sites and the integration time. Both plots show that the uncertainty can be reduced by increasing the number of observation sites and the length of the integration period. The RMSE dramatically decreases to 0.005–0.01 (30 %–50 % in MAPE) when we use four observation sites and 120 min integration. The rate of improvement of CFP, by further increasing the number of sites and integration period, is smaller; the error values slowly decrease until the RMSE and MAPE plateau at 0.002 % and 15 %, respectively. However, establishing more than 10 observational sites in such a small domain is probably impractical. At the SGP site, five Doppler lidar profiling measurements have already been operating over a 90 km×90 km domain. These measurements can be effectively used to estimate cloud fraction without much uncertainty when clouds are homogeneously distributed over the domain.

3.3 Evaluation of a new CFP estimation technique using scanning cloud radar

Forward-radar simulators can be used to evaluate retrieval techniques. We introduce an application for estimating CFP using scanning cloud radar (SCR) measurements presented in Oue et al. (2016). As analyzed in the previous section, profiling radar measurements may produce large uncertainties in CFP estimates. On the other hand, SCRs conduct observations over a much larger domain than zenith-profiling cloud radars such as KAZR (e.g., Lamer et al., 2013; Ewald et al., 2015). Although the SCRs are widely and routinely used to observe 3D cloud fields, the application of SCRs to study shallow cumuli is not straightforward. One of the most significant limitations of the SCR observations is related to the radar sensitivity. Since shallow cumuli over land typically have low reflectivities, the strong reduction in SCR sensitivity with range creates the illusion of the atmosphere being cloudier closer to the radar location (e.g., Lamer and Kollias, 2015). This limitation can introduce uncertainties in the cloud fraction estimates. Oue et al. (2016) addressed uncertainties of radar-estimated CFPs due to the nature of the profiling and scanning radar techniques using CR-SIM-generated observations.

Figure 8 Horizontal cross sections of (a) water content simulated by WRF and (b) Ka-band (35 GHz) Z_hh simulated at 2.4 km above ground level for a LASSO case. In (b), it is assumed that the radar is located at x=0 km and the RHI is scanned along the east–west axis, and the radar sensitivity Z_MIN with $C_{\mathrm{0}}=-\mathrm{50}\mathrm{dBZ}$ was applied. (c) Cloud fraction profiles corresponding to the 10 % cumulative distribution function (CDF) isoline with changing integration time of CWRHI (hence, number of scans). (d) The root-mean-square error (RMSE) from the LES domain-averaged CFP for CDF isolines of 5 % (thin solid line), 10 % (thick solid line), 15 % (dashed line), and 20 % (dashed-dotted line) as a function of integration time. The black dashed line in (c) represents the LES domain-averaged CFP for hydrometeor mixing ratio $\ge \mathrm{0.01}\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}$. (c) and (d) are adapted from Oue et al. (2016).

Figure 8a shows horizontal cross sections of WRF-simulated water content for a shallow convection case (9 June 2015; Oue et al., 2016) from LASSO. Figure 8b shows the CR-SIM simulation of the Ka band (35 GHz) Z_hh accounting for the minimum detectable reflectivity, Z_MIN, of the cross-wind RHI (CWRHI; Kollias et al., 2014) scans from Eq. (5) using $C=-\mathrm{50}\mathrm{dBZ}$. In the CR-SIM analysis, the radar was located along the vertical line in Fig. 8b, and CWRHI scans were performed along the east–west direction, while the clouds were assumed to move along the north–south direction. These figures suggest that the Z_hh from the CWRHI scans cannot capture the clouds with lower water contents that are located far from the radar. This can affect cloud fraction estimates. Because the “true” cloud fraction is estimated from the original model cloud field and thus is known, the CR-SIM runs in different configurations can be used to establish the best method to estimate the cloud fraction while accounting for limitations inherent to the nature of radar measurements. Oue et al. (2016) use the cumulative distribution function (CDF) of the observed Z_hh to define the size of the horizontal domain at each height needed to obtain the best estimate of the domain-averaged CFP. The horizontal domain size as a function of height corresponds to a distance from the radar where Z_MIN is equal to a CDF value of 10 %. Figure 8c shows CFPs using a CDF of 10 % when changing the integration time of the CWRHI, and Fig. 8d shows the RMSE of the estimated CFPs as a function of integration time, adapted from Oue et al. (2016). The figure suggests that the 35 min or longer of CWRHI measurements provide the realistic domain-averaged CFP.

3.4 Investigation of impacts of observation limitations on multi-Doppler radar wind retrievals

Estimation of vertical air motion is essential to understand the dynamics and microphysics of deep convective clouds (e.g., Jorgensen and LeMone, 1989; Wang et al., 2019), evaluate CRM and LES results (e.g., Varble et al., 2014; Fan et al., 2017), and improve convective parameterization in GCMs (e.g., Donner et al., 2001). Multi-Doppler radar techniques have been applied to understand the dynamics and microphysics of the deep convective clouds in different climate regimes (e.g., Friedrich and Hagen, 2004; Collis et al., 2013; Oue et al., 2014). However, the multi-Doppler radar retrievals are not straightforward with potential uncertainties from multiple aspects (e.g., Clark et al., 1980; Bousquet et al., 2008; Potvin et al., 2012). CR-SIM can be used to investigate the impacts of different error sources on the retrieved wind fields.

Oue et al. (2019a) investigated the impacts of the radar VCP for the plan position indicator (PPI) and the observation period on uncertainties in multi-Doppler radar wind retrievals using CR-SIM. They also investigated how the uncertainties attributed to the VCP period can be reduced using the advection-correction technique proposed by Shapiro et al. (2010). The advection correction scheme allows for trajectories of multiple individual clouds, performs smooth grid-box-by-grid-box corrections of cloud locations, and takes into account changes in cloud shape with time by using PPI scans at two times. We summarize their findings, particularly regarding the impacts of radar VCP and period on multi-Doppler radar retrievals.

Figure 9 A diagram of an Observing System Simulation Experiment study to investigate the impacts of radar volume coverage pattern (VCP) on a multi-Doppler radar wind retrieval.

Figure 9 shows a diagram of the analysis process. The input model data are a WRF simulation using the Morrison 2-moment microphysics scheme for an MCS observed on 20 May 2011, during the MC3E field campaign at the ARM SGP site. The horizontal resolution is 500 m; the vertical resolution varies from approximately 30 m near the surface to 260 m at 2 km – above which the resolution remains approximately constant – and the simulation output is saved every 20 s. Measurements from the three X-band scanning ARM precipitation radars (X-SAPR) at the SGP site are simulated using CR-SIM. The CR-SIM-simulated radar reflectivity and Doppler velocity at the model grid are converted into the radar polar coordinates with two different VCPs for each radar: (1) 21 elevation angles ranging from 0.5 to 45^∘ (VCP1, same as the X-SAPR scan strategy during MC3E) and (2) 60 elevation angles ranging from 0.5 to 59.5^∘ with a 1^∘ increment (VCP2). For the both VCPs, the beamwidth is 1^∘, the range-gate spacing is 50 m, and the maximum range is 40 km. The simulated radar reflectivity and Doppler velocity in polar coordinates were used as an input to the 3D-Var multi-Doppler radar wind retrieval algorithm developed by North et al. (2017) to estimate the 3D wind field for a domain of $\mathrm{50}\mathrm{km}\times \mathrm{50}\mathrm{km}\times \mathrm{10}\mathrm{km}$ with horizontal and vertical grid spacings of 0.25 km.

The convective mass flux (MF) is estimated at each height as

$$\begin{array}{}\text{(6)}& \text{MF}=\text{UF}\stackrel{\mathrm{‾}}{w}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{d}}\left[\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}\right],\end{array}$$

where UF is updraft fraction over the horizontal slice of the domain, $\stackrel{\mathrm{‾}}{w}$ is mean vertical velocity over the updraft area, and ${\stackrel{\mathrm{‾}}{\mathit{\rho }}}_{\mathrm{d}}$ is dry air density averaged over the domain. Figure 10 shows comparisons of convective mass flux profiles between simulated multi-Doppler radar retrievals and WRF output for two minimum updraft thresholds of 5 m s⁻¹ (MF₅) and 10 m s⁻¹ (MF₁₀). First, we applied the wind retrieval technique to a snapshot of the forward-model output to bypass the instrument model and examine the uncertainty in the retrieval model (3FullGrid). Figure 10a shows MF profiles from the 3FullGrid simulation (red line) and from the WRF snapshot (black line), with 2 min average (dark gray line) and 5 min average (light gray line). The 3FullGrid MF profile is in good agreement with the WRF output, indicating that the uncertainty in the retrieval model is small; although, it does underestimate the maximum MF for the updraft threshold of 5 m s⁻¹ by 0.05 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (10 % of the true MF) at 5.3 km.

Figure 10 Vertical profiles of convective mass flux with updraft thresholds of 5 m s⁻¹ (solid lines) and 10 m s⁻¹ (dashed lines). Displayed in each panel are different retrieval simulations represented by the colors. The dark gray line in (a) represents the time average of the WRF output over 2 min, and the light gray line in (a) represents the time average of the WRF output over 5 min. The profile from the WRF snapshot is displayed in each panel by a black solid line. Adapted from Oue et al. (2019a).

Figure 10b and c show MF profiles (MF₅ and MF₁₀) obtained from simulated retrievals while considering the effects of VCP (VCP1 and VCP2) and averaging period (“snap” (instantaneous), 2 min, and 5 min averages). For both VCP1 (Fig. 10b) and VCP2 (Fig. 10c), the snapshot and 2 min VCP simulations have similar MF estimates for both sets of MF₅ and MF₁₀ curves, indicating that a 2 min average is sufficient to capture features available from an instantaneous scan. However, the accuracy of these estimates varies with MF profile and VCP. The MF₁₀ estimates for both VCPs systematically underestimate the maximum values occurring between 4.5 and 6.5 km by about 0.5 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (20 %). The performance of the MF₅ estimates for VCP snap and 2 min have strong variations with height. For VCP1 (the less dense scan pattern), MF₅ follows the WRF snapshot below 4.5 km with close agreement between 3 and 4.5 km; however, MF is underestimated around its maximum MF by about 0.075 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ (15 %) and is overestimated below 3 km and above 7 km. The denser scan pattern for VCP2 provides a dramatic improvement around the maximum and above 6 km, but it still shows overestimations below 3 km and above 7 km. Uncertainties are often increased for the VCP simulations when the averaging period is extended to 5 min. For the 5 min VCPs, MF₁₀ estimates for both VCP1 and VCP2 around the maximum are further underestimated, while the MF₅ estimate for VCP2 is further overestimated above 6 km. Other estimates below this height for VCP2 and for all heights for VCP1 are mostly unchanged. These results suggest that the VCP elevation strategy and sampling time extended to 5 min have a significant impact on the updraft properties retrieved at higher altitudes. This is due to the density of data sampled by the VCPs, where greater density particularly improves MF₅ around its maximum, and the deformation of cloud structures within longer sampling periods (exceeding 2 min) that causes uncertainties in the mass continuity assumption.

The rapid volume scan of less than 5 min required in the retrieval of the high-quality vertical velocities is challenging for conventional scanning radars. Most of the improvements required in the sampling strategy of the observing system (higher maximum elevation angle, higher density elevation angles, and rapid VCP time period) can be accomplished using rapid scan radar systems such as the Doppler on Wheels mobile radars (DOW; Wurman, 2001) or phased array radar systems (e.g., Kollias et al., 2018).

3.5 Evaluation of new radar observation strategies

CR-SIM can also be used to examine the performances of new remote-sensing systems and help to choose the most appropriate observation strategy for a new field campaign. Figure 4c and d show the performance of C-band (5.5 GHz) RHI measurements when the radar is located at 24 and 59 km away from the target convective clouds. As expected, the RHI from the greater distance provides the radar observables at lower resolution and includes more attenuation when precipitation clouds are located between the target and the radar. Oue et al. (2019a) investigate the impact of radar data sampling on the multi-Doppler radar wind retrievals for the MCS by an OSSE using CR-SIM. The addition of data from a Doppler radar to form a triple-Doppler radar retrieval, shown in Sect. 3.4, cannot significantly improve the updraft retrievals if the added radar VCP has inferior spatial resolution. Oue et al. (2019a) also show that the updraft retrievals in a limited area around the center of the domain, where data density from the three radars is higher than other areas, produced better results than those in the entire domain. The insights obtained from these OSSEs are beneficial for decision-making regarding radar observation strategies for a field campaign, such as the number of radars required and their locations. For example, Kollias et al. (2018) used CR-SIM to examine how phased array radars improve multi-Doppler radar wind retrievals compared to scanning radars for MCSs.

4 Summary

We present a recently developed comprehensive forward simulator for radar and lidar, CR-SIM, which is suitable for simulating complex, ground-based observational configurations and their synthetic products. CR-SIM can simulate multiwavelength, zenith-pointing and scanning radar observables (radar reflectivity, Doppler velocity, polarimetric fields, radar Doppler spectrum), lidar observables, and multisensor integrated products. The primary idea behind the simulator is to directly compare remote-sensing observations with simulated measurements based on the CRM or LES output, maintaining consistency with the microphysics scheme used in the model. CR-SIM incorporates microphysical and scattering properties independently, so that uncertainties related to microphysical assumptions are separated from uncertainties related to scattering model. This configuration allows CR-SIM to be easily expanded, either by adding microphysical schemes or new scattering classes.

One feature of CR-SIM is that it produces both radar and lidar observables for all the CRM grid boxes while accounting for elevation angles relative to a radar location, similar to Snyder et al. (2017a, b). Other radar simulators also account for radar geometry characteristics such as beamwidth and radar range resolution to simulate scatterers within the radar resolution volume (e.g., Capsoni et al., 2001; Caumont et al., 2006; Cheong et al., 2008). Instead, here, the radar sampling characteristics (e.g., antenna beamwidth, range-bin spacing, total attenuation, sensitivity) are accounted for in our postprocessing instrument model. This feature facilitates configuring any desirable observational setup with a varying number of profiling or scanning sensors from a single model simulation. The CR-SIM multisensor simulations include multiwavelength radars and lidars that allow for simulation of sophisticated virtual products such as ARSCL and 3D-Var multi-Doppler-based wind retrievals. The CR-SIM applications shown in this paper emphasize the value of applying it to high-resolution model output to emulate the sampling by the ground-based observatories. CR-SIM's coupling of CRM microphysical parameterizations with scattering models facilitates improved evaluation of model performance by enabling robust comparisons between model-simulated clouds and observables from radar and lidar while accounting for instrument characteristics and observation limitations.

CR-SIM is easily expanded to include additional microphysical schemes, new scattering classes, scattering calculations, and other applications to simulate multisensor products (e.g., multi-Doppler wind retrievals, ARSCL). At this stage, all ice hydrometeors (e.g., snow, ice, graupel, hail) are modeled as dielectrically dry spheroids. The LUTs of scattering properties incorporated into the current CR-SIM were created using the T-matrix method based on assumptions regarding ice particle composition and shape. More single-scattering properties from other scattering calculation methods can be incorporated by adding LUTs. Moreover, the gaseous attenuation will be considered in the future, as gaseous attenuation can be significant in the millimeter-wavelength radar measurements, and elevation angles will be corrected for the Earth's curvature. The analyses presented here serve as a reference to the CR-SIM package and illustrate numerous applications related to sampling uncertainty, sampling optimization, retrieval uncertainty, and comparison between models and observations.

Appendix A

The radar observables are computed on the basis of the following equations. In the equations below, M is the number of different hydrometeor species coexisting in the same spatial resolution volume of the model, and the subscript i is the index of species; backward and forward scattering amplitudes are denoted as $f_{a,b}^{\left(\mathit{\pi }\right)}$ and $f_{a,b}^{\left(\mathrm{0}\right)}$, respectively. The superscript * in an expression […]^∗ denotes conjugation; R_e[…] and I_m[…] represent the real and imaginary parts of the complex number, respectively; and $\left|\mathrm{\dots }\right|$ refers to the magnitude of the value between the single bars. λ is the radar wavelength in millimeters, and ${\left|K_{\mathrm{w}}\right|}^{\mathrm{2}}$ is the dielectric factor (the value for water =0.92 is used for all hydrometeor species in CR-SIM). The scattering amplitudes are given in millimeters. N_i(D) defines the particle size distribution in terms of the number of particles per unit volume of air and unit bin size (in ${\mathrm{m}}^{-\mathrm{3}}{\mathrm{mm}}^{-\mathrm{1}}$), with the bin equivolume diameter D (in mm). Both the bin fall velocity, V_Fi, and vertical air velocity, w, are given in meters per second. The elevation and azimuth angles are denoted θ and φ, respectively, and u and v are the two components of horizontal wind. The coefficients A_1i–A_7i are the angular moments for one of the three horizontal orientation expressions taken from Ryzhkov et al. (2011).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{Z}_{\mathrm{hh}}={\displaystyle \frac{\mathrm{4}{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{4}}{\left|K_{\mathrm{w}}\right|}^{\mathrm{2}}}}\times \\ {\displaystyle }& {\displaystyle }\sum _{i=\mathrm{1}}^M[\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left\{{\left|f_{bi}^{\left(\mathit{\pi }\right)}\right|}^{\mathrm{2}}-\mathrm{2}{R}_e\left[f_{bi}^{\left(\mathit{\pi }\right)\ast }\left(f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\left(\mathit{\pi }\right)}\right)\right]A_{\mathrm{2}i}\right.\\ {\displaystyle }& {\displaystyle }+\left.{\left|f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\mathit{\pi }}\right|}^{\mathrm{2}}{A}_{\mathrm{4}i}\right\}{N}_i\left(D\right)\mathrm{d}D]\left[{\displaystyle \frac{{\mathrm{mm}}^{\mathrm{6}}}{{\mathrm{m}}^{\mathrm{3}}}}\right],\\ {\displaystyle }& {\displaystyle }{Z}_{\mathrm{vv}}={\displaystyle \frac{\mathrm{4}{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{4}}{\left|K_{\mathrm{w}}\right|}^{\mathrm{2}}}}\times \\ {\displaystyle }& {\displaystyle }\sum _{i=\mathrm{1}}^M[\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left\{{\left|f_{bi}^{\left(\mathit{\pi }\right)}\right|}^{\mathrm{2}}-\mathrm{2}{R}_e\left[f_{bi}^{\left(\mathit{\pi }\right)\ast }\left(f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\left(\mathit{\pi }\right)}\right)\right]A_{\mathrm{1}i}\right.\\ {\displaystyle }& {\displaystyle }+\left.{\left|f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\mathit{\pi }}\right|}^{\mathrm{2}}{A}_{\mathrm{3}i}\right\}{N}_i\left(D\right)\mathrm{d}D]\left[{\displaystyle \frac{{\mathrm{mm}}^{\mathrm{6}}}{{\mathrm{m}}^{\mathrm{3}}}}\right],\\ {\displaystyle }& {\displaystyle }{Z}_{\mathrm{vh}}={\displaystyle \frac{\mathrm{4}{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{4}}{\left|K_{\mathrm{w}}\right|}^{\mathrm{2}}}}\times \\ {\displaystyle }& {\displaystyle }\sum _{i=\mathrm{1}}^M\left[\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left\{{\left|f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\mathit{\pi }}\right|}^{\mathrm{2}}{A}_{\mathrm{5}i}\right\}{N}_i\left(D\right)\mathrm{d}D\right]\left[{\displaystyle \frac{{\mathrm{mm}}^{\mathrm{6}}}{{\mathrm{m}}^{\mathrm{3}}}}\right],\\ {\displaystyle }& {\displaystyle }{Z}_{\text{DR}}={\displaystyle \frac{Z_{\mathrm{hh}}}{Z_{\mathrm{vv}}}}\left[-\right],\\ {\displaystyle }& {\displaystyle }{\text{LDR}}_h={\displaystyle \frac{Z_{\mathrm{vh}}}{Z_{\mathrm{hh}}}}\left[-\right],\\ {\displaystyle }& {\displaystyle }{K}_{\text{DP}}={\displaystyle \frac{\mathrm{180}\mathit{\lambda }}{\mathit{\pi }}}{\mathrm{10}}^{-\mathrm{3}}\\ {\displaystyle }& {\displaystyle }\times \sum _{i=\mathrm{1}}^M\left[\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left\{R_e\left(f_{bi}^{\left(\mathrm{0}\right)}-f_{ai}^{\mathrm{0}}\right)A_{\mathrm{7}i}\right\}{N}_i\left(D\right)\mathrm{d}D\right]\left[{\displaystyle \frac{{}^{\circ }}{\mathrm{km}}}\right],\\ {\displaystyle }& {\displaystyle }\mathit{\delta }={\displaystyle \frac{\mathrm{180}\mathit{\lambda }}{\mathit{\pi }}}\arg \left[\sum _{i=\mathrm{1}}^M\right[\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathit{\{}{\left|f_{bi}^{\left(\mathit{\pi }\right)}\right|}^{\mathrm{2}}+\left[{\left|f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\left(\mathit{\pi }\right)}\right|}^{\mathrm{2}}\right]A_{\mathrm{3}i}\\ {\displaystyle }& {\displaystyle }-\left[f_{bi}^{\left(\mathit{\pi }\right)\ast }\left(f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\left(\mathit{\pi }\right)}\right)\right]A_{\mathrm{1}i}-\left[f_{bi}^{\left(\mathit{\pi }\right)}\left(f_{bi}^{\left(\mathit{\pi }\right)\ast }-f_{ai}^{\left(\mathit{\pi }\right)\ast }\right)\right]A_{\mathrm{2}i}\mathit{\}}\\ {\displaystyle }& {\displaystyle }\times N_i\left(D\right)\mathrm{d}D\left]\right]{[}^{\circ }],\\ {\displaystyle }& {\displaystyle }{A}_{\mathrm{h}}={\displaystyle \frac{\mathrm{10}}{\log \left(\mathrm{10}\right)}}\mathrm{2}\mathit{\lambda }{\mathrm{10}}^{-\mathrm{3}}\sum _{i=\mathrm{1}}^M\left[\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathit{\right\{}{I}_m\left(f_{bi}^{\left(\mathrm{0}\right)}\right)\\ {\displaystyle }& {\displaystyle }-I_m\left(f_{bi}^{\left(\mathrm{0}\right)}-f_{ai}^{\left(\mathrm{0}\right)}\right)A_{\mathrm{2}i}\mathit{\left\}}{N}_i\left(D\right)\mathrm{d}D\right]\left[{\displaystyle \frac{\mathrm{dB}}{\mathrm{km}}}\right],\\ {\displaystyle }& {\displaystyle }{A}_{\mathrm{v}}={\displaystyle \frac{\mathrm{10}}{\log \left(\mathrm{10}\right)}}\mathrm{2}\mathit{\lambda }{\mathrm{10}}^{-\mathrm{3}}\sum _{i=\mathrm{1}}^M\left[\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathit{\right\{}{I}_m\left(f_{bi}^{\left(\mathrm{0}\right)}\right)\\ {\displaystyle }& {\displaystyle }-I_m\left(f_{bi}^{\left(\mathrm{0}\right)}-f_{ai}^{\left(\mathrm{0}\right)}\right)A_{\mathrm{1}i}\mathit{\left\}}{N}_i\left(D\right)\mathrm{d}D\right]\left[{\displaystyle \frac{\mathrm{dB}}{\mathrm{km}}}\right],\\ {\displaystyle }& {\displaystyle }{A}_{\text{DP}}=A_{\mathrm{h}}-A_{\mathrm{v}}\left[{\displaystyle \frac{\mathrm{dB}}{\mathrm{km}}}\right],\\ {\displaystyle }& {\displaystyle }V{Z}_i=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathit{\{}{\left|f_{bi}^{\left(\mathit{\pi }\right)}\right|}^{\mathrm{2}}-\mathrm{2}{R}_e\left[f_{bi}^{\left(\mathit{\pi }\right)\ast }\left(f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\left(\mathit{\pi }\right)}\right)\right]A_{\mathrm{2}i}\\ {\displaystyle }& {\displaystyle }+{\left|f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\mathit{\pi }}\right|}^{\mathrm{2}}{A}_{\mathrm{4}i}\mathit{\}}{V}_{fi}\left(D\right)N_i\left(D\right)\mathrm{d}D\left[{\displaystyle \frac{{\mathrm{mm}}^{\mathrm{6}}}{{\mathrm{m}}^{\mathrm{3}}}}{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}\right],\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{V}_{\mathrm{RW}}={\displaystyle \frac{\frac{\mathrm{4}{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{4}}{\left|K_{\mathrm{w}}\right|}^{\mathrm{2}}}\sum _{i=\mathrm{1}}^M\left[V{Z}_i\right]}{Z_{\mathrm{hh}}}}\left[{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}\right],\\ {\displaystyle }& {\displaystyle }{V}_{\mathrm{D}\mathrm{\_}\mathrm{90}}={\displaystyle \frac{\frac{\mathrm{4}{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{4}}{\left|K_{\mathrm{w}}\right|}^{\mathrm{2}}}\sum _{i=\mathrm{1}}^M\left[w{Z}_{\mathrm{hh}}-V{Z}_i\right]}{Z_{\mathrm{hh}}}}\left[{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}\right],\\ {\displaystyle }& {\displaystyle }{V}_{\mathrm{D}}=\left[u\cos \mathit{\phi }+v\sin \mathit{\phi }\right]\cos \mathit{\theta }+V_{\mathrm{D}\mathrm{\_}\mathrm{90}}\sin \mathit{\theta }\left[{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}\right],\\ {\displaystyle }& {\displaystyle }{\text{SW}}_{\mathrm{H}\mathrm{\_}\mathrm{90}}=\\ {\displaystyle }& {\displaystyle }\sqrt{{\displaystyle \frac{\begin{array}{c}\frac{\mathrm{4}{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{4}}{\left|K_{\mathrm{w}}\right|}^{\mathrm{2}}}{\sum }_{i=\mathrm{1}}^M\left[{\int }_{\mathrm{0}}^{\mathrm{\infty }}\mathit{\right\{}{\left|f_{bi}^{\left(\mathit{\pi }\right)}\right|}^{\mathrm{2}}\\ -\mathrm{2}{R}_e\left[f_{bi}^{\left(\mathit{\pi }\right)\ast }\left(f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\left(\mathit{\pi }\right)}\right)\right]A_{\mathrm{2}i}\\ +{\left|f_{bi}^{\left(\mathit{\pi }\right)}-f_{ai}^{\mathit{\pi }}\right|}^{\mathrm{2}}{A}_{\mathrm{4}i}\mathit{\left\}}{\left(V_{Fi}\left(D\right)-\frac{V{Z}_i}{Z_{\mathrm{hh}}}\right)}^{\mathrm{2}}{N}_i\left(D\right)\mathrm{d}D\right]\end{array}}{Z_{\mathrm{hh}}}}}\\ {\displaystyle }& {\displaystyle }\left[{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}\right],\\ {\displaystyle }& {\displaystyle }{\text{SW}}_{\text{TOT}}=\sqrt{{\text{SW}}_{\mathrm{H}}^{\mathrm{2}}+{\text{SW}}_{\mathrm{T}}^{\mathrm{2}}+{\text{SW}}_{\mathrm{S}}^{\mathrm{2}}+{\text{SW}}_{\mathrm{V}}^{\mathrm{2}}}\left[{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}\right],\end{array}$$

where SW_H, SW_T, SW_S, and SW_V are contributions from different hydrometeor terminal velocity, turbulence, mean wind shear, and cross wind. Tatarevic et al. (2019) describes detailed computations of these contributions.

The lidar observables are calculated as follows. For spherical droplets, using the BHMIE Mie code (Bohren and Huffman, 1998), the single particle extinction σ_α and backscattering cross sections σ_β are computed for a lidar wavelength (i.e., 905, 532, and 353 nm).

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\beta }}_{\text{true}}=\sum _i{\mathit{\sigma }}_{\mathit{\beta }}\left(D_i\right)N\left(D_i\right)\mathrm{\Delta }{D}_i\\ {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\text{ext}}=\sum _i{\mathit{\sigma }}_{\mathit{\alpha }}\left(D_i\right)N\left(D_i\right)\mathrm{\Delta }{D}_i\end{array}$$

The value of refractive index of water used in the calculations is $\mathrm{1.327}+i\mathrm{0.672}\times {\mathrm{10}}^{\mathrm{6}}$ (Hale and Querry, 1973). The attenuated backscatter, β_atten, at a distance z can be written as

$${\displaystyle }{\displaystyle }{\mathit{\beta }}_{\text{atten}}\left(z\right)=\underset{\mathrm{0}}{\overset{z}{\int }}{\mathit{\beta }}_{\text{true}}\left(z\right)\exp \left(-\mathrm{2}{\mathit{\alpha }}_{\text{ext}}\left(z\right)\right)\mathrm{d}z[\mathrm{sr}\mathrm{m}{]}^{-\mathrm{1}},$$

where β_true is the true backscatter at height z, and α_ext is the extinction coefficient:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\beta }}_{\text{true}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\sum _i{\mathit{\sigma }}_{\mathit{\beta }}\left(D_i\right)N\left(D_i\right)\mathrm{\Delta }{D}_i[\mathrm{sr}\mathrm{m}{]}^{-\mathrm{1}},\\ {\displaystyle }& {\displaystyle }{\mathit{\alpha }}_{\text{ext}}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\sum _i{\mathit{\sigma }}_{\mathit{\alpha }}\left(D_i\right)N\left(D_i\right)\mathrm{\Delta }{D}_i[\mathrm{sr}\mathrm{m}{]}^{-\mathrm{1}}.\end{array}$$

Appendix B: List of acronyms

ARM

Atmospheric Radiation Measurement

ARSCL

Active Remote Sensing of Clouds

CDF

Cumulative distribution function

CFP

Cloud fraction profile

CRM

Cloud-resolving model

CR-SIM

Cloud-resolving model Radar SIMulator

C-SAPR

C-band scanning ARM precipitation radar

CWRHI

Cross-wind range–height indicator

GCM

Global climate model

GNU GPL

GNU General Public License

KAZR

Ka-band ARM Zenith-pointing Radar

LASSO

LES ARM Symbiotic Simulation and Observation

LES

Large-eddy simulation

LUT

Look-up table

MAPE

Mean absolute percentage error

MC3E

Midlatitude Continental Convective Clouds Experiment

MCS

Mesoscale convective system

MF

Mass flux

MPL

Micro pulse lidar

OSSE

Observing system simulation experiment

PPI

Plan position indicator

PSD

Particle size distribution

RHI

Range–height indicator

RMSE

Root mean square error

SCR

Scanning cloud radar

SGP

Southern Great Plains

UF

Updraft fraction

VCP

Volume coverage pattern

WRF

Weather Research and Forecasting model

X-SAPR

X-band scanning ARM precipitation radar