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

Abstract

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.

1. 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 ±55° latitude field of view, which as noted is a significant advantage over the use of only ground-based instrumentation. The use of a satellite-based 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.

2. 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 satellite-estimated 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 clear-sky 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.

2.1 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 Earth-atmosphere 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.

2.2 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 1-km 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.

Figure 1.

Daily solar insolation for (a) Puerto Rico, (b) Hispaniola, (c) Jamaica, (d) Cuba, and (e) St. Croix USVI, (f) St. Thomas, and (g) St. John (USVI) on October 16, 2018.

3. 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.

1	ET	Actual evapotranspiration (mm)
2	ETo	Reference evapotranspiration (mm); three methods
3	Tavg	Average air temperature (°C)
4	Tmin	Minimum air Temperature (°C)
5	Tmax	Maximum air Temperature (°C)
6	Td	Dew point temperature (°C)
7	Ts	Effective surface temperature (°C)
8	ea	Actual vapor pressure (kPa)
9	es	Saturated vapor pressure (kPa)
10	RH	Relative humidity (%)
11	u	Wind speed (m s^−1)
12	Rs	Solar radiation (MJ m^−2 day^−1)
13	Rn	Net radiation (MJ m^−2 day^−1)
14	PAR	Photosynthetically active radiation (micromoles m^−2 sec^−1)
15	Ks	Water stress coefficient (unitless)
16	Kc,eff	Effective crop coefficient (unitless)
17	Precip	Rainfall (mm)
18	Precipeff	Effective rainfall (mm)
19	RO	Surface runoff (mm)
20	DP	Deep percolation (mm)
21	θ	Soil moisture content (m^−3 m^−3)
22	Sat	Soil moisture saturation (unitless)
23	rs	Surface resistance (s m^−1)
24	ra	Aerodynamic resistance (s m^−1)
25	LE	Latent heat flux (MJ m^−2 day^−1)
26	H	Sensible heat flux MJ m^−2 day-1)
27	β	Bowen ratio (unitless)

Table 1.

Hydro-agro-climate variables produced daily by GOESWEB.

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)⁻¹ [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

  1. Solar radiation (R_s) is derived from the GOES satellite using the methodology described above.

  2. 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].

  3. Occasionally the satellite-derived solar radiation is not available, in which case the previous days’ R_s values are used.

  4. 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.

• Air temperatures

  1. T_avg, T_min, T_max, and T_d data were obtained from the National Digital Forecast Database (NDFD) website [40] from January 1, 2009, to December 31, 2016.

  2. Temperatures were obtained from the CariCOOS operational gridded Weather Research and Forecasting (WRF) model starting on January 1, 2017.

  3. On occasion, weather parameters from the WRF model are not available. NDFD air temperatures and wind speed are used in those cases. As a final resort, the lapse rate method of [37] is used.

• Wind speed

  1. During the period January 1, 2009, through September 30, 2015, daily average wind speed was obtained from the average of eight NDFD 3-hour values [40].

  2. From October 1, 2015, to the present, daily average wind speed was obtained from the average of 24-hourly wind speed values obtained from the Caribbean Coastal Ocean Observing System (CariCOOS) WRF.

  3. For the reference ET calculation, 10-m wind speeds are adjusted to 2 meters [32].

  4. If wind speed is not available, then the previous day’s data are used.

Model input	ETo	ETa	Source
Solar radiation	X	X	GOES
Air temperature	X	X	NOAA/NDFD, CariCOOS/WRF
Dew point temperature	X	X	NOAA/NDFD, CariCOOS/WRF
Wind speed (2-m height)	X	X	NOAA/NDFD, CariCOOS/WRF
Albedo1, root depth, roughness length, zero-plane displacement		X	ATMET, 2005
Soil texture		X	USDA/SSURGO

Table 2.

Summary table of required GOESWEB input data sources.

¹

Albedo is 0.23 for calculating ET_o.

3.1 Net radiation calculations

Net radiation (R_n) is estimated using the methodology described by [32] and used by [41].

$$R_n=R_{\mathit{ns}}+R_{\mathit{nl}}$$E1

where R_n is net radiation, R_ns is net shortwave radiation, and R_nl is net long wave radiation.

$$R_{\mathit{ns}}=\left(1-\alpha \right)R_s$$E2

where $\alpha$ is albedo and $R_s$ solar radiation. $\alpha$ is defined as 0.23 for estimating RET, and values are derived from a lookup table associated with 32 land cover classes [42] for estimating actual ET (ET_a). R_s is derived from the GOES satellite. The net long wave radiation is estimated from the equation

$${\mathrm{R}}_{\mathrm{nl}}=\mathrm{\sigma }\left[\frac{{{\mathrm{T}}_{max,\mathrm{K}}}^4-{{\mathrm{T}}_{min,\mathrm{K}}}^4}{2}\right]\left(0.34-0.14\sqrt{{\mathrm{e}}_{\mathrm{a}}}\right)\left(1.35\frac{{\mathrm{R}}_{\mathrm{s}}}{{\mathrm{R}}_{\mathrm{so}}}-0.35\right)$$E3

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

$${\mathrm{e}}_{\mathrm{a}}=0.6108exp\frac{\left(17.27{\mathrm{T}}_{\mathrm{d}}\right)}{\left({\mathrm{T}}_{\mathrm{d}}+237.3\right)}$$E4

where T_d is dew point temperature. The calculated clear-sky radiation is estimated by

$${\mathrm{R}}_{\mathrm{so}}=\left(0.75+2{10}^{-5}\mathrm{z}\right){\mathrm{R}}_{\mathrm{a}}$$E5

where z is elevation above mean sea level and R_a is extraterrestrial radiation.

$${\mathrm{R}}_{\mathrm{a}}=\frac{12\left(60\right)}{\mathrm{\pi }}{\mathrm{G}}_{\mathrm{sc}}{\mathrm{d}}_{\mathrm{r}}\left[\left({\mathrm{\omega }}_2-{\mathrm{\omega }}_1\right)sin\left(\mathrm{\phi }\right)sin\left(\mathrm{\delta }\right)+cos\left(\mathrm{\phi }\right)cos\left(\mathrm{\delta }\right)\left(sin\left({\mathrm{\omega }}_2\right)-sin\left({\mathrm{\omega }}_1\right)\right)\right]$$E6

where G_sc is the solar constant = 0.0820 and d_r is the relative distance Earth-Sun, defined as

$${\mathrm{d}}_{\mathrm{r}}=1+0.33cos\left(\frac{2\mathrm{\pi }}{365}\mathrm{J}\right)$$E7

where J is Julian day (e.g., January 1 is 1 and December 31 is 365). ω₁ in Eq. (6) is solar time angle at the beginning of the period and ω₂ solar time angle at end of period, generally expressed as

$${\omega }_s=\frac{\pi }{2}-arctan\left[\frac{-tan\left(\phi \right)tan\left(\delta \right)}{X^{0.5}}\right]$$E8

where φ is latitude and δ solar declination expressed as

$$\delta =0.409sin\left(\frac{2\pi }{365}J-1.39\right)$$E9

and X is defined as

$$X=1-{\left[tan\left(\phi \right)\right]}^2{\left[tan\left(\delta \right)\right]}^2$$E10

and X = 0.00001 if X $\le$ 0.

3.2 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⁻¹, and an aerodynamic resistance equal to 208/u₂, where u₂ is wind speed at 2 m height:

$${\mathit{ET}}_o=\frac{0.408\mathrm{∆}\left(R_n-G\right)+\gamma \left(\frac{900}{T+273}\right)u_2\left(e_s-e_a\right)}{\mathrm{∆}+\gamma \left(1+0.34u_2\right)}$$E11

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.

The second method used to estimate ET_o is the Priestly-Taylor Equation [34], a simplification of the Penman Equation [43, 44]:

$${\mathit{ET}}_o=\alpha \frac{\mathrm{∆}\left(R_n-G\right)}{\mathrm{∆}+\gamma }$$E12

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

$${\mathit{ET}}_o=0.408\left[0.0135R_s\right]\left(T+17.8\right)$$E13

The value 0.0135 is a constant and 0.408 converts the result from MJ m⁻² day to mm (day)⁻¹. 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].

3.3 Energy balance

In GOESWEB, an energy balance approach is used similar to [46]. The basic energy balance equation is given as

$$R_n-\mathit{LE}-H-G=0$$E14

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]:

$$\mathit{LE}=\frac{\rho C_p\left(e_o\left(T_s\right)-e\left(T_a\right)\right)}{\gamma \left(r_a+r_s\right)}$$E15

where $\rho$ 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:

$$H=\frac{\rho C_p\left(T_s-T_a\right)}{r_a}$$E16

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]:

$$r_a=r_{\mathit{ao}}\varphi +r_{\mathit{bh}}$$E17

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

$$r_{\mathit{ao}}=\frac{\mathit{ln}\left[\frac{z-z_{\mathit{\text{disp}}}}{z_o}\right]\mathit{ln}\left[\frac{z-z_{\mathit{\text{disp}}}}{\left(0.1\right)z_o}\right]}{k^{2}u}$$E18

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 model-derived 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 zero-plane 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

$$\varphi =\left[1-\frac{\left[\eta \left(z-z_{\mathit{\text{disp}}}\right)g\left(T_s-T_a\right)\right]}{T_o{u}^2}\right]$$E19

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

$$r_{\mathit{bh}}=\frac{4}{\left(\frac{\mathit{ku}}{\mathit{ln}\left[\left(z-z_{\mathit{\text{disp}}}\right)/z_o\right]}\right)}$$E20

Bulk surface resistance (r_s) is estimated using the equation of [49]:

$$r_s=\frac{\rho C_p\mathit{VPD}}{\mathrm{∆}\left(R_n-G\right)C_f}{\left(\frac{\theta -{\theta }_{\mathit{WP}}}{{\theta }_{\mathit{FC}}-{\theta }_{\mathit{WP}}}\right)}^{-1}$$E21

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.

3.4 Water balance

The water balance is estimated from the equation

$$\mathit{SMD}2=\mathit{\text{Precip}}-{\mathit{ET}}_a-\mathit{RO}-\mathit{DP}+\mathit{SMD}1$$E22

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⁻¹). 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]:

$$\mathit{RO}=\frac{{\left(\mathit{\text{Precip}}-0.2S\right)}^2}{\left(\mathit{\text{Precip}}+0.8S\right)}$$E23

$$S=\left[\left(\frac{25400}{\mathit{CN}}\right)-254\right]$$E24

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.

3.5 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²) for the UPRM and experimental station data were 0.88 and 0.83, respectively. (From [39]).

Figure 2.

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² value for this comparison is 0.88.

Figure 3.

Comparison of remote sensing and pyranometer-measured daily integrated solar radiation at UPRM (From [39]). The r² value for this comparison is 0.83.

Figure 4 shows a comparison of the ET_o computed by the GOESWEB algorithm and from weather station data from the UPR Fortuna Agricultural Experiment Station, near Juana Diaz, PR. The ET_o data covers the period from December 12, 2013 to April 20, 2016 (858 days). Although the vast majority of data pairs fall close to the 1:1 line, indicating close agreement between the two methods, a smaller number of data pairs fall relatively far from the 1:1 line, producing the scatter in the data set. For this comparison, the coefficient of determination (r²) 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⁻¹, 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 4.

Comparison of observed and simulated ET_o at the UPR Fortuna Agricultural Experiment Station, near Juana Diaz, PR. The data cover the period December 12, 2013 through April 20, 2016. The r² value for this comparison is 0.31.

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²), whereas the model represents an area of 1 km² (1,000, 000 m²), and therefore, complete agreement between the two methods would not be expected.

Figure 5.

Comparison of the GOESWEB and weather station soil moisture at the UPR Fortuna Agricultural Experiment Station, near Juana Diaz, PR. The data covers the period January 1, 2014–December 31, 2016. The r² value for this comparison is 0.73.

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² pixel within the watershed. The model does a reasonably good job of simulating monthly stream flow.

Figure 6.

Comparison of observed and simulated stream flow for two watersheds in southwest Puerto Rico. The data cover the 36-month period during 2010–2012. The r² value for this comparison is 0.72 for all data.

4. High-resolution products for Puerto Rico

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 root zone, soil moisture content, actual ET, Rn, LE, and H are included. Approximately 60 mm of rain fell along the northern coast of the island. High values of surface runoff occurred in the rainy area where soil textures have high clay content. High values of percolation below the root zone occurred in small areas where the soil sand content approaches 90%. The soil moisture map indicates a lobe of dry area in Salinas, Cayey, Aibonito, and Cidra. Figure 8 shows the NOAA Drought Monitor for Puerto Rico for October 18, 2018, indicating abnormally dry conditions for Cayey, Aibonito, Cidra, and a portion of Barranquitas. The figure also shows LE 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⁻¹.

Figure 7.

Example of water and energy balance components from the GOESWEB algorithm for November 24, 2015.

Figure 8.

NOAA’s Drought Monitor for Puerto Rico, October 18, 2018. The municipalities of Aibonito, Cayey, and Cidra are experiencing abnormally dry conditions.

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 (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.

Figure 9.

Rainfall over Puerto Rico during the week of Hurricane Maria, during September 2017.

Figure 10.

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.

Figure 11.

Estimated surface runoff over Puerto Rico on September 20, the day that Hurricane Maria made landfall on Puerto Rico.

Figure 12.

Root zone soil moisture saturation for September 19, 2018, 1 day before Hurricane Maria made landfall on Puerto Rico.

5. 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.

Figure 13.

ET_o for October 16, 2018, for (a) Puerto Rico, (b) Hispaniola, (c) Jamaica, (d) Cuba, and (e) St. Croix, (f) St. Thomas, and (g) St. John (USVI).

6. 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.

Acknowledgments

This work was supported by the NOAA-CREST (grant NA06OAR4810162), National Science Foundation (grant 0313747 and 1832576), and the U.S. Department of Agriculture (Hatch Project H402). Any opinions, findings, conclusions, or recommendations expressed in this material are those of the authors and do not necessarily reflect those of the National Oceanic and Atmospheric Administration, NSF, or the USDA. Reference to a commercial product in no way constitutes an endorsement of the product by the authors.

Conflict of interest

None.

Notes/thanks/other declarations

Thanks to Luis Aponte-Bermudez for providing the Matlab® computer code and for reading the daily WRF wind speed data sets and CariCOOS for providing the WRF wind speed, air temperature, and relative humidity data. There are numerous other people that directly or indirectly assisted in this research and to those people we would like to express our appreciation.

Nomenclature

B (unitless)
Cf (dimensionless)
CN (dimensionless)
DP (mm)
dr (dimensionless)
ETa (mm)
ETo (mm)
ea [kPa]
es (kPa)
g (m s−2)
G (MJ m−2 day−1)
Gsc (MJ m−2 min−1)
h (m)
H (MJ m−2 day−1)
J
k (dimensionless)
Kc,eff (unitless)
Ks (unitless)
LE (MJ m−2 day−1)
PAR (micromoles m−2 sec−1)	photosynthetically active radiation
Precip (mm)
Precipeff (mm)
Ra (MJ m−2 day−1)
ra (s m−1)
rao (s m−1)	aerodynamic resistance under conditions of neutral atmospheric stability
rbh (s m−1)
Rdepth (mm)
RH (%)
Rn (MJ m−2 day−1)
Rnl (MJ m−2 day−1)
Rns (MJ m−2 day−1)
RO (mm)
Rs (MJ m−2 day−1)
rs (s m−1)
Rs/Rso
Rso (MJ m−2 day−1)
S (mm)	maximum potential difference between rainfall and runoff at the moment of rainfall initiation
Sat (unitless)
SMD1 (mm)
SMD2 (mm)
Ta (°C)
Tavg [°C]
Td (°C)
Tmax (K)
Tmax (°C)
Tmin (°C)
Tmin (K)
To (°C)
Ts (°C)
U (m s−1)
u2 (m s−1)
VPD (kPa)
z (m)
zdisp (m)
zo (m)
α
δ (rad)
Δ [kPa °C−1]
η (dimensionless)
φ (rad)
σ (MJ K−4 m−2 day−1)
γ (kPa °C−1)
ω1 (rad)
ω2 (rad)
θ (m−3 m−3)
θFC (m3/m3)
θWP (m3/m3)