Annual River Runoff Variations and Trends for the Andes Cordillera

Weanalyzedmodeled river runoff variationswest of theAndesCordillera’s continental divide for 1979/80–2013/14 (35 years). Our foci were annual runoff conditions, runoff origins (rain, snowmelt, and glacier ice), and runoff spatiotemporal variability. Low and high runoff conditions were defined as occurrences that fall outside the 10th (low values) and 90th (high values) percentile values of the period of record. SnowModel and HydroFlow modeling tools were used at 4-km horizontal grid increments and 3-h time intervals. NASA Modern-Era Retrospective Analysis for Research and Applications (MERRA) datasets were used as atmospheric forcing. This modeling system includes evaporation and sublimation from snow-covered surfaces, but it does not take into account evapotranspiration from bare and vegetation-covered soils and from river and lake surfaces. In general for theAndesCordillera, the simulated runoff decreased before 1997 and increased afterward. This could be due to a model precipitation artifact in the MERRA forcing. If so, this artifact would influence the number of years with low runoff values, which decreased over time, while the number of high runoff values increased over time. For latitudes south of ;408S, both the greatest decrease in the number of low runoff values and the greatest increase in high runoff values occurred. High runoff values averaged 84% and 58% higher than low values for nonglacierized and glacierized catchments, respectively. Furthermore, for glacierized catchments, 61%and 62% of the runoff originated from rain-derived runoff during low and high runoff extreme years, respectively; 28% and 30% from snowmelt-derived runoff; and 11% and 8% from glacier-ice-melt-derived runoff. As the results could be MERRA dependent, more work with other precipitation forcings and/or in situ measurements is needed to assess whether these are real runoff behaviors or artifacts.


Introduction
Analyses of observed and modeled weather and hydrological data indicate that extreme events will increase Denotes content that is immediately available upon publication as open access. in frequency, intensity, spatial extent, duration, and timing with a warming climate (e.g., Thodsen 2007;Coumou and Rahmstorf 2012;IPCC 2012, p. 7). Extreme events are defined as ''the occurrence of a value of a variable above (or below) a threshold value near the upper (or lower) ends of the range of observed values of the variable'' (IPCC 2012, p. 5). Gleason et al. (2008) and Lubchenco and Karl (2012) made this definition quantitative, specifying extreme conditions to be defined as occurrences that fall outside the 10th and the 90th percentile values for the record's period.
In the central Chilean Andes (the region between 308 and 388S), hydrological (Favier et al. 2009) and climatological (e.g., Masiokas et al. 2006) studies identified knowledge gaps for nival and glaciological mass balance observations and processes (Rabatel et al. 2011). Snow and glacier monitoring, in combination with processfocused analyses and modeling, are crucial for understanding the hydrometeorological roles of snow and glaciers. Mernild et al. (2017b) revealed that annually ;85% of the Rio Olivares basin (338S) freshwater runoff originated from snowmelt , making snowmelt the dominant contributor of the catchment's runoff. Thus, runoff is dominated by variability in snow precipitation rather than by variability in mean annual air temperature (MAAT), emphasizing a nival runoff regime. Valdés-Pineda et al. (2016) analyzed runoff data series for 107 stations in central-southern Chile and concluded that precipitation was the most dominant factor affecting runoff availability. Peña and Nazarala (1987) simulated runoff from the upper Rio Maipo catchment (338S) and stated that even though runoff from glaciers is small in relation to the total runoff (henceforth, ''total runoff'' is defined as runoff from the catchment outlet), their contribution is significant during dry years (low precipitation years). This is especially the case at the end of the summer season in February and March (austral summer). In their landmark study, Peña and Nazarala (1987) simulated mean monthly glacier-derived runoff for a 4-yr period (1981/82-1982/83 and 1984/85-1985/86) at the hydrometric station El Manzano located near Santiago de Chile (348S). Overall, glacier runoff contributed 2% 6 2% (October; where 6 is the standard deviation), 5% 6 4% (November), 11% 6 8% (December), 16% 6 11% (January), 19% 6 13% (February), and 19% 6 11% (March) of the annual runoff. At the end-of-summer season, the February glacier-derived runoff, as a percentage of total runoff, ranged between 34% (1982, extremely wet year, strong El Niño event) and 67% (1969, extremely dry year, 95% exceedance). Casassa et al. (2015) simulated glacial runoff based on a temperature and shortwave radiation model validated with field-based energy and mass balance data for debris-free, debris-covered, and rock glaciers at Yerba Loca valley (33 8S) and an updated glacier inventory for the Rio Maipo basin. The estimated glacierderived runoff for the Rio Maipo basin at the hydrometric station El Manzano (338S) constituted, on average, 13% of the annual runoff for the year 2014/15 (extremely dry year, similar to 1968/69) and 62% of the runoff in March 2015, although the relevance of glacier runoff strongly decreased downstream (Casassa et al. 2015).
In a warming climate, shrinking glacier area may be crucial for water resources. Paul and Mölg (2014) and Malmros et al. (2016) estimated changes in glacier area over time for glaciers located in the Palena district of Chile (408-448S) and the central Chilean and Argentinean Andes (328-338S), respectively. They showed a mean glacier area shrinkage of 25% (1985-2011) and 30% (1955-2013/14). In the Palena district, 374 glaciers out of 1664 glaciers disappeared (Paul and Mölg 2014). Also, in the semiarid central Chilean and Argentinean Andes, Lliboutry (1998), Leiva (1999), and Pellicciotti et al. (2007Pellicciotti et al. ( , 2014 showed glacier area shrinkage during recent decades. For instance, Georges (2004) estimated glacier area shrinkage to be approximately 20%-30% in the tropical Cordillera Blanca (98S) of Peru; here the glacier area changed from 800-850 km 2 in 1930 to 620 km 2 in 1990. Other glacier studies from catchments in Ecuador (18S; Pouget et al. 2015), Peru (98-128S; Mark et al. 2010;Drenkhan et al. 2015), and Bolivia (168S; Soruco et al. 2015) have shown that glacier area shrinkage produced a reduction in river runoff. For the Patagonia Ice Fields, the maximum potential contribution to sea level rise is 14.7 6 2.9 mm sea level equivalent (Carrivick et al. 2016). Recent studies on the Northern and Southern Patagonian Ice Fields have highlighted post-Little Ice Age glacier area shrinkage on the order of 11%-14% for the period ;1870-2011 (Davies and Glasser 2012;Falaschi et al. 2013;White and Copland 2015). If glaciers along the Andes Cordillera continue to retreat in area, annual glacier runoff will decrease as reductions in glacier size outweigh the effects from melting glaciers (Casassa et al. 2009;Sharp et al. 2011).
In this study, SnowModel and HydroFlow Mernild and Liston 2012) were used to simulate hydrological conditions along the river systems and individual catchments west of the continental divide of the Andes Cordillera. The modeling system does take into account evaporation and sublimation from areas covered by snow and ice, but the model is limited by not including routines for evapotranspiration from river and lake surfaces and from bare and vegetation-covered soils. It is well known that the runoff-rainfall ratio is below 1 in Andes catchments (García et al. 2011), although catchments dominated by snowfall may have relatively high ratios (Berghuijs et al. 2014). Therefore, the simulated runoff may be overestimated compared to real conditions. These simulations focused on the spatiotemporal freshwater river runoff variabilities, including both low and high annual river runoff values from nonglacierized and glacierized catchments. Also, the ratios between total runoff and rainderived runoff, snowmelt-derived runoff, and glacier-icemelt-derived runoff were simulated both for low and high runoff years and for dry years. This is to highlight the importance of snow-derived and glacier-ice-derived runoff during particularly low annual runoff events.
The atmospheric forcings were provided by the NASA MERRA products on an hourly time step and a 2 /38 longitude 3 1 /28 latitude grid. The NASA MERRA atmospheric data were downscaled using meteorological algorithms within MicroMet (Liston and Elder 2006a) to 4-km and 3-h intervals. Observed snow, glacier mass balance, and river runoff time series partly covering the period 1979-2014 for specific nonglacierized and glacierized catchments were used for verifications of SnowModel and HydroFlow simulations. Those tests indicated statistical agreement among snow, glacier mass balance, and river runoff simulations and independent field and remote sensing observations (for further details, see Mernild et al. 2017a,b,c,d).
We simulated and analyzed the annual low and high river runoff values in nonglacierized and glacierized catchments (n 5 4224 corresponds to the number of simulated catchments) along the Andes Cordillera west of the continental divide ( Fig. 1), ranging from Tierra del Fuego (Patagonian subpolar latitudes) to the mountain tropical Andes climate. Only annual values were studied; daily and monthly low and high river runoff values were not analyzed. The analysis covered 1979-2014. This study focused on the following research questions: 1) Does the overall occurrence of low and high river runoff values follow the overall trend in annual runoff? 2) Has the number of annual runoff values (low plus high runoff values) increased since 1979 in connection with the annual runoff trend? 3) What is the spatiotemporal distribution of low and high annual river runoff values along the Andes Cordillera? 4) Is there a link among low and high annual river runoff values and annual rain-derived runoff, snowmeltderived runoff, and glacier-ice-melt-derived runoff? 5) Is runoff from glacierized catchments during dry years mainly from snowmelt-derived runoff or glacier-ice-melt-derived runoff?
This analysis provides us with an understanding of runoff lows and highs, during a period where the future MAAT (toward 2100) is expected to be many degrees Celsius above present-day surface air temperature conditions (Collins et al. 2013).

a. SnowModel and HydroFlow
SnowModel (Liston and Elder 2006b) was used to simulate meteorological, snow, and glacier mass balance conditions for the Andes Cordillera (Mernild et al. 2017a,b,c). HydroFlow was used to simulate individual catchment runoff routing and river runoff to surrounding fjords and seas (Mernild et al. 2017d; Fig. 1), including low and high annual runoff values (presented in this paper). SnowModel is a surface model simulating distributed meteorological, energy balance, snow, ice surface, and freshwater runoff conditions by using several submodels. In this study, five of the six SnowModel submodels were used (DataAssim, a data assimilation submodel, was not applied here): MicroMet is a grid-based meteorological model (Liston and Elder 2006a), Enbal is a surface energy transfer model (Liston 1995;Liston et al. 1999), SnowTran-3D is a blowing snow model Sturm 1998, 2002;Liston et al. 2007), SnowPack-ML is a multiple-layer snowpack evolution model , and HydroFlow is a linear reservoir, gridded hydrologic routing model .
For SnowModel, river runoff was simulated to originate from rain, snowmelt, and/or ice melt, including that from bare glacier surfaces. When a surface is defined as glacier ice, meltwater is assumed to be transferred instantaneously to a linear runoff reservoir. As part of the snow-cover simulation, SnowPack-ML calculates refreezing and internal retention processes and snowpack percolation of snow surface meltwater and rain. These internal snowpack procedures have an influence on the meltwater flow and delay, including how long it takes for the freshwater runoff to reach the coastline . Accounting for these internal snowpack processes in SnowPack-ML improves realism in simulated meltwater transport through the snowpack and across the landscape, yielding improved runoff event timing and values.
In HydroFlow, adjacent cells are linked via a topographically controlled flow network within each individual catchment (Fig. 2). Each grid cell acts as a linear reservoir. Here, water is transferred from a specific grid and from all upslope grids to the connecting downslope grids through an eight compass direction network: north, northeast, east, etc. (Fig. 2). In each grid cell, there are two flow transfer responses: slow and fast . Slow response includes rain, snowmelt, and glacier ice melt within the snow and ice matrices (in the case of glaciers and snow) and within the subsurface (in the case of snow-covered and snow-free land) to the river network. Fast response comprises channelized flow in glaciers and rivers. Liston and Mernild (2012) presented HydroFlow model tests for a specific glacierized and nonglacierized catchment in southeast Greenland. They found that r 2 (square of the linear correlation coefficient) values for simulated runoff compared to observed runoff ranged from 0.63 to 0.77 and that Nash-Sutcliffe coefficient (NSC) values (Nash and Sutcliffe 1970) ranged from 0.60 to 0.61 among individual years.
b. Meteorological forcing, model configuration, and simulation domain The spatiotemporal freshwater river runoff hydrographs, including low and high annual runoff values, were simulated from 1 April 1979 through 31 March 2014 (following the hydrological year from 1 April through 31 March). The NASA MERRA data were used as forcing for the simulations (Bosilovich 2008;Bosilovich et al. 2008Cullather and Bosilovich 2011;Rienecker et al. 2011;Robertson et al. 2011) and were thereafter aggregated to 3-h values, which was the lowest temporal resolution possible to still resolve the diurnal cycle [see Mernild et al. (2014) for further information regarding the NASA MERRA forcing]. MERRA is an atmospheric reanalysis dataset, taking advantage of modern satellite-era datasets FIG. 1. Western part of South America: (a) topography (m; color increment is not linear) and locations of MERRA atmospheric forcing grid points used in the model simulations (black dots); (b) simulated individual drainage basins (n 5 4224; represented by multiple colors) west of the continental divide draining to the Pacific Ocean. The borders are highlighted by straight lines between the countries Colombia (Co), Ecuador (E), Peru (P), and Chile (Ch); and (c) 12 glacierized catchments (red dots with bold black line) where detailed river runoff analyses were conducted (see Fig. 11). (Bosilovich 2008;Bosilovich et al. 2008Liston and Hiemstra 2011). In 1998, MERRA began assimilating the Advanced Microwave Sounding Unit (AMSU) radiance, improving the representation of water cycle processes .
MicroMet downscaled MERRA forcings simulating spatiotemporal fields of atmospheric conditions, following well-known temperature-elevation, wind topography, humidity-cloudiness, and radiation-cloud-topography relationships (Liston and Elder 2006a). In addition, the MERRA fields were spatially distributed using the Barnes objective analysis scheme (Barnes 1964(Barnes , 1973. It is therefore expected that the downscaled trends generated by SnowModel/MicroMet are not significantly different from the trends present in the MERRA dataset. Furthermore, Mernild et al. (2017a) show that spatiotemporal air temperature and precipitation trends are in agreement with observations for the Andes Cordillera.
Topography data from the U.S. Geological Survey's 7.5-arc-s Global Multiresolution Terrain Elevation Data 2010 (GMTED2010) were used (Danielson and Gesch 2011). The topography data were rescaled to a 4-km horizontal grid increment (Figs. 1a, 2). The land-cover file was a hybrid dataset. For nonglacier areas and for glacier areas, the land-cover file was created from a 300-m-resolution 2009 Global Land Cover Map (GlobCover 2.3; http://due.esrin.esa.int/page_globcover.php) and from the Randolph Glacier Inventory version 4.0 (Pfeffer et al. 2014; Fig. 2), respectively. SnowModel does not take into account dynamic glacier changes and subsequent changes in area, hypsometry, and thinning. SnowModel uses a time-invariant glacier area and hypsometry defined from datasets produced during the 2000s and 2010, respectively (Mernild et al. 2017c). The use of time-invariant conditions may result in uncertainties in the simulated glacier surface mass balance FIG. 2. Examples of flow networks calculated from gridded topography and ocean-mask datasets (the gridcell size is 4 km 3 4 km, 16 km 2 ) to illustrate the HydroFlow network configuration over the simulation domain for the individual catchments 1 (11 232 km 2 ), 6 (7984 km 2 ), 7 (15 712 km 2 ), and 9 (3056 km 2 ) (the catchment locations are illustrated in Fig. 1c, and specifications are shown in Table 5). The flow network is marked with black lines, the topographic divide with white lines, the catchment with a transparent gray color, and the ocean with white. Glaciers are represented by light-gray squares and topography (m MSL) by multiple colors. The numbers on the x and y axes are the gridcell numbers, with x 5 1, y 5 1 in the lower left corner of the simulation domain. Each grid cell is 4 km 3 4 km in size.

JULY 2018
M E R N I L D E T A L .
(SMB) budget. This is particularly evident at the beginning of the simulation period, where, for example, simulated runoff based on time-invariant conditions may be underestimated relative to changeable glacier conditions due to changes in glacier area. The use of time-invariant conditions likely underestimates dryseason river runoff at the beginning of the simulation period due to underestimated glacier area. The simulation domain was divided into individual drainage catchments by HydroFlow (Figs. 1b, 2). In total, 4224 catchments were generated for the domain west of the Continental Divide (1 240 432 km 2 ; henceforth, the simulation domain will be defined as the area west of the Continental Divide). Further, the catchments were separated on a country scale (only west of the Continental Divide) belonging to either Colombia (121 952 km 2 ), Ecuador (108 448 km 2 ), Peru (302 384 km 2 ), and/or Chile (707 648 km 2 ), based on administrative country borders and nearby watershed divides, where the watershed divides were not necessarily collocated with the administrative country borders. For Colombia, 214 individual catchments were estimated equal to ;5% of the domain catchments, having a mean catchment area of 570 6 3360 km 2 . For Ecuador, it was 236, ;5%, and 460 6 2770 km 2 ; for Peru, it was 488, ;12%, and 620 6 2030 km 2 ; and for Chile, it was 3286, ;78%, and 215 6 2035 km 2 . These individual catchments are defined as catchments where freshwater runoff is routed through the HydroFlowgenerated routing network and into the surrounding fjords or directly into the Pacific Ocean (Fig. 1).

c. SnowModel and HydroFlow verification
This river runoff study used SnowModel forced with MERRA data to produce simulations for the Andes Cordillera to study climate, snow, glacier mass balance, and hydrological conditions. The SnowModel MERRA simulations were verified against coincident independent field observational datasets, as described by Mernild et al. (2017a,b). Overall, simulated maximum annual snowcover extent during winters 2000-14 agreed with MODIS observations (r 2 5 0.67); 2010-14 snow depth agreement was poorer including scaling of point versus grid cells (r 2 5 0.26, RMSE 5 0.18 m); 1 September 2010-14 snow density errors were 10 kg m 23 (corresponding to a difference of 3% relative to the mean snow density); and the snow line location for 10 January 2001 and 4 May 2001 showed differences of 2.1 6 1.4 km and 1.6 6 1.3 km, respectively. Modeled snow depletion curves for 2000/01-2013/14 correlated with remote sensing observations (r 2 5 0.98) and simulated river runoff with gauged river freshwater runoff during 1979/80-1991/92 (r 2 5 0.51). As described, verification occurred using independent direct field observations (snow depths, snow density, and river freshwater runoff) and remote sensing data (MODIS, MOD10A1 snow-cover extent, snow line, and depletion curves) for specific catchments and regions along the Andes Cordillera. These verifications yielded, on average, acceptable results (no significant differences between simulated values and independent observations). Importantly, scale mismatch between snow point observations and 4-km-resolution model simulations always present obstacles to meaningful comparisons without understanding snow measurement distributions and subgrid variation, especially in rugged remote terrain.
In Mernild et al. (2017d), SnowModel-and HydroFlowsimulated flow network and daily river discharge (hydrograph) time series data were verified against satellite images showing the location of the river network (and the location of the hydrographic station) and independent observed daily river discharge time series (n 5 16, 1979-2014; covering latitudes 258-408S), respectively. On average, the modeled runoff showed an r 2 value of 0.51. For the 16 individual hydrographic stations with the longest observed runoff time series (spanning $98% of the 1979-2014 simulation range), the r 2 values were as large as 0.69. Hydrographic stations used were likely not influenced by upstream anthropogenic influences such as dams (analyzed from Google Earth and maps), providing confidence in our discharge assessment. Similar assessments and outcomes for simulated runoff were described by Beamer et al. (2016) when using SnowModel and HydroFlow for southeastern parts of Alaska, and by Liston and Mernild (2012) for southeast Greenland. Overall, significant verifications were conducted for these parameters: snowcover extent, snow depth, snow density, location of snow line, snow depletion curve, glacier mass balance, and river runoff, yielding agreements between simulations and independent observations (for further detailed information regarding the quantitative validation metrics, see Mernild et al. 2017a,b,c). However, the runoff partitioning between runoff originating from rain, snowmelt, and ice melt was not verified due to insufficient hydrochemistry and water isotope data.

Low and high runoff metrics
For estimating low and high runoff values, this definition was used: the occurrence of a value of a variable above (or below) a threshold value near the upper (or lower) ends of the range of observed values of the variable (IPCC 2012, p. 5). In this analysis of low and high freshwater river runoff conditions and occurrences, we define low and high values to be annual runoff totals that fall outside the 10th and 90th percentile values, respectively, over the period of record. For each individual catchment, the threshold values for the 10th and 90th percentiles were calculated and considered in the analyses. Similar 10th and 90th percentile procedures are well accepted and widely used in atmospheric and climate sciences; for example, these thresholds have been applied in efforts to understand the physical links between climate extremes and climate warming (e.g., Karl et al. 1996;Coumou and Rahmstorf 2012), and the same definition is used here to estimate low and high runoff values.

a. Temporal variability in low and high runoff values
On a country scale for 1979/80-2013/14, annual simulated river runoff appeared to increase (all linear trends were significant) in Colombia, Ecuador, and Chile and decrease in Peru (not significant) [see Fig. 3 and Table 1 for slopes s, r 2 , and p values (p is level of significance; based on linear trends)]. For the entire simulation domain along the Andes Cordillera, annual runoff increased significantly and was especially noticeable after 1997/98. Prior to 1996/97 (for the entire simulation domain), annual runoff appeared to slightly decrease (insignificantly), having a mean value of 1259.3 6 50.3 3 10 9 m 3 yr 21 (equal to 1259.3 6 50.3 km 3 yr 21 ). After 1997, a slight increase (insignificant) having a mean value of 1477.3 6 51.1 3 10 9 m 3 yr 21 was apparent (Fig. 3). This increasing runoff trend for the entire domain is mostly attributable to the increase in annual runoff for Ecuador from 160.9 6 21.5 3 10 9 m 3 yr 21 in the first half to 262.6 6 24.2 3 10 9 m 3 yr 21 in the second half of the time series (Fig. 3). For Colombia the difference in runoff values before and after 1996/97 was less pronounced compared with Ecuador. For Colombia it was 373.4 6 39.6 3 10 9 m 3 yr 21 and 447.4 6 44.8 3 10 9 m 3 yr 21 , respectively. In the second half of the time series, annual simulated runoff increased significantly for Colombia; otherwise, trends were insignificant both in the first half and the second half of the time series (Table 1). Further, in Table 1, the 10th and 90th percentile values are shown indicating different values for the entire period compared to both the first and the second half of the time series. For Chile the runoff in the first half was 492.3 6 23.3 3 10 9 m 3 yr 21 , whereas it changed to 534.6 6 44.9 3 10 9 m 3 yr 21 in the second half of the time series (Fig. 3).
This apparent runoff trend for Ecuador between 58S and 58N (only nonglacierized basins occurred) requires cautious interpretation because of the absence of gauge data. The precipitation increase (significant) over time (Mernild et al. 2017a, their Figs. 4e,f;Mernild et al. 2017d, their Fig. 6) may be a model artifact due, in part, to changes in SnowModel and HydroFlow forcing from MERRA. The addition of AMSU radiance M. Bosilovich 2017, personal communication) in 1998 could explain the precipitation and discharge rise. The artifact appears over this region and is not likely over other regions; in situ runoff observations were available for verification of simulated runoff for other regions. To fully understand this and/or to document potential forcing artifacts over a specific region would require extended analyses with other datasets that are outside the scope of this study.
On the catchment scale, declining individual simulated annual runoff conditions have been reported (Mernild et al. 2017d); this is likely due to shrinking glacier areas outweighing the effect of increased glacierice-derived runoff (Casassa et al. 2009;Sharp et al. 2011), and/or changes in precipitation patterns associated with changes and variations in large-scale atmospheric and oceanic circulation patterns (e.g., Garreaud et al. 2009;Núñez et al. 2013Núñez et al. , 2014. Throughout the simulation period, the occurrence of annual low runoff values (values lower than the 10th percentile) and annual high runoff values (values greater than the 90th percentile) varied through time. For example, 1989/90 was a year with low runoff values for Ecuador, Peru, and Chile, and 2002/03 and 2008/09 were years with high runoff values for Ecuador and Chile, and Colombia and Ecuador, respectively (Table 2). On average, for Colombia, Ecuador, Chile, and the entire domain, the presence of annual low runoff values occurred in the first half of the simulation period (before 1996/97), whereas annual high runoff values occurred in the last half of the period (after 1997/98; Fig. 3, Table 2). For Peru, the occurrence of annual runoff values was different, with both low and high values not clearly divided into specific time periods (Fig. 3, Table 2). Such a distribution in time between low and high runoff values could potentially be due to a model artifact in the MERRA forcing (see the previous discussion and section 4e). To fully understand this and to document potential MERRA forcing artifacts over the Andes Cordillera would require extended analyses with other datasets and/or in situ data to assess if this is a real runoff behavior or artifacts.
Such patterns of low and high rank values for the entire domain and for the individual regions could appear suspicious in cases where low or high runoff values vary by a small percentage (e.g., ,10%) from the mean. However, this percentile threshold approach is well accepted and widely used for time series analyses in atmospheric and climate sciences (e.g., Karl et al. 1996;Coumou and Rahmstorf 2012), just as 80% thresholds are commonly used to assess the presence of drought (e.g., Van Loon and Van Lanen 2012).

b. Latitudinal variability in low and high runoff values
In Fig. 5, the latitudinal distribution for each degree latitude in simulated mean pentadal catchment runoff (1979/80-1983/84, 1984/85-1988/89, etc.) are illustrated (in percentage values), showing both the 90th and 10th percentiles. The 1979/80-1983/84 runoff values outside the 10th and 90th percentile indicate that the greatest clusters of low runoff values occurred in the regions between 58N and 58S and between 438 and 508S. Here, ;30% and ;35% of the low runoff values occurred (Fig. 5a), respectively. The greatest cluster of high runoff values occurred between 18 and 108S, where ;55% of the high runoff values occurred (Fig. 5a). These clusters of both low and high runoff values vary in latitude and in range between the different pentads (Figs. 5a-g), indicating, at a first glance, no systematic pentadal variability in low and high runoff values since 1979/80. In the first and second halves of the time series, values outside the 10th and 90th percentiles indicated that the greatest clusters of low runoff values occurred south of 408S. Regarding high runoff values, no clear geographical distribution occurred for 1980-97, where for 1998-2014 a cluster occurred south of 408S (Figs. 5i,j). However, a systematic pattern occurs (from 1979/80-1983/84 through 2009/10-2013/14) if the pentadal change (linear) in the number of runoff values is analyzed for specific latitude bands (Fig. 5h). For latitudes south of ;408S, the greatest decrease in the number of low runoff values, together with the greatest increase in the number of high runoff values occurred; this follows the overall pattern for the Andes Cordillera toward fewer low-and increased high-runoff values (Fig. 5h). For the latitudes between ;408S and 88N, a reduced change in both runoff low and high values occurred compared to south of ;408S (Fig. 5h). With respect to low runoff values in Colombia, ;40% of the catchments had increasing trends and ;60% had decreasing trends, while those numbers were 20% and 80% in Ecuador, 30% and 70% in Peru, and 0% and 100% in Chile, respectively (Table 3, Fig. 5h). For high runoff values, the values were ;50% and 50%, 0% and 100%, 50% and 50%, and 75% and 25% for Colombia, Ecuador, Peru, and Chile, respectively (Table 3), indicating a variability in low and high runoff trends among the countries.
In comparison, for all simulated basins, the 35-yr mean low and high runoff sums were 1013 3 10 9 m 3 yr 21 (19.6 L s 21 km 22 ) and 1780 3 10 9 m 3 yr 21 (31.5 L s 21 km 22 ), respectively (high runoff values averaged 76% higher than low values; Table 4). The lowest low and high runoff values along the Andes Cordillera occurred from the arid Atacama Desert area and southward to the city of La Serena (188-308S), Chile (Fig. 6a). A similar trend occurred for specific runoff low values. Low and high annual mean runoff values for nonglacierized catchments  increased to the north and south of the Atacama Desert; maximum low and high runoff values in the north, influenced by the rain, and in the south of the domain, influenced by snowmelt, were observed (Mernild et al. 2017c). The greatest differences, for example, in specific runoff values on the individual catchment scale (the difference between high runoff values and low runoff values shown as a percentage difference) occurred south of ;408S, between 108 and 208S, and north of 58N. Here, the differences were greater than ;60% and were likely related to seasonal changes in hydrometeorological conditions. The minimum difference , 10% occurred between the arid Atacama Desert area and southward to La Serena and between 58 and 78S between the cities Piura and Chiclayo in northern Peru (Fig. 6a). These two arid regions have an annual mean precipitation  that spans from almost zero to ;0.25 m w.e. (see Fig. 4d in Mernild et al. 2017a).
In Fig. 7a, the difference between low and high runoff is presented for nonglacierized catchments along the Andes Cordillera. The high runoff was 1304 3 10 9 m 3 yr 21 (scaled to 100%), where the low runoff was ;709 3 10 9 m 3 yr 21 (54%). This indicates that the high runoff was approximately twice the volume compared to the low runoff. Similar runoff conditions occurred for both rain-derived runoff and snowmelt-derived runoff (Fig. 7a). On average, for the period 1979-2014, during high and low runoff years, 95% of the runoff was rain-derived runoff and 5% was snowmelt-derived runoff (Fig. 7b). These values cover variability in runoff along the Andes Cordillera and between low and high runoff values. First, during low runoff events almost all of the catchments located in the northern part of the Andes Cordillera were rain-derived runoff catchments, whereas in the southern part (south of 408S) up to one-third of the runoff was rain-derived runoff and two-thirds was snowmelt-derived runoff. Second, during high runoff events the origin of runoff in the northern part of the Andes Cordillera was insignificantly different from the conditions during low runoff events, whereas in the southern part around one-fourth of the runoff was rain-derived runoff and three-fourths was snowmelt-derived runoff. Overall, for nonglacierized catchments along the Andes Cordillera, snowmelt-derived runoff for both low and high runoff was less important compared to rain-derived runoff (Fig. 7b), indicating that runoff conditions are strongly influenced by rain-derived runoff regimes.

d. Low and high runoff from glacierized catchments
Characteristically, glacierized catchments exhibited runoff (low and high), including relatively high specific runoff values (Fig. 6b), compared with nonglacierized catchments. The specific runoffs from glacierized catchments for low and high runoff events were 63.3 L s 21 km 22 (runoff sum 304 3 10 9 m 3 yr 21 ) and 86.6 L s 21 km 22 (runoff sum 479 3 10 9 m 3 yr 21 ), respectively (high sum runoff values averaged 58% higher than low values; Table 4). In comparison to nonglacierized catchments, the low and high runoff from glacierized catchments were 0.18 3 10 9 and 0.32 3 10 9 m 3 yr 21 , respectively (high values averaged 84% higher than low values; Table 4). A similar overall trend in high and low runoff values was present if only the region south of 408S were analyzed, where the annual mean glacierized catchment runoff was greater (not shown) than the runoff from nonglacierized catchments. The difference between high and low values among glacierized catchments (58%) and nonglacierized catchments (84%) indicate that glacier ice melt provides a baseflow component that reduces the difference between the values. From a spatial perspective along the Andes Cordillera, the greatest percentage difference between annual high and low specific runoff values occurred south of ;408S and between 98 and 178S (Figs. 7a,b). For both of these regions, the difference was .;70%. The difference dropped from .;70% (;408S) to ;20% (308S; Fig. 6b), following the pattern for nonglacierized catchments (Fig. 6a).
For glacierized catchments along the Andes Cordillera, years with both low and high runoff values and dry years have important consequences for water resources (where a dry year is defined as a ''low precipitation year,'' where the annual precipitation falls outside the 10th percentile value). Therefore, there is a need to understand the variability due to rain, snowmelt, and glacier ice runoff components (Figs. 9a-c). Linear regressions between rain-derived runoff and runoff, and snowmelt-derived runoff and runoff, for the 10th and 90th percentiles revealed r 2 values of 0.89 and 0.87, and 0.97 and 0.97, respectively. This indicates that runoff is appropriate for estimating the individual rain-and snowmelt-derived runoff contributions for both low and high values. The glacier-ice-melt-derived runoff and FIG. 6. The latitudinal variability in simulated annual runoff and specific runoff outside the 10th and 90th percentiles from 1979/80 through 2013/14 for each individual catchment, including the percentage difference in specific runoff between the 10th and the 90th percentiles (highlighted in yellow) for (a) nonglacierized catchments (n 5 4037) and (b) glacierized catchments (n 5 187). The x axis is logarithmic. its link to runoff are complex (Fig. 8c); catchmentscale glacier-ice-melt-derived runoff is strongly related to glacier area (r 2 5 0.75 for high runoff values and r 2 5 0.79 for low runoff values; Fig. 8d) and catchment size. Figure 8c shows that annual runoff values , 7.5 3 10 9 m 3 yr 21 (illustrated with a vertical dotted line) are correlated with glacierice-melt-derived runoff for both high (r 2 5 0.60; significant) and low runoff events (r 2 5 0.65; significant; not shown). Annual runoff from glacierized catchments . 7.5 3 10 9 m 3 yr 21 is insignificantly correlated with glacier-ice-melt-derived runoff for both the 10th and 90th percentiles (Fig. 8c); this is likely because runoff from these catchments is influenced not only by glacier-ice-melt-derived runoff, but also by other hydrometeorological processes such as snowmelt and rain. For glacierized catchments along the Andes Cordillera, the difference between low and high runoff is revealed in Fig. 9. The high runoff  was 479 3 10 9 m 3 yr 21 (scaled to 100%), and the low runoff was 304 3 10 9 m 3 yr 21 (63%); therefore, the annual high runoff was less than twice the volume of the low runoff. Similar runoff conditions occurred for rain-derived runoff and snowmeltderived runoff (Fig. 7a), where high runoff was less than twice the volume of the low runoff. Regarding the glacierice-melt-derived runoff, high runoff values were scaled to 100% and low runoff values to 82% (Fig. 9a). Overall, these numbers suggest that high runoff was less than twice the volume of the low runoff. On average, 61% and 62% of the annual runoff originated from rain-derived runoff during low and high runoff years. The relative runoff proportions were 28% and 30% from snowmelt and 11% and 8% from glacier ice melt (Fig. 9b). These values cover variabilities in runoff components between the individual catchments along the Andes Cordillera (Fig. 10). For individual catchments south of 408S, the glacier-ice-melt-derived runoff was up to ;90% of the annual runoff. For individual catchments north of 408S, glacier-ice-meltderived runoff was up to ;10% of the runoff during low runoff events, and ;5% during high runoff events (1979Fig. 10). Casassa et al. (2015) analyzed runoff conditions for a catchment in central Chile for the year 2014/ 15; they showed that glacier-ice-melt-derived runoff was 13% of annual runoff. For dry years, SnowModel-MERRA simulations showed the domain ice-meltderived runoff from glacierized catchments to be 35.5 3 10 9 m 3 yr 21 or ;10% of the runoff from glacierized catchments (not shown); the maximum value for a specific catchment (538S) was 90%. In dry years, the snowmelt-derived runoff from glacierized catchments was 95.8 3 10 9 m 3 yr 21 or ;30% of the overall runoff (not shown).  Runoff sum (m 3 yr 21 ) 10th percentile value 1013 3 10 9 709 3 10 9 304 3 10 9 90th percentile value 1780 3 10 9 (76%) 1304 3 10 9 (84%) 479 3 10 9 (58%) To address the issues of spatiotemporal variability, 12 glacierized catchments along the Andes Cordillera were selected for further analysis (Figs. 1 and 11); these catchments were typically separated by distances between ;200 and 500 km. Low and high annual runoff values occurred throughout the 12 individual catchments, indicating no systematic pattern among the catchments (Fig. 11). For 11 of these 12 catchments, the runoff time series indicated insignificant annual differences in the percentage of rain, snowmelt, and glacier ice runoff (Fig. 11a). The exception was catchment 4 (located in northern Chile), where runoff origin shifts occurred. For catchment 4, the variability in annual runoff ( Fig. 11a) could be due to either domination by a rainderived runoff regime or by a snowmelt-derived runoff regime (Fig. 11b). Overall, runoff seems to be largely driven by rain or snowmelt. A systematic change along the Andes Cordillera occurs, where northern catchments were mainly influenced by rain-derived runoff and snowmelt-derived runoff during low and high annual runoff events, whereas southern catchments mainly were influenced by snowmelt-derived runoff and glacier-ice-melt-derived runoff (Fig. 11c, Table 5). As exemplified for catchments 11 and 12, during annual low and high runoff events the glacier-ice-melt-derived runoff varied between 20% and 40% and between 60% and 80%, respectively, with mean values of 67% (meltwater from both snowmelt and ice melt accounted for 86%) and 75% (88%), respectively (Fig. 11, Table 5).

e. Model limitations and runoff perspectives
SnowModel simulations running on 3-h time steps, and HydroFlow river runoff simulations running on daily time steps, have been used to present a detailed and realistic representation of the spatiotemporal distribution of annual low and high runoff events along the Andes Cordillera. The SnowModel and HydroFlow simulations were based on fully integrated models, including energy balance, blowing snow, multilayer snowpack, and runoff routing submodels, in contrast to other glacier mass balance and runoff models that have unrealistically relied on air temperature as a proxy for the energy available for melt (e.g., Rivera 2004). Accounting for basic physical processes such as fully integrated energy balances and multilayer snow evolution within the snowpack are required components of any credible model used to simulate hydrometeorological processes and conditions, particularly under changing environmental conditions. Also, because incoming solar radiation is the primary source of energy over virtually all snow and glacier ice surfaces, it must be accounted for when calculating snowmelt and icemelt processes; the melt contribution from incoming radiation is an order of magnitude greater than melt from air-temperature-related sensible heat fluxes (e.g., Liston and Hiemstra 2011) in relatively dry environments. Subdiurnal model time steps are also required to account for the subdaily variability in solar radiation and its associated energy-related processes and fluxes.
In Chile, more than 550 individual stream gauges (178-548S) are operated by the Dirección General de Aguas (DGA; http://dgasatel.mop.cl). The majority of these individual stream gauges are located in river catchments where anthropogenic impacts such as damming and irrigation are common. Since these anthropogenic contributions influence the natural variability in the runoff hydrograph, physical modeling of the natural variability during low and high runoff events becomes an important way to assess changes over time, estimate and identify underlying runoff variability, and examine water runoff sources (e.g., Liston and Mernild 2012;Bliss et al. 2014).
SnowModel and HydroFlow have limitations. SnowModel assumes, for example, one-way atmospheric forcing from the downscaled reanalysis MERRA atmospheric data, where the atmospheric forcing is prescribed at each single time step without regard for whether the snow and ice distribution and properties might be different than those in the original MERRA reanalysis. Further, the MERRA forcing after 1998 is influenced by the AMSU radiance , M. Bosilovich 2017. In general, the simulated runoff increased over time and specifically in Ecuador due to an increase in the precipitation forcing (significant; Mernild et al. 2017a, their Figs. 4e,f;Mernild et al. 2017d, their Fig. 6), which may be due to a model artifact in the MERRA forcing. The runoff trend for Ecuador (only nonglacierized basins occurred) requires cautious interpretation in light of absent gauge data. To fully understand this and/or to document potential forcing artifacts over a specific region such as Ecuador would require extended analyses with other datasets and/or in situ data to assess, if these are real runoff behaviors or artifacts. Such an analysis is outside the scope of this study, but it would be highly relevant as a follow-up study to elucidate the necessity to obtain consistent reanalysis precipitation datasets for the Andes Cordillera.
In the natural system, the atmospheric variables would be modified in response to differences and changes in surface conditions (Liston and Hiemstra 2011). HydroFlow assumes that each grid has two flow transfer responses: a slow and a fast response . This is an improvement over other models that assume glacier mass loss and ice-melt-derived runoff contribute instantaneously to sea level rise or hydrographic conditions in fjords and/or the ocean. On the other hand, HydroFlow does not represent processes related to soil moisture and evaporation, and evaporation from the river surface and therefore runoff are not influenced by these processes. Further, HydroFlow does not represent processes related to temporary storage (delay) of water flow in and out of lakes downstream of the flow network. SnowModel and HydroFlow do not take into account the occurrence of jökulhlaups [in the literature also referred to as a sudden glacier lake outburst flood (GLOF); e.g., Dussaillant et al. 2010;Carrivick and Tweed 2016;Emmer et al. 2016]. SnowModel and HydroFlow also neglect glacier mass loss from dynamic activities, melting from internal glacier ice deformation, melting from changes in the internal drainage system, subglacial geothermal melting, and subglacial frictional melting due to basal ice motion. Overall, with their relatively high temporal (3 h) and spatial resolution (4 km), and its physically based representations of the governing processes, the SnowModel and HydroFlow simulations presented herein provide an improvement over previous river runoff model studies in assessing the low and high values of Andes Cordillera river runoff.
This SnowModel and HydroFlow model study and its associated analyses have not addressed how river runoff and low and high runoff events may change in the future. While in a future warming climate, shrinking glacier area may be crucial for river runoff, changes in glacier  Fig. 1c). Also shown are the 90th and 10th mean catchment runoff percentile values (horizontal dotted lines). Note the different scales on the ordinate. (b) The origin and variability of annual runoff from rain, snowmelt, and ice melt and (c) annual mean runoff contributions from rain, snowmelt, and ice melt outside the 10th and 90th percentiles. area have not been considered in this study (SnowModel neglects these changes throughout the simulation period due to the use of a time-invariant glacier area and hypsometry). Therefore, to fully understand the role of glaciers in future Andes runoff processes, magnitudes, and low and high events, the time evolution of the glacier area must be taken into account. Only in this way can potential runoff tipping points be identified, since ultimately the annual glacier runoff amount will decline as reductions in glacier area outweigh the effect of glacier melting (Casassa et al. 2009;Sharp et al. 2011).

Conclusions
This analysis of SnowModel/HydroFlow low and high river runoff provides the following conclusions: 1) For Colombia, Ecuador, and Chile, the presence of annual low runoff events occurred in the first half of the simulation period (before 1996/97), whereas annual high runoff events occurred in the last half of the period (after 1997/98). However, for Peru the occurrence of annual runoff events was different, with both low and high events not clearly divided into specific time periods. 2) Since 1979, the number of annual low runoff values declined (linear; insignificant) for the Andes Cordillera, whereas the number of annual high runoff values increased (linear; significant). For example, the highest number of high runoff events occurred for the pentad 2009-14. On average, the pentadal sums of annual runoff events-low plus high events-increased over time (linear; significant). Also, a relationship between the percentage of catchments with low (10%) and high (90%) runoff events and mean annual catchment runoff is apparent. This highlights how decreasing the catchment runoff can increase the percentage of catchments with low runoff events and how an increase in the catchment runoff increases the percentage of catchments with high runoff events. 3) For Ecuador between 58S and 58N, the increase in simulated runoff was influenced by an increase in simulated precipitation (significant) over time. However, one could argue that a model artifact influences the level of change in Ecuadorian runoff after 1998 (after MERRA began assimilating the AMSU radiance) as long as in situ runoff observations are not available for validation. To fully assess this situation and conclusions 1 and 2, and to document potential MERRA artifacts and their resultant impacts on this work, other precipitation reanalysis datasets and in situ measurements are required to assess whether these runoff behaviors are real. 4) Both low and high annual runoff events varied with respect to latitude (here illustrated at 18 intervals along the Andes Cordillera) and between the different pentads throughout the simulation period. However, the greatest decrease in the number of annual low runoff events and the greatest increase in the number of annual high runoff events occurred south of ;408S. 5) For nonglacierized catchments along the Andes Cordillera, 95% of the runoff originated from rainderived runoff, while 5% originated from snowmeltderived runoff during both high and low runoff events.
In the northern part of the domain, both low and high runoff events were dominated by rain-derived runoff.
In the southern part, during low runoff events, onethird of the runoff was rain-derived runoff and twothirds was snowmelt-derived runoff, and during high TABLE 5. Grid, catchment area, and glacier cover for the 12 catchments together with the origin of runoff from rain, snowmelt, and ice melt for both low and high runoff events. The locations of all 12 catchments are shown in Fig. 1c, where four catchments (marked with *) are illustrated in more detail in Fig. 2. The value in parentheses in the last two columns indicates sum of the runoff from snowmelt and glacier ice (%).
Annual mean runoff from rain, snowmelt, and glacier ice melt (%) Catchment No. and river name Grid, catchment area, and glacier cover (km 2 and %) 10th percentile 90th percentile  (88) runoff events, one-fourth of the runoff was rain-derived runoff and three-fourths was snowmeltderived runoff. Overall, for nonglacierized catchments along the Andes Cordillera, low and high annual runoff was highly influenced by rain-derived runoff regimes. 6) For glacierized catchments, during both annual high and low runoff events, the runoff was proportioned as follows: 61% and 62% from rain-derived runoff, 28% and 30% from snowmelt-derived runoff, and 11% and 8% from glacier-ice-melt-derived runoff, respectively. For individual catchments south of 408S, the ice-meltderived runoff contributed as much as ;90% of the runoff, and north of 408S the ice-melt-derived runoff was up to ;10% of the runoff during low runoff events and ;5% during high events. The ice-meltderived runoff was strongly related to the glacier area. For dry years (low precipitation years), the simulations showed that ice-melt-derived runoff from glacierized catchments averaged 15% of the annual runoff, while the maximum annual value for a specific catchment (538S) was 90%.