Water and Energy Balance Model GOES-PRWEB: Development and Validation

: In 2009, the University of Alabama-Huntsville conﬁgured their GOES satellited-based solar radiation product to include Puerto Rico, the US Virgin Islands (USVI), Dominican Republic, Haiti, Jamaica, and Cuba. The half-hourly and daily integrated data are available at 1 km resolution for Puerto Rico and the USVI and 2 km for Hispaniola, Jamaica, and Cuba. These data made it possible to implement estimates of satellite radiation-based evapotranspiration methods on all of the islands. The use of the solar radiation data in combination with estimates of other climate parameters facilitated the development of a water and energy balance algorithm for Puerto Rico. The purpose of this paper is to describe the theoretical background and technical approach for estimating the components of the daily water and energy balance. The operational water and energy bal-ance model is the ﬁrst of its kind in Puerto Rico. Model validation results are presented for reference and actual evapotranspiration, soil moisture, and streamﬂow. Mean errors for all analyses were less than 7%. The water and energy balance model results can beneﬁt such diverse ﬁelds as agriculture, ecology, coastal water management, human health, renewable energy development, water resources, monitoring, and disaster and emergency management. This research represents a preliminary step in developing a suite of gridded hydro-climate products for the Caribbean Region.


Introduction
Hydrologic water budgets are essential because they provide the information needed to evaluate a region or country's water resources. The regional-scale evaluation helps forecast floods and droughts, maintain sustainable water supplies, and assess the impacts of changing climate and land-use changes on water resources ( [1,2]). Tropical island nations are especially vulnerable to environmental disasters ( [3]). Floods and droughts are common occurrences in the Caribbean region ( [4]). To better understand the hydrologic processes in these islands, tools are needed to estimate the temporal and spatial distribution of the components of the hydrologic cycle, along with other relevant agro-hydro-meteorological variables. Accurate methods are needed to estimate rainfall, evapotranspiration, runoff, where R n is net radiation, the sensible heat flux (H) is where C p is the specific heat [MJ kg −1 • C −1 ] and ρ is the mean air density at constant pressure (kg/m −3 ), T s is surface temperature [ • C], T a is air temperature [ • C], and r a is aerodynamic resistance (sm −1 ). Latent heat flux (LE) is estimated similar to [40], LE = ρ ·C p ·(e(T s ) − e(T a ) γ·(r a + r s ) where r s is surface resistance (sm −1 ), γ is psychrometric constant [kPa • C −1 ]. Vapor pressures were estimated with the equation: e = 0.6108 (exp [(17.27T)/(T + 237. 3)]) [41], where T is T a for the above-canopy vapor pressure, and T is T s for the within-canopy vapor pressure. The soil heat flux (G) is assumed to be zero for the daily analysis. The units of all terms in equation 1 are [MJ m −2 day −1 ]. Equation (1) has only one unknown variable, T s , which is determined by an implicit approach described by [42], which makes no assumption as to the temperature of the evaporating surface. T s is obtained using the recursive root function fzero in MatLab ® (http://www.mathworks.com, accessed on 15 May 2021). The advantage of solving directly for T s is that it eliminates the need to measure surface temperature, which can be hindered by cloud cover ( [43]), and to correct for elevation, slope, and aspect of the surface ( [44]). R n is obtained from the calculation procedure, as presented by [41], which requires R s . In this study, R s is estimated using NOAA's GOES satellite visible imagery. R s is derived from a radiative transfer model first proposed by Gautier et al. [45] (see also [20,24,46,47]), using daily integrated solar radiation data, currently available from the visible channel 2 (0.64 µm) of NOAA's GOES-16 satellite. Before GOES-16, 1 km GOES-12, and -13 visible channel 1 data were used over Puerto Rico. Albedo is obtained from a look-up table [48], assigning the parameters' values to 32 different land covers. The ground-level, 1-km resolution R s product was validated at two locations in Puerto Rico by [25].
Aerodynamic resistance is calculated with the following equation [40]: r a = r ao ·φ + r bh (4) where r ao [s m −1 ] is the aerodynamic resistance under conditions of neutral atmospheric stability [s m −1 ]: In Equation (5), z is the virtual height at which meteorological measurements are taken [m], assumed to be 1.5(z o /0.13) [49]. The 10-m NDFD-or WRF-derived wind speeds are adjusted to the "virtual instrument height," depending on the vegetation's height. The roughness length (z o ) [m] and the zero plane displacement (z disp ) [m] for various land use/vegetation categories were obtained from [48]. z o /0.13 is equal to the canopy height (h) [m]. k is the Von Karman's constant (k = 0.41) [dimensionless]. Wind speed at height z is the wind velocity (u) [m]. The atmospheric stability coefficient [dimensionless] is estimated from [40]: where g is the gravitational constant [m s −2 ], and the coefficient η is commonly taken as 5 [40] [dimensionless]. T o is the average of T s and T a [ • C]. The excess resistance [s m −1 ] in Equation (4) is given by the following equation: The surface or canopy resistance (r s ) is given by [50]: VPD is the vapor pressure deficit [KPa], C f is a coefficient equal to 1 for root depth <1 m and 5 for root depth >1 m in this study [dimensionless], θ is root zone volumetric soil moisture content [m 3 /m 3 ], and θ FC and θ WP are the θ values at field capacity and wilting point [m 3 /m 3 ], respectively. Pedotransfer functions, based on percent sand, silt, and clay, were used to estimate field capacity and wilting point ( [51]). Sand, silt, and clay were obtained from the Soil Survey Geographic Database (SSURGO) of the USDA Natural Resource Conservation Service (NRCS) (https://datagateway.nrcs.usda.gov/, accessed on 15 May 2021).

Water Balance
Actual evapotranspiration is obtained by converting LE for each pixel using the latent heat of vaporization (2.45 MJ kg -1 ). The water balance equation used in this study is presented below ( [20]): SMD2 = P − ET a − RO − DP + SMD1 (9) where SMD1 and SMD2 are the water depth within the soil profile [mm] at the beginning and end of each day, respectively, P is precipitation [mm], RO is surface runoff [mm], and DP is deep percolation [mm]. The water balance is conducted over a depth equal to the root depth (R depth ). Root depth [m] for various land use/vegetation categories was obtained from [48]. Twenty-four-hour rainfall is obtained from NOAA's Advanced Hydrologic Prediction Service (AHPS) website (https://water.weather.gov/precip/, accessed on 15 May 2021). Currently the formulation does not account for soil freezing/thawing and snow storage/melt, although these processes could be added.
RO is estimated using the curve number (CN) method of the soil conservation service [52]: RO = (P − 0.2S) 2 /(P + 0.8S) (10) where S represents the maximum possible difference between P and RO at the moment of rainfall initiation [mm] and CN is a proportion of rainfall converted to runoff [dimensionless] and is adjusted for antecedent rainfall conditions. In this study, CN values were derived for Puerto Rico using the method described by [52], based on hydrologic soil group and land use. Adjustments for antecedent rainfall condition (ARC) were made to CN values based on the following criteria: for 0 ≤ soil saturation (Sat) ≤ 0. 33 [52]; and for ARC II, CN table values are not adjusted (average conditions). In this study we define Sat = 1 for θ = θ FC and Sat = 0 for θ = θ WP . We recognize that the upper limit of Sat is not precisely correct, as it should be equal to the total porosity; however, because the model distributes water to deep percolation at the end of each day, θ FC is effectively the largest value possible in the model. SMD2 is initially estimated as SMD2 i = P − ET a − RO + SMD1. If SMD2 i exceeds the field capacity (FCD), then DP = SMD2 i − FCD, and the value of SMD2 is equal to FCD. If, however, SMD2 i < FCD, then DP = 0, and SMD2 = SMD2 i .

Reference Evapotranspiration
GOES-PRWEB calculates reference evapotranspiration using the Penman-Monteith (PM) equation [41,53]: where ∆ is the slope of the saturated vapor pressure curve [kPa • C −1 ], T is mean daily T a at 2 m height [ • C], u 2 is the wind speed at 2 m height [m s −1 ], e s is the saturated vapor pressure, and e a is the actual vapor pressure [kPa]. Equation (1) applies specifically to a hypothetical reference crop with an assumed crop height of 0.12 m, a fixed r s of 70 s m −1 , and an albedo of 0.23.

GOES-PRWEB Products and Data
Currently, water and energy balance calculations are performed daily for the island of Puerto Rico. Twenty-seven hydro-agro-climate variables are available to the public for download as images (jpg), or in comma-separated values (CSV) and Matlab ® formats, and include: ET a , reference ET (ET o , three methods), minimum, maximum, and average air temperature (T a ), dew point temperature (T d ), effective surface temperature (T s ), saturated and actual vapor pressures (e), relative humidity (RH), wind speed (u), solar radiation (R s ), net radiation (R n ), photosynthetically active radiation (PAR), water stress coefficient (K s ), effective crop coefficient (K c,eff ), rainfall (P), effective rainfall (P eff ), surface runoff (RO), deep percolation (DP), soil moisture (θ), soil saturation (Sat), surface resistance (r s ), aerodynamic resistance (r a ), latent heat flux (LE), sensible heat flux (H), and Bowen Ratio (β). Monthly and annual averages or totals for all variables are available.
Harmsen et al. [54] have described the methodology for estimating ET o . Daily minimum, maximum, and average T a values were calculated from a lapse rate method developed by Goyal et al. [55]. The method works well in PR, where the maximum topographic elevation (MTE) is 1340 m. Other islands with significant topographic relief include Hispaniola (MTE 3098 m) and Jamaica (MTE 1072 m). T d data were assumed to be equal to the minimum daily T a [56]. These simplified methods continue to be used for estimating daily values of ET o for the USVI, Hispaniola, Jamaica, and Cuba, including the use of the worldwide average 2-m wind speed of 2 m s -1 in the Penman-Monteith method [41]. Various modifications have been implemented in the algorithm since 2009 (applicable to Puerto Rico), which are described below.
Two sets of wind speed data were used during the life of the operational model. Data set 1: from1 January 2009 to 30 September 2015. An average of eight 3-h wind speed values obtained from the National Weather Service's National Digital Forecast Database ( [57]). Data set 2: from 1 October 2015 to the present. An average of 24-hourly values of wind speed obtained from the Caribbean Coastal Ocean Observing System (CARICOOS) Weather Research Forecast (WRF) 1-km resolution model. These data are used to estimate the average daily wind speed. The 10-m wind speeds are adjusted to 2 m ( [41]) for the ET o calculation.
Three-hourly minimum, average, maximum T a , and T d data were obtained from the National Digital Forecast Database (NDFD) website [57] from 1 January 2009 to 31 December 2016. Starting on 1 January 2017, hourly weather parameters were obtained from the CARICOOS operational gridded WRF model. Average air temperature, dew point temperature, and wind speed from the two sources were converted to daily averages for input into GOES-PRWEB. Occasionally, estimates of the weather parameters from the WRF model were not available. In those cases, the NDFD T a and wind speed were used. When these data were not available, a lapse rate method was used for T a , and wind speed from the previous day was used.

Study Area
The highly diversified tropical conditions of Puerto Rico represent an interesting test case. The island has climate zones that vary from tropical rainforest (4500 mm rainfall per year) to semi-arid dry forest (750 mm rainfall per year). The region has a wet season from May to November and a dry season from December to April. Capiel and Calvesbert [58] attributed the dry season to a temperature inversion that acts as a general circulation valve, opposing vertical cloud development. There is also a dry period in late June/early July, known as the mid-summertime drought, affecting much of the northern Caribbean region. This drier period during the mid-summer has been attributed to the increase in aerosol concentrations, primarily due to Saharan dust and upper-level vertical wind shear, which tends to inhibit tropical storm formation [59]. Hurricanes and extreme weather events are common on the island (e.g., 2014-2016 drought and Hurricane María in 2017). Wet air carried by the easterly trade wind moves from the northeast to the southwest of the island. Much of this moist air produces orographic rainfall due to mountains that run frp, east to west through the middle of the island (Cordillera Central), with maximum elevation of 1340 m above mean sea level (amsl) at Cerro de Punta; the rain shadow effect results in an arid southern and southwest coastal region. In the west, the sea breeze effect carries wet air from the Mona Channel eastward, converging with the trade wind and resulting in intense convective rainstorms almost every afternoon during the wet season. Jury et al. [60] have described the diurnal cycle and the sea/land breeze's mesoscale features in western Puerto Rico during undisturbed weather conditions (see also [61]).
Puerto Rico has various data sources available for supporting hydrologic analyses. The National Weather Service (NWS) operates the next generation radar (NEXRAD), known as Weather Surveillance Radar 1988 Doppler (WSR−88D), located at Cayey, PR, and maintains 63 rain gauges. The USGS maintains 84 rain gauges, 124 stream gauges, 20 groundwater observation wells and 26 lake/reservoir level gauges. 2021. For the validation study, ET o was estimated using the Penman-Monteith method. Comparisons of estimated maximum and minimum T a (T max and T min , respectively), T d and u 2 (obtained from gridded forecast models), and R s (obtained from satellite) are provided. The remotely sensed R s data were compared with pyranometer data in a previous study by [25] in Puerto Rico, who reported that the uncorrected data produced reasonably accurate results with a maximum 6.22% error between the mean daily estimated and measured R s . In this study, observed and simulated means were evaluated for significant difference at the 0.05 significance level using the two-tailed Student's ttest (www.socscistatistics.com/tests/studentttest/, accessed on 15 May 2021). Test results will be presented as: t (degress of freedom) = the t statistic, p = p value.

Evaluation of Soil Moisture Estimates
GOES-PRWEB produces a single value of θ for the 1-km 2 pixel area. We hypothesized that there was high variability in θ within a 1-km 2 GOES pixel. Notwithstanding the measured variability, we further hypothesized that GOES-PRWEB could estimate the average measured θ for the 1-km 2 area. Undisturbed soil samples were taken to measure the water content at various locations within the studied pixel to test these hypotheses. The study area is located on the University of Puerto Rico-Mayaguez (UPRM) campus in western Puerto Rico. Samples were obtained on 1 October, 13 October, 7 November, and 5 December 2015, and 7 June, 27 June, and 18 August 2016.
Sample locations are presented in Figure 1; some areas of the study pixel in the southwest and the northeast fall outside the UPRM campus. The study area was selected because it corresponded with a GOES visible channel pixel. Soil cores were used to collect soil samples within the study pixel for water content analyses. GPS coordinates were obtained, and organic matter was cleared from the area upon arrival. The metal cores were hammered into the soil with a mallet until the top rim of the core was level with the surface, and a shovel was used to carefully extract the soil-filled core from the ground. Soil samples were quickly removed from the soil core, inserted into a labeled plastic bag, and weighed in the field to record the exact weight of the wet or damp soil. Later, the soil samples were dried with a microwave oven for approximately 20 min, and the dry weight was recorded. The undisturbed θ was estimated from the following equation: θv [cm 3 cm −3 ] = (WWnudist − DWundist)/(Vol ρH2O), where WWundist is undisturbed soil wet weight (gm), DWundist is undisturbed soil dry weight (gm), Vol is the volume of the sample core (cm 3 ), and ρH2O is the density of water (1 gm cm -3 ). Average θ estimates from weather station sensors (Watchdog ET900, Spectrum Technology, Inc., Aurora, IL, USA) located at the Agricul- The undisturbed θ was estimated from the following equation: θ v [cm 3 cm −3 ] = (WW nudist − DW undist )/(Vol ρ H2O ), where WW undist is undisturbed soil wet weight (gm), DW undist is undisturbed soil dry weight (gm), Vol is the volume of the sample core (cm 3 ), and ρ H2O is the density of water (1 gm cm -3 ). Average θ estimates from weather station sensors (Watchdog ET900, Spectrum Technology, Inc., Aurora, IL, USA) located at the Agricultural Engineering Building on the UPRM Campus were also compared with the gravimetric and GOES-PRWEB estimates. The soil moisture sensors were installed at 0.3 m and 0.6 m depths.
Another soil moisture comparison study was conducted over three years at the UPR Agricultural Experimental Station near Juana Díaz, Puerto Rico. The 1-km 2 GOES-PRWEB pixel was assigned a root depth of 1 m, while the station data were based on a weighted average of θ sensors (Time Domain Reflectometry) at depths of 0.0508 m, 0.1016 m, 0.2032 m, 0.508 m, and 1.016 m.

Island-Scale Comparisons
Two island-scale validation studies were conducted:

•
Annual ET a from GOES-PRWEB and the SSEBop model were compared for 2009-2020; • Annual values of the water balance components from GOES-PRWEB and the USGS ( [62]) were compared for 2009-2020.

Basin-Scale Comparisons
Two basin-scale validation studies were conducted: • Comparisons of cumulative monthly streamflow were made for GOES-PRWEB and the US Geological Survey (USGS) stream gage located at the outlet of the Guanajibo watershed in southwest Puerto Rico for the years 2010, 2011, and 2012. Total streamflow was assumed to be the combined flow from surface runoff and deep percolation, where the latter contributes to the stream base flow. Since the model does not explicitly account for groundwater storage, this may be a source of error during short periods. However, over more extended periods, such as a year and during hydrologically normal years, the change in storage was small. • GOES-PRWEB and USGS-based annual ETa estimates were compared for the Guanajibo watershed for 2009-2020. The USGS-based ETa was estimated as P − (RO + DP), where (RO + DP) is considered to be equal to the total measured streamflow.

Pixel Scale Comparison
A pixel-scale validation study was conducted: A comparison of monthly ET a at three locations in Puerto Rico from January 2010 to December 2012 was conducted.

Water Balance Error Analyses
To verify the numerical correctness of the model, island-scale and pixel-scale water balance error analyses are presented.

Results
In this section we present results to support the validation of the model. To better understand the modelestimated ET o , it is necessary to consider the model's climate input data. For the interested reader, graphs of the associated T min , T max , T d , u 2 and R s are presented in Appendix A. The temperature and wind speed data used as input to GOES-PRWEB were derived from gridded model data from NOAA's NDFD website, and R s was derived from the GDM radiative transfer method using GOES-13 satellite data. The T min ( Figure A1) was slightly overestimated relative to the weather station data. On average, T max ( Figure A2) was consistent with the station data; however, gridded data was characterized with more significant variability during specific periods. T d ( Figure A3) was under-estimated relative to the weather station data. u 2 data ( Figure A4) was under-estimated relative to the weather station data, and R s data ( Figure A5) was reasonably accurate relative to the weather station. Figure 2 shows a time series of the estimated ETo at the UPR Agricultural Experimental Station near Juan Diaz, Puerto Rico (latitude 18.033 degrees, longitude −66.533 degrees) from January 2014 to April 2016. The GOES-PRWEB-derived ETo is in good agreement with the ETo from the NRCS SCAN weather station. To better understand the modelestimated ETo, it is necessary to consider the model's climate input data. For the interested reader, graphs of the associated Tmin, Tmax, Td, u2 and Rs are presented in Appendix A. The temperature and wind speed data used as input to GOES-PRWEB were derived from gridded model data from NOAA's NDFD website, and Rs was derived from the GDM radiative transfer method using GOES-13 satellite data. The Tmin ( Figure A1) was slightly overestimated relative to the weather station data. On average, Tmax ( Figure A2) was consistent with the station data; however, gridded data was characterized with more significant variability during specific periods. Td ( Figure A3) was under-estimated relative to the weather station data. u2 data ( Figure A4) was under-estimated relative to the weather station data, and Rs data ( Figure A5) was reasonably accurate relative to the weather station.   The period corresponding with Figure 4 was rainy, and the soil was relatively wet. This is observable because the GOES-PRWEB θ reached the soil field capacity of 0.45 on numerous occasions. During one sampling event on 18 June 2016 ( Figure 5), the measured pixel average θ was significantly lower than GOES-PRWEB (difference approximately 0.08). During the summer of 2016, θ from the weather station sensors frequently dropped below GOES-PRWEB. However, this discrepancy is not of great concern, considering that the weather station water content represents only a single point within the 1-km 2 pixel area. It should be noted that the undisturbed core samples were taken from the top 0.127 m, the two weather stations sensors were installed at 0.3 and 0.6 m, and the GOES-PRWEB assigned root depth was from 0 to 0.87 m. These differences in the three methods' soil control volumes may have contributed to the observed discrepancies.         As another check on the ability of the model to simulate θ, Figure 6 shows a comparison of the estimated and measured (depth-averaged) daily volumetric moisture content at the UPR Agricultural Experiment station near Juana Diaz, Puerto Rico, from 1 January 2014 to 31 December 2016. Soil moisture sensors were located at the following depths: 0.0508 m, 0.1016 m, 0.2032 m, 0.508 m, and 1.016 m. The average GOES-PRWEB θ, average station θ, average percent error, and R 2 for the data presented in Figure 6 were 23.9%, 24.39%, −2.13%, and 0.75, respectively. A Student t-test was performed on 40 randomly selected soil moisture pairs from the 3-year dataset, indicating no significant difference at As another check on the ability of the model to simulate θ, Figure 6 shows a comparison of the estimated and measured (depth-averaged) daily volumetric moisture content at the UPR Agricultural Experiment station near Juana Diaz, Puerto Rico, from 1 January 2014 to 31 December 2016. Soil moisture sensors were located at the following depths: 0.0508 m, 0.1016 m, 0.2032 m, 0.508 m and 1.016 m. The average GOES-PRWEB θ, average station θ, average percent error, and R 2 for the data presented in Figure 6 were 23.9%, 24.39%, −2.13%, and 0.75, respectively. A Student t-test was performed on 40 randomly selected soil moisture pairs from the 3-year dataset, indicating no significant difference at the 0.05 level between the mean observed and simulated values (t(39) = 0.42104, p = 0.674854).

Basin-Scale Streamflow Comparison
The cumulative simulated and observed stream flows from the Guanajibo watershed (USGS gage station at Hormigueros, No. 50138000, catchment area 310.9 km 2 ) were compared. Figure 7 shows the location of the gaged catchment area (red) and the entire watershed (black).

Basin-Scale Streamflow Comparison
The cumulative simulated and observed stream flows from the Guanajibo watershed (USGS gage station at Hormigueros, No. 50138000, catchment area 310.9 km 2 ) were compared. Figure 7 shows the location of the gaged catchment area (red) and the entire watershed (black).

Basin-Scale Streamflow Comparison
The cumulative simulated and observed stream flows from the Guanajibo w (USGS gage station at Hormigueros, No. 50138000, catchment area 310.9 km 2 ) w pared. Figure 7 shows the location of the gaged catchment area (red) and the e tershed (black).  The mean error in the monthly flow for the Guanajibo River for the 36-mo parison is −2.6%, indicating that simulated streamflow was, in general, less tha served. In general, the simulated vs. observed streamflow was in good agreemen ing a two-tailed Student t-test, no significant difference was found between the and simulated mean monthly streamflow at the 0.05 level (t(35) = 0.17268, p = Figure 9 shows the simulated versus the observed stream discharge and Figure the monthly simulated versus observed stream discharge. The results show good agreement; however, there were four months (months 5, 17, 33 and 34) model significantly overestimated the stream discharge. It is also noted that t tended to underestimate stream flow during periods of low flow. The calculat cient of determination for the monthly analysis was r 2 = 0.7136.  The mean error in the monthly flow for the Guanajibo River for the 36-month comparison is −2.6%, indicating that simulated streamflow was, in general, less than the observed. In general, the simulated vs. observed streamflow was in good agreement. Applying a two-tailed Student t-test, no significant difference was found between the observed and simulated mean monthly streamflow at the 0.05 level (t(35) = 0.17268, p = 0.863403). Figure 9 shows the simulated versus the observed stream discharge and Figure 10 shows the monthly simulated versus observed stream discharge. The results show generally good agreement; however, there were four months (months 5, 17, 33 and 34) when the model significantly overestimated the stream discharge. It is also noted that the model tended to underestimate stream flow during periods of low flow. The calculated coefficient of determination for the monthly analysis was r 2 = 0.7136.

Basin-Scale Actual Evapotranspiration Comparison
Independent flux tower data are not available for the period of the study in Puerto Rico. As an alternative, we compared the basin-scale annual average GOES-PRWEB ET a with a USGS-based ET a . In the latter case, ET a = P − (RO + DP), where RO + DP is the total streamflow derived from the USGS stream gage on the Guanajibo River near Hormigueros, Puerto Rico (Figure 7), and P data were obtained from NOAA's AHPS website, which is bias-corrected radar rainfall. Table 2 shows the annual basin analyses for 2009-2020 (2016 was removed; see explanation below). The minimum, maximum, and mean errors in the model ET a for the 12 years were −1.35% (2009), −18.96% (2014), and −6.8%, respectively. Error is equal to 100 (GOES-PRWEB ET a − USGS ET a )/USGS ET a . When applying a two-tailed Student t-test, no significant difference at the 0.05 level was found between the observed and simulated mean annual ET a (t(11) = 1.40728, p = 0.174695). The 2016 ET a estimate from AHPS rainfall and USGS streamflow data was abnormally high and therefore was discarded from the analysis. The cause of the high estimated ET a was the anomalously low value of streamflow, despite that year having the highest recorded rainfall for all the years considered. Figure 11 shows that the 2016 streamflow was an outlier, hence justifying its removal from the dataset. streamflow derived from the USGS stream gage on the Guanajibo River near Hormigueros, Puerto Rico (Figure 7), and P data were obtained from NOAA's AHPS website, which is bias-corrected radar rainfall. Table 2 shows the annual basin analyses for 2009-2020 (2016 was removed; see explanation below). The minimum, maximum, and mean errors in the model ETa for the 12 years were −1.35% (2009), −18.96% (2014), and −6.8%, respectively. Error is equal to 100 (GOES-PRWEB ETa − USGS ETa)/USGS ETa. When applying a two-tailed Student t-test, no significant difference at the 0.05 level was found between the observed and simulated mean annual ETa (t(11) = 1.40728, p = 0.174695).
The 2016 ETa estimate from AHPS rainfall and USGS streamflow data was abnormally high and therefore was discarded from the analysis. The cause of the high estimated ETa was the anomalously low value of streamflow, despite that year having the highest recorded rainfall for all the years considered. Figure 11 shows that the 2016 streamflow was an outlier, hence justifying its removal from the dataset.

Island-Wide ETa Evaluation
Annual actual evapotranspiration from the Simplified Surface Energy Balance (SSE-Bop) model [38] was compared with GOES-PRWEB for 2009-2020 (Table 3)  In most years, the percent difference (error) between the two models was small. The largest errors occurred in 2015 and 2018. Puerto Rico experienced a devastating drought during 2015 (AHPS rainfall 1746 mm). On the other hand, 2018 was a relatively wet year (AHPS rainfall 2205 mm). The overall average error was −0.73 percent, and there was no significant difference at the 0.05 level between the mean SSEBop and GOES-PRWEB annual ET a (Student-t statistics: t(11) = 0.30016, p = 0.767154). Figure 12 shows the GOES-PRWEB and SSEBop annual ETa spatial distribution for 2015. SSEBop produced lower values of ETa in the San Juan Metropolitan area. It has been observed that the SSEBop model produces low values of ET a in coastal urban areas (personal communication, Nick Sepulveda, USGS). On the other hand, GOES-PRWEB may be overestimating ETa in the Metro area. Note that the SSEBop model produced very high values all around the coast (2666 mm) in the pixels adjacent to the ocean. Presumably, SSEBop produces high values of ETa in these pixels because the MODIS thermal image combines ocean and land within the same pixels. Consequently, the estimate is approaching the potential evapotranspiration. GOES-PRWEB does not suffer from this problem because the land surface temperature estimation is based only on the energy balance of the land properties. For the comparison in this study (Table 3), the high SSEBop ETa values along the coasts were removed.

Island-Wide ETa Evaluation
Annual actual evapotranspiration from the Simplified Surface Energy Balance (SSEBop) model [38] was compared with GOES-PRWEB for 2009-2020 (Table 3). SSEBop, a product of the Famine Early Warning Systems Network (FEWS NET), estimates ETa based on a fractional weighting of ETo using Moderate Resolution Imaging Spectrometer (MODIS) thermal imagery, obtained every eight days. The original formulation of the SSEBop model is based on the SEBAL ( [63]) and METRIC ( [35]) algorithms. In most years, the percent difference (error) between the two models was small. The largest errors occurred in 2015 and 2018. Puerto Rico experienced a devastating drought during 2015 (AHPS rainfall 1746 mm). On the other hand, 2018 was a relatively wet year (AHPS rainfall 2205 mm). The overall average error was −0.73 percent, and there was no significant difference at the 0.05 level between the mean SSEBop and GOES-PRWEB annual ETa (Student-t statistics: t(11) = 0.30016, p = 0.767154). Figure 12 shows the GOES-PRWEB and SSEBop annual ETa spatial distribution for 2015. SSEBop produced lower values of ETa in the San Juan Metropolitan area. It has been observed that the SSEBop model produces low values of ETa in coastal urban areas (personal communication, Nick Sepulveda, USGS). On the other hand, GOES-PRWEB may be overestimating ETa in the Metro area. Note that the SSEBop model produced very high values all around the coast (2666 mm) in the pixels adjacent to the ocean. Presumably, SSEBop produces high values of ETa in these pixels because the MODIS thermal image combines ocean and land within the same pixels. Consequently, the estimate is approaching the potential evapotranspiration. GOES-PRWEB does not suffer from this problem because the land surface temperature estimation is based only on the energy balance of the land properties. For the comparison in this study (Table 3), the high SSEBop ETa values along the coasts were removed.

ET a Evaluation at Three Locations
Here, we compare the GOES-PRWEB and SSEBop ET a at three locations in Puerto Rico: Mayagüez, Guánica, and Orocovis. The conditions at each pixel are summarized in Table 4. Figure 13A-C show the GOES-PRWEB versus SSEBop ET a for Mayagüez, Guánica and Orocovis, respectively. It should be noted that the 1-km resolution SSEBop ET a product is not calibrated for PR. Consequently, we will use the term "Deviations" instead of errors. Deviations were estimated as 100(GOES-PRWEB ET a − SSEBop ET a )/SSEBop ET a . The mean monthly deviations were 0.07%, 25.14%, and 9.8% for Mayagüez, Guánica, and Orocovis, respectively.

ETa Evaluation at Three Locations
Here, we compare the GOES-PRWEB and SSEBop ETa at three locations in Puerto Rico: Mayagüez, Guánica, and Orocovis. The conditions at each pixel are summarized in Table 4. Figure 13A-C show the GOES-PRWEB versus SSEBop ETa for Mayagüez, Guánica and Orocovis, respectively. It should be noted that the 1-km resolution SSEBop ETa product is not calibrated for PR. Consequently, we will use the term "Deviations" instead of errors. Deviations were estimated as 100(GOES-PRWEB ETa − SSEBop ETa)/SSEBop ETa. The mean monthly deviations were 0.07%, 25.14%, and 9.8% for Mayagüez, Guánica, and Orocovis, respectively. Paired comparisons of the mean GOES-PRWEB and mean SSEBop for the three locations, using the Student's t-test, were performed. While there was no significant difference at the 0.05 significance level between the Mayaguez means (t(35) = −0.00731, p = 0.99419) and Orocovis means (t(35) = −1.6481, p = 0.103813), there was a significant difference between the means for Guánica (t(35) = −3.29467, p = 0.000774). Because the SSEBop ET a product is not calibrated for PR, it is not possible to determine which model is more accurate. Nevertheless, there is a degree of correlation between the two models ( Figure 13A-C), which provides credibility to the GOES-PRWEB results for the relatively diverse environmental conditions evaluated (i.e., elevation, land cover and climate).

Island-Wide Water Balance Component Comparison
The USGS published estimates of the average water balance components for Puerto Rico ( [62]) as follows: P 1829 mm, ET a 1168 mm, streamflow (RO + base flow [=DP]) 635 mm, groundwater discharge from coastal aquifers to wetlands, estuaries and seabed 25 mm, groundwater withdrawals from coastal aquifers 25 mm, and soil moisture and groundwater storage 25 mm. In terms of rainfall percentages: ET 63.9%, total streamflow (RO + DP) 34.7%, groundwater discharge from coastal aquifers to wetlands, estuaries, and seabed 1.3%, groundwater withdrawals from coastal aquifers 1.3%, and soil moisture and groundwater storage 1.3%.
GOES-PRWEB-estimated water balance components for Puerto Rico from 2009-2020 are presented in Table 5 and compared with the USGS estimates (included in the table). The table shows the water balance components in depths of water in mm and as percentages of rainfall. Average GOES-PRWEB P was 1888.8 mm, ET was 1176.5 mm, and total streamflow (RO + DP) was 728.3 mm.
The GOES-PRWEB and USGS island-wide average ET a , as a percentage of rainfall, were 62.3% and 63.9%, respectively; and the RO + DP were 38.6% and 34.7%, respectively. The island-wide average ET a and RO + DP results from GOES-PRWEB and the USGS are in good agreement. * DP for the USGS includes aquifer recharge, groundwater discharge from coastal aquifers to wetlands, estuaries, seabed, and groundwater withdrawals from coastal aquifers.

Pixel-Scale Water Balance Error Analyses
In this section, the water balance errors are analyzed at the pixel-scale at Mayagüez, Guánica and Orocovis (Table 4). Ideally, P i − ET ai − RO i − DP i + (SM i−1 − SM i ) should equal zero; however, due to numerical errors in the algorithm, water balance errors can occur. Numerical errors may include rounding errors or lack of convergence by the recursive root function employed in the methodology. The subscripts i and i − 1, represent the current and previous time steps. Table 6 presents the average overall mean, minimum, and maximum balance errors for the three locations. The sample size was n = 12 (2009-2020). All values are in millimeters. The overall mean annual errors for the three locations are all less than 0.05 mm and can be considered negligible. The largest overall mean annual negative error was −1.38 mm, while the largest overall mean annual positive error was 21.6 mm. The results in Table 6 are based on mean annual estimated errors presented in Appendix B (Tables A1-A3). Table 7 compares the model water balance error as a percent of rainfall from 2009-2020. Annual change in soil moisture is the difference between the soil moisture on the first and last day of the year. If the model has no water balance error, then the value in column 6 would be zero. The average water balance error is −0.7%. The annual water balance errors are sufficiently small, indicating that the island-wide water balance is acceptable.

Model Applications
The model has been used for various practical applications. For example, [64] reported a web-based GOES-PRWEB method for scheduling irrigation in PR, USVI, Hispaniola, Jamaica, and Cuba. The USGS used GOES-PRWEB to estimate water use by agricultural crops as part of a water withdrawal and use study in Puerto Rico ( [65]). The Scientific Drought Committee of Puerto Rico uses data from the model in their weekly reports [66]. NOAA's National Weather Service (NWS) in San Juan links the near-real-time soil moisture and soil saturation products on their Climate and Drought Information page (https:// www.weather.gov/sju/dss_climo, accessed on 15 May 2021). The U.S. Forest Service and the National Integrated Drought Information System (NISDIS) uses soil moisture, soil saturation, and rainfall deficit information from the model for their bi-monthly reports in Puerto Rico and the USVI. Soil saturation from the model was recently used in a Hurricane María flood modeling study in western Puerto Rico ( [67]). Figure 3 shows considerable variability in θ within the 1-km × 1-km study pixel. This finding illustrates the importance of obtaining a multiple soil samples within a study pixel when attempting to validate the soil water content derived using a satellite method. A validation study in which only a single value of soil moisture is used may be insufficient. The need for a greater number of samples increases with coarser resolution soil moisture remote sensing products, such as the Advanced Microwave Scanning Radiometer 2 (AMSR2) with a 25 km spatial resolution ( [68]) and the soil moisture active passive (SMAP) mission with a 9 km spatial resolution ( [69]).

Discussion of Selected Validation Results
The modeled and measured θ for the three years, presented in Figure 6, are in reasonably good agreement. However, the model did not produce θ values as high as the station immediately after most rainfall events. A possible explanation for the higher observed θ values is that the mean observed θ is weighted towards the shallow sensors because more sensors are near the surface. Another possible reason is that the observed data reflect moisture contents near the total porosity, whereas the maximum possible moisture content in the model is the field capacity. The above explanation may apply to the slight underestimation of θ by the model during wet periods shown in Figures 4 and 5. It is essential to keep in mind precisely what is being compared in Figure 6, namely the average of five θ sensors installed at a single point with an estimate of θ over an area of 1-km 2 . For this reason, it is not reasonable to expect that the model and measured data would ever be in complete agreement.
Ref. [68] compared the AMSR2 soil moisture product with the NRCS SCAN stations in Puerto Rico. They used two down-scaling techniques to improve the comparisons relative to the raw AMSR2 soil moisture data. The mean error for the Juana Diaz Experimental station was −7.5% for the AMSR2 model vs. −2.13% for GOES-PRWEB in this study ( Figure 6). It should be noted that [68] filtered the soil moisture data not using days when the difference between the satellite and station soil moisture values were greater than +/− 0.15. Furthermore, [70] developed an algorithm to estimate soil moisture using a self-organized artificial neural network (ANN) and a stochastic transfer function model for estimating soil moisture at 20 cm depth. Testing the model in Puerto Rico showed that the monthly model had a mean error equal to 2.72%, and the hourly model had mean errors averaging −2.49%. The results from GOES-PRWEB are comparable or better than from the two referenced studies.
For the comparison of monthly stream flow in the Guanajibo Watershed, the mean error for the 36-month period was −2.6%. According to [71], a mean error in stream flow of less than 10% is very good. In addition, [72] calibrated the Vflo model ( [73]) for a one-year period (2003) for the same watershed and obtained a mean monthly error 1.8%. Moreover, [74] calibrated a MIKESHE model ( [75]), which included surface and groundwater flow for the Añasco, Yagüez and Guanajibo watersheds, and obtained a mean monthly error of 11.7%. Finally, [19] compared stream flow from several models, including GOES-PRWEB, with an improved version of the Water Supply Stress Index (WaSSI) water balance model for the El Yunque National Forest in northeastern Puerto Rico. Unfortunately, the authors incorrectly compared GOES-PRWEB RO with their streamflow instead of RO + DP, which consequently invalidates the comparison.
In this study, a lumped method for the surface runoff for the 1-km 2 pixel assumes single values of CN and rainfall. In this case, the CN is based on the predominant land use, soil texture, and hydrologic soil group, all of which may have significant spatial variation in the real world. Similarly, in Puerto Rico, rainfall is highly variable spatially and may not be constant within the pixel area. Furthermore, the lead author's unpublished data collected in Puerto Rico, comparing the AHPS rainfall with rain gage data, indicates that significant rainfall estimation errors can occur. Another possible source of error is the rating curve that the USGS uses to convert the stream depth to stream discharge. Over time, the channel cross-section can change, thus introducing errors into the rating curve. According to the USGS, this was an extensive problem across Puerto Rico, especially after Hurricane María in September 2017 (personal communication, David Hernandez, USGS Hydrologic Data Chief, PR, 16 July 2018).
The DP term in GOES-PRWEB is gross recharge, with a portion potentially going to deep aquifers recharge and a portion, probably the majority, moving laterally through shallow soil or bedrock or to shallow water table aquifers and then discharging to nearby streams. The portion of the DP that flows into the deep aquifers subsequently discharges to the ocean. In the basin-scale stream flow analysis, we assumed that all of DP discharged to the Guanajibo River. This assumption may help explain the overestimation of streamflow for months 5, 17, 33, and 34. On these occasions, it is possible that a larger fraction of the DP went to deep aquifer recharge, and that not all of DP entered the river. It is also noted that the model tended to underestimate streamflow during periods of low flow, which may be because the model does not simulate continuous base flow (i.e., DP); when there is no rainfall, there is no DP or RO. For this reason, the model is not able to estimate daily streamflow.
The basin-scale ET a analysis yielded a mean annual error for the period 2009-2020 of −6.8%. The analysis compared the annual total ET a from GOES-PRWEB with an estimate of AHPS rainfall minus the USGS-measured stream flow. This is consistent with a study in the El Yunque National Forest in northeastern Puerto Rico. Using an improved version of the WaSSI model, [19] obtained a relative error of 7% for GOES-PRWEB annual ET a , relative to their model, for the period 2009-2014. The island-scale annual ET a analysis for the period 2009-2020 yielded a mean difference between SSEBop and GOES-PRWEB of −0.73%. The pixel-scale ET a comparison analysis (same period, same models) for the three locations yielded a mean difference of −1.2%. The method described in this paper uses a lumped composite (one-source) flux surface, as opposed to the two source models [76] and [37], which separates the flux surface into soil and vegetated areas (two-source). The two-source method has been shown to provide more accurate results under certain land covers [76]. Future studies should consider comparing the GOES-PRWEB ET a with the Dis-ALEXI approach. It is also critically important to conduct comparison studies with flux towers in Puerto Rico.
A validation study was also conducted for reference evapotranspiration. ET o estimates were compared with ET o derived from a NRCS SCAN weather station at the UPR Agricultural Experiment Station in Juana Diaz, Puerto Rico. The mean error for the 27-month comparison was 0.96%. Furthermore, [42] compared estimated and weighing lysimetermeasured ET o at Tempe, Arizona for a 26-day period. Estimated ET o was obtained from the generalized Penman-Monteith equation with a constant r s value of 45.6 s m −1 . ET o was measured from a well-watered alfalfa crop. The mean error of the estimation procedure was 1.15%.
The island-wide water budget component analysis showed good agreement between the USGS and the GOES-PRWEB results. The USGS analysis was based on long-term information before 1990, while the GOES-PRWEB analysis was based on conditions for 2009-2020. Molina-Rivera [62] did not provide details associated with the USGS water balance analysis. Table 5 provides estimates of the standard deviations of components of the water balance for the GOES-PRWEB analysis. The standard deviation provides a measure of the variability of the elements of the water balance. Note that there was more variability in the RO+DP data than in the ET a data. During wet years (e.g., 2010, 2011, 2016, and 2017), the RO+DP component increased relative to drier years, whereas the ET a did not increase substantially. The smallest RO + DP occurred in 2009, 2015, and 2019 (446 mm, 407 mm, and 407 mm, respectively). Interestingly, rainfall was significantly greater in 2009 and 2019 than in 2015, when Puerto Rico sustained a severe drought [66]. The result suggests that mean annual rainfall is not necessarily a precise indicator of the amount of RO + DP that we can expect.

Some Model Limitations
In this section, several model limitations are discussed.

•
The model ET o and water and energy balance results are based on weather data (T a , T d , u 2 ) obtained from NDFD or CARICOOS WRF model, Rs from the GOES satellite algorithm, and rainfall from NOAA's AHPS. The main advantage of deriving the input data from these sources is that it is gridded data and is readily assimilated into the model, unlike weather station data which requires interpolation of weather variables. The disadvantage of the gridded data is that it is produced from models and is subject to errors. Furthermore, these data sources are not always available. For example, during Hurricane María, the Doppler Radar (NEXRAD) in Cayey, Puerto Rico, was severely damaged and was not available for nine months. • All water that infiltrates into the soil that exceeds the field capacity becomes DP. This is based on the concept of a "field capacity" and that all water in excess of the field capacity moisture content will percolate past the root zone. The field capacity concept simplifies the soil profile, assuming homogeneous texture and that all of the soil water between the total porosity and field capacity drains within 24 h. Although this may introduce potential errors, as mentioned above (not to mention that an incorrect value of the field capacity could be used), the encouraging results obtained in this study suggest that using the field capacity concept is functionally valid. • Throughout its 12-year life, the model has undergone periodic modifications, ranging from improvements to specific algorithms, changes in input data sources, and adjustment of parameters. Ideally, the model should be reevaluated for ET a , soil moisture, and streamflow after any modifications; however, this is difficult to achieve in practice.

•
The 1-km spatial resolution of the model does not permit estimation for farm-scale conditions. Nevertheless, as shown in [64], the model can be used for specific applications, such as irrigation scheduling, if site-specific information is available.

•
The model is limited to a 1-day time step and, therefore, precludes applications requiring hourly or shorter time steps.
In this section, several model limitations are discussed.
• Daily results should be used with caution. Although daily data are available for use, results for longer time periods will tend to be more accurate because the negative and positive daily errors tend to cancel each other out.

•
In its current form, the model is not capable of estimating snowmelt, which would be required for use at locations in the upper latitudes.

Advantages of GOES-PRWEB
This paper describes a method for estimating the daily water and energy balance (1-km spatial resolution). The operational model is the first of its kind in Puerto Rico. In general, the results indicate that GOES-PRWEB can provide accurate results, which can be used to solve various types of practical problems. The model is already being used extensively for drought monitoring. In agriculture, agronomists can use the weather variables, effective rainfall, crop stress factor, photosynthetically active radiation, and other products to evaluate crop production. Hydrologists and engineers can use the water balance products (ET a , RO, DP, θ) to assess trends in water resources. The solar radiation product can be used to support the design of photovoltaic energy systems. Public health professionals can use the various weather variables to assess climate-related stresses, such as elevated heat indices. Meteorologists can use the energy balance components of the model, including effective surface temperature and Bowen ratio, for analyzing near-surface energy processes. Over time, the model could become a platform for evaluating trends in hydro-meteorological variables related to climate change.
The use of GOES-PRWEB has several advantages over other models. The model includes a suite of over 25 agro-hydro-meteorological variables that are available to the public daily. The authors are unaware of any other model that provides such a wide range of capabilities at the 1-km resolution/daily time scale. The output formats include jpeg images, CSV, and Mathlab ® formats, which are familiar to students, scientists, and professionals. The algorithm runs on a high-performance (gaming) desktop computer, and, therefore, would be practical and economical to use in any region of the world. As noted in Section 4.1, the model data is used by various agencies to support their efforts. It is anticipated that the model will continue to evolve to include more capabilities. Collaboration with other disciplines may yield valuable new products. These products might include indices that help mitigate the severity of forest fires, disease outbreaks from mosquitoes, cases of heat stress during heat waves, drought impacts, pest infestation of crops, etc. It is hoped that the model can be deployed for use in the islands surrounding Puerto Rico in the near term.