Forest Cover and Locality Regulate Response of Watershed Discharge by Puerto Rico Sea Grant

QiongGaoandMeiYu*

Citation: Gao,Q.;Yu,M.ForestCover andLocalityRegulateResponseof WatershedDischargetoRainfall VariabilityinCaribbeanRegion. Forests 2024, 15,154. https:// doi.org/10.3390/f15010154

AcademicEditor:MatthewTherrell

Received:22November2023

Revised:6January2024

Accepted:9January2024

Published:11January2024

Copyright: © 2024bytheauthors. LicenseeMDPI,Basel,Switzerland. Thisarticleisanopenaccessarticle distributedunderthetermsand conditionsoftheCreativeCommons Attribution(CCBY)license(https:// creativecommons.org/licenses/by/ 4.0/). Article ForestCoverandLocalityRegulateResponseofWatershed DischargetoRainfallVariabilityinCaribbeanRegion

DepartmentofEnvironmentalSciences,UniversityofPuertoRico,RioPiedras,SanJuan,PR00925,USA * Correspondence:meiyu@ites.upr.edu

Abstract: Reforestationoftenoccurswhentheeconomyshiftsfromagriculturetoindustryand servicessuchastourism.However,thereisalackofcoherentknowledgeandinvestigationaboutthe impactofreforestationinthetropicsonhydrologicalvariabilityaswellasfloodrisks.Itisunclear howchangesinforestcoverandpatternwillaffectfloodrisksandwatershedresponsetofuture alteredrainfallintensity.ThisstudyusestheSoilWaterAssessmentTool(SWAT+)tosimulatethe impactofreforestation,thelocalityofforest,andtheconcentratedrainfallonthehydrologyofthe largestwatershedinPuertoRico.SWAT+isacomputermodelsimulatingwatershedhydrology drivenbymeteorologicalinputandthecharacteristicsofsoilsandlanduse.Wehypothesizedthat increasedforestcover,especiallyatlowelevationrange,wouldreducefloodriskandthatreduced raindayswhilemaintainingthemeanannualrainfallinvariantwouldincreasestreamdischarge variability.Wefoundthatreforestationsignificantlyreducedlargedischargesbutincreasedsmall discharges;thatforestatlowelevationtendedtoreducelargeandextremedischargesincomparison withforestathighelevation;andthatmoreconcentratedrainfallnotonlyincreasedtherainfall variabilitybutalsoincreasedthedischargevariability.However,boththeimpactofshiftingforest localityandtheresponseofwatershedtoalteredrainfallintensitystronglydependedongeophysical factorssuchasrangesofelevationandslope.Movingforeststolowerelevationinsubbasinswith steeperslopesshowedastrongerreductioninextremedischargesthaninsubbasinswithflatter slopes.Ontheotherhand,subbasinswithsteeperslopestendedtoresponsemorestronglytomore concentratedrainfallwithgreaterincreaseindischargevariabilitythansubbasinswithflatterslopes. Tocopewithfutureincreasedclimatevariability,ourresultsfavorreforestationatlowerelevationfor watershedwithlargeelevationrangesandsteepslopes.

Keywords: hydrologicalvariability;PuertoRico;landusepattern;discharges;rainfallvariability

1.Introduction

Climatevariabilityoftenreferstohowthevaluesoftemperatureorprecipitationdiffer fromthemeanvalues[1].Naturalclimatevariabilityarisesfromperturbation-caused positiveandnegativefeedbackswithintheatmosphericsystem[2].However,themain componentofpresentclimatevariabilityisaresultofhumanactivities[3].Warming increasedtheenergyoftheearth’ssystem,andtheprogressivelossofglaciersandice capsreducedthecoolingofoceanwater.Theseexternalforcingscausedmorefrequent andextremeoceanicandatmosphericeventssuchaslocalizedhighatmosphericpressure interwovenwithmegatropicalcyclones.Consequently,morefrequentextremeweather conditionsimpactcoastalandterrestrialecosystemsaswellashumansocieties[4].Alternatingdroughtandextremelargerainfallbroughtaboutbymajorstormsincreased hydrologicalvariabilityandfloodriskandreducedtheavailabilityofwaterforterrestrial ecosystemsandhumanconsumptions[5].

Notonlydoeshydrologicalvariabilitylinktoclimateprocessesbutitalsoheavilydependsonhowterrestrialecosystemsrespondtotheclimatevariability[6].Thehydrological responsesofterrestrialecosystemtorainfalleventsarecomplex.Soilandvegetationhold

andconsumeaportionofwater,andrunoffisgeneratedwhenrainfallintensityexceeds themaximuminfiltrationrateoflandmatrix.Thedischargeofacatchmentisaresultof accumulatedoverlandandsubsurfaceflows[6].

Terrestrialecosystemoftenactsasabufferorfiltertotherainfallvariability[7].The rootsystemofvegetation,especiallyforest,facilitatestherechargeofsoilandaquifer, whichholdwaterduringlargerainfallandreleasewatergraduallyafter.Surfaceroughness slowsdownoverlandflow,andthehydraulicconductivityofsoilalsogovernsthetiming ofsubsurfaceflow.Ontheotherhand,pavedsurfaceandbarelandlackthisfiltering capacity,andtheconstructionundergroundmayblockthesubsurfaceflow.Thus,the dischargevariabilityofacatchmentdependsnotonlyontherainfallvariabilitybutalso onthecompositionandpatternsoflandcover,especiallytheforestcoverandpattern[8]. Specifically,forestatcriticallocationswithinabasinmayhavemorechancestoregulate waterflowsandtoaffectthevariabilityofstreamdischargesthanatotherlocations[9].

Inthetropics,theislandPuertoRicounderwentprofoundlandusechangesand experiencedincreasingclimatevariability.LocatedattheeastendoftheGreatAntilles andboundedbyAtlanticOceaninthenorthandCaribbeanSeainthesouth,thetropicalislandisgreatlyimpactedbythemorefrequentalternationofseveredroughtand majortropicalstorms[10].TropicalstormsimpacttheislandmostlyduringSeptemberwithextremerainfall[11,12],whereasseveredroughtslimitwatersupplyanddegradetheislandecosystemandharmhumansociety[13].Increasingdroughtsafter2000 (https://www.drought.gov/states/puerto-rico,accessedon1January2024)havebeen observed,andrainydayshavebeenshowntohavedecreased[14]fromthe1950sonward, especiallyinthewetseason.Thepatternsofseveredrought/stormsthatimpactPuerto RicowerefoundcloselyrelatedtoENSO(ElNiñoandSouthernOscillation),NAO(North AtlanticOscillation),andLaNiñaevents[15,16].Majorhurricanesoftenfollowsevere droughts.Forexample,therecentseveredroughtduring2014–2016inPuertoRicoisaligned withthestrongENSOeventinthe2015–2016periodwithapeakwarmingindexof2.2 (https://psl.noaa.gov/enso/mei/,accessedon1January2024).Thisseveredroughtwas followedbyaCategory5hurricaneIrmaandahigh-endCategory4hurricaneMariain 2017,whichdevastatedPuertoRico[10,17].

Industrializationstartedinthe1940sinPuertoRico.Asaresult,urbanexpansionand forestrecoveryfollowedtheabandonmentofagriculture[18].Accordingtoavailableland covermaps[19,20],forestgrewfrom42%ofthelandto55%,andurbanexpandedfrom 14%to17%,inthe10yearperiodfrom2000to2010.Theforestwasmorefragmentedin mountains,especiallyinthewestoftheisland.Thedemandfortourismatthecoastoften competeswithforestsforspace,whichchallengesmangroveforestconservation.Land coverchangenotonlymodifiedthemeandischargebutalsoalteredfloodrisk[21].Itis unclearhowthevariabilityinrainfallislinkedtothedischargevariabilityofriversunder alteredlanduse/cover.

Thisstudyaimedatanoverarchingquestionhowreforestationandforestpatterns modulatethewatershedhydrologyanditsresponsetoalteredrainfallvariability.Specifically:howdoesreforestationaffectthelargedischargeswithfloodrisks;howdoesthe localityorpatternofforestaffectthelargedischarges;andhowdoesincreasingrainintensityalterthedischargevariabilityintropicalwatersheds?Wehypothesizedthatthe increasedforestcoverandforestallocatedatlowelevationtendtoreducethelargedischargesandthatthereducedraindayswithincreasedrainintensitywouldincreasethe variabilitiesofbothrainfallanddischargeinstreams.Weappliedthehydrologicalmodel SWAT+(SoilWaterAssessmentTool)tothewatershedofRioGrandedeLoiza,thebiggest watershedinPuertoRico.SWAT+isacomputermodelaimedtosimulatewatershedhydrologicalprocesseswithinputofspatialpatternsofsoilandlanduseaswellasmanagerial plans.Afterthecalibrationandvalidationofthemodel,wesimulatedthedischargerate withtwolanduse/covermapsof1991and2010.Wethengeneratedtwohypotheticalland usepatternswithforestallocatedtothehighestandthelowestelevation,respectively,and ranSWAT+toinvestigatehowforestlocalityaffectsthedischarges.Asthethirdpartof

thestudy,wegeneratedrainfallwithreducedraindaysbutmaintainedthemeanannual rainfallofthestationsinthewatershedtoinvestigatehowincreasedrainfallintensity wouldaffectthevariabilityofstreamdischarges.

rate with two land use/cover maps of 1991 and 2010. We then generated two hypothetical land use patterns with forest allocated to the highest and the lowest elevation, respectively, and ran SWAT+ to investigate how forest locality aﬀects the discharges. As the third part of the study, we generated rainfall with reduced rain days but maintained the mean annual rainfall of the stations in the watershed to investigate how increased rainfall intensity would aﬀect the variability of stream discharges.

2.MaterialsandMethods

2. Materials and Methods

The Rio Grande de Loiza in Puerto Rico has the greatest drainage basin of 72,688 ha (Figure 1). Located in eastern Puerto Rico, the mountainous part of the watershed intercepts a great amount of moisture carried by the eastly trade wind. The annual rainfall at the rain gauge of RQC00668815 (200 m a.s.l.) is 2402 mm (1981–2010), the greatest among the seven stations within the watershed. The wetland near the river mouth is mostly occupied by mangroves. A reservoir, Lake Loíza, is located at the lower reach of the river down the city of Caguas with the dam built in 1953. The original capacity of the reservoir is 26.8 Mm3 According to the report in 2011, the reservoir water was diverted and supplied to San Juan metropolitan area with a mean rate of 394,000 m3 d 1 [22]. Other water extractions from the lake and streams are unknown.

TheRioGrandedeLoizainPuertoRicohasthegreatestdrainagebasinof72,688ha (Figure 1).LocatedineasternPuertoRico,themountainouspartofthewatershedintercepts agreatamountofmoisturecarriedbytheeastlytradewind.Theannualrainfallattherain gaugeofRQC00668815(200ma.s.l.)is2402mm(1981–2010),thegreatestamongtheseven stationswithinthewatershed.Thewetlandneartherivermouthismostlyoccupiedby mangroves.Areservoir,LakeLoíza,islocatedatthelowerreachoftheriverdownthecity ofCaguaswiththedambuiltin1953.Theoriginalcapacityofthereservoiris26.8mm3 Accordingtothereportin2011,thereservoirwaterwasdivertedandsuppliedtoSanJuan metropolitanareawithameanrateof394,000m3 d 1 [22].Otherwaterextractionsfrom thelakeandstreamsareunknown.

Figure 1. Location of the watershed of Rio Grande de Loiza in the Caribbean (A) and Puerto Rico (B); distribution of DEM, streams, reservoir, and gauges of discharge within the watershed (C); major subbasins and codes for discharge gauges (D). 50055000 (Rio Grande de Loiza at Caguas) and 50057000 (Rio Gurabo at Gurabo) are two upper-reach subbasins, and 50051800 (Rio Grande de Loiza at Highway 183) is a subbasin of 50055000. Maps were created in ArcGIS Pro 3.2 (ESRI, Redlands, CA, USA).

Figure1. LocationofthewatershedofRioGrandedeLoizaintheCaribbean(A)andPuertoRico (B);distributionofDEM,streams,reservoir,andgaugesofdischargewithinthewatershed(C);major subbasinsandcodesfordischargegauges(D).50055000(RioGrandedeLoizaatCaguas)and50057000 (RioGuraboatGurabo)aretwoupper-reachsubbasins,and50051800(RioGrandedeLoizaatHighway 183)isasubbasinof50055000.MapswerecreatedinArcGISPro3.2(ESRI,Redlands,CA,USA).

The Soil Water Assessment Tool (SWAT+) is a hydrological model built as a plugin to QGIS to simulate the hydrological processes of a watershed by computing evapotranspiration, surface runoff, subsurface flow, sediment yield, and nutrient dynamics driven by daily meteorological records [23,24]. Hydrological response unit (HRU), the basic spatial unit of SWAT+, was formed as the result of overlaying the delineated subbasins, soil map, land use map, and slope ranges. Therefore, land use composition and configuration as well as geophysical factors such as elevation and slope interplay to determine watershed stream discharge. Specifically, we used the QSWAT+ 2.0.6 and SWAT+ model of 60.5.2 associated with QGIS 3.16 (QGIS Geographic Information System. Open Source Geospatial Foundation Project. http://qgis.org, accessed on 1 January 2024). The model was set up

TheSoilWaterAssessmentTool(SWAT+)isahydrologicalmodelbuiltasapluginto QGIStosimulatethehydrologicalprocessesofawatershedbycomputingevapotranspiration,surfacerunoff,subsurfaceflow,sedimentyield,andnutrientdynamicsdrivenby dailymeteorologicalrecords[23,24].Hydrologicalresponseunit(HRU),thebasicspatial unitofSWAT+,wasformedastheresultofoverlayingthedelineatedsubbasins,soilmap, landusemap,andsloperanges.Therefore,landusecompositionandconfigurationas wellasgeophysicalfactorssuchaselevationandslopeinterplaytodeterminewatershed streamdischarge.Specifically,weusedtheQSWAT+2.0.6andSWAT+modelof60.5.2 associatedwithQGIS3.16(QGISGeographicInformationSystem.OpenSourceGeospatial FoundationProject. http://qgis.org,accessedon1January2024).Themodelwasset upwithaDEMof20mresolution.ThesoilmapofPuertoRicoSoilSurveyGeographic Database(SSURGO)andlanduse/covermapsof1991and2010wereusedtogenerate HRU.TheresultofthewatersheddelineationandHRUgenerationgave40subbasins, 75channels,and875HRUsbasedonthelanduse/covermapof1991.

TorunSWAT+forRioGrandeLoizaWatershed,weusedeightweatherstationswithin andinclosedistancetothewatershed.ThedataweredownloadedfromNationalClimate DataCenter,whichincludedRQC00661590(CANOVANAS),RQC00669521(TRUJILLO ALTO2SSW),RQC00664276(GURABOSUBSTATION),RQC00661309(GAGUAS1W), RQC00665064(JUNCOS1SE),RQC00668822(SANLORENZOFARM2N),RQC00668815 (SANLORENZO3S),andRQC00664115(GUAVATECAMP).Oftheeightstations, RQC00661590,RQC00669521,RQC00665064,andRQC00664276havedailymaximum andminimumtemperatureobservations.Themissingvaluesofrainfallwerefilledby meansofscalingvaluesfromstationswithmaximumcorrelation,amethodsimilartothat intheR4.2packageof‘fillGap’(RFoundationforStatisticalComputing,Vienna,Austria). Themissingvaluesoftemperaturewerefilledwithregressionontime(withmultiple harmonicwaves)aswellastherainfallintheprevious2days.

Thedailydischargesfrom4USGSgauges(Figure 1)ofRioGrandedeLoizaatHighway183SanLorenzo(USGScode50051800),RioGrandedeLoizaatCaguas(50055000), RioGrandedeLoizabelowLoizadamsite(50059050),andRioGuraboatGurabo(50057000) wereusedforcalibration,validation,andsubsequentanalyses.Thedischargesat50055000 and50057000werefromthetwomajorsubbasinsintheupperreachesofthewatershed (Figure 1).Wehereafterrefertothesesubbasinsusingtheircorrespondinggaugecodes. SWAT+modelmustbecalibratedandvalidatedbeforesimulation,andweusedthe dailyclimateanddischargedatafrom1987to1991tocalibratethemodel.Thecalibration wascarriedoutbymeansofSWAT+Toolbox.Fourparameterswerechosentobeoptimized: cn2,thecurvenumbercontrollingtheproportionofsurfacerunoffduringarainfallevent; k, thesaturatedsoilhydraulicconductivity; flo_min,thethresholdgroundwaterleveltoallow baseflow;and alpha,theparametercontrollingbaseflowtiming.Tovalidatethemodel,we usedmonthlydischargesfrom1981to2010forallfourgauges.

Toinvestigatetheimpactoflandcoverchangesonthevariabilityofstreamdischarges, wesimulatedthedischargesofthethreesubbasinsof50055000,50057000,and50059050 withthelanduse/covermapsof1991[25]andNOAAC-CAPfor2010[19],respectively. Thecalibratedmodelwasrunforeachofthelanduse/covermaps.Weusedthestandard deviationofdailydischargeasthemeasureofdischargevariability,andtheimpactofthe landusechangecanbederivedbycomparingdischargesbetweenlandusemaps.

Twohypotheticallandusescenariosweregeneratedbymanipulatingforestand pasturecoversbasedonthe1991landusemapwhilemaintainingotherlandusetypes unchangedinthetwoupper-reachsubbasinsof50055000and50057000.The highforest scenariowascreatedbymovingtheforesttotopelevationbutpasturetolowelevation, andthe lowforest scenariowasachievedbymovingforesttolowelevationandpastureto topelevation.ThecalibratedSWAT+wasrunforthesetwoscenariostoinvestigatehow theforestlocalityaffectsthestreamdischarges.

Thescenariosofreducedraindaysweregeneratedwithamethodsimilartotherainfall generatorembeddedinSWAT+.Theconditionalprobabilitiesofraintodaygivenyesterday wet,andraintodaygivenyesterdaydrywerecomputedfromtherainfallsequenceforeach ofthe12months.ThedailyrainfallamountwasfittedtoaWeibullminimumdistribution. Theprobabilityofraingivenyesterdaywetwasreducedbydifferentpercentagesand thescaleparameteroftheWeibullminimumdistributionwasincreasedaccordinglyto maintainmeanannualrainfalllargelyunchanged.Thealteredrainfallscenarioswerethen appliedtosimulatewatersheddischarges.QGISwasusedforspatialdatapreparation. Rainfallanalyses,scenariogeneration,anddischargeanalyseswerecarriedoutwithPython 3.10andR4.2[26].

3.Results

3.1.CalibrationandValidationofSWAT+

ThecalibrationofSWAT+yieldedtheNash–Sutcliffemodelefficiencycoefficient (NSE)of0.532fordailydischargesfrom1987to1991atthegauge50059050.Thevalidation ofmodel(Figure 2)wasperformedbycomparingthesimulatedandobservedmonthly

3. Results

3.1. Calibration and Validation of SWAT+

The calibration of SWAT+ yielded the Nash–Sutcliffe model efficiency coefficient (NSE) of 0.532 for daily discharges from 1987 to 1991 at the gauge 50059050. The validation of model (Figure 2) was performed by comparing the simulated and observed monthly discharges from all four gauges during 1981–2010. The NSE varied from 0.58 for 50051800 to 0.83 for 50057000.

dischargesfromallfourgaugesduring1981–2010.TheNSEvariedfrom0.58for50051800 to0.83for50057000.

Figure 2. Validation of SWAT+ for monthly discharges at four stream gauges during 1981–2010. NSE, Nash–Sutcliﬀe eﬃciency.

Figure2. ValidationofSWAT+formonthlydischargesatfourstreamgaugesduring1981–2010.NSE, Nash–Sutcliffeefficiency.

3.2.ImpactofLandUseChangeonDailyDischarges

3.2. Impact of Land Use Change on Daily Discharges

Land use/cover of the three subbasins of 50055000, 50057000, and 50059050 featured reforestation from 1991 to 2010. Forest mostly emerged from the previous grassland. While agricultural land use remained relatively small and constant, urbanization and urban sprawl continued (Table 1).

Landuse/coverofthethreesubbasinsof50055000,50057000,and50059050featured reforestationfrom1991to2010.Forestmostlyemergedfromthepreviousgrassland.While agriculturallanduseremainedrelativelysmallandconstant,urbanizationandurban sprawlcontinued(Table 1).

Table 1. Land cover fractions within three subbasins of the Rio Grande de Loiza watershed.

Table1. LandcoverfractionswithinthreesubbasinsoftheRioGrandedeLoizawatershed. WatershedYearForestPastureAgricultureUrban

The changes in land use altered the hydrology of the watershed. Simulated daily discharges of the three subbasins (Figure 3 upper panels) showed the prevailing patterns of decreased large discharges with the land use of 2010 compared to those with the land use of 1991.

Thechangesinlandusealteredthehydrologyofthewatershed.Simulateddaily dischargesofthethreesubbasins(Figure 3 upperpanels)showedtheprevailingpatternsof decreasedlargedischargeswiththelanduseof2010comparedtothosewiththelanduse of1991.

Bootstrappingwasappliedtoregressthedifferenceinsimulateddailydischarge between2010and1991landusemapsonthelog-transformeddailydischargeswithland useof1991(Table 2).Theslopesoftheregressionaresignificantlysmallerthanzeroasthe 95%confidenceintervalsareinthenegativerange,implyingtheprevailingreducedbig dischargeswiththelanduseof2010.Specifically,themeansloperangesfrom 0.082for 50057000(RioGuraboatGurabo)to 0.349for50055000(RioGrandedeLoizaatCaguas). Themeanslopeof 0.323for50059050(RioGrandedeLoizabelowthedam)islargelya combinationofthetwoupper-reachsubbasins.However,theinterceptsaresignificantly greaterthanzero,i.e.,0.017for50057000to0.133for50059050,implyingtheincreasedsmall discharges(baseflow)withthelanduseof2010.Thesuperpositionoftheincreasedsmall butreducedlargedischargemakesthemeandischargeswiththelanduseof2010slightly decreasedcomparedtothosewiththelanduseof1991(5.24vs.5.29m3 s 1 for50059050).

Figure 3. For three subbasins of the Rio Grande de Loiza watershed, differences in daily discharges between 2010 and 1991 land use versus the log-transformed daily discharges based on 1991 land use (upper panels, land use), and differences in daily discharges between low forest and high forest scenarios versus the log-transformed daily discharges with the high forest scenario (lower panels, forest locality). FL, FH—low forest, high forest scenarios. The horizontal red line shows zero difference in discharge. Points below the lines indicate discharge reduction, whereas those above the lines show increases in discharge.

Figure3. ForthreesubbasinsoftheRioGrandedeLoizawatershed,differencesindailydischarges between2010and1991landuseversusthelog-transformeddailydischargesbasedon1991land use(upperpanels,landuse),anddifferencesindailydischargesbetween lowforest and highforest scenariosversusthelog-transformeddailydischargeswiththe highforest scenario(lowerpanels, forestlocality).FL,FH—lowforest, highforest scenarios.Thehorizontalredlineshowszerodifference indischarge.Pointsbelowthelinesindicatedischargereduction,whereasthoseabovethelinesshow increasesindischarge.

Table2. Coefficientsofbootstrappingregressionofdifferenceindailydischargebetween2010and 1991landuseonlog-transformeddailydischargeswithlanduseof1991(upper),andcoefficientsof bootstrappingregressionofdifferenceindailydischargebetween lowforest and highforest landuse scenariosonlog-transformeddailydischargeswiththe highforest scenario.

SubbasinSlopeIntercept

Bootstrapping was applied to regress the difference in simulated daily discharge between 2010 and 1991 land use maps on the log-transformed daily discharges with land use of 1991 (Table 2). The slopes of the regression are significantly smaller than zero as the 95% confidence intervals are in the negative range, implying the prevailing reduced big discharges with the land use of 2010. Specifically, the mean slope ranges from −0.082 for 50057000 (Rio Gurabo at Gurabo) to −0.349 for 50055000 (Rio Grande de Loiza at Caguas). The mean slope of −0.323 for 50059050 (Rio Grande de Loiza below the dam) is largely a combination of the two upper-reach subbasins. However, the intercepts are significantly greater than zero, i.e., 0.017 for 50057000 to 0.133 for 50059050, implying the increased small discharges (baseflow) with the land use of 2010. The superposition of the increased small but reduced large discharge makes the mean discharges with the land use of 2010 slightly decreased compared to those with the land use of 1991 (5.24 vs. 5.29 m3 s 1 for 50059050).

Differencebetween2010and1991landuse

Table 2. Coefficients of bootstrapping regression of difference in daily discharge between 2010 and 1991 land use on log-transformed daily discharges with land use of 1991 (upper), and coefficients of bootstrapping regression of difference in daily discharge between low forest and high forest land use scenarios on log-transformed daily discharges with the high forest scenario.

Differencebetweenlowforestandhighforest

50055000

Subbasin Slope

Interval Diﬀerence between 2010 and 1991 land use

50055000 −

3.3.ImpactsofForestLocalityonDailyDischarges

Thesimulateddifferenceindailydischargebetween highforest and lowforest scenarios (Figure 3 lowerpanel)showedthattheforestatlowelevationtendstoreducethelarge discharges.Thepatternisespeciallyobviousforthedischargesat50055000and,subsequently,at50059050butisweakandhardlyobservableat50057000.Aftertheremovalof the‘outliers’withcriteriaforthedifferenceindischargesmallerthan 40for50055000and 50057000and 10for50059050,bootstrappingwasappliedtoregressthedifferenceinthe dischargebetweenthetwoforestlocalityscenariosonthelog-transformeddischargewith thehighforestscenario.

Thecoefficientsofthebootstrappingregression(Table 2)showthattheslopefor 50055000is 0.943 ± 0.155,withthe95%confidenceintervalof[ 1.239, 0.624].Since theconfidenceintervalisnegative,wecaninferthattheslopeissignificantlylessthan0. Thesameargumentappliestotheslopeof50059050,forwhichtheslopeof 0.265 ± 0.048 isalsosignificantlylessthan0.Thus,thelargedischargeissignificantlyreducedwith forestatlowelevation.However,thispatternisnottruefor50057000,wheretheslopeof

0.007 ± 0.123with95%confidenceintervalof[ 0.235,0.262]isnotsignificantlydifferent fromzero.

3.4.ResponseofStreamDischargestoAlteredRainfallVariability

TheconditionalprobabilitiesofrainandtheparametersofWeibullminimumdistributionforeachofthe12monthswereshowninTable 3 forthemeteorologicalstation RQC00668815during1981–1990(Table 3).Thesamesetofparameterswerederivedfor eachstationintheperiodsfrom1991to2000andfrom2001to2010.Themeanmonthly rainfallwasregressedontheconditionalprobabilitiesofraininthedaygivenwetordry inthedaybeforeaswellasthescaleofWeibullminimumdistribution.Theequationfor RQC00668815during1981–1990isasthefollowing, P = 281.4 + 24.06λ + 277.24pw + 59.89pd, R2 = 0.990(1) where P isthemeanmonthlyrainfall, λ isthescaleoftheWeibullminimumdistribution, and pw and pd aretheprobabilitiesofraingivenyesterdaywetanddry,respectively.The R2 fortheeightstationswasfoundintherangefrom0.952to0.990withameanof0.980. ThemodelfitwasshowninFigure 4.Ingeneratingrainfallscenarios,wedecreasedthe pw by10to60%atastepof10%.Tokeepthemonthlyrainfallrelativestable,weincreasedthe λ accordingly.Specifically,

where

and

arethecoefficientsbefore pw and λ,respectively.

Table3. TherainfallparametersforRQC00668815during1981–1990. k and λ aretheshapeandscale ofWeibullminimumdistributions,respectively. pw and pd aretheconditionalprobabilitiesofrain todaygivenyesterdaywetanddry,respectively. P isthemeanmonthlyprecipitation.

January1.1914.7410.9360.378106.8 February0.9865.6430.9090.408120.0 March0.9015.3720.8940.277117.7 April0.8855.6320.8990.306126.0

May0.85510.4230.9350.362275.3 June0.9178.6670.9250.413210.2

July0.9679.5680.9730.467260.5

August0.8647.7740.9390.533216.5

September0.8319.4380.9210.351235.0 October0.8779.4030.9350.383243.3

November0.7248.0850.9150.537252.8

December0.8517.0420.9610.423204.3

Figure4. ThefitofmonthlyrainfallonthescaleoftheWeibullminimumdistributionandthe conditionalprobabilityofraininadaygiventheconditionofthedaybeforefortwostationsinthe periodof1981–1990.

Figure 4. The ﬁt of monthly rainfall on the scale of the Weibull minimum distribution and the conditional probability of rain in a day given the condition of the day before for two stations in the period of 1981–1990.

The daily rainfall sequences were generated for each station in the period of 1980 to 2010. Figure 5 shows the rainfall sequence generated for RQC00668815 in the period of three years from 2006 to 2008. The altered rainfall scenarios with the reduction in continuous rain days showed that the rainfall variability increased (Figure 6), whereas the mean annual rainfall was largely unchanged.

daysshowedthattherainfallvariabilityincreased(Figure 6),whereasthemeanannual rainfallwaslargelyunchanged.

Figure 5. Simulated daily rainfall for RQC00668815 in the period of 1 January 2006–31 December 2008 with various reductions in continuous rainfall probability (RCR).

Figure5. SimulateddailyrainfallforRQC00668815intheperiodof1January2006–31December 2008withvariousreductionsincontinuousrainfallprobability(RCR).

Figure6. Variabilityofgeneratedrainfallwiththreelevelsofreductionincontinuousrainfall probability(0,30,and60%)andcorrespondingdischargevariabilityatthetwoupper-reachsubbasins.

Figure 6. Variability of generated rainfall with three levels of reduction in continuous rainfall probability (0, 30, and 60%) and corresponding discharge variability at the two upper-reach subbasins.

SWAT+ simulated discharges at 50055000, 50057000, and 50059050 corresponding to the reduction in continuous rain days (Figures 6 and 7). However, the intercept and slope of the regressions of discharge variability on the reduction in the probability of continuous rain days are quite different among the three subbasins. Rio Gurabo at Gurabo has the greatest intercept (8.54 m3 s 1), whereas Rio Grande de Loiza at Caguas has the steepest slope (0.044). The downstream outlet below the dam showed the smallest intercept and slope.

SWAT+simulateddischargesat50055000,50057000,and50059050corresponding tothereductionincontinuousraindays(Figures 6 and 7).However,theinterceptand slopeoftheregressionsofdischargevariabilityonthereductionintheprobabilityof continuousraindaysarequitedifferentamongthethreesubbasins.RioGuraboatGurabo hasthegreatestintercept(8.54m3 s 1),whereasRioGrandedeLoizaatCaguashasthe steepestslope(0.044).Thedownstreamoutletbelowthedamshowedthesmallestintercept andslope.

Figure 7. Mean discharge variability changes with the reduction in continuous rainfall probability (RCR) for the three subbasins.

Figure7. Meandischargevariabilitychangeswiththereductionincontinuousrainfallprobability (RCR)forthethreesubbasins.

4.Discussion

Landuse/coverchangehasbeenidentifiedasanimportantregulatorofriverdischarges.Vegetation,especiallyforest,tendstoslowdownsurfacerunoff,tendstointercept rainfallvialeavesandstemflows,tendstofacilitatetherechargeofsoils,andtendsto increasebaseflow[27,28].Forestrootsalsoallowsoilstostoremorewaterthanother vegetationsandtendtoreducesurfaceflowbutincreasebaseflow[29],hencereducing thedischargevariability.Ontheotherhand,runoffcoefficientsforimpermeablesurfacesinurbanandtilledagriculturallandaremuchlargerthanthatofforest.Therefore, theimpermeablesurfacesandtilledlandhavelittleorlimitedcapacitytoreducethe hydrologicalvariability.

ReforestationinPuertoRicohaslastedfordecadesaftertheindustrialization[30]. Forestcoverincreasedfrom15.7%in1970sto45.7%in2014[31].TheforestcoverinRio GrandedeLoizawatershedincreasedbymorethan10%overtheperiodfrom1991to2010. Duringthesameperiod,citiesofCaguas,TrujilloAlto,andGurabowithinthewatershed expandedby5~7%despitethestumblingeconomy.Oursimulationshowsthattheeffect ofreforestationdominatesoverthatofurbanexpansionsothatthebigdischargeatthe dailyscaleisreducedforallthethreesubbasins.ThepositiveinterceptsinTable 2 indicate thatthereweregreaterbaseflowswith2010landusethanwith1991landuse.Theability tooffsetthebigdischargediffersamongthreesubbasinsandseemstodependonboth thelandusein1991andthelandusechangefrom1991to2010.RioGrandedeLoizaat Caguas(50055000)hasthelargestnegativemeanslope( 0.349,Table 2);thus,itsabilityto offsetthebigdischargeisthestrongest.Thisisbecausethatthesubbasinhadthehighest forestcoverin1991(46%)andthestrongestincreaseinforestcover(14%)from1991to2010. ThesubbasinofRioGuraboatGurabo(50057000)hastheweakestmeanslope( 0.082) andthustheweakestabilitytooffsetthebigdischargebecauseofitslowestforestcover in1991(34%)andthesmallestincrease(12%)from1991to2010.Moreover,thesubbasin of50057000hadthebiggesturbanexpansionfrom17%to25%.Thedischargefromthe reservoirdamismostlydeterminedbythedischargeatRioGrandedeLoizaatCaguas, whichhasmuchlargerdischargesthanthoseofRioGuraboatGurabo.Thus,theslopeand interceptof50059050(Table 2)areinbetweenandclosetothoseof50055000. Forestcoverisdirectlylinkedtotheconversionofrainfallandrunofftobaseflow.In theinteriortropicalforestofPuertoRico,39to45%ofrainfallwasconvertedintobaseflow instreams[32].Thehigherconversionratioisconnectedtogreaterforestcover.Forest

areahasbeenidentifiedasaone-parametermodelofthevariabilityofstreamflows,and thegreatertheforestarea,themorestablethestreamflows[33].Watershedswithgreater forestcoveralsopossessgreaterbufferingcapacityduringfloodingevents[34].Aprevious simulationstudyofmountainouswatershedsinPuertoRicoshowedthatforestcover significantlydecreasedthebigdischarges[21].Thesynthesisofhydrologicalconsequences ofreforestationincentralAmericaindicatedthatforestcoverregulatedthehydrological processesandecosystemservices[35].Inparticular,increasedforestcoverreducedthe riskofmoderatefloodsandsoilerosion.Thecomparisonofpeakdischargeswithforest landversusnon-forestland[36]showedthatforestdecreasedpeakdischargeby ~1m3 s 1 formoderatestormintensity.Whenstormintensityincreased,thisbufferingcapacity decreased.Oursimulationshowedasimilarreductionindailydischargeswithincreased forestcover.However,thebufferingcapacityofforestinthisstudyseemsmaintainedfor largedischargeevents.

Theimpactofalteringforestlocalityalsovariedamongsubbasins.AtRioGrandede LoizaatCaguas(50055000)andRioGrandedeLoizabelowthedamsite(50059050),placing forestatlowaltitudesignificantlyreduceddischarge,especiallythelargedischarge,in comparisonwithplacingtheforestathighelevation.However,theeffectofshiftingforest localityforthesubbasinofRioGuraboatGurabo(50057000)isnotassignificantasthatfor theothertwosubbasins.

subbasin of 50055000 has higher elevation and steeper slope than 50057000 (Figure 8). In fact, the mean elevation and slope of 50055000 are 265.6 ± 247.8 m and 14.6 ± 14.4°, respectively, in comparison with 172.4 ± 125.2 m and 11.4 ± 9.5° for 50057000. The subbasin of 50059050 has mean elevation and slope as 208.2 ± 143.7 m and 12.9 ± 8.8°, respectively. Greater slopes may imply that the forest roots at higher elevation have less chance to offset subsurface flows before channelization, leaving more chances for forest at lower elevation to intercept runoff and subsurface flow. This explains why the subbasin of 50055000 has the strongest, whereas the subbasin of 50057000 has the weakest, and that of 50059050 has moderate capacity to reduce large discharges by relocating forest to low elevations.

The spatial arrangement, location, and pattern of vegetation or green infrastructure, especially forests and woodlands, have been found to be important for flood risk management. It is important to have green infrastructure located in the path of surface and baseflow to maximize the effect on flood attenuation [9,37,38]. Cross-slope and scattered woody patches are found more effective to smooth the river flow. For example, riparian forest is more functional than upland forest as it has more chances to uptake water and to hinder the surface and subsurface flows [37].

Theabilitytointerceptsurfaceandsubsurfaceflowsdependsalsoonthegeophysical factorssuchasrelativeelevationandslope,aswellaschannelizationofrunoff[9].The subbasinof50055000hashigherelevationandsteeperslopethan50057000(Figure 8). Infact,themeanelevationandslopeof50055000are265.6 ± 247.8mand14.6 ± 14.4◦ , respectively,incomparisonwith172.4 ± 125.2mand11.4 ± 9.5◦ for50057000.Thesubbasin of50059050hasmeanelevationandslopeas208.2 ± 143.7mand12.9 ± 8.8◦,respectively. Greaterslopesmayimplythattheforestrootsathigherelevationhavelesschancetooffset subsurfaceflowsbeforechannelization,leavingmorechancesforforestatlowerelevation tointerceptrunoffandsubsurfaceflow.Thisexplainswhythesubbasinof50055000has thestrongest,whereasthesubbasinof50057000hastheweakest,andthatof50059050has moderatecapacitytoreducelargedischargesbyrelocatingforesttolowelevations.

Figure 8. The elevation and slope distributions within the two upper-reach subbasins.

Figure8. Theelevationandslopedistributionswithinthetwoupper-reachsubbasins.

Land use also aﬀects the response of watersheds to changes in rainfall patterns. More forest cover tends to buﬀer the changes in rainfall variability. However, land use composition and pattern seem to be secondary compared with the geophysical factors: the subbasin of 50055000 had the greatest forest cover, yet its discharge variability has the steepest response to the changes in rainfall pattern (Figure 7). On the other hand, the subbasin of 50057000 has the least forest cover among the three subbasins, and the slope of response to changes in rainfall is not the highest.

Thespatialarrangement,location,andpatternofvegetationorgreeninfrastructure, especiallyforestsandwoodlands,havebeenfoundtobeimportantforfloodriskmanagement.Itisimportanttohavegreeninfrastructurelocatedinthepathofsurfaceand baseflowtomaximizetheeffectonfloodattenuation[9,37,38].Cross-slopeandscattered woodypatchesarefoundmoreeffectivetosmooththeriverflow.Forexample,riparian forestismorefunctionalthanuplandforestasithasmorechancestouptakewaterandto hinderthesurfaceandsubsurfaceflows[37].

Although more forests tend to offset large discharges, steeper slopes tend to generate more runoff as the stronger tangential gravity force tends to pull water down slope, leaving less water infiltrating soil layers [39]. This will provide less chance for baseflow to increase the variability of discharge [40]. Under the scenarios with reduced continuous raining days and more concentrated rainfall, watersheds with steeper slope and greater variation of elevation tend to amplify the effects of increased rainfall variability. In this study, the subbasin of Rio Grande de Loiza at Caguas (50055000) has the greatest and the subbasin of Rio Gurabo at Gurabo (50057000) has the flattest mean slope. Thus, the former has stronger response to decreased rainy days than the latter, despite the former having a larger forest cover. Our results thus imply that the geophysical factors dominate over

Landusealsoaffectstheresponseofwatershedstochangesinrainfallpatterns.More forestcovertendstobufferthechangesinrainfallvariability.However,landusecompositionandpatternseemtobesecondarycomparedwiththegeophysicalfactors:the subbasinof50055000hadthegreatestforestcover,yetitsdischargevariabilityhasthe

steepestresponsetothechangesinrainfallpattern(Figure 7).Ontheotherhand,the subbasinof50057000hastheleastforestcoveramongthethreesubbasins,andtheslopeof responsetochangesinrainfallisnotthehighest.

Althoughmoreforeststendtooffsetlargedischarges,steeperslopestendtogenerate morerunoffasthestrongertangentialgravityforcetendstopullwaterdownslope,leaving lesswaterinfiltratingsoillayers[39].Thiswillprovidelesschanceforbaseflowtoincrease thevariabilityofdischarge[40].Underthescenarioswithreducedcontinuousraining daysandmoreconcentratedrainfall,watershedswithsteeperslopeandgreatervariation ofelevationtendtoamplifytheeffectsofincreasedrainfallvariability.Inthisstudy,the subbasinofRioGrandedeLoizaatCaguas(50055000)hasthegreatestandthesubbasin ofRioGuraboatGurabo(50057000)hastheflattestmeanslope.Thus,theformerhas strongerresponsetodecreasedrainydaysthanthelatter,despitetheformerhavingalarger forestcover.Ourresultsthusimplythatthegeophysicalfactorsdominateovertheland usecompositioninwatershedresponsetothealteredrainfallvariability.Inadditionto theelevationandslopedifferences,biggerwatershedstendtocollectwaterwithdifferent timingsofprecipitation,whichtendtoyieldlessdischargevariability.Thismightapplyto thesubbasinof50059050,whichhastheweakestresponsetothescenarioswithreduced continuousrainingdays.

5.Conclusions

TropicalwatershedsintheCaribbeanregionareunderprogressivelyincreasingseasonalandinterannualclimatevariability,yetthelandusewithinthesewatershedsfeaturesgreatdeforestation/reforestationandurbanexpansion[41].Itisunclearhowthe changesinlandusecompositionandpatternwillaffectthewatershedresponsetoincreased climatevariability.

Inthispaper,weappliedSoilWaterAssessmentTools(SWAT+),acomputermodel forwatershedhydrologygivenlanduse,soil,anddigitalelevationmodel,tosimulate thestreamdischargesofthebiggestwatershedinPuertoRico.Weaimedtoexplorehow reforestationandforestlocalitywouldaffectthelargestreamdischargesandtheresponses ofthewatershedtoalteredrainfallvariability.

Wefoundthatreforestation(12–14%)dominatesoverurbanexpansion(5~8%)to reducelargestreamdischargesandtoincreasebaseflowinthreemajorsubbasins.Shifting foresttolowelevationtendedtoreducelargestreamdischarges.Furthermore,theextent ofchangingforestlocalitydependsstronglyonthegeophysicalfactorssuchasDEMand slope.Astrongerreductioninlargedischargeswithforestshiftedtolowelevationis associatedwithgreaterslopebecausethelow-elevationforesthasmorechancetointercept overlandandsubsurfaceflowsfromsteeperupperslopes.Forsubbasinswithrelativelyflat slopes,shiftingforesthaslimitedeffectsonlargedischarges.Wefoundincreasedrainfall anddischargevariabilitiesareassociatedwithreducedcontinuousraindays.Thechanges indischargevariabilitywithreducedcontinuousraindaysalsodependonDEMandslope sothatthesteepertheslope,thestrongertheincreaseindischargevariability.Henceto copewiththefutureincreasedclimatevariability,reforestationshouldmostlybeplacedat lowelevationrangeforbasinswithsteepslopes.

AuthorContributions: Conceptualization,Q.G.;methodology,Q.G.;formalanalysis,Q.G.andM.Y.; writing—originaldraftpreparation,Q.G.;writing—reviewandediting,Q.G.andM.Y.Allauthors havereadandagreedtothepublishedversionofthemanuscript.

Funding: ThisresearchwasfundedbytheNOAAPuertoRicoSeaGrantNA18OAR4170089.

DataAvailabilityStatement: ThedatasetsusedarepubliclyavailableattheNOAAdataportalat https://www.ncei.noaa.gov/cdo-web/,accessedon1January2024andat https://coast.noaa.gov/ dataviewer/#/,accessedon1January2024andtheUSGSdataportalat https://waterwatch.usgs.gov/, accessedon1January2024.

ConflictsofInterest: Theauthorsdeclarenoconflictsofinterest.Thefundershadnoroleinthedesign ofthestudy;inthecollection,analyses,orinterpretationofdata;inthewritingofthemanuscript;or inthedecisiontopublishtheresults.

References

1. Hernández,A.;Martin-Puertas,C.;Moffa-Sánchez,P.;Moreno-Chamarro,E.;Ortega,P.;Blockley,S.;Cobb,K.M.; Comas-Bru,L.; Giralt,S.;Goosse,H.;etal.Modesofclimatevariability:Synthesisandreviewofproxy-basedreconstructionsthroughthe Holocene. Earth-Sci.Rev. 2020, 209,103286.[CrossRef]

2. Ghil,M.;Lucarini,V.Thephysicsofclimatevariabilityandclimatechange. Rev.Mod.Phys. 2020, 92,035002.[CrossRef]

3. Karl,T.R.;Trenberth,K.E.ModernGlobalClimateChange. Science 2003, 302,1719–1723.[CrossRef][PubMed]

4. Kossin,J.P.;Knapp,K.R.;Olander,T.L.;Velden,C.S.Globalincreaseinmajortropicalcycloneexceedanceprobabilityoverthe pastfourdecades. Proc.Natl.Acad.Sci.USA 2020, 117,11975–11980.[CrossRef]

5. Hansford,M.R.;Plink-Björklund,P.;Jones,E.R.Globalquantitativeanalysesofriverdischargevariabilityandhydrographshape withrespecttoclimatetypes. Earth-Sci.Rev. 2020, 200,102977.[CrossRef]

6. D’Odorico,P.;Laio,F.;Porporato,A.;Ridolfi,L.;Rinaldo,A.;Rodriguez-Iturbe,I.EcohydrologyofTerrestrialEcosystems. BioScience 2010, 60,898–907.[CrossRef]

7. Andrés-Doménech,I.;García-Bartual,R.;Montanari,A.;Marco,J.Climateandhydrologicalvariability:Thecatchmentfiltering role. Hydrol.EarthSyst.Sci. 2015, 19,379–387.[CrossRef]

8. Dwarakish,G.S.;Ganasri,B.P.Impactoflandusechangeonhydrologicalsystems:Areviewofcurrentmodelingapproaches. CogentGeosci. 2015, 1,1115691.[CrossRef]

9. Antolini,F.;Tate,E.LocationMatters:AFrameworktoInvestigatetheSpatialCharacteristicsofDistributedFloodAttenuation. Water 2021, 13,2706.[CrossRef]

10. Yu,M.;Gao,Q.Topography,drainagecapability,andlegacyofdroughtdifferentiatetropicalecosystemresponsetoandrecovery frommajorhurricanes. Environ.Res.Lett. 2020, 15,104046.[CrossRef]

11. Pasch,R.J.;Penny,A.B.;Berg,R. NationalHurricaneCenterTropicalCycloneReport—HurricaneMaria(AL152017)16–30September 2017;NationalHurricaneCenter:Miami,FL,USA,2019;pp.1–48.

12. Pasch,R.J.;Reinhart,B.J.;Alaka,L. NationalHurricaneCenterTropicalCycloneReportHurricaneFiona(AL072022)14–23September 2022;NationalHurricaneCenter:Miami,FL,USA,2023;p.60.

13. Gutiérrez-Fonseca,P.E.;Ramírez,A.;Pringle,C.M.;Torres,P.J.;McDowell,W.H.;Covich,A.;Crowl,T.;Pérez-Reyes,O.Whenthe rainforestdries:DroughteffectsonamontanetropicalstreamecosysteminPuertoRico. Freshw.Sci. 2020, 39,197–212.[CrossRef]

14. Méndez-Lázaro,P.;Nieves-Santiango,A.;Miranda-Bermúdez,J.Trendsintotalrainfall,heavyrainevents,andnumberofdry daysinSanJuan,PuertoRico,1955–2009. Ecol.Soc. 2014, 19,50.[CrossRef]

15. Malmgren,B.A.;Winter,A.;Chen,D.ElNiño–SouthernOscillationandNorthAtlanticOscillationControlofClimateinPuerto Rico. J.Clim. 1998, 11,2713–2717.[CrossRef]

16. Torres-Valcárcel,A.R.TeleconnectionsbetweenENSOandrainfallanddroughtinPuertoRico. Int.J.Climatol. 2018, 38, e1190–e1204.[CrossRef]

17. Yu,M.;Gao,Q.Prolongedcoastalinundationdetectedwithsyntheticapertureradarsignificantlyretardedfunctionalrecoveryof mangrovesaftermajorhurricanes. Landsc.Ecol. 2023, 38,169–183.[CrossRef]

18. Grau,H.R.;Aide,M.;Zimmerman,J.K.;Thomlinson,J.R.;Helmer,E.;Zou,X.M.Theecologicalconsequencesofsocioeconomic andland-useChangesinpostagriculturePuertoRico. BioScience 2003, 53,1159–1168.[CrossRef]

19. OfficeforCoastalManagement. C-CAPLandCover,PuertoRico,2010;OfficeforCoastalManagement:NorthCharleston,SC,USA,2020.

20. Kennaway,T.;Helmer,E.H.TheforesttypesandagesclearedforlanddevelopmentinPuertoRico. GiscienceRemoteSens. 2007, 44,356–382.[CrossRef]

21. Gao,Q.;Yu,M.Reforestation-inducedchangesoflandscapecompositionandconfigurationmodulatefreshwatersupplyand floodingriskoftropicalwatersheds. PLoSONE 2017, 12,e0181315.[CrossRef]

22. Angulo,J.;Ashley,J.;Lebedev,A.;McNeff,N. AssessingtheLifespansofReservoirsinRegion2ofPuertoRico;WorcesterPolytechnic Institute:Worcester,MA,USA,2011.

23. Arnold,J.G.;Srinivasan,R.;Muttiah,R.S.;Williams,J.R.LargeareahydrologicmodelingandassessmentpartI:Modeldevelopment1. JAWRAJ.Am.WaterResour.Assoc. 1998, 34,73–89.[CrossRef]

24. Gassman,P.;Reyes,M.;Green,C.;Arnold,J.TheSoilandWaterAssessmentTool:HistoricalDevelopment,Applications,and FutureResearchDirections. Trans.ASABE 2007, 50,1211–1250.[CrossRef]

25. Helmer,E.;Ramos,O.;Lopez,T.D.M.;Quinones,M.;Diaz,W.MappingtheforesttypeandlandcoverofPuertoRico,acomponent oftheCaribbeanbiodiversityhotspot. Caribb.J.Sci. 2002, 38,165–183.

26. RCoreTerm. R:ALanguageandEnvironmentaforStatisticalComputing;RFoundationforStatisticalComputing:Vienna,Austria,2015.

27. Chapin,F.S.I.;Matson,P.A.;Vitousek,P.M. PrinciplesofTerrestrialEcosystemEcology,2nded.;Springer:NewYork,NY,USA,2011.

28. Francis,J.R.;Wuddivira,M.N.;Farrick,K.K.Exotictropicalpineforestimpactsonrainfallinterception:Canopy,understory,and litter. J.Hydrol. 2022, 609,127765.[CrossRef]

29. Tague,C.L.;Moritz,M.A.PlantAccessibleWaterStorageCapacityandTree-ScaleRootInteractionsDetermineHowForest DensityReductionsAlterForestWaterUseandProductivity. Front.For.Glob.Chang. 2019, 2,36.[CrossRef]

30. Marin-Spiotta,E.;Ostertag,R.;Silver,W.L.Long-termpatternsintropicalreforestation:Plantcommunitycompositionand abovegroundbiomassaccumulation. Ecol.Appl. 2007, 17,828–839.[CrossRef]

31. Yuan,F.;Lopez,J.J.;Arnold,S.;Brand,A.;Klein,J.;Schmidt,M.;Moseman,E.;Michels-Boyce,M.ForestationinPuertoRico,1970 topresent. J.Geogr.Geol. 2017, 9,30–41.[CrossRef]

32. Rodriguez-Martinez,J.;Santiago,M. TheEffectsofForestCoveronBaseFlowofStreamsintheMountainousInteriorofPuertoRico,2010; PuertoRicoDepartmentofNaturalandEnvironmentalResources:Virginia,WV,USA,2017.[CrossRef]

33. Khomsiati,N.;Suryoputro,N.;Yulistyorini,A.;Idfi,G.;Alias,N.Theeffectofforestareachangeintropicalislandstowards baseflowandstreamflow. IOPConf.Ser.EarthEnviron.Sci. 2021, 847,012032.[CrossRef]

34. Ramahaimandimby,Z.;Randriamaherisoa,A.;Vanclooster,M.;Bielders,C.L.DrivingFactorsoftheHydrologicalResponseofa TropicalWatershed:TheAnkaviaRiverBasininMadagascar. Water 2023, 15,2237.[CrossRef]

35. Bonnesoeur,V.;Locatelli,B.;Guariguata,M.R.;Ochoa-Tocachi,B.F.;Vanacker,V.;Mao,Z.;Stokes,A.;Mathez-Stiefel,S.-L.Impacts offorestsandforestationonhydrologicalservicesintheAndes:Asystematicreview. For.Ecol.Manag. 2019, 433,569–584. [CrossRef]

36. Bathurst,J.C.;Amezaga,J.;Cisneros,F.;GaviñoNovillo,M.;Iroumé,A.;Lenzi,M.A.;MinteguiAguirre,J.;Miranda,M.;Urciuolo, A.ForestsandfloodsinLatinAmerica:Science,management,policyandtheEPICFORCEproject. WaterInt. 2010, 35,114–131. [CrossRef]

37. Cooper,M.M.D.;Patil,S.D.;Nisbet,T.R.;Thomas,H.;Smith,A.R.;McDonald,M.A.Roleofforestedlandfornaturalflood managementintheUK:Areview. WIREsWater 2021, 8,e1541.[CrossRef]

38. Tamura,T.ImprovementoftheFlood-ReductionFunctionofForestsBasedonTheirInterceptionEvaporationandSurfaceStorage Capacities.In GreenInfrastructureandClimateChangeAdaptation:Function,ImplementationandGovernance;Nakamura,F.,Ed.; SpringerNatureSingapore:Singapore,2022;pp.93–104.

39. Brooks,K.N.;Ffolliott,P.F.;Magner,J.A. HydrologyandtheManagementofWatersheds,4thed.;Wiley-Blackwell:Danvers,MA,USA, 2013;p.545.

40. Schwyter,A.R.;Vaughan,K.L. IntroductiontoSoilScienceLaboratoryManual;UniversityofWyomingLibraries:Laramie,WY,USA,2020.

41. Yu,M.;Gao,Q.;Gao,C.;Wang,C.ExtentofNightWarmingandSpatiallyHeterogeneousCloudinessDifferentiateTemporal TrendofGreennessinMountainousTropicsintheNewCentury. Sci.Rep. 2017, 7,41256.[CrossRef][PubMed]

Disclaimer/Publisher’sNote: Thestatements,opinionsanddatacontainedinallpublicationsaresolelythoseoftheindividual author(s)andcontributor(s)andnotofMDPIand/ortheeditor(s).MDPIand/ortheeditor(s)disclaimresponsibilityforanyinjuryto peopleorpropertyresultingfromanyideas,methods,instructionsorproductsreferredtointhecontent.