Climate and Land Cover Trends Affecting Freshwater Inputs to a Fjord in Northwestern Patagonia

Introduction

The land-ocean interface in river deltas and fjords is the natural bridge between terrestrial and marine systems. In these areas, nutrients, organic matter and sediments transported by rivers meet with those coming from the coastal and open ocean. Globally, oceans receive more than 36,000 km³ of freshwater and about 20 billion tons of sediments (Milliman and Farnsworth, 2011). About 95% of this volume enters through the rivers. The freshwater input represents important regulating and provisioning ecosystem services, influencing the coastal system’s physical-chemical characteristics, supporting the reproduction of numerous species, many of which sustain important fisheries and subsistence activities (Nixon and Buckley, 2002; Barbier et al., 2011).

The coastal system of Northwestern Patagonia (41–45°S) is comprised by extensive and interconnected fjords, bays and channels that receive important freshwater inputs from the southern Andes Range (Milliman and Farnsworth, 2011; Pantoja et al., 2011; Iriarte et al., 2014). Most of the freshwater inputs come from numerous rivers that are fed by rainfall runoff, snow melt and glacier melt from relatively small basins (< 10,000 km²). Freshwater inputs influence the coastal dynamics at different spatial and temporal scales, producing marked density gradients in the seawater column (Dávila et al., 2002; Saldías et al., 2019), changing nutrient ratios, and bringing high concentrations of organic matter (dissolved and particulate) and silicic acid (Silva et al., 2011; Vargas et al., 2011; Torres et al., 2020). Inorganic nutrients and solar radiation are strongly seasonal, limiting the relatively high primary productivity of estuarine systems in western Patagonia (0.5–3.0 g C m^–2 d^–1) dominated mainly by diatoms (90%) (Iriarte et al., 2007; Jacob et al., 2014; Raymond and Spencer, 2015). Changes in freshwater input could have negative implications on oceanographic processes to fjords in the inner seas of Chilean Patagonia. The ecosystem services that these systems provide depend on the interplay of different water masses, especially including the upper estuarine layer (Sievers and Silva, 2008), which strongly depends on the magnitude, timing, and duration of freshwater input. The depth, extent and chemistry of this layer/water mass also depends on the complexity of processes occurring at the climate-dependent river-ocean interface.

The influence of rivers on the coastal systems in this region has weakened during the last decades (León-Muñoz et al., 2018; Aguayo et al., 2019, 2021) due to a marked decrease in precipitation during summer and autumn (Boisier et al., 2018). The precipitation decline has been attributed to the positive trend of the Southern Annular Mode (SAM), an index that describes the movement of the low-pressure belt that generates westerly winds. The SAM trend has been attributed to stratospheric ozone depletion and increased greenhouse gases concentration, which suggests the effects of a global anthropogenic forcing on the local precipitation regime (Arblaster and Meehl, 2006; Eyring et al., 2013; Boisier et al., 2018). This trend has increased the climate synergy with El Niño-Southern Oscillation (ENSO), promoting severe droughts during austral summers (Garreaud, 2018). Climate projections for the next decades (2020–2070) follow this trend. Recent studies have evaluated the hydrological impacts of climate projections in the region (Puelo River; 41°S) predicting a prolongation of warm and dry conditions, seasonal changes in hydrological regimes and an increase in severe drought frequency (Aguayo et al., 2019, 2021; Pessacg et al., 2020).

Apart from shifts in climate conditions, the magnitude and spatial expansion of land cover changes (LCC) may strongly alter river hydrology (Ellison et al., 2017) and biogeochemistry (Cuevas et al., 2006; Raymond et al., 2008; Aufdenkampe et al., 2011), and thus their influence on the coastal systems. A large proportion of native forest cover has been identified as a factor providing resilience to coastal systems strongly influenced by freshwater against climate change (Desmit et al., 2018; Khoury and Coomes, 2020). Native forest substitution by exotic tree plantations in south-central Chile (34–41°S; Miranda et al., 2017; Uribe et al., 2020) has been identified as one of the causes of lower streamflow during summer (Little et al., 2009; Iroumé and Palacios, 2013; Alvarez-Garreton et al., 2019). Climate change has also increased the recurrence of wildfires in this region (Urrutia-Jalabert et al., 2018), which could potentially alter hydrological processes (Rulli and Rosso, 2007; Boisramé et al., 2019). The impact of LCC on riverine biogeochemistry can affect nutrients and carbon cycling. Dissolved organic matter (DOM) composition is strongly influenced by land cover. For example, different types of vegetation reflected as divergent biomarkers in coastal waters (Hernes and Benner, 2003; Raymond and Spencer, 2015).

Anthropogenic alterations of the natural streamflow regime by way of LCCs make isolating the impacts of climate change on hydrology a daunting task (Yang et al., 2017; Chanapathi and Thatikonda, 2020). In Chile the hydrological impacts of LCCs have been evaluated mainly in south-central Chile (< 41°S; Stehr et al., 2010; Aguayo et al., 2016; Martínez-Retureta et al., 2020; Barría et al., 2021; Galleguillos et al., 2021) and not in tributary basins to the coastal system of northwestern Patagonia. Nevertheless, using rivers like Puelo River as sentinels, several studies have reported anomalous events in fjords of northwestern Patagonia (e.g., surface hypoxia events, harmful algal blooms), mainly during severe droughts (Valle-Levinson et al., 2007; León-Muñoz et al., 2013, 2018; Iriarte et al., 2017; Díaz et al., 2021). Although northern Patagonian basins have not yet been modified by irrigation channels or reservoirs and still maintain extensive little-degraded native forest cover (CONAF and UACH, 2014; Astorga et al., 2018), recent threats such as increased wildfires and their joint effect with climate change have not been evaluated hydrologically to date. Considering this context, here we evaluate the influence of climate and land cover trends on the freshwater inputs of the Puelo River Basin to the Reloncaví Fjord (41.5°S; Figure 1). Specifically, we aim to 1) characterize the water quality of freshwater inputs along the Puelo River network and the Reloncaví Fjord (2018–2019), 2) evaluate trends in natural flow regime (1950–2019), precipitation-streamflow associations (1950–2019), and land cover (2001–2016) for the Puelo River Basin, and 3) estimate the potential future impacts of LCC (2030), mainly caused by forest fires, on the hydrological processes of the Puelo River.

Materials and Methods

Study Area

The study area comprises the Puelo River Basin and the Reloncaví Fjord (Figure 1). The basin (∼9,000 km²; 66% of its drainage area in Argentina and 34% in Chile) is located in the central latitudes covered by the Valdivian Rainforest ecoregion. An important part of the basin is covered by native forest (CONAF and UACH, 2014), with minimal anthropogenic interventions (24% of the basin is located in protected areas). According to Köppen-Geiser’s classification, the basin has temperate and polar climates with tundra characteristics (Beck et al., 2018). The precipitation range follows a strong longitudinal gradient that can vary between 500 mm and 4700 mm per year (Lara et al., 2008; Aguayo et al., 2019). Half of annual precipitation is concentrated during the austral winter. Recently, Aguayo et al. (2021) determined that the Puelo River Basin has a low aridity index (< 1), suggesting that it is energy-limited. The same basin showed, on average, an evaporation index close to 0 during all seasons except in summer, where evapotranspiration represents 37% of precipitation.

The Puelo River has a mean streamflow of 640 m³ s^–1 (2400 mm yr^–1) with a pluvial-nival hydrological regime in the proximity of its mouth, with a maximum peak in winter and a second and lower peak during spring (León-Muñoz et al., 2013). In contrast, the sub-basins located at higher elevations have a nival-pluvial hydrological regime. Puelo River is classified as a high-runoff river (>750 mm yr^–1); it has one of the highest runoff levels in mesoscale basins (Milliman and Farnsworth, 2011). During austral summer and autumn (low-flow season), the streamflow of the Puelo River is significantly correlated with the streamflow of other rivers of northwestern Patagonia (r > 0.4; Lara et al., 2008).

The Puelo River is the main source of freshwater input of the Reloncaví Fjord, the northernmost fjord on Chile’s coast and one of the coastal systems most intensively used for aquaculture in the last decades. Chile is one of the world’s 10 top aquaculture producers, the second-largest exporter of salmon and trout, and the largest exporter of mussel (Mytilus chilensis). The Reloncaví Fjord has historically been used for smoltification and growth of salmonids by industrial salmon farming (Quiñones et al., 2019; Soto et al., 2019). This fjord is also essential for the seed production of Chilean mussel farming (Molinet et al., 2015). The fjord has a J-shape; it is 55 km long and less than 3 km near the head. The bathymetry of the Reloncaví Fjord shows depths >400 m in the area close to the mouth with a deep sill (∼ 150 m depth) located 15 km inland (Valle-Levinson et al., 2007). This fjord has a three-layer circulation pattern with an upper brackish layer flowing into the mouth and is strongly associated with freshwater inputs from the Puelo River and other tributaries (Cochamoì River Q = 100 m³ s^–1; Petrohueì River Q = 250 m³ s^–1). The system’s stratification varies along the fjord and between seasons, decreasing toward the mouth and being maximum in winter. The Reloncaví Fjord has been extensively studied by several research projects, mainly oceanographic, including paleoclimate (Rebolledo et al., 2015), fjord hydrodynamics (Valle-Levinson et al., 2014; Castillo et al., 2016), fjord biogeochemistry (Farías et al., 2017; González et al., 2019; Vergara-Jara et al., 2019; Yevenes et al., 2019) and plankton biology (González et al., 2013; Iriarte et al., 2017).

Data

Hydrology and Climate

Daily instrumental data of precipitation, air temperature, and streamflow were obtained from the Dirección General de Aguas de Chile, Dirección Meteorológica de Chile, Subsecretaría de Recursos Hídricos de Argentina and the Servicio Meteorológico Nacional de Argentina. According to these records, the Puelo River Basin contains one of the largest and most extensive networks of stations in western Patagonia. However, the location of the stations is spatially heterogeneous (Figure 1). Given this limitation, precipitation and temperature monthly data were complemented with CHIRPSv2 (2000–2019; 0.05°; Funk et al., 2015) and MODIS MOD11C3 (2000–2019; 0.05°; Wan et al., 2015) products, respectively. The associated biases between instrumental data and gridded products were corrected with different methodologies detailed in Aguayo et al. (2019). Wind and relative humidity monthly data were directly obtained from the ERA5 atmospheric reanalysis (2000–2019; 0.25°; Hersbach et al., 2020).

Water Quality

Six sampling sites were chosen within the drainage network to characterize the water quality of the Puelo River (Figure 1 and Table 1). Sampling sites were located downstream of Lake Inferior (P1), before and after (P2, P3) the confluences with its main tributaries, the Ventisquero (V1) and Manso (M1) rivers, and downstream of Lake Tagua-Tagua (P4) before its discharge in the Reloncaví Fjord (Figure 1). One additional sampling site was established in the Reloncaví Fjord near the Puelo River mouth (Figure 1 and Table 1). At this site the samples were extracted at 1, 8, and 25 m depth. The 1 m depth sample was intended to represent the upper brackish layer of the fjord (strong freshwater influence). The 25 m depth sample was designed to represent the fjord’s middle layer (below the pycnocline). Finally, the 8 m depth sample aims to reflect the seasonal changes in the system stratification, which, as mentioned, is highly dependent on freshwater inputs and other forcing factors such as tides and wind (Valle-Levinson et al., 2014; Castillo et al., 2016). According to continuous measurements carried out between May and September, 2018 these depths presented mean salinity values of 7, 28 and 30, respectively. Each sampling site was visited six times between 2018 and 2019 (42 site-date combinations). The timing of each field campaign was determined based on the monthly Flow Duration Curve (FDC) of the Puelo River during the period 1990–2019 (Figure 2). Considering the distance between sampling sites and the characteristics of the study area (Figure 1), each campaign was sampled within 3 days, which allowed capturing the variability of short-term events such as floods (Figure 2).

The water samples collected in field campaigns were used to analyze nutrients (N-NH₄, N-NO₃, N-NO₂, N-TOTAL, P-PO₄, P-TOTAL), dissolved organic and inorganic carbon concentration (i.e., DOC and DIC), DOM optical spectroscopy, and dissolved iron. For nutrient analysis, the water samples were collected (0.25 L) and stored at −18°C for further analysis according to standard methods established by American Public Health Association, 2005 in the Institute of Marine and Limnological Sciences, Faculty of Sciences, Universidad Austral de Chile.

N-NO₃ was determined by the cadmium reduction method (APHA 4500-E), whereas N-NO₂ by diazotizing with sulfanilamide and coupling with N-(1-naphthyl)-ethylenediamine dihydrochloride (APHA 4500-NO₂ – B). Ammonium was determined by the phenate method (APHA 4500-NH₃ F), total nitrogen was determined by the sodium hydroxide and persulfate digestion method (APHA 4500-N/C) followed by the determination of total N-NO₃. Total phosphorous (Total-P) was measured using the sodium hydroxide and persulfate digestion method (4500-P B/5) followed by the determination of soluble phosphorus using the method of acid ascorbic - blue indophenol (APHA 4500-PE).

Samples for DOC, DIC and DOM quality were filtered (GFF, nominal pore sizes = 0.7 and 0.22 μm) to analyze DOM optical spectroscopy, and DOC replicates were filtered and fixed by adding 100 μL of fuming HCl (Merck) to stop microbial activity, and transported at 4–7°C to the laboratory for subsequent analysis. Note that samples from the deepest point in the fjord (F25) were not available for DOM. DOC concentrations were measured using high-temperature catalytic oxidation (HighTOC, Elementar Systems) as described in Kamjunke et al. (2017). Fluorescence of DOM was measured using a Varian Cary Eclypse fluorescence spectrometer (Santa Clara, CA, United States) as described in Nimptsch et al. (2014). Briefly, excitation matrices from 240 to 450 nm (5 nm steps) and emission from 300 to 600 nm, inner filter corrections and Raman standardization were done using the FDOMcorr toolbox v1.6 for Matlab (Murphy et al., 2010). Parallel Factor Analysis Model (PARAFAC) was performed using DOMFluor toolbox v1.7 for Matlab (Stedmon and Bro, 2008). Four PARAFAC components (i.e. fluorophores C1 to C5) were identified and split-half validated with 1000 iterations with random starts (Stedmon and Bro, 2008).

For the determination of dissolved iron, water samples (volume: 15 ml) were filtered using a 0.2 μm syringe filter and acidified with suprapure nitric acid (MERCK). The samples were analyzed using total-reflection X-ray fluorescence spectrometry (TRXF) following procedures described in Mages et al. (2003). Acidified water samples (10 μl) were prepared onto quartz carriers and internal Ga standard (5 ng, suspended in Suprapur^® nitric acid for trace analysis—Sigma-Aldrich) was added. After drying on a hot plate (15 min, 60°C), trace elements were determined using a total reflection X-ray fluorescence spectrometer Picotax (BRUKER).

Land Cover Change

Satellite images were used to analyze the land cover change in the Puelo River Basin. The images were obtained from Landsat 5 and 8 for 2001 (December 16) and 2016 (February 9), respectively. These images have a spatial resolution of 30 m. Using the IDRISI GIS Analysis tool¹, selected images were subjected to standard preprocessing procedure as geometric, radiometric and topographic corrections to reduce the effects of atmosphere and shadows on land surface spectral response (e.g., Echeverria et al., 2006; Fuentes et al., 2017).

The supervised land cover classification of each image was performed in ENVI 4.5 software, which uses maximum likelihood statistics based on training points (Segura and Trincado, 2003). Training point selection for digital supervised classification of the 2001 image were taken from CONAF et al. (1999) dataset, a GIS-based data set of thematic maps derived from aerial photographs and satellite imagery from 1998 (Miranda et al., 2018). For the 2016 image classification, 750 training points were obtained from different field campaigns between 2018 and 2019. To increase accuracy of the satellite image classification, the Normalized Difference Vegetation Index (NDVI), Soil Adjusted Vegetation Index (SAVI) and Land Surface Water Index (LSWI) were used as additional spectral bands in both years (Wu et al., 2007).

The spectral information was classified in 10 different land covers categories: primary forests (pristine evergreen forests or those of a natural succession), secondary forest (forests regenerated after a disturbance, whether natural or anthropic), stunted forest (forest under the tree line), exotic tree plantation (commercial plantations mainly of eucalyptus), shrublands, grassland (crop and pastures), urban, water body, snow and ice and bare land. The accuracy assessment consisted of confusion matrices between classified images and independent sampling points collected from aerial photographs and field campaigns. Land cover change analyses were performed using the IDRISI Selva Land Change Modeler module². This analysis consisted of descriptive spatial-temporal statistics of the land covers gains and losses in the time period. The change rates (q) of the different classes were calculated using the equation proposed by Food and Agriculture Organization (1996):

where A₁ and A₂ are the land cover at time t₁ and t₂, respectively.

Wildfires

Past burned areas (>1 ha; 1985–2018) were analyzed as the main potential driver of land cover change in the Puelo River Basin. Burned area reconstruction in the Chilean portion of the basin was performed by applying an algorithm in Google Earth Engine (GEE) (Long et al., 2019). GEE is an open cloud computing platform for geospatial analysis that contains a public catalog of historical satellite images, topography, land cover and other environmental datasets (Gorelick et al., 2017). Taking advantage of the GEE big-data analysis platform, we develop a flexible workflow to reconstruct individual burned area for all fires reported since 1985. This approach processes Landsat images and generates historical burned areas by detecting the multi-temporal reflectance change (before and after the fire) of the Normalized Burn Ratio [NBR = (ρNIR - ρSWIR2)/(ρNIR + ρSWIR2)], which uses the difference between pre- and post-fire reflectance of the near-infrared band (ρNIR) and short-wavelength infrared band (ρSWIR2) extracted from satellite images (Key and Benson, 2003). High values of NBR indicate a burned area that allows the reconstruction of the fire scar in the landscape. All the necessary data are freely available at GEE. This process needs as initial data the ignition point and date, which were compiled by the Chilean National Forestry Corporation (CONAF, 2018). In the Argentinean portion of the basin, the burned areas were obtained directly from MODIS MCD64A1 monthly data (500 m; 2001–2019) and the program of Provincial Agricultural Services (1985–2000).

Statistical Analysis

Water Quality

Nutrient and DOM data were analyzed with a multivariate approach focused on identifying spatial-temporal patterns of coherence between streamflow and water quality variables. Nutrient, DOC, DIC, and DOM data derived from spectroscopic characterization (CDOM (m^–1), color (Pt) and PARAFAC components) were first tested and transformed to meet the normality requirement as appropriate. The relationship between nutrients and DOM variables and between samples from all the sites and campaigns was evaluated using a heatmap and a hierarchical clustering approach on the Euclidean distances of the data matrix (ComplexHeatmap package v1.10 for R; Gu et al., 2016), which identifies clusters among and within samples and variables.

A non-metric dimensional scaling (NMDS) was used on the Manhattan distance matrix of nutrient concentration (N-total, P-total, NH₄, NO₃, NO₂, and PO₄) and DOM quality descriptors (PARAFAC components, CDOM, ColorPt, and DOC concentration). NH₄ and NO₂ were excluded from the analysis as more than 70% of the cases were below the detection limit. The effect of having total P and PO₄ samples below the detection limit in 50 and 30% of the cases was evaluated and the detection limit value was used (Palarea-Albaladejo and Martín-Fernández, 2015). The relationship between the NMDS of nutrients and DOM quality descriptors and streamflow was assessed by fitting streamflow onto the ordination (envfit function) and the result was plotted on the ordination diagram. The significance of the correlation was assessed through a Monte Carlo test (1,000 permutations). Additionally, the grouping of the factors “Campaign” and “Site” was visualized through hulls plotted for each group using the function ordihull and the average of the group scores for each hull (function ordibar). The area of those polygons was then tested against areas produced by randomized groups through a permutation test (1,000 permutations) using the function ordiareatest. The NMDS was generated by MetaMDS function. To assess the effect of the factors campaign and site on DOM quality descriptors, permutation analysis of variance (PERMANOVA; Anderson, 2001) was performed on the corresponding Manhattan dissimilarity matrix (adonis2 function). All the analyses described above were performed using stats v4.1 and vegan v2.5 packages for R.

Hydrological Regime

Trends in the hydrological regime of the Puelo River were examined using the Indicators of Hydrologic Alteration software (IHA; Richter et al., 1996). Based on daily data, this tool calculates 33 hydrological metrics that characterize the intra- and inter-annual variability in streamflow conditions, including the magnitude, frequency, duration, timing, and rate of change of streamflow. The metrics were estimated at the Puelo River fluviometric station located near the mouth (PD in Figure 1). This station has the most extensive records in the basin and possibly the longest in western Patagonia (water years: 1950–2019). This station adequately represents the hydrological behavior of the streamflows recorded by stations located upstream in Chile (Pearson correlation > 0.95 for the Manso and Puelo Rivers before their confluence; Figure 1). Missing daily data were filled with linear regressions with upstream stations (<1% of total data). The magnitude and the significance of the trends were analyzed with the Sen’s slope (Sen, 1968) and Mann-Kendall tests (Mann, 1945) (trend package v1.1 for R).

Cross-wavelet transform (XWT) and coherence (XWC) analyses were applied to identify significant changes in the association between time series of precipitation and streamflow (Grinsted et al., 2004; Cazelles et al., 2008). XWT evaluates when two time series oscillate in common periods and if there is temporal lag between the peaks of the oscillations, and XWC analyzes the strength of this association (in the form of a correlation and a significance level). A high coherence (∼correlation) between two series for a given oscillation frequency and time period implies that the two series oscillate in a common frequency, and also that the lag between the peaks of the oscillations is constant, which suggests an underlying causal mechanism (Grinsted et al., 2004). XWT was calculated from log (streamflow; PD station in Figure 1) or square-root (precipitation; station located near Puerto Montt city 41.5°S), transformed series to attain normality. XWC was calculated in the Undarius High Performance Computer cluster at the Catalan Institute for Water Research, using a Monte Carlo randomization technique as in Grinsted et al. (2004). After XWT and XWC calculations, time series were extracted from correlation and the phase relationship of the most prominent common periods of oscillation shared by the two-time series, and calculated annual means. All data points that did not correspond to statistically significant (p < 0.05) associations as identified during XWT were discarded. The correlation and phase relationship of the same common periods of oscillation were extracted for each season and year. All calculations were performed using the biwavelet package v0.20.19 for R.

Hydrologic Modeling of LCC

Land Cover Scenarios

The hydrological model aims to evaluate the hydrological impact of three different land cover scenarios (SC), projected to the year 2030 based on observed LCC (2001–2016) and historical burned area (1985–2019). The first scenario (SC-1) considered a matrix of transition probability generated from LCC from 2001 to 2016 in the CA-MARKOV Module of IDRISI Selva v17.0 software³. This stochastic model assumes that in a given period pixel persistence or change in either class is dependent only on the immediately previous state, not on historical changes (Zhang and Dai, 2007; Zhang et al., 2011). A cellular automaton model was used to spatialize the Markov model results. This model allowed assessing the spatial changes among pixels, assuming that the changes are dependent on the initial condition and neighboring pixels (Guan et al., 2011; Yang et al., 2012). A standard 5 × 5 contiguity filter type cellular automaton and a 14-year iteration were used to generate the projected land cover for 2030. The second scenario (SC-2) assumed that the zones most likely to be burned are those where wildfires have been previously occurred. Under this assumption, SC-2 replaced SC-1 land cover with burned areas in those places where wildfires have occurred historically (1985–2019). Finally, the third scenario (SC-3) added to SC-2 new burned areas according to a 2 km buffer around the main lakes and urban or touristic spots, since most of the fires have been recorded in the vicinity of these places [see section “Land Cover Change (LCC)”]. The areas affected by wildfires were considered as bare land in the modeling process in order to evaluate the most extreme hydrological condition. This was the most recurrent condition detected during the post-fire analysis of LCC [see section “Land Cover Change (LCC)”]. Finally, the land cover classes were grouped into six categories: forest (primary, secondary, stunted, and exotic tree plantations), shrubland, grassland, bare land (including urban areas), water body, and snow and ice.

Hydrological Model WEAP

The three land cover scenarios for 2030 were simulated in the WEAP hydrological model (Water Evaluation and Planning; Yates et al., 2005), considering the same climatic conditions for each. The WEAP model, previously used in Chile (e.g., Vicuña et al., 2012; Chadwick et al., 2020; Barría et al., 2021), is a semi-distributed model that represents the relevant hydrological processes using empirical functions that describe the distribution of water in two soil water storages (root and deep storage). In this regard, the snow accumulation and snowmelt are based on the degree-day method, and the potential evapotranspiration (PET) was calculated with the Penman-Monteith equation. As a final output, the streamflow composition of each simulated river corresponds to the sum of surface runoff, interflow and baseflow. See Yates et al. (2005) for more details on the equations of the conceptual model.

The WEAP model was previously calibrated and validated for a monthly time step by Aguayo et al. (2021) using the corrected atmospheric products of section “Hydrology and Climate” (period 2000–2019). In this approach, the Puelo River Basin was divided into nine sub-basins (green squares in Figure 1), which were subdivided into different hydrological units according to their elevation and land cover type (baseline year 2001). The model reached a Modified Kling-Gupta Efficiency (Kling et al., 2012) of 0.83 and 0.74 for the calibration (2000–2010) and validation (2011–2019) stages, respectively (sub-basin PD in Figure 1). The model maintained high correlations (r = 0.78 ± 0.07) and dry biases (β = 0.93 ± 0.11) during the validation stage. These biases were mainly associated with high streamflow events when the probability of exceedance was less than 20%. In contrast, minimum annual streamflows were adequately simulated (r = 0.84). The simulated PET near the Reloncaví Fjord was similar to the time series recorded by an evaporation pan located in the same place (2002–2012; R² = 0.87, RMSE = 14 mm). Further details of the calibration/validation process and the parameter estimation can be found in Text S1 and Aguayo et al. (2021).

Results

Water Quality in the Land-Ocean Interface

Nitrogen, phosphorus, dissolved carbon, and iron measured during the six field campaigns in the Puelo River Basin had very low concentrations (Figure 3 and Supplementary Table 2) which are characteristic of well-conserved basins previously reported in southern Chile (Supplementary Table 3). In river samples (P1 to P4, Manso and Ventisquero), total nitrogen, phosphorus, dissolved carbon and iron concentration were mostly < 60 μg N L^–1, < 10 μg P L^–1, < 3.8 mg C L^–1 and < 10 μg Fe L^–1. In contrast, the Reloncaví Fjord had mostly 3–15 times higher mean N, P, C and Fe concentrations (160–400 μg N L^–1, 27–74 μg P L^–1, 6–13 mg C L^–1, 22–26 μg Fe L^–1) which also increased remarkably by a factor of three from surface to 25 m depth (Figure 3 and Supplementary Table 2). High values of total dissolved carbon in the fjord corresponded to the inorganic fraction; DOC concentrations were comparable across samples (Supplementary Table 2). The Ventisquero (V1) and Manso Rivers (M1) presented different mean values from those found by site in Puelo River’s mainstem (Supplementary Table 2). The Ventisquero River showed the lowest mean total N and dissolved C concentration (<30 μg N L^–1 and < 2,14 mg C L^–1), thus diluting nutrient concentrations of the mainstem after their confluence. Manso River transported higher nutrient loads, especially TP, leading to increases in P3 (Figure 3, Supplementary Figure 1, and Supplementary Table 2).

The NMDS separated the samples according to their origin of fjord versus basin and into the different sampling campaigns (Figure 4 and Supplementary Figure 1; stress = 0.159). The first dimension was related, on the positive side, to protein-like DOM components, and with nutrients and DIC on the negative side. Fjord samples (F1 and F8) appeared grouped independently of the sampling campaign, and related to nutrients. The fjord’s surface layer samples (1 m, F1) appear between the deeper fjord samples (8 m, F8) and the basin samples. The second dimension of the NMDS separated the basin samples by their different campaigns based mainly on DOM descriptors. Thus this dimension represents a gradient from colored and humic to protein-like DOM. C2 (Aug. 19) samples presented colored DOM, while those from C1 (May 19) had protein-like DOM. The remaining campaigns appeared ordered in the intermediate values of that second dimension. Streamflow appeared significantly fit to the NMDS ordination (r² = 0.19, p = 0.009), with higher streamflow conditions related to positive values of the second axis. Thus, these conditions (Figure 2) contribute with humic and colored materials (as in August, 2019), while low streamflows corresponded to higher arrival of protein-like materials from the basin (as in May, 2018). Fjord samples from 1 m are grouped together with Puelo River samples for each campaign. The permutation test for the size of the ordination hulls for campaign show that all the campaigns except C3 had smaller areas than random groups, while the test for sites shows that only F8 had a smaller area than random groups (Supplementary Table 4). In agreement with this, the PERMANOVA analysis of nutrients and DOM showed that both campaign and site had a significant effect (Table 2; FCampaign = 5.98, FSite = 10.32, p = 0.001 in both cases), although the effect of site was weaker when nutrients were excluded (Table 2; FSite = 3.26, p = 0.005, FCampaign = 10.12, p = 0.001).

Natural Flow Regime

Puelo River observational records analysis showed trends in magnitude, timing, duration, and change rate of streamflow (Supplementary Table 1 and Figure 5). The streamflow recorded during summer and autumn showed statistically significant decreasing trends (summer in color palette of Figure 5A). The most affected months were January (−6% decade^–1; p < 0.001), February (−6% decade^–1 p < 0.001), March (−4% decade^–1; p < 0.05) and May (−9% decade^–1; p < 0.05). The magnitude of such trends was also reflected in the Flow Duration Curve for the 1950–1979 and 1990–2019 periods (Figure 2). Inter-annual hydrological variability also showed significant changes associated with extreme conditions. Minimum annual streamflows for time windows of 7, 30, and 90 days showed statistically significant trends (Supplementary Table 1 and Figure 5A). The timing of the maximum annual peak streamflow also changed. Only 10% of the years between 1950 and 1990 recorded the peak outside of the autumn or winter. In contrast, since 1990 this percentage has increased to 40% (p < 0.05; Figure 5B), and it has become more frequent to observe them in spring. Furthermore, the duration of the high streamflows registered a significant decrease (0.3 days decade^–1; color palette of Figure 5B). Finally, there were also significant trends in the streamflow fall rate (negative differences between consecutive values; Supplementary Table 1). The anomalies in the natural flow regime have coincided with the SAM trend (Figure 5C). In recent decades, the SAM positive phase has coincided with very warm ENSOs, which has favored extreme drought conditions in Northwestern Patagonia (e.g., 1998 and 2016 in Figure 5). This is not expected, considering that El Niño conditions tend to promote the negative phase of SAM, thus producing a negative correlation between their indices at interannual time-scales (L’Heureux and Thompson, 2006; Ding et al., 2012).

The cross-wavelet transform and coherence analyses suggested shifts in the response of streamflow to precipitation events during the last decades (Figure 6). Beyond the expected year-to-year variability of the correlation (color) and lag (arrows) between streamflow and precipitation at short periods of oscillation (period from weeks to ∼3 months; Figure 6A), the analysis showed a significant common oscillation at 180 days (Figure 6B), corresponding to the relationship of precipitation peak during early winter and the streamflow peaks during early winter and late spring. For the oscillation at a period of 180 days, a significant increase (tau = 0.224, p < 0.05) of the lag between precipitation and streamflow was apparent (Figure 6C), from oscillations mostly around phase (lag ∼ 0, i.e. oscillations are synchronous) during the 1970s, to average lags around 20 days (and up to ∼50 days) during the last decade (i.e. streamflow peak had a lag of ∼20 days with respect to the precipitation peak), which is consistent with Figure 5B. This indicated that years with the prominent streamflow peak during spring are becoming more frequent. The correlation between the time-series at the period of oscillation of 180 days substantially decreased over years, showing longer lags between precipitation and streamflow (Figure 6C). This interpretation is also supported by the correlation between precipitation and streamflow oscillations analyzed at a period of 1 year, which shows consistently high correlation during the last decades (Figure 6B). This is expected if the double streamflow peak during winter-spring simplifies to a single peak, which would correlate better with the prominent precipitation peak during winter.

Land Cover Change (LCC)

The overall accuracy for satellite image classification was 85.3% for 2001 and 88.7% for 2016. The highest accuracy (except for water and bare soil) was achieved in primary forests (2001: 84.6%; 2016: 87.1%). The lowest values of accuracy corresponded to shrubland in 2001 and 2016 (76.6 and 81.8%, respectively). Land cover classification showed that 55% of the Puelo River Basin was covered by native forest in 2001 (Figure 7c and Table 3). Native forest cover rises to 77% without considering the water bodies and high elevated areas of bare land and snow. Other important classes were shrubland (13.5%) and bare land (18.7%) (Tables 1, 3).

Land cover change analysis showed that 91.2% of the total area remained unchanged in the 2001–2016 period (Figure 7a and Table 3). The main changes corresponded to the increase of secondary forest (+3.7% yr^–1), shrubland (+0.2% yr^–1) and bare land (+0.5% yr^–1) (annual rate of change; Table 3). Secondary forest recovery was concentrated in areas with reduced human activity (e.g., Ventisquero River basin; Figure 7a). In contrast, the degradation processes associated with the loss of primary forest were located in the Turbio and Epuyén river valleys, where historical wildfires were concentrated (Figures 7A,B). Although exotic tree plantations and urban areas represented a minimum area in 2016, these land cover classes reported the highest annual rates of change (7.0 and 5.4%, respectively; Table 3).

The reconstruction of the burned area in the Puelo River Basin showed that 649.6 km² (7.1% of the basin) were affected by 74 wildfires between 1985 and 2019 (Figure 7b and Table 1). Of this total, 43 records occurred in Chile (109 km²) and 31 in Argentina (539 km²). Only 15 wildfires accumulated 88% of the total burned area (>10 km²). Of these, only one occurred before 1995 (Table 1). The largest wildfires in Chile and Argentina occurred in the surroundings of Tagua-Tagua Lake (1998; 39 km²) and Puelo Lake (2015; 129 km²), respectively.

Hydrological Response to LCC

The three land cover scenarios projected toward the year 2030 showed LCC of different magnitude. SC-1 will decrease in primary (−1.1%) and stunted forest (−0.4%) (Table 3). SC-2 will add 524 km² of bare land in areas historically affected by wildfires (1985–2019) to SC-1, which affected mainly native forest areas (148 km²), shrubland (197 km²) and grasslands (24 km²) (Figure 7b and Table 3). Finally, SC-3 will add 1,192 km² of possible burned areas to SC-2 according to a 2 km buffer around the principal lakes (Figure 7d and Table 3). Despite the extreme conditions assumed by this scenario, the native forest will decrease from 55.2 to 45.7%, while the bare land will increase from 18.7 to 31.6% (Table 3).

The decrease in vegetation cover shown by the three scenarios (e.g., native forest, shrubland) consistently reduced ET, and altered internal hydrological processes (Figure 8). Reduced vegetation promotes changes in surface hydrological response, reducing infiltration rates and increasing surface runoff. Lower infiltration reduces subsurface flow and percolation into deep storage, which in turn controls base flow. Note that the base flow dominated the streamflow composition during the baseline period (2000–2019) with values of 40 ± 9% (Figure 8).

The magnitude of the hydrological response varied according to the scenarios, from non-significant changes (SC-1) to slight changes (SC-3) (Figure 8). Note that snow accumulation remained constant in Figure 8, since the climatic forcing does not change in the land cover scenarios. Regardless of the scenario evaluated, surface runoff and ET showed the greatest changes (Figure 8). Under the worst-case scenario (SC-3), the winter season recorded the largest projected changes in streamflow composition. In this season, the base flow, interflow and runoff showed changes relative to the baseline period of −2.0, −2.4, and 8.6%, respectively (Figure 8). Despite the seasonal changes in hydrological processes, the same scenario indicated an average annual increase of only 1.1% of freshwater input to the Reloncaví Fjord.

Discussion

This study illustrates the importance of water availability and quality of Puelo River Basin as the main freshwater input to the Reloncaví Fjord system. We show that the hydrology of the Puelo River has changed over time, with major shifts in magnitude (lower flows in summer and autumn), duration (minimum annual streamflows), timing (more floods in spring), and frequency (fewer prolonged floods). Concentrations of nitrogen and phosphorus must be highlighted as some of the lowest recorded for temperate rivers, highlighting the relatively pristine nature of the basin. The analysis of spectroscopic DOM descriptors and DOC concentration showed that streamflow and hydrological moment drive DOM composition in the hydrological network. The freshwater inputs from the Puelo River are fundamental to the haline stratification of the Reloncaví Fjord and the dilution of nutrients in the upper brackish layer. Wildfires have been the main driver of land cover change (LCC) in the past. Future LCC projections will reduce the evapotranspiration and subsurface flow, increasing the surface runoff, which will result in slight seasonal changes in freshwater inputs to the fjord. This study highlights the need for large coastal basins with near-reference conditions as indicators, sentinels or model systems for evaluating how could affect future scenarios of global environmental change affecting freshwater inputs (including volume, seasonal, and water quality) in fjord systems of Patagonia.

Natural Flow Regime and Freshwater Quality

Nutrient concentrations in the Puelo River drainage network across seasons are in the lower range of previously reported values for rivers in southern Chile (Supplementary Table 3; Oyarzun et al., 1997, 2004; Godoy et al., 1999, 2001; Perakis and Hedin, 2002; Oyarzún and Huber, 2003; Little et al., 2008), and elsewhere (de la Crétaz and Barten, 2007). Nutrient concentrations in the Reloncaví Fjord are much higher than those from the Puelo River. Therefore, the low concentration of N and P in these freshwater inputs of the Puelo River may play a significant role in the coastal dynamics, influencing the upper layer of the Reloncaví Fjord. Clearly it can be stated out that the enhanced concentrations of dissolved and total fractions of N, P and Fe reported in the fjord are unlikely to be caused by the Puelo River input, as its nutrient and Fe content were always below the fjord samples. Therefore, it can only be speculated what other sources are responsible for the higher N, P and Fe concentration in the fjord; diffuse input from coastal sides, in-sea contamination from salmon farming and others. The NMDS and PERMANOVA analysis showed that nutrients present stronger differences between the fjord and the basin samples than DOM composition and the fjord’s surface layer (1 m, mainly fresh and brackish water) nutrient concentrations are between the deeper fjord samples (measured at 8 and 25 m depth in the fjord) and the basin samples (Supplementary Table 2). This upper layer can limit the exchanges with deeper layers mostly of oceanic origin (González et al., 2013; Torres et al., 2014). This fjord has likely been more resilient to additional nutrient inputs (salmon farming; Soto et al., 2019) due to both the dilution role of freshwater input and rates of water turnover (Pinilla et al., 2020). Future scenarios of streamflow reduction might reduce water turnover in the fjord and enhance local eutrophication conditions. Hence the pristine conditions of freshwater input could provide important ecosystem services to maintain the resilience of the Reloncaví Fjord ecosystem to global environmental change.

NMDS based on nutrients, spectroscopic DOM descriptors and DOC concentrations (Figure 4) shows that streamflow and hydrological conditions drive DOM composition in freshwater. A closer examination of the origin of DOM from the Puelo River to the Reloncaví Fjord during the low streamflow period illustrates the fragile balance between hydrology and biogeochemical function of the river-fjord system. The PERMANOVA result, with both campaign and site being significant, does not allow us to present hydrology as the only control of DOM quality, but this result together with the NMDS sample clustering leads us to conclude that (i) seasonality and streamflow play a major role in it, (ii) the two major tributaries M1 and V1 contribute with low nutrient concentrations and very variable DOM quality (Figure 4 and Supplementary Table 4).

High streamflow conditions during winter are related to humic and colored DOM, in agreement with behavior as a passive pipe; terrestrial DOM is transferred from soils to streams, and then rapidly transported toward the ocean without significant alteration (Raymond et al., 2016; Casas-Ruíz et al., 2020). In contrast, during low streamflow conditions we found an increase in protein-like DOM sources from degraded peptide material (Fellman et al., 2010) coming from the Puelo River Basin. This is likely due to the disconnection of lateral channel processes in the drainage network resulting in longer water residence time, and a large proportion of DOM from autochthonous sources (Casas-Ruiz et al., 2016; Catalán et al., 2016). However, we found a significant effect of site on DOM quality, and we expect some sections of the Puelo River drainage network to have a disproportionate influence on DOM quality and ultimately on the different layers of the Reloncaví Fjord. In order to constrain the link between hydrology and DOC concentration and DOM quality better, we encourage future studies to trace both the hydrograph and concentration response on a finer temporal scale and especially during extreme hydrologic conditions, as those may lead to especially harmful conditions in the fjord.

While our results help to understand the effect that the interplay between hydrological change and LCCs has on DOM composition, it is still unclear whether the relatively small LCC occurring in the Puelo River Basin have sufficient potential to affect fjord-level processes significantly. Targeting compounds at the molecular level (Coppola et al., 2018) might be key to identify the legacy of wildfires in these aquatic systems further. For example, increases in riverine export of dissolved black carbon have been found decades after wildfires in Brazil’s Atlantic forest (Dittmar et al., 2012). Changes in the composition of riverine organic matter can alter carbon cycling in fjords due to shifts in both sedimentation rates and light penetration (Ward et al., 2017). However, even considering the constraints of our study, based on the high proportion of native forest in the Manso sub-basin, the projected LCC may reduce the DOM inputs from terrestrial vegetation (colored and typically aromatic) in the Reloncaví Fjord. This could lead to reduced flocculation and sedimentation in the coastal waters (Raymond and Spencer, 2015).

Recent hydrological shifts imply that the routing of water in the basin is undergoing profound changes. Considering the pristine nature of the basin, this also may have measurable yet unknown consequences for the export of DOC and nutrients from soils to rivers, and ultimately for the fjord. Although most studies focus on temperate, boreal and artic regions, it is well known that changes in hydrological regimes influence the export of DOC in relatively undisturbed basins (Winterdahl et al., 2016; Baek et al., 2019; Connolly et al., 2020). Unfortunately, our short-term study in terms of water chemistry cannot illuminate this issue, and there are no long-term records of organic matter and nutrient export from Patagonian basins to support claims related to changes in hydrological regimes. However, there is evidence from other regions that DOM export from snow-melt has a different chemical signature than exports from other flow paths (Burns et al., 2016). Furthermore, changes in timing and magnitude of streamflow in the region represent shifts in natural hydrological regimes (Poff et al., 1997; Poff, 2018), and could have additional biological implications in river networks and fjord ecosystems of this region. For example, the predicted shift in the hydrology of the Puelo River might increase the importance of autochthonous sources and decrease the transport of DOM from humic-like terrestrial sources reaching the fjord, except during storm episodes. Considering the relevance of organic carbon and nutrient exports for the ecological dynamics of the fjord and the expected hydrological changes in the future, understanding how this hydrological shift impacts constituent export from the basin should be a research priority.

Wildfires Are the Main Drivers of Land Cover Changes

The native forest cover in the Puelo River Basin (primary, secondary and stunted forests) showed an annual deforestation rate of 0.03%, with more than 90% of the basin area remaining unchanged during the period 2001–2016 (Table 3). As in Gowda et al. (2012), urban areas and exotic tree plantation had a marginal negative net contribution to forest cover throughout the period. However, these land cover classes reported the highest annual rates of change (7 and 5.4%, respectively, Table 3). The low LCC rate demonstrates the good conservation status of the Puelo River basin, differing from the values described for south-central Chile (33–42°S), an area that has been subjected to major LCC and is experiencing an annual rate of forest loss ranging between 0 and 5.8% (Miranda et al., 2017). Also the LCC ranges of the Puelo River basin are in the lower range of northwestern Patagonia, where Echeverría et al. (2012) reported an annual deforestation rate close to 1.0% (1985–1999 vs. 1999–2007).

According to our results, wildfires are main drivers of LCC in the Puelo River Basin. The frequency of fires in northwestern Patagonia peaked in the late nineteenth century, due to widespread burning and clearing of forests by European settlers. Fire frequency declined dramatically around 1910 due to the cessation of intentional fires; it has remained low due to current fire suppression (Kitzberger and Veblen, 1999; Veblen et al., 1999). In this study, the historical reconstruction of burned areas in the Puelo River Basin shows that 649.6 km² were affected by 74 wildfires between 1985 and 2019 (Figure 7 and Table 1). Mundo et al. (2013) reported a strong negative association between the density of fire ignitions and their distance to cities in this region (<5 km). This is consistent with the human origin of the fires in the Puelo River Basin, where most of wildfires have occurred around national parks, lakes and urban areas. Although the annual area burned in the Puelo River Basin has not changed over time, more than 90% of wildfires larger than 10 km² occurred after 1995. Climate trends might facilitate the spread of large wildfires, which added to the high fuel content in the region could increase the frequency of larger wildfires (Kitzberger and Veblen, 1999; González et al., 2018; Urrutia-Jalabert et al., 2018). In comparable mountain systems from western North America there is a strong link between wildfire activity and climate oscillations (Trouet et al., 2006; Holden et al., 2018).

Synergistic Hydrological Effects Between Climate and LCC

The three LCC scenarios showed changes of different magnitude in the Puelo River landscape. If current LCC trends continue (SC-1), native forest will decrease by 1.6% (relative value for 2001–2030), while if LCC intensifies due to potential wildfires (SC-3), native forest could decrease by up to 9.5% by 2030 (Table 3). Under these scenarios, bare land, the most recurrent land cover class after wildfires (Figure 7), would increase by 1.8 and 12.9%, respectively (Table 3). According to the WEAP model, these land cover transitions promote reductions in infiltration rates, which in turn reduce subsurface flow and percolation into deep storage (Figure 8). Reduced percolation can lead to a reduction in base flow, the dominant streamflow component of the Puelo River (40 ± 9%). Consistent with other hydrologic modeling in mountain systems (Havel et al., 2018), other processes observed in the hydrologic response were an increase in surface runoff and a decrease in evapotranspiration. The hydrologic response varied according to the scenario and remained similar between seasons. In the worst scenario for winter (SC-3), the hydrological model predicts changes of −2.0%, −2.4% and 8.6% on the base flow, interflow and runoff, respectively. The magnitude of projected changes is comparable to a deforestation scenario in the Vergara River Basin (Stehr et al., 2010), but lower than what was reported in south-central Chile, where changes in LCC are relatively more drastic and different (intensive forestry; Aguayo et al., 2016; Barrientos et al., 2020; Martínez-Retureta et al., 2020).

The hydrological simulations forced by LCC scenarios were an order of magnitude lower than the previous projections forced by the climate change scenarios. According to dendrochronological reconstructions, recent climate trends are unprecedented in the last centuries (e.g., SAM in Figure 5C; Villalba et al., 2012; Morales et al., 2020), which has been evident in the precipitation and temperature of northwestern Patagonia (Pabón-Caicedo et al., 2020). Recently, Aguayo et al. (2021) showed that temperature in northwestern Patagonia has increased during the whole year, but mainly in summer (0.3 ± 0.1°C decade^–1). In contrast, precipitation has decreased mainly during the autumn (−8 ± 8% decade^–1) and winter (−4 ± 9% decade^–1). Based on these climate trends, recent hydrological modeling in northwestern Patagonia projects an increase in the frequency of severe droughts during summer and autumn months (Aguayo et al., 2019, 2021; Pessacg et al., 2020). Lower streamflow conditions during summer are likely due to higher evapotranspiration (Ellenburg et al., 2018), and lower snow accumulation in snow-dominated regions (Barnett et al., 2005). Timing of the streamflow annual peaks have also changed in the Puelo River Basin, migrating from winter to less prolonged spring peaks. The increasing lag between precipitation and streamflow seasonally (from synchronous peaks to a lag of ∼50 days) suggests a shift in the peak streamflow timing (Figure 6). We hypothesize that the lower precipitation recorded in winter (maximum daily precipitation), together with a greater height of the 0°C isotherm in spring, allows the maximum annual streamflow to occur in spring.

The importance of various hydrological stressors depends on the local climatic, geographic and hydrological conditions, and the respective future trajectories (e.g., Davis et al., 2015; Yang et al., 2017; Chanapathi and Thatikonda, 2020). Hydrological projections for LCC and climate scenarios show opposite directions in the Puelo River Basin. Considering a moderate emissions scenario to the year 2030 (Shared Socio-economic Pathway Two), Aguayo et al. (2021) projected streamflow changes of −15, −18, 8, and 2% to DJF, MAM, JJA and SON, respectively (multi-model mean for 2020–2040 vs. 2000–2019). These values are an order of magnitude larger than those projected by LCC, since the projected climate trend for NP is one of the largest in the world (Cook et al., 2020). For example, Barría et al. (2021) recently determined that the relative impact of the climate factor is more than 10 times larger than the impact of the observed LLC changes in the Aculeo Lake balance (Central Chile; 33°S). Although these results may discourage conservation, Galleguillos et al. (2021) reported that conservative scenarios focused on native forest protection and restoration could partially mitigate the effect of climate change in Mediterranean basins. Future efforts should focus on understanding the hydrological role of different vegetation cover (e.g., Milkovic et al., 2019), under changing climatic conditions in northwestern Patagonia.

Long-Term Hydrological Trends Are Changing the Fjord System

Changes in climate ultimately result in modified trends of surface hydrology. During the last decades, the Puelo River regime showed a clear trend toward longer and more intense hydrological droughts, with major shifts in magnitude (lower flows in summer and autumn), duration (minimum annual streamflows), timing (more floods in spring), and frequency (fewer prolonged floods) (Figure 5). The trends in magnitude were consistent with those previously reported. Recently, Aguayo et al. (2021) showed a decreasing streamflow mainly during autumn (−6% ± 3% decade^–1) and summer (−5% ± 2% decade^–1) for most of the rivers of northwestern Patagonia. As a result of climate change, these results invite understanding the influence of the Puelo River as a non-stationary driver of the Reloncaví Fjord.

Changes in the natural flow regime during extreme climatic conditions, like during ENSO event in summer and autumn (e.g., Valle-Levinson et al., 2007; León-Muñoz et al., 2018), could influence the occurrence of anomalous oceanographic events in the Reloncaví Fjord. The reduction in freshwater input drives the weakening of ocean stratification in the upper layer, allowing vertical advection of saline, less oxygenated and nutrient-rich waters (Valle-Levinson et al., 2007). These mechanisms have the potential to affect the production of mussels (Mytilus chilensis) and salmon (Salmo salar and Onchorhynchus spp.) in aquaculture facilities (Soto et al., 2019). For example, an increase in the salinity of the fjord’s surface layer can alter environmental conditions that minimize the impacts of caligidosis (sea louse) in salmon (Montory et al., 2018). León-Muñoz et al. (2018) reported that the severe drought during 2016 (the lowest streamflow in 70 years) generated the appropriate conditions for a record harmful algae bloom in northwestern Patagonia. As a consequence, the mortality of salmon in a few days was comparable to the mortality that all salmon farming usually experiences in 2 years (US$ 800 million). This extreme event is also a relevant threat to mussel culture (Chávez et al., 2019; Soto et al., 2020). Coastal sites with high input of freshwater seem to ensure better mussel seed provision, survival and capture, apparently because the wild mussel populations providing the seeds are more successful at lower salinities (Molinet et al., 2015, Molinet et al., 2021). Therefore, the decreasing trend in freshwater input to the Reloncaví Fjord during drought years, and the changes in its timing, could be compromising the future of these two industries (Soto et al., 2020).

Conclusion

The Puelo River streamflow showed non-stationary conditions with clear trends toward longer and more intense droughts. Several indicators of the natural flow regime also showed significant trends, such as the timing of the maximum annual peak streamflow, the duration of high streamflows and the streamflow fall rate. Despite the recent increase in wildfires in recent decades, more than 90% of the total area maintained its land cover during the period 2001–2016. Considering these trends, the hydrological simulations based on land cover scenarios toward the year 2030 will slightly influence hydrological processes (e.g., surface runoff, evapotranspiration), and therefore the hydrological regime of the Puelo River will be mostly modulated by climate trends. These climate trends would imply hydrological changes in the opposite direction to those forced by the LCC scenarios. The combined effect of both hydrological stressors and potential decline in water quality (increase in nutrients) is a call for more robust data monitoring and analysis of both drivers. For example, it is essential to advance in the understanding the hydrological role of different vegetation covers under warmer and drier climatic conditions in northwestern Patagonia.

Nutrient concentrations in the drainage network are in the lower range of previously reported values for rivers in southern Chile, and are 3–15 times lower than in the fjord. As result, the fjord′s surface layer limits the exchanges with deeper water layers mostly of oceanic origin. Future decreases in freshwater input might reduce water turnover in the fjord and enhance local eutrophication conditions. Coastal sites with high freshwater input support large-scale salmon aquaculture and ensure better mussel seed survival. Considering the relevance of organic carbon and nutrient exports, understanding how this hydrological shift influences the export of basin constituents warrants further research. For example, to trace both the hydrograph and concentration response on a finer temporal scale during extreme hydrologic conditions, as these may lead to anomalous and extreme events in the fjord (e.g., harmful alga bloom, hypoxia events). This study provides baseline information about present and future threats relevant to the potential of “losing the freshwater-dependent resilience” of the Reloncaví Fjord to future global environmental change. More research is still needed to understand better the influence of water availability and quality on the biogeochemistry, food webs and ecosystem services provided by this river-fjord complex.

Data Availability Statement

The datasets generated for this study can be found in the Zenodo repository (https://doi.org/10.5281/zenodo.4604929). This repository contains hydro-climatic observations, gridded products of precipitation and temperature, land cover for different periods, and the reconstruction of burned areas. The raw data from CHIRPSv2, MODIS MOD11C3, and ERA5 can be found at www.chc.ucsb.edu/data/chirps (last access: June 12, 2021), https://lpdaac.usgs.gov/products/mod11c3v006 (last access: June 12, 2021), and www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5 (last access: June 12, 2021), respectively.

Author Contributions

JL-M: conceptualization, methodology, validation, investigation, formal analysis, and writing—original draft. RA and CC: methodology, investigation, software, formal analysis, data curation, visualization, and writing—reviewing and editing. NC: methodology, investigation, software, formal analysis, and writing—reviewing and editing. RM: software, formal analysis, and writing—reviewing and editing. SW and JN: investigation, formal analysis, and writing—reviewing and editing. DS, IA, and AM: writing—reviewing and editing. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by projects ANID—FONDECYT: N°11170768 “Potential effects of land use change on fjords of western Patagonia under climate change scenarios” and Interdisciplinary Center for Aquaculture Research (FONDAP INCAR 15110027, ANID). RA was supported by ANID PFCHA/DOCTORADO NACIONAL/2019—21190544. RM and NC participated trough the project C-HydroChange (CGL2017-86788-C3-2-P) of the Spanish Ministry of Economy and Competitiveness.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Directorate of Water (DGA), Meteorological Directorate of Chile (DMC), and the Sub-Secretariat of Water Resources of Argentina (SRHA) for their role in collecting data.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.628454/full#supplementary-material

Footnotes

1. ^ www.clarklabs.org/terrset/idrisi-gis/
2. ^ www.clarklabs.org/terrset/land-change-modeler/
3. ^ www.clarklabs.org

References