A High-Precision Sub-Grid Parameterization Scheme for Clear-Sky Direct Solar Radiation in Complex Terrain—Part I: A High-Precision Fast Terrain Occlusion Algorithm

Abstract

In atmospheric modeling, sub-grid parameterization is an important method for studying the topographic effects of solar radiation using high-resolution Digital Elevation Model (DEM) data. For reducing the amount of computation, some approximate methods that can lead to errors are used in existing sub-grid parameterization schemes for clear-sky direct solar radiation (SPS-CSDSR). The lack of a high-precision fast terrain occlusion algorithm (HPFTOA) remains one of the biggest constraints in this field. This study proposed an HPFTOA. It mainly uses two kinds of acceleration algorithms. One method is to use a dynamic, lossless, and fast occlusion search radius. Another way is to use the rectangular grid for calculations within the accuracy of DEM data to avoid coordinate projection conversions. The test results indicate that the HPFTOA can carry out large-scale computation based on DEM data with a resolution of 90 m. Because it rarely uses approximation algorithms and considers the curvature of the Earth, SPS-CSDSR can achieve unprecedented precision. The HPFTOA can also be used in the fields of mountain solar energy assessment, remote sensing, and telemetry, including terrain-obscuring the probe. As computer performance improves and algorithms and execution code are optimized, the application prospects will be very broad.

1. Introduction

It is well known that solar radiation provides all of the Earth’s external heat. It warms the ocean, land, and air, driving weather and climate, and provides almost all of the energy for the ecosystem. Therefore, solar radiation is an important component in surface physical process models for weather, hydrology, climate, and ecological systems.

A mountain can be regarded as consisting of many inclined surfaces. As early as the 1970s, Swift and Kondratyev systematically presented the basic theory of solar radiation on inclined surfaces [1,2].

After computer technology and Digital Elevation Model (DEM) data began to be widely used, Dozier and Frew and other scholars successively proposed mountain solar radiation models based on DEM data [3,4,5,6,7,8]. These methods show that the spatial distribution of clear-sky solar irradiance (CSSI) can be strongly affected by terrain factors (altitude, slope, slope direction, terrain occlusion, etc.)

In order to introduce the terrain effects on solar radiation into the atmospheric model, sub-grid parameterization schemes based on high-resolution DEM data have been developed [9,10,11,12,13,14,15,16]. Based on mountain solar radiation models, these methods calculate the parameterized factors of sub-grid terrain effects on direct solar radiation, diffuse radiation, and reflected radiation. For example, when coupling the terrain effect factor $f_{Mp}$ for direct solar radiation, the downward direct solar radiation flux $E_d$ is instead by ${\displaystyle \frac{E_d}{cos{Z}_M}}{f}_{Mp}$ at the land surface. $E_d$ is calculated by the atmospheric model without terrain. $Z_M$ is the solar zenith angle. The coupling terrain effects on CSSI simulation experiments show that terrain effects on CSSI can alter the simulation results of the energy budget, surface temperature, and precipitation, such as in the Tibetan Plateau and related regions [10,12,15,17,18,19].

However, the above parameterization schemes all use approximate methods.

Table 1. The methods to reduce the amount of computation for SPS-CSDSR.

Index	Method	User
1	Ignoring the cast shadows.	Dubayah (1990) [4]; Zhang (2006) [10]; He (2019) [14]; Gu (2020) [12]
2	Assuming that sub-grid cells in an atmospheric model grid cell have the same slope and the slope direction is evenly distributed in all directions.	Dubayah (1990) [4]; Helbig and Löwe (2013) [14]
3	Using low-resolution DEM data.	Muller and Scherer (2005) [9] used the 30″ (~1 km); Zhang (2006) [10] used the 5′ (~9 km)
4	Using a function to predict the fraction of the surface in shadow in the atmospheric model grid.	Essery and Marks (2007) [9]
5	Ignoring that the slope area is larger than the horizontal projection area.	He (2019) [14]; Gu (2020) [12], Zhang (2022) [15]
6	Calculating terrain effects without the cast shadows on the atmospheric model grid cell by utilizing the solar position factors (altitude, azimuth) and mean value of sub-grid terrain factor (slope, aspect) function values.	He (2019) [14]
7	Based on the method similar to He, using the cast shadowless coverage factor $S{F}_p$ in the atmospheric model grid. $S{F}_p={\displaystyle \frac{1}{n}}{{\displaystyle \sum }}_{k=1}^{N}S{F}_k$. When a sub-grid cell is in cast shadow, $S{F}_k=0$; otherwise, $S{F}_k=1$. N is the number of sub-grid cells in the parent grid.	Zhang (2022) [15]
8	Adjusting $S{F}_p$ to $SF{C}_p$ by the equation $SF{C}_p=1-C_{ad}\left(1-SF{C}_m\right)$, where $C_{ad}=0.1849d{x}^{-1.443}+0.04561$. $dx$ is the east–west grid spacing (km) at this latitude.	Huang (2022) [16]
9	Using an insufficient occlusion search radius.	Huang (2022) [16] used the 27 km; Zhang (2022) [15] used the 9 km
10	Using DEM90 data, within a 27 km search radius and in 100 azimuth directions, the maximum terrain shielding angle was calculated for each sub-grid. These data were converted into a cast shadowless coverage lookup table based on 100 azimuth directions and 100 maximum terrain shielding angle sine levels.	Huang (2022) [16]

lists these approximations of sub-grid parameterization schemes for clear-sky direct solar radiation (SPS-CSDSR). They also parameterize isotropic diffuse radiation and reflection radiation and express the terrain’s shielding of diffuse radiation through sky visibility. In theory, these could all lead to varying degrees of bias. For example, when visibility is less than the distance to the mountain, the shielding of diffuse radiation may be excessive. Because reliable downward solar radiation observational data are lacking in mountainous areas [20], none of their accuracy has been widely verified.

This study will focus on solving the approximate problem about SPS-CSDSR.

The computation for SPS-CSDSR based on high-resolution DEM data (such as 90 m) is very costly. High-resolution DEM data are very dense and typically utilize isometric latitude and longitude coordinates to save storage space. When considering the curvature of the Earth to calculate cast shadows, the data often need to be transformed through projection methods (such as [9]), which can strain storage and memory resources. Another massive computation involves traversing the search terrain occlusion along the solar azimuth. When there is a lack of high-performance computers and efficient algorithms, such expensive computational costs are often unaffordable for a regional atmospheric model.

Early methods ignored the cast shadows and assumed that the terrain has specific distribution characteristics. These methods have facilitated the study of terrain effects on solar radiation, but they are not suitable for practical applications. Other methods used low-resolution DEM data. This does not fully reflect the terrain effect.

In recent years, a new approximation method was proposed by He et al. [14], utilizing the solar position on the atmospheric model grid to parameterize the clear-sky direct solar irradiance without cast shadows. This method was further developed by Zhang et al. [15] and Huang et al. [16]. But processing cast shadows in this way becomes difficult. He et al. [13] ignored cast shadows. Zhang et al. [15] used the cast shadowless coverage $S{F}_p$ in the mode grid; this resulted in negative deviation [16]. Huang et al. adjusted $S{F}_p$ to $SF{C}_p$ [16] (Index 8 in Table 1). However, according to its function, when $S{F}_p$ = 0.001, $SF{C}_p$ will be greater than 0.9 for an atmospheric model with a horizontal resolution of 3 km. This is almost equivalent to ignoring the cast shadows.

In addition, the lookup table (Index 10 in Table 1) for $S{F}_p$ loses some precision [16]. The accuracy of this lookup table is significantly lower than DEM90. Since the horizontal direction only has 100 directions, the horizontal resolution of 90 m can only be achieved within 1432 m. The absolute altitude error of SRTM DEM data is less than 10 m [21]. This accuracy can only be achieved within 2 km because the maximum terrain shielding angle sine is divided into 100 levels. As the distance increases, the accuracy of the lookup table will significantly decrease. For example, at 20 km, one azimuth will represent 1.3 km. The high accuracy of DEM data cannot be fully utilized.

Usually, direct solar radiation occupies the highest proportion during clear skies in the daytime, especially in mountainous areas with high atmospheric transparency. The result of SPS-CSDSR can be used as a parameterized factor for diffuse radiation and reflected radiation. Therefore, the deviation in SPS-CSDSR can reduce the estimation accuracy of total solar radiation, which can then impact the weather and climate prediction results of atmospheric models. Ensuring the accuracy of the calculation for direct solar radiation is crucial.

In summary, the lack of a high-precision fast terrain occlusion algorithm (HPFTOA) is one of the biggest constraints in this field. The purpose of Part I of this study is to optimize algorithms to avoid the use of approximate methods and make the computational cost feasible, thus achieving high precision in SPS-CSDSR.

Based on some key acceleration algorithms, this study presents the HPFTOA. Because it can have a wider range of applications (such as solar energy assessment and remote sensing correction) than SPS-CSDSR and requires to be shown in detail, it is presented as Part I of this study. Part II will focus more on the application of the high-precision SPS-CSDSR in the atmospheric model.

2. Methods

2.1. The Original High-Precision SPS-CSDSR

This paper defines the high-precision SPS-CSDSR needs to use DEM data with a horizontal resolution of at least 90 m. Equations (1)–(5) provided below are basically consistent with those of Müller and Scherer [9], except that they utilized DEM data with a horizontal resolution of 30 arc seconds (~1 km).

The power of direct solar radiation received in complex terrain can be measured by the clear-sky horizontal direct solar irradiance (CSHDSI), which represents the total direct solar radiation received per unit area of a horizontal plane per unit time. In atmospheric models, CSHDSI is equivalent to the direct downward shortwave radiation. This study found that some coupled simulation studies ignored the effect of geometric surface enlargement. It was emphasized by Müller and Scherer [9]. This would lead to an underestimation of total solar radiation. Therefore, in this paper, the acronym “CSHDSI” was used in the equations to prevent confusion and to help clarify this issue.

The CSHDSI in an atmospheric model grid cell can be calculated from sub-grid cells by the following equations.

$$CSHDS{I}_M={\displaystyle \frac{1}{N_s}}{{\displaystyle \sum }}_{k=1}^{N_s}CSHDS{I}_{sk}$$

(1)

with

$$CSHDS{I}_{sk}=DN{I}_{sk}·MAX\left(cos{I}_{sk},0.0\right)·S{F}_{sk}·A_{sk}/A_{sk0}$$

(2)

$$cos I_{sk}=cos {\alpha }_{sk}sin h_{sk}+sin {\alpha }_{sk}cos h_{sk}cos \left({\theta }_{sk}-{\beta }_{sk}\right)$$

(3)

where $M$ and $s$ indicate the atmospheric model grid cell and the sub-grid cell, respectively. $N_s$ and $k$ are the number of sub-grid cells within an atmospheric model grid cell and the sequence number of a sub-grid cell, respectively. $DNI$ is the Direct Normal Irradiance. $I_{sk}$ is the sunlight incidence angle between the plane normal and the solar beam. $S{F}_{sk}$ is the shading factor (0 for the cast shadow, 1 for the no cast shadow). $A_{sk}$ and $A_{sk0}$ are the surface area and the horizontal projected area of the sub-grid cell, respectively. Equation (3) is given by Kondratyev [2], where $I_{sk}$ depends on slope ${\alpha }_{sk}$, aspect (slope direction)${\beta }_{sk}$ (0 for north, increasing clockwise), solar altitude angle $h_{sk}$ between the solar beam and the horizontal plane, and solar azimuth angle ${\theta }_{sk}$ (0 for north, increasing clockwise).

SPS-CSDSR uses the basic assumption that a sub-grid cell and its parent grid cell have the same DNI, i.e., $DN{I}_{sk}\approx DN{I}_M$, thus establishing the following:

$$CSHDS{I}_{Mp}\approx DN{I}_M·f_{Mp}$$

(4)

with

$$f_{Mp}={\displaystyle \frac{1}{N_s}}{{\displaystyle \sum }}_{k=1}^{N_s}(MAX\left(cos{I}_{sk},0.0\right)·S{F}_{sk}·A_{sk}/A_{sk0})$$

(5)

where $p$ indicates parameterization. $DN{I}_M$ can be obtained from the radiation model of the atmospheric model.

Because more computation time is needed based on high-resolution DEM data, $f_{Mp}$ should be calculated offline before the atmospheric model running, such as in the work of Müller and Scherer [9].

The slope ${\alpha }_{sk}$ and aspect ${\beta }_{sk}$ of the sub-grid are calculated by the following equations [22].

$${\alpha }_{sk}=arctan\sqrt{{\left({\displaystyle \frac{\partial H}{\partial x}}\right)}^2+{\left({\displaystyle \frac{\partial H}{\partial y}}\right)}^2}$$

(6)

$${\beta }_{sk}=\{\begin{array}{l}{\displaystyle \frac{3\pi }{2}}-arctan\left({\displaystyle \raisebox{1ex}{$\frac{\partial H}{\partial y}$}\!\left/ \!\raisebox{-1ex}{$\frac{\partial H}{\partial x}$}\right.}\right),if{\displaystyle \frac{\partial H}{\partial x}}>0\\ {\displaystyle \frac{\pi }{2}}-arctan\left({\displaystyle \raisebox{1ex}{$\frac{\partial H}{\partial y}$}\!\left/ \!\raisebox{-1ex}{$\frac{\partial H}{\partial x}$}\right.}\right),if{\displaystyle \frac{\partial H}{\partial x}}<0 \\ 0,if {\displaystyle \frac{\partial H}{\partial x}}=0 {\displaystyle \frac{\partial H}{\partial y}}<0 \\ \pi ,if {\displaystyle \frac{\partial H}{\partial x}}=0 {\displaystyle \frac{\partial H}{\partial y}}>0 \\ undefined, if {\displaystyle \frac{\partial H}{\partial x}}=0 {\displaystyle \frac{\partial H}{\partial y}}=0\end{array}$$

(7)

${\displaystyle \frac{\partial H}{\partial x}}$ and ${\displaystyle \frac{\partial H}{\partial y}}$ are usually calculated by the second-order finite difference algorithm (2FDA) (Equation (8)) or the third-order finite difference algorithm (3FDA) (Equation (9)).

2FDA [3]:

$$\{\begin{array}{l}{\displaystyle \frac{\partial H}{\partial x}}=\left(H_{i+1,j}-H_{i-1,j}\right)/(2·D_x)\\ {\displaystyle \frac{\partial H}{\partial y}}=\left(H_{i,j+1}-H_{i,j-1}\right)/(2·D_y)\end{array}$$

(8)

3FDA [22]:

$$\{\begin{array}{l}{\displaystyle \frac{\partial H}{\partial x}}=\left(H_{i+1,j+1}-H_{i-1,j+1}+H_{i+1,j}-H_{i-1,j}+H_{i+1,j-1}-H_{i-1,j-1}\right)/(6·D_x)\\ {\displaystyle \frac{\partial H}{\partial y}}=\left(H_{i-1,j+1}-H_{i-1,j-1}+H_{i,j+1}-H_{i,j-1}+H_{i+1,j+1}-H_{i+1,j-1}\right)/(6·D_y)\end{array}$$

(9)

where $D_x$ and $D_y$ are the east–west and north–south spacing of the DEM data, respectively.

When using 2FDA and 3FDA for calculating CSHDSI, the equation $A_{sk}/A_{sk0}=1/cos {\alpha }_{sk}$ will be used together in Equations (2) and (5).

2.2. The Two Triangular Planes Algorithm for CSHDSI

The slope smooth effects [23] of 2FDA and 3FDA can cause errors in calculating CSHDSI in complex terrain [20]. Therefore, the two triangular planes algorithm (2TPA) is added in this paper for testing HPFTOA. Analyzing its algorithm shows that the error for $CSHDS{I}_M$ comes only from the occlusion factor, so 2TPA is perfect for the test terrain occlusion algorithm.

In Figure 1, (i, j) represent the grid coordinates of the DEM data. The two triangles are formed by connecting coordinate points. When $tan{\theta }_{sk}\ge 0$, the triangles in a sub-grid cell are positioned in the northwest and southeast (Figure 1 left). When $tan{\theta }_{sk}<0$, they are positioned in the northeast and southwest (Figure 1 right). This avoids internal occlusion as much as possible. Otherwise, excessive internal occlusion is difficult to deal with.

The variables will be calculated at the center (i’, j’) of the DEM grid cell. The local altitude of point (i’, j’) is the average of altitudes at four DEM points surrounding it. When calculating the occlusion, the intersections of DEM gridlines and a line with the solar azimuth from point (i’, j’) are searched. But the DEM gridlines closest to point (i’, j’), which are the edges of the rectangular box in Figure 1, are not included. Otherwise, there would be slightly more cast shadows.

Equations (2) and (3) are changed to Equations (10) and (11), respectively.

$$CSHDS{I}_{sk}^T=DN{I}_{sk}·S{F}_{sk}({{\displaystyle \sum }}_{g=1}^{2}MAX\left(cos{I}_{skg},0.0\right)·A_{skg})/A_{sk0}$$

(10)

with

$$cos I_{skg}=cos {\alpha }_{skg}sin h_{skg}+sin {\alpha }_{skg}cos h_{skg}cos \left({\theta }_{skg}-{\beta }_{skg}\right)$$

(11)

where $T$ indicates the triangulation method, $g$ is a triangular serial number, $I_{skg}$ is the angle between the normal of the triangular plane and the solar beam, and $A_{skg}$ is the triangle area calculated using Heron’s formula.

Equation (5) is changed to Equation (12).

$$f_{Mp}={\displaystyle \frac{1}{N_s}}{{\displaystyle \sum }}_{k=1}^{N_s}S{F}_{sk}\left({{\displaystyle \sum }}_{g=1}^{2}MAX\left(cos{I}_{skg},0.0\right)·A_{skg}\right)/A_{sk0}$$

(12)

${\displaystyle \frac{\partial H}{\partial x}}$ and ${\displaystyle \frac{\partial H}{\partial y}}$ are calculated by coordinates and altitudes at the three points of the triangle. For example, when $tan {\theta }_{sk}<0$, in triangle T₂,

$$\{\begin{array}{l}{\displaystyle \frac{\partial H}{\partial x}}=\left(H_{i+1,j-1}-H_{i,j-1}\right)/D_x\\ {\displaystyle \frac{\partial H}{\partial y}}=\left(H_{i,j}-H_{i,j-1}\right)/D_y \end{array}$$

(13)

2.3. The High-Precision and Fast Terrain Occlusion Algorithm (HPFTOA)

The basic principle of the HPFTOA to judge terrain occlusion based on DEM data is the same as the commonly used method (e.g., Li [24]). However, the HPFTOA takes into account the curvature of the Earth and the curvature of the latitude line. So, it can be fully adapted to a wide range of mountain environments. More importantly, it speeds up the computation in two main ways. One is to use a dynamic, lossless, and fast occlusion search radius. Another is based on a rectangular grid for reducing coordinate conversion when the curvature of the latitude line does not affect the accuracy. Except for the polar regions, it does not need to convert the projection of DEM data coordinates.

Given the critical importance of the terrain occlusion algorithm for SRS-CSDSR, the detailed method of the HPFTOA is provided here. It cannot be implemented in an interpreted programming language such as Python because of the huge amount of computation. Fortran 90 and an Intel compiler were used in this study.

There are five steps to determine whether point A is obscured in Figure 2.

2.3.1. Determining Whether Point A Is Obscured by the Lowest Horizon of the Data Area

The area of CSHDSI is defined as the study area. All the required data areas are defined as a data area, which is larger than the study area because the occlusion can come from outside the study area. The method for determining the size of a data area is given in Section 2.3.2.

Firstly, the lowest solar altitude angle $h_0$ in the study area will be calculated. When the solar altitude angle $h\le h_0$, point A will be obscured by the horizon with an altitude $H_{min0}$. $H_{min0}$ is the minimum altitude in the data area.

In Figure 2, according to geometry, the angle $h_0$ equals the angle AOT, then

$$h_0=-arccos\left((R_0+H_{min0}\right)/\left(R_0+H_A\right))$$

(14)

where $R_0$ is the Earth’s radius.

Because $H_A\ge H_{min0}$, $h_0\le$ 0. This means that the sunlight can be from below the horizon.

The study area can be divided into many small rectangular regions to quickly determine whether these regions are completely obscured by the horizon. A smaller search occlusion radius in the small data area can also be obtained.

When not considering the sunlight from below the horizon, $h_0$ can be directly set to 0.

2.3.2. Obtaining the Maximum Search Occlusion Radius in the Data Area

Geometric analysis can find the farthest distance, where a point with the highest altitude of the data area can occlude point A. This distance will be used as the maximum search occlusion radius in the data area ($SR{A}_{\max })$. In Figure 2, the points C₁, C₂, and C₃ are exactly those three positions of the highest point when the solar altitude angles are $h>0$, $h=0$, and $h<0$. Set $H_{max0}$ represents the maximum altitude in the data area, similar to the principle of analyzing radar terrain masking [25]; Equation (15) for $SR{A}_{max}$ is given as follows.

$$SR{A}_{max}=\{\begin{array}{l}min\left(R_1,R_2\right),if h\ge 0 \\ R_3+R_4,if h_0<h<0\end{array}$$

(15)

with

$$\{\begin{array}{l}{R}_1=\left(H_{max0}-H_A\right)/tan h \\ R_2=\sqrt{\left(H_{max0}-H_A\right)\left(2R_0+H_{max0}+H_A\right)} \\ R_3=\sqrt{\left(H_A-H_{min0}\right)(2R_0+H_A+H_{min0})} \\ R_4=\sqrt{\left(H_{max0}-H_{min0}\right)(2R_0+H_{max0}+H_{min0})}\end{array}$$

Obviously, $SR{A}_{max}$ does not lose accuracy.

When $h<0$, the width of the data area increases by $2R_4$ around the study area. When $h\ge 0$, the increased width of the data area is $R_4$.

2.3.3. Determining Whether Point A Is Obscured by the Horizon in Solar Azimuth

On the DEM coordinate plane, the grid coordinates (i, j; including decimals for point C) of point C are the coordinates of the intersection of the longitude–latitude gridlines and the solar azimuth line. From this, a set of altitudes $H_c$ are obtained within $SR{A}_{max}$. Following the same principles as in Section 2.3.1, it will be determined whether point A is obscured by a new horizon.

Follow the steps to obtain $H_c$.

(1) Using the rectangular grid coordinate within a distance with no loss of DEM data accuracy

$H_c$ is the linear interpolation of the altitudes of two DEM points near point C. Due to the curvature of the latitude line, there is an error in the coordinates of C. When the error just exceeds half of the DEM data spacing, the distance between the coordinates of A and C is defined as $L{C}_{grid}$.

In the actual calculation, a lookup table of $L{C}_{grid}$ at different latitudes and different azimuths was made by the following Equation (16) in advance. An approximate $L{C}_{grid}$ was found from this table based on the latitude and solar azimuth at point A.

Equation (16) is from Aviation Formulary V1.47 [26]. It uses the latitude $latA$ and longitude $lonA$ of point A, solar azimuth ${\theta }_{sk}$, and the distance $Lc$ (arc BD) to calculate the coordinates (latitude $latC$ and longitude $lonC$) of point C. To speed up, Equation (16) is simplified by d instead of $sin d$ based on $SR{A}_{\max }$ being less than 700 km for Mount Everest and sea level.

$$\{\begin{array}{l}d=Lc/R_0 \\ tc=2\pi -{\theta }_{sk} \\ latC=arcsin\left(sin\left(latA\right)cos\left(d\right)+cos\left(latA\right)dcos\left(tc\right)\right) \\ lonC=\{\begin{array}{l}lonA, ifcos\left(latA\right)=0 \\ mod(\left(lonA-arcsin\left(sin\left(tc\right)d/cos\left(latC\right)\right)+\pi ,2\pi \right)-\pi , if cos\left(latA\right)\ne 0\end{array}\end{array}$$

(16)

When the search distance $Lc$ $\le$ $L{C}_{grid}$, the latitude–longitude grid is treated as rectangular. This will speed up computation by avoiding the use of Formula (16).

The following method is used to reduce the calculation time: if $|sin {\theta }_{sk}|\ge 0.5$, select the longitude lines to calculate intersections; if $|sin {\theta }_{sk}|\le 0.8$, select the latitude lines.

(2) Considering the curvature of the latitude line

When $Lc$ > $L{C}_{grid}$, the latitude–longitude coordinates of point C are calculated using Equation (16) with the solar azimuth ${\theta }_{sk}$ and $Lc$. $Lc$ is determined by incrementally increasing the value of $D_{LC}$ on $L{C}_{grid}$. $D_{LC}$ can be 1–5 times as much as $D_y$ when using DEM90. This depends on the balance between the accuracy and speed of calculation.

At high latitudes, a smaller $L{C}_{grid}$ results in slower calculation speed. Because the DEM data become denser in the latitude direction, a $D_{LC}$ that increases based on latitude and azimuth can be applied to speed up.

2.3.4. Obtaining the Search Occlusion Radius in Solar Azimuth

Following the same principles as in Section 2.3.2, the search occlusion radius $SR{S}_{max}$ in solar azimuth can be obtained based on a set of altitudes $H_c$.

2.3.5. Determining Whether A is Obscured within ${\mathrm{SRS}}_{\max }$

In Figure 2, the points B and D are the sea level coordinate positions of point A and point C. The lines AG and BH are horizontal. The lines AB and CH are perpendicular to the lines AG and BH. Point F is the intersection of the solar ray and the line CH. λ is the angle AOC. $H_A$ and $H_c$ are the altitudes of points A and C, respectively. $h$ is the solar altitude angle.

$H_{GC}$ (line GC) is compared with $H_{GF}$ (line GF) to determine whether point A is obscured.

$$S{F}_{sk}=\{\begin{array}{l}0, if H_{GC}\ge H_{GF}\\ 1, if H_{GC}<H_{GF}\end{array}$$

(17)

with

$$\{\begin{array}{l}\lambda =L_{BD}/R_0 \\ L_{DE}=R_0/cos\lambda -R_0 \\ H_{GC}=\left(H_c-L_{DE}\right)·cos\lambda -H_A \\ H_{GF}=\left(R_{0}tan \lambda +\left(H_c-L_{DE})sin \lambda \right)\right)·tan h\end{array}$$

where $L_{BD}$ is the distance (arc BD) between the coordinates of point A and point C. The correction of the Earth’s curvature is considered by $L_{DE}$ and $cos\lambda$.

When $H_{GC}\ge H_{GF}$, point A is obscured and $S{F}_{sk}=0$. Otherwise, $S{F}_k=1$.

Equation (17) is also applicable when the solar altitude angle $h<0$. The analysis diagram is omitted.

2.4. DNI Mode to Be Used for Testing

Equation (18) from the solar radiation model r.sun will be used to simulate $DNI$. In the atmospheric model, it is provided by the radiation model.

$$DNI=G_{0}exp\left(-0.8662T_{LK}m\delta \left(m\right)\right)$$

(18)

with

$$\{\begin{array}{l}{G}_0=I_0\left(1.0+0.03344·cos\left(2\pi N_d/365.25-0.048869\right)\right) \\ m=exp \left(-H/8434.5\right)/\left(sinh+0.50572{\left(h+6.07995\right)}^{-1.6364}\right) \\ {\displaystyle \frac{1}{\delta \left(m\right)}}=\{\begin{array}{l}6.6296+1.7513m-0.1202m^2+0.0065m^3-0.00013m^4, if m\le 20\\ 10.4+0.718m, if m>20 \end{array} \end{array}$$

where $G_0$ is the extra-terrestrial irradiance, $I_0=1367 \mathrm{W}/{\mathrm{m}}^2$ (the solar constant),$N_d$ is the day number starting from 1 January of the year, $T_{LK}$ is the Linke turbidity factor for an air mass equal to 2, $H$ is the altitude, and $h$ is the solar altitude angle (in degrees).

r.sun can calculate direct, diffuse, reflective, and total solar irradiance [5]. It has been implemented in the GRASS^® Geographic Information System environment [27], and its applicability has been demonstrated by many works [28,29,30,31].

Equation (19) will be used to simulate $T_{LK}$ changes with the altitude $H$ on land surface, which borrows from Remund et al. [32]. $T_{LK0}$ is $T_{LK}$ at zero altitude.

$$T_{LK}=exp(ln\left(T_{LK0}\right)\left(1-H/\left(2·8435\right)\right)$$

(19)

3. Materials

3.1. Study Area

Taiwan Island and the eastern part of the Tibetan Plateau were chosen for testing because of the extremely complex topography there.

Taiwan Island is surrounded by seas, and the Central Mountain Range extends from the north to the south of the island. There are 100 peaks above 3 km, and the highest peak reaches 3952 m [33]. This region is ideal to verify the accuracy of the HPFTOA. To eliminate the interference, the altitudes of little islands and the continent surrounding Taiwan Island are set to 0.

The eastern part of the Qinghai–Tibet Plateau is at 84.68°–105.32° E and 24.68°–40.32° N. It includes Mount Everest, the highest peak on Earth, and the eastern part of the Himalayas. The terrain there is choppy.

3.2. Data

In this study, DEM90 is the Shuttle Radar Topography Mission (SRTM) 90 m Digital Elevation Database v4.1 [34] (obtained from https://developers.google.com/earth-engine/datasets/catalog/CGIAR_SRTM90_V4, accessed on 17 June 2023). DEM30 is SRTM DEM V003 ([35], obtained from https://search.earthdata.nasa.gov/, accessed on 22 January 2024).

In this paper, one horizontal resolution of the atmospheric model is about 3 km (Model-3 Km), with each grid cell containing 33 × 33 sub-grid cells of DEM90 data or 97 × 97 sub-grid cells of DEM data with a horizontal resolution of 1 arc second (~30 m) (DEM30). Another horizontal resolution of the atmospheric model is about 9 km (Model-9Km), with each grid cell containing 97 × 97 sub-grid cells of DEM90 data.

3.3. Software

The computer program in this study is written in Fortran 90 and compiled using the ifort compiler (version 2021.5.0) from Intel^® Fortran, which is available at no cost. This compiler is part of the Intel^® oneAPI Toolkits products for Linux. The latest no-cost version can be downloaded by connecting to https://www.intel.com/content/www/us/en/developer/articles/news/free-intel-software-developer-tools.html (accessed on 23 June 2023).

3.4. The Server for Testing

In this study, a very common server, Sugon I620-C30(6132), was used for testing. It was manufactured by Dawning Information Industry Co., Ltd. and purchased in Jinan, China in 2019. It is equipped with two Intel (R) Xeon (R) Gold 6132 CPUs (14 cores, @2.60 GHz, 19.25 MB L3 Cache, Launch Date of the third quarter of 2017), 128 GB of RAM, a Linux operating system, and the Intel (R) Fortran compiler version 2021.5.0. It is referred to as the server with 28 cores (SERVER-28C) in this paper.

4. The Test Plan

According to Equations (4) and (5), the error of $CSHDS{I}_{Mp}$ will mainly come from $DN{I}_M$ and $S{F}_{sk}$. $DN{I}_M$ will come from the atmospheric model when applied in practice, which is not included in this study. The focus of this study is to produce a fast and accurate terrain occlusion algorithm. Therefore, only the accuracy and computation time will be evaluated.

It is difficult to verify the HPFTOA by the observation of cast shadows or CSHDSI. In mountainous areas, the subjective judgment accuracy of using Landsat satellite data to identify cast shadows was only 86% [36]. Reliable downward solar radiation observational data are also lacking in mountainous areas [24]. Downward solar radiation products based on satellite remote sensing also exhibit significant errors in complex terrain [37,38]. The downward solar radiation includes direct, diffuse, and reflected radiation. There is no accurate model to separate them in complex terrain.

This study designs a method to verify the systematic deviation in the HPFTOA by assuming Taiwan Island is in a vacuum atmosphere. In this case, there is no atmospheric absorption or scattering, and DNI ≡ 1367 W/m². By the law of conservation of energy, when the solar altitude angle is greater than zero, the direct solar radiation on the land and the shadows on the sea is the same as when there is no terrain.

When the occlusion algorithm is not appropriate, the cast shadows will be systematically too much or too little.

The computation speed for SPS-CSDSR based on the HPFTOA and DEM90 on the eastern part of the Qinghai–Tibet Plateau will be tested.

5. Testing Results

5.1. Average Error in CSHDSI Based on HPFTOA on Taiwan Island

Figure 3 depicts a topographic map of Taiwan Island and the simulation CSHDSI in Model-3Km on the winter solstice in 2022. The simulation used 2TPA and took $T_{LK0}=1.5$ As can be seen from Figure 3a, the mountains of Taiwan Island run in a north–south direction. Significant differences in direct solar radiation can be observed between mountains and flat surfaces (Figure 3b–d).

Under the assumption of Taiwan Island in a vacuum atmosphere (DNI ≡ 1367 W/m²), using DEM90 and DEM30 data, the systematic bias of the HPFTOA in Model-3Km was verified hourly from 06:50 to 16:50 LST on the winter solstice in 2022. When using DEM30, $D_{LC}=3.0D_y$. When using DEM90, $D_{LC}=1.0D_y$.

The errors in the average CSHDSI based on 2TPA, 2FDA, and 3FDA are shown in Figure 4. The average CSHDSI is the mean value of all grid points within the land and shadows on the sea. The error is the difference between the average CSHDSI with and without terrain. For the perfect model, it should be zero.

Figure 4 shows that the 2TPA method leads to slight positive biases. The HPFTOA has no excessive systematic error based on DEM90 and DEM30 data. According to algorithms, slight biases only come from the occlusion factor. This may be caused by ignoring occlusion between two triangles of a sub-grid cell, despite setting the directions of the triangles to avoid them.

The difference between 2FDA and 3FDA in Figure 4 shows that the multi-order difference algorithm’s slope smoothing effect [23] can lead to systematic negative bias in CSHDSI. It can be mitigated by using a higher resolution of DEM30. According to Equations (8) and (9), the smoothing effect of 3FDA is stronger than that of 2FDA, so its negative bias of 3FDA is more significant. Of course, their partial negative bias may also be related to more occlusion of nearby sub-grid cell. Therefore, the 2TPA method is proposed instead of the multi-order finite difference method to calculate slope, followed by 2FDA.

5.2. Testing Computation Speed

The computation speed for the high-precision SPS-CSDSR based on the HPFTOA and DEM90 on the eastern part of the Qinghai–Tibet Plateau (Figure 5) will be tested.

The plans for testing are listed in

Table 2. The computation time for $f_{Mp}$ based on the HPFTOA and DEM90, and the maximum absolute error of the acceleration plans.

Plan		Computation Time (Minutes)		$\mathbf{MAXAE}\mathbf{in}\mathit{C}\mathit{S}\mathit{H}\mathit{D}\mathit{S}{\mathit{I}}_{\mathit{M}}$ (W/m2)
Index	Content	Summer Solstice	Winter Solstice	Summer Solstice	Winter Solstice
1	$D_{LC}=1.0D_y$	47.23	40.17	/	/
2	$D_{LC}=2.5D_y$	39.49	34.80	5.2	4.8
3	$D_{LC}=5.0D_y$	36.98	33.18	8.5	14.6
4	U27	40.72	24.24	181.6	357.4

. Where $D_{LC}$ is defined in Section 2.3.3 (2), U27 stands for the simultaneous use of three simplified methods: the lowest solar altitude angle $h_0=0$; ignoring the earth’s curvature; and using the fixed search occlusion radius of 27 km. U27 is the usual method, as shown in Index 10 of Table 1. All plans take $T_{LK0}=2.0$ and use 2TPA. The area of Figure 5 is divided into 10 × 8 small regions to speed up. The test uses 26 cores on SERVER-28C.

According to the plans in Table 2, the computing speed for an $f_{Mp}$ interval of 1 h during the daytime of the summer solstice (06:00~21:00 LST) and winter solstice (08:00~19:00 LST) in 2022 was tested. The computation times are listed in Table 2. The maximum absolute errors (the difference with plan 1, MAXAE) in $CSHDS{I}_M$ of plan 2, 3, and U27 were counted and are listed in Table 2.

The test results in Table 2 show that the computation costs for $f_{Mp}$ based on the HPFTOA and DEM90 are acceptable. The plan $D_{LC}=2.5D_y$ is a better scheme that takes into account both accuracy and cost. Its computation time is only about 15% more than the usual method, but its accuracy is guaranteed.

The computation time of plan U27 on the summer solstice is much longer than that on the winter solstice, because there are fewer occlusions around noon in this season, and the searching terrain occlusion must be traversed within the fixed radius of 27 km. Therefore, the acceleration effect of the dynamic lossless fast search radius is significant.

The performance of the CPU has been greatly improved over the SERVER-28C. Next the CPU of Gold 6526Y will be used for evaluation. Its overall performance is twice as good as that of SERVER-28C [39].

Table 3. The estimated wall clock times (in minutes) when using cores of type Gold 6526Y for $f_{Mp}$ a 30 min time interval during the daytime based on DEM90, 2TPA, and the HPFTOA.

Number of Cores	Small-Region		Medium-Region	Large-Region
	10° × 10°	30° × 20°	60° × 40°	80° × 60°
30	9.3~10.6	56.1~63.6	224.2~254.4	448.4~508.8
300	0.9~1.1	5.6~6.4	22.4~25.4	44.8~50.9
1000	0.3~0.3	1.7~1.9	6.7~7.6	13.5~15.3

presents the estimated wall clock time to compute $f_{Mp}$ during the daytime when using varying numbers of cores of type Gold 6526Y, the plan $D_{LC}=2.5D_y$, a 30 min time interval, and different sizes of study regions. The method of estimation is to convert the computation time in Table 2 according to the number of latitude–longitude grid points and the number of computer cores. Table 3 gives the shortest (winter solstice) and longest (summer solstice) wall clock time range during a day.

Table 3 shows that, for a several-day weather simulation, the large-region computation cost of SRS-CSDSR based on the HPFTOA and DEM90 is feasible for a research institution.

Based on the analysis above, DEM30 data are only suitable for a small study region when using the HPFTOA.

Figure 6 (located in the blue rectangle in Figure 5) shows the CSHDSI differences between the three simplified methods of plan U27 and the plan $D_{LC}=1.0D_y$ at 09:00 LST on the winter solstice in 2022. It can be seen that the simplified methods have large deviations in some small regions. The terrain in this region is extremely complex. Mount Everest (86°55′31″ E, 27°59′17″ N) is situated there. Ignoring the effect of Earth’s curvature would significantly underestimate the CSSI at that location. Therefore, the HPFTOA is of great interest for studying glaciers and extreme climates in towering mountains.

6. Conclusions

This paper reviewed the approximate methods used to reduce the computation for SPS-CSDSR and theoretically analyzes the obvious precision reduction they may cause.

Since these approximate methods are mainly used to simplify the processing of terrain occlusion, this study proposes the HPFTOA. It uses many acceleration algorithms on the premise of ensuring accuracy. One important method is to use a dynamic, lossless, and fast occlusion search radius. Another is that when the curvature of the latitude line does not affect accuracy, use the rectangular grid for calculation to avoid coordinate projection conversions. It takes into account the curvature of the Earth and the curvature of the latitude line, so it can be adapted to any mountain environment.

The verification by assuming Taiwan Island in a vacuum atmosphere shows that the HPFTOA has no excessive systematic errors based on DEM90 and DEM30 data.

The computation speed for SPS-CSDSR based on the HPFTOA and DEM90 was tested in the eastern part of the Qinghai–Tibet Plateau, where the terrain is very complex. The results indicate that the HPFTOA can carry out large-scale computation based on DEM90 and is feasible for the regional weather simulation during several days. Those approximation algorithms in Table 1 should no longer be used.

Reliable downward solar radiation observational data are lacking in mountainous areas [24]. So, it is difficult to improve SPS-CSDSR through observation verification. In this case, reducing the use of approximation methods is the main way to ensure the accuracy of SPS-CSDSR.

Because the use of approximation methods is avoided, SPS-CSDSR can be calculated with unprecedented accuracy. High-precision SPS-CSDSR can be built based on the HPFTOA. This will also indirectly benefit the parameterization accuracy of scattered and reflected radiation. From the significant differences in CSHDSI between mountainous and flat surfaces as shown in Figure 3, high-precision SPS-CSDSR will be important for studying the terrain effect in the atmospheric model with a high horizontal resolution of 1 to 3 km.

It is found that the smoothing effect of the multi-order finite difference slope algorithm can lead to negative deviation in direct solar irradiance. Therefore, the 2PTA method is recommended in this paper, followed by 2FDA.

In addition, according to the results of this study, the HPFTOA can have good computational efficiency at low and middle latitudes. At high latitudes, the computational efficiency may decrease significantly. Accuracy and cost need to be balanced based on validation and requirements. This paper has given the preliminary method. For polar regions, a coordinate projection transformation may be required.

Obviously, the application of the HPFTOA will not be limited to SPS-CSDSR. It can also be used in the fields of mountain solar energy assessment, remote sensing, and telemetry, including terrain-obscuring the probe. Because the solar altitude angle is not low when the polar orbit remote sensing satellite passes through, the occlusion search radius is relatively small, and it is does not occur many times during the day, so it will be possible to use higher-resolution data than DEM90.

Excellent algorithms, advanced computing platforms, and efficient code execution are the three factors that determine computing speed. All of this could change rapidly. Therefore, the application prospects of the HPFTOA proposed in this paper will be very broad and inspiring.

Part II of this study will further explore methods to decrease computational costs through day interpolation algorithms.