1 Introduction

Slope and aspect are prominent features in many regions of the world and have important influences on complex terrain. The plain area only accounts for less than 20% of the total land of the world. Some studies (Spreen 1947; Burns 1953) have shown that, together with slope, aspect and surface elevation, these factors account for about 80% of precipitation in some mountainous areas. Therefore, there are a lot of studies on sub-grid slope and aspect. For example, Vico and Porporato (2009) proposed a simple realistic probabilistic description of local slope and aspect as a function of properties of the field of elevation changes applied to assess the effects of topographic features on direct solar radiation mean and standard deviation. Löwe and Helbig (2012) provided a quasi-analytical method to compute the sub-grid slope and aspect over complex terrain. A number of such parameterization schemes have been proposed in the recent years. Results from numerical simulations have demonstrated that including sub-grid heterogeneity parameterization schemes can improve model performance in simulating the details of atmospheric numerical models (ANMs, including regional and global weather and climate models).

Topography has an important influence on global atmospheric circulation (Hughes et al. 2009; de la Torre et al. 2015; Junquas et al. 2015). There has been extensive research into the role of large terrain in numerical weather prediction (NWP) (Dörnbrack and Schumann 1993; Doyle et al. 2013), but there are still a number of deficiencies of sub-grid topography unresolved, especially with regard to the sub-grid slope and aspect. The orographic characteristics parameterization schemes have been extensively explored in previous studies, such as average (Giorgi et al. 2003), standard deviation (Wallace et al. 1983), Laplacian operator (Ma et al. 2016) and characteristic coefficient (Dozier and Frew 1990). The advantage of these methods is that they can describe the terrain features of each sub-grid in detail. The disadvantage is that it requires a large amount of calculation and storage. For global weather and climate models, applying these parameterization schemes on each grid requires too much computation. Taking the range of 100 km as an example, each sub grid is 1 × 1 km2, it needs to calculate 105 times. While the transformation of statistical description is only one calculation toward model simplification. Therefore, taking each sub-grid topography for model calculation is not applicable to global areas. Second, these topographic features have no directionality; that is, they do not consider the factors of slope and aspect, so they are not reasonable for applying for radiation and topographic dynamic uplift.

Based on the above reasons, previous studies have examined slope and aspect sub-grid parameterization for radiation, such as 3D Monte Carlo photon tracing simulation (Chen et al. 2006; Liou et al. 2007; Lee et al. 2013). It is to develop a parameterization scheme to establish a multiple linear expression equation applied to solar fluxes (Lee et al. 2011). This scheme makes the simulation results very close to realistic increased surfaces but depends on observation data, and different multiple linear regression models are required under different terrains. Later, sub-grid topographic effects on solar radiation have been recognized and parameterized in a few regional weather and climate models (Hauge and Hole 2003; Zhang et al. 2006; Arthur et al. 2018; Gu et al. 2012), global climate models (Lee et al. 2015) and land models (Lee et al. 2019; Hao et al. 2021).

The above analysis shows that the effects of sub-grid topography are not explicitly included in ANMs, especially concerning dynamic processes. Nearly all global climate models neglect the sub-grid slope and aspect effects on dynamical lifting. Most ANMs adopt unobstructed horizontal surfaces into dynamic processes. Such simplified parameterizations do not account for sub-grid topographic effects and can lead to large systematic biases in simulating weather and climate processes over complex terrain. In addition to the qualitative influence of the slope and aspect on solar radiation, it has been demonstrated that the impacts on topographic uplift are of great significance (Rontu 2006). For example, California's mountain fires are largely related to the unique terrain and climate conditions in this region. California is located on the Pacific coast west of the Rocky Mountains in the western United States, with a large height difference. When there is a strong cold and high pressure in the mid-west of the United States going south to the Rocky Mountains, due to the strong pressure gradient, a strong northeast wind will blow along the coast of California. This wind falls directly from the Rocky Mountains to the Pacific coast. The wind sinks along the hillside across the mountains, forming a dry adiabatic sinking process, resulting in a significant increase in the temperature and low humidity of the air mass after crossing the mountain and forming a “foehn” (Krick 1933, 1934). There are similar “foehn effects” in almost all mountains around the world, such as the Southern Alps in New Zealand (McGowan and Sturman 1996), the Andes in Chile (Antico et al. 2021) and the Andes in the United States (Hoinka 1985), causing tremendous human and property loss (Richner and Hächler 2013).

In addition, the sub-grid slope and aspect are also of great significance in topographic precipitation (Leung and Ghan 1995). For large mountains, the precipitation increases on the windward side, and rain shadows form on the leeward side. It is due to the slope, which forces the air flow to uplift. In steep topographic slope areas, the vertical velocity of this forced updraft is considerable. Due to the topographic forcing, the airflow is lifted to saturation, and then the saturated air condenses into clouds, thus increasing precipitation. In weather prediction systems, there has been extensive research into topographic precipitation such as rainstorms and local extreme precipitation (Zhong et al. 2015; Chen et al. 2021). If there is a low-level jet in front of the slope, then a large amount of water vapor is transported to the mountain, and there are large water droplets in the convective unstable flow, which can increase the rainfall enhancement effect of the topography. Unfortunately, there is little research on sub-grid topography parameterization schemes in dynamic forcing uplift applied in ANM, especially in general circulation models (GCMs). Only Shen et al. (2007) proposed a parameterization scheme for the dynamic effects of sub-grid topography and its impacts on airflow and rainfall simulation, and applied this into the P–σ regional climate model of Nanjing University.

A simple and efficient parameterization scheme for sub-grid slope and aspect in global climate models needs to account for the effects of dynamic processes, especially the topographic dynamic uplift. However, due to the lack of detailed observation data of vertical velocity, radiation has a large amount of accurate observation data all over the world and the most significant sensibility to topography. Accordingly, we will use radiation for verification. In this study, we propose a sub-grid statistical description method based on Gaussian distribution and apply it to global areas with complex topography. In Sect. 2, we describe an improved method that utilizes the 1-km-resolution topography data and corresponding mathematical transformation of this method. In Sect. 3, we describe the global and regional applicability and reasonability of the method and introduce its application to the Tibetan Plateau. The deviation of surface solar radiation in mountains from unobstructed horizontal surfaces is analyzed in Sect. 4. Conclusions and discussion are given in Sect. 5.

2 Data and method

2.1 Data

The topography data used in this study is from the United States Geological Survey (USGS) DIGITAL ELEVATION MODEL (DEM) with a resolution of 1 km × 1 km. DEM data is derived from field surveys and airborne remote sensing techniques (Holland et al. 2006; Kornus et al. 2006). The observation data is the dataset of high-resolution (3 h, 10 km) global surface solar radiation (Tang et al. 2019) provided by the National Qinghai Tibet Plateau Science Data Center from 1989 to 2018. The result of the correctness verification is an offline calculation.

2.2 Statistical description

As stated in the introduction, slope and aspect have important effects on local radiation and topographic precipitation. Alpert and Shafir (1989) found that orographic rainfall on the small mesoscale is highly predictable with the adiabatic assumption that the uplift is determined by \(\mathrm{V}\cdot \nabla {Z}_{s}.\) In the study of topographic vertical velocity, through the P-coordinate system,

$${\omega }_{s}=\frac{d{p}_{s}}{dt}=\frac{\partial {p}_{s}}{\partial t}+{\overrightarrow{V}}_{s}\cdot \nabla {p}_{s}.$$
(1)

The equation of static equilibrium and deformation formula is as follows:

$$\begin{aligned} \vec{V}_s \cdot \nabla p_s & = - \rho g\vec{V}\cdot \nabla Z_s = - \rho g\left| {\vec{V}} \right|tan\theta_N cos\left( {\theta - \varphi_N } \right) \\ & = - \rho g\sqrt[2]{u^2 + v^2 }\left[ {\left( {tan\theta_N cos\varphi_N } \right)cos\theta + \left( {tan\theta_N sin\varphi_N } \right)sin\theta } \right] \\ \end{aligned}$$
(2)

where \(\uptheta\) is the wind direction, \({\theta }_{N}\) is the slope and \({\varphi }_{N}\) is the aspect.

In addition to dynamics, slope and aspect also play important roles in thermodynamics, that is, radiation. Kondratyev (1977) found that the dependence of |↓\({E}_{\mathrm{s},\mathrm{dir}}\)| on slope angle was formulated with the help of the geometry factor cos(α). Müller and Scherer (2005) found that a correction factor \({f}_{cor}\) is derived for \(\downarrow {E}_{\mathrm{s},\mathrm{dir}}\) combining the effects of slope and aspect, geometric enlargement of a sloping surface, and shadowing. Thus, according to Müller and Scherer (2005), the effective surface solar radiation downward (SSRD) \({\downarrow {E}_{\mathrm{s},\mathrm{dir}}}^{*}\) on an inclined surface can be computed as:

$${\downarrow {E}_{s,dir}}^{*}=\downarrow {E}_{s,dir}\left(\frac{1}{sin{\theta }_{S}}\right)\left(\frac{1}{cos{\theta }_{N}}\right)\left[cos{\theta }_{N}sin{\theta }_{S}+sin{\theta }_{N}cos{\theta }_{S}{cos}\left({\varphi }_{S}-{\varphi }_{N}\right)\right].$$
(3)

Here the ↓\({E}_{\mathrm{s},\mathrm{dir}}\) is the amount of solar radiation (shortwave radiation) reaching the unobstructed horizontal surface. The position of the sun is given by the sun elevation angle \({\theta }_{S}\) and sun azimuth angle \({\varphi }_{S}\). According to Eq. (3) it is possible to compute \({\downarrow {E}_{\mathrm{s},\mathrm{dir}}}^{*}\) by simply multiplying with a correction factor \({f}_{cor}\):

$$\downarrow E_{{\text{s}},{\text{dir}}}^{*} = \downarrow E_{{\text{s}},{\text{dir}}} \times f_{cor} ,$$
(4)

after simplification of Eq. (3),

$$\begin{aligned}{f}_{\text{cor }}&=1+\frac{tan{\theta }_{N}}{tan{\theta }_{S}}cos\left({\varphi }_{S}-{\varphi }_{N}\right)\\ &=1+\frac{(tan{\theta }_{N}cos{\varphi }_{N}){cos}{\varphi }_{S}+({tan}{\theta }_{N}s{in}{\varphi }_{N})sin{\varphi }_{S}}{{tan}{\theta }_{S}}\end{aligned}.$$
(5)

Noted that in this study we only consider the effect of sub-grid slope and aspect on solar direct radiation and do not consider the influence of scattered radiation and other additional radiation. This is because compared with SSRD, diffuse radiation accounts for less than 10%, which is negligible. The sky viewing angle factor is calculated by taking the average value.

It is shown that slope and aspect are used in calculations of topographic vertical velocity and radiation in Eqs. 2 and 5. However, after studying the two orographic factors of slope and aspect, we found that either the slope or aspect PDFs are irregular. After a large number of statistical analyses, it is found that while transforming them into trigonometric functions, tan(slope) × cos(aspect) (\({tan}{\theta }_{N}cos{\varphi }_{N}\) in Eqs. (2) and (5), call for TC, the same below) and tan(slope) × sin(aspect) (\(tan{\theta }_{N}{sin}{\varphi }_{N}\) in Eqs. (2) and (5), call for TS, the same below) conform to the Gaussian distribution. TC and TS are crucial features when calculating radiation and topographic vertical velocity.

In our study, we only consider the sub-grids with slopes > 5°. This is because in practical application, slope is an indicator that determines whether an area has complex topographic features. Plateaus with high altitudes and gentle slopes are defined as high-altitude plains rather than areas with complex topography. Moreover, the weather and climate characteristics of vast plains with very gentle slopes are difficult to distinguish. The dynamic and thermal effects of complex topography are easy to cover up in the regions with slope ≤ 5° and tan(slope) < 0.1. These unfavorable conditions for quantitative analyses and research are difficult to conduct in these regions. Therefore, only the grids with slopes > 5° are discussed.

2.3 Skewness and kurtosis test

The Gaussian distribution test is based on tests for skewness and kurtosis (Bai and Ng 2005). When sample n is limited and large (n ≥ 100), the standard deviation of sample skewness is as follows:

$$\mathrm{s}1= \sqrt{\frac{6(n-2)}{(n+1)(n+3)}} .$$
(6)

The standard deviation of sample kurtosis is as follows:

$$\mathrm{s}2=\sqrt{\frac{24n\left(n-2\right)\left(n-3\right)}{(n+1{)}^{2}(n+3)(n+5)}} .$$
(7)

The k-order central moment is as follows:

$${m}_{k}=\frac{1}{n}\sum_{i=1}^{n} {\left({x}_{i}-\overline{x }\right)}^{k} .$$
(8)

The coefficients of skewness and kurtosis are defined as follows:

$$\tau =\frac{{\mu }_{3}}{{\sigma }^{3}}=\frac{E\left[(x-\mu {)}^{3}\right]}{E{\left[(x-\mu {)}^{2}\right]}^\frac{3}{2}} ,$$
(9)
$$\kappa =\frac{{\mu }_{4}}{{\sigma }^{4}}=\frac{E\left[(x-\mu {)}^{4}\right]}{E{\left[(x-\mu {)}^{2}\right]}^{2}} .$$
(10)

According to the symmetry of the normal distribution, assuming the significance level is α, the final significance level should be \(\frac{\alpha }{2}\). If the skewness coefficient \({\tau }_{\alpha }=\mathrm{s}1\times \kappa -3\) and kurtosis coefficient \({\kappa }_{\alpha }=\mathrm{s}2\times \tau\) meet the normal distribution with significance level α, then the sample can be considered from a normal-distribution population.

2.4 Gaussian distribution

The normal Gaussian distribution uses the standard deviation (σ) or mean value to test, mainly to avoid the occurrence of small probability time. However, this approach lacks flexibility. In this study, this is transformed into a standard normal distribution (standardized transformation) before processing, and a representative value of several sub-grid topography values at the model grid scale is selected and can be easily described and applied.

An important factor of Gaussian distribution is that it is easy to obtain its percentile, which can be calculated according to the probability distribution value P of the normal distribution. The advantages are as followed. Assuming that a series of data × conforms to the Gaussian distribution, it is necessary to find the percentile representing 80% of the values (that is, the representative value of a series of data can represent 80% of the data characteristics). According to the probability P table, the value \({Z}_{p}\) corresponding to Z during Z-Scores (\(\mathrm{z}= \frac{x-\mu }{\sigma }\)), where \(\mu\) is the average and \(\sigma\) is the standard deviation) transformation, can be found. The representative value y is as follows:

$${y}_{p}=\mu +{Z}_{p}*\sigma ,$$
(11)

where \({Z}_{p}\) can be seen by table that the closest value when p = 0.8 is 0.7995, which corresponds to \({Z}_{p}\) = 0.84. In this way, the representative value (such as the average value) of the series of data can no longer be calculated according to the range of value range, which simplifies the problem.

3 Analysis the statistical description of global sub-grid topography

3.1 Statistical description of global sub-grid topography

According to the above method, the spatial distribution of the global sub-grid topography after trigonometric function conversion is shown in Fig. 1. It can be seen that the value of TC and TS do not change significantly with the increasing of horizontal resolution, especially in areas with complex topography around the world such as the Tibetan Plateau, Rocky Mountains and Andes Mountains It is indicated that the resolution has little impact on the characteristics of this complex orography areas, and the steepness of the topography and orientation of the mountains can not be covered up in high resolution.

Fig. 1
figure 1

Spatial distribution of global terrain after transformation of trigonometric function at different resolutions a TS and b TC at the 1° × 1° resolution, c, d As in a, b, but at the 0.5° × 0.5° resolution, e, f As in a, b, but at the 0.25° × 0.25° resolution

Full size image

Table 1 shows the proportions of the grids with the Gaussian distribution of TC and TS at different resolutions. The grids with the Gaussian distribution of TS at 1° resolution show the largest (98.3%) proportion, while those at a resolution of 0.1° show the smallest proportion of 95.92%. Those of TC at 0.5° resolution account for the largest proportion at 98.74%. Similar to TS, at 0.1° resolution, the proportion for TC is the lowest at approximately 96.22%. From 1° to 0.5° resolution, the total grid number increases by 1.3 × 104, while the grids that do not conform to the Gaussian distribution increase from approximately 100 to more than 200, accounting for approximately 1.5% of the total grids. At 0.1° resolution, only 1.5 × 105 out of 3.68 × 106 grids (accounting for ~ 5.1%) do not conform to the Gaussian distribution. Although the number of grids that do not conform to the Gaussian distribution increases, the increasing proportion is small. Table 1 shows that the slope and aspect of almost all topographies conform to the Gaussian distribution after the trigonometric function transformation. Although the proportion of grids conforming to the Gaussian distribution decreases slightly with the increasing of resolution, the reduction has little influence on the global scale (all grids in Table 1 refer to the grids with slope > 5. If the proportion of land in the unit grid is greater than 10%, then the grid is defined as land). The above analysis shows that with the increasing of resolution, the Gaussian distribution characteristics are stable, and the method is feasible. Therefore, we can confirm that the global sub-grid topography conforms to the Gaussian distribution after the statistical description. It is more obviously indicated from Fig. 2 that the Gaussian distribution conforms to the distribution characteristics of global complex topography. The spatial distributions are similar at different resolutions.

Table 1 The proportion of TC and TS conforming to Gaussian distribution at different resolutions
Full size table
Fig. 2
figure 2

The spatial distribution of global terrain conforming (red dot) and unconforming (blue dot) the Gaussian distribution at different resolution, a TC and e TS at the 1° × 1° resolution; b, f, As in a, e, but at the 0.5° × 0.5° resolution; c, g, As in a–e, but at the 0.25° × 0.25° resolution, d, h, As in a, e, but at the 0.1° × 0.1° resolution

Full size image

3.2 Characteristics of sub-grid statistical description in different regions

The spatial distributions of the proportion of TC and TS Gaussian distribution over the world are given above. Since different continents have different topographic and regional characteristics, the spatial distributions of elevation and topography in different continents are provided. The standard deviation of topography is commonly used in climate models such as CESM1.0.4 to represent the steepness of topography. The turbulent drag from sub-grid-scale topography is represented by the Turbulent Mountain Stress (TMS) scheme using USGS topographic data (Liang et al. 2017). Figure 3 shows the spatial distribution of the standard deviation of sub-grid topographic height and the grids conform to Gaussian distribution on the four continents. The average standard deviation of geography in Asia is the largest, reaching 1925 m while in Africa is the smallest at only 86.9 m. It is implied that the complex topography in Asia accounts for a large proportion, while the overall topography in Africa is relatively gentle. However, the standard deviation cannot represent the direction of the topography, and it is necessary to use slope and aspect for further analysis.

Fig. 3
figure 3

Maps of standard deviation of orography (shaded, units: m) and regional grids conforming to the Gaussian distribution (red dots), a TC and b TS over East Asia; c, d As in a, b but over North America; e, f As in a, b but over Africa; g, h As in a, b but over South America

Full size image

It can be seen that the Gaussian distribution in East Asia (Fig. 3a, b) accounts for the largest proportion. As the standard deviation of geography increases, the grids of the Gaussian distribution gradually increase, such as in Siberia and Southwest China, especially over the Tibetan Plateau. Except for the peaks of the mountains, the topographic statistical characteristics of most areas over the entire Tibetan Plateau conform to the Gaussian distribution. In addition, the Gaussian distribution in South Korea and North Korea accounts for a larger proportion, which may be related to the topography of Halla Mountain and the large mountain ranges in northeastern North Korea. However, the proportion of the Gaussian distribution in Africa (Fig. 3e, f) is very small. Although the altitude in Southeast Africa is higher than 1000 m and there is even an Ethiopian plateau called the “Roof of Africa,” we find that only a few grids conform to the Gaussian distribution, accounting for only approximately 5% of all land grids. Most other African regions do not conform to this feature because in those regions the topographic gradients are relatively gentle and less complex. The relationship between complex topography and altitude is not necessary or sufficient. Areas with complex topography may have high altitudes, but not all areas with high altitudes have complex orography. In North America and South America (Fig. 3c, d, g, h), the Rocky Mountains and Andes Mountains are the hotspots in previous studies. In these two regions, the topography basically conforms to the Gaussian distribution. However, on the Brazilian Plateau, the second largest plateau in the world, the land surface is relatively gentle, and the topography is less complex. Only a few grids on the east coast conform to the Gaussian distribution. In addition, there are some complex terrains near Mount Rey and Mount Roraima in Venezuela and the east coast of Canada that conform to the distribution characteristics of Gaussian statistics. Figure 3 shows the standard deviations of topography in Asia, North America, Africa and South America. The correlations between the Gaussian distribution and the complex-orography distribution in Asia, North America, Africa and South America are 91.1%, 90.6%, 85.6% and 89.9%, respectively. This result proves the reasonability of the sub-grid geostatistical description method. In addition, several mountainous regions with typical complex climate characteristics are selected for further analysis.

As shown in Fig. 4, most of the grids that conform to the Gaussian distribution are in mountainous and plateau areas. We selected the Tibetan Plateau, Rocky Mountains and Andes Mountains for detailed analyses. Figure 4a, b show the proportions of grids conforming to the Gaussian distribution over the entire world and three mountains for all land grids (slope > 0°). Worldwide and for the three typical mountainous landscapes, the proportion of grids with the Gaussian distribution at 1° resolution is the largest among the different resolutions. With the increasing of resolution, the proportions of grids with the Gaussian distribution of TC and TS decrease. This may be related to a decrease in effective grids and needs further analysis. Figure 4c, d show the proportion of the effective grids (slope > 5°) with Gaussian distribution over global and three mountainous landscapes. Table 1 clearly shows that the proportion of the global Gaussian distribution is close to 100%, and the proportion changes little as the resolution improves. For the three mountainous landscapes, in general, the Gaussian distribution of TC is slightly larger than that of TS, which is approximately 8% higher. Both TC and TS show the largest proportions of the Gaussian distribution at the resolution of 0.5°, followed by 0.25°, and that at 0.1° resolution is the smallest. From a regional perspective, the Gaussian distribution in the Rocky Mountains accounts for the largest proportion, followed by the Andes Mountains, and the Tibetan Plateau shows the smallest proportion, which is approximately 10% less than the Rocky Mountains on average. This indicates some differences in the topographic features of the above three mountain landscapes.

Fig. 4
figure 4

Proportion of the grids conforming the Gaussian distribution over the globe and three typical mountainous landscapes at different resolutions. a TC averaged over all land grids, b TS over all land grids, c TC over the complex terrain area, d TS over the complex terrain area

Full size image

4 Results of radiation deviation improvement

The above analysis gives the statistical characteristics of slope and aspect after mathematical transformations over three mountainous landscapes. We further take the radiation over the entire Tibetan Plateau area as an example to verify the reasonability of this method. A region of 25 °N–40 °N and 70 °E–105 °E was selected to cover the Tibetan Plateau. It has the most complex topography surrounded by the highest mountains in the world such as the Himalayas, Pamir, Kunlun Shan, Hengduan Mountains and Kailas Range (Chen et al. 2003; Wu et al. 2007). The slope is steep, and it is sensitive to radiation to radiation. This study uses the radiation scheme in CESM (Community Earth System Model) Version 1.2.1 developed by NCAR (National Center for Atmospheric Research). The SSRD of unobstructed horizontal surfaces (CTL), considering the sub-grid topography after radiometric correction according to Eqs. (5) and (6), are calculated offline (Rad) and compared with the observation data (OBS). The results are shown in Fig. 5. It indicates that the maximum SSRD of CTL in spring (Fig. 5a) is concentrated on the southwest side of the Tibetan Plateau and is approximately 280 W/m2. The minimum is located on the southeast side, with the spatial distribution characteristics indicating a decrease from southwest to southeast. The maximum SSRD of Rad and OBS (Fig. 5b, c) is located on the top mountains of the Tibetan Plateau at ~ 280 W/m2, with the spatial distribution characteristics showing a decrease from the top of the mountain to the surroundings. The SSRD is the largest in summer, which may be caused by the solar direct point near the Tropic of Cancer, the closest latitude distance to the Tibetan Plateau and the highest solar elevation angle. In this season, the maximum of the CTL SSRD (Fig. 5d) is located on the west side of the Tibetan Plateau at more than 350 W/m2, but the SSRD received by the slope at the south side of the Plateau is very small at only approximately 100 W/m2. At this time, the distribution of Rad and OBS SSRD (Fig. 5e, f) is still similar to that in spring, and the maximum is located over the Tibetan Plateau, decreasing from the center to all surroundings. Although the maximum SSRD received by CTL (Fig. 5g) in autumn is located on the Tibetan Plateau, which is similar to the results of OBS and Rad (Fig. 5h, i), the magnitude of SSRD is still much smaller, and its maximum is only about 205 W/m2, while the maximums of OBS and Rad are approximately 230 W/m2. Since the solar direct point in winter is in the southern hemisphere, the SSRD over the Tibetan Plateau is very small, and we will not analyze it here.

Fig. 5
figure 5

Spatial distribution of SSRD during boreal spring (top), summer (middle) and autumn (bottom) under applying the sub-grid slope and aspect scheme (left), unobstructed horizontal surface scheme (middle) and observation data (right), respectively (unit: W/m2). Black contour indicates elevation above 3000 m

Full size image

We further quantitatively analyze the radiation results of the Rad method. Table 2 gives the advantages of the Rad method when evaluated in terms of root mean square error (RMSE). It indicates that in spring, summer and winter, the RMSE of Rad is smaller than that of CTL. In the results of CTL, the RMSE is largest in autumn and smallest in summer; while in Rad, the RMSE is largest in spring and smallest in autumn. The most improved season of RMSE is autumn, with an improvement of approximately 18.2 W/m2. This improvement ratio reaches 38.4%, which may be due to the largest deviation between CTL and OBS in autumn that produces the best improvement. The improvement effect in spring is close to that in summer, at approximately 21%.

Table 2 The RMSE in the control and Rad test during different seasons
Full size table

Figure 6 shows the spatial distribution of the deviation between Rad, CTL and OBS. Overall, the deviation between the SSRD of CTL and OBS is concentrated on the south side of the Tibetan Plateau and windward slope (where the wind blows upward along the slope) in summer and autumn. The SSRD in spring (Fig. 5b) is higher than that in autumn (Fig. 5h), which leads to the difference between CTL and OBS in spring (Fig. 6a) is less than that in autumn (Fig. 6c). The same as in summer (Fig. 6b). Thus, it leads to a strong positive deviation between OBS and CTL on the windward slope of Tibetan Plateau. However, there is a significant improvement after considering the sub-grid slope and aspect. From the spatial distribution, it can be found that compared with OBS, the SSRD of CTL in autumn is generally lower (Fig. 6c), about 10–80 W/m2 lower than OBS. The largest deviations are located at the south and southeast foothills of the Tibetan Plateau, reaching 85 W/m2. After considering the slope and aspect (Fig. 6f), the deviation between Rad and OBS is significantly reduced, especially in the complex terrain area at the foot of the Tibetan Plateau. This result is consistent with the fact that the best improvement for RMSE occurs in autumn. The deviation of SSRD between CTL and OBS in summer (Fig. 6b) is about – 40 to 110 W/m2. The area with the largest difference is located at the southern foot of the Tibetan Plateau, reaching 105 W/m2. After considering the topography, the deviation is almost completely eliminated (Fig. 6e). In spring, except near the southern foothills of the Tibetan Plateau, the SSRD of CTL is significantly lower than that of OBS in most areas, with a deviation of about − 10 to 70 W/m2 (Fig. 6a). After considering the effect of sub-grid slope and aspect, the deviation between Rad and OBS is also reduced (Fig. 6d), but the improvement effect in spring is insignificant, and the reason needs to be further researched. A significant negative deviation occurs on the southern edge of Tibetan Plateau in spring after considering slope and aspect due primarily to smaller solar incident angles, which are blocked by surrounding mountains and/or shaded by the slope.

Fig. 6
figure 6

The spatial distribution of SSRD difference between the observation and control run (top, OBS-CTL) during a spring, b summer and, c autumn; d–f, As in a–c, but for the observation and Rad run (OBS-Rad, unit: W/m2)

Full size image

It can be seen from the above analysis that in spring, summer and autumn, the SSRD of CTL over the Tibetan Plateau is almost lower than that of OBS, and Rad effectively increases the SSRD in this region, but the specific areas need to be further analyzed. Figure 7 shows the spatial distribution of the difference between the SSRDs of Rad and CTL. The difference means \({\downarrow {E}_{\mathrm{s},\mathrm{dir}}}^{*}-\downarrow {E}_{\mathrm{s},\mathrm{dir}}\) and the improvement ratio is (\({\downarrow {E}_{\mathrm{s},\mathrm{dir}}}^{*}-\downarrow {E}_{\mathrm{s},\mathrm{dir}}\)) /\(\downarrow {E}_{\mathrm{s},\mathrm{dir}}\)). At heights above 3000 m, compared with CTL, the SSRD of Rad almost increased in spring, summer and autumn. The best improvement occurred in autumn, with a range of 130 W/m2 located at the southeastern edge of the Tibetan Plateau, with an improvement ratio of more than 130%. The entire windward slope improved by 80–100%, and the improvement on the Tibetan Plateau mountains is less at approximately 10–40%. The area with the best improvement effect in summer (Fig. 7b, e) is the south side of the plateau, with an improvement range of 125 W/m2 or approximately 105%, while the SSRD on the north side decreases by approximately 40 W/m2, accounting for ~ 10%. The improvement in spring (Fig. 7a, d) is concentrated on the east and southeast foothills of the Plateau, with an improvement range of 105 W/m2 and improvement ratio of ~ 80–100%. The results are consistent with those of previous studies (Lee et al. 2013, 2019). Significantly, more detailed topographic factors such as the sky view factor and surface albedo were discussed in previous studies, and only a specific time of a day was considered. In this study, we mainly consider the direct solar shortwave and do not consider factors such as shading, longwave and shortwave scattering. Thus, the reflected fluxes of radiation were improved more accurately by previous researchers.

Fig. 7
figure 7

The spatial distribution of SSRD difference between the Rad run and control run (Rad-CTL, units: W/m2) during a spring, b summer and, c autumn; d–f As in a–c, but for the relative difference ((Rad-CTL)/CTL)

Full size image

In order to focus on the improvement of the statistical method on the SSRD, we set up a control test using the grid-averaged topography (called the Mean test) to calculate the radiation, and an experiment using the statistical description method to analyze the difference in the maximum and regional average SSRD over the Tibetan Plateau with different slopes. Figure 8 shows the difference between the Mean and Rad tests. These results suggest that both the Mean test and the Rad test have improved SSRD. Regional maximum (Fig. 8a) and average (Fig. 8b) of SSRD over the Tibetan Plateau increase with increasing slope. Obviously, the results of the Rad test are better than those of the Mean test. The difference in the maximum SSRD between the Rad and Mean tests in spring is about 27–43 W/m2, and that of the regional average SSRD is about 25–30 W/m2. The divergence of the maximum SSRD in summer is about 27–57 W/m2, and that of the regional average SSRD is about 13–22 W/m2. The divergence of the maximum SSRD in autumn is about 19–34 W/m2, and that of the regional average SSRD is about 28–43 W/m2. The above analysis indicates that in summer, the divergence of the maximum SSRD is the largest, while it is the second largest in autumn and the lowest in spring. The difference in the regional average SSRD over the Tibetan Plateau between the Rad and Mean tests is the largest in autumn, while it is second in spring and the least in spring. This result is consistent with the previous conclusion that RMSE improved most in autumn, while the improvement effect in spring is not significant. Compared with spring and autumn, the SSRD over the Tibetan Plateau is largest in summer, so its maximum improved the most.

Fig. 8
figure 8

SSRD of a Maximum and b regional average with different slopes over Tibetan Plateau in spring (bule lines), summer (red lines) and autumn (green lines). Solid (dotted) lines indicate Rad (Mean) tests (units: W/m2)

Full size image

5 Conclusions and discussion

In this study, a statistical description method of sub-grid slope and aspect based on the trigonometric function is proposed for ANM over complex topography. We found that the PDF of this method conforms to the Gaussian distribution in global areas over complex topography and that the distribution feature is applicable to typical local mountainous landscapes. This method has been applied to calculate SSRD over the Tibetan Plateau. The results indicated that the statistical method performs better in simulating the maximum and regional average SSRD when compared to unobstructed horizontal surfaces and the grid-average method. The improvement is the most significant in autumn, and the divergence of the regional average SSRD ranges from 28 to 43 W/m2 over the Tibetan Plateau.

On the basis of the above results, we only analyzed the topography in terms of SSRD over the Tibetan Plateau. No research in other regions such as Africa. More research is needed to determine whether it is more reasonable to calculate sub-grid orography for each sub-grid because there are so few grids in Africa that follow the Gaussian distribution. On the other hand, the topography does not take into account the characteristics of the underlying surface, and we did not consider the albedo. Instead, we focused primarily on how sub-grid slope and aspect affect radiation variability and dynamic uplift. However, it is unclear why Africa lacks obviously terrain sensibility on radiation.

Besides, less attention has been paid to the sub-grid slope and aspect in numerical weather forecast models or climate models. It is unknown what kind of probability distribution function the sub-grid slope and aspect followed. Therefore, we provided a feasible method for numerical models from the perspective of statistics. The Gaussian distribution is applicable at the global scale and remains stable as the resolution increases. It can be observed that our method provided a new perspective on how sub-grid topography parameterization could be used in numerical models. Meanwhile, we analyzed the topographic features of regional landscapes and found that they still conform to the Gaussian distribution. It demonstrates the regional applicability of our methodology. There could be other more effective schemes for numerical models. Regarding our approach, it not only can be applied to both thermal processes and dynamic uplift, but it also makes calculations more effective.

Moreover, our method highlights how the statistical method can be applied to ANMs because it is independent of observational data and is applicable to both thermodynamics and dynamics. We established its viability in low-resolution radiation. However, the resolution of topography data is inadequate when compared to earlier statistical methods (Lee et al. 2013), and the statistical method only takes into account the effect of slope and aspect on SSRD. The difference between Rad and OBS is caused by the fact that this approach does not take into account how solar radiation is affected by variables like diffuse reflection, direct reflection, coupled reflection, and other aspects. Future studies based on Eqs. (4) and (5) can be applied to ANMs for testing based on this statistical method, and it will be taken into consideration whether variables like diffuse-reflected and diffuse flux can be amended similarly.

For GCMs in thermodynamics and dynamics, we presented a universal sub-grid topography statistical method, but it has not been applied to ANMs for numerical experiments. Only considering the sub-grid slope and aspect of SSRD over the Tibetan Plateau is verified offline in theory. In addition, sub-grid slope and aspect are of great significance in dynamics, especially for topographic vertical motion, such as in the vertical uplift of topography and the foehn effect. Few studies have taken slope and aspect into account in the sub-grid topographic parameterization scheme to simulate complex terrain precipitation, but improving precipitation simulation over complex terrain areas by the parameterization method has always been the main work for NWP (Fonseca et al. 2015). Previous studies only incorporated the parameterization scheme of slope and aspect in GCMs to represent surface radiation processes, which indirectly affected the ability of simulating precipitation (Gu et al. 2020). Note that we not only considered the thermodynamic and dynamic effects of the method in GCM but also further explored the effect of the method on the topographic precipitation simulation ability under the combined influence of thermodynamic and dynamic conditions in the future. In addition, further analyses will involve a wider test that introduces this method into AGCM considering the impact of adding topographic vertical motion on atmospheric circulation, temperature and precipitation.