The total mass and spatial–temporal variability of aerial cryosphere over the Tibetan Plateau from 2003 to 2020

Changes in snow, ice, and ecological system over the Tibetan Plateau (TP) are extremely sensitive to local precipitation and radiation budget, which are largely modulated by the atmospheric ice. However, how much ice is there in the atmosphere over the TP and how it is distributed are still unclear. The total mass, spatial distributions, and long-term trends of atmospheric ice over the TP were evaluated by using four sets of satellite retrieval data (Aqua, Terra, the Suomi National Polar-orbiting Partnership (Suomi NPP), and NOAA-20) and ERA5 reanalysis data from 2003 to 2020. Based on the estimations using multiple satellite datasets, we concluded that the total mass of atmospheric ice could be up to 0.26 ± 0.03 Gt over the TP from 2013 to 2020. The spatial distributions of atmospheric ice derived from various datasets were highly consistent. In general, the southwest and northeast areas of the TP were the low-concentration areas (0.05 kg/m2 in average), while the southeast area was the high-concentration area (0.09 kg/m2 in average), and this spatial pattern was most evident in summer. The high values around (0.15 kg/m2) were centered over Linzhi and its surrounding areas. The plentiful water vapor transported by southwest summer monsoon and steep topography jointly led to rapid growth of atmospheric ice in Southeast Tibet, which was the dominant reason for the higher ice concentration in this area.


Introduction
As a sensitive area of climate change, the air temperature of the TP had increased much faster than that of the other areas at the same latitude (Liu and Chen, 2000), and the mass and energy balance of ice and snow in this area are greatly influenced by radiation and precipitation associated with the aerial cryosphere. The aerial cryosphere is composed of all ice bodies in the atmosphere, including those in clouds (ice clouds and ice crystals) and those under clouds (snowflake and hail) (Qin et al., 2018). It is an element with the largest coverage in the cryosphere, which has an important impact on the radiation budget and water cycle in the climate system (Dou et al., 2020). Previous studies have indicated that the atmospheric ice dominated the distribution and lifetime of ice clouds and mixedphase clouds over the TP (Long, 2016;Li et al., 2018;Zhao et al., 2021), thus profoundly affecting the surface radiation budget. Therefore, it is of great significance to investigate the spatial-temporal variability of the atmospheric ice over the TP.
With the rapid progress of satellite remote sensing technology and the improvement of data assimilation, the remote sensing data and reanalysis data were increasingly used to explore the pattern of different phases of water in the atmosphere in the past few decades (Liljegren et al., 2001;Horváth and Davies, 2007;O'Dell et al., 2008;Lebsock and Su, 2014;Khanal et al., 2020;Qi et al., 2022). Some studies suggested that the cloud water path (CWP) increased in most regions of the TP from the 1980s to the beginning of this century (Li et al., 2008;Yue et al., 2016). The annual average of the ice water path (IWP) had been relatively stable from 2000 to 2015 without any obvious trend (Li et al., 2018). In contrast, on the seasonal scale, the IWP in the premonsoon season (March to May) and winter showed a significant upward trend in the western TP, while in monsoon season (June to September), a significant downward trend in the western TP was observed (Zhao et al., 2021).
Earlier studies estimated the total mass (Dou et al., 2020;Xu et al., 2022) and investigated the spatial-temporal distribution of atmospheric ice on a global scale using multiple sets of satellite remote sensing data and high-quality reanalysis data (Duncan and Eriksson, 2018;Dou et al., 2020;Xu et al., 2022). However, the total mass and spatial-temporal distribution of atmospheric ice over the TP are poorly understood, and there is also a lack of uncertainty analysis based on the multi-source data. Here, we used four sets of satellite retrieval products, namely, Aqua, Terra, the Suomi NPP, and NOAA-20 in CERES_SSF_1deg-Month L3 data and ERA5 reanalysis to carry out a quantitative study of atmospheric ice over the TP and reveal the main controlling factors of its spatial pattern. Note that acronyms and abbreviations used in this study are given in Table 1 The boundary data on the TP in this study were downloaded from the National Tibetan Plateau Scientific Data Center (TPDC). The TP is about 3,360 km from east to west and 1,560 km from south to north, ranging from 25°59′30″N to   (Zhang et al., 2002;Zhang et al., 2021a;Zhang et al., 2021b).

Topography data
The topography data we used were produced by The General Bathymetric Chart of the Oceans (GEBCO), which showed a global terrain model for the ocean and land, providing elevation data in meters, on a 15 arc-second interval grid (Group, 2019).

Remote sensing retrieval data
The remote sensing products were retrieved by the Terra, Aqua, Suomi NPP, and NOAA-20 satellites from the CERES_SSF_1deg-Month L3 data. Clouds and the Earth's Radiant Energy System (CERES) was a NASA project that looked into the role of cloud and radiation feedback in the climate system (Wielicki et al., 1996). When coupled with the global climate model, these data contributed to a better understanding of how clouds and radiation interacted with atmospheric thermodynamics. There were currently five satellites carrying CERES instruments, and four satellites being employed in this study were sun-synchronous. The Terra satellite crosses the equator from south to north, and because it does so, at about 10:30 a.m., it is also known as the "Morning Satellite." It has two CERES instruments, Flight Model 1 (FM1) and Flight Model 2 (FM2), as well as a Moderate Resolution Imaging Spectroradiometer (MODIS). A MODIS has 36 bands with a wavelength range of 0.412-14.45 μm and a spatial resolution of 250-1,000 m, allowing it to cover the entire globe in 2 days (Barnes et al., 2003). The Aqua satellite crosses the equator from north to south at about 1:30 p.m., and so it is known as the "Afternoon Satellite." It has two CERES instruments, FM3 and FM4, and a MODIS. The "Afternoon Satellite" Suomi NPP is also equipped with the CERES instrument FM5 and the Visible Infrared Imaging Radiometer Suite (VIIRS). The VIIRS contains 22 imaging and radiation bands with a wavelength range from 0.41 to 12.5 μm and a spatial resolution of 370-750 m (Cao et al., 2014). NOAA-20, the successor of the Suomi NPP, which is also an "Afternoon Satellite," is equipped with the FM6 CERES instrument and VIIRS (Cao et al., 2018). The four satellites mentioned previously were equipped with advanced spectral imagers so as to provide 24-h day-time and night-time IWP monthly average data CERES_SSF_1deg-Month L3 data, with a spatial resolution of 1°×1°. These products had the advantages of long time series, wide application range, and high space coverage. We chose 24-h day-time and night-time IWP data because it better represented the whole IWP in the research area even though the retrieval algorithms used during the day-time and night-time were different. Meanwhile, because the time series of the Suomi NPP and NOAA-20 satellite products were short, the payloads carried by them were the same, and their orbits were similar, so they were merged into one new product with the time series from 2013 to 2020.

Reanalysis data
The reanalysis data used here were ERA5, which was the fifth-generation reanalysis data product made by the European Centre for Medium-Range Weather Forecasts (ECMWF). A new data assimilation system, known as the ECMWF Integrated Forecasting System "Cy47r3" 4D-var, was applied for ERA5 by ECMWF. It was based on numerical calculation of the model and assimilation of a large amount of observation data. Due to these advantages, ERA5 contain detailed records of the global atmosphere, surface, and ocean waves since 1950. It had a significant improvement in the spatial resolution and numerical accuracy of variables when compared to the previous reanalysis data product ERA-Interim (Hersbach et al., 2020). We selected the monthly averaged ERA5 on single levels from 1979 to the present as the research data and employed the combination of two variables, total column cloud ice water and total column snow water, as the IWP to characterize the total ice content.

Methodology
IWP (g/m 2 ) is a primary indicator for calculating the atmospheric ice concentration. It was defined as the vertical integral of the ice water content (IWC) (g/m 3 ) in the atmospheric column (Holl et al., 2014).
The CERES_SSF_1deg-Month L3 data assumed that the atmosphere was plane parallel, and the whole cloud column was composed entirely of ice. Hence, the IWP value was the product of the effective size and optical thickness of ice crystals. The specific algorithm was as follows: where ρ is the density of ice (−0.9 g cm −3 ), R e is the effective radius, τ is the optical depth, and Q is the extinction efficiency, which depends on R e and ranges from 2.01 to 2.11 (Tian et al., 2018).

Total mass
The three satellite retrieval products (Aqua, Terra, and the Suomi NPP-NOAA-20) during 2013 and 2020 were used to estimate the total mass of atmospheric ice over the TP by considering the area weight because the Suomi NPP and Frontiers in Earth Science frontiersin.org NOAA-20 products were only produced since 2013 ( Table 2). The total mass of atmospheric ice calculated based on Aqua and Suomi NPP-NOAA-20 was very close (~0.28 Gt) during this period, which is larger than that from Terra (0.23 Gt). The reason for this condition might be the difference in their orbits. The Aqua and Suomi-NPP were both A-train satellites. Combined with the estimation results of all satellite datasets, we concluded that the multi-year average total mass of atmospheric ice was 0.26 ± 0.03 Gt (~2.6 ×10 8 tons) over the TP, which accounted for approximately 5‰ of the global aerial cryosphere. The total mass of atmospheric ice had a significant seasonal variation, with a maximum of 0.32 ± 0.05 Gt in summer and a minimum of 0.22 ± 0.01 Gt in winter. In general, the result of ERA5 (total mass 0.24 Gt) was comparable with that of Terra and smaller than those of Aqua and Suomi NPP-NOAA-20. This is also applied on a global scale, and the mass of atmospheric ice given by the reanalysis data was also on the low side. The main reason was that the current atmospheric models did not perform well enough in describing the ice formation mechanism, ignoring the secondary ice formation mechanism and underestimating the number concentration of ice condensation nuclei (ECMWF, 2021).

Spatial distribution and seasonal variability
The ice concentration ranged from 0.05 to 0.15 kg/m 2 over the TP, with an average of 0.08 kg/m 2 . As shown in Figure 1, the spatial variability of atmospheric ice over the TP was obvious. Three sets of satellite retrieval products consistently revealed the spatial pattern of high ice concentration in the eastern TP and low ice concentration in the western TP ( Figures 1A,C,D). ERA5 generally underestimated the spatial heterogeneity of atmospheric ice content, and there are only a few high-value areas in the southeast edge of the TP ( Figure 1B). All the four datasets observed the high-value center of the ice concentration around Linzhi City in southeast Tibet, where the value was up to 0.15 kg/m 2 . Additionally, the IWP over the southern slope of the TP was relatively higher than those of other areas. Overall, the spatial distributions derived from various satellite datasets were consistent, particularly between Aqua ( Figure 1A) and Suomi NPP-NOAA-20 ( Figure 1D).
In order to clarify the reasons for the spatial difference of atmospheric ice over the TP, the vertical profile of wind field, relative humidity at 95°E, 30°N, and the IWP of corresponding grids were analyzed from 2003 to 2020 combined with the topography (Figure 2). A large amount of water vapor was   Frontiers in Earth Science frontiersin.org 05 transported to the southeast of the TP (around 27°N) through the Bay of Bengal channel, with the relative humidity exceeding to 65% at 850 hPa in this area, which provided a sufficient water vapor condition for the formation of atmospheric ice. In addition, strong updraft carried abundant water vapor to high altitudes, producing lots of atmospheric ice through condensing. The previous two

Frontiers in Earth Science
frontiersin.org 06 factors are the main reasons for the high-value center of atmospheric ice in southeast Tibet. The terrain blocking prevented water vapor from northward penetration, resulting in the low ice concentration in northern TP.
Affected by southwest monsoon, the content of atmospheric ice in summer and autumn was generally higher than that in the other two seasons, especially for the eastern TP (Figure 3). This led to higher annual mean ice concentrations in the eastern TP. The seasonal distributions of the ice content derived from various satellite retrieval products were also consistent, only with slight differences in the central and western TP inland areas (Figure 3, Supplementary Figure S1, S2). For example, the IWP of the Aqua satellite in winter was slightly higher than those of Terra and Suomi NPP-NOAA-20 in this region during the same period ( Figure 3D; Supplementary Figures S1D, S2D).
ERA5 reanalysis could basically capture the seasonal variations of the atmospheric ice content derived from satellite retrieval products, although there were still some discrepancies in the amount of ice. In summer, the ice content tended to decrease from the east to west but that from the ERA5 reanalysis tended to decrease from southeast to northwest, and the IWP was overestimated by −0.02 kg/m2 in the southeast Tibet ( Figure 2B, Figure 3B, Supplementary Figures S1B, S2B). However, in winter, ERA5 underestimated the IWP value in most areas of the TP, except for the northwest TP ( Figure 2D, Figure 3D, Supplementary Figures S1D, S2D).

Annual variability and long-term trend
As shown in Figure 5, the variability of the atmospheric ice content over the TP within a year presented a distinct single peak structure, with an upward tendency in the first half of the year, peaked in summer, and then began to decline. Among them, the growth of atmospheric ice was more "twists and turns" than its reduction process (a tiny maximum would arise around March-May), and this pattern was noticeable across all satellite products. For the annual maximum value of the ice content, Aqua was close to Suomi

Frontiers in Earth Science
frontiersin.org NPP-NOAA-20 and ERA5, and higher than Terra. As for the minimum value, the three sets of satellite data were relatively close and significantly higher than ERA5 ( Figure 5), indicating that the underestimation of the atmospheric ice content over the TP by ERA5 was mainly caused by the simulation deviation of ice growth in winter due to the underestimation of water vapor over the TP during this season. The atmospheric ice content showed no obvious trend over the TP before 2017 whether for satellite products or reanalysis data, while the annual minimum ice content has increased in recent years ( Figure 5). Similarly, the annual average ice content also increased from 2017, except for the result of Aqua ( Figure 5). Dou et al. (2020)pointed out that the total mass of the global atmospheric ice showed an increasing trend with climate warming in the past decades. However, this study indicated that the mass of atmospheric ice in the TP had not completely followed this law but showed certain regional characteristics, although it seemed to have increased over the last few years.
From the time series of the atmospheric ice content in different seasons, there were great differences between various datasets, which were not only reflected in the IWP value but were also reflected in the change trend. In general, Aqua was close to Suomi NPP-NOAA-20 and higher than the other two datasets. Terra was close to ERA5 in spring and summer, while it is significantly higher than ERA5 in autumn and winter ( Figure 6). The interannual variability of various datasets was relatively consistent in autumn and winter, but there were obvious discrepancies in the other seasons. Additionally, all sets of data showed a consistent increasing trend in winter in recent years, although there were some differences in the ice mass. In spring and autumn, Suomi NPP-NOAA-20 and ERA5 showed a consistent increasing trend, while the other two datasets showed a weak decreasing trend ( Figure 6).
In view of the increasing trend of the atmospheric ice content in winter in the past few years, we further analyzed the spatial manifestation and the possible causes of this trend. By comparing the spatial distribution of the atmospheric ice content in 2017 (low-value year) and 2019 (high-value year), it can be seen that the increase of atmospheric ice was reflected in most areas except the northeast TP ( Figure 7). Specifically, the increase in the ice content based on Suomi NPP-NOAA-20 is the most significant, with an increase of 0.08 kg/m 2 in the central and southwest regions of the TP. It can be seen from the moisture flux in winter that the water vapor convergence was clearly intensified over the TP since 2017 (Supplementary Figure. S3), which provided sufficient water vapor for the growth of ice crystals, resulting in the overall increase of the atmospheric ice content in winter during recent years. This study applied four sets (Suomi NPP was merged with NOAA-20, and there are actually three datasets) of satellite retrieval products and ERA5 reanalysis to estimate the total mass, spatial-temporal distribution, and long-term trend of atmospheric ice over the TP, analyzed the causes of spatial differences, and discussed the uncertainties of the results. This study indicates for the first time that the total mass of atmospheric ice over the TP was to be 0.26 ± 0.03 Gt, with the maximum in summer (0.32 ± 0.05 Gt) and the lowest in winter (0.22 ± 0.01 Gt). This result was obtained from three sets of satellite data, and there were some differences among them. For two datasets using the MODIS instruments, according to the quality summary, the Aqua satellite did not have calibration anomalies, but Terra's band 1 (central wavelength of 0.65 μm) had two calibration anomalies in 2003 and 2009 and one in band 5 (central wavelength of 1.2 μm) in 2002 (CERES, 2021). Coupled with the difference of their orbits, these factors will result in the difference of retrieval values of the IWP. In comparison with ground-based, ship-based, and airborne observations and other high-quality satellite data (e.g., CALIPSO data), Aqua-MODIS satellite data showed the highest performance, and Suomi NPP-VIIRS data were overall comparable to Aqua-MODIS and would gradually replace Aqua and Terra MODIS retrieval data completely in the future (Minnis et al., 2016;CERES, 2021). In conclusion, although there are some differences among various datasets, the total mass of atmospheric ice over the TP estimated based on them is close and thus, has high reliability.
Under the influence of monsoon, the water vapor on the TP is more abundant in summer and autumn, especially in southeast Tibet, resulting in a higher atmospheric ice content in these two seasons and a spatial pattern of low concentration in the west and high concentration in the east. The high-value center (0.15 ± 0.035 kg/m 2 ) was located in Linzhi City and its surrounding area, which could be well captured by ERA5 reanalysis. In terms of spatial distribution, the results of different satellite data and ERA5 reanalysis are consistent, although the latter slightly underestimates the spatial heterogeneity. For the ice content, the value given by Aqua is comparable with that of Suomi NPP-NOAA-20 and is larger than those of the other two datasets in the eastern TP. ERA5 is close to Terra in the description of the atmospheric ice content.
From the time series of atmospheric ice over the TP, it can be seen that the annual variation is unimodal. Since 2017, the annual minimum increased, and the annual amplitude declined. This might affect the frequency and distribution of ice clouds over the TP, which needed further study. There was no evident long-term trend before 2017, while after that, apart from Aqua, the other datasets showed that the mass of atmospheric ice has increased over the TP in recent years. It is to be noted that, in winter, all four datasets show a consistent increasing trend since 2017, and this trend is reflected in most areas except the small area of the northeast TP. Since 2017, there is a significant increase in the moisture convergence over the TP, which promotes the growth of atmospheric ice, leading to an increasing trend of the ice content during this period. The largest uncertainty of this study comes from the discrepancies among different satellite products. Next, it is necessary to verify the reliability of the present satellite products by using the in situ measurement, ground radar observations, or other high-quality satellite data. It is also necessary to use the new generation of weather radar, terahertz sensors, and lidar to retrieve more advanced atmospheric ice products.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author. The boundary data of the TP was downloaded from the National Tibetan Plateau Scientific Data Center (TPDC) (http://data.tpdc. ac.cn/zh-hans/), and the topographic data was downloaded from the website of GEBCO Project (https://download.gebco.net/). Remote sensing data was downloaded from the NASA Langley Research Center CERES ordering tool at (https:// ceres.larc.nasa.gov/data/). ERA5 reanalysis data was downloaded from Copernicus Climate Change Service (C3S) Climate Date Store (https://cds.climate.copernicus.eu/cdsapp#!/ search).

Author contributions
TD and CX jointly conceived the study; YY, GX, and TD performed the analyses; all of the authors discussed the results, contributed to interpretations, and wrote the manuscript.

Funding
This study was funded by the National Key Research and Development Program of China (2018YFC1406104) and the National Nature Science Foundation of China (NSFC 41401079).