Evaluation of Performance of Polar WRF Model in Simulating Precipitation over Qinghai-Tibet Plateau

The Qinghai-Tibet Plateau (TP), an area covering 2.5 million km² and an average elevation exceeding 4000 m, is the largest and the highest plateau on the earth (Bibi et al. ^[1]; Chen et al. ^[2]; Gao et al. ^[3]; Kang et al. ^[4]; Li et al. ^[5]; Tong et al. ^[6]). To some extent, the plateau is like a huge high mountain system, which is composed of mountains, plateaus, and rivers and interspersed with extensive broad valleys, basins, and lakes. Besides the extremely complex and various terrain, there are also a large number of glaciers covered by vast amounts of permafrost and snow, which rank third except for the Antarctic and Arctic (Kang et al. ^[4]; Yang et al. ^[7]; Yao et al. ^[8]; You et al. ^[9]). Hence, the Qinghai-Tibet Plateau was also known as the "Third Pole" (Bibi et al. ^[1]; Chen et al. ^[2]; Gao et al. ^[3]; Kang et al. ^[4]; Li et al. ^[5]; Tong et al. ^[6]; You et al. ^[9]; Li et al. ^[10]; Sun et al. ^[11]; Yao et al. ^[12]; Yao et al. ^[13]). TP plays a crucial role in the earth's climate system and acts on East Asia, even global weather and climate. Many excellent studies have confirmed that the forcing of the TP's significant dynamic and thermal effects change the difference between land and sea, affect atmospheric circulation, modulate the East Asian Monsoon, South Asian Monsoon, and the mid-latitude westerly and their synergistic effects (Kang et al. ^[4]; You et al. ^[9]; Chen et al. ^[14]; Chen et al. ^[15]; Chen et al. ^[16]; Duan et al. ^[17]; Duan and Wu ^[18]; Jiang et al. ^[19]; Kuang et al. ^[20]; Lai et al. ^[21]; Liu et al. ^[22]; Liu et al. ^[23]; Wang et al. ^[24]; Wang et al. ^[25]; Wang et al. ^[26]; Wu et al. ^[27]; Wu et al. ^[28]; Ye ^[29]; Zhang et al. ^[30]; Zhou et al. ^[31]). Moreover, the significant processes of ice/snow and water cycle on TP can exert strong coupling processes, especially the coupling interaction of cryosphere-hydrosphere-atmosphere with other spheres in the Earth system. It is believed to have a huge influence on regional and global weather climate and ecological environment (Bibi et al. ^[1]; Kang et al. ^[4]; Yao et al. ^[8]; You et al. ^[9]; Li et al. ^[10]; Yao et al. ^[12]; Chen et al. ^[14]; Liu et al. ^[22]; Zhang et al. ^[30]; Chan et al. ^[32]; Li et al. ^[33]).

It is very highly confident that the TP has been warming and wetting since the 1950s, and this trend beganto accelerate after the 1980s, which was indicated by multiple facts (Bibi et al. ^[1]; Kang et al. ^[4]; Li et al. ^[10]; Yao et al. ^[12]; Chen et al. ^[14]; Kuang et al. ^[20]; Zhang et al. ^[30]; Chan et al. ^[32]; Li et al. ^[34]; Wu et al. ^[35]; Yang et al. ^[36]; You et al. ^[37]). To better understand these evident changes, the precipitation variations over TP were one of the key concerns for the research front (Kang et al. ^[4]; Yang et al. ^[7]; Yao et al. ^[8]; Yao et al. ^[12]; Wang et al. ^[24]; Wang et al. ^[26]; Arias et al. ^[38]; Chen et al. ^[39]; Cuo et al. ^[40]; Ning et al. ^[41]; Wang et al. ^[42]; Zhu et al. ^[43]; Zhu ^[44]). Numerical models have been widely used to investigate precipitation variations. Therefore, many studies have been devoted to improve the performance of precipitation simulation. The terrain/land use, the cumulus parameterization scheme (CUs), the planetary boundary layer parameterization scheme (PBLs), and cloud microphysical parameterization schemes (MPs) in numerical simulations, including the sensitivity to the accuracy of the precipitation have also been investigated (Gao et al. ^[3]; Li et al. ^[5]; Wang et al. ^[24]; Chen et al. ^[45]; Chen et al. ^[46]; Fu et al. ^[47]; Gao et al. ^[48]; You et al. ^[49]; Li et al. ^[50]; Lin et al. ^[51]; Norris et al. ^[52]; Ou et al. ^[53]; Tong et al. ^[54]; Wang et al. ^[55]; Wu et al. ^[56]; Xu et al. ^[57]; Yu et al. ^[58]). However, there are still few studies on the comprehensive optimization of multiple parameters with the same model version for a large area over the TP and the surrounding regions. In addition, due to the special characteristics, the precipitation mechanism and its evolution process in this region appear more complicated, especially under the rapid retreat of the cryosphere that has huge effects and feedback on surface energy balance, water vapor flux, clouds, etc. (Kang et al. ^[4]; Yang et al. ^[7]; Yao et al. ^[8]; You et al. ^[9]; Yao et al. ^[12]; Niu et al. ^[59]; Pepin et al. ^[60]). It still faces challenges for numerical models and reanalysis to describe the precipitation over TP (e.g., Tong et al. ^[6]; Chen et al. ^[45]; Bai et al. ^[61]; Liu et al. ^{[62, 63]}; Wang et al.). One factor contributing to these uncertainties is associated with the coupled processes of the special topographic forcing and the interaction of strong ice/snow-atmospheric feedbacks.

To the sophisticated description of "polar environmental characteristics" over TP, such as better heat transfer for ice/snow and more suitable model settings. Here, the polar-optimized WRF model V4.1.1 (Polar WRF) is used for downscaling simulations over the TP and its surrounding regions. This study aims to (1) investigate the impacts of multiple options to study the "better-performing" configuration suite enhancing the accuracy and reliability representation of precipitation for using Polar WRF; (2) generally evaluate whether precipitation can be reasonably simulated with this polar-optimized model over the TP by a full calendar year and thus are more statistically robust than sensitivity test and case studies. This investigation is also an extension of the utilization and test of the Polar WRF (V4.1.1) beyond the Antarctic and Arctic (Xue et al. ^{[65, 66]}). It follows the series of verification studies documenting ongoing Polar WRF comprehensive evaluations and focuses on the "hottest topics" of modeling precipitation over the TP. The next sections are organized as follows: Section 2 concentrates on the Polar WRF model description. Section 3 depicts the experiments, model configuration, and data used in this study. The results of the simulation are performed in Section 4. Finally, the conclusion and discussion are provided in Section 5.

The polar version of the regional WRF model named "Polar WRF" is being developed, tested, and applied to climate variability and change problems in polar regions (e.g., Xue et al. ^[65]; Bromwich et al. ^{[67, 68, 69]}; Hines ^[70]; Hines and Bromwich ^[71]; Wilson et al. ^{[72, 73]}). Polar WRF is based on the standard WRF (ARW) and has been improved by the Polar Meteorology Group at the Byrd Polar and Climate Research Center of Ohio State University. The primary task of the modifications is to make the model keep state-of-the-art skills in operational weather forecasting and research simulation over polar regions (Xue et al. ^[66]; Bromwich et al. ^[67]; Wilson et al. ^{[72, 73]}; Kumar et al. ^[74]; Powers et al. ^{[75, 76]}; Walsh et al. ^[77]). Hence, it focuses on the "polar environmental characteristics" such as optimizing heat transfer for ice sheets and revised surface energy balance in the Land Surface Model (LSM), comprehensively describing the ice, snow process, and optimal treatment of the ice/snow-air feedback, improving and testing cloud microphysics for polar conditions. The improvements are present in different Polar WRF versions, and some are also absorbed by the rolling update of standard WRF (Xue et al. ^[65]; Bromwich et al. ^{[68, 69]}; Hines et al. ^{[71, 70][78, 79]}; Hines and Bromwich ^[80]).

Upon the release of WRF 4.X, the Polar WRF was pushed to V4.1.1 in 2019 with the significant new features modified Noah-MP and the improvement in the downward shortwave radiation implemented that compared with the Polar WRF V3.X. This version of the model serves good performance both in the Antarctic and Arctic. The results of the benchmark over the Antarctic and the Southern Ocean indicate that the model has good skills for variables, including surface and vertical of the atmosphere, for intra-annual variations in large-scale atmospheric circulation (Xue et al. ^[65]). For the polar mesoscale system, the model was applied over the Arctic region by simulating polar lows that demonstrated good performance (Xue et al. ^[66]). Considering the model performance under multiple polar conditions and the better utilization experience needed to investigate further, the Polar WRF V4.1.1 is continuously utilized in this study.

Topography, land cover, and vegetation can have significant impacts on processes such as orographic drag, albedo, radiation, surface roughness, soil temperature, and humidity, which play crucial roles in the heat and moisture fluxes conduction and exchange between the land surface and atmosphere (Skamarock et al. ^[81]). Finer represents the land surface information such as terrain, land use, and vegetation in regional models is highly confident of improving model skills (Cao et al. ^[82]; Li et al. ^[83]; Li et al. ^[84]; Li et al. ^[85]; Nahian et al. ^[86]; Schicker et al. ^[87]; Teklay et al. ^[88]; Yin et al. ^[89]). To better capture the characteristics over TP, the new terrain and land use are used to replace the default static data that was released with the model. The SRTM (Shuttle Radar Topography Mission) (Farr et al. ^[90]) 3 arc-second digital elevation data (the spatial resolution of the data is ~90 m×90 m) is applied (Fig. 1) to replace the SRTM 30 arc-second terrain data (spatial resolution of ~900 m×900 m) (De Meij et al. ^[91]). Even on a grid with a horizontal resolution of 10 km, there are significant differences (exceeding±300m) between the two in steep mountain areas of the TP. More refined land use data of TP from the National Tibet Data Center (TPDC) (Xu ^[92]), with a spatial resolution of 300 m×300 m, to optimize the model's default land use (Schicker et al. ^[87]). Compared with the original land use, the classification of permafrost and tundra has been added to the land surface data, and the main distribution areas are the Altun Mountains and the Kunlun Mountains in the northern part of TP. In addition, the representation of the land surface, such as the water body (lakes), ice, and snow on the plateau, also can be more vivid (figure not shown).

The model sets a single domain with 360×300 grid points on the Lambert projection (centered at 32.5ºN, 87.5ºE) and 10 km horizontal resolution, which covers the TP and its surrounding regions in western China (Fig. 1). In the vertical direction, the model utilizes HVC hydrostatic pressure coordinates with 71 vertical levels and reaches the top of 3 hPa, as it has been noticed that a higher model top provides better treatment of upward propagating gravity waves (e.g., Xue et al. ^[65]; Wilson et al. ^[73]; Bromwich et al. ^[93]). The optimized Noah-MP (Xue et al. ^[65]) is used for the LSM to sufficiently perform land-atmosphere interaction. The Rapid Radiative Transfer Model (RRTMG) (Iacono et al. ^[94]) is selected for longwave and shortwave radiation that relies upon a wide range of previous tests (e.g., Xue et al. ^[65]; Bromwich et al. ^{[67, 68, 69]}; Hines and Bromwich ^[71]; Wilson et al. ^[72]). There are more variations for other physical processes and closely related to sensitive tests. The details will be discussed in the next section.

Figure 1.  Model domain overlaid SRTM 3 arc-second terrain data at 10 km resolution. The black dots indicate the location of the in situ observations from TPDC. The gray dash lines along at 34ºN and 95ºE are used to divide the domain into four sub-regions (NW, Northwest; NW, Northeast; SE, Southeast; SW, Southwest).

Here, the monthly running method is carried out to better take into account the impacts of external forcing on the regional climate (Xue et al. ^[65]). However, to pursue a more economical way, the spin-up time is shortened to 5 days; that is, the model runs continuously for 33–36 days per month and outputs 3 hours at 0000, 0300, 0600, 0009, 0012, 0015, 0018, 0021 UTC for each day. The first 5 days at 0000 UTC before the 1st of each month are used for the mode spin-up, and the rest are the simulation results. The ERA-Interim (Dee et al. ^[95]) 60 model levels reanalysis is employed to initialize and drive model runs as the higher vertical resolution driving data showed better performance than those of less such as ERA5 (Hersbach et al. ^[96]) with 37 pressure vertical levels (Xue et al. ^[65]). As the reanalysis did not resolve the lakes over the TP, it is often recommended to solve the "hot lake" issue by the average surface air temperature used to adjust the lake surface temperature instead of the default that was extrapolated from the nearest resolved water bodies (e.g., Gao et al. ^[97]; Heath et al. ^[98]).

One target of this study is to investigate the impacts of multiple options on the precipitation over TP and the surrounding regions. Considering configurations are highly correlated with each other; therefore, based on land static data optimization and the optimized Noah-MP, 4 categories of 19 sensitive experiments are comprehensively examined to further discuss the model performance, including the nudging setting, 5 CUs, 2 PBLs, and 6 MPs. All these numerical experiments are carried out in 2015 because more in situ observation and the new land use over TP from the TPDC are available during this period. Table 1 shows the summary for each experiment.

Nudging can force the model simulation towards a series of analyses grid-point by grid-point (Skamarock et al. ^[81]), which is used to reduce model drift (e.g., Xue et al. ^[65]; Bromwich et al. ^[67]; Hines et al. ^[79]; Glisan et al. ^[99]), particularly for simulations of extended periods (e.g., months). Inherited from the WRF, Polar WRF can afford to nudge u, v, t, and q from driving fields with spectral and grid nudge methods. Previous studies have shown that model skills are very sensitive to nudge settings when long-term runs or dynamic aspects need to be addressed (Xue et al. ^[66]; Cha et al. ^[100]; Liu et al. ^[101]; Otte et al. ^[102]; von Storch et al. ^[103]). Therefore, the 6 experiments (nudging category in Table 1) for methods and nudging coefficients are investigated to increase the knowledge of nudging settings.

Several excellent investigations have identified that CUs schemes have a considerable impact on the simulation of precipitation over TP and nearby areas (Gao et al. ^[3]; Gao et al. ^[48]; You et al. ^[49]; Lin et al. ^[51]; Ou et al. ^[53]; Wu et al. ^[56]; Yu et al. ^[58]; Gao et al. ^[97]; Chen et al. ^[104]; Gao et al. ^[105]; Qian et al. ^[106]; Zheng et al. ^[107]). Referring to these findings, here we mainly study the influence of CUs settings by the following 3 typical kinds of simulation (a total of 5 experiments). (1) The Grell-Devenyi Ensemble Scheme (GD) and its improved version Grell-3 scheme (Grell et al. ^[108]). The GD combines multiple schemes and variants run within each grid point and gives an average feedback result. Currently, each grid box has 144 members in this scheme. Grell-3 (G3) is based on GD and is especially suitable for high-resolution simulations that should not be coarser than 10 km (Skamarock et al. ^[81]). (2) Kain-Fritsch schemes, including the Modified Kain-Fritsch (KF) (Kain ^[109]) and Multi-Scale Kain-Fritsch (MKF) (Zheng et al. ^[107]). KF using the mass flux with downdrafts and CAPE removal time scale for deep and shallow convection sub-grid scheme. The MKF stemmed from the original KF scheme and introduced new functions, such as resolution-dependent that adapt to the time scale of the convective process. (3) The NoCU represents the CUs that are not used, which also achieves good results at the "gray zone" grid spacing (~4–10 km resolution) (Ou et al. ^[53]; Qian et al. ^[106]).

For the PBLs, the Mellor-Yamada-Nakanishi-Niino (MYNN) Levels 2.5 (MYNN2.5) (Nakanishi et al. ^[110]) and the Yonsei University schemes (YSU) (Hong et al. ^[111]) are examined here. MYNN2.5 is tested by the Polar Meteorology Group, which is outstanding for polar environmental features and orographic influences (e.g., Xue et al. ^[65]; Bromwich et al. ^[67]; Hines et al. ^{[78, 79]}). Moreover, the improvement of the reduction in the downward shortwave radiation bias has been verified as a significant improvement of Polar WRF V 4.1.1 (Xue et al. ^[65]; Olson et al. ^[112]). For the TP, YSU also has been frequently used (Ou et al. ^[53]; Gao et al. ^[97]; Chen et al. ^[104]; Gao et al. ^[105]; Qian et al. ^[106]). More importantly, when the MKF scheme is selected for CUs, the PBLs matched with YSU; that is, only the combination of MKF +YSU can be used in the model right now (Skamarock et al. ^[81])

MPs dealing with cloud and precipitation particles must be parameterized in the numerical weather and climate models. The parameterizing of microphysics remains highly challenging because of the complexity of the underlying physics and the lack of understanding of these coupling processes (Bauer et al. ^[113]; Benjamin et al. ^[114]; Randall et al. ^[115]). Hence, a reasonable selection and modification of MPs are expected to achieve better simulation results of clouds and precipitation. Here, MPs are studied based on the research advances that have been extensively tested in the Arctic, Antarctica, and TP (Gao et al. ^[48]; You et al. ^[49]; Xue et al. ^{[65, 66]}; Bromwich et al. ^[68]; Hines et al. ^[78]; Hines and Bromwich ^[80]; Gao et al. ^[97]; Qian et al. ^[106]; Zhu et al. ^[116]; Li et al. ^[117]; Ma et al. ^[118]; Shen et al. ^[119]; Yang et al. ^[120]; Zhang et al. ^[121]). Modify the liquid water droplet concentration in Morrison 2-mom to 50 cm^-3 which has been verified to significantly improve the performance of polar clouds such as supercooled water in clouds (Xue et al. ^[66]; Hines et al. ^[78]; Hines and Bromwich ^[80]). TP, similar to the Antarctic and Arctic, is a colder region, and the liquid water of clouds more easily exists in the form of supercooled water. For the TP, the Morrison 2-mom specified droplet concentration of 100 cm^–3 (referred to as M100) because of a greater amount of could particles (probably twice for maximum) and more prone to supercooled water than the Arctic (Yan et al. ^[122]). P3 (Morrison et al. ^[123]) used a new approach to predict the transformation of ice particles, named "one ice category, " which represents a combination of ice, snow, and possibly graupel. The setting of double-moment for rain, ice, and cloud water and supersaturation-dependent activation is applied in this study (P51). Thompson aerosol-aware (Thompson et al. ^[124]) has been modified for mid-latitude convective, orographic, and snowfall conditions and considers water- and ice-friendly aerosols. WRF Double-Moment 6-class scheme is the double moment for warm-rain, and other processes are basically the same as WRF Single-Moment 6-class scheme (WSM6) (Lim et al. ^[125]). The condensation nuclei (CNN), such as cloud water and rain, adopt the "prediction method." It performs reasonably well for high-resolution simulation studies in China (Zhu et al. ^[116]; Li et al. ^[117]; Shen et al. ^[119]; Zhang et al. ^[121]). Compared with the two-moment schemes, single-moment schemes such as WRF Single-Moment 5-class scheme(WSM5) and WSM6 are often used as they can maintain a good balance between model skills and computational resources. And the WSM6 is also used and tested by the GRAPES_Meso operational weather forecasting model (Shen et al. ^[119]; Zhang et al. ^[121]; Huang et al. ^[126]; Ma et al. ^[127]). Table 1 shows a summary of the microphysics mentioned above, the hydrometeors, and whether ice-phase/mixed-phase processes are also included.

Considering the limited performance of observation and reanalysis over TP, multiple datasets are used to evaluate the simulations. The in situ observation obtains the daily precipitation dataset located in the TP and surrounding areas supplied by the TPDC (1979–2015, about 184 stations Fig. 1) (China Meteorological Data Network ^[128]). Satellite precipitation estimation products used the Tropical Rainfall Measuring Mission (TRMM) 3B42 V7 dataset with 3-hourly and 0.25° × 0.25° resolution (Huffman et al. ^[129]). CRA-RA/Land, the first generation global atmosphere and land reanalysis released by the China Meteorological Administration. Probably assimilation of various conventional and satellite observations over East Asia, CRA-RA/Land is in good agreement with and even outperforms other reanalyses such as ERA5 datasets in near-surface variables over China (Li et al. ^[130]; Shen et al. ^[131]; Yang et al. ^[132]; Yu et al. ^[133]; Zhang et al. ^[134]). Precipitation estimates from CRA-Land with a horizontal resolution of ~34 km and a 6-hour temporal interval. The ERA-Interim precipitation has been investigated as a benchmark and is also used in this study. The data is provided every 3 hours with a horizontal resolution of 0.75°×0.75°.

Due to the different spatial resolutions, the grid data are extracted and horizontally interpolated onto the domain by the WRF Preprocessing System (WPS) (e.g., Xue et al. ^[65]; Bromwich et al. ^[67]) with the same re-grid method. It is beneficial for comparisons using the same resolution (10 km) and grid points through the utilization of bilinear interpolation to interpolate grid data to station locations when comparing data between grid and in situ.

In this section, the multiple aspects of the model option for precipitation over the TP nearby areas are examined by the sensitivity experiments. Considering the main targets of this study and successor work, the variations of monthly and daily precipitation between multiple model settings are evaluated against in situ, satellite observations, and reanalyses. It focuses on the spatiotemporal comparison of large-scale and average conditions and does not analyze the specific rainfall processes, such as comparing rainfall forecasts by the Threat Scores (Jiang et al. ^[135]; Mesinger ^[136]; Yang and Tung ^[137]). Here, January and July, the two months with the largest difference in a year are used to investigate. Since the precipitation is mainly occupied in the summertime, the analysis of the precipitation characteristics of each experiment in July 2015 is discussed below.

Figure 2 shows the results of precipitation in July 2015 under different nudging settings and the references with ERA-Interim (EI), CRA-Land (CRA), and TRMM 3B42 (TRMM), respectively. Overall, the model well captures the main characteristics of precipitation distribution over the domain. Precipitation decreases from southeast to northwest, and there is a clear borderline along the edge of the plateau (such as an altitude of 3000 m). The Bay of Bengal has the most precipitation amount, with a maximum value of more than 800 mm, whereas the Tarim Basin is the region with the least, with most areas less than 10 mm. At the same time, it can also be seen that EI has the largest precipitation, with most areas on the plateau having more than 50 mm during this month. CRA and TRMM are roughly similar, and the covering of precipitation greater than 50mm is significantly smaller than EI. For TRMM, we also noticed that the precipitation in the northern part of the domain (north of 35ºN, such as the Tianshan area) is insufficient. It is likely to be related to the low coverage of TRMM data in extratropical regions.

Figure 2.  Precipitation of July 2015 over the model domain. (a)-(c) the model results with specified nudging settings. (a) Spectral nudging wave number was specified as 17, (b) the wave number was specified as 7, and (c) grid nudging; (d)-(f) ERA-Interim (EI), CRA-Land (CRA), TRMM 3B42 (TRMM). The black solid line indicates the altitude of 3000 m.

Further analysis shows that the results of spectral nudging (Fig. 2 (a) and (b)) located a very weak anomalous precipitation belt along the southern side of the Himalayas. This belt is also accompanied by an abnormally strong surface short-wave down and weak long-wave radiation (not shown), which indicates that the low cloud cover and the convective process might be restricted. In addition, the tests also found that the intensity of this abnormal precipitation belt is related to the wave number used. For instance, when the wavenumbers are selected as 17, 7, and 3 (the cutoff wavelength is roughly equal to ~350 km, ~800 km, and ~1800 km, respectively), the belt gradually blurred with the wavenumber decreases. In other words, the model has more freedom as the wave number decreases, and the impacts of nudging on medium scales are also weakened. That means the modulation for overestimation precipitation is suppressed. The grid nudging result is different from spectral nudging, which is consistent with EI, CRA, and TRMM maintaining a strong precipitation belt on the south side of the Himalayas. Some possible reasons can be inferred as follows. The spectral nudging skills are determined by the nudge scale, which is related to the model domain and the scale of driving fields (Skamarock et al. ^[81]; Liu et al. ^[101]; Miguez-Macho et al. ^[138]). The southern foothills of the Himalayas are more sensitive to scale due to the complex local topographical and convective variations. Thus, the skill of spectral nudging really penalizes here. However, the option of grid nudging can avoid "scale selection" and good performance with precipitation and moisture adjustment (e.g., Bowden et al. ^[139]; Wootten et al. ^[140]). Therefore, it shows better performance.

Nudging toward the water vapor mixing ratio is believed to be beneficial for the simulation of precipitation and moisture (e.g., Mai et al. ^[141]; Spero et al. ^[142]). Therefore, the impacts of nudge q and its strength are examined. First of all, sensitive test results indicate nudging q can significantly improve precipitation on the plateau. In the central and northern parts of the plateau, where precipitation is overestimated (such as exceeding 50 mm), there is a significant decrease (not shown). More details are shown in Fig. 3. It gives the difference in monthly precipitation between model results with specified grid nudging settings and EI, CRA, and TRMM, respectively. Although the simulation precipitation minus each dataset shows more variation, it is consistent that the superfluous precipitation can be reduced by nudging q, and the greater the intensity of the nudging coefficient, the more significant the effect on the reduction of precipitation. Compared with EI, the model mainly shows a negative difference, especially along the western and southern mountain margins, where there is a precipitation belt with negative 80–100 mm. This is probably caused by the overestimation of precipitation over these areas from EI. Both the differences in the model minus CAR and TRMM indicate the negative values in the west while positive in the central and eastern parts of the plateau. The maximum positive region is located at the east edge of the plateau and the Hengduan Mountains. It is consistent with the general inability of current models to simulate precipitation in this region (Chen et al. ^[45]). It is worth noting that with increasing the intensity coefficient of nudging q, the positive trend on the central and eastern TP will decrease, and the changes are intensifying correspondingly in the western plateau. Thence, trade-offs for the obvious regional features possessed by nudging q. The intensity coefficient is 0.0002 for q and other settings, including nudging u, v, t, and ph, using the same value of 0.0003. All the nudging applied above the vertical 35 layers higher than the top of the PBL reasonably improved the model performance at the upper layers (Xue et al. ^[65]; Mai et al. ^[141]).

Figure 3.  Precipitation difference in July 2015 with specified grid nudging settings. Model minus ERA-Interim (EI), CRA-Land (CRA), TRMM 3B42 (TRMM), respectively. The columns from left to right are turned off nudging q, nudging q coefficients as 0.0001, 0.0002 and 0.0003, respectively. The black solid line indicates the altitude of 3000 m.

Since EI significantly overestimates precipitation in the southern TP and TRMM underestimates in the north, Fig. 4 only shows the monthly precipitation difference of the model minus CRA for July 2015. At present, MKF can only be combined with YSU. Therefore, the suite of cumulus with MYNN2.5 and YSU are tested, respectively. It is clear that the model simulated with CUs usually produced too much precipitation, but NoCUs seems precipitation is not enough conversely. Furthermore, the impacts of different CUs on precipitation are also evident. Fig. 4a-c shows that G3 has the most significant positive difference in precipitation, and KF produces too much precipitation covering large areas such as the east and south of TP (the Hengduan Mountains). We also noticed that the G3 scheme produced more precipitation than GD (Fig. 4b-c), which had been investigated by the simulation analysis of the precipitation over the southern slopes of the TP (Wu et al. ^[56]). The MKF performed much superior to KF (Fig. 4d-e). It means the scale adaptive cumulus scheme shows better skills for precipitation over the TP. In other words, it is that making use of scale-aware performs superior for solving interaction/uncertainty problems, especially for simulations forced by complex terrain and low confidence in initial and lateral boundary conditions. Therefore, MKF is suggested for priority use.

Figure 4.  Precipitation difference in July 2015 from the model with specified CUs and PBLs. Model minus CRA-Land. The black solid line indicates the altitude of 3000 m.

Figure 4 also remarkably indicates that GD +MYNN2.5 and MKF+YSU are better at simulating precipitation. And MKF+YSU seems more sophisticated as over the eastern side of the TP (such as west of Sichuan Basin), the central part (the Three-River headwater region), and northern India (figure omitted), the precipitation performed by GD+MYNN2.5 more superfluous than that of MKF+YSU. Due to the different setting suites, KF+YSU and KF+MYNN2.5 served as supplementary to determine MKF and YSU, which one improved model skills more significantly. Compared with MKF+YSU (Fig. 4d), KF+YSU (Fig. 4e) overestimated precipitation at central and eastern TP. It indicates that MKF has better performance than KF. Correspondingly, the results of KF+MYNN2.5 (Fig. 4a) and KF+YSU (Fig. 4e) revealed that MYNN 2.5 shows the ability to improve precipitation slightly compared to TP. However, as the MKF performs greater modeling skills, the MKF+YSU performs much better. In other words, compared with MKF, both YSU and MYNN2.5 could have a limited impact on precipitation simulation.

Based on the optimized settings for nudging, CUs, and PBL, microphysics options can be further examined to expand the knowledge for a "better-performing" configuration suite. Fig. 5 compares the spatial distribution of the precipitation difference from the model with specified microphysics minus CRA. Significant differences indicate that the model configured with different microphysics offers a wide variety of capabilities or limitations. Among the four double-moment schemes, M100, P51, and TPA carried similar results. Compared with CRA, precipitation is overestimated in the eastern slope of TP (such as the western Sichuan Basin), the southern foot of the Nyainqentanglha Mountains, and the western slope of the Hengduan Mountains, etc. On the contrary, in the Tianshan Mountains and other places in the west of the TP (west of ~85ºE), the modeling precipitation is less than that of the CRA. Specifically, M100 is close to P51. However, the negative values of M100 in the west of the plateau are smaller, and it is much closer to CRA. TPA is slightly poor, and too much precipitation is produced in the southeast, such as Yunnan and Guizhou. For WDM6, it is evident precipitation is underestimated in the east of TP while slightly overestimated in places such as the Kunlun Mountains in the northwest of TP. In the two single-moment schemes, WSM5 is similar to WDM6 but with significantly less precipitation. WSM6, which extended the WSM5 by adding additional processes related to graupel (Hong ^[143]), performs much better than WSM5 and is roughly comparable to M100 and P51.

Figure 5.  Precipitation difference in July 2015 from the model with specified MPs. Model minus CRA-Land. The black solid line indicates the altitude of 3000 m.

To evaluate the finer differences, the RMSE (Root Mean Square Error) and precipitation difference from model data against in situ observations were calculated. Fig. 6 shows the RMSE value and spatial distribution from the model with specified MPs accompanied by EI, CRA, and TRMM that are used as references. First, RMSE exhibits increasing from the northwest (values are usually less than 4 mm) to the southeast (most are greater than 10 mm). In particular, the maximum RMSE is mainly distributed in the Hengduan Mountains and on the edge of the Himalayas. One aspect is likely related to the higher absolute precipitation that produces larger biases. On the other hand, it is still challenging for the model to simulate precipitation under the complex terrain and strong water vapor transport conditions over the TP and surroundings. Compared with different MPs, M100 and WDM6 show less RMSE than the other two double-moment schemes at the southeast of TP. The single-moment schemes indicate a similar spatial pattern to that of double-moment schemes, and WSM5 shows less RMSE than WSM6. Furthermore, we analysed the differences in monthly precipitation. Although WDM6 and WSM6 perform better RMSE, precipitation is systematically underpredicted compared to in situ observation (Fig. 7), which is consistent with the previous findings compared with CRA (Fig. 5). By comparison with in situ dataset, Fig. 6 and Fig. 7 also verify that the CRA shows the best performance in the 3 kinds of grid precipitation products.

Figure 6.  Precipitation RMSE (colored dots) in July 2015 with from the model with specified MPs and ERA-Interim (EI), CRA-Land (CRA), and TRMM 3B42 (TRMM), respectively. Model/Reanalysis minus in situ observation. The black solid line indicates the altitude of 3000 m.

Figure 7.  Precipitation difference (colored dots) in July 2015 from the model with specified MPs and ERA-Interim (EI), CRA-Land (CRA), and TRMM 3B42 (TRMM), respectively. Model/Reanalysis minus in situ observation. The black solid line indicates the altitude of 3000 m.

Considering the significant spatial variation of precipitation differences through the TP, the region is divided into four sub-regions along 34ºN (Li ^[144]) and 95ºE (Gao et al. ^[3]). That is, Northwest (NW), Northeast (NE), Southeast (SE), and Southwest (SW), and the number of in situ stations in each sub-region are 24, 53, 68, and 29, respectively. Fig. 8 gives the whole TP and its sub-regions of station-averaged precipitation in July 2015 from the model with specified MPs. It can be seen that among the four sub-regions, the precipitation in the SE is the largest while the NW is the smallest, which is about one-tenth of SE. NE and SW are equivalent, probably half of the SE. These characteristics are consistent with the spatial distribution discussed above that the precipitation decreases from southeast to northwest over TP. The difference is mainly in the southern TP (south of 34°N), especially in the SE and SW, where MPs perform significant impacts on precipitation. For instance, in SE where is the most significant precipitation difference, the maximum (TPA 201 mm) is more than 2.5 times larger than the minimum (WSM5 77 mm). From Fig. 8, the overarching findings are that the M100, WSM6, P51, and TPA are relatively overestimated precipitation, whereas WSM5 and WDM6 are underestimated precipitation. It also can be seen that M100 is superior to others and is closest to CRA and in situ observation, respectively.

Figure 8.  Station-averaged precipitation in July 2015 from the model with specified MPs, in situ observation (OBS), ERA-Interim (EI), CRA-Land (CRA) and TRMM 3B42 (TRMM), respectively (units: mm).

Standing on the "better-performing" configuration investigated above, the performance of precipitation simulation using Polar WRF can be confidently further discussed in the following sections. The refined model parametrizations are as follows: grid nudging is used, where the q intensity coefficient is 0.0002 while the rest are 0.0003; the MKF cumulus parameter scheme; the YUS boundary layer parameter scheme and the M100 as the microphysics. Here, as mentioned in the introduction, a full calendar year simulation and compared with multiple reanalyses and observations are carried out.

Figure 9 shows the annual precipitation over TP and the surrounding regions from simulation, EI, CRA, TRMM, and in situ observations. All the results show an evident similar pattern in that precipitation decreases from southeast to northwest, as discussed above. Further comparison with station observations, EI exhibits more precipitation over the domain while TRMM does not have enough precipitation north of ~35ºN. For CRA, most areas are very close to the observation except for some stations on the edge of 3000 m of the southeast. Although the precipitation on the western slope of the Hengduan Mountains and the eastern slope of TP is relatively heavy, the distribution of M100 performs consistently with observations. By examining the annual precipitation frequency, the abnormally high frequency of precipitation can be seen in the EI, CRA, and M100. But M100 shows much better performance than the two reanalyses, such as the western slopes of the Hengduan Mountains, the eastern slopes of TP, and the Tianshan Mountains (Fig. 10). Corresponding to these features, the annual precipitation intensity from models and reanalyses is also weaker than observations. M100 presents superior to the two reanalysis data, especially in the southeast of TP (Fig. 11).

Figure 9.  The distribution of annual precipitation over the Qinghai-Tibet Plateau and the surrounding regions. The black solid line indicates the altitude of 3000 m.

Figure 10.  The distribution of annual precipitation frequency over the Qinghai-Tibet Plateau and the surrounding regions. The black solid line indicates the altitude of 3000 m.

Figure 11.  The distribution of annual precipitation intensity over the Qinghai-Tibet Plateau and the surrounding regions. The black solid line indicates the altitude of 3000 m.

Monthly and daily precipitation can describe more spatial and temporal variations and have been used to examine precipitation performance (e.g., Bai et al. ^[61]; Li et al. ^[130]; Jeworrek et al. ^[145]). Fig. 12 shows monthly station-averaged precipitation, including the whole TP and its sub-regions from January to December 2015. Precipitation has significant annual temporal variation and regional spatial characteristics. The precipitation is mainly concentrated in the summer half year. The amounts from May to September account for 70%–80% of the annual precipitation. Among them, in situ observation (OBS) is 80.0%, M100 is 74.4%, EI is 72.2%, CRA is 79.5%, and TRMM is 80.3%. Moreover, the precipitation from April to October contributed more than 90% of the annual precipitation (only EI accounted for 87.4%, and the rest ranged from 90.2% to 92.9%). This intra-annual variability is in step with the outbreak, maintenance, and retreat of the Asian summer monsoon.

Figure 12.  Monthly station-averaged precipitation from January to December in 2015. (a) NW, 24 stations; (b) NE, 57 stations; (c) SE 68, stations; (d) SW, 29 stations; (e) TP, 178 stations.

For the three precipitation products, CRA and TRMM are relatively close to those observed by the station, and CRA is slightly larger. The precipitation from TRMM in the northern TP (NW, NE) is somewhat underestimated, whereas that in the southern TP (SE, SW) is slightly overestimated. This is also the reason why TRMM has the smallest statistical precipitation difference between the whole TP and station observations. Obviously, EI shows too much precipitation, especially in the sub-regions of SE and SW. Two reasons might be related to this. One, EI has relatively large precipitation in these areas, so the ability to represent precipitation in complex terrain areas such as the south is limited. The two horizontal resolution of EI is the lowest in the four grid data used here; it may produce larger errors than others after interpolation.

The model's monthly precipitation is evidently better than EI and generally close to CRA and TRMM, even particularly slightly superior to CRA and TRMM in the NW and SW sub-regions. However, over the SE, the monthly precipitation trend is quite well depicted by the model, but the simulated precipitation is still too overestimated. This is probably due to the sparse observations and the limited direct observation data available for data assimilation at the NW and SW. CRA failed to show a better advantage, while the downscaling of the model has certain advantages. In particular, the model optimizes the ice-snow-atmosphere interaction and performs better in the area with large ice and snow coverage in the west of TP. On the contrary, observations in the eastern region are relatively abundant, and CRA performs better in this region as richer observations have been assimilated (Li et al. ^[5]; Li et al. ^[130]). In addition, since the initial and lateral boundary conditions derived from EI, it will also lead to large positive abnormalities where precipitation is overestimated dramatically, such as in the southeast of TP. These findings also suggest the importance of data assimilation (improvement of initial conditions) and model-driving data for the regional model.

Unifying the mode results in the daily temporal resolution as the in situ, Fig. 13 displays the time series of daily station-averaged precipitation in 2015. The model shows an excellent ability to capture the daily precipitation evolution during the annual cycle. It skillfully describes the precipitation characteristics over TP and nearby areas. In the view of temporal distribution, precipitation is mainly produced from April to October (90–300 days on the X-axis in Fig. 13), which is consistent with the above finding. The daily precipitation also shows good agreement between the model and observations. The disadvantage is that the peak of daily precipitation simulated by the model is smaller than the observation. It indicates that the model is limited in simulating some heavy/extreme precipitation, which also caused the shortcoming of weak intensity for precipitation mentioned in 4.2.1. In addition, it can also be seen clearly that the daily precipitation in the SE from April to July (90–210 days on the X-axis in Fig. 13) is greater than the observation. Which further verifies the results discussed above.

Figure 13.  The time series of daily station-averaged precipitation in 2015. The grey solid line and black dashed line indicate observation and model output, respectively. (a) NW, 24 stations; (b) NE, 57 stations; (c) SE, 68 stations; (d) SW, 29 stations; (e) TP, 178 stations.

To combine multiple data into a comprehensive evaluation, Fig. 14 shows the daily precipitation of the Taylor diagram for M100, CRA, and TRMM, respectively. The daily precipitation simulated by the model (black dot) is at the forefront among the results, and the correlation of each region is in the range of 0.75–0.86 (Table 2). which is better than others. Also, the distribution of each region exhibited in the figure is relatively concentrated, which indicates the model skills in various regions of TP are relatively stable. In general, the model performance is comparable to that of CRA, especially in the NW and SW, which is better than the former. This is also supported by a quantitative comparison of the correlation, bias, and RMSE that are shown in Table 2. Although the monthly precipitation from TRMM is the closest to the observation, the analysis of daily precipitation shows that it is not superior to others. The results generally show a smaller correlation and larger RMSE (the blue box in Fig. 14). This suggests that TRMM is limited for characterizing the evolution of the refined precipitation. EI (the red triangle in Fig. 14) performs worst, and the ability to describe the daily precipitation is extremely limited. Moreover, the distribution of each point is too widespread to draw the SW whose results are beyond the scope of the illustration

Figure 14.  Daily precipitation of Taylor diagram for Model (M100), ERA-Interim (EI), CRA-Land (CRA) and TRMM precipitation (TRMM) compared with what was observed in 2015. The direction of the azimuth cosine is the correlation, the distance to the azimuth is the normalized variance ratio, and the distance to the reference point (OBS) represents the RMSE of the observation. The poor results of EI in the SW area are not shown in the figure; that is, only 4 values are displayed.

Considering it is still facing challenges for numerical models and reanalysis to describe the precipitation over the TP and its surrounding region, this study extended the testing and applied the Polar WRF model to expand the knowledge of modeling precipitation.

Firstly, based on the model ground surface properties upgrade and the optimized Noah-MP, the impacts of nudging, CUs, PBLs, and MPs on summer precipitation are comprehensively examined by 4 categories of 19 groups of sensitive experiments with the 10 km horizontal resolution. The "better-performing" options suite is as follows.

(1) Nudging helps to quantify forecasting and modeling uncertainties. The spectral nudging skills are determined by the nudge scale which is more sensitive to the complex local topographical and convective variations over the TP and nearby areas. The option of grid nudging can avoid "scale selection" and good performance with precipitation and moisture adjustment. That is, grid nudging is suggested. By trade-offs, keep the balance of model "drift" and "freedom, " the intensity coefficient is 0.0002 for q and other settings, including nudging u, v, t, and ph using the common values of 0.0003.

(2) For CUs, the MKF is suggested for priority use because the scale-aware scheme performs superior in solving interaction/uncertainty for simulations forced by complex terrain and low confidence in initial and lateral boundary conditions. The ensemble scheme displayed the more complex situation. GD and G3 show different simulation skills, and GD is much better. However, due to the limited work here, further research is absolutely needed to verify why the G3 scheme causes significantly more precipitation on the east side of the TP.

(3) Compared with the CUs, both YSU and MYNN2.5 could have a limited impact on precipitation simulation. Although MYNN 2.5 also shows the ability to improve precipitation. However, because only the combination of MKF+YSU can be used in the model right now, YSU is selected for PBLs at the end.

(4) A wide variety of capabilities or limitations are shown in MPs' options. The overarching findings are that the M100, WSM6, P51, and TPA are relatively overestimated precipitation, whereas WSM5 and WDM6 are underestimated precipitation. It also can be verified that M100 performs more reasonably than others and is closest to CRA and in situ observation, respectively.

Then, precipitation simulation for a full calendar year and compared with multiple reanalyses and observations are used for further evaluation. The downscaling simulations indicate that Polar WRF can successfully capture the general features over TP and performs more sophisticated than EI and TRMM. For the monthly and daily precipitation, it is equal to or even slightly better than CRA, especially in the western part of the plateau where the observations are sparse and the ice-snow-atmosphere interaction is more significant. In addition, It is worth noting that the performance of multiple precipitation products works quite differently over the TP and the surrounding regions. Overall, CRA is better than TRMM, and TRMM is superior to EI. Therefore, multi-sources of blended precipitation should be considered as a reference in the future.

This study also extends the findings of model defects. Although the model exhibits reasonable skills, it still underperforms in simulating precipitation under the complex terrain and strong water vapor transport conditions over the TP and the surroundings. The shortcomings include precipitation overestimation and precipitation with higher frequency and weaker intensity, which certainly need to be further improved. Furthermore, due to the key purpose of this study, only limited testing for the model based on knowledge of using Polar WRF and WRF model. More model option suites and refined evaluations are yet to be carried out, such as multi-year simulations, higher convection-permitting resolution, and multiple physics ensembles. Also, the precipitation diurnal variation, the rainfall peak, and the duration of rainfall events should be taken into account in future work.

In addition, better-driving data should be further considered. As the higher vertical resolution driving data is beneficial for the Polar WRF model (Xue et al. ^[65]), the ERA-Interim 60 model levels reanalysis is employed to initialize and drive model runs. However, ERA-Interim has been verified to overestimate precipitation over the TP, especially in the sub-region of SE and SW. It also probably contributes to the overestimation of model precipitation in these regions, with the ERA-Interim being used as the global model field used for nudging.

To summarize, this study comprehensively investigates the impacts of multiple options. It provides valuable information about the simulation of precipitation across a variety of Polar WRF model metrics over the TP and its surrounding regions. Despite the inevitable limitations studies entail, hopefully, our findings can help Polar WRF users and developers with their model configurations and better understand the model characteristics beyond the Antarctic and Arctic.