severity along the Usumacinta River, Mexico : Identifying the anthropogenic signature of tropical forest conversion

Anthropogenic activities are altering flood frequency-magnitude distributions along many of the world ’ s large rivers. Yet isolating the impact of any single factor amongst the multitudes of competing anthropogenic drivers is a persistent challenge. The Usumacinta River in southeastern Mexico provides an opportunity to study the anthropogenic driver of tropical forest conversion in isolation, as the long meteorological and discharge records capture the river ’ s response to large-scale agricultural expansion without interference from development activities such as dams or channel modifications. We analyse continuous daily time series of precipitation, temperature, and discharge to identify long-term trends, and employ a novel approach to disentangle the signal of deforestation by normalising daily discharges by 90-day mean precipitation volumes from the contributing area in order to account for climatic variability. We also identify an anthropogenic signature of tropical forest conversion at the intra-annual scale, reproduce this signal using a distributed hydrological model (VMOD), and demonstrate that the continued conversion of tropical forest to agricultural land use will further exacerbate large-scale flooding. We find statistically significant increasing trends in annual minimum, mean, and maximum discharges that are not evident in either precipitation or temperature records, with mean monthly discharges increasing between 7% and 75% in the past decades. Model results demonstrate that forest cover loss is responsible for raising the 10-year return peak discharge by 25%, while the total conversion of forest to agricultural use would result in an additional 18% rise. These findings highlight the need for an integrated basin-wide approach to land management that considers the impacts of agricultural expansion on increased flood prevalence, and the economic and social costs involved.


Introduction
Global climate change and anthropogenic activities are disrupting flood frequency-magnitude distributions along many of the world's large rivers, posing threats to human populations and critical infrastructure.Many competing drivers contribute to modify a river's longterm discharge patternland-cover and land-use change through deforestation and agricultural expansion, urban expansion, water abstraction for human consumption and irrigation, channel modifications and infrastructural developments, as well as shifting precipitation patterns due to climate change (Bruijnzeel, 2004;Júnior et al., 2015;Malmer, 1992).Each of these drivers contributes its own signal to a river's hydrograph by altering the timing, duration, and magnitude of peak flows, yet evaluating the relative impact of any one of these drivers in isolation is a persistent challenge due to their interconnectedness, and constantly varying climatic conditions (Rogger et al., 2016;van Dijk et al., 2009).
The signal derived from forest conversion is particularly difficult to identify and quantify, as large-scale land-cover change historically takes place over long periods, amid many competing and overlapping developments, and often during times with spatially sparse climate and discharge data (Bruijnzeel, 1990(Bruijnzeel, , 1993)).Forest removal changes the pathways of precipitation through altering the processes of interception (the proportion reaching the surface), infiltration (the partition of surface and subsurface), and retention (the proportion reaching the rivernetwork after infiltration).Replacing dense tropical forest canopy with short vegetation (grass, shrubs, or crops) significantly reduces interception losses (Spracklen et al., 2018), compacts soil and removes organic material reducing infiltration rates (Germer et al., 2010;Muñoz-Villers and McDonnell, 2013), and dramatically decreases evapotranspiration rates, particularly in the humid tropics (Spracklen et al., 2018;Zhang et al., 2001).By reducing interception losses, evapotranspiration, and the infiltration rate of soils, the removal of forest increases the proportion of rainfall entering the river network, and the rate at which it reaches the network, thereby increasing the potential for large-scale flooding (Brown et al., 2005;Bruijnzeel, 1989;Fritsch, 1993).Whilst this has long been understood on the conceptual level (Blackie and Edwards, 1979;Clark, 1987;Gilmour et al., 1987), and many studies that relate tropical forest conversion to increases in mean river discharge exist at the small experimental scale (<1 km 2 ) and lower-mesoscale (<10 km 2 ), there are relatively few studies that examine the effects at the meso-or large-scale (>100 km 2 ) (Bruijnzeel, 1990(Bruijnzeel, , 1997;;Costa et al., 2003;Fritsch, 1993).In addition, due to poor data coverage, the relative inaccessibility of the study areas, and the temporal and spatial heterogeneity of competing factors, a number of early large-scale studies reported inconclusive or contradictory findings (Dyhr-Nielsen, 1986;Gentry and Lopez-Parodi, 1980;Qian, 1983;Richey et al., 1989;Wilk and Plermkamon, 2001).
Despite the complexities and complications that large-scale studies of this type have to contend with, there is a growing body of research that examines the hydrological impact of tropical forest conversion at the river-basin scale (Brown et al., 2005;Dias et al., 2015;Gao et al., 2020;Gerold and Leemhuis, 2008;Júnior, Tomasella, and Rodriguez, 2015;Karamage et al., 2017;Molina et al., 2012;Recha et al., 2012;Robinet et al., 2018;Sahin and Hall, 1996;Van der Weert, 1994).The overwhelming majority of these studies conclude that large-scale deforestation increases annual water yields and mean discharges, yet the large-scale impact on storm-flow generation, storm-flow-pathways, and how these pertain to flood generation is less clear (Bruijnzeel, 2004;Robinet et al., 2018).At the global scale, Bradshaw et al. (2007) undertook an analysis of flood severity and found that deforestation increases flood risk in tropical regions, yet a subsequent analysis suggested that 83% of the variation may be accounted for by population density as a root cause rather than a direct causal link with removal of trees (van Dijk et al., 2009).It is this persistent challenge of disentangling landscape responses to tropical forest conversion from competing anthropogenic and climate signals that has prohibited the systematic quantification of the impact forest conversion has on longterm flood risk (Bruijnzeel, 2004;van Dijk et al., 2009).
Here we attempt to fill this gap by clearly demonstrating the impact of forest conversion on peak discharge generation in the absence of competing anthropogenic influences, and accounting for climatic variation between our reference time-periods to significantly clarify the role tropical forest plays in large-scale flooding.This is an important contribution, as lacking a clear understanding of the role of forest cover in preventing future flood disasters has complicated the planning and implementation of effective land management strategies together with feasible flood mitigation policies (Calder and Aylward, 2006;FAO and CIFOR, 2005).

Data and methods
Here we analyse 55 years  of climate and discharge data from the Usumacinta River sub-basin in southeastern Mexico, which provides an opportunity to study the hydrological impact of large-scale tropical forest conversion in isolation.The past 20 years have seen an increase in the severity of flooding in the States of Tabasco and Chiapas in southeastern Mexico (Atreya et al., 2017;Gama et al., 2010Gama et al., , 2011)), which encompass the connected river basin of the Grijalva and Usumacinta Rivers.Unlike the Grijalva sub-basin, which has several large hydropower dams along the main channel, the flow of the Usumacinta sub-basin is unobstructed for the entirety of its course from the Guatemalan highlands to the Gulf of Mexico.The absence of large urban settlements and the relative inaccessibility of the terrain stifled development and large-scale forest conversion until the 1990 ′ s, meaning the 55 year discharge record from 1959 to 2014 captures the hydrological impact of deforestation on an otherwise consistent landscape.To identify and quantify an anthropogenic signal of tropical forest conversion in the discharge record of the Usumacinta River, we analysed continuous daily time series of precipitation, temperature, and discharge to compare long-term trends, and accounted for climatic variability between reference periods by examining the ratio of daily discharge to 90-day precipitation totals from across the contributing area.We then successfully reproduced this signal using a distributed hydrological model to simulate the response of the study area to widespread forest conversion (VMod: Lauri et al., 2006), and evaluated both the historical and future impact agricultural expansion has on the flood frequency-magnitude distribution along the Usumacinta River.In this paper, we employ a novel approach to disentangle climatic variability from within the daily discharge record at the intra-annual scale, and contribute to the wider, global scale debate concerning the role of forests in controlling large flood events, especially in the humid tropics.

Study area
The Grijalva and Usumacinta are the largest rivers in Mexico (Day et al., 2003); both begin in the highlands of Guatemala and traverse the States of Chiapas and Tabasco before converging in the low-lying floodplains just 50 km from the Gulf of Mexico.The combined Grijalva-Usumacinta river-basin covers some 130,000 km 2 and produces mean monthly discharges between 3000 and 6000 m 3 /s, equivalent to around 30% of the total surface runoff of the country (Areu-Rangel et al., 2019).The Grijalva and Usumacinta sub-basins are separated both hydrographically and socio-economically, as the Grijalva sub-basin was the historical focus of development within the region, leaving the Usumacinta sub-basin relatively undisturbed until the 1980s (Villela and Martínez, 2018).As such, the main population centres, areas of industrial oil and natural-gas extraction, projects of irrigated agriculture, and four hydropower dams, are all concentrated within the Grijalva subbasin.However, an intensification of agricultural development projects since the 1990s has seen rapid large-scale conversion of natural forest cover across much of the Usumacinta sub-basin, driven by agricultural expansion that extends into many of the headwater regions.Our study focuses on the Usumacinta sub-basin, measured at the Boca del Cerro gauging station located on the main branch of the Usumacinta River outside Tenosique (Fig. 1), which is the main urban centre along the river course, home to around 50,000 inhabitants (INE, 2005).We identify four periods within the discharge and climate records of the Usumacinta sub-basin that correspond to distinct stages of development, and draw comparisons between these land-cover signatures.

Discharge data
Our analysis of the Usumacinta flow regime is based on discharge data collected by the Servicio Meteorológico Nacional (SMN), under the Comisión Nacional del Agua (CONAGUA) of Mexico.We use data from the Boca del Cerro gauging station (30019) (Fig. 1), which is the only station within the Usumacinta sub-basin to provide a continuous record of average daily discharges (m 3 /s) across the entire study period .

Climate data
Observations of daily precipitation, maximum temperature, and minimum temperature were available at 15 weather stations that collectively form a continuous spatially robust dataset spanning the period 1959-1992.These data were provided by the Global Historical Climatology Network (GHCN) (Menne et al., 2012), andGarcía (1977).We infilled a spatial data gap located within Guatemala using additional data points taken from the Princeton University Global Meteorological Forcing (PGF) v.1 (Sheffield, Goteti, and Wood, 2006), formed of a suite of global observation-based datasets with the NCEP/NCAR reanalysis.We applied a bias correction to these data against ground data using a multi-variable scaling method (Santander Meteorology Group, 2015; Wilcke et al., 2013).After 1992, there are few GHCN stations within the Usumacinta sub-basin that record daily climate observations, therefore we supplemented these with a combination of the Tropical Rainfall Measuring Mission (TRMM) v7 satellite derived precipitation data, and the Climate Prediction Centre (CPC) Global Daily Temperature data for the period 1999-2018.To preserve model calibration across all timeperiods, we sampled the gridded data at points corresponding to the observation stations used in earlier iterations.We used daily composites of the rain gauge-adjusted, 3-hourly, 0.25-degree TRMM product (3B42) (Kummerow et al., 1998) that has been shown to reliably reproduce rainfall in the humid tropics, and has been used extensively in climate analyses and model forcing (Ferreira et al., 2012;Ji, 2006;Lauri et al., 2014;Shrivastava, 2014;Tapiador, 2017;Wu et al., 2015).A comparison of TRMM data with gauged observations taken within the Grijalva River sub-basin showed a strong correlation, particularly at monthly timescales (NSE: 0.6-0.82 with 30 day moving average).CPC daily temperatures combine the GHCN observation dataset and the Climate Anomaly Monitoring System (CAMS) dataset and interpolate them across a 0.5-degree grid (Fan and Dool, 2008), which has proven a reliable forcing for climate models (Nashwan, 2019).Comparing measurements of max and min temperature to ground observations across the wider Grijalva-Usumacinta river basin for the period 1998-2003, we found a consistent negative bias by the CPC dataset, which we again corrected using a scaling factor (Santander Meteorology Group, 2015; Wilcke et al., 2013).For our analysis of climate trends across timeperiods, we used spatially averaged data interpolated across the study area from point data used as our model inputs.

Model description and set up
The Environmental Impact Assessment Centre of Finland's (EIA) Integrated Water Resources Management modelling tool (IWRM VMod), is a physically based hydrological model distributed across a square grid representation that couples sub-models resolving energy and water balances at the grid scale, with a 1-D river-channel network model that routes outflows between cells.VMod first constructs a grid mesh over layered raster inputs representing elevation (m), soil type, flow direction, and land-cover class.It then interpolates daily climate data for each  (1959-1973), LC1992 (1978-1992), LC2004 (1999-2007), and LC2014 (2008-2014).
grid cell from input data at discrete points (max and min temperature ( • C), precipitation (mm)), and calculates potential evapotranspiration (PET) using the Hargreaves-Samani method (Hargreaves and Samani, 1982), before solving energy and mass balances across two subsurface soil layers and surface-atmosphere transfers following Dingman (1994).Runoff generated from each cell is then routed through a 1-D river channel network to give discharge outputs, which are calibrated against historical records.For a detailed description of the model construction and the governing equations, see Lauri et al. (2006).
For elevation data, we used SRTM 90 m (Jarvis et al., 2008), from which the model inferred flow direction data and the river channel network, which we adjusted to ensure alignment with satellite imagery.We prepared soil data from the FAO world soil map (FAO, 2009) by reclassifying the original classifications into six classes with default parameterisations, but later amended these as part of the calibration process.We then defined four periods each with a distinct land cover signature ranging from LC1970 with almost total forest cover, to LC2014 with just 42% dense forest cover across the study area.LC1970 land cover map was inferred from the International Satellite Land-Surface Climatology Project's (ISLSCP II, Ramankutty et al., 2010) historical land cover maps for 1950 and 1970 that showed little deforestation outside of the area around Tenosique, which corresponds with accounts of land use described in Tudela (1989).LC1992 land cover map corresponds to land use classifications described in the Central American Vegetation/Land Cover Classification andConservation Status (1992 -1993) (CCAD, 1998).LC2004 and LC2014 land cover maps were derived from MODIS land cover classifications (MCD12Q1) from 2004 and 2014 respectively (Friedl and Sulla-Menashe, 2019).Each of the land cover maps were reclassified from the original interpretations into five classes: water, forest, rain-fed cropland, pastureland, and urban.We then aggregated each of the raster inputs to match the model's grid sizing, which we set to a resolution of 2.5 × 2.5 km.

Model calibration, validation, and testing
As the focus of this study is to assess the impact of forest conversion on the hydrological regime of the Usumacinta River, we calibrated our model using 4 periods with distinct land cover signatures to distinguish between behaviours driven by soil type characteristics, and those driven by vegetation dynamics.Our initial calibration was against discharge data for the period 1978-1985, assuming a ubiquitous forest cover.Although this assumption contradicts our land cover map LC1992, this period has the most abundant and reliable climate data needed for a robust calibration of the soil parameters controlling the timing of runoff.We used the period of 1968-1973 as a validation of the initial calibration in a cyclical process to identify the interactions of soil and vegetation effects to best approximate the parameters for each soil type and the 'forest' land-cover class across both periods, which we then tested against the period 1959-1966.Having thus calibrated the soil type and forest land-cover parameters, we applied the model to the period 2008-2014, where variations in model performance stem entirely from forest land-cover alterations (assuming consistent soil parameters across time-periods).Whilst maintaining the model calibration from the initial periods , we introduced two additional land cover classes (rain-fed cropland and pastureland) to correct the model deficiencies using 2003-2007 as a validation period for a repetition of the cyclical calibration process.This final model calibration, including the nonforest land-cover classes, was tested against the periods 1986-1992 and 1999-2003.The initial calibration phase (1968)(1969)(1970)(1971)(1972)(1973)(1974)(1975)(1976)(1977)(1978)(1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)) required all major parameters to be adjusted, but primarily focused on hydraulic conductivities (horizontal and vertical directions), soil layer depths, storage capacities, weather interpolation values, surface runoff coefficients, and computational grid values.The second phase (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014) focused solely on defining the vegetation characteristics of the non-forest types, primarily the evapotranspiration and interception parameters, as well as the surface model components pertaining to vegetation differences.
For the calibration phases, we used the Nash-Sutcliffe efficiency coefficient (NSE; Nash and Sutcliffe, 1970) as the objective function.We then assessed the overall model performance against observed discharge by comparing relative biases of the total annual flow, low-flow, and high-flow indices (i.e., the 95th and 5th percentile of the discharge record respectively), as well as comparing the mean monthly discharges and distributions of annual maximum, minimum, and mean discharges across the entire discharge record .

Assessing hydrological and climatic changes between periods
To assess changes in the hydrological regime of the Usumacinta River, we first analysed trends in the long-term annual maxima, minima, and mean discharges for the duration of the discharge record.We repeated these analyses for mean annual temperature data, as well as total annual and seasonal precipitation data, where we defined the wettest season as June -November, and the drier season as December -May.To examine changes in the intra-annual flow regime, we divided the discharge and climate records into four periods for comparison, each representative of a distinct land-cover signature: LC1970 (used for period 1959-1973), LC1992 (1978-1992), LC2004 (1999-2007), and LC2014 (2008-2014).We compared the mean day-of-year (doy) discharges, temperatures, and precipitation (30-day totals) for each period against the long-term means across the period of the entire discharge record.Finally, to determine whether variances in mean discharges between LC periods are attributable to alterations in water availability, we normalised mean day of year (doy) discharges by average 90-day precipitation totals scaled by contributing area.In the case of LC1970, which has a ubiquitous covering of vegetation, the discharge recorded at the Boca del Cerro gauging station will be directly proportional to the amount of precipitation fallen within the study area, less interception/ evapotranspiration and changes to groundwater storage (and alterations to flow due to extraction or damsnot applicable in the Usumacinta).Interception/evapotranspiration is a function of temperature controlled by vegetation characteristics, and changes to groundwater storage can be assumed negligible when summed over multiple years.Therefore, normalising monthly discharge by precipitation totals scaled by area should reveal a consistent proportionality that reduces intra-annual variation, such that: where Nd i is the normalised discharge for the i th day of the year, n is the number of years in the observation record, d i,n is the discharge (m 3 /s) on the i th day of the n th year, A is the contributing area (m 2 ), and p j,n is the precipitation total (mm/day) on the j th day of the n th year.Nd i then represents the dimensionless (after unit conversion factors) proportion of discharge to the average amount of water fallen over the entire study area as precipitation in the previous 90-days.We find that 90-day precipitation totals are most suitable for removing intra-annual variation to consistent proportionalities in this study area.
Were the periods 1978-1992(LC1992), 1999-2007(LC2004), and 2008-2014 (LC2014) to maintain the same ubiquitous forest cover, then the discharge record at the Boca del Cerro gauging station should display the same proportionality to precipitation as displayed in LC1970 (assuming a consistent temperature distribution across time-periods).Any variation from the intra-annual normalised discharge pattern displayed in LC1970 will be the result of alterations to the vegetation effects controlling interception/evapotranspiration, and thus represents an anthropogenic signature of forest conversion to agricultural expansion.

Forest conversion scenarios
To assess the impact of potential future forest conversion on the hydrological regime of the Usumacinta River, we developed scenarios that progressively removed forest area according to observed historical patterns projected into the future.LC1970, LC1992, LC2004, and LC2014 represent 98%, 87.3%, 73.3%, and 42.1% forest cover respectively.We randomly converted forested pixels from the initial land cover map (LC1970) sampled from areas later converted to either crop agriculture or pastureland (as displayed in LC1992, LC2004, and LC2014), maintaining the observed proportions of each (as displayed in LC2004the most reliable partition of land classes), to infill the proportion of non-forest cover for 100-50% forest cover scenarios, after which the forest conversion was entirely randomised for 25% and 0% forest cover scenarios.We then used these progressive land cover maps to run forestconversion model scenarios across the most complete record of climate forcings (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018).This allowed us to draw direct comparison of the mean doy discharges and hydrological extremes under different projections of forest conversion, and to assess the likely impact of continued agricultural expansion on severe flooding along the Usumacinta River, and by extension, the Grijalva River.

Long term climate trends and hydrological signals
Using interpolations of daily climate data, and discharge data spanning the entire 55-year study period , we found statistically significant positive trends for the mean annual temperature (pvalue < 0.01, Fig. 2a), and annual maximum, mean, and 10th-percentile (low flow) discharges (p-values 0.1, 0.04, and <0.01 respectively, Fig. 2c).We found no statistically significant trends in neither the total annual precipitation (Fig. 2b), the total drier season precipitation, nor the total wet season precipitation across the years.
From a comparison of the mean daily discharges taken for each of the land-cover classification periods (LC1970-LC2014), there are considerable increases in the first wet season peak (Jun-Aug) and again in the second peak (Sep-Nov) (Fig. 3a).Whilst the drier-season flows look comparatively stable across the periods, the proportional increases from the historical base discharge show statistically significant drier season gains (positive 7-60%, Fig. 4a), which is not evident in the drier season precipitation between time-periods (Fig. 4c).The largest proportional change comes at the onset of the wet-season in June, where mean monthly discharges are >75% larger in LC2014 compared to LC1970 (Fig. 4a), whilst in terms of magnitude, discharge increases in October are equivalent to those in June, with LC2014 exhibiting a mean monthly discharge ~900 m 3 /s larger than that for LC1970.To characterise the different landscape responses under each of the land-cover classifications independently of the varying climatic conditions, we compared the ratio of daily discharge totals at the mouth of the study area to the mean daily precipitation to have fallen across the entire study area for the previous 90 days.The results show an average proportionality between 0 and 1 that is relatively consistent throughout the year (Fig. 3d) when compared to the variation displayed in the precipitation totals (Fig. 3c).However, this proportionality alters dramatically between early and later land-cover classifications, with LC2004 and LC2014 showing changes of up to 58% (May) from those of LC1970 (Fig. 4c), predominantly in the drier season months where the ratio of mean daily discharge to precipitation has risen from ~0.45 to consistently above 0.7 (Fig. 3c).

Hydrological model calibration and validation
Both the initial calibration phase that concentrated on soil characteristics and forest cover parameterisation , and the later calibration phase that focused on defining the non-forest land-cover types (2008-14) display good agreements with observed data, each obtaining an NSE of >0.80 (Table 1).The test phases, which were not included as part of the calibration procedure, display NSEs of 0.64 and 0.74 respectively.The slightly poorer fit to the older test period may in part be due to the reliability of the forcing climate data, which was sparser and contained a number of data gaps making the interpolation less consistent.Overall, the modelled discharge series displayed a robust agreement to the observed data set, with an NSE of 0.76 (Table 1).
In addition to consistently performing better than the sample mean (indicated by the NSE), an important component of hydrological modelling is the faithful reproduction of key aspects of the regime.Comparing the flow duration curve (Fig. 5b), the distribution of annual means (Fig. 5c middle lines), and the mean daily discharges (Fig. 5d), the model performs well -reproducing the distribution of flows characteristic of the Usumacinta River.However, as both the measures of bias (Table 1) and the distributions of maximum and minimum discharges (Fig. 5c top and bottom) attest, the model tends to marginally underestimate high-flows (Q5), and overestimate low-flows (Q95).
For our purposes, the most important component of the model is the representation of hydrological processes affected by forest cover, and   the impact of forest conversion to agricultural land use on the hydrological regime.As the differences in discharge records between landcover classification periods are driven by the combined effects of climatic variability and landscape dynamics, normalising the discharge record by variations in climate should reveal a signal of forest conversion to agriculture.This signal is present in the observed discharge record, where LC2004 and LC2014 display fundamentally different behaviour with respect to the proportion of water reaching the study area outflow (Fig. 4d).This same signal is present in the model results (Fig. 6a), which suggests that the model represents the impact of vegetation cover on the hydrological cycle adequately, and the effect of forest conversion on discharge.Re-running the model with a homogeneous forest land-cover class, whilst maintaining the original climate input data, produces results that display a more uniform proportionalitywhere the ratio of daily discharge to mean precipitation across the study area more closely resembles that of the historical base case LC1970 (Fig. 6b).

Hydrological analysis under future forest conversion scenarios
To investigate the impact of forest conversion to agricultural landuse on the hydrological regime and flood magnitude-frequency distribution of the Usumacinta river, we ran the calibrated model using climate data from 1999 to 2018 under different scenarios of forest cover, representing 100%, 50%, 25%, and 0% forest cover, as well as the 2018 land cover classification.We found that each successive forest cover scenario shows a clear increase in discharge throughout the year (Fig. 7a), though the increase is not uniformly proportional (Fig. 7b).Drier season flows (Dec-May) display a larger proportional increase with decreasing forest cover compared to the wet season (Jun-Nov).The 2014 land cover scenario has a forest cover of 42%, yet exhibits discharge patterns that more closely resemble FC25 (25% forest cover scenario) than FC50 in the drier season months.This is most likely due to the proportion of cropland compared to pastureland represented, as we maintained the cropland: pastureland ratio observed in LC2004 throughout forest cover scenarios, while the observations of LC2014  identify a smaller proportion of cropland.Lastly, we fitted a generalized extreme value (GEV) distribution function to the annual maxima for each of the forest cover scenarios (including the additional FC625 scenariowith 62.5% forest cover) to ascertain to what extent agricultural expansion has increased the expected return flood historically, and to what extent it is likely to increase it in the future.We found that the 10-year return flood under the current land cover classification (LC2014) has increased 25% compared to FC100, and that the continued expansion of agricultural land-use could increase it a further 10% under the FC25 projection, and 18% under the FC0 (Fig. 8).This means that the return period for a record high flood during the study period, i.e. approximately the 2008 peak discharge (Fig. 5a), would fall from the current estimate of 22 years under LC2014 to just 8 years under the total forest conversion scenario (FC0).

Discussion
The hydrological changes evident in our scenarios of forest conversion within the Usumacinta river sub-basin clearly suggest that forests play an important role in controlling the frequency and magnitude of floods in the humid tropics of southeastern Mexico.In addition, these results may offer broader insights into the hydrological functioning of forests in similar climate conditions, and support the assertion that forests play an important role in controlling floods globally.Here we discuss the implications and processes underlying these results, as well as additional drivers and externalitiesboth biophysical and linked socio-political issues.

Impact of forest conversion on hydrological processes
The conversion of forest to short vegetation is typically associated with a large reduction in interception losses over longer time-periods (Spracklen et al., 2018).However, interception losses vary with precipitation intensity -the highest losses are associated with lower intensity events, whereas interception losses are unlikely to significantly affect higher intensity precipitation events, as the canopy capacity is  rapidly exceeded (Bandeira et al., 2018;Fleischbein et al., 2005;van Dijk et al., 2009).Evapotranspiration is another mechanism by which forests can reduce the proportion of precipitation reaching the rivernetwork compared to grass and cropland.Forests have rates of evapotranspiration between 20% and 80% greater than tropical grasslands (Schlesinger and Jasechko, 2014;Spracklen et al., 2018;von Randow et al., 2004;Zhang et al., 2001).As with interception losses, the difference in evapotranspiration between forest and shorter vegetation will have a pronounced effect on base-flow characteristics and small-to midsized flood peaks resulting from prolonged lower-intensity precipitation events, than upon floods produced by short duration extreme rainfall (Brown et al., 2005;Bathurst et al., 2011).The third mechanism by which forest conversion affects the delivery of water to the rivernetwork is by altering the permeability of soils, and thus the partitioning of surface to subsurface flow.The exposure of bare soil to intense rainfall (Lal, 1987(Lal, , 1996)), the compaction of topsoil by machinery or grazing (Gilmour et al., 1987;Kamaruzaman, 1991), and the removal of roots and organic intrusions (Lal, 1983;Aina, 1984), all contribute to reduce soil permeability and rainfall infiltration after forest conversion (Bruijnzeel, 2004;Germer et al., 2010;Moraes et al., 2006;Muñoz-Villers and McDonnell, 2013).Unlike the previous two mechanisms, a reduction in soil permeability (and the associated increase in runoff generation) will significantly affect the magnitude and timing of flood peaks -particularly those resulting from the most intense precipitation events (Kamaruzaman, 1991;Robinet et al., 2018;Van der Plas and Bruijnzeel, 1993).
The combined effect of reduced interception losses and evapotranspiration following forest conversion should be most evident in the drier season, when the volume of water returned to the atmosphere is a significant proportion of the total volume that falls across the study area.Our results clearly support these hypotheses, as the ratio of discharge at the study area outflow to the average precipitation across the study area shows a dramatic increase in the drier season months for land cover classifications LC2004 and LC2014 (Fig. 4d).The difference between interception and evapotranspiration losses following forest conversion are less significant in the wetter months, as evidenced in the more consistent ratio of discharge to precipitation across land cover classes between June and December (Fig. 4d).Whilst the total volume of water reaching the river-network after forest clearing during the wetter season may not alter as significantly as during the drier season, the timing and magnitude of peak discharges may shift as a result of increased runoff generation due to reduced soil infiltration rates.This may account for the observed increases in discharge during Aug -October, despite there being no significant alteration in precipitation totals (Fig. 4).
An analysis of extreme value distributions requires sufficiently long time-series data to encompass a range of extremes; therefore, a direct comparison of historical discharges across each of the land-cover classification periods does not yield robust results, as the time-periods only cover 7-15 years.However, using our calibrated model, we simulated the expected range of annual extremes across a consistent 20-year period (1999-2018) of climate data for each scenario of forest cover (FC100-FC0).This allowed us to explore the impact of both historical and future forest conversion on the flood magnitude-frequency distribution of the Usumacinta, and to characterise the role that forests play in mediating large-scale flood events.Our results indicate that the largescale conversion of forests to agriculture has intensified discharge extremes, and that continued conversion is likely to exacerbate fluvial flooding in the future.Whilst the quantification of this impact may only be valid within the Usumacinta context, the general patterns observed here are indicative of alterations that forest conversion makes to the underlying processes that mediate large-scale flooding.Due to the interconnectedness of the Usumacinta and Grijalva sub-basins, and the resemblance of their topographic and climatic characteristics, we can reasonably assert that our results will hold true across the wider Grijalva-Usumacinta basin.Whilst the current flood regime of the Grijalva differs significantly from the Usumacinta, due to the large-scale infrastructural developments, including four hydropower dams along its course, our results nevertheless have important implications for the management of reservoirs if the conversion of forests to agriculture and pastureland in the Grijalva upper-basin continues.Reassessing the reservoirs' operating levels may be necessary to take into account the shift in landscape response to intense rainfall events, and to accommodate more frequent higher-magnitude discharge events during the wettest periods.

Additional externalities
In addition to forest conversion, the major biophysical factors that will affect the flood frequency-magnitude distribution in the Grijalva-Usumacinta river basin in the coming decades are climatic changes, major shifts in agricultural production, urban expansion, and infrastructural development.Of these, climate change has the most potential to radically alter the hydrological regime of southeastern Mexico.The expected impact of global warming on future patterns of precipitation across southeastern Mexico is a reduction of annual totals and protraction of drier periods, with the possibility of increased extreme events during the wetter season (Fuentes et al., 2015;Imbach et al., 2018;Karmalkar et al., 2011).Less frequent, more intense rainfall interspersed with longer dry periods will likely result in a more rapid conveyance of rainfall to the river network.Increased soil degradation and compaction during dry spells may further reduce the infiltration capacity of soils (Batey, 2009;Bruijnzeel, 2004), followed by an intensification of precipitation extremes that will generate a higher proportion of runoffpotentially overloading the river network capacity and causing widespread flooding.Each of the biophysical factors affecting the flood regime of the Grijalva-Usumacinta has the potential to either exacerbate or counteract the impact of future climate change, but this will depend upon the multifaceted interplay of sociopolitical drivers that shape the land-use patterns, rural development strategies, and plans for urbanization in the region.
Since the 1980s, the major agricultural activity that has superseded forest conversion in the Grijalva-Usumacinta river-basin is extensive cattle raising, which is often practiced not only for the economic value of the herd but also as an indication of land ownership and a form of land speculation (Kaimowitz and Angelsen, 2008;Tudela, 1989).At an average of one head of cattle per hectare, the profitability of extensive cattle raising in the humid tropics is considerably lower per unit of area than in many alternative production systems, including cacao, citrus, or intensive agrosilvopastoral systems (Nahde-Toral et al., 2013).The reason why many farmers prefer extensive cattle raising is that it requires relatively few investments and the labour costs are relatively low (Tudela, 1989).Several governmental agricultural development programmes in the region have also promoted cattle raising, with subsidies, low-interest loans, and technical assistance.
A major shift in agricultural and environmental policies that promotes the reforestation of large portions of the upper-basins might mitigate some of the negative impacts of forest removal on flood generation (Bruijnzeel, 2004).At present, there is an initiative to reforest 1 million hectares of pastures and croplands with fruit and timber trees across the Grijalva-Usumacinta riverbasin by the federal government; however, the tendency seems to be land conversion to largescale sugar cane, African oil palm, and gmelina plantations (El Heraldo de Tabasco, 3 of April 2020).If this one-million-hectare reforestation programme succeeded in its goals, this could shift the current land cover scenario (LC2014) towards the FC625 land cover scenario (see Fig. 8) in terms of high vegetation cover.Whilst reforestation's capacity to entirely rehabilitate the hydrological functioning of converted land is uncertain (Bruijnzeel, 1997(Bruijnzeel, , 2004)), restoring comparable levels of interception losses and evapotranspiration to pre-conversion levels could be expected within 3-5 years (Malmer, 1992;Brown et al., 1997), accompanied by a significant reduction in peak-and storm-flows (Bruijnzeel, 1989).A transition from the flow regime associated with LC2014 to FC625 would reduce the occurrence of the expected 20-year return peak discharge to once every 36 years (Fig. 8).Such a reduction in flood frequency would have enormous economic and social ramifications, as the 2007 flood alone has been estimated to have incurred losses to the state of Tabasco between 800 and 3000 million USD (Kauffer, 2013;Räsänen et al., 2017;Ishizawa et al., 2017).A thorough investigation of environmental, social and economic impacts of basin-wide reforestation programmes should be undertaken to assess the long-term viability of balancing the cost of reforestation programmes against the associated reductions in flood damage expenses.
As a result of a severe flood disaster in 1999, the federal government of Mexico initiated an Integrated Flood Control Program (PICI) in the State of Tabasco (Aparicio et al., 2009).However, due to a series of delays and budgetary ambiguities, many of these projects were abandoned before completion and few were operational when the 2007 flood event occurred (Perevochtchikova and Lezama de la Torre, 2010).A renewed effort following the devastating 2007 flood saw a number of infrastructural projects completed, which appear to have reduced the socioeconomic impacts of the 2010 flood, despite its magnitude being larger than that of 2007 (Ishizawa et al., 2017).However, an issue with many of these projects is that they divert the flood hazard from one place to another, rather than eliminating the risk of serious flooding.Through levees, gate structures, embankments, and water channels, floodwater has been redirected from the more affluent urban areas to socioeconomically vulnerable rural areas and indigenous communities, who are under-represented in the political decision-making, lack the resources to protect themselves from the negative consequences of flooding and to recover quickly after a disaster (Nygren, 2016).In recent years, the long-term sustainability of technocentric flood-control measures has been questioned, and there has been a shift towards more integrated flood-management approaches, with basin-wide land-use strategies to "make room for the river", and to enhance residents' social resilience to flooding (Butler and Pidgeon, 2011;Nygren, 2016;Sletto and Nygren, 2015;Rijke et al., 2012;Räsänen et al., 2017).To avoid future catastrophes in southeastern Mexico, a broader and more integrated approach to flood management is needed with strategies to reduce the strain on the river-network capacity, and shift from merely water resources management towards approaches that consider the complex interactions and overlappings between river basin management, land-use changes, hydrological infrastructure, coastal zone restoration, and flood-risk prevention.

Limitations of the study and opportunities for future research
We focused our study on the Usumacinta River because its channel has remained free from infrastructural development, and agricultural expansion happened at time when the landscape response to widespread forest conversion was captured in satellite imagery and gauging station records.This allowed us to characterise the anthropogenic signal within the discharge record generated from agricultural expansion that would otherwise have been lost amongst competing signals.However, the benefits derived from the study area's relative remoteness come with certain disadvantages and challenges.The main constraint in this research has been data availability and coverage, a limitation common to these kinds of studies, but particularly difficult in this instance due to the geographical remoteness, and the transboundary nature of the Usumacinta River.The upper reaches of the Usumacinta sub-basin are located within the boundary of Guatemala, where there are very few records of climate data and no records of discharge.Even in the relatively data-rich areas within Mexico, data sets lacked consistency and continuity.We acknowledge the limitations and uncertainty that these data constraints place on our work, nevertheless; supplementing observed data with satellite-derived data we were able to consistently replicate discharge records across the entire study period.
The methods used in this study to explore the processes that shape the landscape's hydrological response to widespread forest conversion would be directly transferable to river-basins with similar characteristics across the tropics.However, our assessment of future flood patterns is not directly transferable to the wider Grijalva-Usumacinta river-basin, as these will also depend upon the operation of the reservoirs along the Grijalva River, which issue was outside the scope of this study.To comprehensively account for future flood risk within this area, a detailed analysis of future trends in discharge that incorporates the current reservoir operating rules is required.Such a study should also incorporate scenarios of future changes to precipitation and temperature patterns due to climate change.To implement a fully integrated flood risk-management strategy, a more thorough understanding of the changing flood-management policies and land-uses patterns, and their linkages to residents' socio-economically differentiated vulnerabilities and capabilities of flood resilience, would also be needed.

Conclusion
This study is one of the few comprehensive assessments that quantifies the impact of widespread tropical forest conversion on river discharge and flood-magnitude conducted at the large scale.By analysing 55 years of climate and discharge data, we identified a signal of forest conversion within the discharge record of the Usumacinta River, southeastern Mexico.Comparing the proportion of water falling as precipitation to that reaching the study area outflow between stages of deforestation extent, we found tropical forest conversion significantly alters the hydrological functioning of the landscape.
Between 2010 and 2020, the net loss of global forests is estimated at 4.7 million ha per year, with the large majority of that loss occurring in the tropical regions of Latin America and Africa (FAO, 2020).In addition to the large carbon source that this loss represents, there are a number of ecosystem services essential to humanity's well-being that can no longer function, such as; weather regulation (Spracklen et al., 2018), air purification, biodiversity protection, and pest control (Chivian, 2002;Diaz et al., 2006;Laurance and Williamson, 2001).Yet, despite this, governments are reluctant to conserve tropical forests as the economic valuation of these ecosystem services is difficult to quantify, and the global demand for timber and agroindustrial products incentivises forest conversion for agriculture and the commercial extraction of forests to rapidly generate capital (Rudel et al., 2009;Hosonuma et al., 2012).Our findings bring to light the potential for tropical forests to play a key role in the mitigation of large flood events, and the impact continued deforestation can have on the magnitude and frequency of future flood events across the tropics.Due to the socio-economic costs and environmental impacts these increases in flood magnitude represent, such findings may contribute to a holistic evaluation of the benefits derived from conserving forest cover, and promote the implementation of integrated flood-management approaches that include comprehensive river-basin-wide land-use and resource-management practices.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 4 .
Fig. 4. Proportional change from the historic base case of LC1970 (1959-1973) for discharge (a), temperature (b), precipitation (c), and ratio of discharge to precipitation (d) at Boca del Cerro station of Usumacinta sub-basin.

Fig. 5 .
Fig. 5. A) Observed (blue) and modelled (red) discharges for entire study period at Boca del Cerro station of Usumacinta sub-basin.B) Flow duration curve showing distribution of observed and modelled discharges.C) Distribution of maximum (top), mean (middle), and minimum (bottom) discharges for observed and modelled discharges.D) Mean day of year discharges taken across the whole study period for observations (blue) and modelled results (red) shaded between the 5th and 95th percentiles, with max and min points marked.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Proportional change from the historic base case of 1959-1973 for the ratio of discharge to precipitation, at Boca del Cerro station of Usumacinta sub-basin, using observed land cover classifications (a), and assigning a uniform forest land cover classification for the same periods (b).

Fig. 7 .Fig. 8 .
Fig. 7. Impact of forest conversion at Boca del Cerro station of Usumacinta sub-basin.A) Mean day of year discharges for forest conversion scenarios of 100% (green), 50% (blue), 25% (yellow), and 0% (red), as well as the land cover classification for 2014 (dotted black).B) Proportional change of the mean monthly discharge from the 100% forest cover scenario.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1
Model verification indicators at Boca del Cerro station of Usumacinta sub-basin (orange: calibration, blue: validation, green: test), LC refers to periods of distinct land cover classifications.