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

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 water 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 atmosphe-ric 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 large-scale 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 observations. 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 one-second 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 verti-cal 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, Zhou et al. (2010) 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 analyzed 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 numeri-cal 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 quan-tify 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.

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.

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 Precipitation 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⁻¹ s⁻¹. 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.

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 (Stephens et al., 2008). 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 CloudSat 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⁸ grid points.

Figure 1.  Schematic diagram showing tracks of the CloudSat satellite passing over China during 2007–2010. The southwest–northeast (northwest–southeast) slanted lines denote the tracks around 0500 (1700) UTC. The distance between any two neighboring tracks is about 2.5° and the tracks are repeated every 16 days on average.

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⁻³, 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⁻³. 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).

Figure 2.  Vertical profiles of liquid water content (LWC; mg m⁻³) and ice water content (IWC; mg m⁻³) for stratiform clouds (St) and convective clouds (Dc) averaged over China. St-IWC and St-LWC denote IWC and LWC for stratifum clouds while Dc-IWC and Dc-LWC denote IWC and LWC for convective clouds. The results are based on the statistics of cloud products from CloudSat/CALIPSO passing over China during 2007–2010.

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.

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.

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¹² km² over 34°–44°N, 108°–121°E.

Figure 3.  Model domain for CWR assessment over North China in April and August 2017 (d02 area shown by the black box is the region selected for analysis in this study).

In order to validate the two quantification methods, rationality and consistency of the results obtained by the two methods are compared and discussed below.

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.

According to Table 2 below, the sum of income terms 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.

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⁻² 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⁻². 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 large-value centers of the simulated MM_v are consistent with those in the reanalysis water vapor product.

Figure 4.  Comparison of (a) simulated MM_h with (b) CWP retrievals from CERES (unit: g m⁻²) over North China in April 2017.

Figures 5 and 6 present the comparison of simulated precipitation with CMPA-Hourly v1.0 [China Hourly Merged Precipitation Analysis combining observations from automatic weather stations with CMORPH on 0.1° × 0.1° grid; Shen et al., 2014; CMORPH denotes the US Climate Prediction Center (CPC) morphing technique satellite precipitation estimates]. The CMPA-Hourly v1.0 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.

Figure 5.  Distributions of monthly accumulated precipitation from (a, c) CPEFS model simulations and (b, d) CMPA-Hourly v1.0 observations in (a, b) April and (c, d) August 2017.

Figure 6.  Evolution of daily precipitation from the CPEFS model simulation (red line) and CMPA-Hourly v1.0 observation (black line) in (a) April and (b) August 2017.

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.

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.

Figure 7.  Multi-year (2001–2009) averaged distributions of (a) diagnosed TCC and (b) CERES observed TCC (unit: %).

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⁻², and the CERES CWP is within 100–450 g m⁻². 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 3-mm 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.

Figure 8.  Multi-year (2001–2009) averaged distributions of (a) diagnosed MM_h and (b) CERES retrieved CWP (unit: g m⁻²).

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.

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 diagnos-tic 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 equations in the atmosphere are relatively reasonable in the CWR-DQ method. It is thus feasible to approximate E_s as the residual term of the atmospheric water substance balance equation, and approximate C_vh − C_hv to the residual of the atmospheric hydrometeors balance equation.

Comparisons of other CWR contributors are also listed in Table 4. Since the state variables change with time, MM_h and MM_v are compared instead. In general, daily values of MM_h, Q_hi, Q_ho, P_s, MM_v, Q_vi, and Q_vo obtained by the two methods are well correlated. Except for Q_vo in April, R for each individual contributor in both months exceeds 0.78. Among the contributors, P_s and MM_v have the highest correlation with R exceeding 0.95. The fitting equations and curves of daily values (figures omitted) indicate that diagnostic and numerical quantification results of MM_v are close, while the other 6 variables have higher diagnostic value than their numerical simulated counterparts.

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 numeri-cal 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.

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 atmosphe-ric 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 rea-sonable. 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 diagnos-tic 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 long-term 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.

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 cloud-resolvingmodel (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 medium-term CWR forecasting and climatic CWR prediction; while the CWR-DQ method is more suitable for long-term 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.