The Use of Visible Geostationary Operational Meteorological Satellite Imagery in Mapping the Water Balance over Puerto Rico for Water Resource Management

A solar insolation satellite remote sensing product for Puerto Rico, the US Virgin Islands (USVI), Dominican Republic, Haiti, Jamaica, and Cuba became available in 2009 through a collaboration between the University of Puerto Rico-Mayagüez Campus and the University of Alabama in Huntsville. Solar insolation data are available at 1 km resolution for Puerto Rico and the USVI and 2 km resolution for the other islands, as derived from 500 m resolution GOES-16 visible imagery. The insolation data demonstrate the powerful utility of satellite-derived fields for water resource applications, specifically the routine production of potential and reference evapotranspiration. This chapter describes the theoretical background and technical approach for estimating components of the daily water and energy balance in Puerto Rico. Useful information can be obtained from the model, which benefits disaster and emergency management, agriculture, human health, ecology, coastal water management, and renewable energy development at the island scale.


Introduction
Estimates of incoming solar radiation (also known as "insolation") have been made from geostationary satellite data for many years, since the early to mid-1970s [1]. Related to the present effort, Geostationary Operational Environmental Satellite (GOES) visible channel ($0.64 μm) data have been processed within a scalable and flexible insolation model, which is well documented and described in detail below. For ongoing water management support over Puerto Rico and the broader Caribbean, the Diak-Gautier insolation model [2] has been specifically structured to provide daily integrated, gridded solar insolation at 1-2 km spatial resolution. The insolation model has been rigorously tested and validated and operates on GOES imagery from GOES-4/-5 through the present day GOES-16/-17. Geostationary satellites are optimal for providing spatially and temporally continuous fields across all regions in their AE55°latitude field of view, which as noted is a significant advantage over the use of only ground-based instrumentation. The use of a satellitebased insolation algorithm also ensures that a consistent algorithm is applied across an entire region, one which relies on data from only one instrument, specifically, the GOES Imager.
Over Puerto Rico (PR) and the Caribbean, as well as in other subtropical and tropical regions, evapotranspiration (ET) is a critical variable for water management, both in hydrologic flow simulations involving potential ET (PET) and water allocation and agricultural water use involving reference ET (RET or ET o ). Importantly, solar insolation is a large, yet often unknown, determinant for temporal variation in PET and RET. Solar insolation is a primary determinant of spatial variation, particularly in areas with heterogeneous cloud cover, as common to subtropical and tropical regions where small cumulus clouds dominate the regional cloud climatology.
For an ET product to be desirable, it must be spatially continuous, rather than consisting of only point values derived from local weather station networks. Thus, mapping of ET is greatly facilitated by satellite-derived estimates that contain the actual spatial variability and distribution of solar insolation. Prior to 2009, regions across Puerto Rico and the Caribbean did not had access to a consistent, spatially continuous method of computing RET and PET. The original motivation for development of the Geostationary Operational Environmental Satellite-Water and Energy Balance (GOESWEB) model was to develop a robust insolation calibration framework coupled to a satellite-based insolation model, to provide a key radiative dataset that can grow over time toward 10-year and longer timeframes, thus forming an ET climatology that can be extended indefinitely.
The GOES-based insolation datasets are used in conjunction with other information, including net radiation (R n ), air temperature, relative humidity, wind speed, and land cover information, in the formulation of daily, 1-and 2-km estimates of RET across the Caribbean. RET is valuable for farm-and city-based water management, as well as irrigation scheduling; PET can be used as input into surface and groundwater hydrological models, whereas the solar insolation data themselves may be used as data input in certain ecosystem models.

GOES solar insolation data
The use of geostationary satellite visible data has been used for estimating solar insolation for over 30 years. The main methods used for such estimation range from statistical-empirical relationships, such as [3], to varying complex physical models [2,[4][5][6][7][8][9][10][11][12]. Studies such as [13,14] proved the utility and feasibility of satelliteestimated solar insolation methods, demonstrating that fairly accurate results can be produced from such models; hourly insolation estimates obtained from the most current models are within 5-10% of ground-based pyranometer data, during clearsky conditions (15-30% for all sky conditions), while daily estimates are found to be within 10-15% [15]. Studies by [16,17] have further highlighted the overall utility of these methods.
The main advantages of using satellite-estimated insolation, over those collected by pyranometer networks, include wide-area spatial coverage, high spatial resolution (1-2 km), and the ability to produce useful data in remote, inaccessible, or in potentially hazardous areas, over large water bodies and oceans (e.g., [18]), and in locations where the installation of a ground-based pyranometer network is prohibitive. As an alternative to the methods used in this study, [19,20] describe the use of the Global Energy and Water Cycle Experiment (GEWEX) Surface Radiation Budget (SRB) downward solar flux [11], as used within the North American Land Data Assimilation project. Error statistics for the SRB product are comparable to those shown in [21], as used in this study, yet SRB resolutions are at best 0.5°and 3 hourly [22].
Related to the PRWEB applications to be developed here, the solar insolation product is derived from National Oceanic and Atmospheric Administration (NOAA) GOES-East satellite visible (0.64 μm) imagery. These data were processed using [4] methods to produce daily integrated solar insolation throughout Puerto Rico at 1-km horizontal spatial resolution. This 1-km resolution is chosen as it provides solar insolation observations between cumulus clouds, which comprise a significant component of the cloud climatology in subtropical regions.

Details of the GOES solar insolation model
The GOES solar insolation model is developed by [4], which was later modified by [2] and updated by [23], and most recently by [24], which is the 2017 version of the solar insolation model employed in this study. This model will be referred to as the "GD" model from this point forward and employs a simple physical model that represents cloud and atmosphere radiative processes. The GD model was shown to perform even better than more complex solar insolation methods over a variety of land-surface and climatic conditions [5,17,18,23,26]. When comparing with pyranometer data, these prior studies list root mean square errors in hourly and daily insolation estimates as a percentage of the mean pyranometer observed value, which range from 17-28% to 9-10%, respectively. In [24,25] the higher magnitudes of these errors were reported ($28 and $10%, respectively) in a study over northern central Florida using GOES-12 data. However, the GD model has been proven to be valuable in operational use of near-real-time, regional-, and continental-scale insolation estimates for several main applications, including land-surface carbon and water flux assessments [27][28][29], the generation of agricultural forecasting products [30,31], and subsurface hydrologic modeling.
The GD model is based on conservation of radiant energy in the Earthatmosphere column, with two modes for estimating solar insolation received at Earth's surface: (1) clear and (2) cloudy conditions. These modes are determined based on satellite-derived, visible channel surface albedo data. A reference albedo grid representative of clear-sky conditions per satellite pixel is developed within the GD algorithm, which captures the temporal changes in land-surface characteristics over time and season. This running 2-week minimum of this albedo data, reassessed at solar noon daily, is stored for each GOES satellite visible data pixel. This approach is considered representative of the true land-surface albedo, which is more accurate than using the daily estimated value as the latter may be corrupted by high albedo values when even low-cloud amounts are present during a given day. Note that this minimum albedo is wavelength-specific, is unique to the GOES Imager visible sensor (which includes some near-infrared contribution), and is not a true surface albedo.
As the GD algorithm runs across a series of GOES images per day, the digital brightness at each image pixel is compared to that of the stored clear-sky reference 2-week minimum albedo for that pixel. If the brightness exceeds that threshold, the pixel is deemed partly or completely cloudy. Based on this determination per GOES pixel, either the clear or cloudy model of atmospheric radiation processes (within the GD model) is used to calculate surface solar insolation received. Both clear and cloudy models incorporate parameterizations for Raleigh scattering, ozone absorption, and water vapor absorption within the atmospheric column, using simple bulk relationships, such as fixed ozone and aerosol contents. This rough parameterization works because these produce secondary sources of error to the instantaneous surface solar insolation. The cloudy GD component estimates a cloud-top albedo and separately accounts for atmospheric effects above and below the cloud.
For the water vapor absorption parameterization, a fixed, approximate annual median value of precipitable water (PW) of 3.0 cm was used, which is considered appropriate for Puerto Rico. This annual median value helps to estimate atmospheric column-integrated PW during the initial processing. [PW is defined as the amount of water that would precipitate out of a vertical column of the atmosphere if all the water vapor were condensed into liquid]. PW data are used to calculate the slantwise path and subsequently the absorption coefficients [4]. Real-time PW data from numerical forecast model output may also be used in the GD model, versus setting a constant value.

GOES data processing and quality control
The GOES-East series of satellites (the most recent additions being GOES-13 and -16) are in geostationary orbit above the Earth's equator at À75°W, which provides continuous, 5-15-minute resolution observations in visible and infrared radiation channels at high spatial (500 m to 1 km). GOES data are thus ideal for high-resolution estimates of solar insolation as used in GOESWEB, to be described below. Although the GOES visible sensors have a nadir (the point directly below the satellite) spatial resolution of 1 km (GOES-13 and prior) or 500 m (GOES-16), this resolution decreases the further from nadir the instrument scans: for Puerto Rico, the highest resolution attainable is about 1.25 km and 525 m, respectively, for GOES-13 and -16. All solar insolation data used for this study were provided at 1km resolution. A simple method for computing sunrise and sunset times per pixel across the domain was used, as a means of determining daytime conditions. Potential significant GOES data issues that may impact the error in the solar insolation product include (1) sensor degradation with time and (2) sun glint effects. The effects of the latter are small. In general, GOES satellite data are available on a continual basis with high reliability (>99%). As an example, Figure 1 shows the daily integrated solar radiation for October 16, 2018, for Puerto Rico, the USVI, Hispaniola, Jamaica, and Cuba.

The GOESWEB modeling framework
GOESWEB performs daily water and energy balance calculations for the island of Puerto Rico. Twenty-seven hydro-agro-climate variables are available to the public for download ( Table 1). Downloadable formats are available as images (jpg) or in comma-separated values (csv) and Matlab® formats. The variables in Table 1 are also available as monthly and annual averages or totals. Simplified versions of the algorithm have been developed for estimating reference ET on the islands of the USVI, Hispaniola, Jamaica, and Cuba.
ET o is estimated by three methods: Penman-Monteith [32,33], Priestley-Taylor [34], and Hargreaves-Samani [35]. In [36], they described the methodology used to estimate ET o in the earliest version of the algorithm. T avg , T min , and T max values were estimated from a lapse rate method developed by [37]. T d was assumed to be equal to the minimum daily T min [38]. Wind speed was assumed to be the worldwide average 2-m wind speed of 2 m(s) À1 [32]. The algorithms for Hispaniola, Jamaica, Cuba, and the USVI continue to use these simplified methods for estimating daily values of ET o .
Water and energy balances were added to the algorithm for Puerto Rico. The daily meteorological data used are described below. Table 2 summarizes the GOESWEB input data sources.
• Solar radiation i. Solar radiation (R s ) is derived from the GOES satellite using the methodology described above.
ii. The ground level, 1-km resolution R s product became available in Puerto Rico in March of 2009 and has been validated at two locations in Puerto Rico by [39].
iii. Occasionally the satellite-derived solar radiation is not available, in which case the previous days' R s values are used.
iv. Prior to GOES-16, 1 km GOES-12 and -13 visible channel 1 data were used over Puerto Rico and the USVI and 2 km data over the other islands.   iii. For the reference ET calculation, 10-m wind speeds are adjusted to 2 meters [32].
iv. If wind speed is not available, then the previous day's data are used.

Net radiation calculations
Net radiation (R n ) is estimated using the methodology described by [32] and used by [41].
where R n is net radiation, R ns is net shortwave radiation, and R nl is net long wave radiation.
where α is albedo and R s solar radiation. α is defined as 0. 23 where σ is the Stefan-Boltzmann constant, T max is maximum absolute temperate during the 24-hour period, T min is minimum absolute temperature during the 24-hour period, e a is actual vapor pressure, R s /R so is relative shortwave radiation (limited to ≤1.0), and R so calculated clear-sky radiation. Actual vapor pressure is estimated by e a ¼ 0:6108 exp where T d is dew point temperature. The calculated clear-sky radiation is estimated by where z is elevation above mean sea level and R a is extraterrestrial radiation.
where G sc is the solar constant = 0.0820 and d r is the relative distance Earth-Sun, defined as where J is Julian day (e.g., January 1 is 1 and December 31 is 365). ω 1 in Eq. (6) is solar time angle at the beginning of the period and ω 2 solar time angle at end of period, generally expressed as where φ is latitude and δ solar declination expressed as δ ¼ 0:409 sin 2π 365 J À 1:39 (9) and X is defined as and X = 0.00001 if X ≤ 0.

Reference evapotranspiration estimates
The Penman-Monteith (PM) equation is given by Eq. 1 [32], which applies specifically to a hypothetical reference crop with an assumed crop height of 0.12 m, an albedo of 0.23, a fixed surface resistance of 70 sec m À1 , and an aerodynamic resistance equal to 208/u 2 , where u 2 is wind speed at 2 m height: where Δ is the slope of the vapor pressure curve, G is soil heat flux, γ is the psychrometric constant, T is mean daily temperature at 2 m height, e s is the saturation vapor pressure, and e a is the actual vapor pressure.
where α is the Priestly-Taylor constant. Values in the literature for α range from 1.26 [34] to 1.32 [45]. In this study we use a value of α equal to 1.3. The third method used to estimate ET o is the Hargreaves-Samani ET o Equation [35] given by The value 0.0135 is a constant and 0.408 converts the result from MJ m À2 day to mm (day) À1 . In [38] they showed that this method produces comparable results with the PM method in PR.
The PM method is considered superior to the other two methods because it accounts for the major variables that control ET (R n , T, VPD and u), and the PM method has been rigorously validated [33].

Energy balance
In GOESWEB, an energy balance approach is used similar to [46]. The basic energy balance equation is given as R n is obtained from the calculation procedure presented above. Albedo, which is used in the R n calculation, is obtained from a lookup table [42], which assigns values of the parameters to 32 different land covers.
LE, H, and G are the latent, sensible, and soil heat fluxes, respectively. LE is estimated using the following Equation [47]: where ρ is mean air density, C p is specific heat, r a is aerodynamic resistance, and r s is surface resistance. G is the soil heat flux, assumed to be zero for the daily analysis. H is estimated using the following equation: The effective surface temperature is difficult to obtain from remote sensing under cloudy conditions. Therefore, T s is obtained by an implicit approach similar to that described by [48]. When Eq. (14) is expanded using Eqs. (1), (15), and (16), T s is the only unknown variable, which is obtained using the recursive root function fzero in MatLab® (http://www.mathworks.com).
The aerodynamic resistance (r a ) is calculated with the following Equation [46]: where r ao is the aerodynamic resistance under conditions of neutral atmospheric stability and r bh is the excessive resistance. r ao is expressed as where z is the virtual height at which meteorological measurements are taken. In this study z is assumed to be within the inertial sublayer and equal to 1.5(z o /0.13) [47], which is equivalent to the canopy height (h). The NDFD or WRF modelderived wind speeds at 10 m height are adjusted to the "virtual instrument height," depending on the height of the vegetation. Roughness length (z o ) and the zeroplane displacement (z disp ) are derived from a lookup table for various land use/ vegetation categories [42]. k is Von Karman's constant (k = 0.41). u is the wind velocity at height z. From [46], the atmospheric stability coefficient is where g is the gravitational constant and the coefficient η is taken as 5 [46]. The temperature, T o , is the average of the values of T s and T a . Other variables and parameters were previously defined. The excess resistance in Eq. (17) is given by the equation Bulk surface resistance (r s ) is estimated using the equation of [49]: where VPD is the vapor pressure deficit, C f is a calibration coefficient equal to 1 for root depth <1 m and 5 for root depth >1 m, and θ FC and θ WP are the volumetric soil moisture content (θ) at field capacity and wilting point, respectively. Field capacity and wilting point were obtained from regression equations of [50] based on percent sand, silt, and clay. Soil properties for sand, silt, and clay for Puerto Rico were obtained from the USDA Natural Resource Conservation Service (NRCS) and Soil Survey Geographic (SSURGO) database.

Water balance
The water balance is estimated from the equation where SMD1 and SMD2 are the depths of soil moisture in the root zone (Rdepth) at times 1 and 2, respectively. In GOESWEB the time step is 1 day. Precip is rainfall, RO is surface runoff, and DP is deep percolation below the root zone. The daily ET a is obtained by converting LE to an equivalent depth of water by dividing by the latent heat of vaporization (2.45 MJ kg À1 ). Root depths for various land use/vegetation categories are obtained from [42] lookup table. Twenty-four-hour rainfall is obtained from NOAA's Advanced Hydrologic Prediction Service (AHPS). In PR, AHPS rainfall is bias-corrected radar rainfall using rain gauge data.
Surface runoff is estimated using the curve number (CN) method of the NRCS [51]: where S is the maximum potential difference between rainfall and runoff at the moment of rainfall initiation and CN is a proportion of rainfall converted to runoff, adjusted for antecedent rainfall conditions. CN values were derived for Puerto Rico using the method described by [51], based on land use, hydrologic soil group, and antecedent rainfall conditions. To estimate DP, the following procedure is followed: SMD2 i = Precip -ET a -RO + SMD1. If the value of SMD2 i is larger than the depth of water in the soil profile at field capacity (FCD), then DP = SMD2 i -FCD, and the value of SMD2 is equal to FCD. If SMD2 i < FCD, then DP = 0 and SMD2 = SMD2 i .

GOESWEB model accuracy and validation
In this section accuracy and validation data are presented for remotely sensed solar radiation, RET, soil moisture, and stream flow. Solar radiation is a critically important variable in the estimation of ET. Figures 2 and 3 show comparisons of the daily integrated solar radiation at the University of Puerto Rico (UPR) Fortuna Agricultural Experimental near Juana Diaz, PR, and the UPR-Mayaguez Campus (UPRM) in Mayaguez, PR, respectively [39]. The figures show a high degree of correlation between the remote sensing solar radiation and the measured solar radiation. The coefficients of determination (r 2 ) for the UPRM and experimental station data were 0.88 and 0.83, respectively. (From [39]).  Comparison of remote sensing and pyranometer-measured daily integrated solar radiation at the UPR Fortuna Agricultural Experiment Station, near Juana Diaz, PR (From [39]). The r 2 value for this comparison is 0.88.

11
data set. For this comparison, the coefficient of determination (r 2 ) was 0.31. The average GOESWEB and weather station ET o were 4.6 mm and 4.14 mm, respectively, and the average calculated error was 11.2%. It should be noted that the weather station at this location does not comply with the required "reference conditions" for computing ET o . Reference conditions refer to a grass-type vegetation with an approximate height of 0.12 m, an albedo of 0.23, and a fixed surface resistance of 70 sec m À1 , receiving adequate water. The climate of southern PR is semiarid, and there are frequent times when there was no vegetation at all on the ground surrounding the weather station. Figure 5 shows a time series comparison of soil moisture from GOESWEB and soil moisture from a weather station located at the UPR Fortuna Agricultural Experiment Station. The weather station soil moisture is an average of five sensors positioned at depths of 0.0508 m (2 in.), 0.1016 m (4 in.), 0.2032 m (8 in.), 0.508 m (20 in.), and 1.016 m (40 in.). Immediately after rainfalls the weather station soil moisture tended to rise to higher soil moisture values than the soil moisture from the model. It is important to know that maximum soil moisture values in GOESWEB are limited to the field capacity, as excess water is routed below the root zone as deep percolation. Furthermore, the sensor soil moisture represent a single point (approximately 1 m 2 ), whereas the model represents an  area of 1 km 2 (1,000, 000 m 2 ), and therefore, complete agreement between the two methods would not be expected. Figure 6 compares the monthly stream flow values for two watersheds in southwest Puerto Rico. Observed stream flow values were obtained from the US Geological Survey (USGS). The results are presented as a depth of water in millimeters (i.e., monthly stream volume/watershed area). The total stream flow for the model was assumed to be the surface runoff plus the deep percolation (or aquifer recharge). The latter term represents the stream base flow. To obtain the monthly value of stream flow in the model, the surface runoff and deep percolation were averaged for every 1 km 2 pixel within the watershed. The model does a reasonably good job of simulating monthly stream flow. Figure 7 shows an example of selected water and energy balance components for Puerto Rico on October 16, 2018. Rainfall, surface runoff, percolation below the    and H fluxes, which sum to the Rn (i.e., Eq. (14)). ETa is the LE flux divided the latent heat of vaporization constant equal to 2.45 MJ kg À1 . Figure 9 shows the rainfall during the week of September 17, 2017, the same week Hurricane Maria occurred. The maximum rainfall for the week was nearly 1300 mm (51 in.) in southeast Puerto Rico. The rainfall data were derived from rain gauge data, since the Doppler radar in Cayey, PR, was severely damaged during the hurricane. The National Weather Service (NWS) combined the gauge rainfall for September 20 and 21. The maximum rainfall during the 2-day period was 950 mm  Rainfall over Puerto Rico on September 20, the day that Hurricane Maria made landfall on Puerto Rico. The gauge rainfall reported by the NWS was for the 20th and 21st; therefore the rainfall for September 20 was assumed to be half the amount.    (37.5 in.) in southeast Puerto Rico. To simulate the daily hydrology, rainfall was evenly divided between the 2 days. Figures 9 and 10 show the rainfall and surface runoff for September 20, respectively. Note that the surface runoff is almost identical to the rainfall, as seen in Figure 11. Nearly 100% of the rainfall was converted to surface runoff because the soils were already saturated the day before Hurricane Maria arrived (September 19), as shown in Figure 12.

High-resolution ETo products across the Caribbean
GOESWEB provides daily values of ET o for Puerto Rico, the USVI, Hispaniola, Jamaica, and Cuba. As an example, the ET o for each of the islands for October 16,2018, are presented in Figure 13. In the study by [52], they describe a web-based method for determining irrigation requirements using the GOESWEB ET o maps.

Conclusions
The above study demonstrates the operational utility of incorporating spatially continuous, high spatial resolution (1 km) GOES-16-derived solar insolation, using the model described by [24], into the water balance model GOESWEB, to then estimate the complete water budget. In this demonstration, applications of water balance were performed over the US territory of Puerto Rico, a subtropical location that is very sensitive to high rates of ET, relative to various crop types and vegetation characteristics, and that also receives high amounts of rainfall. High rainfall causes significant runoff, for which the GOESWEB water balance model can help identify related to actual rainfall events. Expanding GOESWEB to other island regions would be a future avenue for the research and algorithm development activities described here.

B (unitless)
Bowen ratio C f (dimensionless) calibration coefficient CN (dimensionless) curve number, proportion of rainfall converted to runoff DP (mm) deep percolation or the soil water that passes below the root zone d r (dimensionless) relative distance Earth-Sun ET a (mm) actual evapotranspiration ET o (mm) reference evapotranspiration e a [kPa] actual vapor pressure e s (kPa) saturated vapor pressure g (m s À2 ) gravitational constant G (MJ m À2 day À1 ) soil heat flux density G sc (MJ m À2 min À1 ) solar constant = 0.0820 h (m) canopy height H (MJ m À2 day À1 ) sensible heat flux J Julian day (e.g., January 1 is 1 and December 31 is 365) k (dimensionless) Von Karman's constant (0.41) K c,eff (unitless) effective crop coefficient K s (unitless) water stress coefficient LE (MJ m À2 day À1 ) latent heat flux PAR (micromoles m À2 sec À1 ) photosynthetically active radiation Precip (mm) rainfall precipitation Precip eff (mm) effective rainfall R a (MJ m À2 day À1 ) extraterrestrial radiation r a (s m À1 ) aerodynamic resistance r ao (s m À1 ) aerodynamic resistance under conditions of neutral atmospheric stability r bh (s m À1 ) excess resistance R depth (mm) root depth RH (%) relative humidity R n (MJ m À2 day À1 ) net radiation R nl (MJ m À2 day À1 ) net long wave radiation R ns (MJ m À2 day À1 ) net shortwave radiation RO (mm) surface runoff R s (MJ m À2 day À1 ) solar radiation r s (s m À1 ) surface resistance R s /R so relative shortwave radiation R so (MJ m À2 day À1 ) clear-sky radiation S (mm) maximum potential difference between rainfall and runoff at the moment of rainfall initiation Sat (unitless) soil moisture saturation SMD1 (mm) depths of water in the soil profile beginning of the day (24 hours) SMD2 (mm) depths of water in the soil profile end of the day (24 hours) T a (°C) air temperature T avg [°C] mean daily air temperature T d (°C) dew point temperature