the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
ZJUAERO V0.5: an Accurate and Efficient Radar Operator designed for CMAGFS/MESO with the capability to simulate nonspherical hydrometeors
Hejun Xie
In this study, we present a new forward polarimetric radar operator called the Accurate and Efficient Radar Operator designed by ZheJiang University (ZJUAERO). This operator was designed to interface with the numerical weather prediction (NWP) model of the global forecast system/regional mesoscale model of the China Meteorological Administration (CMAGFS/MESO). The main objective of developing this observation operator was to assimilate observations from the precipitation measurement radar (PMR). It is also capable of simulating the groundbased radar's polarimetric radar variables, excluding the Doppler variables such as radial velocity and spectrum width. To calculate the hydrometeor optical properties of ZJUAERO, we utilize the invariantimbedding Tmatrix (IITM) method, which can handle nonspherical and inhomogeneous hydrometeor particles in the atmosphere. The optical database of ZJUAERO was designed with a multilayered architecture to ensure the flexibility in hydrometeor morphology and orientation specifications while maintaining operational efficiency. Specifically, three levels of databases are created that store the singlescattering properties for different shapes at discrete sizes for various fixed orientations, integrated singlescattering properties over shapes and orientations, and bulkscattering properties incorporating the size average, respectively. In this work, we elaborate on the design concepts, physical basis, and hydrometeor specifications of ZJUAERO. Additionally, we present a case study demonstrating the application of ZJUAERO in simulating the radar observations of Typhoon Haishen.
 Article
(11180 KB)  Fulltext XML
 BibTeX
 EndNote
The development of regional models with finer horizontal resolutions, such as the Chinese operational regional numerical weather prediction (NWP) model, known as the regional mesoscale model of the China Meteorological Administration (CMAMESO) (Chen et al., 2008; Shen et al., 2023), necessitates the acquisition of more convectivescale information about the atmosphere to improve quantitative precipitation forecasts. Fortunately, the measurements from spaceborne and groundbased weather radars provide valuable sources of threedimensional, kilometerscale volume data with high temporal resolutions. However, weather radar can only observe the amplitude and phase of electromagnetic waves echoed from meteorological objects, specifically various types of hydrometeors. Except for the Doppler radar variables, such as radial velocity (beyond the scope of this study), it is challenging to establish a connection between the prognostic variables simulated by the NWP model and the observable polarimetric radar variables, which are inferred from the statistical moments of voltage time series collected by the receiver of the weather radar electronics system (Zhang, 2016).
The software package introduced in this work is referred to as a “forward radar operator”, designed to transform the model prognostic variables into the observation space, resulting in equivalent simulated synthetic radar variables. Utilizing a unified forward radar operator for assimilations and retrievals is believed to be superior to employing an ensemble of retrieval relationships along with preprocessing procedures and corrections for different frequencies and platforms. In essence, using radar data in the observation space is preferred over the model space due to the highly nonlinear and nonunique nature of the processes that observational operator of polarimetric radar describes.
Several radar operators have been developed and published over the past several decades. For instance, Jung et al. (2008) implemented a polarimetric radar simulator known as the Polarimetric Radar data Simulator developed by the Center for Analysis and Prediction of Storms (CAPSPRS) at the University of Oklahoma. This simulator uses spheroids to characterize hydrometeors and computes optical properties using online Rayleigh approximations or offline lookup tables (LUTs) constructed by the extended boundary condition method (EBCM) described in Mishchenko and Travis (1994). This simulator has been applied to lowfrequency bands, such as the S, C, and X bands. Ryzhkov et al. (2011) described another radar operator for research purposes, specifically tailored for spectraresolving cloud microphysics models, although it is more computationally expensive. Zeng et al. (2016) described an efficient radar operator that is onlinecoupled to the Consortium for Smallscale Modeling (COSMO) and Icosahedral Nonhydrostatic (ICON) weather and climate model, makes use of Mie/Tmatrix scattering lookup tables of solid, liquid, and melting (mixedphased) hydrometeors, and is named the Efficient Modular VOlume scanning RADar Operator (EMVORADO). While early versions of EMVORADO focus on nonpolarimetric radar variables, later developments on EMVORADO have enabled its capability to simulate dualpolarization variables and conduct sufficient evaluations (Trömel et al., 2021; Shrestha et al., 2022). In addition to the above operators, Wolfensberger and Berne (2018) reported a crossplatform polarimetric radar operator termed a POLarimetric forward radar operator for the COSMO model (COSMOPOL). This operator was designed for the COSMONWP model and can simulate melting particles. The optical database of the COSMOPOL was constructed also using the EBCM, characterizing all hydrometeor particles as homogeneous spheroids. However, in the COSMOPOL, the hydrometeor orientations and shape parameter settings are fixed at the level of the optical database, which limits sensitivity testing and finetuning. Wang and Liu (2019) reported a forward reflectivity observation operator (together with its tangent linear and adjoint operator) with the simulation capability of frozen hydrometeors designed to assimilate the data of groundbased radar in the Weather Research and Forecasting Model (WRF) in which the simulated reflectivity is parameterized as fastpolynomial relationships with respect to the mixing ratios of hydrometeors. More recently, Oue et al. (2020) developed the Cloudresolving model Radar SIMulator (CRSIM), which can simulate polarimetric radar and light detection and ranging (lidar) observations for various cloudresolving models (CRMs), including the WRF, ICON, the Regional Atmospheric Modeling System (RAMS), and the advective statistical forecast model (SAM). CRSIM has a unique capability for explicitly representing the rimming procedure by coupling it with the prognostic variable known as the rimming ratio in the Predicted Particle Properties (P3) microphysics package (Morrison and Milbrandt, 2015). However, CRSIM is currently limited to groundbased platforms and offers no explicit treatment for melting particles. Moreover, the fastparameterized forward radar operator developed by Zhang et al. (2021) contains a melting scheme module for targeting data assimilation purposes.
This work aims to design a crossband and crossplatform radar operator for research purposes, such as microphysics package validation, and operational data assimilation use in the global forecast system/regional mesoscale model of the China Meteorological Administration (CMAGFS/MESO). The software prototype of this radar operator, hereafter referred to as the Accurate and Efficient Radar Operator designed by ZheJiang University (ZJUAERO), essentially addresses the scattering computation of hydrometeors and the construction of an optical properties database as the key aspects in the evolution of the radar operator. We utilize a semianalytical scattering computation approach of a Tmatrix to ensure accuracy and feature a multilayered optical database that includes singlescattering properties at discrete sizes and orientations, integrated singlescattering properties over shapes and orientations, and bulkscattering properties incorporating the size average. Additionally, ZJUAERO allows flexibility in the particle orientation and shape probability distribution tuning. Notably, this software has also inherited established techniques, such as subbeam sampling, used for simulating the effects of beambending, beambroadening, or beamshielding (Ryzhkov, 2007).
The development of ZJUAERO was primarily motivated by the future data assimilation purpose of the precipitation measurement radar (PMR) aboard the FengYun3 Rain Measurement satellite (FY3 RM) (Zhang et al., 2023). The parameters of the instrument FY3 RM/PMR are comparable with those of the Global Precipitation Measurement/Dualfrequency Precipitation Radar (GPM/DPR). The DPR aboard the GPM was designed to obtain the storm structure, rainfall rates, drop size distributions (DSDs), pathintegrated attenuation, and other useful information that is not available from passive sensor observations (Iguchi et al., 2003; Iguchi, 2020). Both PMR and DPR have two bands (the Ku and Ka bands), while PMR is designed with a swath of 303 km, horizontal resolution of 5 km at nadir, and radial resolution of 250 m.
This paper is organized as follows. Section 2 provides readers with an overview of the general concepts of ZJUAERO. The variables (matrices) that describe the scattering properties of hydrometeors are presented here to clarify the notation convention used in this study and eliminate ambiguity in the context. Details of the implementation of hydrometeor settings, including dielectric models, aspect ratio models, orientation preference models, and particle size distributions (PSDs) are also listed in Sect. 2. Section 3 elaborates on the flexible architecture of the stateoftheart nonspherical scattering property database, which distinguishes the ZJUAERO from its predecessors. The characteristics of the multilayered optical database are illustrated using a nonspheroid hydrometeor model, namely the Chebyshevshaped raindrop. Moreover, Sect. 4 presents a case study of observations of a tropical cyclone given by the spaceborne radar GPM/DPR and compared with simulations given by the ZJUAERO. Sensitivity tests of PSDs and nonsphericity are also described in this section. Section 5 summarizes this study and describes the development plans for the ZJUAERO.
The ZJUAERO was developed to simulate the observable polarimetric radar variables for radar systems aboard on various platforms, including groundbased, spaceborne, and potentially airborne radars in the future. The conceptual graph of the multiplatform radar operator is shown in Fig. 1. While the physical principles of the weather radar detection process are universal, certain factors, such as beambroadening, beambending, and beamblocking among others (as indicated along the beam trajectory in Fig. 1), which are critical to groundbased radars are not equally important across platforms. For example, the beambending effect is typically negligible for spaceborne radar due to large absolute elevation angles (usually 70–90°, as illustrated in Fig. 1 for spaceborne radar) and shorter beam trajectories (usually <20 km) below the model top, as compared to the groundbased radar (the trajectory can reach up to about 250 km). However, the simulations of spaceborne and groundbased radars share consistent hydrometeor setting entries for snow/graupel, melting snow/graupel, and rain, such as the dielectric model, particle morphology (distribution), size distribution, orientation preference, and others.
2.1 Flowchart and concepts
Figure 2 provides an overview of the ZJUAERO simulation procedure for a singleradar scan. ZJUAERO consists of five modules represented by green boxes in Fig. 2. These modules are the (A) NWP interface submodule, (B) beam submodule (decomposed into B1 and B2 for groundbased and spaceborne radars, respectively, (C) interpolation submodules, (D) hydrometeor submodule, and (E) core submodule. The workflow of ZJUAERO can be summarized as follows:
 1.
The NWP interface submodule (A) reads NWP prognostic variables from external storage files and interpolates the data defined on the original model grid (such as horizontal ArakawaC grid and vertical Charney–Phillips grid in CMA–GFS/MESO; Chen et al., 2008; Shen et al., 2023) to a regular grid on which all variables are colocated and evenly spaced horizontally (in the space of projection). The prognostic variables include hydrometeor variables (Q_{c}, Q_{i}, Q_{r}, Q_{s}, and Q_{g}, representing the mixing ratios of cloud water, cloud ice, rain, snow, and graupel, respectively) and dynamic variables (U, V, W, T, P, and Q, representing the zonal, meridional, vertical wind, temperature, pressure, and water vapor mixing ratio, respectively). Additionally, static model information such as orography data defined on model grids is required for simulating partial beam blocking. This step prepares for a quick and convenient secondround interpolation from the regular model grid to the radar beam trajectories (in step 3). It is worth pointing out that ZJUAERO can also interface with the output of WRF NWP model (Skamarock et al., 2019). Enabling ZJUAERO to interface with grid data from another NWP model involves only technical adjustments that require basic information about that NWP model's horizontal mesh (projections), vertical grid, and variable mapping table.
 2.
The beam trajectory submodule (B) calculates the propagation trajectories of the radar beam. For groundbased radar (B1), users have the option of using an online trajectory solver that uses temperature and humidity profiles above the radar site. Specifically, the atmosphere refractive index Na are determined from atmosphere temperature T, pressure P, and water vapor mixing ratio Q. The trajectory is determined by a raytracing ordinary differential equation (ODE) solver (Zeng et al., 2014). Additionally, ZJUAERO offers an alternative option of using an offline $\mathrm{4}/\mathrm{3}{R}_{\mathrm{E}}$ solver for groundbased radar in ZJUAERO. Multiple subtrajectories (determined by the horizontal quadrature number N × vertical quadrature number M) are sampled within the 3 dB beamwidth of the main lobe for the N × M subtrajectories. The observable radar variables are calculated by integrating bulkscattering properties over the antenna patterns (as described in step 5) to obtain the final results (beambroadening and beamblocking are considered in this way). We applied the subbeam sampling and averaging methods described in Zeng et al. (2016). Nevertheless, for spaceborne radar (B2), the beam trajectory is computed using the geometry shown in Fig. A1, and subtrajectory averaging is not supported at this time. For more details on the beam trajectory module of spaceborne radar (B2), please refer to Appendix A.
 3.
The interpolation submodule (C) uses a trilinear interpolation algorithm to interpolate the NWP prognostic variables to the gates of radar beam trajectories. The quick trilinear interpolations are performed in a twostep manner: (1) vertically interpolating the data defined on the eight vertices of the cube containing the radar gates to the four vertices of the horizontal box surrounding the radar gate and (2) performing bilinear interpolation; we adopted the approach described in Appendix A of Wolfensberger and Berne (2018) and reimplemented it as a Cython extension.
 4.
The hydrometeor submodule (D) specifies the properties of hydrometeors in each radar gate along each trajectory, usually by loading presets of microphysicsconsistent constants of hydrometeors. This includes the orientation preference, probability distribution of particle morphology, particle size distribution (PSD), and other parameters. Those presets of hydrometeor properties can also be modified by users to perform “forced” (i.e., inconsistent with microphysics) simulations for research purposes. The PSD parameters of hydrometeors are solved in this step from prognostic bulkhydrometeor variables in NWP models (mass concentrations and number concentrations). For more details, see Sect. 2.3.
 5.
The core module (E) finally calculates the polarimetric radar variables.
 5.1
The bulkscattering properties of particles in each radar gate are computed by integrating the singlescattering properties across the size spectrum of each hydrometeor type and summing over hydrometeor types, which can be conducted either online (E1; research mode) or offline (E2; operational mode) in ZJUAERO. The scattering property LUTs of ZJUAERO are consulted in this step, which is composed of three levels (level A, level B, and level C). The multilayered architecture will be introduced in detail in Sect. 3.1. The research mode is more flexible since it accesses the level B database for singlescattering property, while the operational mode is more efficient by accessing the level C database for bulkscattering properties straightforward. We provide users with tool scripts for level A to B and level B to C conversions (integration parameters are alterable in YAML configure files).
 5.2
Once the bulkscattering properties on each subtrajectory grid point within each radar gate are available, the antenna pattern integration involves integrating the bulkscattering properties within the scanning volume of each radar gate.
 5.3
Then, the core module calculates the intrinsic polarimetric radar variables on each radar gate, based on the bulkscattering properties incorporating the antenna pattern integration presented in step 5.2. These radar variables include singlepolarization reflectivities (Z_{H} for horizontal reflectivities), dualpolarization variables (Z_{DR}, K_{DP}, δ_{hv}, Φ_{DP}, and ρ_{hv} for differential reflectivity, differential phase shift, backscatter differential phase, and copolar correlation coefficient, respectively), and attenuation variables (a_{H} and a_{V} for horizontal and vertical attenuation coefficients, respectively). The definitions of these variables can be found in Appendix D. For a detailed explanation of the intrinsic radar variables, please refer to Zhang (2016).
 5.4
In the final step of core module, the observable radar variables are obtained by taking into account the attenuation and phase shift accumulated along the beam trajectories.
 5.1
The above procedures of steps 2–5 are independent for each single beam, which guarantees the toplevel scalability of the forward operator. Therefore, we are using the sharedmemory Python parallel library (multiprocessing) to speed up the simulation.
In addition, for those computationally intensive tasks (such as the trilinear interpolation in step 3), we are employing the technique of mixed programming (building C/Fortran extensions that interface with Python scripts) to further accelerate the computation.
The performance of the forward operator is generally satisfactory; ZJUAERO can complete a groundbased station volume scan with nine plan position indicator (PPI) sweeps in 2 min on a modern laptop CPU with a sixcore processor (i710750H) if the online size distribution integration is performed and the operator takes the level B singlescattering property database as input. If the level C bulkscattering property database is used, such a volume PPI scan can be completed even faster (in 30 s). Such efficiency can support data assimilation purposes while also preserving flexibility for research purposes.
In this paper, we do not elaborate on the algorithm details of certain issues, such as (a) trilinear interpolation, (b) subbeam sampling and antennapatternweighted averaging, and (c) a raytracing trajectory solver. For trilinear interpolation, we follow the approach described in Wolfensberger and Berne (2018) to interpolate the model grid data to the radar gates of each subtrajectory grid point. Regarding subbeam sampling and antennapatternweighted averaging, we utilized Gauss–Hermite quadrature, as outlined in Caumont et al. (2006). As the raytracing trajectory solver, we offer users both an offline beam trajectory solver ($\mathrm{4}/\mathrm{3}{R}_{\mathrm{E}}$) and an online beam trajectory solver (Zeng et al., 2014). All of these methods are reimplemented in Python using either efficient NumPy/SciPy API or customized Cython extensions.
2.2 Physical basis
At this point, we can specify the formulation convention used in the radar operator and move on to the nonspherical optical database design of ZJUAERO.
The amplitudescattering matrix S_{FSA} of a single particle is defined as follows:
In Eq. (1), E indicates the vector electric field, and the superscripts “sca” and “inc” represent scattering and incident waves, respectively. The subscripts “h” and “v” indicate two decomposed components of the vector electric field in the horizontal and vertical directions, respectively. Specifically, the horizontal and vertical unit vectors are defined as unit vectors, $\widehat{\mathit{\varphi}}$ and $\widehat{\mathit{\theta}}$, of the spherical coordinate system under the forward scattering alignment (FSA) convention. The differences between FSA and backwardscattering alignment (BSA) are described in Appendix C. k_{0} is the wave number in free space, and r is the distance from the particle center.
The scattering matrix relates the incident electric field to the scattered electric field, and it must have the dimension of L (L is the dimension symbol of length). In practice, the amplitudescattering matrix is usually expressed in the units of millimeters. The amplitudescattering matrix of a single hydrometeor particle is obtained from scattering computations, specifically using the Tmatrix method in this study.
Apart from the 2 × 2 complex amplitudescattering matrix S defined on complex electric field vector bases, one can define the 4 × 4 real matrix, known as the Mueller matrix, Z, and extinction matrix, K, which describe polarimetric lightscattering and extinction properties of particles on the real Stokes vector bases, respectively. We kept the definitions of Z and K consistent with those mentioned in a study by Mishchenko (2014).
The Mueller and extinction matrices can be derived from the amplitudescattering matrix (Mishchenko, 2014). Specifically, the forward amplitudescattering matrix shows a linear relationship with K. For example, the formulas of the matrix elements of K used in this study are shown below:
Here, the notations of ℑ and ℜ indicate the real and imaginary parts of a complex number, respectively. On the other hand, the backscatteringamplitude matrix can be used to calculate Z in backscattering geometry. For example, the formulas of the matrix elements of Z used in this study are shown below:
The complete set of equations from the elements of the amplitude matrix S to each element of the Mueller matrix Z or the extinction matrix K can be found in Mishchenko (2014). The elements of Z and K are both in the dimension of L^{2} (i.e., they are usually expressed in the units of mm^{2}).
We could compute the bulkscattering properties 〈X〉^{1} used in the expression of several polarimetric radar variables by performing size distribution integrations over elements of Z and K matrices:
Here, N(D) represents the number concentration distribution function (in units of mm^{−1} m^{−3}) over the particle spectrum, and D is the diameter^{2} of the hydrometeor (in units of mm). The elements of bulk matrices 〈Z〉 and 〈K〉 were usually expressed in units (of mm^{2} m^{−3}).
Note that we can only apply the particle ensemble integration over the complex amplitudescattering matrix, S, in the forward scattering geometry. This is because integrating over extinction matrix elements is proportional to integrating corresponding forward amplitudescattering matrix elements. For example:
Equation (5) is derived by performing ensemble mean on Eq. (2).
Therefore, we use Mueller and extinction matrices to represent radar variables because ensemble means can be performed directly on them, as is not the case for amplitude matrix. Also, they have a unified dimension of L^{2}. The LUTs in ZJUAERO store Mueller and extinction matrices instead of the amplitude matrix (see Appendix B, which will be further described in Sect. 3.1). The equations for radar variables based on Mueller and extinction matrices can be found in Appendix D.
2.3 Hydrometeor specifications
Table 1 summarizes the hydrometeor specifications in ZJUAERO, with the following important notes:

During the earlydevelopment stage of ZJUAERO, we initially used the dielectric model for water proposed by Liebe et al. (1991). However, we transitioned to a more accurate and contemporary dielectric model developed by Ellison (2007). Nevertheless, we retained the outdated option and optical property lookup table from the old dielectric model in our archive for future reference and comparison (see the column of the refractive index model in Table 1).

In principle, it is encouraged to use PSD schemes compatible with the microphysics package in the NWP model to ensure consistent hydrometeor settings in simulations. Therefore, ZJUAERO, which is designed for CMA–GFS/MESO, provides PSD options that are compatible with the singlemoment microphysics scheme WSM6 (Hong and Lim, 2006) and Thompson (Thompson et al., 2008), which are widely used in the global and regional operational models of CMA. For instance, option C1 for rain and option B1 for graupel are compatible with the WSM6 package, while option C2 for rain and option B2 for graupel are compatible with the Thompson package (see the column of particle size distribution in Table 1). However, for the snow category, we implemented the PSD scheme of Field et al. (2005) as the only option since it is the widely acknowledged as the best globally applicable temperaturedependent PSD model for solid precipitation. Additionally, we have provided users with some additional PSD schemes for sensitivity assessment, such as options C3 (Wang et al., 2016) and C4 (Abel and Boutle, 2012) for the rain category. Those PSD schemes in the ZJUAERO that are incompatible with the microphysics package used in the NWP model are referred to as forced PSD schemes.

Solid hydrometeor categories, such as snow and graupel, are known to have relatively larger uncertainties associated with parameterizations of aspect ratios and orientation preference. To address these uncertainties, a field research study conducted by Garrett et al. (2012) used an in situ observation instrument called the multiangle snowflake camera (MASC). This instrument was used to measure the aspect ratios and orientations of over 30 000 solid particles in the eastern Swiss Alps. The particles were then classified into aggregates (corresponding to the snow category in this study), rimed particles, and graupel, as described in Garrett et al. (2015). The fitted model from this study was used as a priori knowledge of hydrometeor shape specifications in the ZJUAERO (see the column of shape parameters in Table 1):
$$\begin{array}{}\text{(6)}& p(\mathit{\gamma};{D}_{\mathrm{max}})\propto {\left(\mathit{\gamma}\mathrm{1}\right)}^{K\left({D}_{\mathrm{max}}\right)\mathrm{1}}\mathrm{exp}\left({\displaystyle \frac{\mathit{\gamma}\mathrm{1}}{\mathrm{\Theta}\left({D}_{\mathrm{max}}\right)}}\right).\end{array}$$Equation (6) provides the probability distribution function of the aspect ratio in which γ is the reciprocal of the aspect ratio (minor axis/major axis, always less than unity). It is assumed to follow a gamma distribution with an offset coefficient of 1 (i.e., γ>1). The shape coefficient, K(D_{max}), and a scale coefficient, Θ(D_{max}), are functions of the particle maximum dimension. We used the power law relationships of K(D_{max}) and Θ(D_{max}) that were fitted by Wolfensberger and Berne (2018) as follows:
$$\begin{array}{}\text{(7)}& \begin{array}{l}{K}_{\mathrm{snow}}\left({D}_{\mathrm{max}}\right)=\mathrm{8.42}{D}_{\mathrm{max}}^{\mathrm{0.57}};{\mathrm{\Theta}}_{\mathrm{snow}}\left({D}_{\mathrm{max}}\right)\\ \phantom{\rule{1em}{0ex}}=\mathrm{0.053}{D}_{\mathrm{max}}^{\mathrm{0.79}}.\\ {K}_{\mathrm{graupel}}\left({D}_{\mathrm{max}}\right)=\mathrm{3.2}{D}_{\mathrm{max}}^{\mathrm{0.42}};{\mathrm{\Theta}}_{\mathrm{graupel}}\left({D}_{\mathrm{max}}\right)\\ \phantom{\rule{1em}{0ex}}=\mathrm{0.074}{D}_{\mathrm{max}}^{\mathrm{0.67}}.\end{array}\end{array}$$ 
Since solving the PSD parameters from NWP prognostic hydrometeor mass concentrations (for bulk microphysics) requires knowledge of the mass of various particle sizes, the mass–diameter relationship is crucial in determining the PSD (see the column of mass–diameter relationship in Table 1). In this study, all the hydrometeor categories follow the gamma distribution (the widely used exponential distribution is just a special case of gamma distribution when μ=1),
$$\begin{array}{}\text{(8)}& N\left(D\right)={N}_{\mathrm{0}}{D}^{\mathit{\mu}}{e}^{\mathrm{\Lambda}D},\end{array}$$in which N_{0} is the intercept, Λ is the slope, and μ is the shape coefficient of that gamma distribution. As is often the case (for all PSD options in ZJUAERO; see Table 1), μ is a prescribed constant in the microphysics package, while N_{0} either equals a prescribed constant or relates to Λ by a power law in which x_{1} and x_{2} are parameters fitted by drop size distribution (DSD) observations (see Sect. 3.4):
$$\begin{array}{}\text{(9)}& {N}_{\mathrm{0}}={x}_{\mathrm{1}}{\mathrm{\Lambda}}^{{x}_{\mathrm{2}}}.\end{array}$$If a hydrometeor category is of a singlemoment microphysics scheme, given the mass concentration of arbitrary hydrometeor category “x” Q_{x} (in units of kg m^{−3}),
$$\begin{array}{}\text{(10)}& {Q}_{x}=\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}\stackrel{\mathrm{\u203e}}{m}\left(D\right)\cdot N\left(D\right)\mathrm{d}D.\end{array}$$If the mass–diameter relationship can be approximated as a power law form $\stackrel{\mathrm{\u203e}}{m}\left(D\right)=a{D}^{b}$ in which $\stackrel{\mathrm{\u203e}}{m}\left(D\right)$ is the average mass of a given size of a particle (considering that some hydrometeor categories have a probability distribution over shape parameters such as aspect ratio). Then we can solve the unknown parameter Λ and determine all relevant PSD information that pertains to that radar gate analytically by plugging Eq. (8) into Eq. (10) as follows:
$$\begin{array}{}\text{(11)}& \mathrm{\Lambda}={\left[{\displaystyle \frac{a{x}_{\mathrm{1}}\mathrm{\Gamma}(\mathrm{1}+b+\mathit{\mu})}{{Q}_{x}}}\right]}^{{\scriptscriptstyle \frac{\mathrm{1}}{\mathrm{1}+b+\mathit{\mu}{x}_{\mathrm{2}}}}}.\end{array}$$Again, there is a microphysicsconsistent mass–diameter relationship $\stackrel{\mathrm{\u203e}}{m}\left({D}_{\mathrm{max}}\right)=\frac{\mathit{\pi}}{\mathrm{6}}{\mathit{\rho}}_{\mathrm{sp}}{D}_{\mathrm{max}}^{\mathrm{3}}({\mathit{\rho}}_{\mathrm{sp}}$ is the overall density of the solid precipitation sphere particle) for snow and graupel. Many microphysics schemes, such as WSM6, simply treated solid precipitation categories as spheres with different ice–air mixture ratios and hence different overall densities ρ_{sp}. However, this practice can result in a huge inconsistency between the mass concentration represented by radar operators and the microphysics schemes since the actual average mass of a given size bin is represented by the following probabilityweighted averaging over the aspect ratio for solid hydrometeors shown below:
$$\begin{array}{}\text{(12)}& \stackrel{\mathrm{\u203e}}{m}\left({D}_{\mathrm{max}}\right)=\underset{\mathrm{1}}{\overset{\mathrm{\infty}}{\int}}p(\mathit{\gamma};{D}_{\mathrm{max}})m(\mathit{\gamma};{D}_{\mathrm{max}})\mathrm{d}\mathit{\gamma}.\end{array}$$It turns out that fitting Eq. (12) into power laws is troublesome when the probability distribution function p(γ;D_{max}) varied dramatically with diameter. This problem will become even worse when we introduce other nonspherical particles, such as snowflakes.
To resolve this error when using a traditional mass–diameter relationship, we implemented a benchmark PSD solver employing a numerical method (i.e., Newton iteration) for particles with complicated morphology specifications. The simplest case, a PSD represented by an exponential distribution with a fixed intercept parameter N_{0}, can be used as an example. Here the equation for the unknown PSD parameter Λ can be expressed as follows:
$$\begin{array}{}\text{(13)}& \begin{array}{rl}{Q}_{x}& =\sum _{{D}_{i}={D}_{\mathrm{l}}}^{{D}_{i}={D}_{\mathrm{u}}}\stackrel{\mathrm{\u203e}}{m}\left({D}_{i}\right)N\left({D}_{i}\right)\mathrm{\Delta}D\\ & =\sum _{{D}_{i}={D}_{\mathrm{l}}}^{{D}_{i}={D}_{\mathrm{u}}}\stackrel{\mathrm{\u203e}}{m}\left({D}_{i}\right){N}_{\mathrm{0}}{e}^{\mathrm{\Lambda}{D}_{i}}\mathrm{\Delta}D,\end{array}\end{array}$$in which the integration in Eq. (10) is truncated within a specific diameter range of [D_{l},D_{u}] and further discretized as a summation. The expression of the exponential distribution is then substituted in for the second equality. $\stackrel{\mathrm{\u203e}}{m}\left({D}_{i}\right)$ at discretized size bins is precomputed by Eq. (12) and treated as a constant. The mass m(γ;D_{i}) for each morphology specification can be calculated as the density of the pure ice ρ_{ice} multiplied by the volume occupied by the nonspherical hydrometeor particle model V(γ;D_{i}) (calculated from mathematical formulas of geometrical bodies).
We then define the function F(Λ) and its derivative F^{′}(Λ) in the Newton iteration as follows:
$$\begin{array}{}\text{(14)}& \begin{array}{l}F\left(\mathrm{\Lambda}\right)=\sum _{{D}_{i}={D}_{\mathrm{l}}}^{{D}_{i}={D}_{\mathrm{u}}}\stackrel{\mathrm{\u203e}}{m}\left({D}_{i}\right){N}_{\mathrm{0}}{e}^{\mathrm{\Lambda}{D}_{i}}\mathrm{\Delta}D{Q}_{x},\\ {F}^{\prime}\left(\mathrm{\Lambda}\right)=\sum _{{D}_{i}={D}_{\mathrm{l}}}^{{D}_{i}={D}_{\mathrm{u}}}\stackrel{\mathrm{\u203e}}{m}\left({D}_{i}\right){N}_{\mathrm{0}}{e}^{\mathrm{\Lambda}{D}_{i}}{D}_{i}\mathrm{\Delta}D.\end{array}\end{array}$$Based on the above formulation, the iteration relationship to derive the (n+1)th guess Λ_{n+1} from the nth guess Λ_{n} can be expressed as shown below:
$$\begin{array}{}\text{(15)}& {\mathrm{\Lambda}}_{n+\mathrm{1}}={\mathrm{\Lambda}}_{n}{\displaystyle \frac{F\left({\mathrm{\Lambda}}_{n}\right)}{{F}^{\prime}\left({\mathrm{\Lambda}}_{n}\right)}}.\end{array}$$While performing iterations online (the benchmark PSD solver) may lead to a decrease in the efficiency of the ZJUAERO, this problem can be resolved using bulkscattering property (BSP) LUTs instead of singlescattering property (SSP) LUTs (see Sect. 3.5).
In Sect. 2, we provide a comprehensive description of the design of ZJUAERO. Specifically, we emphasize that the hydrometeor optical property database includes the elements of Z and K (in units of mm^{2}) for singlescattering properties and those of 〈Z〉 and 〈K〉 (in units of mm^{2} m^{−3}) for bulkscattering properties, both in the FSA convention. In the first two subsections, we will delve into the design of the database (LUT) in ZJUAERO with more details.
Furthermore, ZJUAERO currently encompasses three types of hydrometeors: (1) rain, (2) snow, and (3) graupel. Among these hydrometeors, the rain category offers a nonspherical shape parameter option, known as the Chebyshev shape. This shape differs from the traditional spheroid shape commonly used in other radar observation operators. Therefore, we will evaluate the database contents using the Chebyshev raindrop as an example in the last three subsections.
3.1 The multilayered architecture
The ZJUAERO optical property database is designed with a multilayered architecture consisting of three levels: (1) the raw singlescattering property (RSSP) database (level A), (2) the applied singlescattering property (ASSP) database (level B), and (3) the BSP database (level C). These levels are described in detail in Table 2.
The RSSP level A database contains the optical properties of individual particles without any averaging or integration over the shape parameters and orientations, which normally consumed a significant amount of memory resources (∼ 10^{1} GB). However, using the RSSP in ZJUAERO would require the online integration of orientations and shape parameters, leading to a significant slowdown in radar operator performance. On the other hand, if shape and orientation averaging were applied during the creation of the database, and the raw optical data (RSSP) are discarded, then the resulting database would lack the flexibility needed for modifying the shape and orientation distributions. Additionally, future enhancements to the ZJUAERO database may involve incorporating more sophisticated shape parameters for the nonspherical hydrometeors. Therefore, retaining the RSSP level A database is essential to accommodate uncertainties and increase the convenience in experiments and simulations related to shape and orientation parameters.
Building on the RSSP level A database, the ASSP level B database improves the computational efficiency by carrying out averaging or orientation over shape parameters and orientations offline. Finally, the BSP level C optical database integrates the optical properties stored in the ASSP level B database over PSD and enables fast batch running in ZJUAERO for operational use. In summary, the multilayered architecture of the lookup table (stored in netCDF4 format files separately for different database levels) is to strike a balance between the flexibility of the database and the computational efficiency required by radar operators. Future releases of ZJUAERO will provide software tools for flexible conversions between the database levels, allowing users to easily configure integration parameters.
3.2 Scattering geometry
Figure 3 depicts the Cartesian coordinate system used to determine the scattering geometry of a hexagonal plate particle. The laboratory coordinate system, denoted as OX_{L}Y_{L}Z_{L}, is aligned by vertically positioning its OZ_{L} axis and placing its OX_{L} axis in the vertical plane specified by the incident radar beam and OZ_{L}. This alignment sets the azimuthal angle ϕ_{inc} of the radar beam to zero. In Fig. 3a, the direction of the radar beam is shown, which is practically specified using the elevation angle $e\in \left[\mathit{\pi}/\mathrm{2},\mathit{\pi}/\mathrm{2}\right]$. For groundbased radar, this angle is positive, while for spaceborne radar, it is negative. In this context, the polar angle of the radar beam, denoted as θ_{inc}, is related to the elevation angle through ${\mathit{\theta}}_{\mathrm{inc}}=\mathit{\pi}/\mathrm{2}e$. The shape of the particle is defined in the particle coordinate system OX_{P}Y_{P}Z_{P}. In the specific example shown in Fig. 3d, the OZ_{P} axis is set perpendicular to the basal face of a hexagonal plate. The origin “O” is placed at the mass center of the particle, and the OY_{P} axis intersected two opposite vertices of the hexagonal basal face. Using the ZYZ convention of the Euler angles specified by α, β, and γ (rotations were performed with respect to the OZ_{P}, OY_{P}, and OZ_{P} axes, respectively), any arbitrary orientation of a given particle can uniquely be determined (see Fig. 3b–d).
With a specified scattering geometry, we can now outline the procedures for computing the scattering properties of particles with fixed orientations (steps 1–3) and then integrate over them for scattering properties with specific orientation preference (step 4):

We used the Tmatrix method to compute the Tmatrix only once for a given particle and a radar beam wavelength. The guidelines for selecting a scattering computation approach are as follows: for both axially symmetric (i.e., the dielectric constant distribution of electromagnetic medium in spherical coordinates ε(rθφ) is irrelevant with azimuth angle φ) and homogenous particles, we used the EBCM Tmatrix code (Mishchenko and Travis, 1994), while for particles without axial symmetry or those that are inhomogeneous (or shapes for which the EBCM approach suffers from numerical stability issues), we applied the invariantimbedding Tmatrix (IITM) code (Bi et al., 2013; Bi and Yang, 2014; Bi et al., 2022; Wang et al., 2023).

With the Tmatrix computed in Step 1, we efficiently calculated the forward and backwardamplitudescattering matrix S_{FSA}in Eq. (1) for tuples of scattering geometry (α, β, γ, and e), using the method of Mishchenko (2000).

The forward and backwardamplitudescattering matrices S_{FSA} were then converted into the backscattering Mueller and extinction matrices, Z and K, respectively (Mishchenko, 2014).

When averaging over α and γ, we considered that the atmosphere is generally horizontally isotropic on the scale of hydrometeor particles. It is believed that they should have no preference for Euler angles α and γ, except for extreme conditions such as the lighteninginduced reorientation of ice crystals (Hubbert et al., 2010). Therefore, we performed internal averaging over α and γ for elements of Z and K in the scatteringcomputation code before generating the RSSP level A database:
$$\begin{array}{}\text{(16)}& \stackrel{\mathrm{\u203e}}{\mathbf{X}}(e,\mathit{\beta})={\displaystyle \frac{\mathrm{1}}{\mathrm{4}{\mathit{\pi}}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi}}{\int}}\mathrm{d}\mathit{\alpha}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi}}{\int}}\mathrm{d}\mathit{\gamma}\mathbf{X}(e,\mathit{\alpha},\mathit{\beta},\mathit{\gamma}),\mathbf{X}=\mathbf{Z},\mathbf{K}.\end{array}$$^{3}Hence, the level A database only has two residual scattering geometry dimensions: (1) e and (2) β, allowing for a reasonable volume for archiving.
The integration in the level A to level B database conversion tool can be formalized as the integration over the canting angle β as follows:
$$\begin{array}{}\text{(17)}& \stackrel{\mathrm{\u203e}}{\stackrel{\mathrm{\u203e}}{\mathit{X}}}\left(e\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\pi}}}\underset{\mathrm{0}}{\overset{\mathit{\pi}}{\int}}\mathrm{sin}\left(\mathit{\beta}\right)\mathrm{d}\mathit{\beta}p\left(\mathit{\beta}\right)\stackrel{\mathrm{\u203e}}{\mathit{X}}(e,\mathit{\beta}),\mathit{X}=\mathbf{Z},\mathbf{K}.\end{array}$$^{4}Here p(β) is the probability distribution of canting angle β. A Gaussian distribution in the polar angle is often used to approximate the probability distribution of canting angle,
$$\begin{array}{}\text{(18)}& p\left(\mathit{\beta}\right)=\mathrm{exp}\left({\displaystyle \frac{\mathit{\beta}}{\mathrm{2}{\mathit{\sigma}}_{\mathit{\beta}}^{\mathrm{2}}}}\right).\end{array}$$Here the standard deviation of canting angle, σ_{β}, can be found in the orientation preference column in Table 1.
Particle symmetry can be considered in the orientation averaging of Eqs. (16) and (17) to avoid a redundant evaluation of Z and K. For example, in the case of a hexagonal plate with an 6fold azimuthal symmetry and xyplane reflectance symmetry, Z and K only need to be evaluated and averaged over $\mathit{\alpha}\in \left[\mathrm{0},\mathit{\pi}/\mathrm{6}\right]$ and $\mathit{\beta}\in \left[\mathrm{0},\mathit{\pi}/\mathrm{2}\right]$.
3.3 Raw singlescattering properties: level A
To improve the accuracy of raindrop shape representation in the forward radar operator ZJUAERO, it is essential to apply an established raindrop model that accounts for various effects such as the surface tension, hydrostatic and aerodynamic pressures, and static electric forces. By incorporating these factors, a more accurate geometric representation of observed raindrop shapes can be achieved compared to the most commonly used spheroid model, as proved by photographs of water drops falling at the terminal velocity in air (Pruppacher et al., 1998). In a study conducted by Chuang and Beard (1990), an equilibrium differential equation was utilized to iteratively determine the shape of a falling raindrop with a given mass which corresponds to a specific diameter D_{eq}. The obtained results were then fitted to Chebyshev polynomials as follows (Chuang and Beard, 1990):
This equation involves the distance between the raindrop's surface to its mass center, denoted as χ(θ), where θ is the zenith angle in spherical coordinates (see Fig. 4d). Chuang and Beard (1990) provided Chebyshev coefficients a_{n}(D_{eq}) on diameter grids. Therefore, for a given D_{eq} size, the corresponding Chebyshev coefficients can be obtained through interpolation.
The Chebyshev model of raindrops deviates from xyplane reflectance symmetry, resulting in a flatter bottom surface compared to the spheroid. Conversely, the top surface of the raindrop described by the Chebyshev model is sharper. Note that the disparity between the two models becomes more pronounced with increasing D_{eq}, as shown in Fig. 4a–c. It can be noted from Fig. 4c that larger raindrops are more prone to aerodynamic effects and therefore have a flatter base compared to the reference “spheroid”.
Comparing the optical properties of two shapes with significantly different aspect ratios cannot reveal the effects resulting from introducing Chebyshev shapes. To address this, we defined the aspect ratio of the Chebyshev model in Fig. 4d, which represents the vertical maximum dimension b^{′} versus the horizontal maximum dimension a^{′}. Figure 5 illustrates a comparison between the aspect ratio of the Chebyshev model, as defined above, and the aspect ratio of the commonly used spheroid raindrop model (Thurai et al., 2007; Brandes et al., 2002). It is apparent that for common raindrops with D_{eq}<8 mm, the aspect ratios of the two models are approximately equal. Therefore, we can infer that the differences in optical properties between the spheroid and Chebyshev models arise from higherorder shape specifications rather than the firstorder aspect ratio parameter (also proved by comparisons between two shapes with identical aspect ratio; figure not shown).
A detailed examination of the optical properties of Chebyshev model rain droplets and their sensitivity against traditional spheroid model were conducted by Ekelund et al. (2020). To compare the radar variables between the spheroid and Chebyshev models, they used an extended precision version of the EBCM Tmatrix code (Mishchenko and Travis, 1994). However, it is important to note that the EBCM may encounter numerical instability issues due to the extremely high imaginary part of refractive indices for liquid water around the K band (10–40 GHz), where the Ku and Ka bands are located. To ensure integrity and accuracy, this study presents an optical property database (with a userfriendly radar operator interface) for the Chebyshev raindrop model at the Ku and Ka bands (13.6 and 35.5 GHz, respectively) using the IITM code (Bi et al., 2013).
To better understand the impact of singlescattering properties on radar variables, it is necessary to analyze the RSSP–level A database. Additionally, we choose to introduce intermediate quantities called the “SSP factors for radar variables”, which are illustrative and facilitate our understanding of how SSPs of a given particle size affect radar variables.
The SSP factor “[z_{h}]” for the horizontal reflectivity “z_{h}” is defined as
Here, the quantity enclosed in the square brackets [z_{h}] indicates the SSP factor of the radar variable z_{h}. We can perform a particle ensemble mean on it, which is often indicated by angle brackets, as follows (Zhang, 2016):
Here, the last equals sign is the definition of the horizontal reflectivity (see Eq. D1 in Appendix D).
Similarly, the SSP factor “[a_{h}]” for the horizontal attenuation coefficient “a_{h}” is defined as follows (the extinction cross section by essence):
If we perform particle ensemble mean on it, then we find that it is proportional to the definition of the horizontal attenuation coefficient (see Eq. D5 in Appendix D):
As shown by Eqs. (21) and (23), these factors ([z_{h}] and [a_{h}]) are the equivalent radar parameters of radar reflectivity z_{h} and specific attenuation a_{h}, respectively, for a single particle of size D_{eq}. For further SSP factors for a level A database diagnoses, please refer to Fig. 8.
It is not surprising to observe that the scattering properties exhibited symmetry with respect to the radar beam elevation angle, e=0°, for the spheroid model. This is because spheroids possess reflectance symmetry in the xy plane. However, this symmetry does not hold true for the Chebyshev model. We find that the SSP factor [z_{h}] has a more pronounced peak than spheroids near the nadir (specifically, e=90°, which is frequently encountered for the verticalprofiling cloud radar), as shown in Fig. 6a. Conversely, the peak near the zenith (i.e., $e=\mathrm{90}\mathit{\xb0}$ for the spaceborne radar) is weaker. This phenomenon is also found for the Chebyshev model at 94.1 GHz (Ekelund et al., 2020), which can be attributed to its flatter bottom surface and sharper top surface geometries in Ekelund et al. (2020), as depicted in Fig. 4. Additionally, we find that the [z_{h}] factor for groundbased scattering geometries (e=0–20°) is close to the values obtained for the spheroid model. However, deviations from the spheroid model were typically much smaller for lowerfrequency bands, such as the Ku band (13.6 GHz) (figure not shown here).
However, the forward scattering properties such as [a_{h}] of the Chebyshev model still maintain their symmetry with respect to the beam elevation angle e=0°, which can be easily understood based on the reciprocal theory (Van de Hulst, 1981). It is important to note that the attenuation effects of Chebyshev raindrops are slightly stronger compared to spheroid ones for nearnadir and zenithscattering geometries.
Figure 6 also demonstrates the stability of the deviations between the Chebyshev model and the spheroid model (CmS hereafter) in terms of the temperature dimension. However, the finding is different for [z_{h}] against the orientation preference dimension β (see Fig. 7a). It can be interpreted that a varying orientation presents a significantly different “view” of raindrops to the observer, while the varying temperature essentially keeps the same view of the raindrop with different dielectric constants. It was found that the positive and negative CmSs of [z_{h}] at nadir and zenith, respectively, hold true only if β<20°. In cases where particle groups have larger canting angles (β>20°), the CmS for [z_{h}] can produce a reversal of their signs (Fig. 7a). Since the column “orientation preference” in Table 1 shows that the standard deviation of β=7° for raindrops under normal turbulence conditions does not exceed this threshold for raindrops, we can infer that the CmS deviations found for β=0° only show a minor decrease if orientation averaging is performed.
As for [a_{h}], we learned that the conclusions of CmS at β=0° always hold true when β is sufficiently large (Fig. 7b).
We also examined the results across the entire PSD range of raindrops shown in Fig. 8. It was found that the aforementioned CmS results for [z_{h}] were significant when D_{eq}>3 mm (Fig. 8g), while those for [a_{h}] were significant only for larger raindrops (i.e., D_{eq}>5 mm; see Fig. 8j).
It is worth mentioning that raindrops had significantly higher [z_{h}] for larger absolute values of the beam elevation angle e (Fig. 8a). This can be easily understood since the zenith or nadir observation geometries tend to capture a larger cross section of oblatespheroidlike models.
As for other SSP factors of polarimetric radar variables, the CmS deviations are either not stable against size spectra (Fig. 8h and i) or too small in terms of its base magnitude (Fig. 8l), making them not significant for a particle group. However, the CmS effects of SSP factor for differential attenuation ([a_{h}] − [a_{v}]) in Fig. 8k is significant at low absolute elevation angles, indicating that the SSP factor of the vertical polarization attenuation [a_{v}] is significantly stronger for the Chebyshev model (considering [a_{h}] is close for two shapes at low absolute elevation angles in Fig. 8j).
The canting angle of raindrops is generally small, so the CmS deviations of SSP factors found in Figs. 6 and 8 only show a minor decrease in the level B ASSP database when compared to the curves of β=0° in the level A RSSP database. Hence, we omit any further analysis of the orientation averaging of raindrops and SSP factors of the level B database here.
3.4 PSD options and analysis of level B to level C LUT conversion
In this subsection, we will focus on describing the PSD options for raindrops and the methods used to analyze their impacts on the conversion of SSP factors of radar variables to intrinsic radar variables (i.e., the ASSP to BSP conversion).
In ZJUAERO, a total of six options for PSD schemes for raindrops were available, all of which are listed in Table 3. All schemes were designed for rain modeled by a singlemoment microphysics scheme with exponential distribution assumptions. Each PSD schemes can be visualized as a N_{0}–Λ curve, as shown in Fig. 9.
The constant N_{0} parameterization originally proposed by Marshall and Palmer (1948) and used in microphysics packages, such as Hong and Lim (2006), provides a rough approximation for various liquid precipitation scenarios. However, modern PSD schemes for rain, such as the one by Thompson et al. (2008), aim to capture the observed transition from high concentrations of drizzlesized drops (D_{eq}<0.5 mm) in stratocumulus clouds to PSDs dominated by large millimetersized raindrops in heavy precipitation.
We classified the six PSD schemes into two groups by how the intercept parameter N_{0} is determined. Group A is characterized by a power law N_{0}–Λ relationship, which is represented as a straight line in the duallogscale diagram in Fig. 9. On the other hand, the intercept parameter, N_{0}, for group B schemes follows a formalization described by Eq. (24), transitioning from N_{2} to N_{1} as the water mass concentration, Q_{R}, increases. This can be visualized as the tanhlike curves shown in Fig. 9:
The values of the reference rainwater mass concentration Q_{R0} and the dynamic range [N_{2}, N_{1}] of the intercept parameter are displayed in the group B in Table 3.
^{*} The ThompsonTuned PSD scheme is for numerical tests only; therefore, it is not a user option in Table 1.
Figure 9 demonstrates the large variability in the intercept parameter N_{0} for a specified rainwater content (RWC) among different PSD schemes. That variability can be shown by measuring the vertical coordinate difference in the intersection points between a given isoRWC line and different N_{0}–Λ lines of PSD schemes. For PSD schemes expressed by an exponential distribution, a larger intercept parameter N_{0} means more smaller drops when the RWC is fixed. Therefore, the Thompson scheme (Thompson et al., 2008) prioritizes smaller drops, while the Wang scheme (Wang et al., 2016) places more emphasis on larger drops.
According to global aircraft in situ measurements carried out by Abel and Boutle (2012), the six singlemoment (SM) microphysics assumptions in Table 3 cover the natural variability in the raindrop PSDs ranging from continental precipitation to maritime precipitation. Typical continental PSD with intensive large drops, such as in Wang et al. (2016), and typical maritime PSD, such as in Thompson et al. (2008), with a large population of drizzlesized raindrops are both present. That is why we have chosen the six SM microphysics schemes for the PSD sensitivity test.
It should be noted that although SM microphysics has been used for years and will continue to be used, it has obvious shortcomings with respect to reproducing the polarimetric features observed in convective and stratiform events. Since doublemoment (DM) microphysics in CMAMESO is still experimental, the consistent DM scheme in ZJUAERO is still in the development stage. However, the natural variability in the raindrop PSDs for DM microphysics schemes is still within what is displayed in Fig. 9; therefore, it is safe to use those SM assumptions for a sensitivity test.
To explore the effects of different PSD schemes on the radar variables, a new intermediate quantity called “backscatagrand” is analyzed in Fig. 10, which is a spectra of particle size defined as the product of the horizontal backscattering cross section σ_{bsca,h}(D) and size distribution function N(D). It helps us identify particles for which the size range dominates the energy of the backscattering.
The idea of analyzing backscatagrand was inspired by the definition of “extagrand” used in the analysis of the RTTOVSCATT bulkscattering tables (Geer et al., 2021), which helps with diagnosing the singleparticle contribution to the extinction coefficients in the radiative transfer simulations of passive microwave sounders. The concept of extagrand is based on the insight that the total extinction of an ensemble of particles is primarily influenced by a fraction of particles within a narrow diameter range. Therefore, by analyzing the SSP for particles within that size range, we can infer the BSP of the entire ensemble of particles. The extagrand σ_{ext}(D)N(D) (in units of mm m^{−3}) is factorized as the product of the mass distributions m(D)N(D) (in units of kg mm^{−1} m^{−3}) and the extinction cross section per unit volume ${\mathit{\sigma}}_{\mathrm{ext}}\left(D\right)/m\left(D\right)$ (in units of mm^{2} kg^{−1}). These quantities are determined by PSDs and SSPs, respectively. The extagrand quantities are appreciable only for size ranges in which both the mass distribution m(D)N(D) and mass extinction efficiency ${\mathit{\sigma}}_{\mathrm{ext}}\left(D\right)/m\left(D\right)$ are large enough to produce significant products of extagrand.
However, for weather radar applications, the principal quantity is backscattering rather than the extinction, as is the case for passive instruments. Hence, we try to analyze quantities of the horizontal polarization backscattering cross section per volume ${\mathit{\sigma}}_{\mathrm{bsca},\mathrm{h}}\left(D\right)/m\left(D\right)$ (parallel definition of mass extinction efficiency, which we refer to as mass backscattering efficiency hereafter) and the backscatagrand σ_{bsca,h}(D)N(D) (parallel definition of extagrand) in which σ_{bsca,h}(D) indicates the horizontal polarization backscattering cross section.
In Fig. 10, we can find that the mass distributions and backscattering contribution spectrum (i.e., the socalled backscatagrand) exhibit large uncertainties due to the natural variability in the PSD schemes. The mass backscattering efficiency has multiple oscillations caused by the resonance effect, particularly in highfrequency bands such as the Ka band. Specifically, two important peaks of mass backscattering efficiency at D_{eq}=2.5 mm and D_{eq}=6.5 mm and one important dip at D_{eq}=5.0 mm exist (solid blue line in Fig. 10d). For extremely heavy precipitation (${Q}_{R}\sim {\mathrm{10}}^{\mathrm{2}}$ kg m^{−3}), the peak of the mass distributions might coincide with the dip of the mass backscattering efficiency at D_{eq}=5.0 mm. This phenomenon is unique to PSD schemes that contain more large drops (such as Wang2016 and AB2012), which leads to a loss of the bulk reflectivity if the total mass of the raindrop remains constant.
Under all but extremely heavyprecipitation conditions (Q_{R}<10^{−2} kg m^{−3}; see Fig. 10c and g), the Thompson2008 PSD scheme stands out as an outlier. Its extremely large N_{0} compared to other schemes (see Fig. 9) leads to a significantly smaller peak in mass distributions as small as D_{eq}<1.0 mm, while the peak of other PSD schemes is approaching the first peak of the mass backscattering efficiency at D_{eq}=2.5 mm. Accordingly, the total reflectivity computed with the Thompson2008 scheme must be much smaller at ${Q}_{R}={\mathrm{10}}^{\mathrm{3}}$ kg m^{−3} than those computed with other PSDs.
The relative importance of particles in the entire spectrum can be diagnosed with the curves of backscatagrand, which is the ultimate goal to introduce this concept. For example, we can learn from Fig. 10h that the contribution of backscattering is dominated by particles with the diameter D_{eq} at around 2.5 mm for the MP1948 scheme, while the contribution from particles of 2.5 mm size and 6.5 mm size are almost equally important for the Wang2016 and AB2012 schemes. To go further, we can infer that the CmS deviations mentioned in Sect. 3.3 can only affect the bulkscattering properties computed with the Wang2016 and AB2012 schemes, since the CmS deviations of backscattering are only significant for particles with D_{eq}>3 mm (Fig. 8g).
3.5 Bulkscattering properties: level C
Figure 11 shows the intrinsic radar variables for the Chebyshev model which were computed using the BSP database and the spheroid model at the Ka band. In Sect. 3.4, we hypothesized based on the backscatagrand plot, and we can now confirm these hypotheses individually.

The application of PSD schemes, such as Wang2016, can result in a significant reduction in the reflectivity Z_{H} (by approximately 5 dBZ) for extremely heavyprecipitation scenarios compared to the default PSD scheme MP1948 (refer to Fig. 11a; ${Q}_{R}\sim {\mathrm{10}}^{\mathrm{2}}$ kg m^{−3}).

The reflectivity Z_{H} computed by the Thompson2008 scheme for moderately heavyprecipitation scenarios is considerably lower (by over 10 dBZ) compared to other PSD schemes in group A in Table 1 (refer to Fig. 11a' ${Q}_{R}\sim {\mathrm{10}}^{\mathrm{3}}$ kg m^{−3}).

The CmS deviations are only significant (reducing Z_{H} by about 2 dBZ) for extremely heavyprecipitation scenarios and PSD schemes that emphasize larger drops, such as Wang2016 (refer to Fig. 11a, ${Q}_{R}\sim {\mathrm{10}}^{\mathrm{2}}$ kg m^{−3}).

The CmS deviations of attenuations, α_{H}, are never significant for regular mass concentrations of rain Q_{R} (see Fig. 11b). This finding can be explained by the fact that the CmS deviations in RSSP (Fig. 8j) are significant only if D_{eq}>7 mm, which represents a group of drops sharing a small fraction of the mass distribution, even for extremely heavy precipitation.
From Fig. 11, it can be concluded that the horizontal polarization intrinsic radar variables Z_{H} and α_{H} at the Ka band (actually also for the Ku band; not shown here) are more sensitive to the uncertainty in the PSD schemes than the CmS optical property deviations for spaceborne radar observation geometries. The CmS deviations are only notable for extremely heavyprecipitation scenarios in which large drops are prominent. Those conclusions focusing on the Ka band (35.6 GHz) generally agree well with what was found by Ekelund et al. (2020).
The sensitivities of intrinsic polarimetric radar variables with respect to CmS deviations for groundbased radar bands (S, C, or X band) and viewing geometries were found to be negligible. Therefore, they are not shown here.
To assess the simulation capabilities of ZJUAERO, we conducted case studies using realworld data and investigated the sensitivities of the PSD schemes and the new Chebyshev raindrop model. For this study, we chose the GPM core satellite (referred to as GPM hereafter) and specifically analyzed the overpass of Typhoon Haishen, the first supertyphoon of the 2020 northwest Pacific typhoon season, on 5 September 2020 at 09:21 UTC. We used the simulation and observation data from the dualfrequency precipitation radar (DPR) aboard the GPM for our analysis.
To conduct the ZJUAERO simulations, we utilized model grid data generated by the operational run of CMAMESO, which was initialized at 00:00 UTC on 5 September 2020. The CMAMESO grid data had a horizontal resolution of 3 km and consisted of 50 vertical layers. Note that the WSM6 microphysics package was selected for the CMAMESO operational run. However, any forced PSD schemes could be applied in the simulation of the radar operator, as mentioned in Sect. 2.3.
To demonstrate the reliability of the forward radar operator for groundbased polarimetric radars, we have also performed a case study of a mesoscale convective system (MCS) which can be found in the user manual (see the “Code and data availability” section). The results are reasonable but relatively trivial compared to previous radar operators, so we will not display them in this section.
4.1 Simulation of a tropical cyclone case
According to Iguchi (2020), the GPM/DPR's Kuband radar basically follows the instrument characteristics of the Tropical Rainfall Measuring Mission (TRMM) Precipitation Radar (PR), while the new Kaband radar is sensitive to light rain and snow. By combining data from two channels, more accurate estimates of DSD parameters can be obtained (Rose and Chandrasekar, 2006; Liao et al., 2014). The Kaband radar has two modes (Iguchi et al., 2010), namely (1) a highsensitivity mode for light rain and snow (highsensitivity beam) and (2) a matchedbeam mode in which the sampling volumes of the Ka and Kuband radar are collocated (matched beam). However, since 21 May 2018, the GPM/DPR has switched its scan pattern so that now a full swath can be considered the matched beam in the evaluation of a forward radar operator. Therefore, in this case study, we used data from the Ka band in the matchedbeam mode to estimate DSD and drop morphology parameters.
The dualfrequency ratio (DFR) is defined as the difference between the logspacemeasured reflectivity at two channels (Ku and Ka bands). Previous studies have shown that the DFR can be used to distinguish between stratiform and convective rain (Le and Chandrasekar, 2012).
Figure 12 displays the observation and simulation of the Kuband radar reflectivity at different levels, with the last row showing the mismatch between them (i.e., the OmB of reflectivity). Figure 13 presents the cross section of radar reflectivity between points A and B in Fig. 12 separately for the Ku and Ka bands. Additionally, the last column shows the DFR as defined above. The radar operator applied the ThompsonTuned PSD option (for developers only) and the Chebyshev morphology options A2 (see Table 3) for the category rain in the simulations. For snow, we use default option A2. For graupel, we use default option A2 and B1. The reflectivities in Fig. 12a–d are masked by a flag called “flagPrecip” (available at the ground) offered in the L2A swath data of GPM/DPR, while the simulation of reflectivity in Fig. 12e–h is masked by the sensitivity threshold of 12 dBZ to keep it the same with sensitivity of GPM/DPR observation. We applied the attenuation in simulations while using the “zFactorMeasured” product of GPM/DPR (no attenuation correction applied). As for the calculation of OmB reflectivity, the radar gates for which the reflectivity is undetectable (below the instrument sensitivity) either in the simulation or observation are filled with a “background” reflectivity of 0 dBZ to generate the OmB reflectivity.
Based on the comparison between observations and simulations in Figs. 12 and 13, it can be found that the regional model of CMAMESO is able to capture some of storm structures, such as the cyclone eye and eyewall. However, the structures of outer spiral rain bands in the simulation (Figs. 12e and 13d) appear more contiguous and vague compared to the isolated towering bands depicted by the GPM/DPR measurements (Figs. 12a and 13a). Although the NWP model could not accurately predict the cloud and precipitation timing and position, the probability distribution of the simulated radar reflectivity should be unbiased when compared with observations (detailed examinations on that will be conducted in Sect. 4.2); otherwise, a systematic bias in the NWP model or the radar operator might be identified. Since (a) the precipitation forecast of CMAMESO model can be frequently calibrated against a rain gauge network (Bárdossy and Das, 2008; Cattoën et al., 2020), and (b) a forecast lead time of 9 h is beyond the 6 h spinup time of the rainwater content in CMANWP (Ma et al., 2021), we assume that the NWP model CMAMESO has no significant bias in this case.
Notably, the freezing level in this case was found within the altitudes ranging from 4 to 5 km (see the 0 °C NWP model background isothermal line in Fig. 13), which is believed to contain abundant melting particles. Actually, the GPM/DPR observations revealed a weak brightband (BB) signature below the freezing level, which can be recognized at the Ku band (with an average reflectivity enhancement of approximately 3 dB; see Fig. 13a) but is not so clear at Ka band (see Fig. 13b). As reported by longterm radar observation statistics of melting layers documented in Fabry and Zawadzki (1995), the magnitude of the BB signature in the melting layer is significantly weaker in deepconvection regions. Therefore, the weak BB signatures in the melting layer of this tropical cyclone case can probably be attributed to the large riming rates (Zawadzki et al., 2005) in the convective precipitation of tropical cyclone eyewall and rain bands. However, the simulations we conducted did not show a BB signature, which was attributed to not considering melting or mixedphased particles and the shortcoming of the microphysics scheme. Considering the lack of melting schemes in ZJUAERO, it is expected that the melting layer would exhibit a positive OmB signature due to lack of melting hydrometeors in ZJUAERO. However, in Fig. 13g and h, we encountered difficulties in identifying a continuous positive mean bias of the OmB for both the Ku and Ka bands in the melting layer. This challenge arose due to the large mislocation errors in the precipitation predicted by the NWP model used in this study. To address this issue, future analysis could employ horizontal averaging or examine the probability distribution function of reflectivities in the melting layer. Implementation of melting particle models and their relevant validations will be a subject for upcoming publications.
Moreover, extremely large DFR values (>30 dB) were found in simulations (Fig. 13f) but not in observations (Fig. 13c) in the 0 to 2 km layer, which could be attributed to an overestimation of the attenuation in the Kaband simulation for heavy precipitation (Fig. 13e). This hypothesis is supported by two facts: (1) many weak reflectivity (∼ 10 dBZ) regions are found right beneath the strong reflectivity gates aloft at around 4 km for the Ka band (Fig. 13e), and (2) the weak reflectivity regions of the Ka band collocate with the strong reflectivity region (∼ 40 dBZ) of the Ku band (Fig. 13d).
The OmB plots, as shown in Figs. 12 and 13, are useful tools for verifying and calibrating new observation operators and identifying their deficiencies. Overall, the results are generally reasonable and demonstrate the capability of ZJUAERO to simulate the reflectivity of dualfrequency spaceborne radar such as GPM/DPR. More analysis on the bias of the probability distribution of reflectivities/DFR will be presented in the next subsection.
4.2 Sensitivity assessments on hydrometeor PSDs and morphologies
We also performed sensitivity assessments for the PSD and nonspherical morphology options of the rain hydrometeor category. Figure 14 shows the observation and simulation distributions of radar reflectivity at the Ku and Ka bands and the dualfrequency ratio (DFR) at two altitude layers (0–2 and 2–4 km). Since we were tuning the liquid hydrometeor, the melting and solid layers in the last two columns (4–6 and 6–8 km) did not concern us. The radar operator applied the Chebyshev morphology in all simulations except the group with the label “ThompsonTuned(Spheroid)”, which was simulated with the spheroid raindrop option. Other options of ZJUAERO are identical with those in Fig. 12.
Based on statistical analysis (see Fig. 14g and h), we found that the mass concentration of rain in the CMAMESO at the storm rain bands was primarily dominated by moderately heavy rain (${Q}_{R}\sim {\mathrm{10}}^{\mathrm{3}}$ kg m^{−3}) rather than extremely heavy rain (${Q}_{R}\sim {\mathrm{10}}^{\mathrm{2}}$ kg m^{−3}). According to Sect. 3.5 and Fig. 11, applying the Thompson2008 PSD leads to a significant reduction in the simulated reflectivity in conditions of moderately heavy rain (${Q}_{R}\sim {\mathrm{10}}^{\mathrm{3}}$ kg m^{−3}), which is consistent with our findings in Fig. 14a, b, c, and d for the Ku and Ka bands of liquid layers.
Due to the strong wind shear in the rain bands of the tropical cyclone, large drops can be broken apart (Radhakrishna, 2022), causing the rain PSD in such conditions to behave irregularly and deviate from the prevailing parameterizations. Figure 14a, b, c, and d demonstrate that almost all the PSD schemes significantly overestimated the reflectivity for both Ku and Ka bands, except for the Thompson2008 scheme, which underestimated the reflectivity.
To address the discrepancy, we designed a new PSD scheme, referred to as ThompsonTuned, by tuning the parameters in Eq. (24). The optimization procedure can be formulated as follows.

We use a scoring method to quantify the match between observation and simulation histograms, as proposed by Geer and Baordo (2014):
$$\begin{array}{}\text{(25)}& s=\sum _{j=\mathrm{1}}^{\mathrm{nbands}}\sum _{i=\mathrm{1}}^{\mathrm{nbins}}\left\mathrm{log}\left({\displaystyle \frac{{N}_{\mathrm{sim}}^{j}\left(i\right)}{{N}_{\mathrm{obs}}^{j}\left(i\right)}}\right)\right.\end{array}$$Here, s represents the final score, with N_{sim} and N_{obs} indicating the counts in the ith bins of the simulation and observation histograms, respectively. The superscript “j” denotes jth band. To prevent infinite values in the summation, we assume a count of 0.1 in bins with empty values in either the simulation or observation histograms.

Next, we establish a grid of free parameters for optimization which includes three parameters (N_{1}, N_{2}, and Q_{R0}). Due to the broad range of these free parameters, the grid is set up in a quasilogarithmic scale.

Subsequently, simulations are conducted, and the score from Eq. (25) is evaluated for each grid point. By identifying an optimal region in the parameter space with the best score, we refine the grid in that region to pinpoint a more precise parameter subregion. This iterative process is carried out to finetune the parameters.

Finally, we identified an optimal point at which the black lines in Fig. 14a, b, c, and d closely matched the observed distributions. The parameters are listed in Table 3 and plotted in the N_{0}–Λ diagram in Fig. 9. The tuned PSD scheme is reasonable as it falls between the MP1948 and Thompson2008 schemes in the N_{0}–Λ diagram, with an emphasis on smaller particles compared to MP1948. This may be attributed to the unusual DSD in tropical cyclones.
While no in situ DSD observations are available to support the tuned PSD schemes used in this study, the implications of the tuning experiments are interesting, considering that the matchedbeam observation of the Ku and Ka bands were designed for the DSD estimation.
As suggested by Sect. 3, the CmS effects were negligible in this case, as it is difficult to discriminate the solid and dashed black lines in Fig. 14. This is likely due to the moderately heavy rainfall in tropical cyclones with strong wind shears which prevent the large drops from highlighting the CmS deviations. However, the CmS effects could be significant in cases of extremely heavy rain in the supercell storms, according to our observation system simulation experiment (OSSE). An OSSE of GPM/DPR overpassing a mesoscale convective system (recorded with extreme heavyprecipitation events) reported CmS decrease effects of more than 1 dB at the Ku band and more than 2 dB at the Ka band (figures not shown). This sensitivity test is consistent with what we found in Fig. 11a and demonstrates the value of introducing the Chebyshevshaped raindrop model in certain scenarios (e.g., vertical pointing cloud radar, airborne radar, and spaceborne radar working at highfrequency bands). As for polarimetric radar variables such as Z_{DR} and K_{DP} for groundbased radar at sideviewing geometry, the CmS effects are generally negligible.
In summary, Sect. 2 of our study introduced the basic formulations and concepts of the design in the ZJUAERO. These concepts included the general procedure of the software, radarvariable calculations, and the available hydrometeor settings (shape parameters, dielectric constant models, canting angles, and particle size distributions). Formulations of polarimetric radar variables are derived starting from singleparticle backscattering Mueller matrix Z and extinction matrix K.
In Sect. 3, we highlighted the unique features of ZJUAERO, specifically its multilayered design for the optical database of nonspherical particles. We demonstrated this by displaying the scattering properties using the example of the Chebyshevshaped raindrop particle model and comparing it to the properties of traditional spheroid raindrops. We conducted LUT demonstrations for two database layers, namely level A (raw singlescattering properties database) and level C (bulkscattering property database). We also introduced a new intermediate quantity named backscatagrand to diagnose the PSD integrations of optical properties. We concluded that the Chebyshevshaped raindrop model shows noticeable differences in the bulkscattering properties (compared to spheroid model) only at zenith and nadirviewing geometries and for millimeterwavelength radar bands. This difference is more prominent for continental DSD models (e.g., Wang2016) with larger drops. These deviations can reach up to 2–3 dB at the Ka band for spaceborne radars in heavy continental precipitation regions where large drops dominate. Given the lower uncertainties in simulating the reflectivity of the liquid phase compared to the solid and melting phases, such a difference deserves attention in specific applications such as comparing data from groundbased and spaceborne radar observations (Warren et al., 2018).
Furthermore, in Sect. 4, we validated the simulation reliability and capability of ZJUAERO by analyzing a case study of a tropical cyclone using input from the CMAMESO for simulating spaceborne radar observations. We found that ZJUAERO provides reasonable simulation results, except for the brightband signature at the melting layer, which can be attributed to the current version of ZJUAERO and does not consider melting or mixedphased particles. We also performed sensitivity assessments of PSDs and morphology options for rain in ZJUAERO and found that the ThompsonTuned singlemoment PSD scheme provides the best fit of the reflectivity histogram in the simulation to the GPM/DPR observation. However, using either the Chebyshevshaped raindrop particle model or the spheroid model makes little difference to the simulation results since the tropical cyclone precipitation has a maritime DSD where small drops dominate.
Currently, ZJUAERO is an efficient forward radar operator that has the advantage of handling the complexities of nonspheroid particle models. Therefore, it is a powerful research tool for studying the sensitivities of polarimetric radar observations with respect to the nonsphericity of hydrometeor particles. It also applies parallel acceleration techniques to boost its performance, allowing operational applications of data assimilation in NWP models employing singlemoment (SM) microphysics (such as CMAGFS/MESO). A ZJUAERO forward radar operator can be applied in data assimilation (DA) studies using indirect DA methods such as the Bayesian approach (Caumont et al., 2010) and direct DA methods such as ensemble Kalman filter (EnKF) that both require no tangent linear (TL) or adjoint (AD) versions of the forward radar operator.
ZJUAERO also has some limitations. For example, it currently cannot be applied in DA research based on the variational method. Simplification and linearization works are involved to obtain a TL/AD version of the forward operator. Moreover, PSD solvers for DM microphysics schemes have already been implemented in ZJUAERO for the experimental CMAMESO DM microphysics, but there are many validation and evaluation works to be done. Also, unlike the EMVORADO, which use a distributedmemory parallel design and interface with the COSMONWP model online, ZJUAERO applies a sharedmemory parallel design, and the NWP model input should be stored in external files.
Overall, the results were satisfactory, and the ZJUAERO operator is ready for experimental and operational usage. However, several aspects of ZJUAERO still need improvement, and we have listed the ongoing development tasks as follows:

Improve frozen hydrometeor modeling by introducing irregularshaped aggregates and riming snow.

Improve the modeling of melting particles (including melting snow, melting graupel, and melting hail) using layered inhomogeneous modeling rather than an effective medium approximation for the mixture matrix.

Compute the optical properties of melting particles using the IITM code and extend the SSP/BSP database with a new dimension of the water fraction.

Develop more concrete models for singlecrystal, aggregated, and rimed snow to replace the “soft spheroid” model and create a corresponding SSP/BSP database.

Model hail as nonspheroid and inhomogeneous particles.

Include cloud ice, which plays an significant role in the highfrequency bands of spaceborne radar.

Do more tests for doublemoment microphysics schemes in CMAMESO.

Conduct more case studies, particularly with measurements from the spaceborne radar FY3 RM/PMR. Notably, the L1 product of spaceborne radar FY3 RM/PMR has been accessible online (released by the National Satellite Meteorological Center (NSMC) on 22 November 2023). We have already implemented an external I/O (input/output) module to interface with the data format of FY3 RM/PMR in ZJUAERO, but the time coverage of observation data is too limited for us to find a good demonstration case. Therefore, more case studies and finetuning should be conducted with future measurements from FY3RMPMR.
In conclusion, ZJUAERO is an observation operator that facilitates the exploitation of measurement data from both spaceborne and groundbased radars. Its versatility and effectiveness make it a valuable tool for data assimilation in CMAGFS/MESO. Moreover, ZJUAERO has the potential to be applied in various other studies within a wide range of contexts.
For spaceborne radars aboard rain measurement satellite platforms, such as the Tropical Rainfall Measuring Mission/Global Precipitation Measurement/FengYun3 Rain Measurement (TRMM/GPM/FY3 RM), the trajectory of the radar beam can be treated without considering the beambending effects while still maintaining precision. The WGS84 coordinates of satellites, denoted as A (scLat, scLon, scAlt), in addition to the center of the footprint, A_{1} (lat, lon), can be obtained through the satellite radar L1/L2 products (hereafter referred to as swath files). These coordinates can then be used to calculate the local elevation angle of a given radar gate, C, using the knowledge of trigonometry in Fig. A1b. The length of segments H and R_{E} can be computed by converting the (latitude, longitude, and altitude) WGS84 coordinates to the Earthcentered, Earthfixed (ECEF) coordinates. The range of radar gate AC is provided by the spaceborne radar observation system in the L1 product of GPM/DPR known as “scRangeEllipsoid” (Iguchi et al., 2010). When neglecting the beambending phenomenon in the spaceborne radar detection, the local elevation angle, e^{′}, can be expressed as ${e}^{\prime}=e\mathit{\alpha}$ in which e is the elevation angle on the satellite. The angle α could be determined using trigonometry, given that AC represents the range of the radar gate C.
The grids of LUT dimensions are presented in Table B1. Only the Z and K matrix elements that appear in the formulas of polarimetric radar variables in Eqs. (D1)–(D8) are stored in the database. For temperatures, we set the minimum temperature for supercooled liquid particles as 233 K, considering that homogeneous freezing starts at lower temperatures. For solid hydrometeors, the lowest temperature in the database is 203 K, and lower environment temperatures encountered are taken as 203 K. We provided LUTs for six bands that are widely used for groundbased and spaceborne weather radars. The ranges of diameters and aspect ratios for solid hydrometeors are suggested by observations of field research (Garrett et al., 2015), while the range of diameters and aspect ratios of the liquid hydrometeor have already been discussed in Sect. 3.3. The number of diameter bins is sufficient to produce a reasonable size spectrum of hydrometeors for various hydrometeor mass concentrations, as shown in Fig. 10. The range of mass concentration grids for BSP LUTs is the same as that of Geer et al. (2021). Other external LUTs, such as Eriksson et al. (2018), need to undergo format conversions and interpolations (e.g., for different diameter grids) before it can be applied in ZJUAERO.
To alleviate the truncation and quadrature error in the particle size integration, a renormalization technique is applied after the particle number in each size bin is calculated. We first calculate the renormalization factor κ as follows:
Here, m(D_{i}) is the mass of hydrometeor particle with a diameter of D_{i}, and N(D_{i}) is the particle number in that size bin, while ΔD is the diameter step of the numerical integration. Then (κ<1) gives the ratio between the mass concentration presented by our PSD and the mass concentration of hydrometeors given by NWP model Q_{m}. Then we can calculate the corrected particle number distribution N^{corr}(D_{i}):
Although the optical properties of hydrometeors are consistently computed and stored in the FSA convention, this convention is not used for monostatic radar applications. Figure C1 displays the differences between the backscattering alignment (BSA) convention and the FSA convention (Chandrasekar, 2001). The incident light propagates along −X_{L} direction and comes into contact with the particle at the origin O. $\widehat{h}$, $\widehat{v}$, $\widehat{k}$ are used to represent unit vectors of horizontal polarization, vertical polarization, and the propagation direction of electromagnetic waves, respectively. If we use the definition of $\widehat{h}$ and $\widehat{v}$ unit vectors under the FSA convention, then the horizontal unit vector of backscattering light is inconsistent with the incident light (${\widehat{h}}_{\mathrm{FSA}}^{\mathrm{inc}}$, ${\widehat{v}}_{\mathrm{FSA}}^{\mathrm{inc}}$, ${\widehat{k}}_{\mathrm{FSA}}^{\mathrm{inc}}$) = (${\widehat{h}}_{\mathrm{FSA}}^{\mathrm{bsca}}$,${\widehat{v}}_{\mathrm{FSA}}^{\mathrm{bsca}}$, ${\widehat{k}}_{\mathrm{FSA}}^{\mathrm{bsca}}$), where “bsca” indicates backscattering light. To resolve this convention conflict between the polarimetric radar observation system and scattering computation, BSA forces the following relationship by reversing the direction of the backscattering horizontal unit vector definition: (${\widehat{h}}_{\mathrm{FSA}}^{\mathrm{bsca}}$, ${\widehat{v}}_{\mathrm{FSA}}^{\mathrm{bsca}}$, ${\widehat{k}}_{\mathrm{FSA}}^{\mathrm{bsca}}$) = (${\widehat{h}}_{\mathrm{BSA}}^{\mathrm{bsca}}$, ${\widehat{v}}_{\mathrm{BSA}}^{\mathrm{bsca}}$, ${\widehat{k}}_{\mathrm{BSA}}^{\mathrm{bsca}}$). Hence, the amplitudescattering matrix in the BSA convention is related to that of FSA as follows:
The following amplitudescattering elements used for polarimetric radar variable computation are represented in the BSA convention, unless otherwise stated (i.e., a conversion is needed in the core submodule before obtaining radar variables).
Next we describe how to perform intrinsic radar variable calculations using bulk matrices 〈Z〉 and 〈K〉:

Horizontal/vertical reflectivity ${z}_{\mathrm{h}/\mathrm{v}}$ (in units of mm^{6} m^{−3}):
$$\begin{array}{}\text{(D1)}& {\displaystyle}\begin{array}{l}{z}_{\mathrm{h}}=\frac{{\mathit{\lambda}}^{\mathrm{4}}}{{\mathit{\pi}}^{\mathrm{5}}\mathrm{}{K}_{\mathrm{w}}{\mathrm{}}^{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{\mathit{\sigma}}_{\mathrm{bsca},\mathrm{h}}N\left(D\right)\mathrm{d}D,\\ =\frac{\mathrm{4}\mathit{\pi}{\mathit{\lambda}}^{\mathrm{4}}}{{\mathit{\pi}}^{\mathrm{5}}\mathrm{}{K}_{\mathrm{w}}{\mathrm{}}^{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{\left{S}_{\mathrm{hh}}\right}^{\mathrm{2}}N\left(D\right)\mathrm{d}D,\\ =\frac{\mathrm{2}{\mathit{\lambda}}^{\mathrm{4}}}{{\mathit{\pi}}^{\mathrm{4}}\mathrm{}{K}_{\mathrm{w}}{\mathrm{}}^{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}({Z}_{\mathrm{11}}{Z}_{\mathrm{12}}{Z}_{\mathrm{12}}+{Z}_{\mathrm{22}})N\left(D\right)\mathrm{d}D,\\ =\frac{\mathrm{2}{\mathit{\lambda}}^{\mathrm{4}}}{{\mathit{\pi}}^{\mathrm{4}}\mathrm{}{K}_{\mathrm{w}}{\mathrm{}}^{\mathrm{2}}}\left(\u2329{Z}_{\mathrm{11}}\u232a\u2329{Z}_{\mathrm{12}}\u232a\u2329{Z}_{\mathrm{21}}\u232a+\u2329{Z}_{\mathrm{22}}\u232a\right),\end{array}\text{(D2)}& {\displaystyle}{z}_{\mathrm{v}}={\displaystyle \frac{\mathrm{2}{\mathit{\lambda}}^{\mathrm{4}}}{{\mathit{\pi}}^{\mathrm{4}}\mathrm{}{K}_{\mathrm{w}}{\mathrm{}}^{\mathrm{2}}}}\left(\u2329{Z}_{\mathrm{11}}\u232a+\u2329{Z}_{\mathrm{12}}\u232a+\u2329{Z}_{\mathrm{21}}\u232a+\u2329{Z}_{\mathrm{22}}\u232a\right).\end{array}$$In Eq. (D1), σ_{bsca,h} indicates the horizontal backscattering cross section of a particle (i.e., ${\mathit{\sigma}}_{\mathrm{bsca},\mathrm{h}}=\mathrm{4}\mathit{\pi}{\left{S}_{\mathrm{hh}}\right}^{\mathrm{2}}$). λ is the wavelength of radar, ${K}_{\mathrm{w}}=\left\left({\mathit{\epsilon}}_{\mathrm{w}}\mathrm{1}\right)/\left({\mathit{\epsilon}}_{\mathrm{w}}+\mathrm{2}\right)\right$ is the dielectric factor (ε_{w} is the dielectric constant of water at the wavelength of radar for a fixed temperature of 283.15 K). Similarly, vertical reflectivity can be derived in Eq. (D2).

Differential reflectivity z_{dr} (dimensionless):
$$\begin{array}{}\text{(D3)}& {z}_{\mathrm{dr}}={\displaystyle \frac{{z}_{\mathrm{h}}}{{z}_{\mathrm{v}}}}={\displaystyle \frac{\u2329{Z}_{\mathrm{11}}\u232a\u2329{Z}_{\mathrm{12}}\u232a\u2329{Z}_{\mathrm{21}}\u232a+\u2329{Z}_{\mathrm{22}}\u232a}{\u2329{Z}_{\mathrm{11}}\u232a+\u2329{Z}_{\mathrm{12}}\u232a+\u2329{Z}_{\mathrm{21}}\u232a+\u2329{Z}_{\mathrm{22}}\u232a}}.\end{array}$$ 
The specific differentialphase shift (° km^{−1}) upon propagation, K_{DP}, is defined by
$$\begin{array}{}\text{(D4)}& \begin{array}{l}{K}_{\mathrm{DP}}={\mathrm{10}}^{\mathrm{3}}\cdot \frac{\mathrm{180}}{\mathit{\pi}}\cdot \mathit{\lambda}\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}\mathfrak{R}({S}_{\mathrm{hh}}^{\mathrm{fwd}}{S}_{\mathrm{vv}}^{\mathrm{fwd}})N\left(D\right)\mathrm{d}D\\ ={\mathrm{10}}^{\mathrm{3}}\cdot \frac{\mathrm{180}}{\mathit{\pi}}\cdot \underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}\frac{\mathrm{2}\mathit{\pi}}{{k}_{\mathrm{0}}}\mathfrak{R}({S}_{\mathrm{hh}}^{\mathrm{fwd}}{S}_{\mathrm{vv}}^{\mathrm{fwd}})N\left(D\right)\mathrm{d}D\\ ={\mathrm{10}}^{\mathrm{3}}\cdot \frac{\mathrm{180}}{\mathit{\pi}}\cdot \u2329{K}_{\mathrm{34}}\u232a.\end{array}\end{array}$$Here, 10^{−3} in Eq. (D4) represents the coefficient in the unit conversion (from mm^{2} m^{−3} to km^{−1}), and the coefficient $\mathrm{180}/\mathit{\pi}$ is used to convert the unit of the result (from radii km^{−1} to ° km^{−1}). Here, the superscripts of “fwd” over the matrix elements indicate that these are elements of the forward scattering matrix.

The oneway linearscale attenuation coefficient^{5} of horizontal/vertical polarization ${a}_{\mathrm{h}/\mathrm{v}}$ (in units of km^{−1}):
$$\begin{array}{}\text{(D5)}& {\displaystyle}\begin{array}{l}{a}_{\mathrm{h}}={\mathrm{10}}^{\mathrm{3}}\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{\mathit{\sigma}}_{\mathrm{ext},\mathrm{h}}N\left(D\right)\mathrm{d}D\\ ={\mathrm{10}}^{\mathrm{3}}\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}\frac{\mathrm{4}\mathit{\pi}}{{k}_{\mathrm{0}}}\mathfrak{I}\left({S}_{\mathrm{hh}}^{\mathrm{fwd}}\right)N\left(D\right)\mathrm{d}D\\ \\ ={\mathrm{10}}^{\mathrm{3}}\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}({K}_{\mathrm{11}}{K}_{\mathrm{12}})N\left(D\right)\mathrm{d}D\\ ={\mathrm{10}}^{\mathrm{3}}\left(\u2329{K}_{\mathrm{11}}\u232a\u2329{K}_{\mathrm{12}}\u232a\right)\end{array}\text{(D6)}& {\displaystyle}{a}_{\mathrm{v}}={\mathrm{10}}^{\mathrm{3}}\left(\u2329{K}_{\mathrm{11}}\u232a+\u2329{K}_{\mathrm{12}}\u232a\right).\end{array}$$Again, 10^{−3} in Eq. (D5) serves as the coefficient (in unit the conversion from mm^{2} m^{−3} to km^{−1}). σ_{ext,h} indicates horizontal extinction cross section of a given particle (i.e., ${\mathit{\sigma}}_{\mathrm{ext},\mathrm{h}}=\frac{\mathrm{4}\mathit{\pi}}{{k}_{\mathrm{0}}}\mathfrak{I}\left({S}_{\mathrm{hh}}^{\mathrm{fwd}}\right))$. Similarly, vertical attenuation coefficient can be derived in Eq. (D6).

The total differentialphase shift upon backscattering δ_{hv} in units of degrees (hereafter deg.) is represented as shown below:
$$\begin{array}{}\text{(D7)}& \begin{array}{l}{\mathit{\delta}}_{\mathrm{hv}}=\frac{\mathrm{180}}{\mathit{\pi}}\mathrm{\angle}\left(\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{S}_{\mathrm{hh}}{S}_{\mathrm{vv}}^{\ast}N\left(D\right)\mathrm{d}D\right)\\ =\frac{\mathrm{180}}{\mathit{\pi}}\mathrm{\angle}\left(\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}\left[\mathrm{0.5}({Z}_{\mathrm{33}}+{Z}_{\mathrm{44}})+\mathrm{0.5}({Z}_{\mathrm{43}}{Z}_{\mathrm{34}})i\right]N\left(D\right)\mathrm{d}D\right)\\ =\frac{\mathrm{180}}{\mathit{\pi}}\mathrm{\angle}\left(\left[\u2329{Z}_{\mathrm{33}}\u232a+\u2329{Z}_{\mathrm{44}}\u232a\right]+\left[\u2329{Z}_{\mathrm{43}}\u232a\u2329{Z}_{\mathrm{34}}\u232a\right]i\right)\end{array}.\end{array}$$The notation “∠” in Eq. (D7) indicates the phase of the complex value. The coefficient $\mathrm{180}/\mathit{\pi}$ is used to convert the unit of the result from radii to deg.

Copolar correlation coefficient ρ_{hv} (dimensionless):
$$\begin{array}{}\text{(D8)}& \begin{array}{l}{\mathit{\rho}}_{\mathrm{hv}}=\frac{\left\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{S}_{\mathrm{hh}}{S}_{\mathrm{vv}}^{\ast}N\left(D\right)\mathrm{d}D\right}{\sqrt{\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{\left{S}_{\mathrm{hh}}\right}^{\mathrm{2}}N\left(D\right)\mathrm{d}D\cdot \underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{\left{S}_{\mathrm{vv}}\right}^{\mathrm{2}}N\left(D\right)\mathrm{d}D}}\\ =\frac{\left\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}\left[\mathrm{0.5}({Z}_{\mathrm{33}}+{Z}_{\mathrm{44}})+\mathrm{0.5}({Z}_{\mathrm{43}}{Z}_{\mathrm{34}})i\right]N\left(D\right)\mathrm{d}D\right}{\sqrt{{\scriptscriptstyle \begin{array}{c}{\int}_{\mathrm{0}}^{\mathrm{\infty}}\mathrm{0.5}({Z}_{\mathrm{11}}{Z}_{\mathrm{12}}{Z}_{\mathrm{21}}+{Z}_{\mathrm{22}})N\left(D\right)\mathrm{d}D\\ \cdot {\int}_{\mathrm{0}}^{\mathrm{\infty}}\mathrm{0.5}({Z}_{\mathrm{11}}+{Z}_{\mathrm{12}}+{Z}_{\mathrm{21}}+{Z}_{\mathrm{22}})N\left(D\right)\mathrm{d}D\end{array}}}}\\ =\frac{\left\left(\u2329{Z}_{\mathrm{33}}\u232a+\u2329{Z}_{\mathrm{44}}\u232a\right)+\left(\u2329{Z}_{\mathrm{43}}\u232a\u2329{Z}_{\mathrm{34}}\u232a\right)i\right}{\sqrt{{\scriptscriptstyle \begin{array}{c}\left(\u2329{Z}_{\mathrm{11}}\u232a\u2329{Z}_{\mathrm{12}}\u232a\u2329{Z}_{\mathrm{21}}\u232a+\u2329{Z}_{\mathrm{22}}\u232a\right)\\ \cdot \left(\u2329{Z}_{\mathrm{11}}\u232a+\u2329{Z}_{\mathrm{12}}\u232a+\u2329{Z}_{\mathrm{21}}\u232a+\u2329{Z}_{\mathrm{22}}\u232a\right)\end{array}}}}\end{array}.\end{array}$$In Eq. (D8), the copolar correlation coefficient ρ_{hv} is the amplitude of complex copolar correlation coefficient $\stackrel{\mathrm{\u0303}}{{\mathit{\rho}}_{\mathrm{hv}}}$, whose phase is the total differentialphase shift upon backscattering in Eq. (D7).
Please note that the summation over hydrometeor types was omitted in Eqs. (D1)–(D8) for the sake of clarity. Readers can easily obtain the more complicated, real expressions of radar variables by applying extra summations over hydrometeor types for 〈Z〉 and 〈K〉 elements before carrying out the calculations.
The aforementioned radar variables ${z}_{\mathrm{h}/\mathrm{v}}$, ${a}_{\mathrm{h}/\mathrm{v}}\left({\mathit{\alpha}}_{\mathrm{h}/\mathrm{v}}\right)$, z_{dr}, K_{DP}, δ_{hv}, and ρ_{hv} are often referred to as intrinsic radar variables determined by the atmosphere and hydrometeor conditions in local radar gates. Under the assumption of firstorder multiplescattering model, although the wave is scattered only once before it is received, the twoway propagation effects such as attenuation and phase shift are taken into account using the wave number of the effective medium composed of air and hydrometeors along the beam trajectory (Zhang, 2016). We can derive the observable radar variables as shown below:
In Eqs. (D9) and (D10), r_{g} is the range of the radar gate, and the variable of integration r is the range along the radar beam trajectory. ${z}_{\mathrm{h}/\mathrm{v}}^{\prime}$ marked with a prime is the horizontal/vertical observable reflectivity of the radar gate attenuated on the way from the transmitter to the particle and on the way back from the particle to the receiver. Φ_{DP} is the total phase shift interpreted as the differentialphase shift upon twoway propagation plus the differential phase shift upon backscattering.
Equations (D11) and (D12) give the observable horizontal/vertical reflectivity ${Z}_{\mathrm{H}/\mathrm{V}}^{\prime}$ (in units of dBZ) and differential reflectivity ${Z}_{\mathrm{DR}}^{\prime}$ (in units of dB).
Codes of the forward radar operator ZJUAERO V0.5 and the packaged Conda environment and the user manual are available on Zenodo (https://doi.org/10.5281/zenodo.11307123; Xie et al., 2024a). The database of the scattering properties (i.e., the LUTs) is also released with the software package.
Two cases of this forward radar operator are presented in the user manual of ZJUAERO to demonstrate its usage for spaceborne and groundbased radar, respectively. The NWP model grid data and the radar observation products for those two demonstration cases are also available on Zenodo (https://doi.org/10.5281/zenodo.11307206; Xie et al., 2024b).
HX performed the coding and visualization of the LUTs and designed the case study experiments. LB and WH supervised this study. All authors contributed to the writing of the paper.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This research has been supported by the National Natural Science Foundation of China (grant nos. U2342213 and 42022038) and National Key Research and Development Program of China (grant no. 2022YFC3004004). We are grateful to NASA and JAXA for providing the L2A data of the dualfrequency spaceborne precipitation radar GPM/DPR, which serve as a good validation dataset for the present spaceborne forward radar operator.
This research has been supported by the National Natural Science Foundation of China (grant no. 42022038) and the National Key Research and Development Program of China (grant no. 2022YFC3004004).
This paper was edited by Yuefei Zeng and reviewed by two anonymous referees.
Abel, S. and Boutle, I.: An improved representation of the raindrop size distribution for singlemoment microphysics schemes, Q. J. Roy. Meteor. Soc., 138, 2151–2162, 2012.
Bárdossy, A. and Das, T.: Influence of rainfall observation network on model calibration and application, Hydrol. Earth Syst. Sci., 12, 77–89, https://doi.org/10.5194/hess12772008, 2008.
Bi, L. and Yang, P.: Accurate simulation of the optical properties of atmospheric ice crystals with the invariant imbedding Tmatrix method, J. Quant. Spectrosc. Ra., 138, 17–35, 2014.
Bi, L., Yang, P., Kattawar, G. W., and Mishchenko, M. I.: Efficient implementation of the invariant imbedding Tmatrix method and the separation of variables method applied to large nonspherical inhomogeneous particles, J. Quant. Spectrosc. Ra., 116, 169–183, 2013.
Bi, L., Wang, Z., Han, W., Li, W., and Zhang, X.: Computation of optical properties of coreshell superspheroids using a GPU implementation of the invariant imbedding Tmatrix method, Frontiers in Remote Sensing, 3, 35, https://doi.org/10.3389/frsen.2022.903312, 2022.
Brandes, E. A., Zhang, G., and Vivekanandan, J.: Experiments in rainfall estimation with a polarimetric radar in a subtropical environment, J. Appl. Meteorol., 41, 674–685, 2002.
Cattoën, C., Robertson, D., Bennett, J., Wang, Q., and CareySmith, T.: Calibrating hourly precipitation forecasts with daily observations, J. Hydrometeorol., 21, 1655–1673, 2020.
Caumont, O., Ducrocq, V., Delrieu, G., Gosset, M., Pinty, J.P., Du Chatelet, J. P., Andrieu, H., Lemaître, Y., and Scialom, G.: A radar simulator for highresolution nonhydrostatic models, J. Atmos. Ocean. Tech., 23, 1049–1067, 2006.
Caumont, O., Ducrocq, V., Wattrelot, É., Jaubert, G., and PradierVabre, S.: 1D+3DVar assimilation of radar reflectivity data: A proof of concept, Tellus A, 62, 173–187, 2010.
Chandrasekar, V.: Polarimetric Doppler weather radar [electronic resource]: principles and applications, Cambridge University Press, ISBN 9780511016073, 2001.
Chen, D., Xue, J., Yang, X., Zhang, H., Shen, X., Hu, J., Wang, Y., Ji, L., and Chen, J.: New generation of multiscale NWP system (GRAPES): general scientific design, Chinese Sci. Bull., 53, 3433–3445, https://doi.org/10.1007/s114340080494z, 2008.
Chuang, C. C. and Beard, K. V.: A numerical model for the equilibrium shape of electrified raindrops, J. Atmos. Sci., 47, 1374–1389, 1990.
Ekelund, R., Eriksson, P., and Kahnert, M.: Microwave singlescattering properties of nonspheroidal raindrops, Atmos. Meas. Tech., 13, 6933–6944, https://doi.org/10.5194/amt1369332020, 2020.
Ellison, W.: Permittivity of pure water, at standard atmospheric pressure, over the frequency range −25 THz and the temperature range −100 °C, J. Phys. Chem. Ref. Data, 36, 1–18, https://doi.org/10.1063/1.2360986, 2007.
Eriksson, P., Ekelund, R., Mendrok, J., Brath, M., Lemke, O., and Buehler, S. A.: A general database of hydrometeor single scattering properties at microwave and submillimetre wavelengths, Earth Syst. Sci. Data, 10, 1301–1326, https://doi.org/10.5194/essd1013012018, 2018.
Fabry, F. and Zawadzki, I.: Longterm radar observations of the melting layer of precipitation and their interpretation, J. Atmos. Sci., 52, 838–851, 1995.
Field, P., Hogan, R., Brown, P., Illingworth, A., Choularton, T., and Cotton, R.: Parametrization of iceparticle size distributions for midlatitude stratiform cloud, Q. J. Roy. Meteor. Soc., 131, 1997–2017, 2005.
Garrett, T. J., Fallgatter, C., Shkurko, K., and Howlett, D.: Fall speed measurement and highresolution multiangle photography of hydrometeors in free fall, Atmos. Meas. Tech., 5, 2625–2633, https://doi.org/10.5194/amt526252012, 2012.
Garrett, T. J., Yuter, S. E., Fallgatter, C., Shkurko, K., Rhodes, S. R., and Endries, J. L.: Orientations and aspect ratios of falling snow, Geophys. Res. Lett., 42, 4617–4622, 2015.
Geer, A. J. and Baordo, F.: Improved scattering radiative transfer for frozen hydrometeors at microwave frequencies, Atmos. Meas. Tech., 7, 1839–1860, https://doi.org/10.5194/amt718392014, 2014.
Geer, A. J., Bauer, P., Lonitz, K., Barlakas, V., Eriksson, P., Mendrok, J., Doherty, A., Hocking, J., and Chambon, P.: Bulk hydrometeor optical properties for microwave and submillimetre radiative transfer in RTTOVSCATT v13.0, Geosci. Model Dev., 14, 7497–7526, https://doi.org/10.5194/gmd1474972021, 2021.
Hong, S.Y. and Lim, J.O. J.: The WRF singlemoment 6class microphysics scheme (WSM6), AsiaPac. J. Atmos. Sci., 42, 129–151, 2006.
Hubbert, J., Ellis, S., Dixon, M., and Meymaris, G.: Modeling, error analysis, and evaluation of dualpolarization variables obtained from simultaneous horizontal and vertical polarization transmit radar. Part I: Modeling and antenna errors, J. Atmos. Ocean. Tech., 27, 1583–1598, 2010.
Iguchi, T.: Dualfrequency precipitation radar (DPR) on the global precipitation measurement (GPM) mission's core observatory, Satellite Precipitation Measurement, 1, 183–192, ISBN 9783030245689, https://link.springer.com/chapter/10.1007/9783030245689_11 (last access: 8 September 2023), 2020.
Iguchi, T., Hanado, H., Takahashi, N., Kobayashi, S., and Satoh, S.: The dualfrequency precipitation radar for the GPM core satellite, IGARSS 2003, in: 2003 IEEE International Geoscience and Remote Sensing Symposium. Proceedings (IEEE Cat. No. 03CH37477), Toulouse, France, 21–25 July 2003, 1698–1700, https://doi.org/10.1109/IGARSS.2003.1294221, 2003.
Iguchi, T., Seto, S., Meneghini, R., Yoshida, N., Awaka, J., Le, M., Chandrasekar, V., and Kubota, T.: GPM/DPR level2 algorithm theoretical basis document, NASA Goddard Space Flight Center, https://gpm.nasa.gov/resources/documents/gpmdprlevel2algorithmtheoreticalbasisdocumentatbd (last access: 10 April 2023), 2010.
Jung, Y., Zhang, G., and Xue, M.: Assimilation of simulated polarimetric radar data for a convective storm using the ensemble Kalman filter. Part I: Observation operators for reflectivity and polarimetric variables, Mon. Weather Rev., 136, 2228–2245, 2008.
Kobayashi, T. and Adachi, A.: Measurements of raindrop breakup by using UHF wind profilers, Geophys. Res. Lett., 28, 4071–4074, 2001.
Le, M. and Chandrasekar, V.: Precipitation type classification method for dualfrequency precipitation radar (DPR) onboard the GPM, IEEE T. Geosci. Remote, 51, 1784–1790, 2012.
Liao, L., Meneghini, R., and Tokay, A.: Uncertainties of GPM DPR rain estimates caused by DSD parameterizations, J. Appl. Meteorol. Clim., 53, 2524–2537, 2014.
Liebe, H. J., Hufford, G. A., and Manabe, T.: A model for the complex permittivity of water at frequencies below 1 THz, Int. J. Infrared Milli., 12, 659–675, 1991.
Ma, Z., Zhao, C., Gong, J., Zhang, J., Li, Z., Sun, J., Liu, Y., Chen, J., and Jiang, Q.: Spinup characteristics with three types of initial fields and the restart effects on forecast accuracy in the GRAPES global forecast system, Geosci. Model Dev., 14, 205–221, https://doi.org/10.5194/gmd142052021, 2021.
Marshall, J. S. and Palmer, W. M. K.: The distribution of raindrops with size, J. Meteorol., 5, 165–166, 1948.
Mishchenko, M. I.: Calculation of the amplitude matrix for a nonspherical particle in a fixed orientation, Appl. Optics, 39, 1026–1031, 2000.
Mishchenko, M. I.: Electromagnetic scattering by particles and particle groups: an introduction, Cambridge University Press, ISBN 9781139019064, 2014.
Mishchenko, M. I. and Travis, L. D.: Tmatrix computations of light scattering by large spheroidal particles, Opt. Commun., 109, 16–21, 1994.
Morrison, H. and Milbrandt, J. A.: Parameterization of cloud microphysics based on the prediction of bulk ice particle properties. Part I: Scheme description and idealized tests, J. Atmos. Sci., 72, 287–311, 2015.
Oue, M., Tatarevic, A., Kollias, P., Wang, D., Yu, K., and Vogelmann, A. M.: The Cloudresolving model Radar SIMulator (CRSIM) Version 3.3: description and applications of a virtual observatory, Geosci. Model Dev., 13, 1975–1998, https://doi.org/10.5194/gmd1319752020, 2020.
Pruppacher, H. R., Klett, J. D., and Wang, P. K.: Microphysics of Clouds and Precipitation, Aerosol Sci. Tech., 28, 381–382, https://doi.org/10.1080/02786829808965531, 1998.
Radhakrishna, B.: Raindrop size distribution (DSD) during the passage of tropical cyclone Nivar: effect of measuring principle and wind on DSDs and retrieved rain integral and polarimetric parameters from impact and laser disdrometers, Atmos. Meas. Tech., 15, 6705–6722, https://doi.org/10.5194/amt1567052022, 2022.
Rose, C. and Chandrasekar, V.: A GPM dualfrequency retrieval algorithm: DSD profileoptimization method, J. Atmos. Ocean. Tech., 23, 1372–1383, 2006.
Ryzhkov, A., Pinsky, M., Pokrovsky, A., and Khain, A.: Polarimetric radar observation operator for a cloud model with spectral microphysics, J. Appl. Meteorol. Clim., 50, 873–894, 2011.
Ryzhkov, A. V.: The impact of beam broadening on the quality of radar polarimetric data, J. Atmos. Ocean. Tech., 24, 729–744, 2007.
Shen, X., Su, Y., Zhang, H., and Hu, J.: New Version of the CMAGFS Dynamical Core Based on the Predictor–Corrector Time Integration Scheme, Journal of Meteorological Research, 37, 273–285, 2023.
Shrestha, P., Mendrok, J., and Brunner, D.: Aerosol characteristics and polarimetric signatures for a deep convective storm over the northwestern part of Europe – modeling and observations, Atmos. Chem. Phys., 22, 14095–14117, https://doi.org/10.5194/acp22140952022, 2022.
Skamarock, W. C., Klemp, J. B., Dudhia, J., Gill, D. O., Liu, Z., Berner, J., Wang, W., Powers, J. G., Duda, M. G., and Barker, D. M.: A Description of the Advanced Research WRF Model Version 4.3, No. NCAR/TN556+STR), https://doi.org/10.5065/1dfh6p97, 2019.
Thompson, G., Field, P. R., Rasmussen, R. M., and Hall, W. D.: Explicit forecasts of winter precipitation using an improved bulk microphysics scheme. Part II: Implementation of a new snow parameterization, Mon. Weather Rev., 136, 5095–5115, 2008.
Thurai, M., Huang, G., Bringi, V., Randeu, W., and Schönhuber, M.: Drop shapes, model comparisons, and calculations of polarimetric radar parameters in rain, J. Atmos. Ocean. Tech., 24, 1019–1032, 2007.
Trömel, S., Simmer, C., Blahak, U., Blanke, A., Doktorowski, S., Ewald, F., Frech, M., Gergely, M., Hagen, M., Janjic, T., KalesseLos, H., Kneifel, S., Knote, C., Mendrok, J., Moser, M., Köcher, G., Mühlbauer, K., Myagkov, A., Pejcic, V., Seifert, P., Shrestha, P., Teisseire, A., von Terzi, L., Tetoni, E., Vogl, T., Voigt, C., Zeng, Y., Zinner, T., and Quaas, J.: Overview: Fusion of radar polarimetry and numerical atmospheric modelling towards an improved understanding of cloud and precipitation processes, Atmos. Chem. Phys., 21, 17291–17314, https://doi.org/10.5194/acp21172912021, 2021.
Van de Hulst, H.: Light scattering by small particles, Dover Publications, ISBN 9780486642284, 1981.
Walters, D. N., Best, M. J., Bushell, A. C., Copsey, D., Edwards, J. M., Falloon, P. D., Harris, C. M., Lock, A. P., Manners, J. C., Morcrette, C. J., Roberts, M. J., Stratton, R. A., Webster, S., Wilkinson, J. M., Willett, M. R., Boutle, I. A., Earnshaw, P. D., Hill, P. G., MacLachlan, C., Martin, G. M., MoufoumaOkia, W., Palmer, M. D., Petch, J. C., Rooney, G. G., Scaife, A. A., and Williams, K. D.: The Met Office Unified Model Global Atmosphere 3.0/3.1 and JULES Global Land 3.0/3.1 configurations, Geosci. Model Dev., 4, 919–941, https://doi.org/10.5194/gmd49192011, 2011.
Wang, J., Dong, X., Xi, B., and Heymsfield, A. J.: Investigation of liquid cloud microphysical properties of deep convective systems: 1. Parameterization raindrop size distribution and its application for stratiform rain estimation, J. Geophys. Res.Atmos., 121, 11637–11651, https://doi.org/10.1029/2018JD028727, 2016.
Wang, S. and Liu, Z.: A radar reflectivity operator with icephase hydrometeors for variational data assimilation (version 1.0) and its evaluation with real radar data, Geosci. Model Dev., 12, 4031–4051, https://doi.org/10.5194/gmd1240312019, 2019.
Wang, Z., Bi, L., and Kong, S.: Flexible implementation of the particle shape and internal inhomogeneity in the invariant imbedding Tmatrix method, Opt. Express, 31, 29427–29439, 2023.
Warren, R. A., Protat, A., Siems, S. T., Ramsay, H. A., Louf, V., Manton, M. J., and Kane, T. A.: Calibrating groundbased radars against TRMM and GPM, J. Atmos. Ocean. Tech., 35, 323–346, 2018.
Wolfensberger, D. and Berne, A.: From model to radar variables: a new forward polarimetric radar operator for COSMO, Atmos. Meas. Tech., 11, 3883–3916, https://doi.org/10.5194/amt1138832018, 2018.
Xie, H., Bi, L., Wang, Z., and Han, W.: An Accurate and Efficient Radar Operator Designed for CMAGFS/MESO with Capability of Simulating Nonspherical Hydrometeors, Zenodo [code], https://doi.org/10.5281/zenodo.11307123, 2024a.
Xie, H., Bi, L., and Han, W.: Example Case Data of Forward Radar Operator ZJUAERO Release V0.6.1, Zenodo [data set], https://doi.org/10.5281/zenodo.11307206, 2024b.
Zawadzki, I., Szyrmer, W., Bell, C., and Fabry, F.: Modeling of the melting layer. Part III: The density effect, J. Atmos. Sci., 62, 3705–3723, 2005.
Zeng, Y., Blahak, U., Neuper, M., and Jerger, D.: Radar beam tracing methods based on atmospheric refractive index, J. Atmos. Ocean. Tech., 31, 2650–2670, 2014.
Zeng, Y., Blahak, U., and Jerger, D.: An efficient modular volumescanning radar forward operator for NWP models: description and coupling to the COSMO model, Q. J. Roy. Meteor. Soc., 142, 3234–3256, 2016.
Zhang, G.: Weather radar polarimetry, Crc Press, ISBN 9781315374666, https://doi.org/10.1201/9781315374666, 2016.
Zhang, G., Gao, J., and Du, M.: Parameterized forward operators for simulation and assimilation of polarimetric radar data with numerical weather predictions, Adv. Atmos. Sci., 38, 737–754, 2021.
Zhang, P., Gu, S., Chen, L., Shang, J., Lin, M., Zhu, A., Yin, H., Wu, Q., Shou, Y., and Sun, F.: FY3G satellite instruments and precipitation products: first report of China's Fengyun rainfall mission inorbit, Journal of Remote Sensing, 3, 0097, https://doi.org/10.34133/remotesensing.0097, 2023.
Angle brackets are only used to indicate the integration over the particle size distribution in this study. Integration over azimuthal orientation and shape distributions are indicated by overlines.
In the weather radar community, the particle diameter D usually refers to equalvolumesphere diameter for a liquid hydrometeor, while D is often regarded as maximum dimension when describing the solid and melting types of hydrometeors.
Overlines over elements of Z and K are omitted for brevity for level A database elements elsewhere in this paper.
Two overlines over elements of Z,K, 〈Z〉, and 〈K〉 are omitted for brevity for levels B/C database elements elsewhere in this paper and in formulas of polarimetric radar variables.
The oneway specific attenuation (in dB scale) ${\mathit{\alpha}}_{\mathrm{h}/\mathrm{v}}$ (in units of dB km^{−1}), which relates to the linearscale attenuation coefficient by ${\mathit{\alpha}}_{\mathrm{h}/\mathrm{v}}=\mathrm{10}{\mathrm{log}}_{\mathrm{10}}\left(e\right)\cdot {a}_{\mathrm{h}/\mathrm{v}}=\mathrm{4.343}\cdot {a}_{\mathrm{h}/\mathrm{v}}$, is also used in many studies of weather radar that consider the change in the logarithm base from e to 10.
 Abstract
 Introduction
 General descriptions of ZJUAERO
 Database of hydrometeor optical properties
 Case studies
 Summary and ongoing tasks
 Appendix A: Specifications on spaceborne radar trajectory solver
 Appendix B: Appendix B: grids of LUT dimensions in ZJUAERO
 Appendix C: Forward scattering alignment (FSA) and backwardscattering alignment (BSA)
 Appendix D: Calculations of radar variables
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 General descriptions of ZJUAERO
 Database of hydrometeor optical properties
 Case studies
 Summary and ongoing tasks
 Appendix A: Specifications on spaceborne radar trajectory solver
 Appendix B: Appendix B: grids of LUT dimensions in ZJUAERO
 Appendix C: Forward scattering alignment (FSA) and backwardscattering alignment (BSA)
 Appendix D: Calculations of radar variables
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References