Modeling of Diffuse Solar Radiation and Impact of Complex Terrain over Pakistan Using RS / GIS

Diffuse solar radiation is subject to the combined influence of ground and sky factors, such as topography, geography of the area and cloud cover. This study attempts to quantify the impacts of topography, sky factors and the cloud cover on the distribution of diffuse solar radiation over Pakistan. Distributed modeling approach by considering anisotropy scattering mechanism was adopted. Digital elevation model and observed data are used to derive average monthly diffuse solar radiation values over the rugged terrains of Pakistan. Extraterrestrial solar radiation model, sky view factor model (openness model) and digital elevation model (DEM) are applied to investigate the impacts of ground factors, while diffuse solar radiation model for horizontal surface was considered for sky factors. Furthermore, corrected MODIS cloud fraction data are incorporated using GIS plat form. Results show that the highest amount of diffused solar radiation occurs during the monsoon months along the eastern side of the River Indus, when the sky is covered by clouds of various heights and densities. The variation due to topography is evident in mountainous areas, particularly in the North Pakistan and over the Baluchistan Plateau.


Introduction
Pakistan is the 33 rd largest country (803,940 km 2 ) with more than twice the size of Japan featuring heterogeneous terrain from 0 m altitude in coastal belts to 8000 m in Himalaya ranges [1].The whole country has only 5 meteorological stations measuring solar radiation and none of them have diffuse solar radiation measurements.It is, therefore, imperative to find other ways to estimate the solar radiation over such complex terrain.
Terrain is an important factor in describing the distribution of energy fluxes.Global irradiance on a slope is the sum of its direct (beam), diffuse and ground-reflected components.After the solar radiation enters the earth's atmosphere, it is partially scattered and partially absorbed.The scattered radiation is called diffuse radiation.Estimation of that is variable from day to day.Diffuse cloud radiation appears to only contribute a minor part to radiation energy from above the mid-visible to the infrared spectrum, but can contribute up to 40% of the radiation energy from the mid-visible through mid-ultraviolet spectrum [2].
It is obvious that the relative proportion of direct to diffuse radiation depends on the geographical location, season of the year, elevation from the mean sea level, and time of the day.On a clear-sky day, the diffuse component of solar radiation will be about 10% -20% of the total radiation, but during an overcast day, it may reach up to 100%.This implies, practically, that in the solar radiation and energy calculations, weather, meteorological and climatological conditions implications must be taken into consideration in addition to the astronomical.The instantaneous total radiation can vary considerably through the day depending on the cloud cover, dust concentration, humidity, etc. [3].That is why monthly mean values are often considered.
Diffuse sky irradiance under cloud-free conditions may be estimated by assuming an isotropic sky and calculating the proportion of the sky seen from a point that is using the equivalent of the view-shed operation in GIS [4].Under cloudy or partly cloudy conditions, diffuse radiation is anisotropic and may be explicitly modeled.The applications of digital elevation model (DEM) data provide a new means for diffuse solar radiation (DSR) calculation over the complex terrain [4]- [6].

Study Area
Pakistan lies between 24˚ to 37˚N latitude and 62˚ to 75˚E longitude, covering a total land area of 796,096 km 2 (Pakistan official website).The topography of Pakistan is a combination of landscapes varying from plains to deserts, forests, hills, and plateaus ranging from the coastal areas of the Arabian Sea in the south to the lofty mountains including Hindu-Kush, Karakoram and Himalaya ranges in the north.Pakistan is divided into three major geographic areas: the northern highlands, the Indus River plain, and the Baluchistan Plateau.The northern highlands of Pakistan contain the Karakoram, Hindu Kush and Pamir mountain ranges, which have some of the world's highest peaks, including K2 (8611 m/28,251 ft) and Nanga Parbat (8126 m/26,660 ft) [7].The Baluchistan Plateau lies to the West, and the Thar Desert in the East.The vast alluvial plains were laid down through ages by the River Indus and its tributaries in Punjab and Sindh.The 1609 km (1000 mi) long.The Indus River and its tributaries flow through the country from the Kashmir region to the Arabian Sea [8] [9].

Meteorological Data
The daily global solar radiation data from 1978 to 2000 and the monthly cloud cover and monthly sunshine duration data from 1979 to 2008 were obtained from Pakistan Meteorological Department (PMD).The data for cloud cover and sunshine duration are available from 21 stations, while data for solar radiation are obtained from 6 meteorological observatories of Pakistan (only Global Solar Radiation) and 5 Meteorological stations of China (Diffuse and Global).Data quality from Pakistan remains a pronounced concern of apprehension due to lack of funds and modern infrastructure with trained personals.The continuity of solar radiation data can be categorized into two groups, one from 1979 to 2000 and other from 2001 to 2008.The data in group one are used for simulation due to relatively more consistency in terms of availability.

Satellite and DEM Data
DEM data for Pakistan are obtained from Shuttle Radar Topography Mission (SRTM) (ftp://e0mss21u.ecs.nasa.gov/srtm/).The data of 3 Arc-Second, almost with 900 m cell size, are employed in this research.The cloud-fraction data were retrieved from Moderate Resolution Imaging Spectroradiometer (MODIS) Level-2 product to estimate cloud cover over Pakistan at are solution of about 1 km.

Methodology
From the models of Hay (1979) [10], Gueymard (1987) [11], Skartveit et al. (1985) [12] and Perez et al. (1986) [13], we can calculate the global and diffuse radiation on an inclined (slope) surface.Hay's model distinguishes from the three other models by its simple structure and the small number of input values.
The Hay's model is structured as follows where , which is also known as the conversion factor; V is topographic openness.
If the sky is overcast, theoretically the direct irradiation will be zero, therefore, DSR is an important component of the surface-received global solar radiation so that both are closely related [14].Referring to the researches of decomposition models, we present the function for simulating d H , as fol- lows, ( ) ( ) There, a, b and c are empirical coefficients, s the percentage of sunshine and H the global solar radiation over a horizontal surface and d H is the diffuse solar radiation at horizontal plane.Equation ( 2) is of explicit physical implication.For entirely overcast weather (s = 0), there is no direct incoming solar radiation such that solar radiation striking a horizontal surface is composed completely of diffused solar radiation, i.e., d H H = .On a fine day ( ) the horizontal surface received solar radiation is all made up mainly of direct radiation, with DSR component reaching a minimum [7].
The empirical coefficients a, b and c of Equation ( 2) are determined by measurements [15].For investigating the time dependent characteristics of empirical coefficients, we take advantage of two kinds of data sets for establishing separately, a unified and a monthly model.
The unified model (Table 2) is formulated with all monthly d H data from all stations concerned as one sample.The monthly model (Table 1) is constructed on the data of the same respective month from all the stations.
From Table 2, by considering the changing characteristics of empirical coefficients with time, monthly models can improve the accuracy of   Here n is the sample length.
diffuse solar radiation.As we know, there are a limited number of meteorological observatories having data of horizontal global radiation.Conventional stations keep only percentage of sunshine measurements in a consequence, horizontal global radiation must be acquired before fitting DSR by Equation (2).
A closed model for horizontal solar radiation simulation is constructed using Equation ( 2) combined with the model for horizontal direct radiation simulation [16].This guarantees that the sum of fitted horizontal direct and diffuse radiation equals to the fitted global radiation on the horizontal surface.

Various studies show that global solar radiation ( )
H and sunshine percentage ( ) s are closely related [17] [18]; typically, a linear model is used for their relation; ( ) There and c c a b are empirical coefficients [16].
Through combination of Equation (2) with Equation (3), we get the following: Using Equation ( 4) we need only data of sunshine percentage as input to simulate diffuses solar radiation

Model of Topographic Openness
Any particular location over the complex terrain may be obstructed by topography, which may reduce the dif-fuse irradiance from its corresponding sky directions.This obstruction can either from self-shadowing by the slope (shading) itself or from adjacent terrain (shadowing).The topographic openness model for sky view factor V can be calculated to give the ratio of diffuse sky irradiance at a point to that on an unobstructed horizontal surface [4].
In rugged terrains, slope-received diffuse solar radiation from particles in the sky is associated with shielding of the surrounding terrain.In previous studies, typically, the topographic openness V for a slope is as follow- where α is the grade of slope.
The model (Equation ( 6)) gives just the self-shielding of a single slope with an infinite length of itself.However, the topographic openness is also dependent on the inter-shielding of terrains around.Thus, the shielding effect requires numerical integration within 2π azimuth.The topographic openness model for a point P over the rugged surface is given below, 1) With ∆Φ (degrees) as the step length of the azimuth, we divide 2π circumference into n azimuths where ( ) int is a function which gets the integer part of the value 360 ∆Φ .
2) With south as the start direction, i.e. 0 0 Φ = , and ∆Φ as step length of azimuth, in a clockwise direction, the n azimuths can be given by 0 0,1, 2, , 1 3) Calculation of the openness i V at the direction i Φ , Fu (1983) [20] proved that openness of point P at the direction of i Φ is where i α denotes the largest elevation angle of P at the direction of i Φ , i.e., the largest shielding angle caused by terrain at the direction of i Φ .Starting from P and drawing a straight line i L in the direction of i Φ , then, we can compare elevation an- gles of P at the points along i L to find the maximum.In calculation, terrain elevations are obtained from the DEM that consists of grids with constant length and width.In our proposed model, a space step length L ∆ is used in determining the elevation angles along i L .The minimal grid length and width are used as L ∆ which is denoted as: where sizex is the resolution of DEM in the -axis x , while sizey the resolution of DEM in the -axis y . With L ∆ as the space step length, the elevation of point j is ( ) where L ω ∆ and y L ∆ are the increase of the space step lengths in the and x y direction, respectively, with ( ) ( ) sin and cos Therefore, for point P , the shielding angle of point j is as follows; ( ) ( ) For point P , the greatest shielding angle caused by terrains in the direction of The length of i L is not essentially infinite, but taking a certain shielding radius R is enough to meet the needs of calculation and N in Equation (11) and Equation ( 13) stands for the number of calculations, depend- ing on the screening radius R and the step length L ∆ .
Since tangential function increases monotonously in the bound ( ) , the following expressions are used to increase efficiency of calculation.

(
) ( ) , Then Equation ( 9) is employed to compute the openness i V of the point P at the direction of i Φ .
The impossibility to ensure and in integers leads to the that ( ) has to be obtained with resample method.Bilinear interpolation method [21] is used in the calculation.4) Calculation of topographic openness of point P : we make averaging of openness i V at i Φ along 2π circumference, resulting in the topographic openness of P , that is 1 0

Diffuse Solar Radiation Distributed over the Irregular Topography
By putting the values in Equation ( 1), we have the spatial distribution pattern of diffuse solar radiation over the rugged topography.Figure 1 shows annual diffuse solar radiation spatial distribution.The DSR quantity randomly decreases from north-east to south-west and has a gradual decrease towards south with minimum values lyingin Kharan desert (Baluchistan).The same tendency has also been observed in Thar (Sindh), Cholistan and Thal (Punjab) deserts and maximum values have been observed along the eastern side of Indus basin due to maximum cloud-cover associated with Indian monsoon weather system over the area.There is an uneven variation in mountains of Gilgit Baltistan and Kashmir regions in northern Pakistan due to varying slopes in term of direction and height, i.e. aspect ratios and complex topographic effects on weather.

Seasonal Variations
The seasonal (winter and summer) variations have been shown in Figure 2. The color ramping shows a clear difference between winter and summer.The values of DSR are higher during the summer monsoon than winter season.During summer, DSR has an increasing trend from west to east, while in winter there is no significant variation in the plains of Pakistan.

Topographic Impact on Diffuse Solar Radiation
To analyze the effect of topography on diffuse solar radiation over the region of Pakistan, we calculate the difference between diffuse solar radiations over rugged terrain and diffuse solar radiation on flat land as shown in Figure 3.The result shows a prominent difference in the northern areas including Gilgit Baltistan, Potwar and Baluchistan Plateaus (purple color); on the contrary, there is less difference in the plains, particularly in the areas comprised of deserts

Conclusions and Limitations
There is a significant impact of topography on the distribution of diffuse solar radiation over the uneven terrain, particularly in Gilgit Baltistan, Potwar Plateau and Baluchistan Plateau (Ref: Figure 1).Figure 4 represents the annual sum of diffuse solar radiation over Pakistan.The central and south eastern part of Pakistan have more homogenous distribution of diffuse solar radiation as compared to North and North-West part due to the tomography orientation and aspects.
Since there is no observed data available for diffuse solar radiation in the study area, therefore, the general  idea about the distribution of temperature and precipitation (as shown in Figure 5) could help us to have better understanding and some grounds for the diffuse radiation pattern and its spatial distribution validity.It is evident from the temperature distribution that, it follows the terrain patterns and displays latitudinal variations.It is worth to be noted that temperature is relatively lower on the western side of the Indus River as compared to

HH
αβ is the diffuse solar radiation quantity on rugged terrain; d H is the diffuse solar radiation quantity reaching the flat surface (horizontal plane); H is the global solar radiation quantity on a horizontal plane; b H is the direct solar radiation quantity on a horizontal plane; 0 is extraterrestrial solar radiation quantity on a horizontal surface; b R is the ratio of extraterrestrial solar radiation quantity on rugged terrain ( ) d H αβ to that on a horizontal surface ( ) 0 H , i.e.

3 . 2 . 2 .H
model approaches to the circumsolar model.Considering the topographic openness model (sky view factor) V and conversion factor b R , Equation (1) is a universal model for calculating diffuse solar radiation (DSR) quantity on rugged terrain.Model of Fitting Horizontal Surface Diffuse Solar Radiation ( ) b Even in mountainous area, meteorological stations usually situate at open flat having no obstacles in certain slope.Therefore, observations show the horizontal-surface characteristics.Our fitting model is based on station d H observations.

dH
simulation effectively.Hence for the locations having data of horizontal global radiation ( ) H , we use monthly model with conjunction of percentage sunshine duration to simulate the percentage of sunshine duration is being calculated through incorporating digital elevation data to possible sunshine duration data[19], i.e. coefficients, being calculated through MODIS cloud products corrected by observed cloud data.This process is done using ArcGIS.By considering monthly percentage of sunshine measurements (1978 to 2008) from the 21 stations throughout Pakistan, calculations of Equation (4) lead to the fittings of d H on a monthly basis for each of the stations from 1978 to 2008.Then, using inverse distance weighting (IDW) interpolation method, spatial distribution of monthly d H quantity with resolution of 1km by 1kmfor Pakistan has been generated.

Figure 4 .
Figure 4. Annual distribution of diffuse solar radiation over the rugged terrain of Pakistan.

Figure 5 .
Figure 5. Mean annual temperature and precipitation during 1971-2000 (courtesy of Pakistan Meteorological Department).easternside[22].The relatively high temperature causes more convective cloud formation and may feed by moisture from the Arabian Sea and the Bay of Bengal during monsoon season in particularly.Furthermore, the North part of Pakistan is relatively cooler due to the latitudinal location, topographic feature and cold wind from west and Siberian high.Similarly this area gets maximum precipitation due to the moisture being feed along the foothills of Himalaya from Bay of Bengal and Arabian Sea which could get cooling for condensation either from Siberian High or so called western disturbance[23] [24].Now if we describe the distribution of diffuse solar radiation in conjunction with temperature and precipitation patterns, we can conclude that during the monsoon season, diffuse radiation has increasing trend from western side of the Indus River to Eastern side of the Indus River due to distribution of cloud cover and or ographic patterns.The maximum values of diffuse radiation lays in the upper Khyber-Pakhtunkhwa ( ) 2 2628 -2867 M J m − ⋅ ⋅ , while minimum values ( ) 2 1068 -1562 M J m − ⋅ ⋅ lays in Baluchistan Plateau.

Table 1 .
Statistics of monthly model for diffuse solar radiation.

Table 2 .
Statistics of unified Model for diffuse solar radiation.