Quantifying the Cloud Water Resource: Methods Based on Observational Diagnosis and Cloud Model Simulation

Based on the concepts of cloud water resource (CWR) and related variables proposed in the first part of this study, this paper provides details of two methods to quantify the CWR. One is diagnostic quantification (CWR-DQ) based on satellite observations, precipitation products, and atmospheric reanalysis data; and the other is numerical quantification (CWR-NQ) based on a cloud resolving model developed at the Chinese Academy of Meteorological Sciences (CAMS). The two methods are applied to quantify the CWR in April and August 2017 over North China, and the results are evaluated against all available observations. Main results are as follows. (1) For the CWR-DQ approach, reference cloud profiles are firstly derived based on the CloudSat/CALIPSO joint satellite observations for 2007–2010. The NCEP/NCAR reanalysis data in 2000–2017 are then employed to produce three-dimensional cloud fields. The budget/balance equations of atmospheric water substance are lastly used, together with precipitation observations, to retrieve CWR and related variables. It is found that the distribution and vertical structure of clouds obtained by the diagnostic method are consistent with observations. (2) For the CWR-NQ approach, it assumes that the cloud resolving model is able to describe the cloud microphysical processes completely and precisely, from which four-dimensional distributions of atmospheric water vapor, hydrometeors, and wind fields can be obtained. The data are then employed to quantify the CWR and related terms/quantities. After one-month continuous integration, the mass of atmospheric water substance becomes conserved, and the tempospatial distributions of water vapor, hydrometeors/cloud water, and precipitation are consistent with observations. (3) Diagnostic values of the difference in the transition between hydrometeors and water vapor (Cvh − Chv) and the surface evaporation (Es) are well consistent with their numerical values. (4) Correlation and bias analyses show that the diagnostic CWR contributors are well correlated with observations, and match their numerical counterparts as well, indicating that the CWR-NQ and CWR-DQ methods are reasonable. (5) Underestimation of water vapor converted from hydrometeors (Chv) is a shortcoming of the CWR-DQ method, which may be rectified by numerical quantification results or by use of advanced observations on higher spatiotemporal resolutions.


Introduction
China faces extreme shortage of freshwater resources, and the per capita water resources are less than a quarter of the world average. In the change of water substance in the atmosphere, hydrometeors (cloud water) play an important role. The cloud water resource (CWR) refers to the hydrometeors that participate in the atmospheric wa-ter change in a certain area and time period but have not formed ground precipitation, and thus can be utilized/exploited (Cai, 2013;Zhou et al., 2020a). At present, artificial precipitation (snow) enhancement is the main technique to exploit the CWR. Scientific quantification and understanding of spatiotemporal distribution characteristics of the CWR are very important for weather modification planning and reasonable exploitation of the CWR.
In the water change process in the atmosphere, hydrometeors are inseparable from water vapor. The water mass in the atmosphere must remain conserved. Zhou et al. (2020a) defined 16 CWR contributors (including 12 independent ones) and 12 characteristic variables associated with precipitation and CWR based on the water substance balance equations in the atmosphere, and proposed the algorithms for the calculation of these physical variables. Among them, the CWR contributors include the mass of hydrometeors at a specific time, advection of hydrometeors, and conversion between water vapor and hydrometeors in a certain period, as well as surface evaporation and precipitation. Based on the above quantities, the CWR and related characteristic variables, such as the gross mass of hydrometeors in the atmosphere, precipitation efficiency, renewal time, etc. can be quantitatively calculated. Same items can be quantified for water vapor, as well as the atmospheric water substance (which is the summation of values for hydrometeors and water vapor).
From the perspective of observations, in order to diagnose and quantify the CWR in a specific area during a certain time period, the data that are most needed are the three-dimensional (3D) time-varying mass of atmospheric water substance and wind fields. The atmospheric water vapor and winds are continuously distributed in space and time, which can be obtained from observation-based atmospheric reanalysis data (such as the NCEP/NCAR reanalysis products). Surface precipitation can also be directly observed. To obtain the state variables of atmospheric hydrometeors and their advection requires observations of 3D time-varying clouds (including spatial distribution of cloud cover and cloud water content). However, clouds are spatially and temporally discontinuous, and systematic cloud observations are not available. Therefore, how to obtain 3D cloud data is a critical issue in the CWR diagnostic quantification.
The spatial distribution of clouds is unstable and discontinuous. Ground-based observation of cloud shape and cloud cover is the main source of cloud information in earlier times. Du et al. (2012) investigated the characteristics of local CWR based on cloud cover observations at local stations. Zhang (2002) examined research papers and sorted out relevant data used for cloud water content calculation, including water content in various types of cloud, cloud cover area, cloud depth, cloud volume, total cloud water content, etc. Using this dataset, Li et al. (2010) calculated the gross mass of cloud water in Jiangxi Province of China.
Since the 1980s, development of meteorological satellites has provided an effective method for large-scale cloud observations. At present, various cloud products released by the International Satellite Cloud Climate Program (ISCCP; Rossow and Schiffer, 1991), Clouds and the Earth's Radiant Energy System (CERES; Zhang et al., 2012), Moderate Resolution Imaging Spectroradiometer (MODIS; Ackerman, 2008;Wang et al., 2011), and so on, have been widely applied. These cloud products mainly include cloud cover, cloud water path, cloud optical thickness, etc. Many researchers in China have used the above satellite products to study spatiotemporal distribution characteristics of cloud cover and cloud water path in different regions (Chen et al., 2005;Zhang et al., 2006;Li X. Y. et al., 2008;Peng et al., 2010;Li M. et al., 2015;Bo et al., 2016). In these studies, the cloud water path, i.e., the integrated mass of atmospheric hydrometeors, is generally regarded as the CWR.
According to the definition of CWR proposed by Zhou et al. (2020a), the state quantity of hydrometeors in the previous studies is only a small part of the CWR. The calculation of the CWR should also consider advection or horizontal transport of atmospheric hydrometeors and other sources/sinks terms newly generated/transformed by physical processes such as cloud condensation and cloud evaporation. However, due to the lack of largescale temporally and spatially continuous 3D cloud observations, few observational studies have been conducted on the characteristics of advection/transport of atmospheric hydrometeors. In addition, the difficulty in observation of cloud condensation and evaporation makes it even harder for the CWR study.
Due to the imperative need to understand the vertical structure of clouds, researchers have been exploring how to use temperature and humidity soundings to identify cloudy areas (Essenwanger and Haggard, 1962;Poore et al., 1995). Theoretically, the relative humidity (RH) is around 100% in a stable cloudy area; specifically, it is closer to 100% within liquid-phase clouds and greater than 100% in mixed-phrase or ice clouds. Therefore, RH can be used to determine cloud cover area. However, observed RH in actual clouds can seldom reach 100%, which is related to errors caused by observational instruments. Current data of 3D humidity distributions are mainly obtained from soundings, and retrievals based on remote sensing are verified against sounding observa-tions. Therefore, the sounding observation errors affect the global and regional humidity observations and reanalysis products, and instrument errors must be considered in the application of such data.
Based on multi-case analyses of sounding data at onesecond intervals and CloudSat-borne cloud radar observations, Zhou and Ou (2010) proposed an RH threshold (RH c ) that can be used to analyze vertical structure of clouds. With the L-band sounding data over China for two years and the spatiotemporally matched CloudSat satellite observations, Cai et al. (2014) calculated skill scores of BS (Bias Score) and TS (Threat Score) to evaluate RH c at different altitudes, and proposed a set of algorithms to identify single-point cloud vertical structure mainly based on L-band second-interval sounding data. On this basis, Cai et al. (2015) established a statistical relationship between CloudSat observed clouds and the RH extracted from ECMWF products, obtained RH c at various altitudes for cloud identification, and eventually retrieved the 3D spatial distribution of cloud cover. This is a preliminary exploration of 3D cloud cover diagnosis.
According to Winker and Pelon (2003) and Winker et al. (2007), the CALIPSO satellite can observe both weak water vapor condensation layers at high altitude and ice clouds with small optical thickness. By combining CALIPSO with CloudSat, a more complete cloud vertical structure can be detected. Compared with CloudSat observation only, the CloudSat/CALIPSO joint observations have a higher accuracy in identifying terrestrial cloud areas, and the results are closer to CERES observations (Tang, 2018). Apparently, our diagnostic method for cloud cover needs further optimization. Meanwhile, it is also necessary to develop diagnostic methods and products for 3D distributions of cloud water content for estimation of the state quantity and advection of atmospheric hydrometeors, which, combined with wind field and precipitation observation, can be used to calculate the CWR and its contributors as well as related characteristic variables.
Relevant numerical case studies have also been conducted in China. For example,  investigated the atmospheric water substance budget and precipitation generation mechanism for a stratiform cloud system over Henan Province of central China using the MM5 modeling system. They analyzed the water vapor transport and evaluated the budget of water vapor and hydrometeors over the study area. Chen et al. (2011) used the GRAPES (Global/Regional Assimilation and PrEdiction System) model and simulated a precipitation process that occurred in Chongqing. They presented the distributions of water vapor and hydrometeors, and ana-lyzed the cloud microphysical processes and precipitation efficiency. Using the Weather Research and Forecasting (WRF) model dynamic core coupled with the Chinese Academy of Meteorological Science (CAMS) cloud microphysics scheme, Tao et al. (2015) simulated the structure and moisture budget of a stratocumulus cloud system over Beijing. They estimated the moisture budget, and calculated water vapor condensation rate, sublimation rate, and PE for hydrometeors and atmospheric water vapor. Nonetheless, the above studies are all based on single-case simulations. Neither the properties of cloud water over long periods were discussed, nor the quantification of CWR was conducted.
According to Zhou et al. (2020a), the CWR cannot be observed directly, but can be quantified through an observation-based diagnostic method and through numerical simulations. In this study, a diagnostic method and a numerical simulation method are firstly established for the calculation of CWR based on the concepts and algorithms of CWR proposed by Zhou et al. (2020a). The diagnostic method relies firstly on satellite observations and atmospheric reanalysis data to retrieve cloud profiles and 3D cloud fields, then on the water substance budget/balance equations in the atmosphere to derive the CWR and related quantities. The numerical quantification of CWR is based on outputs from a cloud resolving model with a comprehensive explicit cloud microphysics scheme developed by the CAMS. In this paper, the two methods are introduced and then implemented to quantify the CWR in North China in two typical months, i.e., April and August 2017. Finally, the results from the two sets of methods are compared with each other and with observations to evaluate the performance of the methods.

Theoretical basis and the key issue
According to the definition and calculation algorithms of the CWR proposed by Zhou et al. (2020a), the basis of the CWR quantification for any period and region is the water substance balance equations in the atmosphere. Hydrometeors (or cloud water), water vapor, and the water substance in the atmosphere comply with the water budget/balance equations: where the 16 quantities/terms (including 12 independent terms) related to the hydrometeors, water vapor, and water substance (sum of hydrometeors and water vapor) in the atmosphere, indicated by the subscript h, v, and w, respectively, are termed as the contributors of the CWR. The quantities M x1 and M x2 (the subscript x can be replaced by h, v, and w) are the mass of hydrometeors, or water vapor, or water substance in the atmosphere at the initial and end time of the study period, termed as instant state variables; Q xi and Q xo are the inflow and outflow of hydrometeors, or water vapor, or water substance in the atmosphere during the study period, termed as advection variables; and C vh is the mass of hydrometeors converted from water vapor through condensation (deposition included) while C hv is the mass of water vapor converted from hydrometeors through evaporation (sublimation included) during the study period. Together with surface evaporation (E s ) and precipitation (P s ), C vh and C hv are termed as source/sink terms. In Eq. (1), the income terms that maintain or bring in water substance in a specific area over a certain time period are listed on the left-hand side while the expenditure terms are on the right-hand side. The sum of the income terms is equal to that of the expenditure terms, which is defined as the gross mass of the hydrometeors, water vapor, and water substance (GM h , GM v , and GM w ). Among the GM h , those hydrometeors that have not formed surface precipitation but remain in the air are the CWR. According to Eq. (1), the CWR contributors in any period and region can be quantitatively calculated based on appropriate data, and the CWR is derived as, Furthermore, the CWR-related characteristic variables such as precipitation efficiency of hydrometeors, water vapor, and total water substance (i.e., PE h , PE v , and PE w ), condensation efficiency of water vapor (CE v ), mean mass of hydrometeors and water vapor (MM h and MM v ), as well as their renewal time (RT h and RT v ), are also derived. Together with GM h , GM v , and GM w , these variables represent the characteristics of the CWR in a specific region over a certain time period. For detailed formulas, refer to Eqs. (6)-(13) in Zhou et al. (2020a). As indicated in Zhou et al. (2020a), the size of the study area and the duration of the time integral/average affect the calculation results of CWR and related variables. The CWR in a large area mainly depends on cloud physical processes, followed by advection along the boundary.
Fundamental data needed for the calculation of CWR and related variables are surface precipitation, surface evaporation, and four-dimensional cloud water, water vapor, and wind fields, as well as the conversion between water vapor and hydrometeors. These data can be easily obtained from appropriate cloud-resolving models. On the other hand, the water vapor and wind fields can be obtained from reanalysis products (such as NCEP/NCAR reanalysis data) and precipitation can be observed directly, but no systematic cloud observations are available. Thus, how to obtain observational 3D time-varying cloud information (including cloud cover and cloud water content) is a key issue that needs to be solved.
In this paper, we propose two quantification methods for CWR. One is the numerical simulation quantification method based on a cloud-resolving model (CWR-NQ), and the other is the diagnostic quantification method based on observational data (CWR-DQ). We decompose the study area into grid meshes with the same horizontal resolution as the basic data, i.e., 3 km × 3 km grids for the CWR-NQ method and 1° × 1° grids for the CWR-DQ method, respectively. The CWR contributors and characteristic variables are firstly computed at each grid, and then integrated over space and time to obtain the overall results. The state variables are instantaneous, which cannot be accumulated over time; the advection variables can be integrated along the regional boundary, which cannot be simply accumulated in space; and the sink/source terms can be accumulated both temporally and spatially.

Numerical quantification
The key issue of the CWR-NQ method is the availability of a cloud-resolving model that can provide precise descriptions of cloud microphysical processes and cloud distribution. The microphysics scheme developed by Hu and Yan (1986) of the Chinese Academy of Meteorological Sciences (CAMS) includes the water content of cloud (q c ), rain (q r ), ice crystals (q i ), snow (q s ), and graupel (q g ); number concentrations of rain droplets, ice crystals, snow, and graupel particles; and cloud droplet spectral width (F c ). A total of 31 cloud microphysical processes including collision, deposition, automatic conversion between rain and cloud, nucleation, proliferation, freezing, and so on, are considered (Hu and Yan, 1986;Tao et al., 2015). Meanwhile, a quasi-implicit scheme is implemented to calculate these microphysical variables to ensure that the calculation is stable, positive definite, and conservative. The CAMS microphysics scheme was coupled into the WRF version 3.1 by Gao et al. (2011). To meet the need for operational cloud forecast and weather modification, the CAMS scheme was coupled with the WRF version 3.5 at the Weather Modification Center of the China Meteorological Administration, forming the WRF_CAMS mesoscale cloud-resolving model (Liu et al., 2016), also called the Cloud and Pre-cipitation Explicit Forecast System (CPEFS). After passing the initial examination for cloud forecast and operational service, the CPEFS has been applied for operational service in 8 regions of China since 2016, and the cloud products over China are released twice a day at CMA intranet (Sun et al., 2019). Verification of the CPEFS outputs indicates that the model can well simulate various types of cloud system in China. Consequently, the CWR-NQ method is implemented based on this model.
In the simulation, double nested domains are adopted with horizontal resolutions of 9 and 3 km in the parent and child domains (hereafter d01 and d02), respectively. There are 34 uneven levels in the vertical direction, with the model top at 50 hPa. Taking North China as an example (see Section 3), the model domain is centered at 39.9°N, 116.3°E, and the number of gird points for d01 and d02 is 254 × 182 and 406 × 352, respectively. The initial and boundary conditions of the model are from the NCEP/NCAR reanalysis data, and the lateral boundary condition is updated every 6 h. The model outputs for d01 and d02 are at 6-and 1-h intervals, respectively. The time step of integration is 60 seconds. Consecutive simulations for one month or several months can be achieved as needed.
In the numerical quantification of CWR contributors, q v is a direct model output, while q h is the sum of q c , q r , q i , q s , and q g . In addition, the model directly outputs condensation rate (r con ) and evaporation rate (r eva ). The value equals hsvc + hsvr + hsvi + hsvs + hsvg, where hsvc and hsvr are the condensation and evaporation rates of water vapor during its transformation to cloud water and rain, and hsvi, hsvs, and hsvg are sublimation and evaporation rates of water vapor transformation to ice, snow, and hail. When the value is positive, it is defined as the condensation rate; otherwise, it is defined as the evaporation rate. The unit of r con and r eva is kg kg −1 s −1 . The surface evaporation rate (e) and rain intensity (I) are calculated in the land surface process and cloud microphysics scheme of the model and are direct model outputs. Based on these parameters, the CWR contributors can be quantified according to the corresponding algorithms in Zhou et al. (2020a), and the CWR-related characteristic variables can also be obtained eventually.

Satellite data
The CloudSat/CALIPSO products used in this study include 2B-GEOPROF, 2B-CWC-RO, 2B-GEOPROF-LIDAR, 2B-CLDCLASS-LIDAR, and ECMWF-AUX (http://www.cloudsat.cira.colostate.edu). 2B-GEOPROF is cloud radar profiles, which mainly provides effective values of cloud mask to identify cloud boundary and cloud top. 2B-CWC-RO provides the profile information of cloud liquid water content and cloud ice water content . 2B-GEOPROF-LIDAR is cloud geometric profiles observed by cloud radar and lidar, which provides cloud fraction detected at each vertical grid point. 2B-CLDCLASS-LIDAR provides eight types of cloud products (Winker et al., 2007), including high cloud (Ci), altostratus (As), altocumulus (Ac), stratus (St), stratocumulus (Sc), nimbus (Ns), cumulus cloud (Cu), and deep convective cloud (Dp), with the first six types of cloud collectively called stratiform clouds (St) while the rest two convective clouds (Dc) in this paper. ECMWF-AUX is an intermediate product that contains the atmosphere state variables provided by ECMWF through the interpolation to each CloudSat cloud profiling radar (CPR) scanning point.
The CloudSat/CALIPSO satellite observation orbits pass over China twice a day around 0500 and 1700 UTC. The data used in this study cover the period from January 2007 to December 2010. Figure 1 shows the Cloud-Sat tracks over China. The satellite-borne radar has a resolution of 2.5 km along the orbit and a resolution of 1.4 km perpendicular to the trajectory. Each scan profile contains 125 grid points with a vertical resolution of 500 m. Under the condition of cloud mask ≥ 20; or cloud mask < 20 and cloud fraction ≥ 99%, this grid is considered to have clouds. If cloud mask = 0 and cloud fraction = 0, this grid is clear and cloudless (Kahn et al., 2008). For those grids that do not meet the above two conditions, to determine whether the grid has clouds or not is difficult due to possible instrument errors, hence such grids are not considered in our study. During the study period 2007-2010, there are 5994 tracks passing over China, corresponding to 7.33 × 10 8 grid points.

Cloud region detection and cloud water content statistics
In this study, we use the CloudSat/CALIPSO joint observations during 2007-2010 to first determine a climatological relationship between cloud existence and atmospheric relative humidity and temperature. Following Cai et al. (2015) for cloud profile identification, the threat score (TS) method is employed to obtain various relative humidity thresholds (RH c ) at different temperature ranges. Specifically, the CloudSat/CALIPSO joint observations are taken as true values; the TS scores for forecast or detection of cloudy and clear sky conditions (TS cloud and TS clear , respectively), corresponding to various relative humidity thresholds, are statistically calculated; and the accuracy ratio (AR), false alarm ratio (FAR), and missing alarm ratio (MAR) for the cloud forecast/detection, are also calculated. By analyzing the TS scores, the RH c standards for identifying cloud regions are derived and listed in Table 1.
Overall, the RH c for cloud identification first decreases and then increases with decreasing temperature. Within the temperature range of −5 to −25°C, RH c for cloud formation is less than 70%; when temperature is higher than 5°C and lower than −45°C, RH c rises to 80%. The TS cloud is generally higher than 0.4 for various ranges of temperature lower than 15°C, and the AR is within 62%-76%. The highest TS cloud and AR occur between −5 and −10°C, with the values reaching 0.52 and 75.48%, respectively. Meanwhile, FAR and MAR are relatively low. Despite the temperature change, TS clear is typically higher than 0.76.
By using 2B-CWC-RO and 2B-CLDCLASS-LIDAR products from 2007 to 2010, the vertical profiles of average liquid water content (LWC) and ice water content (IWC) for stratiform clouds and convective clouds with temperature over China are obtained (Fig. 2). The average LWC is typically less than 450 mg m −3 , decreasing with descending temperature. LWC of convective clouds (Dc-LWC) is greater than that of stratiform clouds (St-LWC). IWC is significantly smaller than LWC, with an average value less than 150 mg m −3 . As the temperature descends, both Dc-IWC and St-IWC increase first and then decrease, with the peak values appearing within −15 to −20°C. The vertical profiles of cloud water content for different types of cloud are obtained by adding LWC and IWC within each temperature range (Fig. 2).

Deriving the 3D cloud diagnosis products
The NCEP/NCAR atmospheric reanalysis data for temperature (T), relative humidity (RH), water vapor content (q v ), and wind (u, v) are used, with a horizontal resolution of 1° × 1° and 26 vertical levels from 1000 to 10 hPa. The data are downloaded from https://rda.ucar. edu/datasets/ds083.2/. The data are available at 0000, 0600, 1200, and 1800 UTC, with 6-h intervals.
Based on the RH c table (Table 1), the proxy climatology of LWC and IWC (Fig. 2), and NCEP/NCAR reanalysis temperature and RH fields, 3D cloud products including cloud detection (i.e., cloudy or clear diagnosis; CLD), cloud base height (CBH), cloud water content (CWC), and cloud water path (CWP; integrated cloud water content, namely M h ), are obtained over China from 2000 to 2017. The diagnostic cloud dataset has the same resolution as the reanalysis data, i.e., horizontal resolution of 1° × 1° and time interval of 6 h. In the vertical direction, CLD and CWC are located at 26 levels from 1000 to 10 hPa, while CBH and CWP are single-layer data. The 3D cloud diagnosis products (Zhou et al., 2020b) are then generated, which provide the data basis for the following diagnostic quantification of CWR. These products can be directly applied to CWR estimation in different regions of China; alternatively, localized RH c and vertical profiles/3D fields of cloud liquid/ice water content can also be derived with the above approach based on the CloudSat/CALIPSO observations, and then used to estimate the regional CWR.

Quantification of CWR and related variables
The 1° x 1° daily precipitation dataset (1DD; https:// www.ncei.noaa.gov/data/global-precipitation-climatology-project-gpcp-daily/access/) (Huffman et al., 2001;Zi et al., 2007) of Global Precipitation Climatology Project (GPCP) is used to calculate surface precipitation. Using the 3D diagnostic cloud products obtained in subsection 2.3.3, together with water vapor and wind fields extracted from NCEP/NCAR reanalysis and GPCP products, the CWR contributors are partially quantified, including the 6 state variables such as M h1 and M h2 , 6 advection items such as Q hi and Q ho , as well as P s . The source/sink terms C vh , C hv , and E s are estimated from the budget equation of atmospheric hydrometeors and total water substance. Note that we have the following formulas that can be transformed from Eq. (1), In Eq. (3), the result is treated as C vh when the value is larger than 0, otherwise it is treated as C hv . The net condensation (C vh − C hv ) calculated by the above diagnostic method is accumulative for any extended study period and area, and the bias is small. However, C vh and C hv offset each other during individual time periods and at individual grids, which results in lower diagnostic values of C vh and C hv than their actual values. Because of this, GM h and CWR are underestimated, while PE h is overestimated. As the resolution of the reanalysis data increases, C vh and C hv are calculated in grids with a higher resolution. Correspondingly, biases of GM h , CWR, and PE h estimates will decrease.

Assessment and verification of the quantification results
Using both the numerical (CWR-NQ) and diagnostic (CWR-DQ) quantification methods, we take North China in April and August 2017 as an example to demonstrate the CWR calculation results. Gridded values of CWR and related variables are quantified first, and the monthly values over the entire region are obtained then. The double-nested model domain is displayed in Fig. 3, in which the black box (d02) is the study region for both methods, which covers an area of 1.34 × 10 12 km 2 over 34°-44°N, 108°-121°E.
In order to validate the two quantification methods, rationality and consistency of the results obtained by the two methods are compared and discussed below.

Verification of the numerical quantification results
The challenging issue in the CWR-NQ method is to ensure that the atmospheric water mass remains conserved and realistic after continuous simulation over a long period of time. The results of numerical simulations are verified against cloud and precipitation observations and evaluated from the perspective of atmospheric water mass conservation. 3.1.1 Evaluation of the atmospheric water mass conservation According to Table 2  for atmospheric hydrometeors (M h1 + Q hi + C vh ) is 53.36 billion ton in April 2017, and the sum of expenditure terms (M h2 + Q ho + C hv + P s ) is 53.31 billion ton. The deviation is 0.05 billion ton, equivalent to 0.09% of GM h . For atmospheric water substance, the sum of income and expenditure terms, (M w1 + Q wi + E s ) and (M w2 + Q wo + P s ), are 402.87 and 404.15 billion ton, with a difference of −1.28 billion ton, which is equivalent to 0.13% of GM w . The deviations as percentages of GM h and GM w in August are 0.94% and 0.13%, respectively. Thus, after the CPEFS model runs continuously over a month, the water substance in the atmosphere remains approximately conserved.

Evaluation of simulated cloud and precipitation
To verify whether the spatial distributions of numerically quantified MM h are close to their observations, the simulated MM h is compared with CWP retrieved from the CERES satellite data in April 2017 (Fig. 4). Note that MM h bears the same physical meaning as the CWP retrieved from satellite products; for easy comparison, MM h in original unit of kg is converted to column hydrometeors content in water depth over a unit area (g m −2 or μm). The simulated MM h has a good consistency with the CWP retrievals from the CERES, and both show a gradual decrease from south to north, with a maximum value of 350 g m −2 . The simulated large-value center deviates slightly from that retrieved from satellite observations. In addition, the simulated MM v is compared with that derived from the NCEP/NCAR reanaly-sis (figures omitted). The magnitude, spatial distribution, and largevalue centers of the simulated MM v are consistent with those in the reanalysis water vapor product.   Fig. 4. Comparison of (a) simulated MM h with (b) CWP retrievals from CERES (unit: g m −2 ) over North China in April 2017.
DECEMBER 2020 dataset is obtained from the CMA intranet. The simulated and observed precipitation agree well in their monthly spatial distributions (Fig. 5) and daily variations (Fig. 6), and the values of simulated daily precipitation are highly consistent with the observations. The CPEFS model can well reproduce the multiple precipitation    events that occurred in April and August 2017. The observed rainbelt is stronger in the south than in the north and in the east than in the west, which is reasonably simulated by the model. The maximum accumulated precipitation is 100 and 270 mm for April and August 2017, respectively, and the values are consistent between model simulations and CMPA-Hourly observations. Apparently, the atmospheric water mass remains conservative in the numerically quantified results and are consistent with observations, indicating that the CWR-NQ method and results are reasonable.

Evaluation of the diagnostic 3D cloud fields
Satellite remote sensing is an effective approach for large-scale cloud observations. The CERES cloud products have a high spatial resolution. Therefore, horizontal distributions of the diagnosed clouds are compared with the CERES observations. Monthly mean total cloud cover (TCC) and CWP for the period 2001-2009 from CERES_SYN1deg_Ed4.1 (Young et al., 1997;Wong et al., 2000;Chambers et al., 2002;Minnis et al., 2020), which is downloaded at https://ceres-tool.larc.nasa. gov/ord-tool/jsp/SYN1degEd41Selection.jsp, are used to verify the diagnostic TCC and MM h over China for the same period. The diagnostic TCC is derived from the cloud identification product. For the 1° × 1° column in any time, when cloud is identified in no less than 1 grid point in the vertical direction, this column is considered to have clouds at this time. In a certain period, the ratio of cloud occurrence times to the total number of times involved in the statistics is the TCC of the column.
Comparison of spatial distributions of TCC between the diagnostic results and the CERES is shown in Fig. 7. In general, the diagnosed 2001-2009 averaged TCC (Fig.  7a) agree well with that retrieved from the CERES (Fig.  7b) in their spatial distribution. The high values of TCC are mainly located over the Sichuan basin, the Tianshan Mountain, and the Himalayas, while the desert areas in Xinjiang Region, the southern Qinghai-Tibetan Plateau, the Inner Mongolia, and North China are low TCC areas. In the northeastern part of Heilongjiang Province, however, the diagnosed TCC is slightly higher than that from the CERES. Sun-Mack et al. (2007) compared global TCC observed by different satellites in April 2007, and found that the global average TCC from the CERES is about 0.14 lower than that observed by CALIPSO. We thus infer that it is reasonable for the diagnosed TCC being slightly higher than that from the CERES.
Comparison between the diagnosed MM h and observed CWP is displayed in Fig. 8. Overall, the annual average distribution of the diagnosed MM h is consistent with the CWP products retrieved from the CERES satellite observations. MM h gradually decreases from south to north, presenting a distinct zonal distribution pattern. High values of MM h are found over the Sichuan basin and southeast of Zhejiang, Jiangxi, and Hunan provinces. Low values of MM h are largely located in western Xinjiang Region, Qinghai Province, Tibetan Region, and some other places of Northwest China (except for the west Tianshan Mountain in Xinjiang). The above features are consistent with the distributions of CWP in China shown in the study of Li X. Y. et al. (2008) based on the ISCCP data. The diagnosed MM h is within the range of 200-800 g m −2 , and the CERES CWP is within 100-450 g m −2 . The diagnostic value is higher than that retrieved from the CERES, which may be attributed to the fact that the diagnosed TCC is higher than the retrieved TCC. On the other hand, the CERES CWP may be underestimated. Comparison of ground-based radar and lidar observations and satellite data on overcast stratus cases between March 2000 and December 2004 showed that, the differences between the satellite and  surface values of LWP are 17.2 ± 32.2% (Dong et al., 2008). Another comparison of marine boundary layer cloud LWP from CERES Edition 4 and ground-based 3mm cloud radar at the Azores indicated that, satellite-retrieved LWP is smaller than that from the cloud radar retrieval (Xi et al., 2014). In addition, Tian et al. (2018) found that CERES IWP is lower than the GOES (Geostationary Operational Environmental Satellite) retrieval by 12%. Thus, the CERES CWP may be in general underestimated.
In addition, we also compared the diagnosed cloud base height (CBH), a parameter indicating the cloud vertical structure, with that retrieved from a laser ceilometer (Viasala CL51) installed at Southern Beijing Observatory (39.80°N, 116.47°E) for the entire year of 2015 with 4-time daily data, and they are generally accordant. The above comparisons demonstrate that the horizontal distribution and vertical structure of the cloud fields obtained by the diagnostic method are generally consistent with the corresponding observations. The 3D cloud diagnosis scheme proposed in this paper is reasonable, and can provide 4-time daily data of 3D spatial distributions of cloud cover and cloud water content, especially the vertical structure of large-scale clouds that cannot be obtained from satellite remote sensing and ground-based observations.

CWR contributors
In the CWR-DQ method, the two terms C vh − C hv and E s , which are key to the rationality of the diagnostic results, are derived from the balance equations of water substance and hydrometeors in the atmosphere [see Eqs. (3) and (4) in subsection 2.3.4]. Their direct observations cannot be obtained, and thus they are compared with numerical simulations. Table 3 presents the results of sink/source terms of the two methods with their difference (ΔV in percentage) calculated by where V DQ and V NQ refers to the diagnostic value and numerical value, respectively. According to Table 3, the diagnosed E s is higher than the simulated one by about 29.37% (April) and 9.4% (August), respectively; however, the simulated P s is also smaller than the observed by 24.09% and 13.34% in the two months, respectively. Therefore, the diagnosed E s is reasonable. The diagnosed C vh is slightly lower than the simulated with differences of 4.94% (April) and 16.12% (August); the diagnosed C hv is significantly lower than the simulated by 79.56% in April, and the difference can reach up to 1098.14% in August. This may be because in the diagnostic method, C vh − C hv is firstly estimated based on the hydrometeors equation and the result is then attributed to C vh or C hv depending on its signs. For a small area and during a short period, C vh competes with C hv ; as the area and time period increase, they will offset each other. Consequently, the difference of C vh − C hv between the two methods could be small, while diagnostic C vh and C hv values are smaller than the simulation values. Table 4 lists the fitting equations and the correlation coefficient (R) of quantified daily E s , C vh − C hv , and C vh obtained by the two methods. The diagnostic daily value of the three parameters is very close to their simulated values with R in the range of 0.69-0.97. This indicates that the residual terms of the water budget/balance equa-  Table 5 lists values of the CWR-related characteristic variables obtained from the two methods. Since the diagnostic C vh and C hv are lower than their numerical values, the diagnostic CWR is also smaller than its numerical value. However, the difference in GM h between the two methods are smaller than that in CWR. In addition, the higher diagnostic P s results in larger diagnostic PE h than its numerical result. Values of GM v and PE v calculated from the CWR-DQ method are also slightly higher than those from the CWR-NQ method with differences overall less than 20%. The results of CE v , PE v , and RT v are relatively close, with a smaller difference than that of PE h and RT h .
The above comparison and analysis indicate that the results obtained by the two quantification methods are consistent, and they can mutually support and confirm each other. However, the underestimation of C hv leads to the underestimation of CWR and overestimation of PE h , which is the main deficiency of the CWR-DQ method.

Discussion
Two quantification methods, the CWR-NQ method and the CWR-DQ method, have been developed and applied for CWR study over North China in April and August 2017. The results are compared and evaluated to demonstrate the feasibility of the two methods in quantifying the CWR.
To obtain objective results, the CWR-NQ method requires that the mass of total water substance in the atmosphere remain conserved even after long-term continuous model integration, and the simulated cloud fields and precipitation be as realistic to the observations as possible. In this paper, the simulation of CWR in North China in April and August 2017 shows that after one month of continuous integration, the mass of atmospheric water substance is basically conserved, and the temporal and spatial distributions of water vapor, cloud water, and precipitation are consistent with the observations. These results indicate that the CWR-NQ method is reasonable. The obtained CWR results with high spatiotemporal resolution are the main advantage of this method. In the future, with increasingly advantageous computational resources and constantly improved cloud resolving models, the CWR-NQ method is expected to play an important role in the mid-and long-term refined forecast of regional CWR and provide guidance in the design of operational schemes for artificial precipitation enhancement and CWR exploitation. The difficulty of the CWR-DQ method lies in the fact that C vh − C hv and E s cannot be measured directly. They are, respectively, the residual term of the atmospheric hydrometeors balance equation and atmospheric water substance balance equation. Comparison of quantification results obtained by diagnostic and numerical methods for April and August 2017 in North China shows that the discrepancies of C vh − C hv and E s between the two methods are, respectively, 22.94% and 29.37% in April, and 14.04% and 9.4% in August. This indicates that the two independent methods are in general consistent, and it is feasible to implement the diagnosis method to approximate C vh − C hv and E s from the atmospheric water substance balance equations. Correlation analysis of daily quantification results shows that the diagnostic quantities of CWR contributors and characteristic variables are correlated as well and consistent with the numerical results.
It is worth noting that C vh and C hv are underestimated by the CWR-DQ method, which in turn leads to the underestimation of CWR and overestimation of PE h , especially in summer when precipitation is strong. However, the diagnostic C vh − C hv is well correlated and consistent with its numerical result, suggesting that it may be feasible to use numerical quantification results to correct biases of diagnostic C vh and C hv . By improving the spatiotemporal accuracy of observations and increasing vertical resolution in the calculation, the biases of diagnostic C hv are expected to reduce to achieve better accuracy of diagnosed CWR, PE h , and so on in the future.
Based on the CWR-DQ method, combined with longterm satellite cloud observations, atmospheric reanalysis products and precipitation observations, we can realize climatic quantification of atmospheric hydrometeors content and CWR in China and around the globe. Furthermore, investigation of spatiotemporal distribution characteristics of CWR, and associated features and trends in various regions will provide a scientific basis for CWR exploitation plan and implementation.

Conclusions
The primary purpose of this study is to establish two quantification methods for cloud water resources based on numerical simulation (CWR-NQ) and observational diagnosis (CWR-DQ). By using the two methods, CWR quantification in North China are then achieved in April and August 2017, and the quantification results obtained from the two methods are compared to the observation and between each other. Major conclusions are as follows.
(1) Based on the CloudSat/CALIPSO joint cloud observations, a 3D cloud diagnosis method is optimized, which solves the key difficult issue of the CWR-DQ method. The diagnosed cloud products are compared with observations of satellites and laser ceilometers, revealing that the horizontal and vertical structures of the diagnosed clouds are in good agreement with observations.
(2) For the CWR-DQ method, the diagnostic CWR contributors and related characteristic variables are consistent and well correlated with those obtained from the CWR-NQ method. This suggests that, although the diagnosed C vh and C hv are smaller than their numerical values, the CWR-DQ method is still reasonable. In the future, by improving the spatiotemporal resolution and accuracy of observations and calculations, the biases of diagnostic C vh and C hv will be reduced to further improve the accuracy of the diagnostic CWR and related variable.
(3) For the CWR-NQ method, after one month of continuous integration of the WRF_CAMS mesoscale cloudresolving model (i.e., CPEFS), the mass of water substance in the atmosphere is basically conserved, and the temporal and spatial distributions of water vapor, cloud water path, and precipitation are consistent with observations, indicating that the CWR-NQ method is reasonable.
In general, the CWR-NQ method and the CWR-DQ method proposed in this study are basically acceptable. However, in practical applications, more in-depth investigations are needed to optimize variables such as E s , C vh , C hv , and so on. The two CWR quantification methods could have broad application prospects in the future. The CWR-NQ method can be used for short-and mediumterm CWR forecasting and climatic CWR prediction; while the CWR-DQ method is more suitable for longterm worldwide CWR research. Analysis of the spatial and temporal scales of variables associated with CWR, as well as their distribution and variation characteristics can provide valuable information for the exploitation of CWR in various regions.