## 1 Introduction
### Reference Parameters and Datasets  
To achieve the work's objectives, we evaluated evapotranspiration (ET), land surface temperature (LST), normalized difference vegetation index (NDVI), and gross primary production (GGP) using remote sensing datasets from 2012 to 2021. ArcGIS was used to perform most of the geoprocessing; however, data analysis was also performed using RStudio. Attribute tables for each map were created to build correlation graphs, and each table was created as a spreadsheet using the coordinates of the datapoints within the study area for it to be compatible with RStudio. Spreadsheets were created using latitude and longitude coordinates along with datapoints (approximately 2500 points). We then used the Create Fishnet function on ArcMap to place/plot the datapoints onto the study area. Then, the values of the points were extracted using the Extract Multi Values to Points tool. For each attribute table, we used a 50 \(\times\) 50 fishnet. When \"clipped\" to the study area, a total of approximately 1349 points was used. For a more accurate estimation of ET and LST values, adequate attribute tables were built, as gaps within these datasets had to be amended using other sources that over- or underestimated values. To estimate the resilience of the ecosystem to atmospheric events, we compared and correlated the estimations for ET, NDVI, LST, and GPP, calculated from Landsat and MODIS. We used Landsat 8-9 OLI TRIS for LST, and MODIS products were used for ET, GPP, and NDVI. To compare the data from MODIS products, multiday values were modified into daily values. Gap-filled GPP maps were converted from an 8-day composite to daily values, and the same process was performed for MOD16A2 for ET estimations. The coordinate system for MODIS products is undefined; therefore, we defined the coordinate system using WGS 1984 UTM Zone 20N. Figure 2: Spatial distribution of imperviousness within the Culebrinas River Watershed in Western Puerto Rico. The urban areas can be identified in red clusters as part of the main towns of the municipalities of Aguadilla, Aguada, Moca, and San Sebastian. To calculate ET, we used MOD16A2 at a 500 m, 16-bit spatial resolution [20]. It has a valid range from \(-3276.1\) to \(3276.7\):  
\[\text{ET}=\text{MOD16A2-ET}^{*}0.1 \tag{1}\]  
where MOD16A-ET represents the values obtained from MODIS products for ET estimations. The 0.1 is a scale factor provided by the Land Processes Distributed Active Archive Center (LP DAAC) for ET estimations. To obtain the daily values of ET, we converted the datasets from an 8-day composite to daily ET estimations. During atmospheric events, heavily clouded zones formed in our area of study, generating uncertainty within the datasets. For ET estimations, we used GOES-PR data to amend the MOD16A2 datasets. The unit kg/m\({}^{2}\)/day is equivalent to mm/day. Thus, MODIS-derived ET is comparable/compatible with GOES-PR data. As GOES-PR has a 1 km resolution, the points plotted from GOES-PR did not fully cover the gaps. Thus, we used inverse distance weighted interpolation (IDW) to finalize amending the zones with cloud interference. IDW is a deterministic spatial interpolation that assumes that the closer the value of the interpolation is to the initial point, the more accurate it is. In general, the interpolated values from GOES-PR underestimated ET values as compared with the values obtained from MOD16A2. [PERSON] and [PERSON] [21] presented a procedure to estimate LST from Landsat 8-9 using emissivity estimation (Equation (8)). Emissivity estimations for temperature generate a grade of uncertainty to the final LST estimation, and as per [PERSON] [22], emissivity estimations from Landsat 8-9 present good estimations for LST. For this procedure, Band 10 was used to determine brightness temperature, and Bands 4 and 5 were used for the NDVI. Other values were taken from the metadata of the maps (Table 1).