Assimilating remote sensing data into GIS-based all sky solar radiation modeling for mountain terrain

Remote Sensing of Environment 231 (2019) 111239

Contents lists available at ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Assimilating remote sensing data into GIS-based all sky solar radiation modeling for mountain terrain

T

⁎

Shuhua Zhanga, Xingong Lib, , Jiangfeng Shea, Xiaomin Penga a b

Jiangsu Provincial Key Laboratory of Geographic Information Science and Technology, School of Geography and Ocean Science, Nanjing University, Nanjing, China Department of Geography and Atmospheric Science, University of Kansas, Lawrence, KS 66045, United States of America

A R T I C LE I N FO

A B S T R A C T

Edited by Jing M. Chen

Solar radiation is the ultimate energy resource of earth surface energy balance and the main driving force of atmospheric, ecological and hydrological processes. Solar radiation over complex terrain has large spatial and temporal variation because of terrain shading and high cloud heterogeneity. While most existing GIS-based solar radiation models only work under clear sky condition, this research presents a solar radiation model which considers both terrain shading and anisotropic cloud attenuation and diffuse radiation using MODIS atmospheric products. Specifically, we use skyshed map, sky cloud maps, and sky weight map to represent angular distribution of sky obstruction, anisotropic cloud properties, and diffuse radiance over hemispherical sky, respectively. Combining skyshed map, sky weight map and sky cloud maps, we develop a solar radiation model where 3D geometrical relationships among sun, cloud, and terrain are considered and anisotropic diffuse radiance and cloud attenuation are modeled. Model results are evaluated using field observations in the Kunlun Mountains of western China. At Terra and Aqua overpass time, our model performs well with a mean relative bias (MRB) of −0.2%. It underestimates in clear and partly cloudy sky with a MRB of −5.86% and −4.79% and has a mean absolute relative bias (MARB) of 8.11% and 21.59% respectively. It overestimates under overcast sky with a MRB of 1.68% and has a MARB of 31.71%. Our model performs better when compared with existing instantaneous solar radiation products. For daily solar radiation, our model shows good performance with a MRB of 1.43% and MARB of 17.02%. Our model shows significant spatial variation of solar radiation within the study area with the influences from terrain and cloud. Our research provides a novel and improved approach to assimilating remote sensing data into GIS-based solar radiation modeling in mountainous terrain where observations are sparse and difficult to obtain.

Keywords: Terrain shadow Cloud attenuation Sky cloud map Solar radiation

1. Introduction Solar radiation powers many atmospheric, hydrological, and biological processes on earth surface (Tian et al., 2001; Marsh et al., 2012; Zhang et al., 2013; Chen et al., 2014). In mountainous terrain, spatial and temporal heterogeneity of solar radiation determines local variation and dynamics of eco-hydrological factors, such as air temperature, evapotranspiration, and snow/glacier melting, which are key inputs for ecological and hydrological models (Funk and Hoelzle, 1992; Li and Williams, 2008; Allen et al., 2006; Aguilar et al., 2010). Many different methods have been developed to model solar radiation over mountainous area. Those methods can be classified into three categories. The first and most accurate method is deterministic or stochastic 3D radiative transfer models (Helbig et al., 2009). Deterministic models, such as the spherical harmonics discrete ordinate method

⁎

(SHDOM), use an iterative process to compute the source function term of radiative transfer equations on a grid of points in space and the diffuse part of the source function is represented with a spherical harmonic expansion. Stochastic Monte Carlo (MC) method computes radiation by tracing stochastic trajectories of photons through the atmosphere and all the interactions of the photons with atmosphere and terrain surface in 3D space can be modeled (Chen et al., 2006; Mayer, 2009; Helbig et al., 2009). However, both deterministic and stochastic 3D radiative transfer models are computationally very expensive. The second category of methods uses remote sensing observations to retrieve solar radiation. Those methods include physical models (Miesch et al., 1999) and hybrid methods (Wang et al., 2018). Physical models have the same shortcoming of significant computational expense as 3D radiative transfer models do. Hybrid methods, as shown in the most recent study of Wang et al. (2018), first simulate solar

Corresponding author. E-mail address: [email protected] (X. Li).

https://doi.org/10.1016/j.rse.2019.111239 Received 22 October 2018; Received in revised form 26 May 2019; Accepted 1 June 2019 0034-4257/ © 2019 Elsevier Inc. All rights reserved.

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

radiation received on a flat surface and at the top of atmosphere (TOA) using a radiative transfer model with different atmospheric conditions. A simpler model, in Wang et al. (2018)’s case an artificial neural network (ANN), is then built to estimate solar radiation received on a flat surface using remote sensing atmospheric products and TOA radiance as inputs. The result is further corrected for topographic effects using cosine of incident angle, cast shadow and sky view factor. Their approach, however, is questionable because the trained ANN model expects the TOA from a flat surface but the remote sensing TOA is a result of various terrain effects. In addition, their approach only works for clear sky as it does not quantify cloud effects at all. The third category of solar radiation models is based on digital elevation models and considers terrain shading effect, and some of them are implemented in GIS software (Dozier and Frew, 1990; Dubayah and Rich, 1995; Kumar et al., 1997; Fu and Rich, 1999; Li et al., 1999; Corripio, 2003). These models can be generally categorized as either solar based or ground based (Liu et al., 2012; Fu and Rich, 1999; S. Zhang et al., 2015). Solar based models follow sun's position in a day and detect terrain shading for all the locations on a DEM at specific solar positions. Those models are limited to a certain geographical extent. Ground based models, on the other hand, calculate sun's track (sunmap) and visible sky (skyshed map) on the hemispherical sky at each location on a DEM. Shading from surrounding terrain in direct and diffuse radiation can be obtained by overlaying sunmap or sky map (discretized hemispherical sky map) with skyshed. Solar based models are more accurate for direct radiation calculation (S. Zhang et al., 2015) while ground based models have some advantages for modeling diffuse radiation where sky map and skyshed map can be used to represent radiance distribution and surrounding terrain shading effect over hemispherical sky, respectively (Li et al., 2016). DEM-based GIS solar radiation models are efficient and they can accurately model terrain shading effect, including self-shading, shading from surrounding terrain and sky view. However, atmosphere effects in these models are based on simple empirical attenuation equations and most of them are just for clear sky condition. Cloud has significant attenuation on solar radiation in mountainous terrain as it occurs frequently due to terrain uplift and topographical variation (Fitzpatrick et al., 2004; Fitzpatrick and Warren, 2005; Chen et al., 2012; Zhang et al., 2013; Šúri and Hofierka, 2004). Terrain shading and cloud attenuation are the two major challenges in modeling solar radiation in mountain terrain (Šúri and Hofierka, 2004). Remote sensing atmospheric products provide an opportunity to empower GIS-based solar radiation models to account for both terrain shading effect and cloud attenuation. There are some studies that have tried this approach using MODIS (MODerate resolution Imageine Spectroradiometer) and other remote sensing atmospheric products (Van Laake and Sanchez-Azofeifa, 2004, 2005; Durr and Zelenka, 2009; Roupioz et al., 2018; Houborg et al., 2007). The main limitations of these studies are that terrain shading, cloud attenuation and their coupling effects on direct and diffuse radiation over mountainous areas are not fully considered. For example, while Houborg et al. (2007) considered pixel-sun geometry, they only considered terrain selfshading effect and the isotropic diffuse radiation component in their model is not enough for cloudy sky, especially for partly cloudy sky because of significant anisotropic radiance distribution. In addition, most existing studies modeled diffuse radiation using sky view factor (Dozier and Frew, 1990), which assumes an isotropic cloud distribution and diffuse radiation. Wang et al. (2005), however, found that substantial error may occur if the anisotropic nature of diffuse radiation is not considered. In this research, we proposed a new modeling approach for GISbased solar radiation modeling in mountain terrain, which takes the advantage of remote sensing atmosphere and cloud products and considers both terrain effect and anisotropic diffuse radiation in the model. We used sky cloud maps and sky weight map to represent anisotropic cloud and radiance distribution in the hemispherical sky. By overlaying

skyshed map, sky weight map, and sky cloud maps, our model fully describes the geometry relationships among sun, cloud and terrain and represents anisotropic effects of cloud and sky radiance.

2. Remote sensing data and ground observations Our model uses the 30-m ASTER Global Digital Elevation Model (http://www.gscloud.cn) and MODIS remote sensing data (https:// modis.gsfc.nasa.gov/). Atmosphere attenuation under clear sky condition is modeled using aerosol optical thickness (AOD), water vapor, ozone optical depth, and surface pressure data from MODIS atmosphere products MOD04/MYD04, MOD05/MYD05 and MOD07/MYD07. Because AOD and water vapor are missing for cloudy pixels, we used corresponding MERRA-2 datasets to fill the gaps (Gelaro et al., 2017). We used cloud optical thickness (COT), cloud water path (CWP), cloud effective radius (CER), cloud phase optical properties (COP) and cloud top height (CTH) data from MODIS (Collection 6.1) cloud products MOD06/MYD06 to determine cloud transmittance and cloud shadow for cloudy sky. Surface albedo data is extracted from MCD43 product. Sensor scan time in MOD06/MYD06 are also used to determine sensor's overpass time. Sensor zenith, azimuth angle, cloud top height, and DEM are also used to obtain cloud position in 3D space. Spatial resolutions of the remote sensing datasets are different (Table 1). AOD and water vapor data have the resolution of 10 km and 1 km respectively. Both ozone optical depth and surface pressure have the resolution of 5 km. All the cloud properties datasets have the resolution of 1 km. Although the atmosphere products have relatively coarse resolutions, the atmosphere is relatively stable and the coarse resolutions are usually acceptable. The resolution of the albedo dataset is 500 m. The time period of the remote sensing datasets used span between June 21, 2012 and December 22, 2014 during which we have field observations. Because the terrain and atmosphere products have different resolutions, we resampled the DEM from 30 m into 1 km to match the resolution of the MODIS cloud products to correct cloud position. For modeling solar radiation, we used 30 m DEM to simulate terrain effects at the three stations. When calculating the solar radiation for the entire watershed, we used 1 km DEM because of computation expense of using 30 m DEM. Field solar radiation measurements are collected between 2012 and 2014 at three automatic weather stations (AWS) in the Tizinafu mountain watershed in the Kunlun Mountains in the arid region of western China (Fig. 1). The three stations were installed in the same year but in different month (XHX in June 2012 and KD and MMK in December 2012). Campbell CS 300 Pyranometer (accuracy ± 5%) was installed to measure solar radiation. Global solar radiation is recorded every 10 min at the stations. The stations have different terrain shadow and their elevation ranges from 1927 to 2972 m (Fig. 1). Table 1 List of remote sensing products and their spatial and temporal resolutions used in our model. Product

Dataset

Spatial resolution

Temporal resolution

MOD06/MYD06

COT, CWP, CER, COP, CTH AOD Water vapor Ozone optical depth, Surface pressure Surface albedo AOD, Surface pressure, Water vapor，Ozone optical depth

1 km

Instantaneous

10 km 1 km 5 km

Instantaneous Instantaneous Instantaneous

500 m 0.5° × 0.625°

Daily 3 h or 1 h

MOD04/MYD04 MOD05/MYD05 MOD07/MYD07 MCD43 MERRA-2

2

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

Fig. 1. The location and skyshed of three AWSs (MMK, KD, and XHX) in the Tizinafu mountain watershed of western China.

3. All sky instantaneous global solar radiation

seen from the satellite is off its true location above the ground (Fig. 3). To model cloud attenuation in direct and diffuse solar radiation, we must know cloud position in 3D space. Cloud position from sensor view projection (A in Fig. 3) needs to be corrected to its orthogonal projection (B in Fig. 3) based on cloud top height and sensor view geometry. Cloud top elevation in MODIS cloud top (geopotential) height (CTH) data (MOD06/MYD06) is relative to mean sea level. Cloud top height relative to terrain surface (CTHr) is calculated as

Our model calculates broadband solar radiation separately for clear and cloudy sky conditions. We use sky cloud maps to represent cloud property distributions over hemisphere sky. Those sky cloud maps are then combined with sky weight map and skyshed map in an anisotropic diffuse radiation module where anisotropic radiance distribution, cloud attenuation and terrain shading effects are considered. For clear sky, direct radiation and anisotropic diffuse radiation are calculated separately. For cloudy sky, which means there is cloud along the path of solar ray in our model, direct radiation is neglected and only anisotropic diffuse radiation is calculated. The atmosphere is assumed to have two layers: a clear layer and a cloudy layer. Iqbal's C model (Iqbal, 1983) and Xie et al. (2016)’s model are used to calculate clear and cloudy sky atmospheric transmittance (more details in the Appendix A.1 and A.2), respectively. Fig. 2 provides a flowchart showing the major processes in the model, which are discussed below.

CTHr = CTH − DEM

(1)

where CTH is MODIS cloud top height and DEM is the ASTER DEM resampled to 1 km to match with CTH's resolution. Cloud location under orthogonal projection (i.e., B in Fig. 3) can be calculated by

π ⎞ ⎛ ⎞ X = x + ⎜CTHr /tan ⎜⎛ − SZA⎟ sin(SAA) ⎟/ R ⎝2 ⎝ ⎠ ⎠

3.1. Cloud data correction

π ⎞ ⎛ ⎞ Y = y − ⎜CTHr /tan ⎛⎜ − SZA⎟ cos(SAA) ⎟/ R 2 ⎝ ⎝ ⎠ ⎠

Fig. 3 shows the parallax effect on cloud from MODIS cloud product. The oblique viewing angle causes the apparent position of the cloud as

where X and Y are cloud column and row number under orthogonal projection and x and y are cloud column and row number in sensor 3

(2)

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

Fig. 2. Model flowchart.

maps which is discussed next.

3.2. Skyshed and sky cloud map Terrain shading and cloud attenuation are the two most important factors in modeling solar radiation in mountainous terrain. Skyshed, i.e., projected visible hemisphere sky, which is called viewshed in ESRI© ArcGIS' Solar Analyst, can be used to represent shadow from surrounding terrain and hemispherical diffuse radiation received at a point on the terrain (Fig. 5a). A sky cloud map is similar to an upward-looking hemispherical (or fisheye) cloud photography at a location. The dimensions of the sky cloud map are made same as those of the skyshed for consistence and overlaying purposes. In the model, we used a dimension of 200*200 for skyshed map and sky cloud maps. Each grid cell center on a sky cloud map corresponds to a zenith and azimuth angle on the hemisphere. A ray can be formed using the azimuth and zenith angle at the center of a cell and it can be determined whether the ray intersects with any cloud cells on the position corrected CTH. In the process, we use a 3*3 window around the cell because of gaps in the corrected cloud map. If there is no cloud in the window, then the cell on the sky cloud map is clear (i.e., no cloud attenuation occurs in the cell). If there are one or more cloud cells in the window, the average cloud property of the intersecting cloud cells is calculated and saved on the sky cloud map. An example COT sky cloud map and its overlay with the skyshed map at the XHX station at Terra overpass time on July 29, 2012 are shown in Fig. 5b and c, respectively. With corrected cloud data, a sky cloud maps for each MODIS cloud property (i.e., CTH, COT, CWP, CER and COP) is generated for use in our model.

Fig. 3. Geometric relationships in cloud position correction. A in the figure is cloud position in sensor view projection and B is cloud position in orthogonal projection.

viewing projection respectively. SZA and SAA are sensor zenith and azimuth angle (relative to the south), respectively. R is the cell size of the raster data. Using the above equations, MODIS datasets CTH, COT, CWP and CER are corrected to orthogonal projection. Fig. 4 shows a comparison between the original and corrected CTH cloud map. The sensor azimuth angle for the original map is −1.48 rad (relative to the South and east is negative and west is positive) and sensor zenith angle is 0.95 rad. Comparing Fig. 4a and b, cloud moves toward south-east after position correction. Fig. 4b also shows many gap pixels within corrected cloud. While Lu (2016) and Chen et al. (2012) choose to fill these gaps before using them in solar radiation modeling, we handle the gaps in the process of constructing sky cloud 4

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

Fig. 4. The effect of cloud position correction. (a) Original cloud map; (b) Position corrected cloud map. The sensor azimuth angle for the original map is −1.48 rad (relative to the South and east is negative and west is positive) and sensor zenith angle is 0.95 rad.

3.3. Anisotropic diffuse radiation

I (θs , a3) =

Diffuse solar radiation comes from every direction in the hemisphere sky that is not blocked by surrounding terrain. In most cases, diffuse radiation is anisotropically distributed over the hemispheric sky as circumsolar radiance decreases rapidly with angular distance from solar disk, and because of near horizon brightening or darkening radiance due to air mass effect and a uniform background radiance for cloudless sky. In addition, non-uniform cloud distribution over hemisphere sky also contributes to anisotropic diffuse radiation. In our model, radiation from the clear layer is anisotropically redistributed and anisotropic cloud distribution and attenuation in the cloud layer are considered using sky weight map and sky cloud map. We divide the sky into a series of sectors (sky sectors hereafter) along a certain number of azimuth and zenith directions. The sky sectors, when projected onto the same grid used in the skyshed map, become a sky map. The sky map with different weight in every sky sector, become a sky weight map. Anisotropic diffuse radiance at the top of the cloudy layer is modeled, depending on sky condition, by distributing either the diffuse radiation (clear sky) or the total radiation (cloudy sky) at the bottom of the clear layer to the sky sectors through a weighting scheme as

Ii = Ia ∗ wi

[2θs sin θs − 0.02πsin (2θs )] ⎫ ⎬ ⎭

where Ia is, depending on sky condition, the diffuse radiation or total radiation at the bottom of the clear layer. Ii and wi are the diffuse radiation and weight of the ith sector, respectively. Sector weights are calculated based on Brunger and Hooper (1993)’s radiance distribution formulation (Eqs. (4)–(6)) where a sector's weight is the proportion of diffuse radiation from the sector relative to all the sectors.

(4)

ψ = acos (sin θsinθs cos (φ − φs ) + cos θcosθs )

(5)

(6)

To calculate the weight for a sky sector, we let Gd, diffuse radiation over a horizontal surface, be 1. L(θ, φ) is the radiance in the sky sector with zenith angle of θ and azimuth angle of φ and (θs, φs) is sun's position; ψ is the angular distance in radians between the direction of solar disk (θs, φs) and the sky sector at (θ, φ). a0, a1, a2 and a3 are the parameters that depend on cloud properties and we use four sets of the parameters based on the COT in the visible sky. Sky weight maps with four different COT ranges are shown in Fig. 6. The maps show more anisotropic weight distributions with smaller COTs but also a similar spatial pattern where weights are higher near sun's position but decline further away from the sun. The weights are nearly isotropic for an overcast sky with a COT of > 30. To consider anisotropic cloud attenuation on diffuse radiation in the cloudy layer, the sky weight map is then overlaid with sky cloud maps (Fig. 7b). Mean COT, CWP, COP, CER within each sky sector are calculated. Cloud gap, i.e., the fraction of a sky sector that is covered by clouds, is also calculated as the ratio of the number of COT pixels over the total number of pixels in a sky sector (Fig. 7d). We then further overlay the result with the skyshed map (Fig. 7c) to consider terrain shading effect. A sky sector or parts of a sky sector may be blocked by surrounding terrain. Sky gap, which is the visible fraction of a sky sector, is calculated for every sector. Diffuse radiance from a sky sector is then calculated as

(3)

a + a1 cos (θ) + a2 exp (−a3 ψ) ⎤ L (θ, φ) = Gd ⎡ 0 ⎢ ⎦ ⎣ π (a 0 + 2a1/3) + 2a2 I (θs , a3 ) ⎥

2(1 − exp (−a3 π )) ⎤ 1 + exp (−a3 π /3) ⎧ × π − ⎡1 − × ⎢ ⎨ πa3 (1 + exp (−a3 π /2)) ⎥ a32 + 4 ⎣ ⎦ ⎩

Isi = Ii ∗ cos (INAi ) ∗ (τic ∗ CGi + (1 − CGi )) ∗ SGi

(7)

where Ii is the diffuse radiance in the ith sector at clear layer (Eq. (3)). τic is the cloud transmittance in the visible part of the ith sector and it's calculated based on Xie et al. (2016) with more details in the Appendix A.2. INAi is the incidence angle at the centroid of the ith sky sector. SGi

Fig. 5. A skyshed map and a COT sky cloud map at the XHX station: (a) Skyshed map; (b) a COT sky cloud map at Terra overpass time on July 29, 2012; (c) overlay of the skyshed map and COT sky cloud map. 5

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

Fig. 6. Sky weight maps (with a unit of 1/100) for a solar position with solar azimuth angle of −1.08 rad and solar zenith angle of 0.3 rad under different cloud conditions. Note that left side in the maps is east and west side is right which is the convention used in the skyshed map.

the window is completely clear, it is a clear sky. Otherwise it is a cloudy sky. We use a small window instead of just the cell where the sun is located, primarily because of the concerns on the errors in CTH and the propagation of the errors in the subsequent processing on CTH. Fig. 8b shows the accuracy of sky condition determined using MODIS cloud data with different window size when compared with conditions determined from field observations. To determine sky condition from field observations, we examined 144 overpass solar radiation measurements in January, April, July and October and compared them with the maximum overpass solar radiation in a month. We decided that when overpass solar radiation is > 80% of the maximum overpass solar radiation in the month it is considered as a clear sky and when overpass solar radiation is between 20% and 60% of the maximum overpass solar radiation it is considered as cloudy sky. Based on Fig. 8b, the overall accuracy almost doesn't change with window size larger than 7*7. We finally selected 11*11 cells as the window size to determine sky condition as it gives the highest accuracy for our model. If the sky is clear, direct radiation is the main component and diffuse radiation is relatively small. Direct radiation at the bottom of the clear layer is further corrected by incidence angle and surrounding terrain shading, which is decided by comparing solar altitude angle and the horizon angle in solar direction, and is calculated as

and CGi are the sky gap and cloud gap of the ith sky sector respectively. Total diffuse radiation is the sum of the radiation from all the sky sectors, i.e., Isd = ∑iIsi. 3.4. Instantaneous solar radiation As mentioned before, we divide the atmosphere into two layer: a clear layer and a cloudy layer in our model. At the bottom of the clear layer, direct and diffuse radiation are calculated using Eqs. (8) and (9):

Icb = I0 ∗ τb

(8)

Icd = I0 ∗ τd

(9)

where I0 is extraterrestrial radiation, τb and τd are beam transmittance for direct and diffuse radiation, respectively, and are calculated based on Iqbal (1983) with details in the Appendix A.1. Total solar radiation at the bottom of the clear layer is

Ict = Icb + Icd

(10)

For the cloudy layer, we first determine the sky condition at MODIS overpass time. A small window (the yellow box in Fig. 8a), centered at the cell where the sun is located, is used to determine whether there is cloud along the path of the incoming sun ray based on sky cloud map. If

Fig. 7. Overlaying a sky weight map, a sky cloud map and a skyshed map for calculating diffuse radiation at the XHX station. (a) A sky weight map with 16 divisions in both azimuth zenith directions with a different weight in each sector to represent anisotropic diffuse radiance; (b) overlay of the sky weight map and a sky cloud map at Terra overpass time on July 29, 2012; (c) overlay of the skyshed map, sky weight map, and sky cloud map; (d) a cloud gap map based on COT. 6

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

Fig. 8. Determine sky condition. (a) A window around sun's position is used to detect sky condition; (b) accuracy of sky condition determined using MODIS cloud data with different window size when compared with field observations.

4. Model evaluations

Table 2 Model error statistics under all sky condition at all overpass time. Terra

XHX KD MMK All stations

MRB(%)

MARB(%)

Number of samples

MRB(%)

MARB(%)

Number of samples

3.79 −10.95 7.24 0.02

27.14 22.07 21.35 23.52

761 665 431 1857

−4.13 −2.95 5.87 −0.41

30.76 30.57 22.87 28.07

782 677 439 1898

Isb = Icb ∗ cos (θs ) ∗ sf

Model results are compared with field observations at both Terra and Aqua overpass time under different sky conditions characterized by a sky clearness index. Relative bias (RB) is calculated at each overpass time as

Aqua

RB =

Is 1 − ρa ρg

(13)

where IM is modeled solar radiation and Io is observed solar radiation at the stations. Mean relative bias (MRB) and mean absolute relative bias (MARB) are then calculated for all overpass time to estimate model's systematic and absolute error respectively.

(11)

4.1. Model performance under all sky condition

where Icb is the direct radiation at the bottom of the clear layer, θs is solar incidence angle. If cos(θs) is negative, the surface is self-shaded, otherwise it is not. sf is a surrounding shading factor and is 0 if it is under shadow and otherwise is 1. sf is calculated by comparing solar altitude angle and horizon angle which is calculated based on the 30-m DEM. Diffuse solar radiation is calculated using the anisotropic diffuse radiation module discussed in Section 3.3 where Ia in Eq. (3) is the diffuse radiation at the bottom of the clear layer calculated using Eq. (9). If it is cloudy, we assume direct radiation is negligible and solar radiation is mainly from diffuse radiation which is also calculated based on the anisotropic diffuse radiation model discussed in Section 3.3. The total solar radiation (i.e., Ict from Eq. (10)), instead of the diffuse radiation, at the bottom of the clear layer is anisotropically distributed over the hemisphere for modeling anisotropic cloud attenuation within the cloudy layer. Finally, considering the multiple reflectance between land surface and atmosphere or cloud, total solar radiation at the surface is calculated as

It =

IM − Io × 100 Io

Table 2 shows the MRBs and MARBs at all the stations from Jun. 21, 2012 to Dec. 22, 2014 under all sky conditions. The mean relative bias (MRB) is −0.20%, indicating a good model performance with a slight underestimation. A possible reason for the underestimation is that the reflected radiation from surrounding terrain is neglected in our model. The mean absolute relative bias (MARB) is 25.79%. The correlation coefficients between observed and modeled solar radiation are 0.79 and 0.75 (both with a P < 0.01) for Terra and Aqua, respectively (Fig. 9). Most of the data points fall around the 1:1 line. Data points in the red circles in Fig. 9 are the overpasses whose sky conditions are misjudged in our model. Table 2 also shows that the model underestimates at the XHX and KD stations but overestimates at the MMK station though it performs better at the MMK station than at the XHX and KD stations. Seasonal MRBs and MARBs are shown in Table 3. The model performs better at Terra overpass time than at Aqua overpass time in all seasons except for the winter. With a MRB of −12.84%, the model underestimates more in summer than any other seasons. In addition, the model performs better in spring and autumn than it does in summer and winter. The possible reason is there are more clear sky days in spring and autumn than in summer and winter (Bo et al., 2016).

(12)

where ρa is atmosphere spherical albedo for clear or cloudy sky and is calculated based on Qin et al. (2015) and ρg is surface albedo from MODIS. Is is downward solar radiation which is (Isb + Isd) for clear sky and only Isd for cloudy sky. For high albedo surfaces (e. g., dry new snow and ice), Eq. (12) is problematic as multiple reflectance between cloud and surface is spectral dependent.

4.2. Model performance under different sky conditions We use sky clearness index (CI), which is the ratio of observed and extraterrestrial solar radiation, to examine model performance under clear (CI ≥ 0.75), partly cloudy (0.35 < CI < 0.75) and overcast (CI < 0.35) sky conditions (Iqbal, 1983). The percent of clear sky, 7

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

Fig. 9. Scatter plots of modeled and observed instantaneous solar radiation at (a) Terra and (b) Aqua overpass time. Data points in the red circle are the overpasses whose sky conditions are misjudged in the model. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

that our model neglects direct radiation under cloudy sky. The model shows overestimation for overcast sky with a MRB of 1.27% and 2.09% at Terra and Aqua overpass time respectively, and has a relatively large MARB of 31.71%. In Table 4, we also see that the underestimation is higher under clear sky than in cloudy sky which may be caused by the neglecting of reflected radiation in our model.

Table 3 Seasonal error statistics under all sky conditions. Terra

Spring Summer Autumn Winter

Aqua

Average

MRB(%)

MARB(%)

MRB(%)

MARB(%)

MRB(%)

MARB(%)

−1.15 −12.44 −1.65 3.20

22.00 25.45 20.77 30.12

−0.33 −13.25 5.46 −10.05

27.95 31.44 27.15 27.60

−0.74 −12.84 1.90 −3.43

24.97 28.45 23.96 28.86

4.3. Comparison with other studies To provide a context to assess our model, we compared our model results with the products and accuracies reported in existing studies (Table 5). Lu (2016) uses a Monte Carlo 3D radiative transfer model with 3D cloud constructed from MODIS cloud data for cloudy sky and gives a similar overall MARB of 24.6% which is obtained at one station at a relatively flat terrain at Lasha in Tibet. We also compared our results with the global 5-km instantaneous incident shortwave radiation product from Breathing Earth System Simulator (BESS) (Ryu et al., 2018). Since the BESS dataset only considered elevation influence but not terrain shading effects, we corrected BESS with self-shading effect for comparison with our results. The BESS dataset shows overestimation and has a MARB of 35.84% and 36.55% at Terra and Aqua overpass time respectively. For partly cloudy sky and overcast sky, our model shows significant accuracy improvement over BESS with a MARB reduction of 7.32% (from 28.92% to 21.59%) and 18.29% (from 50% to 31.71%) at Terra and Aqua overpass time, respectively. The BESS dataset without terrain self-shading correction has a much larger MARB of 47.68% indicating the correction is necessary in mountainous area. We also compared our model results with downscaled CERES SYN and CMSAF METEOSAT East solar radiation products estimated from satellite observations. CERES SYN datasets have a relatively larger spatial resolution of 1° and CMSAF datasets have a spatial resolution of 5 km. We downscaled the solar radiation datasets using the method in Antonanzas-Torres et al. (2014, 2018) with 30 m DEM. Details of the downscaling method are described in Appendix A.4. The MARB of downscaled CERES SYN and CMSAF solar radiation products at Terra

Table 4 Model performance under different sky conditions. Note that for clear and overcast sky, we only used the days that the sky condition was determined accurately based on MODIS cloud data. Terra

Clear sky Partly cloudy sky Overcast

Aqua

MRB(%)

MARB(%)

MRB(%)

MARB(%)

−4.45 −4.41 1.27

6.31 19.59 30.45

−7.28 −5.17 2.09

9.90 23.60 32.97

partly cloudy sky and overcast sky is 15%, 60% and 25%, respectively. The percentage of cloudy sky, including partly cloudy and overcast sky, at Terra and Aqua overpass time is 84.93% and 85.56% respectively, which indicates the frequent occurrence of clouds in the study area and a slightly more cloud occurrence in the afternoon (Aqua overpass time) than in the morning (Terra overpass time). From Table 4, we can conclude that under clear sky condition, the model performs very well and slightly underestimates solar radiation with a MRB of −4.45% and −7.28% and a MARB of 6.31% and 9.90% at Terra and Aqua overpass time, respectively. Under partly cloudy sky the model underestimates solar radiation with a MRB of −4.41% and −5.17% at Terra and Aqua overpass time, respectively. A possible reason for the underestimation is

Table 5 Comparison of our model results with existing studies. Note that we take the error statistics directly from Lu (2016) and we obtained the BESS (Ryu et al., 2018) data products from the authors and calculated error statistics by ourselves. MARB (%)

Anisotropic diffusion

Isotropic diffusion

BESS

BESS without correction

Terra overpass Aqua overpass Overall

23.52 28.07 25.79

26.29 31.26 28.78

35.78 36.55 36.16

48.9 46.67 47.68

8

Lu (2016)

Downscaled CERES SYN

Downscaled CMSAF

24.6

38.89 45.55 42.22

36.82 42.8 39.81

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

Table 6 Error statistic for daily solar radiation.

XHX KD MMK Average

MRB(%)

MARB(%)

4.80 −11.03 10.51 1.43

16.57 17.67 16.80 17.02

Fig. 10. Scatter plots of measured and estimated daily solar radiation at (a) XHX; (b) KD; (c) MMK stations.

and Aqua overpass times are 42.22% and 39.81%, which are larger than our model (Table 5). We also downscaled the CERES SYN datasets using the method proposed by Bair et al. (2016) which shows larger errors than those using the downscaling method by Antonanzas-Torres et al. (2014, 2018). Hinkelman et al. (2015) evaluated 3-h CERES SYN dataset using averaged field measurements within the hours and had smaller errors than the 1-h CERES SYN datasets used in our evaluation. Overall, our model has significant improvement on accuracy when compared with existing studies under cloudy sky condition and downscaled satellite solar radiation products. Considering terrain effects, including self-shading and surrounding terrain shading, are necessary for solar radiation modeling over mountainous area. We also examined the effectiveness of modeling anisotropic diffuse radiation by comparing it with an isotropic radiation model which is described in detail in the Appendix A.3. When using an isotropic radiation model, our model gives an overall MARB of 28.78% which is 3% larger than that (25.79%) of using the anisotropic model. Table 5 also shows the MARBs at Terra and Aqua overpass time using isotropic and anisotropic diffuse radiation model. Compared with the isotropic model, the anisotropic model performs better at Aqua overpass time than at Terra overpass time as it has more cloud in the afternoon than in the morning in the Kunlun Mountains (Liu et al., 2002). Overall, the anisotropic model performs better than the isotropic model. 4.4. Daily solar radiation model Many methods have been developed to upscale modeled instantaneous solar radiation to daily radiation. We use adjusted sinusoidal interpolation (SIN) to derive daily solar radiation based on instantaneous solar radiation estimated at Terra and Aqua overpass time (Wang et al., 2010). Given the instantaneous solar radiation I(T0) at overpass time T0, the instantaneous solar radiation I(t) at time t can be interpolated as:

sin I (t ) = I (T0 )

sin

(t − Tr ) π Ts − Tr (T0 − Tr ) π Ts − Tr

(14)

where Tr and Ts are the time of sunrise and sunset. With the instantaneous solar radiation estimated at Terra and Aqua overpass time, we have two functions, ITerra(t) and IAqua(t), based on the above equation. We assume that the relationship between cloud and surface elevation remains the same within the interpolation time period. We use ITerra(t) and IAqua(t) to integrate solar radiation from sunrise to the time of Terra overpass and from the time of Aqua overpass to sunset, respectively. Between the overpass time of Terra and Aqua, we calculate weighted mean of two sinusoidally interpolated solar radiation:

IT − A (t ) =

TAqua − t TAqua − TTerra

ITerra (t ) +

t − TTerra IAqua (t ) TAqua − TTerra

(15)

The overall function to calculate daily-integrated solar radiation is:

Isin =

∫T

TTerra

r

ITerra (t ) dt +

∫T

TAqua

Terra

IT − A (t ) dt +

∫T

Ts

Aqua

IAqua (t ) dt

(16)

As shown in Section 4.2, there are many partly cloudy sky in our study area and cloud may change quickly. As such, we calculated clear sky daily solar radiation using the same atmosphere attenuation model discussed before and modeled terrain effects including both self9

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

Fig. 11. (a) Instantaneous solar radiation in the Tizinafu watershed at Terra overpass on June 23, 2012; (b) Original cloud top height at the overpass time; (c) Corrected cloud top height at the overpass time; (d) hillshade map for interpreting (e); (e) difference in the annual sum of solar radiation in 2013 with and without considering topographic effect; (f) instantaneous solar radiation for clear sky. The three horizontal zones in the map are due to missing AOD data in MODIS filled using the data from MERRA-2.

4.5. Solar radiation in space

shading and surrounding terrain shading with a time step of 20 min. We then used clear sky radiation and cloud fraction in the visible sky to adjust the interpolated daily solar radiation as

IDaily = IClear ∗ (1 − cf ) + cf ∗ Isin

An instantaneous solar radiation map for the Tizinafu watershed at Terra overpass is shown in Fig. 11a. Solar radiation was modeled at a resolution of 1 km. Terrain self-shading was also calculated using the slope and aspect derived from the DEM with a resolution of 1 km. Cloud covers most of the watershed at the overpass time (Fig. 11b). Here we use CTH to show cloud distribution in the watershed as all other cloud datasets are corrected based on CTH. Comparing the solar radiation map (Fig. 11a) with the original cloud map (Fig. 11b) in the watershed, although their spatial patterns look similar, they are not always in consistence. For example, while there is few clouds for the part of the watershed in the dashed red box on the original cloud map (Fig. 11b), solar radiation in the same area is patchy (Fig. 11a). Comparing the cloud in the dashed red boxes in Fig. 11b and c, cloud moves toward southeast after position correction because the sensor looks at the ground from south-east with sensor azimuth angle of −1.48 rad (relative to south and east is negative) and sensor zenith

(17)

where IClear is clear sky daily solar radiation and cf is the maximum cloud fraction which is calculated as the ratio of the number of COT pixels and the number of all visible pixels at both Terra and Aqua overpass time. Table 6 shows the error statistics of daily radiation at the three stations. It shows an overestimation with a MRB of 1.43% and has a MARB of 17.02% indicating a good performance. While the model shows overestimation at the XHX and MMK stations it underestimates at the KD station. Fig. 10 shows the scatter plots and the correlation coefficients between estimated and measured daily solar radiation which are all > 0.9 at the stations.

10

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

skyshed map for shading determination, DEM data errors can still propagate into the calculation of horizon angle (S. Zhang et al., 2015 ). This is especially the case when horizon angle is very close to solar altitude angle where error in horizon angle can lead to large error in solar radiation modeling. As an example, at one Terra overpass time, the solar altitude angle (0.481 rad) is very close to the horizon angle (0.478 rad) and our modeled solar radiation (520 W/m^2), which includes direct solar radiation, is several times larger than the observed (37 W/m^2), which is primarily diffuse radiation, at that moment. For diffuse radiation, skyshed and sky weight map determine how much diffuse radiation is available. Sky weight map is discretized as sky sectors for diffuse radiation modeling. The resolution of the sky sectors is based on the number of zenith and azimuth divisions. Li et al. (2016) indicate that more divisions lead to more accurate sky view factor for diffuse radiation under clear sky though the factor is more sensitive to zenith divisions. We calculated MRB and MARB with different sky divisions at Terra overpass time under overcast sky condition as solar radiation is nearly uniform and is more sensitive to sector skygap and incidence angle. We find that the MARB almost doesn't change with > 16*16 sky sectors, but MRB tends to decrease and leads to negative with more sky divisions (Fig. 13). We suggest that more sky divisions should be used if there are more overcast days. We used a 16*16-division on sky sectors in our model as more divisions significantly increase computing time. MODIS atmosphere products are sensitive to the presence of clouds in the field of view, especially for ozone optical depth and AOD (Liu et al., 2015; Borbas et al., 2011; Sayer et al., 2014). Those are potential error sources for modeling solar radiation under cloudy sky condition. We used MERRA-2 datasets to fill the gaps when MODIS atmosphere datasets are missing though MERRA-2 datasets have a coarse resolution. In addition, MODIS cloud product CTH is used to estimate cloud position in 3D space, which further influences sky cloud map generation. CTH is retrieved based on cloud top pressure and gridded meteorological products from climate models that provide temperature profiles. Although improved algorithms are used to retrieve cloud top pressure in MODIS Collection 6.1, gridded meteorological products from the Global Data Assimilation Model (GDAS) model are often unreliable and have a spatial resolution that is too coarse for regions of highly variable topography (Menzel et al., 2015). In addition, because MODIS cannot detect multiple cloud layers in a vertical column, the process of correcting cloud products to the orthogonal projection introduce errors in cloud properties' map. After corrections, the cloud properties' map may have gaps (Fig. 4) as clouds seldom exist only at a single layer with similar height (Lu, 2016). The MODIS cloud position and data errors can propagate into sky cloud maps which are used to detect sky condition and calculate diffuse radiation. As an evidence, we found that only 83% of all sky conditions

angle of 0.95 rad. Although there are many gap pixels in Fig. 11c, the solar radiation map (Fig. 11a) rarely shows any gaps, indicating the gap filling method in the process of generating sky cloud maps works well. Comparing the solar radiation and cloud distribution maps in the dashed red boxes in Fig. 11a and c, cloud influence on solar radiation is now more consistent in space. There is also an area of high solar radiation along the southwest edge of the watershed in the dashed red oval in Fig. 11a. Although the area is cloudy on the original cloud map in Fig. 11b, we see that a linear area of clear sky outside the watershed in the dashed red oval in Fig. 11b is moved into the watershed after correction in the dashed red oval in Fig. 11c. All those indicate that considering the geometry relationships between sun, cloud and terrain and correcting cloud position are important for cloudy sky, especially for partly cloudy sky. We also calculated and compared the annual sum of solar radiation from our model which considers topographic effect and a model without considering topographic effect over Tizinafu watershed in 2013. Fig. 11e shows the difference (without topographic effect minus with topographic effect) between the annual sum of solar radiation. Most locations in the watershed have positive difference except for the areas of high altitude close to the southern edge of the watershed where snow and cloud exists all year round. The possible reason for the negative difference in Fig. 11e is how cloud information is used in models. In our model, cloud in solar path is used for direct radiation and cloud over hemispherical sky is used for diffuse radiation. In the model without considering topographic effect, only the cloud pixels right above at the stations are used. Referring to the hillshade map Fig. 11d, we can see that large positive difference in Fig. 11e is mainly located in the areas with complex terrain while small positive difference is in the areas of relatively flat terrain, especially in the area close to watershed outlet. Because of a cloudy sky on June 23, 2012, terrain effect is not significant at the overpass time. So, we calculated the solar radiation on a clear sky day for the entire watershed (Fig. 11f) with 30 m DEM, where the terrain effect is obvious and significant. 5. Discussions In mountainous terrain, terrain shading is an important factor that affects both direct and diffuse solar radiation. Horizon angle derived from a DEM is the key for detecting whether there is direct solar radiation which is the major component in clear sky and can change abruptly because of terrain shading (Fig. 12). Although using horizon angle directly derived from a DEM is more accurate than using a

Fig. 12. Clear and cloudy sky solar radiation observed at the XHX station on two days in January of 2013. Because of terrain shading, solar radiation changes abruptly on clear sky day.

Fig. 13. MRB and MARB with different sky divisions under overcast sky. 11

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

may also be caused by neglecting the reflected radiation from surrounding terrain in the model. Our model indicates that position correction on MODIS cloud data and considering the geometric relationships between sun, cloud and terrain are necessary for modeling solar radiation in mountain terrain under cloudy sky, especially partly cloudy sky. The relatively large errors under cloudy sky condition may be caused by the errors and uncertainty in MODIS cloud datasets used in our model. We also note that the ground stations are located in extreme mountainous terrain in the northwestern Tibetan Plateau with elevation ranges between 1575 and 6234 m and experience strong terrain effects (Fig. 1). The retrieval of cloud properties in such a terrain is more complex and less accurate, which lead to the difficulty of modeling solar radiation under cloudy sky condition in the area. In the future, we would like to evaluate other remote sensing cloud products, especially high temporal resolution remote sensing cloud products, such as Pathfinder Atmospheres Extended (PATMOS-x) data or Himawari-8 data (Yu et al. 2019) with more overpass times on a day, to improve daily solar radiation estimation. In addition to considering physical-process based narrow-band radiative transfer models for modeling atmospheric attenuation, we also would like to add more measurements at the stations to verify atmosphere products and various solar radiation components, and to include the modeling of reflected solar radiation. In addition, sky cloud map generation is time consuming process in our model as it runs pixel by pixel for each overpass time. Accelerated computing, such parallel computing, is necessary to speed up the process.

can be correctly determined using MODIS cloud data with a window size of 11 × 11 cells (Fig. 8b). For the overpass time with correct sky condition determination from MODIS cloud data our model shows a smaller MARB of 19.97% (compare with 23.52% in Table 2) and 24.58% (compare with 28.07% in Table 2) for Terra and Aqua respectively. 6. Conclusions In this research, we present an approach to model all sky solar radiation for mountainous terrain using primarily MODIS atmospheric products. We used sky cloud maps and sky weight maps to represent anisotropic distribution of cloud properties, diffuse radiance in hemispherical sky respectively. The geometric relationships between sun, cloud and terrain are obtained by overlaying skyshed map, sky weight map, and sky cloud map and both terrain and cloud effects on direct and diffuse radiation are fully modeled. Our model is evaluated with field observations at three stations with different terrain shading effects. Under clear sky condition, our model performs very well and only slightly underestimates solar radiation with a MRB of −5.86% and has a MARB of 8.1%. Under cloudy sky, the model has a MRB of −4.79% and 1.68%, and a MARB of 21.59% and 31.71% for partly cloudy sky and overcast sky, respectively. The all sky MRB of −0.2% shows a small model bias. The all sky MARB of 25.79% is relatively high as about 85% of the overpass time is under cloudy sky condition in the study area. Compared with BESS shortwave radiation products, our model reduced MARB by 10.37% for all sky. The most significant improvement over BESS is the 7.32% and 18.29% reduction in MARB under partly cloudy sky and overcast sky respectively, indicating our model provides a better approach for solar radiation modeling in mountain terrain when compared with existing studies. Daily solar radiation also performs well with a MARB of 17.02%. Instantaneous solar radiation map is also generated for the Tizinafu watershed. The map is spatially consistent with corrected cloud maps in the watershed. Terrain effect is more significant for clear sky which

Acknowledgements This work is supported by the National Natural Science Foundation of China, Grant/Award Number: 41701442, Project funded by China Postdoctoral Science Foundation (2018M632281) and CAS Strategic Priority Research Program (XDA20020100).

Appendix A A.1. Calculating clear sky transmittance (τb and τd) We compared several different broadband radiative transfer models and decided to combine Iqbal's C model (Iqbal, 1983) and the aerosol attenuation model in Yang et al. (2001) to calculate clear sky transmittance. Below are the various transmittances used in our model:

τr = exp (−0.0903ma0.84 (1 + ma − ma1.01)) τo = 1 − (0.1611U3 (1 + 139.48U3

)−0.3035

(1)

− 0.002715U3 (1 + 0.044U3 + 0.0003U3

τg = exp (−0.0127ma0.26)

(2) (3)

τw = 1 − 2.4959U1 ((1 + 79.034U1

)0.6828

τa = e (−mrβ (0.6777 + 0.1464mrβ − 0.00626(mrβ τaa = 1 − 0.1(1 − ma +

2 )−1)

ma1.06)(1

+

6.385U1)−1

(4)

2 −1.3

))

)

(5) (6)

− τa)

ma1.02))

(7)

τda = 0.79τo τg τw τaa Fc (1 − τas )/(1 − ma + ma1.02))

(8)

ma = mr (exp (−z /8430))

(9)

τdr = 0.79τo τg τw τaa (0.5(1 − τr )/(1 − ma +

mr = (cos (θz ) + 0.15(93.885 − θz )−1.253)−1

(10)

U1 = wmr

(11)

U3 = lmr

(12)

τas = τa/ τaa

(13)

where τr, τo, τg, τw, τa are the broadband transmittances of Rayleigh scattering, ozone, permanent gas absorption, water vapor and aerosol extinction, respectively. τaa is the transmittance of direct radiation due to aerosol absorptance. τdr and τda are the Rayleigh-scattered and aerosol scattered diffuse 12

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

transmittance after the first pass through the atmosphere, respectively. mr is air mass. ma is pressure-corrected air mass. θz is solar zenith angle. z is the altitude above sea level in meters. w and l, both derived from MODIS atmosphere products, are precipitable water and ozone layer thickness, respectively. β is the Angstrom turbidity coefficient which can be calculated using MODIS AOD data and the method of Y. Zhang et al. (2015). Under clear sky, beam transmittances for direct radiation (τb) and diffuse radiation (τd) are:

τb = 0.975τr τo τg τw τa

(14)

τd = τdr + τda

(15)

A.2. Cloudy sky atmospheric transmittance Various cloud attenuation models, from simple parameterizations to complex radiative transfer models have been developed. Using radiative transfer models may be the most accurate but input data requirements and computation efficiency makes them less attractive for practical use, especially in mountainous areas. Simple parameterizations based on sound cloud attenuation physics and using few cloud properties are more suitable for mountainous terrain. Here we use the cloud transmittance in Xie et al. (2016) which has shown good performance. (16)

τc = τcb + τcd −τ/μ0

, where τ is the cloud optical thickness in a sky sector and μ0 is where τcb is cloud extinction and is based on the Beer-Bouguer-Lambert law τcb = e the cosine of the incidence angle at the centroid of the sky sector. For cloudy sky, diffuse transmittance is the key component and involves different processes for water cloud and ice cloud. As such, the parameterization in Xie et al. (2016) is adopted for water and ice cloud, respectively.

(

τ ⎛ − log10 τp (De , μ 0 ) τcd = f (De , τ , μ0 ) Tdu, p (De , μ0 ) exp ⎜ D (μ0 ) ⎜ ⎝

) ⎞⎟ ⎟ ⎠

(17)

where, f(De, τ, μ0), Tdu, p(De, μ0), τp(De, μ0) and D(μ0) are parameterized for water and ice cloud, respectively. For water clouds, they are parameterized as:

τp (De , μ0 ) = [2.885 + 0.002(De − 60)] μ0 − 0.007347

Tdu, p (De , μ0 ) = 0.7846[1 + 0.0002(De − 60)] μ00.1605

D (μ0 ) = −0.644531μ0 + 1.20017 + 0.129807μ0−1 − 0.00121096μ0−2 + 1.52587 × 10−7μ0−3

f (De , τ , μ0 ) = 1 + sinh[0.012μ0 (τ − τp (De , μ0 ))] and for ice cloud they are parameterized as:

2.8487μ0 − 0.0029(De ≤ 26) τp (De , μ0 ) = ⎧ [2.8335 + 0.006(100 − De )] μ0 − 0.00612(De > 26) ⎨ ⎩ Tdu, p (De , μ0 ) = 0.756μ00.0883

D (μ0 ) = −0.0549531μ0 + 0.617632 + 0.17876μ0−1 − 0.002174μ0−2 f (De , τ , μ0 ) = 1 + sinh[0.01μ0 (τ − τp (De , μ0 ))] A.3. Isotropic model For isotropic model, diffuse radiation is assumed isotropically distributed over the hemispherical sky. Direct radiation is calculated as

Isb = I0 ∗ τb ∗ τcb ∗ cos (θs ) ∗ sf

(18)

where I0 is extraterrestrial radiation, τb is beam transmittances for direct radiation in clear sky. τcb is cloud extinction and is based on the BeerBouguer-Lambert law τcb = e−τ/μ0, and θs is incidence angle. sf is surrounding terrain shading factor and it is 0 if there is shadow, otherwise it is 1. For diffuse radiation, if it is clear sky, the diffuse radiation is calculated using our diffuse model using an isotropic weight. If it is cloudy sky, diffuse radiation from a sky sector is calculated as:

Isi = Id ∗ wi ∗ cos (INAi ) ∗ SGi

(19)

Id = I0 τb τcd

(20)

where Id is the diffuse radiation on horizontal surface, INAi is the incidence angle at the centroid of the ith sky sector. wi is the isotropic weight in the ith sky sector which can be calculated based on Fu and Rich (1999). SGi is the fraction of the visible part of the ith sky sector and is calculated by overlaying skyshed and sky map. τb is direct transmittance in clear sky and τcd is diffuse transmittance for cloudy sky which is calculated based on Xie et al. (2016) and the cloud properties in small window around sun position (Fig. 8a) are used. A.4. Downscaling methods for CERES SYN and CM SAF We downscaled CERES SYN and CM SAF datasets using the method developed by Antonanzas-Torres et al. (2018, 2014) except we didn't use their 13

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al.

error interpolation process. Global (GH) and direct (DH) horizontal radiation are provided in the two datasets. Diffuse horizontal (DFH) radiation is obtained as the difference between the global and direct radiation. The diffuse horizontal radiation is further divided into two components: isotropic diffuse radiation (DFHiso) and anisotropic diffuse radiation (DFHani) based on the anisotropic index (ka) which is defined as the ratio of direct radiance Id and extra-terrestrial irradiance Ie as: (21)

DFH = DFHiso + DFHani

ka =

Id Ie

(22)

Based on the anisotropy index ka, isotropic diffuse radiation is calculated as:

DFHiso = DFH ∗ (1 − ka)

(23)

DFHani also referred to circumsolar diffuse radiation, considers the incoming portion from the circumsolar disk and can be calculated as: (24)

DFHani = DFH ∗ ka

Topographical effects include terrain shading effect on both direct radiation and diffuse radiation which are modeled using a binary terrain shading factor (Shadow) and the sky view factor (SVF) in the Solar Analyst module of ArcGIS, respectively. Downscaled global solar radiation (Gd) is calculated as: (25)

Gd = DH ∗ Shadow + DFHiso ∗ SVF + DFHani ∗ Shadow

occurrences of mountain permafrost. Permafr. Periglac. Process. 3 (2), 139–142. https://doi.org/10.1002/(ISSN)1099-1530. Gelaro, R., McCarty, W., Suarez, M., et al., 2017. The modern-era retrospective analysis for research and applications, version 2 (MERRA-2). J. Clim. 30, 5419–5454. https:// doi.org/10.1175/JCLI-D-16-0758.1. Helbig, N., Lowe, H., Lehing, M., 2009. Radiosity approach for the shortwave surface radiation balance in complex terrain. Journal of The Atmosphere Sciences 66, 2900–2912. Hinkelman, L., Lapo, K., Cristea, N., Lundquist, J., 2015. Using CERES SYN surface irradiance data as forcing for snowmelt simulation in complex terrain. J. Hydrometeorol. 16, 2133–2152. https://doi.org/10.1175/JHM-D-14-0179.1. Houborg, R., Soegaard, H., Emmerich, W., Moran, S., 2007. Inferences of all-sky solar irradiance using Terra and Aqua MODIS satellite data. Int. J. Remote Sens. 28 (20), 4509–4535. Iqbal, M., 1983. An Introduction to Solar Radiation. Academic Press, Toronto, ON. Kumar, L., Skidmore, A., Knowles, E., 1997. Modelling topographic variation in solar radiation in a GIS environment. Int. J. Geogr. Inf. Sci. 11 (5), 475–497. https://doi. org/10.1080/136588197242266. Li, X., Williams, M.W., 2008. Snowmelt runoff modelling in an arid mountain watershed, Tarim Basin, China. Hydrol. Process. 22 (19), 3931–3940. https://doi.org/10.1002/ hyp.7098. Li, X., Cheng, G., Chen, X., Lu, L., 1999. Modification of solar radiation model over rugged terrain. Chin. Sci. Bull. 44 (15), 1345–1350. Li, X., Zhang, S., Chen, Y., 2016. Error assessment of grid-based diffuse solar radiation models. Int. J. Geogr. Inf. Sci. 30 (10), 2032–2049. Liu, R., Liu, Y., Du, B., 2002. Cloud climatic characteristics of the Tibetan Plateau from ISCCP data. Journal of Nanjing Institute of Meteorology 25 (2), 226–234. Liu, M., Bardossy, A., Li, J., Jiang, Y., 2012. GIS-based modelling of topography-induced solar radiation variability in complex terrain for data sparse region. Int. J. Geogr. Inf. Sci. 26 (7), 1281–1308. https://doi.org/10.1080/13658816.2011.641969. Liu, H., Tang, S.H., Zhang, S.L., Hu, J.Y., 2015. Evaluation of MODIS water vapour products over China using radiosonde data. Int. J. Remote Sens. 36 (2), 680–690. Lu, X., 2016. Estimation of the Instantaneous Downward Surface Shortwave Radiation Using MODIS Data in Lhasa for All-Sky Conditions. International Development, Community and Environment (IDCE), pp. 99. Marsh, C.B., Pomeroy, J., Spiteri, R.J., 2012. Implications of mountain shading on calculating energy for snowmelt using unstructured triangular meshes. Hydrol. Process. 26, 1767–1778. https://doi.org/10.1002/hyp.9329. Mayer, B., 2009. Radiative transfer in the cloudy atmosphere. The European Physical Journal Conferences 1, 75–99. Menzel, W.P., Frey, R.A., Baum, B.A., 2015. Cloud top properties and cloud phase algorithm theoretical basis document. available at. http://modis-atmos.gsfc.nasa.gov/_ docs/MOD06-ATBD_2015_05_01.pdf. Miesch, C., Briottet, X., Kerr, Y.H., Cabot, F., 1999. Monte Carlo approach for solving the radiative transfer equation over mountainous and heterogeneous areas. Appl. Opt. 38 (36), 7419–7430. Qin, J., Tang, W., Yang, K., Lu, N., Niu, X., Liang, S., 2015. An efficient physically based parameterization to derive surface solar irradiance based on satellite atmospheric products. Journal Geophysical Research: Atmospheres 120, 4975–4988. https://doi. org/10.1002/2015JD023097. Roupioz, L., Jia, L., Nerry, F., Menenti, M., 2018. Estimation of daily solar radiation budget at kilometer resolution over the Tibetan Plateau by integrating MODIS data products and a DEM. Remote Sens. 8, 504. Ryu, Y., Jiang, C., Kobayashie, H., Dettof, M., 2018. MODIS-derived global land products of shortwave radiation and diffuse and total photosynthetically active radiation at 5 km resolution from 2000. Remote Sens. Environ. 204, 812:825. Sayer, A.M., Munchak, L.A., Hsu, N.C., Levy, R.C., Bettenhausen, C., Jeong, M.J., 2014. MODIS collection 6 aerosol products: comparison between Aqua's e-Deep Blue, Dark

References Aguilar, C., Herrero, J., Polo, M.J., 2010. Topographic effects on solar radiation distribution in mountainous watersheds and their influence on reference evapotranspiration estimates at watershed scale. Hydrol. Earth Syst. Sci. 14, 2479–2494. https://doi.org/10.5194/hess-14-2479-2010. Allen, R.G., Trezza, R., Tasumi, M., 2006. Analytical integrated functions for daily solar radiation on slopes. Agric. For. Meteorol. 139 (1–2), 55–73. https://doi.org/10.1016/ j.agrformet.2006.05.012. Antonanzas-Torres, F., Martinez-de-Pison, F., Antonanzas, J., Perpinan, O., 2014. Downscaling of global solar irradiation in complex areas in R. Journal of Renewable and Sustainable Energy 6, 063105. Antonanzas-Torres, F., Martínez-de-Pisón, F., Antonanzas, J., Perpinan, O., 2018. Downscaling of Global Solar Irradiation in R. Bair, E., Rittger, K., Davis, R., Painter, T., Dozier, J., 2016. Validating reconstruction of snow water equivalent in California's Sierra Nevada using measurements from the NASA Airborne Snow Observatory. Water Resour. Res. 52, 8437–8460. https://doi. org/10.1002/2016WR018704. Bo, Y., Wang, Y., Li, J., Wang, C., 2016. Temporal and spatial variation features of cloud water and its relation to precipitation over the Tibetan Plateau. J. Glaciol. Geocryol. 38 (6), 1679–1690 (In Chinese). Borbas, E., Seemann, S.W., Kern, A., Moy, L., Li, J., Gumley, L., Menzel, W.P., 2011. MODIS atmospheric profile retrieval algorithm theoretical basis document collection 6. Available at: http://modis-atmos.gsfc.nasa.gov/_docs//MOD07_atbd_v7_ April2011.pdf. Brunger, A., Hooper, F., 1993. Anisotropic sky radiance model based on narrow field of view measurements of shortwave radiance. Sol. Energy 51 (1), 53–64. Chen, Y., Hall, A., Liou, K.N., 2006. Application of three-dimensional solar radiative transfer to mountains. J. Geophys. Res. 111, D21111. Chen, L., Yan, G.J., Wang, T.X., Ren, H.Z., Calbo, J., Zhao, J., Mckenzie, R., 2012. Estimation of surface shortwave radiation components under all sky conditions: modeling and sensitivity analysis. Remote Sens. Environ. 123, 457–469. Chen, M., Zhuang, Q.L., He, Y.J., 2014. An efficient method of estimating downward solar radiation based on the MODIS observation for the use of land surface modeling. Remote Sens. 6, 7136–7157. https://doi.org/10.3390/rs6087136. Corripio, J.G., 2003. Vectorial algebra algorithms for calculating terrain parameters from DEMs and solar radiation modelling in mountainous terrain. Int. J. Geogr. Inf. Sci. 17 (1), 1–23. https://doi.org/10.1080/713811744. Dozier, J., Frew, J., 1990. Rapid calculation of terrain parameters for radiation modeling from digital elevation data. IEEE Trans. Geosci. Remote Sens. 28 (5), 963–969. https://doi.org/10.1109/36.58986. Dubayah, R., Rich, P.M., 1995. Topographic solar-radiation models for GIS. Int. J. Geogr. Inf. Syst. 9 (4), 405–419. https://doi.org/10.1080/02693799508902046. Durr, B., Zelenka, A., 2009. Deriving surface global irradiance over the Alpine region from METEOSAT Second Generation data by supplementing the HELIOSAT method. Int. J. Remote Sens. 30 (22), 5821–5841. Fitzpatrick, M.F., Warren, S.G., 2005. Transmission of solar radiation by clouds over snow and ice surfaces. Part II: cloud optical depth and shortwave radiative forcing from pyranometer measurements in the southern ocean. J. Clim. 18, 4637–4648. Fitzpatrick, M.F., Brandt, R.E., Warren, S.G., 2004. Transmission of solar radiation by clouds over snow and ice surfaces: a parameterization in terms of optical depth, solar zenith angle, and surface albedo. J. Clim. 17, 266–275. Fu, P., Rich, P.M., 1999. Design and implementation of the Solar Analyst: An ArcView extension for modeling solar radiation at landscape scales. In: Paper Presented at 19th Annual ESRI User Conference, San Diego, CA, Available from: http:// proceedings.esri.com/library/userconf/proc99/proceed/papers/pap867/p867.htm. Funk, M., Hoelzle, M., 1992. A model of potential direct solar radiation for investigating

14

Remote Sensing of Environment 231 (2019) 111239

S. Zhang, et al. Target, and “merged” data sets, and usage recommendations. Journal of Geophysical Research: Atmospheres 119 (24), 13965–13989. Šúri, M., Hofierka, J., 2004. A new GIS-based solar radiation model and its application to photovoltaic assessments. Trans. GIS 8, 175–190. https://doi.org/10.1080/ 13658810802022806. Tian, Y.Q., Davies-Colley, R.J., Gong, P., Thorrold, B.W., 2001. Estimating solar radiation on slopes of arbitrary aspect. Short communication. Agric. For. Meteorol. 109, 67–74. Van Laake, P.E., Sanchez-Azofeifa, G.A., 2004. Simplified atmospheric radiative transfer modelling for estimating incident PAR using MODIS atmosphere products. Remote Sens. Environ. 91, 98–113. Van Laake, P.E., Sanchez-Azofeifa, G.A., 2005. Mapping PAR using MODIS atmosphere products. Remote Sens. Environ. 94, 554–563. Wang, Q., Tenhunen, J., Schmidt, M., Otieno, D., Kolcun, O., Droesler, M., 2005. Diffuse PAR irradiance under clear skies in complex alpine terrain. Agriculture Forest. Meteorology 128, 1–15. Wang, D., Liang, S., Liu, R., Zheng, T., 2010. Estimation of daily-integrated PAR from sparse satellite observations: comparison of temporal scaling methods. Int. J. Remote Sens. 31 (6), 1661–1677. Wang, T., Yan, G., Mu, X., Jiao, Z., Chen, L., Chu, Q., 2018. Toward operational shortwave radiation modeling and retrieval over rugged terrain. Remote Sens. Environ.

205, 419–433. Xie, Y., Sengupta, M., Dudhia, J., 2016. A fast all-sky radiation model for solar applications (FARMS): algorithm and performance evaluation. Sol. Energy 135, 435–445. Yang, K., Huang, G.W., Tamai, N., 2001. A hybrid model for estimating global solar radiation. Sol. Energy 70, 13–22. Yu, Y.C., Shi, J.C., Wang, T.X., Letu, H.S., Yuan, H.S., Zhou, W., Hu, L., 2019. Evaluation of the Himawari-8 Shortwave Downward Radiation (SWDR) product and its comparison with the CERES-SYN, MERRA-2, and ERA-Interim datasets. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 70, 519–532. https://doi.org/10.1109/JSTARS.2018.2851965. Zhang, H.L., Xin, X.Z., Li, L., Liu, Q.H., 2013. An improved parametric model for simulating cloudy sky daily direct solar radiation on titled surface. IEEE Journal of selected topics in applied earth observation and remote sensing 6 (1), 180–187. Zhang, Y., Li, X., Bai, Y., 2015. An integrated approach to estimate shortwave solar radiation on clear-sky days in rugged terrain using MODIS atmospheric products. Sol. Energy 113, 347–357. Zhang, S., Li, X., Chen, Y., 2015. Error assessment of grid-based direct solar radiation models. Int. J. Geogr. Inf. Sci. 30 (10), 2032–2049. https://doi.org/10.1080/ 13658816.2015.1055273.

15