Solar Irradiance and Pan Evaporation Estimation from Meteorological Satellite Data

Knowledge about spatial and temporal variations in surface global solar radiation (GSR) and evaporative water loss from the ground are important issues to many researches and applications. In this study empirical relationships suitable for Taiwan were established for GSR retrieval from geostationary satellite images using the Heliosat method for the period from 2011 2013. The derived GSR data has been used to generate consecutive maps of 10-day averaged pan evaporation (Epan) as the basis to produce regional ET estimation using a strategy that does not require remote sensed land surface temperatures (LST). The retrieved daily GSR and the derived 10-day averaged Epan were validated against pyranometer and class-A pan measurements at selected Central Weather Bureau (CWB) stations spread across various climatic regions in Taiwan. Compared with the CWB observed data the overall relative mean bias deviations (MBD%) and root mean square differences (RMSD%) in daily solar irradiance retrieval were about 5 and 15%, respectively. Seasonally, the largest MBD% and RMSD% of retrieved daily solar irradiance occur in spring (9.5 and 21.3% on average), while the least MBD% (-0.3% on average) and RMSD% (9.7% on average) occur in autumn and winter, respectively. For 10-day averaged Epan estimation, the mean MBD% and RMSD% for stations located in the coastal plain areas were 0.1 and 16.9%, respectively. However, in mountainous areas the mean MBD% and RMSD% increased to 30.2 and 34.5%, respectively. This overestimation was due mainly to the large differences in surrounding micro-environments between the mountainous and plain areas.


IntroDuctIon
Solar radiation is the driving force in the transportation and exchange of mass, energy and momentum at the interface between the Earth's surface and the atmosphere, which makes it an important input parameter for many meteorological, climatic and hydrological related models (e.g., Liang et al. 1994;Berbery et al. 1999).As the energy source of photosynthesis reaction, solar radiation serves as a critical factor in determining the gross primary production of terrestrial and marine ecosystems (e.g., Monteith 1972;Roebeling et al. 2004).Knowledge of daily irradiation received at a particular site is also essential for many energy related applications such as solar energy utilization, architectural design and building environment planning (e.g., Oliver and Jackson 2001;Kumar and Umanand 2005).Accurate data on the spatial and temporal variations in surface global solar radiation (GSR) are therefore absolutely necessary for many basic researches and engineering applications.
Over flat plains and on clear days incoming solar radiation may be calculated from a sophisticated radiative transfer model which considers the effects of scattering and absorption on solar radiation due to water vapor, O 3 , O 2 , and CO 2 , as well as other air molecules during its passage through the atmosphere together with information on the solar constant, the Earth-Sun distance, and the solar zenith angle (e.g., Chou and Lee 1996).The calculated solar irradiance is strongly affected by the amount of clouds and aerosol loadings at a given location and time, which makes the importance of on-site GSR measurements irreplaceable.Over complex terrains the significance of terrain shading on the GSR increases with decreasing scales of spatial resolution, particularly in the winter season and at scales smaller than 5 km (Lai et al. 2010).
Traditionally, the solar radiation received at a place is interpolated and/or extrapolated from a pyranometer network.Although measurement error for a properly calibrated and maintained pyranometer is generally accepted to be 3 -5%, the irradiance error estimated from nearby ground-measuring stations increases with the distance between stations.Studies by Perez et al. (1997) and Zelenka et al. (1999) indicated that, within only 20 -30 km, the relative root mean square differences (RMSD%) already reach 10 -15 and 20 -25% for daily and hourly GSR estimates, respectively.The micro-variability in the irradiation field in conjunction with measurement uncertainties causes a sharp rise in root mean square difference (RMSD) within the first 10 km.Hence, extrapolation, even from a nearby station, may not be satisfactory for time intervals such as hourly and daily.
The averaged distance between Central Weather Bureau (CWB) pyranometer stations in Taiwan is about 20 km in the northern region, and about 60 km in the southern and eastern regions.Most of these pyranometer stations are located in flat urban areas along the coast.There are only three CWB pyranometer stations, all concentrated in Central Taiwan, measuring incident solar radiation at mountainous areas that cover almost 2/3 of Taiwan.Considering the spacing between pyranometer stations, the geographical distribution of mountainous areas on the island, and the spatio-temporal heterogeneity of cloud cover caused by the combined effects of prevailing winds, diurnal circulation and orographic lifting, it is doubtful that the current pyranometer network in Taiwan can provide good estimates of hourly and daily GSR.Furthermore, according to Lin et al. (2004), most installed pyranometers in Taiwan are not regularly calibrated.The accumulated errors of solar irradiance measurements from 1986 -2002 were systematically lower by 15 -30%.Therefore, a high-quality short-duration solar radiation atlas is not yet available for the above mentioned basic researches and engineering applications.
Since cloud cover plays a major role in determining the transmitted solar irradiance, geostationary meteorological satellites have become a valuable tool for estimating GSR on large scales because of their capabilities in providing high spatial and temporal resolution cloud information.The reported RMSD% of GSR estimates derived from satellites were typically in the range of 20 -25% and 10 -15% at hourly and daily intervals, respectively (Beyer et al. 1996;Zelenka et al. 1999).These errors are comparable to the accuracies obtained by interpolation from a pyranometer network with stations 20 -30 km apart.It has been shown that GSR estimates made from satellite data offer relatively equal quality within the observed area and the accuracy of the estimates does not depend on the geographical area (Zelenka et al. 1999).Therefore, GSR estimates from geostationary satellites, though inherently less accurate than ground-based measurements, may be more suitable to generate site/time specific data at any arbitrary locations and times.
Currently most of the operational calculation schemes to retrieve solar irradiance from geostationary satellite images are semi-empirical and based on statistical methods.The Heliosat algorithm, originally proposed by Cano et al. (1986) to estimate GSR from visible band of meteorological satellites images, has now become the most popular methodology (Pagola et al. 2014).This method has been modified and improved with 4 major versions and many minor modifications.Up to Heliosat-3 the basic idea is that the amount of GSR over an area is statistically related to the cloud cover above.Therefore, a cloud index at each pixel in the image should first be calculated based on the measured albedo and prior determined reference albedos under clear and thick cloud conditions specific to that particular pixel.The GSR is then computed from a set of linear relationships describing changes in incident extraterrestrial solar radiation transmitted to the ground with the derived cloud cover index.The cloud index is basically a relative measure of the reflectivity to the darkest and brightest points on the satellite image.In Heliosat-2 version a histogram technique is applied to determine the ground reflectivity for each slot (images acquired at the same UTC-time of day belong to the same ''slot'') and month to maintain the sun-ground-satellite angle fairly constant in order to improve the accuracy of the derived cloud cover index (Rigollier et al. 2004).In Heliosat-3 version the ground reflectivity is parameterized as a function of the coscattering angle, i.e., the angle between the directions to the sun and satellite as seen from the ground, to replace the time consuming histogram technique which has to be applied separately to each slot and month (Dagestad and Olseth 2007).However, Heliosat-4 version is based on direct modeling of solar radiation propagation through the atmosphere and was developed jointly by MINES ParisTech and the German Aerospace Center (DLR).The required inputs for Heliosat-4 are atmospheric properties for the clear sky, the ground albedo, and cloud properties that are derived from an appropriate processing of images taken by the Meteosat Second Generation (MSG) satellite (Wald 2014).
Evapotranspiration (ET) describes the rate of water loss to the atmosphere through evaporation and transpiration processes, and it is a key component in the hydrological cycle and surface energy balance computation (Brutsaert 1982).On a global basis the mean ET from land surface accounts for more than half of the total precipitation (Chahine 1992;Oki and Kanae 2006).It is therefore crucial to have reliable ET information for applications such as water resources management, irrigation planning, numerical modeling of atmospheric, and hydrological processes, etc. Factors affecting ET rates include solar radiation, wind speed, air temperature and humidity, soil water content, vegetation type, growth stage, planting density, and management practices (Rosenberg et al. 1983;Monteith and Unsworth 1990).Many ground-based techniques, such as weighing lysimeter, Bowen ratio, and eddy covariance (Dugas et al. 1991;Wright 1991;Fritschen and Fritschen 2005;Meyers and Baldocchi 2005), were developed to measure ET from a given site.However, these point measurements cannot represent ET at regional scales when considering the complexities existing over a heterogeneous large area.Deployment of a ground-based measurement network is often costly, time consuming, and laborintensive.Malfunctions in instruments and accessibility to difficult terrain also generate many data collection problems for long term observation in remote areas.
Remote sensing techniques have emerged as a very useful tool to map regional ET because of its capabilities in providing repetitive and spatially contiguous images and without site accessibility issues in areas to be surveyed (Jiménez et al. 2011).A substantial amount of algorithms, varied greatly in input requirements and main assumptions, have been proposed over the past few decades to provide ET estimations at scales varying from diurnal to annual and from local to global (Kalma et al. 2008;Li et al. 2009).Those remote sensed ET models can be categorized as empirical or physically-based approaches.The empirical approach, generally for weekly to seasonally estimates, employ empirical relationships and input data directly derived from remote sensing observations (e.g., Seguin and Itier 1983;Menenti and Choudhury 1993;Gillies et al. 1997;Yuan and Bauer 2007).The physically-based approach, generally developed for instantaneous estimates and involved physical processes of varying complexity, usually requiring auxiliary ground and atmospheric observations that cannot readily be measured through remote sensing techniques (e.g., Anderson et al. 1997;Bastiaanssen et al. 1998;Mecikalski et al. 1999;Su 2002;Allen et al. 2007).For both approaches, land surface temperature (LST) is the most essential input.Despite considerable progress having been achieved in land surface variables retrieval from remote sensing data (e.g., Liang 2004;Liu et al. 2008Liu et al. , 2014;;Trigo et al. 2008;Li et al. 2013), accuracy in retrieving LST under cloudy conditions still needs to be improved (Jiang et al. 2004;Jiménez-Muñoz et al. 2009;Yu et al. 2009;Benmecheta et al. 2013;Li et al. 2013).Therefore, remote sensing strategies that do not require remote sensed LST are worth exploring.
To derive ET from crop fields the Food and Agriculture Agency (FAO) of The United Nations recommends the following procedures (Allen et al. 1998).First, a potential daily ET from a hypothetical reference crop surface (ETo) should be computed by the recommended Penman-Monteith equation.Actual evapotranspiration rates of crops from fields well watered and under optimal agronomic conditions (ETc) are then determined by multiplying ETo with the corresponding crop coefficients (Kc), i.e., ETc = Kc ETo.For crops under water and/or other environmental stresses, the evapotranspiration rates of crops are further adjusted (ETc adj ) by incor-porating corresponding stress coefficients (Ks) and adjusted crop coefficients (Kc adj ) as ETc adj = Ks Kc adj ETo.
However, FAO Penman-Monteith equation calculation requires daily data on maximum and minimum air temperatures, vapor pressure, net radiation, and wind speed at 2 m above ground level.All of these required inputs are still not easily accessed by current remote sensing techniques.Alternatively, the use of pan evaporation (Epan) measurements to predict ETo for periods of 10 days or longer is warranted (Allen et al. 1998).The Epan can be related to the reference evapotranspiration (ETo) by empirically derived pan coefficients (Kp) as ETo = Kp Epan.The evapotranspiration process is determined by the amount of energy available to vaporize water.Solar radiation is the primary driving force of evapotranspiration from ground surfaces.Therefore, a study to relate the estimated GSR from satellite images to Epan may provide the basis of algorithms for regional ET estimation which do not require remote sensed LST as an input.
Therefore, the objectives are to establish the empirical relationships over the Taiwan region for (1) deriving reference albedos at various Sun-satellite configurations to estimate instantaneous GSR by Heliosat method using visible images from geostationary satellites and (2) estimating Epan at 10-day intervals using the retrieved GSR data with the intention to develop a different remote sensing methodology for regional ET estimation.Materials including CWB weather stations selected for validation and model development, satellite used, image preprocessing, and methodologies in developing the empirical relationships are described in section 2. Results from the proposed empirical relationships as well as error analysis discussions with ground measurements are presented in section 3. The paper closes with conclusions and recommendations in section 4.

cWb Weather Stations and Ground Measurements
Observations from 12 CWB weather stations distributed at different regions in Taiwan were used in this study (Fig. 1).SML, ALM, and JM are stations located in mountainous areas, while the remaining stations are located in plains areas along the coast.The elevations of selected stations vary from 10 -3998 m.Incident GSR and air temperature were measured using the Eppley Precision Spectral Pyranometer and platinum resistance thermometer, respectively.Measurements are given in the format of hourly values, e.g., the values at 12:00 means measurements between 11:30 and 12:30.The Epan data were measured by Class A pan and reported on daily basis.
Hourly GSR and air temperature data and daily Epan data from 1 st January 2011 to 31 st December 2013 (the study period) at these CWB stations were retrieved from databases archived and maintained by the Data Processing Branch of CWB.

Satellite and Image Preprocessing
Imageries acquired by the Multi-functional Transport Satellite (MTSAT) geostationary satellites, owned and operated by Japan Meteorological Agency (JMA) and Japanese Ministry of Land, Infrastructure, Transport and Tourism (MLITT) to observe Asia-Pacific region, have long been used by the CWB for its daily weather forecasting along with other satellites.MTSAT-2, replaced MTSAT1R for routine operation in 2010, is the second satellite in the series.MTSAT-2 sits in geostationary orbit 35800 km above the equator at around 145 degrees east longitude and scans the full-disk of the earth once every half-hour.It provides images with one visible (VIS, 0.55 -0.80 μm) and four infrared channels (IR1, 10.3 -11.3 μm; IR2, 11.5 -12.5 μm; IR3, 6.5 -7.0 μm; IR4, 3.5 -4.0 μm).At Taiwan region, ground resolution of the VIS band is approximately 4 km latitude by 2.5 km longitude.
Archived MTSAT-2 images (named as YYYY-MM-DD_HHMM.IRX.LCC.2byte.gz) of three consecutive years (2011 -2013), provided by CWB, were used in this study.In order to reduce the storage requirement and increase the computation speed, subsets of the Taiwan region were first extracted from archived images based on the latitude and longitude information from a grid data file (lat_lon_LCC.dat.gz) provided by the Meteorological Satellite Center (MSC) of CWB.During the extraction process the reflectivity of visible channel at each pixel was also converted from raw digital counts based on a calibration table provided by MSC.All retrieved subset images were then resampled using an area-weighted average scheme to produce apparent reflectivity maps, at grid spacing of 5 × 5 km, within the Taiwan Region for following studies.There were a total of 12029 daytime images available in this study.In general the daytime length varies from 5AM to 7PM in summer and from 7AM to 5PM in winter.

Estimation of Instantaneous Solar Irradiation
The Heliosat-3 approach was applied in this study.In order to produce a more representative ground albedo reference ( ground t , representing clear conditions) and thick cloud conditions ( cloud t ) for the various sun-satellite configurations, the apparent albedo (t) at every pixel from all daytime images within the studied period (2011 -2013) were first plotted by month against the co-scattering angle (W), i.e., the angle between the directions to the sun and satellite as seen from the center of that particular pixel.As discussed by Dagestad and Olseth (2007), this scheme avoided the difficulties of having few clear situations available for each month/hour combination and also reduced the changing sun-ground-satellite geometry effects with dates of a year and time of a day.Third-order polynomials were then fitted to those 4-percentile and 98-percentile values within each 10° bin to derive the required functions of ground t and cloud t with W for each pixel, respectively, as shown in Eqs. ( 1a) and (1b).
where A s and B s are fitted coefficients.Once the relationships have been derived values for ground t and cloud t could easily be computed for a given location at specific W calculated for a particular date and time.These values were then used as the lower and upper limits in cloud index (n) computation using Eq. ( 2).
Ideally, the n values should be within the range from 0 -1, representing the cloud cover fractions.Because ground t and cloud t were calculated from functions fitted to the 4-percentile and 98-percentile reflectance values [Eq.( 1)], the value of n can become negative under very clear conditions (i.e., when t < ground t ).On the other hand, the value of n can also be greater than 1 under heavy overcast conditions (i.e., ).A clear sky index (k) was further calculated from the computed cloud index (n) and the empirical relationships proposed by Rigollier and Wald (1998) This formulation ensures that instantaneous global irradiance (G) equals clear sky global irradiance (G clear ) when n = 0 and not less than 5% of G clear even under the thickest clouds (n > 1.1).The k = 1 -n relation assumes that atmospheric absorptivity does not change from clear to overcast conditions.When the atmosphere is clearer than the "reference atmosphere" for which the clear sky model applies, the negative cloud index would give a global irradiance higher than the clear sky model.The G was computed through multiplying G clear by k as shown in Eq. ( 4), The clear sky global irradiance (G clear ) is modeled as Cano et al. (1986) by Eq. ( 5).
where a x is the atmospheric transmittance, I 0 is solar constant, d r is the correction for distance between the Sun and Earth, θ is the solar elevation angle.In this study, a x was set as 0.7 (Cano et al. 1986), I 0 was set as 1367 W m -2 (Fröhlich and Brusa 1981), d r and θ were computed based on the formulas given by Michalsky (1988).Compared to the typical range (0.66 -0.68) of a x used by Lai et al. (2010) for clear days, the G clear calculated using Eq. ( 5) may differ at most by 6%.

Development of Empirical relations for Epan Estimation
The performance of six radiation-based empirical models commonly used for ETo estimation (Table 1) were first evaluated by comparing with Epan data measured over each 10-day period at 11 selected weather stations (Epan data were not available for station TN).To compare the daily global irradiance (Rs) performance estimated from different sources, integrated pyranometer measurement values at weather stations and those derived from satellite images were used, separately, in computation for the 10-day averaged daily evaporation values.The latent heat vaporization (λ), slope of saturation vapor pressure curve (Δ), psychro-metric constant (γ), and atmospheric pressure (P) quantities were calculated using formulas given by Allen et al. (1998).
where T and z are 10-day averaged daily temperatures (in °C) and elevation (in meter) above sea level at corresponding station, respectively.There were a total of 1138 observations from all stations over the entire study period (2011 -2013).To develop the empirical formula for Epan mapping all the available 10-day averaged daily Epan measurements at each selected CWB station were first arranged chronologically and then split in two parts based on the odd or even number in its sequencing.Using the split-half method the model was developed on one half and tested on the other and vice versa.As explained later in section 3.6, Epan data from mountain stations (i.e., SML, ALM, and JM) were left out from the training dataset in developing the model, but were kept in the validation test in order to check the developed model performance.The odd numbered half contained 439 data from plain stations and 135 data from mountain stations, while the even numbered half had 432 and 132 data from stations located at plain and mountain areas, respectively.

Error Analysis
In order to evaluate the performance of the derived ( ) ground t W and ( ) cloud t W functions, distribution statistics for the computed cloud indices at pixels corresponding to selected CWB station locations within the study period were computed under contrasting timeframes, such as morning (from Sunrise to 12PM) and afternoon (from 12PM to Sunset) or four seasons (Spring: March to May, Summer: June to August, Autumn: September to November, and Winter: December to February).The rationality and consistency of the computed cloud indices were used to as the basis for performance judgements.
The accuracy in estimating the daily global radiation and 10-day averaged daily Epan were evaluated using statistical indices including: mean bias deviations and root mean square differences (MBD and RMSD) as well as their relative contributions in percent of the mean observed values (MBD% and RMSD%).
where N is the number of days compared, X i and X i l are the estimated and observed values, respectively.For the GSR evaluation the observed data are the daily global irradiances derived by summing the hourly reported data from each selected station while the estimated data are those derived by integrating the retrieved instantaneous irradiations at pixels corresponding to the CWB station location with time.In order to explore the seasonal effects the evaluations were conducted based on the seasons (Spring: March to May, Summer: June to August, Autumn: September to November, and Winter: December to February).For the ETo model performance comparison, the observed data were pan measurements at each station while the estimated data were computed from selected empirical models.For the Epan validation the observed and estimated data were those from the validation dataset described above.

changes of co-Scattering Angle (}) with Days and time
Ground reflectivity changes with the sun-ground-satellite geometry can be represented using the co-scattering angle.The geometry changes with the time in a day and with seasons, and also changes somewhat even for the same time of day during a month.Figure 2 showed co-scattering angle variations at Chiayi station.Other locations in Taiwan showed similar changes.The co-scattering angle periods less than 60 degrees are from early morning to about 3PM and 4PM on winter and summer solstice days, respectively.In the late afternoon the co-scattering angles are generally larger than 60 degrees.

changes in Apparent reflectivity with co-Scattering Angle (})
The lower end of apparent reflectivity ( ground t ) showed a decreasing trend with increasing co-scattering angle (}) because the amount of shadows seen on the ground increased with } (Fig. 3).The reflectivity would reach peak value when } approached zero because no shadows should be seen from the satellite, which is also called the opposition effect (Shkuratov and Helfenstein 2001).Not considering the opposition effect, a fixed ground t value, e.g., a mean reflectivity of all } bins, would underestimate the global irradiance when } changes toward zero because a higher apparent reflectivity could give a higher cloud index using Eq. ( 2).This would cause global irradiance underestimation when the } values were small, e.g., in the morning (Fig. 2).Conversely, the global irradiance would be overestimated in late afternoon when } changes toward the higher end.As shown in Fig. 3 the third-order polynomial fitted to the 4-percentile value within each 10° bin of W for each month could represent the lower bound of reflectivity.Therefore, the errors in ground t estimation due to opposition effects could be minimized using the polynomial fitting.

Model
In theory thick clouds were close to a Lambertian reflector and should give rather homogeneous back scattering, which would produce similar reflectivity for all W bins at the higher end.However, a rapid decrease in apparent reflectivity at the higher end ( cloud t ) with the co-scattering angle (}) was observed (Fig. 3).This phenomenon might result from the apparent albedo at visible band computing scheme by CWB, which is presumably the ratio of reflected solar radiance received by the satellite to Solar Constant (1367 W m -2 ).As shown in Fig. 2a large } values generally occurred when the solar elevation angle was low.Therefore, the rapid decline in apparent reflectivity at the higher end with increasing } angles was due to the rapid decrease in incident solar irradiance following the cosine law.The cloud index would be significantly underestimated if a constant value [e.g., 0.8 by Dagestad and Olseth (2007)] was used for all } angles to expedite the computation speed.As shown in Fig. 3, the thirdorder polynomial fitted to the 98-percentile value within each 10° bin of W for each month could represent the higher bound of reflectivity.Therefore, possible errors resulting from cloud t estimation in cloud index computation using Eq. ( 2) were also minimized by adopting the polynomial fitting scheme.

cloud Index Statistics
The distributions of all computable n values in the morning and afternoon periods are shown in Fig. 4. The n values distribution is more positively skewed and with sharper peak at the lower end, i.e., having smaller n values in the morning than in the afternoon.As shown in the figure the medians in the morning were less than those in the afternoon, which indicated that the percent of cloud cover was less in the morning than in the afternoon.Negative n values (very clear conditions) occurred most often before 9AM while n values greater than 1 (heavy overcast conditions) occurred most often in the period between 11AM to 3PM (Fig. 5).This is in agreement with Kerns et al. (2010) that skies are generally clearer in the morning than in the afternoon.On clear mornings when the sun is able to heat the Earth's surface the conditions for cumulus formation are favorable.After reaching a maximum in the afternoon the cumulus formation activity decreases and the clouds dissolve in the late afternoon or early evening.
The distribution patterns of the n values also varied among different seasons and regions.For stations situated along the east coast of Taiwan (IL, HL, and TD), the n values in winter and spring were distinctively higher than those in summer and autumn, and also higher than other regions (Table 2).For stations situated in the northern region of Taiwan (TP and HC), n values in summer and autumn were relatively smaller than those in winter and spring.For stations situated in the western and southern regions of Taiwan (TC, CY, TN, and KS), there were no clear distinction in n values between seasons but the n values in winter and spring were distinctively smaller than other regions.For stations situated in mountainous areas (SML, ALM, and JM), the distinctions in n values between seasons were not clear.
The above described spatial and temporal variations in n values are explainable by the circulation of the East Asia monsoon and topographical characteristics in Taiwan.The northeast monsoon during the winter half-year and the southwest monsoon from the equatorial ocean during the summer half-year are the two major prevailing monsoons.The Central Mountain Range, which has over 100 peaks higher than 3000 m, runs from north to south and divides the island into eastern hills and western plains.Therefore, during the northeast monsoon prevailing seasons (winter and spring), the cloud cover in the eastern and northern regions of Taiwan are high resulting from orographic lifting on the windward side, while the western and southern regions have low cloud cover from being located on the leeward sides of mountainous terrain.During the summer local vertical mixing enhanced by stronger solar heating in conjunction with humid air from the warm ocean favors the formation of broken convective cumulus clouds on land, which results in a similar cloud cover over the southwestern plain area.Cloud formation over the mountainous area results from the combined effects of solar heating, orographic lifting, blocking of prevailing winds, and their interactions with diurnal wind cycle.Therefore, the amount of clouds is more evenly distributed over the mountainous area and is higher than those over the western coastal plains.

comparisons between Measured and Estimated GSr
Estimated from Eq. ( 5) the maximum possible (assuming n = 0 all day long) daily GSR in the Taiwan region was 27.0, 27.1, 24.6, and 20.8 MJ m -2 day -1 for spring, summer, autumn, and winter, respectively, while the corresponding  minimum (assuming n = 1 all day long) were 1.0, 1.2, 0.8, and 0.7 MJ m -2 day -1 , respectively.As shown in Table 3 the minimums and maximums for the measured daily GSR exceeded the possible ranges (marked in bold) in many stations, which is an indication of the lack of proper pyranometer calibration.At the JM station the highest median value occurred in autumn.This was inconsistent with the seasonal changes at other stations where the highest median value occurred in summer instead.Compared with the measured values, the distribution patterns of daily GSR retrieved from MTSAT images by this study were much more reasonable either in the ranges of solar irradiation received or the seasonal consistency (Table 4).However, linear relations with high R 2 values existed between the retrieved and the observed values (Table 5).The somewhat smaller R 2 values for those mountainous stations were probably the results of large spatial variations in cloud cover due to complex interactions between the terrain and atmospheric circulation.This indicated that the pyranometer linear responses at each station to the intensity of incident solar radiation were still good.However, the relations were not close to 1:1 lines.The positive intercepts indicated that the solar irradiance observations were systematically underestimated, which is in agreement with the findings of Lin et al. (2004).The slopes of values less than 1 indicated that the observed GSR were underestimated under overcast condition, but overestimated for clear conditions, which may be the causes of those extraneous observed minimums and maximums at many stations.
The overall relative MBD% and RMSD% averaged  across all stations and four seasons was about 5 and 15%, respectively, but the MBD% and RMSD% varied very differently among stations and seasons.Seasonally, the largest MBD% and RMSD% both occurred in spring (9.5 and 21.3% on average), while the least MBD% and RMSD% occurred autumn (-0.3% on average) and in winter (9.7% on average), respectively.Station TP had the largest MBD% (40.9%) in winter and RMSD% (36.7%) in spring.Therefore, before the pyranometers are carefully calibrated the values estimated from satellite images may be the best surrogates if accurate GSR data are required.

Performance of the Six Empirical Formulas using Measured and Estimated GSr as Inputs
The goodness of fit between the 10-day averaged ETo values computed from the six empirical formulas with corresponding Epan measurements at selected CWB stations are listed in Table 6.Judging from the overall MBD and RMSD averages and ranges, the Turc, Hargreaves, and Hansen models were considered having better performance than the other three models.The Makkink model was not selected because it had higher estimation errors for stations located in plain areas.This indicated that an empirical formula with a format similar to the Turc, Hargreaves, and Hansen models may be more suitable for estimating ETpan within the island of Taiwan.
The Turc and Hargreaves models require T and Rs as explicit inputs but the Hansen model uses only Rs as the explicit input.When comparing the MBD% and RMSD% using different Rs sources as the model inputs it is clear that using Rs retrieved from satellite images would produce smaller MBD% and RMSD% for stations located in plain areas (Fig. 6).However, in mountainous areas the retrieved Rs tend to produce larger MBD% and RMSD%.As discussed above the retrieved Rs were intended for a 5 × 5 km area and thus were not good surrogates for measured values at any particular site in the mountains where the cloud cover and amount may change very rapidly.
However, no significant benefits were observed by explicitly including T as an extra explicit input.Therefore, for input data simplicity it may be worthwhile to use Rs retrieved from satellite images as the sole input for Epan estimation provided the temperature effects can be included in computing the latent vaporization heat (λ), saturation vapor pressure curve slope (Δ), and psychrometric constant (γ).The possibility of using Rs as sole explicit input in computing Epan has another advantage because it is still difficult, if not impossible, to retrieve LST under cloudy conditions using remote sensing techniques.Table 5. Linear relations (Y = a + bX) between the estimated (Y) and observed (X) daily global radiation at selected CWB weather stations during years from 2011 -2013 and discrepancy analysis.
Note: 1: a, Y, and X are in unit of MJ m -2 day -1

Development of Empirical Formulas for Estimating Epan using Derived GSr
As discussed above a Hansen type model might be suitable for Epan estimation in the Taiwan region.To overcome the difficulties in retrieving accurate LST by remote sensing techniques, the empirical formulas provided by Chiu et al. (2014) for monthly temperature lapse rate (TLR) patterns at four geographical climatic regions in Taiwan were tested first to estimate the temperature at each studied CWB station (Fig. 1).The assigned regions considered both the influence of the prevailing winter monsoon on the TLR and the goodness of fit of the empirical equations to be used.Good 1:1 linear relations between the estimated and measured temperatures indicated that the regional arrangements and the empirical formulas selected have properly corrected the geography and seasonal change effects on the estimated temperatures (Fig. 7).
Based on the previous six selected model performance discussion, a Hansen type model was developed using the training dataset described in section 2.4.As error analysis discussed above (Fig. 6), Epan measurements at mountain stations were not closely related to the corresponding GSRs retrieved from images because of mismatching in representative areas.Therefore, data pairs in training dataset from mountain stations were left out in deriving the model.Values of λ, Δ, and γ were computed using temperature estimated by corresponding TLR formula and height above sea level.The derived empirical model for Epan estimation was listed as Eq. ( 12). .
The validation test results for the developed Epan estimation model at selected CWB stations are shown in Fig. 8.The very small MBD% (0.1%) and RMSD% (16.9%) for data pairs located in plain areas indicated that the proposed empirical model performed well in coastal plain areas.However, for data pairs located in mountainous areas, the large differences in surrounding microenvironments (e.g., percent vegetation cover, air humidity, ground surface roughness, etc.) between mountainous and plain areas might have resulted in the Epan overestimation and thus gave a very large MBD% (30.2%) and the corresponding RMSD% (34.5%).However, with proper site-specific pan coefficient (Kp), the estimated Epan may be converted into ETo required for further applications.

Spatio-temporal Variation of cloud Index, GSr, and Epan in taiwan region
From the analyzed 12029 MTSAT-2 daytime images within the study period (2011 -2013), the monthly clear skies frequency occurrence (n ≤ 0.2) (Fig. 9) and cloudy skies (n ≥ 0.8) (Fig. 10), and monthly averaged daily GSR (Fig. 11) and Epan (Fig. 12) atlases for the Taiwan region were compiled.As shown in Fig. 9 the clear skies time over land was less than that over the ocean during the summer half-year (May to October).This phenomenon is particularly apparent in the summer season (June to August) because of stronger atmospheric stability over the colder ocean surface.Resulting from the prevailing southwest monsoon during the summer season, the clear skies time decreased with increasing height in topography was also vary apparent in the western half of Taiwan in this season.The prevailing northeast monsoon during the winter season (December to February) significantly decreased the clear skies time in the north and northeastern regions.The coastal plain area in the southwest region experienced more clear skies all year long than any other places in Taiwan.The chance for clear skies could reach 60% in July.Very cloudy skies occurred most often over land in the north and northeastern regions and the surrounding ocean during the winter half-year resulting from prevailing northeast monsoon and in the mountainous area due to orographic lifting (Fig. 10).In April and May, frontal systems might cause very cloudy skies island wide.In August, strong convective lifting on land in conjunction with humid air from the ocean from the southwest monsoon might be the main causes for increased chance of very cloudy skies over the southwest region of Taiwan.Mountains on the windward side of the northeast monsoon experienced more overcast skies than any other places in Taiwan.The chance for cloudy skies could reach 30% in winter.In July and September, the chance for experiencing overcast skies was less than 10% island wide except for a few places in the mountains.The region near Yuan Yang Lake (24°34'37"N, 121°24'09"E) was the place that constantly experienced overcast skies more than 10% year round in Taiwan.
As affected by the spatial and temporal variations in cloud cover disclosed above, the amounts of incident GSR were larger on the western half of Taiwan than on the eastern half, and the mountainous area received less solar radiation than the plains area (Fig. 11).The coastal plain in the southwestern region received the largest amount of solar radiation every month.The lowest monthly averaged daily solar radiation was 9 MJ m -2 day -1 in December while the highest was 21 MJ m -2 day -1 in June and July.The total annual solar irradiation in Taiwan ranged from 1120 -1725 kWh m -2 year -1 .About 6% of Taiwan received solar irradiation less than 1200 kWh m -2 year -1 , 71% between 1200 and 1500 kWh m -2 year -1 and 23% more than 1500 kWh m -2 year -1 .The averaged total annual solar irradiation received over the island within the study period was about 4.7 × 10 13 kWh year -1 , which was 200 times more than the total annual electric power generated.The daily Epan values were also larger in the western half of Taiwan than those in the eastern half, and the values decreased with increasing heights in mountainous areas (Fig. 12).The coastal plain in the southwestern region showed the largest value in every month where the lowest monthly averaged daily Epan was 2.0 mm day -1 in January while the highest was 5.0 mm day -1 in June and July.The total annual Epan in Taiwan ranged from 750 -1300 mm year -1 .About 11% of Taiwan had annual Epan less than 850 mm year -1 , 42% between 850 and 1000 mm year -1 , 38% between 1000 and 1200 mm year -1 , and 9% more than 1200 mm year -1 .The total annual potential evaporative loss over the island was about 3.4 × 10 10 m 3 year -1 , which was about 36% of the total annual precipitation received.

concluSIonS
The empirical relationships for the reference albedos to be used in computing cloud indices for Heliosat method application in Taiwan were developed using three years' MT-SAT-2 daytime images.The computed cloud indices were reasonable and consistent with knowledge regarding the spatial and temporal variations in cloud cover over Taiwan.The estimated daily solar irradiation derived by integrating instantaneous solar radiation retrieved from MTSAT-2 images was evaluated against the measured data from twelve CWB stations during years from 2011 -2013.The overall MBD% and RMSD% for daily GSR retrieval were about 5 and 15%, respectively.The results clearly showed that the MTSAT-2 data can be used for mapping global solar irradiation over Taiwan with a ground resolution of 5 × 5 km.This accuracy is a great achievement when considering the effects of very sparse radiometric network in mountainous area on producing a high quality solar radiation atlas.
A Hansen type empirical model for estimating 10-day averaged daily Epan was developed and validated.The empirical model used solar radiation retrieved from MTSAT-2 images as input in conjunction with latent vaporization heat (λ), saturation vapor pressure curve slope (Δ), and psychrometric constant (γ) modified by surface temperatures estimated from TLR relationships.The model performed quite well in coastal plain areas, but overestimated in mountainous areas.The results showed that the solar radiation retrieved from MTSAT-2 images can be used for mapping Epan over Taiwan with a ground resolution of 5 × 5 km.The results provided insights into the choice of the most suitable empirical equation when mapping ET was prevented either by lacking adequate ground measurement density required by FAO recommended procedures or the unavailability of good quality LST data retrieved by remote sensing techniques under clouds.
The achievements in this study pave the road to produce a more reliable solar irradiation atlas over Taiwan since dates when geostationary meteorological satellites images were available and also presents ground breaking work to provide more reliable ET estimation at regional scale for better water resources management and planning.Though coefficients to convert Epan into ET have been proposed for several pan deployment setups and major crops (Allen et al. 1998), to further convert from Epan maps into ET maps, site-specific pan coefficient (Kp) for converting Epan into ETo as well as coefficients similar to crop coefficients (Kc) for converting ETo to ET are still needed for every major land cover type in Taiwan.ET flux measurements using the eddy correlation method over extensive areas may serve as ground truth data for testing the accuracy of estimated ET.Studies on estimating proper site-specific Kp and Kc coefficients as outlined above have already began.

Fig. 1 .
Fig. 1.Topography and locations of selected weather stations distributed within the study region.The red lines represent borderlines of geographical climatic regions for land surface temperature (LST) estimation using the empirical formulas of temperature lapse rate (TLR) provide by Chiu et al. (2014).NWr: Northwest region, WWr: Windward region, LWr: Leeward region, and Er: East region.(Color online only)

Fig. 2 .
Fig. 2. Changes of co-scattering angles with (a) time and (b) solar elevation angle on days of equinoxes and solstices at Chiayi, Taiwan.

Fig. 3 .
Fig. 3. Changes in apparent reflectivity with co-scattering angle at Chiayi in (a) March, (b) June, (c) September, and (d) December.Lines at top and bottom are the fitted third order polynomial curves used for cloud t and ground t

Fig. 5 .
Fig. 5. Counts of the computed cloud indices had values (a) less than zero and (b) greater than one at selected CWB weather stations during years from 2011 -2013.

Fig. 7 .
Fig. 7. Comparisons of observed air temperatures and temperatures estimated using the empirical formulas of temperature lapse rate (TLR) provide by Chiu et al. (2014).Blue circle data pairs are for stations located at plain area, while red cross pairs are for mountainous stations.(Color online only)

Table 2 .
Descriptive statistics of the computed cloud indices, separated into four seasons (Spring: March to May, Summer: June to August, Autumn: September to November, and Winter: December to February), at selected CWB weather stations during years from 2011 -2013.

Table 3 .
Descriptive statistics of daily GSR (MJ m -2 day -1 ) observed at selected CWB weather stations during years from 2011 -2013.

Table 4 .
Descriptive statistics for daily GSR (MJ m -2 day -1 ) retrieved from images at pixels corresponding to locations from selected CWB weather stations during years from 2011 -2013.