Clarifying regional hydrologic controls of the Marañón River, Peru through rapid assessment to inform system-wide basin planning approaches

In response to rising energy and water demands of an increasing population, agriculture, and consumptive industries, rivers are being dammed and diverted in ever more remote and pristine parts of the world. Rivers provide an upstream-downstream continuum that is critical for maintaining healthy environmental, social, and economic conditions along the corridor, and for connecting distinct and diverse ecosystems that are geographically disparate, but which share a common river passage (Vannote et al., 1980). Connectivity across ecosystems is central for a natural system’s ability to respond to change (Groves et al., 2012). Hydropower potentially introduces a twofold disruption on a river system’s resilience by simultaneously imposing major disturbance and physically disconnecting the river continuum.

Conserving the functioning of the river continuum is intrinsic to protecting highly valued environmental and social aspects of the river. Recent calls for a system-wide, holistic approach to river development are gaining momentum to more accurately balance perceived economic benefits of hydropower with impacts to biodiversity, ecosystem services, and local communities (Winemiller et al., 2016). Dams compromise historic biotic and abiotic regimes locally and regionally (Bruno and Siviglia, 2012; Konrad et al., 2012), so why not plan for them over the same scale?

Basin-scale planning, however, requires up-to-date, spatially explicit, accessible data for inventorying environmental and social landscapes (Hartmann et al., 2013). Limited data is a common challenge for scientists and planners studying remote river systems subject to current development pressures. Remote sensing holds promise for data acquisition in remote regions but has its own limitations with cloud cover, thick forests (Chang et al., 1996; Foster et al., 1997), and deep snowpacks (Chang et al., 1987). Filling gaps in our knowledge of core data sets for key environmental characteristics, such as streamflow, sediment transport, and biodiversity, is required to understand the upstream-downstream linkages across regional basins (Groves et al., 2012).

The Marañón River, which drains the eastern flanks of the Peruvian Cordillera Blanca, is one such poorly studied, yet highly important, river corridor under strong development pressure. Since the 1700s the Marañón River, when compared to the longer, but lower discharge, Ucayali River, has been considered a major headwater source of the Amazon River (Lee et al., 2014), and it is one of the last major free-flowing connections from the glaciated crest of the Peruvian Andes to the Amazon jungle. As many as 81 dams have been proposed for the Marañón River and its tributaries (Finer and Jenkins, 2012), with plans for construction of 20 large dams on the mainstem itself.

Understanding the main contributors to river discharge is important to system scale planning because discharge patterns provide the context for ecosystem geography, sediment dispersal, water quality, and human-development patterns along the river. This paper focuses on a dry-season transect to help fill a hydrologic data gap for the upper Marañón River by identifying key hydrologic processes and transitions in the context of hydropower development. We address these goals:

1. Map intra-basin geographical variations in precipitation inputs and cryospheric resources to better understand the role that new water plays in the basin’s water resources.

2. Utilize hydrochemical and isotopic data from samples to describe dominant hydrologic input processes over a 1,745 m elevation gradient and spanning key ecologic transitions from the arid alpine (2,100 m ASL) to the humid jungle (355 m ASL).

3. Clarify the relative source water contributions to the Marañón’s discharge and consider potential changes to the river’s long-term discharge due to climatically sensitive ice and snow resources found in the corridor.

More broadly, to respond to data gaps on similar rivers in remote regions, we present this case study as an example of a combined approach using field observations and remote-sensing image-processing techniques to efficiently characterize mountain-basin hydrologic controls at regional scales.

The Marañón River originates at the glacierized crest of the Cordillera Blanca, an Andean sub-range with the highest concentration of tropical glaciers in the world (Kaser and Osmaston, 2002). The river flows south to north from its mountainous headwaters before it turns eastward to join the Ucayali River to form the upper Amazon River near Iquitos (Figure 1 inset).

Our study entailed sampling at key transitions across the alpine headwaters to the humid, lush, lower-lying jungle (Figures 1, 2b, c). For the purposes of defining and understanding these transitions, we have categorized the study domain into three reaches that flow through three zones: alpine, transitional, and jungle zones (Figure 2). The alpine reach flows through a series of narrow, cascading and highly eroded canyons with an average gradient of 4.1 m/km, whereas the transitional reach has a shallower fall (average gradient 2.2 m/km) through a combination of broader canyons and wide, agricultural valleys. The jungle zone is characterized by lush vegetation, with the river passing through constrictions (pongos) in rugged moderate-elevation mountains, where high discharge and shallow gradient creates massive whirlpools and roller-coaster-like rapids (average gradient 0.94 m/km).

The upper Marañón River corridor runs between two geologically distinct Andean sub-ranges (Figure S1). Eastern (right bank) tributaries drain the Marañón complex, comprising crystalline plutonic and strongly metamorphosed rock. Western (left bank) tributaries have geology that is dominated by a mix of younger marine sedimentary rocks (limestone, sandstone and fluvial conglomerates), younger volcanics, and some ore deposits.

In general, geologic classifications run parallel with the river alignment until the Marañón River’s confluence with the Chamaya River. At this point the channel bend to the northeast cuts across the grain of the Andes. The result of this cross-cutting is that a similar combination of mostly sedimentary geologies occur in the right and left bank subbasins downstream of the Chamaya, notably the high discharge Chinchipe and Utcubamba subbasins (Figure S1).

Vegetation along the study site varies dramatically. The highest elevations are bare rock, ice and snow. The puna ecosystem (local term for the more general ‘páramo’), is sandwiched between snowline and treeline (Young et al., 2017) and is a tropical Andean high-elevation land cover in the headwaters of the Marañón and its tributaries. Puna is generally a “grassland matrix” including tussock, graminoids and forbs. Wetlands, peatlands, fens and bogs are common in alpine regions of the Cordillera Blanca, with ponded water in the wet season and still-saturated slowly draining soils during the dry season (Polk, 2016). Puna is also prized for its water regulation and its role as an efficient water reservoir (Hall et al., 2015; Poulenard et al., 2001) thereby providing important ecosystem services. Alpine ecosystems have been modified over time by both people and climate, with changes of vegetation character and extent in the Cordillera Blanca being driven by agriculture, temperature fluctuations, and reductions of glacial and snow meltwater inputs (Polk, 2016).

Below the puna biome, moderate elevation arid zones (above 2,100 m) are dominated by shrublands (locally called matorrales) with some open montane forests. In contrast, the lusher forest flanking the river as it runs below 400 m is dense tropical forest as seen by the deep greens in satellite images and in Figure 2a.

Glacier surface coverage in the Marañón basin is restricted to those areas that are sufficiently cold year-round to allow for persistent perennial snow fields. Tropical glacier systems are highly susceptible to even small perturbations in temperature and humidity (Kaser, 1999; Wagnon et al., 1999) as well as larger climate systems (Schauwecker et al., 2014). With increasing ambient temperatures recorded from the 1970s to the present, the glacial extent in the Cordillera Blanca has decreased rapidly over the last few decades with glacier termini shifting to higher elevations especially on eastern aspects (Marañón side) of the range (Racoviteanu et al., 2008). Previous Andean studies (Baraer et al., 2012; Bury et al., 2013) show that many basins adjacent to the Marañón headwaters have undergone, or are about to undergo “peak water” (Jansson et al., 2003), a term used to describe the short-term increase in glacial melt-water generation during the final loss of the frozen-water reserves as ice masses wane in a warming climate.

Glacier meltwater is responsible for the majority of streamflow in some Cordillera Blanca headwater catchments (Mark et al., 2005) and can be especially critical for dry-season discharge (Baraer et al., 2012). The percent glacierized area of headwater catchments is significantly related to the specific discharge of that basin (Mark and Mckenzie, 2007), with ice melt’s role diminishing for larger basins where rainfall and snowmelt dominate. This is also demonstrated in the Himalaya (Racoviteanu et al., 2013; Wilson et al., 2016).

We used boats to collect a Lagrangian hydrochemistry data set, by following water downstream as it mixes and evolves (Moody, 1993). We sampled for 30 days, covering 620 kilometers of the upper Marañón River in July and August 2015, during which there were no major mainstem storms that had a noticeable impact on discharge. This data-collection design is ideal for evaluating the environmental response to landscape and discharge as it changes with distance downriver (Konrad et al., 2012). Together with the use of remotely sensed imagery and mapping techniques, we used hydrochemistry and isotopic tracers from water samples to address study goals and to estimate tributary contributions using two-end-member mixing models. Other research activities during this study included sampling for presence and abundance of freshwater invertebrates, beach pebble counts and aerial filming to use for structure-from-motion sediment quantification. This work represents the first baseline environmental data set on the upper Marañón River.

The extensive study domain necessitated prioritization of sampling locations in order to capture the diversity of subbasins within the system in a reasonable time frame and budget. Sampling at larger tributaries was prioritized because these represent an aggregated hydrochemical signature of significant landscapes of the upper Marañón basin. The final sampling plan aimed to achieve a balance between tributaries meeting a basin-area threshold, practicality of sampling, sample count, spatial representation over the 620 km reach, and representativeness of landscape features anticipated to have a hydrochemical impact (e.g., unique geologic deposits, odd uplands, faults, geothermal source waters, etc.). Sampled subbasins make up 60% (38,773 km²) of the entire upper Marañón basin. Attributes of sampled subbasins are summarized in the supplementary material (ST1).

River-water samples (16 mainstem, 16 tributaries) were mostly collected in pairs to allow for rigorous mixing calculations of river discharge (Figure 3, ST2). Samples were collected upstream of the targeted confluences at locations that ensured little to no influence of hyporheic flow. Three additional samples were collected based on observed anomalies in the field: one downstream of the Vijus mine wastewater pipe (km 225), one in a low volume tufa limestone deposit on Rio Magdalena (km 437), and one at Rio Conjun (km 458) to increase representation of right bank subbasins in the dataset.

Care was taken to access well-mixed channels for all sampling. Samples were stored in 60-mL polypropylene Nalgene¹ bottles that were triple rinsed with deionized water, air dried on a sterile paper towel, and capped until used for sampling. A 350-mL sampling container was used for sample collection. Sample bottles and the sampling container were rinsed three times with sample water by vigorously shaking while capped and then emptying a partially filled bottle. The sampling container was then filled mid-channel by plunging the opening upstream at arm’s depth from the boat or from a wading individual, whichever was feasible given the river depth at that point. Ultraviolet treatment (Steripen¹, 1 L sterilization) was used to disinfect the water from live biota. Sample water was filtered using a pressure syringe filtering system and 47-mm diameter Gelman¹ A/E glass fiber filters with an approximate 1-μm pore size. Three times the sample volume (3 × 60-mL) from the sampling container (350-mL) was pushed through the filter prior to filling the sterile 60-mL sample bottles. A positive meniscus was attained prior to capping to eliminate any air headspace in the sample. The sterilized, filtered samples were stored in a dark cooler for long-term transport on boats for the remainder of the expedition.

All water samples were transported to the University of Colorado for analysis. Instruments, methods, precisions, and detection limits for all analytes are provided in the supplementary material (ST3). The University of Colorado’s Geological Sciences Department Laboratory for Environmental and Geological Sciences analyzed samples for dissolved silica, copper, iron and cations (Calcium – Ca²⁺, Magnesium – Mg²⁺, Sodium – Na⁺, Potassium – K⁺) by means of Inductively Coupled Plasma Mass Spectrometer (ICP-MS). The Kiowa wet chemistry lab at the Institute for Arctic and Alpine Research analyzed samples for pH, ANC via manual gran titration, and anions (Chloride – Cl^–, Nitrate – NO₃^–, Sulfate – SO4^2–) via ion chromatography. We measured ionic concentrations, but these ions will be hereafter noted by their elemental or compound abbreviation.

The Kiowa lab also analyzed samples for isotopes of oxygen (δ¹⁸O) and hydrogen (δD) using a Picarro L2130-i ultra high-precision isotopic water analyzer. Stable water isotopes are reported as a δ (per mil, or ‰) ratio of the sample to Vienna Standard Mean Ocean Water (VSMOW), calculated as:

where R is the ratio of stable isotopes of oxygen or hydrogen, O¹⁸/O¹⁶ or H²/H¹ respectively.

Two-end-member mixing models were used to estimate the size of each sampled tributary discharge relative to the Marañón mainstem. This work extends a historical legacy of mixing models, because one of the oldest two-end-member mixing models ever implemented was used by Antonio Raimondi in 1879 (Raimondi, 1879) to estimate the relative discharges of the Marañón and Ucayali Rivers.

In the present study, we use mixing models to estimate water inputs along a river reach. In our mixing models the two ‘end members’ are considered to be tributary discharge and the Marañón River discharge upstream of the tributary confluence with the mainstem (ST2 in supplementary material). Using conservative tracer concentrations, discharge contributions from tributaries relative to the mainstem downstream of the tributary confluence (percent flow) were calculated using equation 2.

Subscripts in equation 2 relate to the sampling configuration in Figure 3, where a refers to the tributary sample, b refers to the mainstem sample upstream of the confluence, and c refers to the mainstem sample downstream of the confluence. The tracer concentration in the sample is referred to as C.

Because our approach uses a targeted but minimal sampling program, we are forced to accept that the chemistry of each mainstem sample (C_c in equation 2) may represent water contributions from other tributaries of a similar composition in the same reach. Because of this assumption the tributaries sampled define the reaches we were able to quantify. The uneven subbasin sampling between left (west) and right (east) bank tributaries was primarily a practical issue due to inaccessibility to right bank streams, and as mentioned above, this defines the reaches we were able to assess. We have a sound understanding of the chemistry of the right bank tributaries from the three right bank samples collected in the Cordillera corridor. We expect a similar chemical composition in un-sampled right bank tributaries waters because all have a similar geology to the rivers that were sampled – dominated by the Maranon Complex (Supplementary figure S1) – and similar topography (Supplementary figure S2).

Discharge contributions were calculated for the suite of conservative tracers (ions and isotopes) analyzed in water samples yielding a range of discharge contribution fractions for each tributary. An error-weighted mean (EWM) of the mixing models for each constituent was determined for each tributary and was used as the modeled discharge contribution for that tributary. Tributary and mainstem discharges were calculated at each subbasin using the discharge record at the Balsas gage on August 3, 2015 (196 m³s^–1), the day we passed the gage. Discharge calculations were propagated upstream and downstream from the Balsas gage using the EWM percent contribution solution for each subbasin.

Error-weighted mean is preferred over other metrics because it considers uncertainty, w, in the mixing model results due to propagated analytical error as calculated by equation 3.

In equation 3, f refers to the fraction of discharge contribution calculated from the mixing model, C is the constituent concentration, w_C is the analytical uncertainty in laboratory analysis specific to the constituent analysis, and w_f is the analytical uncertainty propagated through to the mixing model fraction. This study did not replicate samples or perform duplicate laboratory analysis on samples due to limited sample volumes, therefore w_C is based on the annual laboratory analytical precision (percent relative standard deviation) over all runs in a year (between 4,000–6,000 samples) for a given analyte and a given upstream, downstream or tributary sample, where w_Csample = w_c × Csample or the method detection limit (ST3), whichever is greater. Accordingly, w_C is the same across end members for each constituent and these values are provided in table ST3 in supplementary material.

As is evident from equation 3, mixing model calculations involving two end members of similar concentrations have larger associated uncertainties than models with end members that have disparate concentrations. Error-weighted means (EWM, equation 4) allow for the lower-uncertainty mixing models to dominate the tributary discharge calculation while still considering all model results across analytes.

To understand the most influential constituents in the EWM calculation, we calculated an “influence metric” using equation 5 for each constituent. This provides additional guidance for hydro-chemical interpretation for results.

Total error to the final EWM mixing fraction integrates error across constituents and was calculated as follows (equation 6).

To evaluate mixing model performance for all mixed mainstem samples concentration reconstructions were computed using the EWM solution. Normalized error for each sample was calculated to compare reconstructed concentration with observed concentrations (equation 7) and to examine model bias, if any.

Discharge gages at three locations on the upper Marañón River recorded data over varying historical periods (Figures 1, 5a): 1) Balsas (855 m ASL) from 2015-present, 2) downstream of the Chamaya River confluence at Corral Quemado (CQ, 410 m ASL) from 1977–1981, and 3) downstream of the Chinchipe and Utcubamba confluences at Rentema (355 m ASL) from 1977–1981.

The Balsas gage recorded discharge during the 2015 study period however there are no on-going discharge records downstream of Balsas. To understand the relationship between Balsas discharge and discharge records at CQ and Rentema, we examine multi-year average historical discharges for 1977–1981 at CQ and Rentema as compared to a smoothed historical-average discharge at Balsas. The smoothed Balsas hydrograph is reconstructed from graphs presented in the 2014–2015 Peru Ministry of the Environment National Meteorological and Hydrologic Service (SENAMHI) hydrologic monitoring report (Servicio Nacional de Meteorolgia e Hidrologia del Peru (SENAMHI, 2015).

The 620-km, south-north, linear Marañón River flows through varying hydro-climatology regimes related to elevation gradients and river distance from headwaters. Climate monitoring stations in the Peruvian Andes are sparse, and those stations with precipitation data cover time periods that often do not overlap. To better characterize precipitation distribution over the Marañón basin we used the WorldClim data set (version 1.4) (Hijmans et al., 2005). WorldClim utilizes an interpolated spline algorithm to produce global maps of terrestrial climate from weather station data that adequately represent spatial climate variations, even over remote data-scarce regions. Monthly average precipitation inputs over the length of the dataset (1960–1990) were calculated at 30-arcsecond resolution (~0.93 km at the equator). These gridded monthly means were averaged at the subbasin scale to assess the range of precipitation inputs to major basins across the Upper Marañón basin as well as for the sampled subbasins.

Rain samples were not collected during the study; however, the Global Network of Isotopes in Precipitation (GNIP) (IAEA/WMO, 2016) monitoring network provides two sites that are not within the study domain but are expected to be representative of the alpine and jungle precipitation found in the Marañón.

The closest alpine site is located at 4477 m ASL at Marcapomacocha, 130 km south from the source of the Marañón River, and was used to provide representative rain isotopic values for precipitation inputs to the upper elevations of the Marañón basin (Figure 1A). Similarly, no GNIP record exists for the lower elevations within the study domain, however the Puerto Almendras at Iquitos GNIP precipitation record near the eastern extreme of the Marañón basin (98 m ASL) provides representative isotope values for lower elevation, Amazonian precipitation systems likely to affect the lower subbasins of the Marañón.

The Randolph Glacier Inventory (RGI) (Pfeffer et al., 2014) was used to map and quantify glacier surface areas in the Marañón Basin using glacier area data for the Cordillera Blanca (Cogley et al., 2015; Racoviteanu, 2007) and geographic information system (GIS) software. RGI is a processed subset of the Global Land Ice Measurements from Space (GLIMS) (Raup et al., 2007) glacier monitoring project.

Snow maps were derived from satellite remote sensing. Standard snow maps that exploit band difference ratios of snow surfaces (Dozier, 1989) classify a pixel as either snow covered or snow free. Their performance is degraded in alpine environments such as the Peruvian Andes where mixed land cover pixels are common at NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS) satellite (Justice et al., 1998) resolution (500 m pixel size). Spectral un-mixing algorithms that calculate the snow cover fraction of each pixel (Nolin et al., 1993; Rosenthal and Dozier, 1996) yield snow maps that are more accurate and exhibit less sensitivity to different land cover types (Rittger et al., 2013). Snow cover maps were created using the MODIS Snow Covered Area and Grain Size (MODSCAG) algorithm (Painter et al., 2009) to estimate fractional snow cover. We quantify the probability of snow cover in a given pixel by calculating the number of positive snow-covered observations divided by the total number of clear-sky observations over the 2001–2014 period. This metric estimates the probability that a pixel will have snow on any given day over the 13-year record, indicating the relative snow importance in that pixel. Glacierized areas have a snow probability of 1 because this technique cannot reliably differentiate between snow and ice surfaces. Either seasonal snow cover or ice is present year-round over the glacier footprint.

Separating snow from clouds with MODIS data presents persistent challenges (Rittger et al., 2013), especially in the lower-lying subbasins (jungle zone) where moisture-rich air masses first meet the Andean foothills. Cloud masking algorithms were not able to consistently differentiate between snow and cloud in the jungle zone. To address this shortfall, we utilized the WorldClim temperature dataset (1970–2000) (Hijmans et al., 2005) to map the mean minimum temperature for the coldest month (MMTCM) at approximately 1-km² or 30-arcsecond resolution. We used a 2°C MMTCM upper temperature threshold to indicate where precipitation could plausibly fall as snow. Pixels where MMTCM occasionally fall below the 2°C threshold are not likely to remain persistently cold, resulting in a quick melt response and little influence to snow cover probability.

Our rain-snow temperature threshold is considered conservative as compared to the 1.2°C optimized rain-snow air temperature threshold based on 17.8 million precipitation phase observations (Jennings et al., 2018). Additionally, Jennings and others (2018) find lower elevation and humid locations similar to the lower Maranon basin have cooler thresholds than this optimized temperature, that is, the minimum snowline elevation is probably higher than we identify ensuring we do not erroneously disregard areas where snow may occur.

Using the MMTCM temperature map we mask those areas that are too warm for precipitation to fall as snow, eliminating any snow probability error introduced by clouds in the lower elevation, cloud-prone basins. The final snow probability map is validated with Landsat imagery. The temperature-thresholding data analysis process is demonstrated visually in Supplementary Figure S2.

The puna biome, a non-wetland tropical alpine vegetation ecosystem unique to the Andean cordillera, is mostly characterized by tussock and bunch grasses but also includes shrublands and woodlands (Polk, 2016; Young et al., 2017). We map the puna biome between 3,500 m and 4,800 m elevation (Sanchez-Vega and Dillon, 2006) and the latitudinal range that bounds the ecosystem following Ochoa-Tocachi and others (2016).

Although puna classification can include somewhat steep slopes, the water storage functions of punas are of particular interest to the hydrologic questions addressed here. The gravity-driven drainage of water on steeper slopes limits their ability to contribute to water reservoirs serving dry season baseflow in the Marañón. Accordingly, to more accurately map the puna areas with greater water storage potential, we filter out areas steeper than 20° in our puna biome mapping procedure.

The Marañón basin undergoes pronounced seasonal swings in precipitation, melt patterns and river discharge. During the dry season (June–September) World Clim monthly mean precipitation decreases in all areas, but lower elevation subbasins are noticeably wetter (47–122 mm rain/month) than upper subbasins (10–28 mm rain/month) (Figure 4a, b, c). During the wet season (October-March) the difference in magnitude of precipitation between lower and higher elevation basins is not as large as the dry season. Peak mean monthly precipitation occurs in March (118–178 mm/month) for all subbasins (Figure 4a).

Moisture masses generally move into the Maranon basin via the river corridor from the Amazon basin to the east. This results in the northern and eastern basins bearing the brunt of orographic rain, with an increasingly dry progression from the northeastern major subbasins to those in the south, as visually demonstrated by Figure 4b.

Discharge recorded at gages (Figure 5a) exhibits a ‘flashy’ or rapid response to individual storms in hydrographs across the catchment (Figure 5b). Rentema discharge peaks appear to be strongly affected by localized storms in the low-elevation, high-discharge tributaries, Chinchipe and Utcubamba rivers (blue and cyan in Figure 4b, respectively). The Chamaya basin (orange in Figure 4b) measured at Corral Quemado (CQ), exhibits a subdued flashy behavior compared to Rentema that is consistent with the lower rainfall there. The three large lower basins, Chamaya, Chinchipe and Utcubamba, are hereafter referred to as the lower-lying basins (LLB). Dry season in the LLB is in August and September, with rain events and resultant discharge peaks that are smaller than at other times of the year. The Chinchipe and Utcubamba rivers drain parallel to the Andean ranges in contrast to all the other tributaries which cut across the ranges.

In contrast, high-elevation catchments demonstrate few dry-season discharge peaks due to a general lack of storms over this period (Figure 5b). Despite little precipitation during the June through September dry season, the higher-elevation reaches of the Marañón have an impressive 200 m³s^–1 baseflow.

The smoothed monthly average historical discharge record at Balsas indicates the discharge during the study period was within its normal range and that discharge observed at this location is considerably less than the discharges of CQ and Rentema (Figure 5). The limited discharge records (1978–1982 at CQ and Rentema, record commenced 2015 at Balsas) do not allow for evaluating if the Marañón has experienced “peak water.” However, Baraer et al. (2012) find that melt water-influenced dry season flow likely experienced “peak water” in the 1970s in the neighboring glacierized Santa River basin at La Balsa (approximately 1850 m elevation) downstream of Huaraz.

The fractional contribution of meltwaters from snow and ice at the subbasin scale in the Marañón corridor is controlled by elevation with only the highest altitude topography having cold enough temperatures to receive substantial snow (Figure 6a, Supplementary Figure S2). Similarly, only the three highest elevation subbasins have glacial sources: Yanomayo, Putchka, and Marañón headwaters. The percent (area) of these subbasins that is glacier covered ranges from 0.4% (27.6 km²) in the Marañón headwaters to 4.7% (105.8 km²) in the Yanomayo basin (Figure 6a, Table 1). Total glacier-covered area in the study domain is 168.5 km². Glacier mass balance budgets are dominated by exchanges on the glacier surface (snowfall and melt) (Cuffey and Paterson, 2010) suggesting that the low seasonal snow occurrence in the high-elevation glacierized subbasins is inadequate to reverse glacial ablation trends in the basin.

As defined by the mapping criteria specified in Methods above, the puna biome covers 14.8% of the total study domain (sampled plus unsampled subbasins) (Figure 6b). A majority (55%) of the Headwater subbasin is identified as puna biome representing nearly 3400 km². In contrast to its small glacier area presence (0.43%), this result suggests the Headwater subbasin’s water sample may be an appropriate characterization of the puna water end member.

Due to the elevation control on the puna ecosystem, the higher elevation alpine basins have a larger proportion of puna area, generally around 30%, than subbasins in the transition or jungle zones (Table 1). Notably the lower-lying, large volume Chinchipe, Utcubamba, and Chamaya subbasins have little-to-no puna presence, emphasizing the importance of rainfall generated runoff in contrast to the groundwater reservoirs provided by puna.

The Marcapomacocha GNIP station data demonstrate the temporal variability of precipitation isotopic composition throughout the year and is also used to determine the local meteoric water line (LMWL) for the Marañón headwaters (Figures 7, 8). This station is located at 4,477 m ASL and 130 km south from the headwaters of the Marañón in a similarly arid alpine environment. At the Marcapomacocha site, the wet season (October–May, δO¹⁸ ranging –18.9‰ to –15.2‰) is more isotopically depleted than the dry season (June–September, δO¹⁸ ranging –10.7‰ to –7.1‰) (Figure 7a, left colored bars). The general trends of precipitation isotopes observed at Marcapomacocha agree with the isotopic trends reported for monthly cumulative precipitation over 2006–2007 at 3060 m in Huaraz (Baraer et al., 2015), approximately 80 km west of the Marañón corridor. Huaraz reports more enriched precipitation isotopes than Marcapomacocha, as expected given its lower elevation and the dominant isotope-elevation effect we observed.

No major variation in isotopes is expected across the basin width within the Cordillera corridor (that is, from right ridge to left ridge) given the moisture movement up the river valley from the Amazon. Isotope values from any spillover from the east onto the eastern Cordillera corridor ridge is expected to be dictated by elevation, leading to little variation since left and right bank ridges are of similar altitude.

Precipitation isotopes in the jungle are more enriched than in the alpine as demonstrated by the Puerto Almendras-Iquitos GNIP record at 98 m ASL (Figure 7a, right colored bars). No published low-elevation precipitation isotopes exist for the LLBs preventing a study domain-specific comparison to the Puerto Almendras-Iquitos GNIP record. Prevailing air masses during the June–September dry season come from south easterly trade winds that move δO¹⁸-enriched recycled moisture from the Amazon Basin into the Andes (Windhorst et al., 2013). The dry-season isotopic values in the jungle LLBs are consistent with receiving moisture from the Amazon lowlands to the east. Considering the Amazon origin of moisture masses moving upstream from the east, we consider the Puerto Almendras-Iquitos to be an appropriate indicator of precipitation isotopes for our downstream study extent.

Marañón basin surface-water isotopic values fall within the range of GNIP precipitation at the alpine and jungle site (Figure 7a), with high-elevation mainstem and tributary isotope values more closely related to wet-season alpine precipitation than dry-season precipitation. This indicates the marginal role of dry-season rain or snow in the high-elevation water balance and the likelihood that wet season precipitation is stored in natural reservoirs for dry-season release.

The isotopic character of the subbasins’ surface waters can be grouped into two regimes: the middle and upper subbasins have confining slopes adjacent to the river (hereafter referred to as the Cordillera corridor), and the LLBs are rain-fed jungle tributaries with catchments more directly accessed by Amazon moisture systems. The Cordillera corridor tributaries show a gradual δO¹⁸ isotopic enrichment with decreasing elevation from –14.67‰ to –10.31‰. Enrichment of δO¹⁸ in the mainstem is observed with mixing of progressively enriched tributary inputs from below the confluence of the Putchka River at 1,850 m (–14.25‰) to the final sampling point at Montenegro (–7.41‰). The smaller tributaries in the upper reaches result in a gradual, incremental enrichment of the mainstem Marañón. In contrast, the abrupt isotopic enrichment of Marañón water downstream of the Chamaya confluence at 430 m ASL (Figure 7a) is consistent with large-discharge, Amazonian sourced inputs from the isotopically enriched Chamaya, Chinchipe and presumably Utcubamba tributaries (LLBs).

Downstream enrichment through the Cordillera corridor suggests the elevation effect (Dansgaard, 1964) is a contributing factor to the isotopic gradient of Marañón waters (Figure 7a). Precipitation isotopic signatures affected by altitude are the result of squeezing out of moisture from an air mass as air is adiabatically cooled while rising to greater elevations (Gat, 1996). In precipitation, heavier isotopes (e.g., O¹⁸) are removed preferentially over lighter isotopes (e.g., O¹⁶), leaving behind water vapor that is progressively more depleted in heavier isotopes as it moves up in elevation.

Previous Andean studies have shown spatial variability of isotopes is dominated by the elevation effect, with a range of depletion rates reported from –0.17‰ per 100 m (Garcia et al., 1998) to –0.24‰ per100 m (Gonfiantini et al., 2001). Indeed these rates are similar to the globally established lapse rate of –0.28‰/100 m (Poage and Chamberlain, 2001). In the Cordillera corridor of the Marañón the δO¹⁸ depletion rate with elevation can be described by a linear fit correlating to an isotopic depletion rate of –0.37‰ per 100 m:

Using the above Andean depletion rates, over the 1,520 m elevation change in the Cordillera corridor from the Putchka River (highest tributary at 2,100 m ASL) to the Silaco River (580 m ASL) the elevation change results in a depletion of –2.58 to –3.64‰ δO¹⁸. The actual change observed between these tributaries, –5.87‰ δO¹⁸, suggests that additional factors, probably related to hypsographic effects within individual basins, contribute to the isotopic transition observed over the study domain. In addition to topographic effects, vegetation transpiration, evaporation or sublimation effects on alpine meltwaters and standing water storage in the punas may affect isotopic signatures of surface waters.

The increase in Deuterium excess (D-excess) in precipitation (indicated as the y-intercept of δO¹⁸-δDeuterium regressions in Figure 8) of the LMWLs relative to the Global Meteoric Water Line (GMWL) reflects moisture recycling in Amazonia as systems travel across several thousand kilometers of humid, forested regions (Friedman, 1977; Salati et al., 1979). The Amazonian moisture system influence on precipitation D-excess is mirrored in surface waters receiving these rain inputs, as demonstrated by the rapid D-excess increase in LLB tributaries (Figure 7b). Transpiration, not post-deposition snowpack processes, as the driver of high D-excess is logical given the far greater magnitude of water recirculation through the Amazon moisture systems as compared to the snow volumes mapped across the study domain (Figure 6).

The δO¹⁸-δDeuterium bi-plot (Figure 8) also sheds light on the impact of evaporation or sublimation fluxes to the watershed because oxygen isotopes fractionate at a higher rate during evaporation and sublimation than hydrogen isotopes thereby changing the slope of the relation. The influence of evaporation and sublimation on water surfaces in the tropical Andes appears to be highly site-specific, with some studies reporting a high evaporation signal (Williams et al., 2001) while others do not (Baraer et al., 2015; Stichler et al., 2001). Some evaporative effect is observed in the alpine samples (m = 7.39) relative to the LMWL (m = 8.2), suggesting possible importance to surface water of snow, ice, or ponded (puna) water in the alpine reach, whereas in the alpine zone effects from transpiration are expected to be minimal due to limited vegetation.

Ion concentrations in the mainstem Marañón generally increase as the river progresses downstream through the Cordillera corridor (Figure 9) but are highly variable across subbasins in the study domain indicating somewhat diverse geology of the area (Figure S1). Tributaries Yanomayo, Shuve, Conjun, Chamaya and Chinchipe are lower in Ca than the mainstem, while all other tributaries are more concentrated, some by nearly a factor of two (e.g., Rupac = 3,206 µeq/L Ca; Cresnejas = 3,371 µeq/L Ca; Magdalena = 3,289 µeq/L Ca). A marked dilution of the mainstem Marañón occurs with the inputs of the Chamaya and Chinchipe/Utcubamba tributaries (ST2), demonstrating lower solute concentrations combined with considerable volume as compared to the mainstem.

Alkalinity, SO₄, and Ca dominate the constituents of Marañón waters (Figure 9, supplementary Figure S3). The right bank tributaries generally exhibit lower concentrations of SO₄ and Ca than left bank tributaries (Figure 9). High levels of SO₄ observed in the left bank tributaries are attributable to the presence of minor evaporates in the Peruvian Andes marine sediments such as the Late Jurassic marine deposits prominent in the left bank tributaries (Stallard and Edmond, 1983) (Figure S1). Jungle tributaries have a higher dominance of carbonate hardness than alpine or transitional tributaries that are balanced by other anions other than the carbonates (Piper, 1944). Highly enriched stream waters are observed in the tufa stream at Rio Magdalena, which is also a distinct outlier in the Piper diagram (supplementary Figure S3).

PCA consolidates multi-variable geochemical data into a data representation, capturing most of the variation of the original data with fewer degrees of freedom (the principal components) (Christophersen and Hooper, 1992). PCA can be used to identify hydrochemical controls of the system (Figure 10). Two principal components explain 74.2% of the variability in the data with the third PC adding only 8.7% additional explanation. PC1 explains 47.2% of the variability and is associated with chemical weathering products and an additional influence from evaporite deposits (Stallard and Edmond, 1983). The PC2 explains 31.0% of variability and is associated with the isotopic variables and silica. Water samples plotted on the PCA bi-plot show that controls of hydrochemical signatures vary over the different zones (Figures 2c, 10). The 95% PCA confidence ellipses for the mean of each zone show that the alpine and jungle zones are well discriminated.

As is evident, the alpine zone is highly dependent on PC1 which is related to geochemical weathering products of the catchment, whereas it is rather independent of PC2 that includes isotopic variables. This is not surprising given the disproportionately high contribution of suspended sediments and geochemical constituents contributed to Amazon waters by the high relief areas of the Amazon basin (Gibbs, 1967; Stallard and Edmond, 1983). Geochemical weathering controls in the alpine may also indicate the importance of groundwater for sustaining dry-season baseflow. The alpine’s independence to isotopic variables suggests a uniform, consistent source water throughout the zone. In contrast, the confidence ellipse for the jungle samples incorporating the LLBs plots along the 45-degree axis indicating this zone is affected by both components and thus weathering products, and a defined isotopic source-water signature (associated with the influence of Amazon moisture systems as discussed above), are both influential. The transition zone lies in between.

Using a two end-member mixing model, relative tributary volumes are calculated using a suite of conservative tracers. A core assumption in this approach is that the chemistry of each tributary may represent inputs from additional tributaries of a similar composition in the same reach.

The error-weighted means (EWM) of mixing model results for each tributary are shown in Figure 11 and Table 2, with less than 2% total mixing fraction error for all tributaries (0.84% mean error). Of the Cordillera corridor tributaries, the Putchka, Yanomayo and San Miguel rivers contribute the most discharge to the Marañón relative to mainstem discharge at their confluence (35%, 22%, 20%, respectively). This is somewhat expected given that the Putchka and Yanomayo subbasins make up a relatively large proportion of total basin area given their location high in the catchment. Also, the Putchka and Yanomayo subbasins are glacierized whereas the San Miguel receives year-round water from Laguna Piaz, a substantial lake 3km upstream of the Marañón confluence. The natural reservoirs (glacier ice and lake) in these subbasins provide dry-season flow inputs absent in other tributaries.

An especially significant mixing fraction result is that the Chamaya and Chinchipe/Uctubamba rivers are major dry season inputs to the mainstem Marañón, with the discharge of each tributary nearly equaling or exceeding the discharge of the Marañón mainstem, upstream of the confluence (Figure 11). The Utcubamba River, flowing in from the right bank 2 km upstream of the Chinchipe confluence, was not sampled for logistical reasons. The discharge contribution calculated between the upstream-of-Chinchipe and downstream-of-Chinchipe sampling points more than doubles the mainstem flow and includes the volume contribution of the Utcubamba on the assumption that the geology (refer site description above and Supplementary figure S1), topography (Supplementary figure S2, right panel), vegetation, and geographic placement within the basin is similar to the Chinchipe. A rainfall runoff regime similar to the Chinchipe is expected in the Utcubamba subbasin given precipitation mapping results (Figure 4). Refer hatched pattern subbasin shown in Figures 1, 2c and 5a for Utcubamba basin outline.

Discharge volumes for tributaries and the mainstem at each sampled subbasin (Table 2) demonstrate the similar order of magnitude of inputs provided from the Cordillera corridor subbasins in contrast to the LLBs. Mainstem discharge reconstructions at the historical gage locations of Corral Quemado (516 m³s^–1) and Rentema (1546 m³s^–1) are similar to the historical average flows at those locations in July which are estimated to be 600 m³s^–1 and 1100 m³s^–1, respectively (Figure 5). The higher flow calculated at the Rentema gage site may be due to the timing of Chinchipe sampling, which occurred shortly after a precipitation event in the Chinchipe headwaters. This would result in a temporary pulse in Chinchipe discharge and a higher overall mainstem discharge at the Rentema gage site as compared to the historical July average. Notably, the mainstem discharge increases by a factor of 31 across the domain, demonstrating the scale of hydrologic variation observed in this study.

Constituents that most influenced the EWM result varied depending on the tributary. SO₄, Ca, Mg and Cl were common influential ions among most subbasins as calculated by the influence metric (Table 2a). Isotopic variables, δO¹⁸ and δD, were important to the EWM calculation outside of the four high-elevation subbasins.

Mixing-model performance was evaluated by way of comparing the mainstem concentrations using the EWM mixing fraction result with the observed concentrations in water samples (Figure 12a). The modeled versus observed comparison demonstrates very good model performance (m = 1.02, b = 4.96, R² = 0.99). The sample below the Cajas River confluence is an obvious Acid Neutralizing Capacity (ANC) outlier as marked in Figure 12a. In this instance the mixing model result is affected by an unsolvable 2-component mixing model problem wherein both the Cajas tributary and upstream Cajas ANC observations (the two end members) are smaller than the downstream Cajas sample (mixed sample result). This physically infeasible mixing situation suggests an additional but unmeasured chemical input to the Marañón between the mainstem samples collected upstream and downstream of the Cajas confluence. The EWM approach somewhat ameliorates the effect of this outlier on the final mixing fraction result by reducing the influence of highly uncertain mixing models like this one on the EWM.

Normalized errors for all analytes are shown in Figure 12b. The median error is near zero for all analytes, with a tight interquartile spread for all but Cl^– and ANC which are highlighted as red and green points, respectively, in Figure 12a. All other constituent values are shown as black points.

The mismatch between modeled and observed concentrations for the ‘Below Cajas’ ANC sample increases the magnitude of normalized error for ANC. Consideration of the overall magnitude of Cl^– concentrations is important for interpreting chlorine’s larger normalized error values. Because Cl^– concentrations are low (i.e., generally on the order of 100–200 µeq/L), a small difference between modeled and observed values results in a high normalized error as compared to more concentrated constituents with concentrations >1,000 µeq/L.

Of note, the modeled isotope (δO¹⁸, δD) performance is especially good despite their low magnitude values making them prone to being influenced by small analytic laboratory error. This strong result indicates that they are reliable tracers for use in mixing models.

As mentioned above, this limited data mixing model approach acknowledges that the chemistry of each mainstem sample may represent inputs from the sampled tributary as well as additional tributaries of a similar composition in the same reach. Errors arise when non-sampled inputs of markedly different composition enter the mainstem between the tributary confluence and the downstream sample point assumed to represent a fully mixed mainstem sample (Figure 3). An example of this effect is discussed in relation the Cajas River mixing model solution in the section above. Drawbacks of using the minimal-sampling approach for efficient characterization are also demonstrated in the spread in mixing fraction results across analytes for a given tributary (Table 3b). This problem, which we deem minor, would have been further reduced if we had had the capacity to take additional samples once the river became well mixed downstream closer to tributary confluences.

Unmeasured inputs with unusual chemical concentrations, for example from point sources or geologically anomalous subbasins, can sometimes be anticipated. For example, in the case of San Miguel River, several notable inputs occur between the upstream (end member) and downstream (fully mixed) sample including wastewater byproducts associated with the town of Chagual, riverside Vijus mine’s discharge, and major in-channel earthworks associated with a bridge construction site. More constrained mixing model results with lower variability across constituents (Table 3b, e.g., Chamaya, Chinchipe tributaries) have few potentially anomalous inputs, with upstream and downstream sampling sites located much closer together.

The EWM approach to mixing model fractions constrains the spread observed in mixing model results, collapsing the suite of results into an average value influenced proportionally by constituent-specific uncertainties. This approach provides a metric for identifying influential constituents in the EWM calculation, those with lower associated uncertainties. As is evident by equation 3, uncertainty values increase rapidly as the difference in concentrations, ΔC, between end members (tributary and upstream sample, in this case) approach 0. Conversely, end members with high ΔC have smaller uncertainties. The constituents that are found to be heavily influential (low uncertainty) in this study (Table 3a, e.g., SO₄), are those with concentrations that oscillate from high to low, providing continually large ΔC as the river progresses down river. This suggests that the presence of SO₄ is highly variable across the domain, whereas Na or K are more constrained, consistently present at high concentrations yielding higher uncertainty in mixing model calculations and thus less influence to the EWM.

The suite of constituents that influence the EWM changes across the domain. Isotopic variables are non-influential to EWM in the higher basins when isotopic levels across subbasins are comparatively homogenous. In contrast, isotopes consistently influence EWM at mid- and lower-basins when isotopic composition changes rapidly with elevation (refer slope changes in Figure 7a). In summary, the ‘influence metric’ highlights constituent variability within the system, and high constituent variability generally relates to lower 2-component mixing model fraction uncertainty yielding greater influence on the EWM. Table 3 demonstrates this result visually. The greyed values in Table 3a, b and c are those with less than 5% influence on the error weighted mean calculation. These low-influence constituents (Table 3a) in most cases correlate to the outlier mixing model results (Table 3b) as is expected since constituents with high uncertainty (low influence) produce mixing fractions with higher deviation. Tables 3a and 3b combine in Table 3c to determine each constituent’s contribution to the EWM. The small total error (Table 2) is a result of having relatively few dominant constituents influencing the EWM.

We had initially assumed, based on naming and length expectations, that the largest water contribution would be from the highest elevation Marañón headwaters and its tributaries, that is, the river has a mainstem dominance. Discharge in the Marañón River is not driven by the headwater hydrologic system. Instead, our mixing model results show that much of the discharge is derived from the low lying Chamaya, Chinchipe and Utcubamba basins (Figure 11). Our finding that the Marañón mainstem dry-season discharge is decidedly dominated by the high-runoff LLB subbasins (Figure 11) contradicts a general expectation that rivers are named for the tributary that controls annual discharge. Moreover, this finding shows another example that river length may not be a critical variable from which to anticipate relative discharge of continental-scale river systems.

At these subbasin confluences, Marañón discharges are controlled by high-precipitation inputs to the LLBs with more direct access to Amazon moisture systems. Moisture recycling (re-evaporation cycles) also leads to increases in D-excess as observed in the Chamaya (D-excess = 14.34) and Chinchipe/Utcubamba (D-excess = 16.05) tributaries, and this is reiterated in our PCA results showing moisture recycling variables having strong influence on jungle water hydrochemistry.

Given the Amazon moisture-movement patterns, the LLBs are the first subbasins of the study domain to intercept the wet air masses from the east. As these masses hit the topography of the Andes, considerable new precipitation falls in the LLBs even during the dry winter season (June–September), as evidenced by greener vegetation (Figure 2a) and year-round flashy hydrographs (Figure 5b).

In contrast, subbasins in the Cordillera corridor are flanked on both banks by high peaks and lie in a rain shadow imposed from both the east and west. The bulk of the precipitation in the Cordillera corridor comes in March and April, synchronized with the end of the Amazonian rainy season when occasional Antarctic storms reach the Amazon (Friedman, 1977). These fronts can be especially moisture rich and able to overcome the initial orographic rain-out effect. Additionally, a subtle east-west snow cover trend is observed in the Cordillera corridor with snow cover probabilities higher on the eastern ridges (Figure 6a), suggesting some Amazonian moisture spillover to the eastern flanks of the Marañón corridor during the colder winter months.

The consistent 200 m³s^–1 dry season baseflow observed at Balsas, however, implicates a reliable and massive reservoir to support this substantial flow over the winter period when the Cordillera corridor receives few, if any, meteoric inputs (Figure 4). The water storage of approximately 1.5 km³ is required to sustain this baseflow for the three-month dry season, equating to approximately a three-year water supply for Lima’s 10 million person population.

Precipitation falling as snow plays a minor role in the water balance throughout the system because those areas of the Marañón basin that are cold enough for snow lack major winter precipitation inputs. Only the Yanomayo basin has an appreciable glacierized area (5%), and the low snow cover to the upper basins implies low probability for glacier regeneration, supporting observations of rapidly wasting glacier systems across the Cordillera Blanca region (Racoviteanu et al., 2008).

The majority of precipitation to the high-elevation catchments comes during the warmer wet season months of January through March when most precipitation falls as rain, as demonstrated by the flashy wet season hydrograph at Balsas (Figure 5b). The isotope values for the alpine baseflow are bracketed by fall and spring alpine precipitation isotopes, and they are dissimilar to winter precipitation isotopes. After infiltration to groundwater, fractionation largely stops thereby preserving the isotopic signature of summer rain. These data suggests dry season baseflow waters are likely sourced from stored wet season rain.

Without the natural storage and delayed release functions of snowpacks and glaciers, the question revolves around other environmental systems with the capacity to store the summer season rains for delayed winter release. Our puna ecosystem mapping (Figure 6) shows the Cordillera corridor subbasins have 29–55% of their basin area identified to be shallow gradient (<20°) puna biome. These environments are well documented to hold vast water storage reservoirs that in turn provide ecosystem services throughout the Andes. The PCA analysis (Figure 10) supports the routing of alpine surface waters through the subsurface, with hydrochemical variables being the primary influence in differentiating alpine water samples but consistent isotope values in alpine waters implying the source of baseflow is regenerated during a confined season.

Discharge amount and geochemical data (Figure 10) indicate that the substantial baseflow that provides dry season water to the Cordillera corridor – the site of many of the proposed hydropower projects – during the dry winter season is sourced from puna reservoirs above 3500 m ASL. Punas and puna-like ecosystems called páramo modulate the release of water received during extreme weather patterns (prolonged wet periods) through significant soil storage capacity, releasing water gradually during dry seasons to provide sustained baseflow to rivers (Balslev and Luteyn, 1992; Buytaert et al., 2007). Puna water storage functions are well known to provide large volume water storage as ecosystem services, including serving the majority of municipal water supplies to major centers including Quito and Bogata (Borrelli et al., 2015; Buytaert et al., 2006). Indeed, basins with water reservoirs have been shown to produce substantially more runoff than basins without these ecosystems at least partially due to low levels of evapotranspiration and the high storage capacity afforded by highly organic soils (Buytaert et al., 2007, 2006). These distinct natural wetland reservoirs likely serve as critical wet season rain storage role for the Marañón. Ponded wet season rain present in puna ecosystems also more likely explains the increased fractionation signal observed in baseflow.

A key question for this research revolved around the role that snow and ice play in the hydrology of the Marañón because of their high sensitivity to a warming climate (Barnett et al., 2005). This in turn raises questions about long-term discharge sustainability, an important consideration for viable hydropower development and basin-wide planning. The direct role of snow and ice melt as a dry-season river source (i.e., meltwater’s direct translation to open channel flow) appears to be limited due to (a) little snow and limited glacier coverage in the headwaters (Figure 6a) and (b) rain’s control on flow in the LLBs (Figure 4).

However, groundwater-surface water interactions in this region are complex and shifts in these interactions are possible due to secondary effects from glacial recession (Gordon et al., 2015). Upward vegetation expansion into melt-out areas (Young et al., 2017) may increase evapo-transpiration and decrease available alpine groundwater reserves, although this effect may not be noticeable given its magnitude relative to the apparently massive groundwater volumes present in the Marañón basin. Increasing fragmentation and changes to the spatial distribution of water-controlled ecosystems such as punas and wetlands are thought to be caused by compromises to subsurface connectivity due to lessening meltwater inputs (Polk, 2016; Young, 2015). In this way glacier recession and resulting changes to meltwater generation may indirectly affect baseflow by way of impacting groundwater-surface water interactions and subsurface water storage reservoirs in the Marañón, eventually affecting the character of baseflow.

In summary, water resource sensitivity to climate is complex in the Marañón basin. The natural, gravity-fed subsurface puna reservoirs appear to provide ample storage of wet season rain to service substantial dry season baseflow through the Cordillera corridor. The puna’s hydrologic regulation function is compromised by anthropogenic disturbance (Ochoa-Tocachi et al., 2016), and conservation strategies are needed to protect these central ecosystem services. Long-term climate effects on baseflow may be more far-reaching and convoluted, and require further study. In contrast, in the high-runoff LLBs, long-term discharge variation due to climate changes would likely be caused by changes in precipitation patterns (timing and amount).

This case study demonstrates the need for a multipronged approach to the characterization of remote rivers such as the Marañón. Baseline information is critical to provide a quantitative and objective evaluation mechanism for development projects. Yet, mobilizing and conducting studies over large regions and in compressed timeframes is challenging. Processing of remotely sensed imagery and mapping can yield targeted, efficient field studies that can be used to represent hydrologic settings in ways that are not obvious in the field or from a literature review. In a sense, chemical analyses of stream-water samples effectively represent remote sensing of this upstream watershed and provides information not available through satellite imagery alone.

Dry-season discharge volume in the Marañón River, before its turn east to enter the Amazon lowlands, is dominated by low-elevation tributaries receiving large precipitation inputs from Amazonian moisture masses, and is not strongly influenced by discharge volume in the mainstem upper Marañón. In the upper basin the substantial dry season baseflow in the Cordillera corridor appears to rely on the slow release of wet season rain from massive natural storage reservoirs present in the higher elevation puna biome. Together, these high- and low-elevation subbasins define the dry-season response that the Marañón River will have to major hydropower development in terms of the long-term reliability of water supply.

Numerous remote mountain river systems, especially those in the Andes and across the Hindu Kush Himalaya, are poorly understood but are targeted for large-scale development as global energy demands rise. Mountain rivers serve as the connecting conduit across elevations and ecosystems that maintain the continuum of material and nutrient transport, with far-reaching environmental, social and economic consequences. There is an urgent need for the environmental science community to use creative methods to efficiently characterize large-scale systems under development pressure. Understanding hydrologic controls and baseline conditions allow evaluation of the upstream-downstream effects of development decisions before actions are taken, possibly resulting in irreversible changes to some of the most important and bio-diverse ecosystems on the planet.

Water sample hydro chemistry and isotope data unique to this study is uploaded as online supporting information. The following data used here is available publicly online:

• Web reference for public download:http://worldclim.org/version2

• Timespan of data used: gridded average monthly from 1960–1990 (full record)

• Web reference for public download:https://snow.jpl.nasa.gov/portal/

• Timespan of data used: 2001–2014

• Web reference for public download:http://nsidc.org/data/nsidc-0272

• Timespan of data used: latest GLIMS version updated 2012

• Sites/regions: refer data submitter references below for site reference identification numbers

• Web reference for public download:http://www-naweb.iaea.org/napc/ih/IHS_resources_isohis.html

• Timespan of data used: monthly data over 2006–2012

• Sites: Pto Almendres-Iquitos (WMO code 8437701) and Marcapomacocha (WMO code 8457001)

The success of the Marañón Project is due in no small part to the 16-person team of river runners who paddled 620km of the Marañón River and assisted with the study presented here. Especially sincere thanks to James “Rocky” Contos, Tomas Binimelis, Pedro Peña, Ben Griffin, Jen-ai Stokesbary, Nate Mack, Nancy Moore and Christian Martin for tireless coordination, risk management and logistics to get this month-long remote river trip on the water and down the river safely. Thank you also to Henry Worobec for his recording efforts, and to the suite of scientists who shared the motivation to initiate this baseline data set: Jaime Ruth Goode, Natalie Kramer Anderson, Jorge Peralta, Stephanie Gustavson, Diana Silva and Florencia Trama. Thank you to Lily Jackson for navigating geologic data sources, guidance on geologic mapping and grouping of rock types. Much gratitude goes to Holly Miller and her staff in the INSTAAR wet chemistry lab and Fred Luiszer at the University of Colorado Geological Sciences Department’s Laboratory for Environmental and Geological Sciences (LEGS) for the high-quality lab analysis of all water samples. The final version of this manuscript was greatly improved thanks to Deborah Martin’s in-house U.S. Geological Survey professional review and thoughtful and thorough peer reviews provided by Katalyn Voss and Oliver Chadwick.

The Maranon Project was partially funded by National Geographic Committee for Science and Exploration grant # 9752-15 awarded to Natalie Kramer Anderson (PI), and by generous friends and family who supported The Marañón Project crowd-funding campaign. Robert F. Stallard has been funded by the National Research Program for Water at the U.S. Geological Survey. Publication of this article was funded by the University of Colorado Boulder Libraries Open Access Fund.

The authors have no competing interests to declare.

• Contributed to conception and design: AH, RS

• Contributed to acquisition of data: AH, RS, KR

• Contributed to analysis and interpretation of data: AH, RS, KR

• Drafted and revised the article: AH, RS, KR

• Approved the submitted version for publication: AH, RS, KR