A Changing Hydrological Regime: Trends in Magnitude and Timing of Glacier Ice Melt and Glacier Runoff in a High Latitude Coastal Watershed

With a unique biogeophysical signature relative to other freshwater sources, meltwater from glaciers plays a crucial role in the hydrological and ecological regime of high latitude coastal areas. Today, as glaciers worldwide exhibit persistent negative mass balance, glacier runoff is changing in both magnitude and timing, with potential downstream impacts on infrastructure, ecosystems, and ecosystem resources. However, runoff trends may be difficult to detect in coastal systems with large precipitation variability. Here, we use the coupled energy balance and water routing model SnowModel‐HydroFlow to examine changes in timing and magnitude of runoff from the western Juneau Icefield in Southeast Alaska between 1980 and 2016. We find that under sustained glacier mass loss (−0.57 ± 0.12 m w. e. a−1), several hydrological variables related to runoff show increasing trends. This includes annual and spring glacier ice melt volumes (+10% and +16% decade−1) which, because of higher proportions of precipitation, translate to smaller increases in glacier runoff (+3% and +7% decade−1) and total watershed runoff (+1.4% and +3% decade−1). These results suggest that the western Juneau Icefield watersheds are still in an increasing glacier runoff period prior to reaching “peak water.” In terms of timing, we find that maximum glacier ice melt is occurring earlier (2.5 days decade−1), indicating a change in the source and quality of freshwater being delivered downstream in the early summer. Our findings highlight that even in maritime climates with large precipitation variability, high latitude coastal watersheds are experiencing hydrological regime change driven by ongoing glacier mass loss.

freshwater reservoirs, with the ability to temporarily store water over diurnal, seasonal, and long-term (decadal to millennial) time scales (Jansson et al., 2003). Watersheds containing even as little as 5% glacier cover exhibit modified flow patterns compared to their ice-free equivalents, with lower annual and monthly variability, and with a maximum seasonal flow in mid-summer (Fountain & Tangborn, 1985). These differences arise because while runoff from non-glacierized watersheds is dominated by precipitation, runoff in glacierized basins is energy balance dominated (Lang, 1986).
Watersheds downstream of glaciers with persistent negative net mass balance display a distinctive longterm streamflow pattern. This pattern is characterized initially by increasing discharge due to higher rates of glacier mass loss up until a maximum (often referred to as "peak water" (Gleick & Palaniappan, 2010)), followed by decreasing discharge due to shrinking glacier area and volume (Jansson et al., 2003). Whether or not a particular glacierized basin or region has passed peak water is linked to several factors. Huss and Hock (2018) found through a global glacier mass balance modeling study that percent ice cover and absolute glacier size exhibit controls over the timing of peak water. Similarly, Moore et al. (2009) identified for Western North America that northern basins with larger glaciers still show increasing runoff while basins with smaller glaciers further south have already passed peak water. On the other hand, a study by Carnahan et al. (2018) identified through glacier flow modeling that glacier dynamics (characterized by glacier response times, linked primarily to climate and slope) and landscape evolution (i.e., vegetation succession after deglaciation) had a roughly equal impact on basin runoff in response to glacier retreat. Together, these findings indicate that peak water is likely to occur at different times in different regions.
Knowing whether a glacierized watershed is pre-or post-peak water is crucial information for downstream concerns such as infrastructure, ecosystems, and ecosystem resources (Moore et al., 2009). In a study that forecast glacier streamflow to 2100, the large glaciers of the Gulf of Alaska were predicted to experience peak water the latest (between 2060 and 2070) of all regions globally (Huss & Hock, 2018). However, the fate of individual watersheds was less certain due to large intrabasin variability and calibration to regional glacier mass balance observations rather than local runoff measurements.
From an ecological perspective, freshwater from glaciers-whether from melted glacier ice, melted firn, or terrestrial water that has passed through a glacier system-carries a unique biogeochemical signature relative to other freshwater sources. For example, glacier runoff has been found to control fluxes of limiting nutrients crucial for primary productivity in riverine and marine environments. A study on streams discharging the Juneau Icefield found that glaciers serve as an important source of phosphorus and nitrogen in those streams (Hood & Scott, 2008), while nearby glacierized watersheds in Alaska such as the Copper River have proven a critical source of iron to the Gulf of Alaska (Crusius et al., 2011). Glacier meltwater also serves as a major source of bioavailable organic carbon to both riverine (Fellman et al., 2015) and near-shore marine ecosystems Lawson et al., 2014). Moreover, glacier runoff possesses physical properties that are distinct from other terrestrial water sources. Hood and Berner (2009) found that both summer stream turbidity and water temperature can be predicted by the percentage of glacier cover within a basin. These physical conditions are in turn critical for biological productivity at all trophic levels.
The Juneau Icefield, a coastal temperate icefield adjacent to the Gulf of Alaska, is one of the largest icefields in North America. This area experiences extreme amounts of precipitation (Pelto et al., 2013) and among the highest variability in precipitation of any climatic zone in Alaska (Bieniek et al., 2014), both of which may obscure runoff trend detection. In addition to downstream riverine and nearshore marine ecosystems, the icefield lies in close proximity to the coastal city of Juneau, Alaska, such that community infrastructure, including bridges over glacial rivers and residential areas prone to flooding from glacial outburst floods (Kienholz et al., 2020), is also closely linked to physical changes in the icefield.
To assess these changes, several recent studies have evaluated glacier mass balance of the Juneau Icefield. These have primarily relied on geodetic approaches (e.g., digital elevation model differencing) that determine bulk volume loss between two known dates. Despite sourcing imagery from different satellite sensors and covering different time spans, all studies calculated negative glacier-wide mass balance rates over the investigated periods between 1962 and 2016 (Berthier et al., 2010(Berthier et al., , 2018Larsen et al., 2007;Melkonian et al., 2014). A recent study has also modeled future glacier mass balance for the icefield under different climate scenarios, projecting a volume loss of 58%-68% of the icefield by 2100 (Ziemen et al., 2016). This estimate falls on the upper end of regional projections of a 32%-58% loss of Gulf of Alaska glaciers as a whole (Hock & Huss, 2015).
Given the close coupling to surrounding ecosystems and infrastructure, and its persistent state of negative mass balance, the purpose of this study is to examine how components of runoff from the western Juneau Icefield have changed over the past several decades. In particular, we leverage a distributed, high-resolution model to evaluate: (a) trends in the annual or seasonal volume of total runoff, glacier runoff, and glacier ice melt; (b) shifts in timing of the onset or end of glacier runoff and/or ice melt season; (c) shifts in winter glacier runoff events or volume, and (d) changes in timing or magnitude of total runoff, glacier runoff, and glacier ice melt. This study is the first to examine recent changes in timing and magnitude of different hydrological cycle variables in this region and, in turn, to assess whether trends of increasing or decreasing runoff can be detected in a high latitude maritime environment. These findings provide key information for socio-ecological systems downstream, and leave us better poised to project future changes in ongoing climate change.

Study Area
Bordered by mountain ranges spanning from sea level to >5,000 m a.s.l., and with a maritime climate that delivers an average of 2 m w. e. and a peak of 7 m w. e. of precipitation per year (Daly et al., 2008), the Gulf of Alaska coastline is characterized by both extensive glacier cover and extreme volumes of freshwater runoff. Unlike other major watersheds in North America that are dominated by large rivers, 78% of runoff into the Gulf of Alaska is delivered from the steep topography to the coast via short (∼10 km average), small drainages (Neal et al., 2010). In coastal Alaska, glacier termini often lie below treeline, placing glacier ice directly adjacent to the mixed forest of the northern Pacific temperature rainforest. Together, these qualities set up a tight coupling between ice and snowmelt from alpine terrain and the nearshore marine ecosystems downstream.
The Juneau Icefield (Figure 1), centered at 58.9°N and 134.2°W, spans the coast mountains between Southeast Alaska, USA, and Northwestern British Columbia, Canada. It is the third largest icefield in North America with an area of >3,700 km 2 and elevations ranging from sea level to ∼2,300 m a.s.l (Kienholz et al., 2015). All outlet glaciers are currently lake-or land-terminating although, as it finishes a tidewater glacier cycle advance (Truffer et al., 2009), the large (∼725 km 2 ) Taku Glacier is ∼60% protected by a shoal moraine with the remaining portion of the terminus abutting a proglacial lake and short river.
Although the highest elevations receive snowfall throughout the year, C-band synthetic-aperture radar reveals that snow and/or ice melt occurs over the entire icefield during July and August (Ramage et al., 2000). Moreover, because temperatures frequently hover near the freezing point on the coast, low elevations may see ice melt and rain throughout the year. In addition to typical patterns of increasing precipitation with elevation, the icefield also experiences a strong decreasing precipitation gradient from southwest to northeast (i.e., with increasing distance from the coast) due to the prevalence of southwesterly weather systems moving inland from the Gulf of Alaska (Royer, 1998;Stabeno et al., 2004). These patterns are evidenced both in measurements (Pelto et al., 2013) and mass balance modeling studies (Roth et al., 2018;Ziemen et al., 2016).
The spatial domain in this study comprises all terrain draining the western portion of the Juneau Icefield directly to the coast (Figure 1). Though we calculate glacier mass balance for the entire icefield for purposes of calibration, we focus our analysis of runoff on those watersheds of the icefield that supply direct runoff to marine ecosystems. This amounts to a spatial domain of 6,405 km 2 , of which 2,860 km 2 or 44% is glacier ice covered.

Data & Methods
In remote and rugged settings with scarce ground observations, glacio-hydrological models can help fill knowledge gaps about the hydrological regime at high spatial and temporal resolution. To estimate glacier mass balance and total runoff for 1980 to 2016, we use the energy and mass balance model SnowModel (Liston & Elder, 2006a), coupled with the SoilBal routine for calculating evapotranspiration over ice-free areas YOUNG ET AL. 10.1029/2020WR027404 (Beamer et al., 2016), and the linear reservoir runoff routing routine HydroFlow  ( Figure 2). These model routines are described below, as are the data and approaches used for initialization, calibration, and validation.

SnowModel
SnowModel is a distributed surface energy and mass balance model for simulating snow distribution and evolution in terrain where snow and ice are present (Liston & Elder, 2006a). It uses meteorological, elevation, and surface type data as inputs, and accounts for all first-order processes involved in snowpack evolution, including: snow accumulation; forest canopy interception, unloading, and sublimation; snow-density evolution; and snowpack and ice melt. SnowModel is comprised of several routines: (a) MicroMet, (b) En-Bal, and (c) SnowPack ( Figure 2).
MicroMet is a quasi-physically-based data assimilation and interpolation routine that distributes coarse-resolution meteorological forcing over high-resolution topography (Liston & Elder, 2006b). MicroMet adjusts coarse-resolution climate data in two ways: (a) all available data are spatially interpolated over the domain, and (b) physical submodels are applied to generate more realistic values at each grid cell and time step. MicroMet also estimates solar and incoming longwave radiation based on topography, location, and time of year, as well as cloud cover based on relative humidity and temperature. Limitations of the routine include a focus on simple topographic relationships (e.g., precipitation interpolation between stations does not YOUNG ET AL.  include an orographic component), and a lack of two-way feedback (near-surface atmospheric conditions do not change in response to e.g., new snow deposition) (Liston & Elder, 2006b).
EnBal performs surface energy balance calculations in response to atmospheric conditions generated in MicroMet. The energy balance is calculated by: where α is the surface albedo (for either melting snow below a forest canopy, melting snow in a non-forested area, fresh snow, or glacier ice), Q si is solar radiation reaching the terrain surface, Q li is incoming longwave radiation, Q le is emitted longwave radiation, Q h is the turbulent exchange of sensible heat, Q e is the turbulent exchange of latent heat, Q c is conductive energy transport, and Q m is the energy available for melt (Liston & Elder, 2006a). Energy terms are added at the snow-or ice-atmosphere interfaces, and any surplus energy is assumed to be available for snowmelt or for glacier ice melt once overlying snow has been removed (Mernild et al., 2006).
SnowPack, described in Liston and Elder (2006a), simulates snow depth and snow water equivalent evolution within the snowpack based on precipitation and melt energy. Compaction-based snow densification follows the formulation of Anderson (1976), wherein density changes in response to snow temperature and the weight of overlying snow. Density changes also occur by snow melting and rain-on-snow events, which redistribute water through the snowpack and increase density until a maximum is reached. Here, we use SnowPack as a single-layer model, given our focus on bulk snow water equivalent. After the removal of the snowpack through melting, SnowModel initiates melting at the glacier ice surface if ice is present.
SnowModel has some limitations. To avoid infinite snow accumulation at high elevations during multi-year simulations, each year's end-of-summer snowpack over glacier cells is reset to zero under the assumption that residual snow is converted to glacier ice. This prevents the formation of firn, which can play a distinct YOUNG ET AL.
10.1029/2020WR027404 5 of 30 role in glacier mass balance. Moreover, SnowModel is a surface model and does not include a glacier flow model to redistribute ice mass under climate forcing. SnowModel also does not account for changes in either glacier extent by retreat or area-altitude distribution by thinning or ice flow and instead keeps a constant surface and extent representing conditions during a reference year/period (Section 3.2.1). Finally, SnowModel does not include snow and ice mass loss due to dynamic processes, such as frictional melting from viscous heating (internal deformation of the ice) or sliding at the glacier bed (Mernild et al., 2014). We discuss these limitations further in Section 5.1.4.
SnowModel has been applied in a number of glaciology investigations at similar spatial scales as our study, including in Alaska and Greenland (Liston & Hiemstra, 2011;Liston & Sturm, 2002;Mernild et al., 2006Mernild et al., , 2007Mernild et al., , 2010Mernild et al., , 2015Mernild et al., , 2017. Recently, SnowModel has also been applied along with the SoilBal and HydroFlow routines to model freshwater discharge from 1980 to 2014 for the full Gulf of Alaska watershed (Beamer et al., 2016), a study which informs several of our model configuration choices.

SoilBal
SoilBal, a soil moisture routine, was developed by Beamer et al. (2016) to introduce evapotranspiration (ET) into SnowModel-HydroFlow for ice-free and vegetated landscapes. SoilBal first calculates potential evapotranspiration (PET) using the Priestley-Taylor equation (Priestley & Taylor, 1972), which uses only daily air temperature and top-of-canopy net radiation as input data. The Priestley-Taylor equation has been applied to many forested landscapes (see Komatsu (2005) for a review of studies) and has been found to outperform more complex formulations for a mixed temperate mountainous forest (Shi et al., 2008). After PET is calculated, a soil water balance (Hoogeveen et al., 2015) is solved using inputs of PET, runoff from SnowModel, and gridded soil water storage. SoilBal produces daily grids of actual evapotranspiration, surface, and base flow runoff. The latter two are summed and drive the hydrological routine HydroFlow.

HydroFlow
Using this surface and base flow runoff from SoilBal, the HydroFlow routine next simulates routing to downslope areas   (Figure 2). In HydroFlow, each grid cell acts as a linear reservoir (i.e., a reservoir with discharge linearly proportional to water input) that transfers water from itself and upslope cells to the downslope cell, creating a linked flow network. HydroFlow assumes that within each cell there are two transfer functions with two time scales, each associated with different routing mechanisms. Runoff enters first into the slow-response reservoir, which accounts for the time it takes for water transport through the snow, ice, and soil matrices. Water is then routed via the fast-response reservoir, which generally represents some form of channel flow, such as glacier runoff routing or streamflow. Residence time coefficients for each reservoir in each grid cell are a function of many elements, including: surface slope; snow, ice, and soil porosity; snow temperature (cold content); density of glacier crevasses and moulins; hydrostatic water pressure; and soil and land-cover characteristics. HydroFlow therefore assigns residence time coefficients and velocities for four dominant surface types that account broadly for these processes: snow-covered ice, snow-free ice, snow-covered land, and snow-free land. A coupled system of equations solves for slow-and fast-response flow, yielding a discharge hydrograph for each grid cell. A full description of HydroFlow is available in Liston and Mernild (2012).

Model Configuration
Our model simulations cover the water years between October 1, 1980 and September 30, 2016 and are run using a daily time step and grid cell size of 200 × 200 m. Figure 1 shows our model spatial domain, which encompasses the full extent of all observational datasets used for calibration and validation (described below). For this study, unless otherwise specified, reported findings on glacier mass balance include ice within the red outline of the Juneau Icefield, in order to match recent estimates from both Berthier et al. (2018) and Ziemen et al. (2016). However, when discussing freshwater runoff, our spatial domain encompasses all western Juneau Icefield terrain that drains directly to the coast ("coastal watershed domain" in Figure 1). We do not include terrain that routes freshwater into large interior rivers (i.e., the Taku or Yukon rivers, with drainage areas of 17,000 km 2 and 850,000 km 2 ) given the size disparity (the Mendenhall River is the largest coastal drainage at 289 km 2 ) and different climatological regimes (Bieniek et al., 2012). We focus our analysis instead on the unique hydrological regime of the short and steep coastal drainages, particularly given their relevance for downstream estuaries and their prevalence throughout high latitude coastal regions in Alaska (e.g., Glacier Bay, Prince William Sound) and beyond (e.g., Patagonia, New Zealand, Norway).
To evolve the snowpack and route water through the landscape, SnowModel-HydroFlow requires the following data sets.

Elevation, Land Cover, and Soil Type
We use a digital elevation model (DEM) from the United States Geological Survey (USGS) National Elevation Data set (available at https://nationalmap.gov/elevation.html), representing surface topography from the early 2010s as measured by Interferometric Synthetic Aperture Radar data. Data are available at a 1 arcsec (∼30 m) resolution over ∼95% of the domain, and 2 arcsecs (∼60 m) over portions of Canada. The DEM is hydrologically corrected (i.e., depressionless) and we resample to 200 m resolution. We do not modify glacier surface elevations or extents through the 1980 to 2016 model period given that DEMs and glacier outlines for the 1980s are not available for the full icefield.
Land cover classes are obtained from the 2011 North American Land Change Monitoring System (NAL-CMS) and are available at a 30 m resolution (Homer et al., 2015). We resample to 200 m and reclassify to the surface classes defined in Liston and Elder (2006a), including such categories as coniferous forests, alpine meadow, bare rock, and glacier ice. Classes primarily impact model parameter space via albedo (see Equation 1), and control surface processes after snow is removed (e.g., water routing through soils or over bare rock, or continued melting of glacier ice). To more accurately delineate glacierized areas, we also modify the NALCMS grid using higher precision glacier outlines from the mid-2000s from the Randolph Glacier Inventory (RGI) v6.0, available at https://www.glims.org/RGI/rgi60_dl.html (Kienholz et al., 2015;Pfeffer et al., 2014). We do not update surface type information related to for example, vegetation succession after deglaciation during the model period, due to a lack of information dating back to the 1980s.
To classify soil types, we use the gridded Harmonized World Soil Data set (HWSD) version 1.2 (available at http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/) (Fischer et al., 2008), which we resample from 1 km resolution to 200 m. Soil textures are similarly obtained from the HWSD, corresponding to USDA soil moisture texture classes. We use a uniform soil depth of 1,000 mm, the common reference depth used in HWSD and in line with the value used by Beamer et al. (2016) for a Gulf of Alaska SnowModel-HydroFlow modeling study.

Meteorological Data
SnowModel requires forcing of daily temperature, relative humidity, wind speed and direction, and precipitation. We use reanalysis data from NASA's Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) (Gelaro et al., 2017), available at http://gmao.gsfc.nasa.gov/reanalysis/MER-RA-2/. We choose this product given that in their modeling study on freshwater runoff to the Gulf of Alaska, Beamer et al. (2016) found that Version 1 of MERRA (Rienecker et al., 2011) performed best in reproducing measurements of point glacier mass balance and local domain streamflow out of several reanalysis options. MERRA-2 has in turn been found to perform better in North America than MERRA-1 for precipitation, and snow amounts in particular have a lower bias and better correlation to reference data in neighboring parts of Canada (Reichle et al., 2017).
We compare the product to observational meteorological records within our domain and discuss the outcomes in Section 4.2.

Model Calibration Datasets
To constrain our estimates of glacier mass change and freshwater runoff for the Juneau icefield, we use multiple calibration datasets including: snow water equivalent observations, a geodetic glacier mass balance estimate, and streamflow measurements.

Snow Water Equivalent
Point observations of snow water equivalent (SWE) used for model calibration ( Figure 2) are obtained from field campaigns conducted as part of this study in collaboration with USGS. In late April of each 2013, 2014, YOUNG ET AL.

10.1029/2020WR027404
and 2015, we carried out SWE observations at six locations along the Gilkey Glacier centerline between 300 and 1,900 m a.s.l and on the Taku Glacier between 1,000 and 1,870 m a.s.l. SWE values were derived using measured density profiles obtained from snow core samples, representing stratigraphic balances, and converted to SWE following standard glaciological protocols (Østrem & Brugman, 1991). Data for the Gilkey Glacier are available at Young (2019) and for the Taku Glacier at McNeil, Campbell, et al. (2019).
We also incorporate helicopter-borne ground-penetrating radar (GPR) observations collected by USGS along the Taku Glacier and Gilkey Glacier centerlines in spring 2014 and 2015, in collaboration with the above-described field campaigns. Raw GPR data were sourced from O' Neel et al. (2018), and were processed by USGS and converted to snow depths using the methods described in McGrath et al. (2015). Density data were sourced from six contemporaneous snow cores measured along each corresponding flight path, where densities were linearly interpolated between locations by the increment 1/n, where n are equally spaced observations between core sites. This data set is equivalent to ∼121,000 and ∼39,000 SWE point observations in 2014 and 2015, that we averaged to single annual values within each model grid cell.

Geodetic Glacier Mass Balance
Several studies have derived geodetic bulk volume loss estimates for the Juneau Icefield, including Larsen et al. (2007) 2018)). The mass balance result from Berthier et al. (2018), calculated from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) imagery, agrees closely with laser altimetry and is therefore the current best estimate overlapping with our study interval.
In our calibration process, we aim to reproduce the mean annual glacier-wide mass balance rate from Berthier et al. (2018) for the same spatial domain (i.e., the glacier outline for the Juneau Icefield obtained from the Randolph Glacier Inventory v5.0, which is identical to v6.0). Because the early and late ASTER scenes used in Berthier et al. (2018) represent mosaics of different acquisition dates, the authors listed their geodetic estimate as generally spanning 2000 to 2016, without citing specific start or end dates. For comparison to the model, we select start and end dates as the beginning and end of the associated water years, that is, October 1, 2000 and September 30, 2016.

Streamflow Measurements
Semi-continuous time series of discharge data (Q) are available for four stream gauges in the Juneau area ( Figure 1), including three streams instrumented by the USGS (Mendenhall River, Lemon Creek, and Montana Creek; data available at https://waterdata.usgs.gov/nwis/rt) and one (Cowee Creek) monitored by researchers at the University of Alaska Southeast. Data are available for different time periods for each. The four instrumented basins represent a range of size above the gauge, percent glacier cover, elevation range, and distance between glacier outflow and gauge (Table 1). This range increases our ability to test model performance across different flow regimes.

Calibration Approach
To correctly characterize glacier mass change and freshwater discharge, we adopt a two-stage calibration approach. The first stage is automated within SnowModel, and leverages the built-in data assimilation routine SnowAssim ( Figure 2). SnowAssim is used to compile and interpolate all available snow water equivalent data (Liston & Hiemstra, 2008). The scheme uses a Barnes objective analysis algorithm that applies a Gaussian weighting function whereby the weight of each measurement on a grid point's value decreases with increasing distance from the measurement (Barnes, 1964(Barnes, , 1973Koch et al., 1983). It does not account for uncertainties in the snow observations. SnowAssim is run prior to regular SnowModel simulations. It optimizes interpolation by calculating the differences between observed and modeled snow values and retroactively applying multiplicative factors to precipitation values or melt factors to create improved fields prior to the assimilated observations. SnowModel is then run again using the new precipitation fields as input. This early, automated form of calibration generates more accurate spatial distribution of snow depth and SWE throughout the season rather than only at the time of observation.
For the second calibration stage, we begin by identifying which of the SnowModel-HydroFlow parameters to treat as tuning parameters and which can be prescribed. SnowModel-HydroFlow has an extensive suite of parameters, many of which have been constrained empirically or from modeling experiments. Based on a review of other SnowModel-HydroFlow studies in glacierized terrain, we initially select seven parameters for tuning: glacier albedo, fresh snow albedo, melting (non-forested) snow albedo, monthly precipitation lapse rates, monthly temperature lapse rates, and factors for modifying each the slow and fast reservoir velocities in the HydroFlow routing routine (acting to increase or decrease fluid residence time). We adopt a traditional grid search approach to tuning model parameters, beginning with a broad search across the parameter space, as guided by the literature (Table 2). Preliminary simulations indicate that results are relatively insensitive to values of fresh snow albedo and the factor for slow reservoir velocities. We therefore focus our calibration on the remaining five parameters. We zoom in on narrower ranges of physically realistic values with a finer grid. All other SnowModel parameters are set to default SnowModel values, a select list of which is also shown in Table 2.
We next establish calibration datasets and appropriate metrics to evaluate model performance. We first prioritize achieving a match between our modeled glacier mass balance and the long-term geodetic estimate from Berthier et al. (2018) for the same 2000 to 2016 period. We aim to minimize is the annually averaged glacier-wide mass change rate for 2000 to 2016 from the model and  geo B is −0.68 ± 0.15 m w. e. a −1 . We next compare HydroFlow output of discharge (Q) to streamflow data for the four local drainages, aiming to maximize r 2 and to obtain Nash-Sutcliffe Efficiency (NSE) nearest to 1 (Nash & Sutcliffe, 1970). We generate separate statistics for each instrumented basin, but prioritize matching those with the highest percent glacier cover (Mendenhall River, 56%, and Lemon Creek, 46%).
In summary, we prioritize our performance metrics in the following order: nearest to 0 for each simulation's 2000 to 2016 annually averaged glacier-wide mass balance rate; (b) NSE nearest to 1 for streamflow discharge, prioritizing the statistics for more glacierized basins first. For our final time series analysis, we identify out of our 215 simulations all those that yield: (a) an annually averaged glacier mass balance rate for the full icefield that falls within the error bounds of the  geo B goal value for October 1, 2000 to September 30, 2016 (priority (a), and; (b)) a glacierized basin NSE of ≥0.75 (priority 2). This yields an ensemble of 16 simulations. Among this ensemble is a midpoint member with a mass balance rate most closely matching the goal value, that is, with  diff B = 0, as well as two ensemble end members whose mass balance rates correspond to the upper and lower limit of the Berthier et al. (2018) estimate error bars. We use the midpoint simulation for most of the analyses to follow, and use the end member simulations to derive upper and lower estimates of uncertainty.

Model Validation
For independent validation, we compare model results to meteorological data from the nearest long-term automated weather stations, point ablation observations, and a time series of terrestrial water changes for the Juneau Icefield area derived from gravimetry data.

Meteorological Variables
To assess the performance of the MicroMet routine, we compare daily MicroMet-interpolated MERRA-2 air temperature fields to observations from two nearby National Oceanic and Atmospheric Administration (NOAA) airport weather stations (Figure 1). Data are available for the Juneau station from https://www. ncdc.noaa.gov/cdo-web/datasets/GHCND/stations/GHCND:USW00025309/detail, and for the Skagway station from https://www.ncdc.noaa.gov/cdo-web/datasets/GHCND/stations/GHCND:USW00025335/detail (Table 3). We compare mean monthly modeled fields in order to identify possible seasonal biases. We use this data set for validation rather than apply bias corrections to the model fields for temperature or precipitation because both weather stations are biased to low elevations, and we have no additional information YOUNG ET AL. 1.26 Note. That we also list a selection of prescribed parameters that are not varied.  for spatially distributing corrections across the complex topography between stations. We therefore assume that any biases are accommodated for by adjustment to the tuning parameter suite.

Ablation Observations
We also compare model results to point snow and ice ablation observations at stake sites from several sources. We glean observations for Taku Glacier from Criscitiello et al. (2010) and McNeil, Campbell, et al. (2019), from Lemon Creek Glacier from Criscitiello et al. (2010) and McNeil, Sass, et al. (2019), for Mendenhall Glacier from Motyka et al. (2002) and Boyce et al. (2007), and from our field campaigns on the Gilkey Glacier in 2013-2015 (Young, 2019). Snowmelt values were calculated by subtracting SWE equivalent values between snowpacks at known start and end dates. Ice melt values used exposed stake height changes multiplied by an assumed glacier ice density of 900 kg m −3 . All ablation observations are compared to model output extracted for the same location and time span.

GRACE Gravimetry Data
Finally, to independently validate model results for the full water balance in this region, we leverage gravimetry data from the Gravity Recovery and Climate Experiment (GRACE) tandem satellites, which measured gravity changes of all Earth system components between 2003 and 2016. We choose GRACE data from NASA Goddard Space Flight Center Geodesy Laboratory's high resolution v2.4 mass concentration (mascon, i.e., grid cell) solution, which provides mass change estimates at ∼30-day intervals and 1 ° × 1 ° (∼12,390 km 2 ) resolution . Data are available at https://earth.gsfc.nasa.gov/geo/data/grace-mascons. This solution represents the full terrestrial water budget after other signals such as atmospheric loading and Earth and ocean tides have been removed. The solution includes snowfall, rain, and runoff from nonglacierized and glacierized terrain, including glacier ice melt. We choose this GRACE product because it is one of few that explicitly corrects for local mass increases from post-Little Ice Age disintegration of the Glacier Bay icefield (Larsen et al., 2005), as estimated using the ICE-5G glacial isostatic adjustment model (Peltier, 2004). This GRACE product also compares well with regional-scale mass balance model simulations for the Gulf of Alaska (Beamer et al., 2016;Hill et al., 2015) and to mass loss estimates from NASA's Ice, Cloud, and Land Elevation Satellite (ICESat) . Moreover, this solution is among the first to provide information for constructing 95% confidence intervals on mass changes for individual mascons based on estimates of noise and leakage, as detailed in Loomis et al. (2019).
For comparison of our model results to the GRACE time series, we include in our model validation domain all terrain within two adjacent GRACE mascons surrounding the icefield (shown together in one box as "GRACE domain" in Figure 1). We extract this spatial domain and select mass change estimates at dates corresponding with the mid-points of the GRACE time series monthly averages. We calculate the long-term mass loss trend by fitting an annual sinusoid to data using a least squares approximation. Individual annual amplitudes are calculated by subtracting annual minima from maxima, an approach deemed appropriate for the Gulf of Alaska region due to its clean seasonal signal relative to noise .

Water Balance, Glacier Mass Balance, and Runoff Calculations
Using SnowModel-HydroFlow as described above, the water balance for our domain is calculated by: where S is the volume of water stored within the seasonal snowpack, glacier ice, or top 1 m of soil; P is precipitation input (rain or snow); R is runoff (defined as the water immediately available for routing to downslope areas); ET is evapotranspiration; and SU is sublimation at the snow surface. Dot notation indicates that all quantities are taken to be rates. Note that because none of the glaciers within the domain are ocean-terminating, we do not include marine iceberg calving or submarine melt within Equation 2.
Although several glaciers are lake-terminating, previous studies on the Mendenhall Glacier-the largest of the lake-terminating glaciers-revealed that iceberg calving represents only 4%-6% the amount of ice lost through surface melt (Boyce et al., 2007;Motyka et al., 2002). Similar to Ziemen et al. (2016), we therefore consider ice discharge into lakes to be a small component of Juneau Icefield glacier mass balance, and an even smaller part of water balance of the coastal watershed.
In SnowModel, runoff R is all water that is, immediately available to be routed downstream, and is calculated by: where R gim is glacier ice melt, R sm is snowmelt that does not refreeze or fill pore space within the snowpack, R rbs is rain on bare surfaces (i.e., rain that does not fall onto snow or soil substrates), and rss n R and rss o R are rain on already-saturated snow and soil.
We note that the term "glacier runoff" is used ambiguously within the literature and often represents different physical quantities (O'Neel et al., 2014;Radić & Hock, 2014). For our purposes, we define glacier runoff as all runoff in Equation 3 produced over glacierized cells. This formulation is identical to two studies that modeled runoff for the Gulf of Alaska (Beamer et al., 2016;Neal et al., 2010) as well as to the quantity defined conceptually in O'Neel et al. (2014) as total runoff from the glacier surface (concept 5). We use the term "glacier ice melt" separately, to denote meltwater from the glacier surface only after snow cover has been removed (i.e., it is one component of glacier runoff). We calculate both quantities throughout the study.
We calculate the area-averaged glacier mass balance using Equation 2 over glacierized grid cells only (noting that evapotransporation (ET) goes to zero over glacier surfaces). Glacier mass balance therefore represents a portion of the full spatial domain's water storage S. The contribution of non-glacierized cells makes up the remaining portion. All comparisons of model output to stream gauge instruments are comparisons to: that is, discharge Q (a flux) is all runoff that has been routed to a known gauge location, after evapotransporation ET has been taken into account.
Finally, comparisons of model output to GRACE data are to water storage S, given that the GRACE satellites measure all changes in water mass distribution over Earth's surface.

Trend Analyses
We evaluate trends in magnitude and timing of hydrological variables (e.g., total runoff, glacier runoff, glacier ice melt, and water balance), integrated over the full spatial domain draining west to the coast. For trends in magnitude, we examine spatially and temporally integrated quantities including annual volumes of total runoff, precipitation, glacier runoff (the sum of ice melt, snowmelt, and rain on the glacier surface), glacier ice melt (i.e., melt at the glacier surface after snow has been removed), and water balance. We also identify maximum and minimum daily values for each year. Further, we examine volumes of glacier runoff and ice melt for spring and summer seasons, where each season's start and end dates are defined by the maximum, minimum, and inflection points of the domain-and temporally averaged annual air temperature climatology from the MicroMet-interpolated climate input data. Here, "winter" falls between December 24 to April 6, 'spring' is April 7 to July 17, 'summer' is July 18 to October 11, and "fall" is October 12 to December 23. Finally, we assess cold season volumes of glacier runoff and glacier ice melt. Here, the cold season is defined as the period between late-fall termination and spring onset of glacier runoff and ice melt, which correspond to the latest and earliest dates that respectively follow or precede a period of at least two weeks of glacier runoff/ice melt below a near-zero threshold. This two-week criteria was chosen out of several algorithms for best reproducing manually selected dates.
For additional trends in timing, we use the complete time series to test for trends in: day of year of minimum daily volumes of total runoff and water balance; day of year of glacier runoff and glacier ice melt onset and end, as well as the length of the season in between; and number of non-zero days of cold season glacier runoff and ice melt. For trends in the timing of peak flows (i.e., maximum daily volumes of total runoff, water balance, glacier runoff, and glacier ice melt), we test for day of year trends in a time series smoothed with a 14-day running mean in order to minimize the effect of extremes.
Trends are detected using the Mann-Kendall test for significance, a non-parametric test (i.e., data do not have to meet the assumption of normality). Trends are calculated using the Theil-Sen estimator, a non-parametric approach that is, more robust against outliers than simple linear regression, making it well-suited to, and commonly used in, hydrological applications (Helsel & Hirsch, 2002). To identify the statistical significance of each trend, we report a harmonic mean p-value, a formulation for combining p-values from tests that cannot be guaranteed to be independent (Wilson, 2019), for example, model simulations with the same input data and physics but variation in parameter values. We calculate a harmonic mean p-value for every trend by equally weighing our midpoint and two end member simulation p-values.
In reporting our findings, we take an approach that extends beyond the traditional method of judging results as meaningful solely by the p-value ≤0.05 criteria. This has been challenged in recent years, citing limitations such as variation in p-value statistics across replicate studies (Halsey et al., 2015) and difficulty in interpreting results when the p-value is high and the null hypothesis cannot be rejected (Cohen, 2016). We turn instead to recommendations from Halsey (2019) and Tomczak and Tomczak (2014) to include in our analysis a measure of effect size (which in our case is the trend itself) as well as 95% confidence intervals around the trend, to provide additional insight into the range of possibilities that are reasonably likely. We also heed advice from Amrhein et al. (2019) that including factors such as background evidence, data quality, and understanding of underlying mechanisms can contribute to meaningful interpretation of statistical results. As such, we include as an interpretive tool for the reader a qualitative assessment of our confidence that a positive trend should be detected, in the context of our full suite of results and a priori current knowledge from the literature for each climatological and hydrological variable.

Model Initialization and Calibration Process
We evaluate the impact of the initial calibration routine SnowAssim by comparing SnowModel on-glacier point SWE estimates to observations from glacier mass balance field and airborne campaigns (Figure 3). We observe that model reproduction improved markedly from r 2 = 0.45 to r 2 = 0.90 and RMSE = 0.45 m w. e. to RMSE = 0.18 m w. e (Figure 2). This highlights that the SnowAssim routine produces more realistic SWE fields irrespective of location or duration between observations.
In the second calibration phase, we succeed in tuning parameters to reproduce the geodetic mass balance rate from Berthier et al. (2018)  simulations that meet this criteria, we focus our primary analysis on the midpoint simulation with a mass balance rate of exactly −0.68 m w. e. a −1 , and consider the ensemble end members-whose mass balance rates are nearest the upper and lower error bounds from Berthier et al. (2018) -to be the limits of our uncertainty. Best-fit parameter values are shown in Table 2. This step of calibrating to a long-term mass balance rate is crucial for correctly characterizing glacio-hydrological systems. Had we not undertaken this step, our initial simulations using SnowModel default parameter values would have yielded a mass balance rate of +0.08 m w. e. a −1 .
Our ability to reproduce observations from stream gauge records on the four instrumented basins varies by the amount of glacier cover (Figure 4). For the two glacierized basins with the largest percent cover, comparison of modeled to observed monthly discharge yields stronger agreement: for Mendenhall River (56% glacier cover), we obtain NSE = 0.84 and r 2 = 0.88, and for Lemon Creek (46% glacier cover), we find NSE = 0.76 and r 2 = 0.82. The model, however, is unable to reproduce many of the large peaks in the daily Mendenhall discharge record, several of which are associated with recent (2011 and on) glacier lake outburst floods from an upstream tributary basin (Kienholz et al., 2020). The model does not include a mechanism to generate these impulsive events. For the two basins that are predominantly forested, modeled to observed agreement is weaker: for Montana Creek (2% glacier cover), we find NSE = −1.37 and r 2 = 0.45, and for Cowee Creek (11% glacier cover), we obtain NSE = −0.81 and r 2 = 0.47. We also note that the Mendenhall River and Lemon Creek watersheds show evidence of seasonal biases between modeled and observed quantities, with the model generally overproducing runoff in summer and under-producing in fall. We discuss this, and provide possible reasons for the modeled-to-observed discrepancy in less-glacierized basins, in Section 5.1.2. Altogether, weighing all four basins according to both above-gauge basin area as well as length of observational record, we calculate a weighted NSE = 0.21 and weighted r 2 = 0.73. We believe this performance is acceptable given that, rather than any one process in isolation, streamflow represents an integration of all glacio-hydrological processes in the watershed, and thereby has the potential to integrate any sources of error with input data as well as model physics into a single metric. Because our model performs well in reproducing other calibration datasets, particularly in glacierized watersheds (e.g., our estimate for the 2000 to 2016 mass balance rate for the Mendenhall Glacier alone is −0.73 m w. e. a −1 , which matches the estimate of −0.73 ± 0.13 m w. e. a −1 from Berthier et al. (2018)), we are confident in the calibrated model performance.

Meteorological Conditions
We find strong correlation between MicroMet-interpolated MERRA-2 air temperature fields and observations from the NOAA airport weather stations in Juneau and Skagway (r = 0.96 and p < 0.001, and r = 0.94 and p < 0.001). However, we find systematic biases between modeled and observed temperatures, when averaged monthly, with lower than observed temperatures in winter months (as large as −2.1°C in Juneau and −5.5°C in Skagway) and higher than observed temperatures in summer months (as large as 2.0°C in YOUNG ET AL.   Table 1. Solid black lines represent 1:1 lines. Note the differing axis scales. Juneau and 2.8°C in Skagway). For daily precipitation, modeled and observed volumes were moderately correlated in both Juneau (r = 0.72 and p < 0.001) and Skagway (r = 0.63 and p < 0.001). Mean monthly modeled fields also overproduced precipitation, particularly in fall and early winter months, with biases between 1.3 and 4.7 mm w. e. d −1 for Juneau and 0.8-2.3 mm w. e. d −1 for Skagway.

Snow and Ice Ablation
The model reproduces independent point snow and ice ablation observations with an r 2 = 0.79 and RMSE = 1.63 m w. e ( Figure 5). The larger RMSE values are not unexpected given the predominance of ablation measurements at lower elevations in the ablation area (60% of the observations are at <800 m a.s.l.), which on large glaciers with undulating surface topography often display substantial local variability that may not be well-captured by the model (e.g., Young et al., 2018). However, we note that the model appears to underpredict melt for more negative point mass balances, which may be related to the aforementioned higher-than-observed precipitation in early fall and winter.

Glacier Mass Balance
Our modeled, tuned annual glacier-wide mass balance rate for the Juneau Icefield is −0.68 m w. e [−0.15, +0.11] a −1 for 2000 to 2016, where values in the square brackets correspond to the asymmetric lower and upper uncertainty bounds from our simulation ensemble end members. Extending to the full model period of October 1, 1980 to September 30, 2016, we calculate a rate of −0.57 [−0.11, +0.12] m w. e. a −1 for the icefield, suggesting an acceleration in recent decades. Finally, for all ice contained within the domain draining to the coast, our model estimates a mass balance rate of −0.81 [−0.08, +0.11] m w. e. a −1 for 1980 to 2016, suggesting that the ice nearest the coast (i.e., to the west of the topographic divide) experiences greater rates of mass loss than the more interior glaciers. Cumulative glacier-wide specific mass balance for the full model period is shown in Figure 6. Annual glacier mass balance over this time period and domain is comprised of, on average, 3.07 ± 0.01 m w. e. a −1 of precipitation, 3.85 [−0.08, +0.10] m w. e. a −1 of glacier runoff, and 0.03 ± 0.01 m w. e. a −1 of sublimation from the snow surface.

Freshwater Runoff
For the watershed encompassing the portion of the Juneau Icefield that drains to the coast ("coastal watershed domain" in Figure 1), we estimate mean annual freshwater runoff of 20.0 [+0.5, −0.4] km 3 a −1 for 1980 to 2016. Of this, 11.0 ± 0.3 km 3 a −1 (or 55 ± 1%) is glacier runoff (i.e., runoff sourced from the glacier surface). The water balance volume we calculate is, on average, −2.1 [+0.4, −0.3] km 3 a −1 , though as we discuss below in Section 5.1 this is believed to be an underestimate of the long-term water storage loss. For ice-only cells, we calculate water storage losses (i.e., glacier volume loss) of 2.4 [−0.3, +0.2] km 3 a −1 for the same time period, which means that glacier volume loss (the percentage of runoff due to the persistent negative mass balance trend, rather than seasonal magnitudes of glacier runoff) comprises 12 ± 1% of total To better understand the linkages between individual water balance components, we assess the correlation between different modeled quantities. We find that annual volumes of glacier runoff and total runoff for the domain are highly correlated (r = 0.95, p < 0.001), while glacier runoff and glacier ice melt are less so (r = 0.82, p < 0.001). Glacier ice melt is also weakly correlated with total runoff (r = 0.67, p < 0.001).
YOUNG ET AL.  Overall, this suggests that glacier runoff is controlled more by precipitation than by glacier ice melt, which we discuss further in Section 5.2.2.

Water Balance and Comparison With GRACE
For the 2003 to 2016 period overlapping with GRACE data availability, we calculate a glacier-wide mass balance rate for all ice in the GRACE domain of −0.51 [−0.18, +0.13] m w. e. a −1 (or −2.5 [−0.9, +0.6] km 3 a −1 ), in close agreement with the GRACE-derived negative trend estimate of −0.55 m w. e. a −1 (−2.7 km 3 a −1 ), as shown in Figure 8a. Correlation between these two time series is robust, with r = 0.95 and p < 0.001 ( Figure 8b). These results showcase the model's ability to reproduce the climatic conditions over the ice-covered portions of the domain that are driving sub-and interannual water storage changes.
However, in comparing GRACE to modeled results for ice and land cells together, we observe that correlation is less strong (r = 0.60, p < 0.001). This discrepancy can be seen in the SnowModel land + ice time series in Figure 8a primarily as a lack of agreement in the overall trend, which is not sufficiently negative at −0.002 m w. e. a −1 . We discuss this further in Section 5.1.3. Nonetheless, our full SnowModel land + ice water balance produces seasonal amplitudes (mean annual accumulation = 25.8 km 3 a −1 , ablation = −26.6 km 3 a −1 ) that are more in line with those from GRACE (18.1 and −21.5 km 3 a −1 ) than those from ice cells alone (9.0 and −12.1 km 3 a −1 ). This result is encouraging given that the GRACE solution we use measures the change of the total of all terrestrial water balance components.

Trends in Magnitude and Timing
We next assess trends in the timing and magnitude of different hydrological variables, and summarize results of trend detection tests in   guide a column with a qualitative assessment of our confidence that a positive trend should indeed be present in each specific variable, given the trend result in context with our full suite of results as well as a priori information.
To help interpret our model output results, we first assess trends in the principal input variables of precipitation and mean air temperature. We find no reliable trend in annual precipitation volume, but do find an increase in mean air temperature of 0.1°C per decade (Table 4). This is consistent with recent analyses of air temperature trends in Alaska, including Bieniek et al. (2014) who found a 0.2°C increase in the northern portion of the Juneau Icefield between 1980 and 2012.
Of all variables tested, the most statistically robust (p ≤ 0.05) trends are related to shifts in timing of the peaks of the 14-day smoothed glacier ice melt curve (occurring 2.5 days earlier per decade) and glacier runoff curve (occurring 4.4 days later per decade) (Figure 9). The day of year of the water balance minimum is also found to be occurring 3.5 days earlier per decade.
From a seasonal perspective, the most statistically robust trends with the largest effect sizes occur in our hydrological variables in the spring season ( Figure 10). We also observe an increase in glacier ice melt in summer.
Among the different hydrological variables examined, the most robust trends are related to glacier ice melt (Table 4). These include the volume of spring glacier ice melt (increasing by 16.5% per decade) and, with slightly less statistical strength, the annual volume of glacier ice melt (9.6% per decade), both of which are visible in Figure 11. Our results also suggest an increase in the magnitude of the maximum daily volume of glacier ice melt (10.2% per decade).
The extreme amounts and large interannual variability in precipitation in this domain (Bieniek et al., 2014) increasingly act to obscure trend detection as the proportion of non-glacier ice grid cells grows in a particular hydrological variable. In other words, when examining volumes, we observe the pattern that trends for glacier ice melt, glacier runoff, and total runoff exhibit respectively smaller proportion change with less robust statistical significance. For example, in spring months, we calculate p-values of 0.05, 0.11, and 0.25, and respective trends of 16.5%, 6.8%, and 2.7% per decade for those three variables (Table 4). This pattern holds true for each spring, summer (not shown in table), and annual periods, and disappears during fall and winter months when glacier ice melt ceases almost entirely.
YOUNG ET AL.

10.1029/2020WR027404
17 of 30 Finally, our results also suggest trends for variables associated with colder months, including an increase in the number of days of non-zero glacier runoff during the cold season (2.4 days per decade), but a decrease in the volume of glacier runoff during winter months (−5.8% per decade) ( Note. Here all variables are defined as positive (e.g., glacier ice melt is positive even though it represents a loss), such that positive/negative trends correspond to increasing/decreasing quantities in all cases. p-values are given by the harmonic mean of individual Mann-Kendall tests for the midpoint, upper, and lower end member simulations, and bold indicates the trends that are statistically strongest. Trends are given by the Theil-Sen slope and a 95% confidence interval is provided for each. The percent change per decade is indicated for the mean trend (column 3) relative to the 1980 to 1989 period. Finally, the last column shows our qualitative assessment of confidence that a positive trend should be present, given our results and in context with the literature (1 = highly confident, 2 = confident, 3 = somewhat confident, -= not confident).  Of the remaining variables tested, none show trends we believe to be significant according to our methods. Of these, fall season volumes show the lowest p-values of any season for all hydrological variables, followed by the winter season. Maximum and minimum daily volumes of the different variables do not exhibit changes in either volume or timing. We also do not detect changes in the frequency of cold season glacier ice melt events. Finally, we do not detect reliable trends in the onset and end of glacier ice melt or glacier runoff, nor in the length of the melt season in between. Any of these variables, however, may show significant trends in future years.
Finally, we observe a significant shift in spatially distributed volumes of freshwater from the beginning and end periods of our model interval ( Figure 12). Precipitation decreases consistently across the full domain, particularly over the central portion of the icefield. This also appears to be reflected in increased glacier melt over the same central icefield region, although those increases are generally more widespread and exhibit a greater negative deviation from the 1980 to 2016 means. Runoff increases over the full icefield between the beginning and end periods, and a decrease in runoff over non-ice terrain aligns with a corresponding decrease in precipitation. Overall, the water balance over the icefield shows a shift from slightly positive during 1980-1990 to strongly negative in 2010-2016.

Model Performance
Overall, our model calibration approach achieves robust agreement with calibrating datasets of snow water equivalent observations, long-term geodetic glacier-wide mass balance, and discharge in highly glacierized basins. These results highlight our ability to effectively combine the suite of different physically based sub-routines needed to reproduce accumulation, ablation, and hydrological processes in these complex, glacierized basins.
YOUNG ET AL.

Parameter Tuning-System Dominated by Ice and Snow Albedo
During the calibration process, we find that the annual glacier mass balance rate-which we aim to match to the estimate from Berthier et al. (2018) as our primary calibration criteria-is most sensitive to the values for glacier ice albedo and melting snow albedo in clearings (i.e., non-forested areas, including over glaciers). We tune both parameters to values on the low end of typical ranges seen in the literature (i.e., 0.30 to 0.40 for glacier ice albedo and 0.40 to 0.50 for melting snow albedo in clearings). The lower values may be explained by the presence of both snow algae (documented on another coastal icefield in Alaska in Ganey et al. (2017), and observed by the first author in the field) as well as dust and black carbon (Nagorski et al., 2019). Both of these light-absorbing impurities contribute to an amplifying feedback process by lowering albedo and increasing melt rates, which in turn consolidates material on the snow surface and further increases melt rates. On the Juneau Icefield, dust and black carbon concentrations in surface snow increase later in the melt season Nagorski et al. (2019), suggesting that snowpack "aging" should be taken into consideration in future melt modeling efforts. Incorporating this process by allowing for monthly varying albedo values would likely improve our SnowModel-Hy-droFlow simulations of late-summer freshwater discharge by increasing glacier ice melt and snowmelt during those months. Modeled glacier mass balance rates were insensitive to the value of fresh/dry snow albedo, consistent with the fact that the coastal Juneau Icefield is dominated by aged or wet snow during the runoff season.
Among the range of precipitation lapse rates we employed, the smallest lapse rates performed best (0.20/0.05 km −1 for January/June) ( Table 2). At the scale of the icefield, this can be explained physically by increases in precipitation with elevation being largely canceled out by decreasing precipitation with distance from the coast Roth et al. (2018). As SnowModel only applies a single lapse rate over the entire domain, we effectively combine these two effects into a small lapse rate value. This pattern in precipitation lapse rates may be equally important in other coastal regions with extreme topography rising steeply from sea level and lying along a strong coastal-to-continental gradient. We also find that normal to shallow temperature lapse rates (3.9/7.7°C km −1 for January/June) perform the best overall, in agreement with well-established findings that glaciers can impose a dampening effect on local atmospheric lapse rates ( Our hydrological simulations reveal that model discharge results are relatively insensitive to the slow reservoir velocity parameter, indicating that most runoff is routed through creeks and streams or over fast-flow terrain such as glacier ice and bare rock. This is supported by the shallow soil reference depth cited in the Harmonized World Soil Data set (Fischer et al., 2008), and by the modest fraction of forest coverage within the model domain (17% forest, 14% shrubland/grasses/meadows).

Challenges With Reproducing Stream Gauge Records
While our model adequately reproduces gauge observations in the two basins with high glacier cover (≥45%), gauge-matching results in the two lesser glacierized basins (≤15%) are weaker. This mismatch is evident as an overproduction of discharge in spring, an underproduction in summer, and an underproduction in winter (Figure 4). These patterns are similar in the more glacierized basins, but to a lesser extent. Spring and summer discharge discrepancies may be explained by our finding that MicroMet-interpolated MERRA-2 air temperature fields are generally higher in spring and lower in summer compared to observations (see Section 4.2), and may therefore generate too much early snowmelt in spring, and too little glacier ice melt in summer. This is consistent with a comparative study of reanalysis products for hydrological ap-YOUNG ET AL. plications that found that MERRA-2 does not maintain snow in mountainous terrain for long enough into spring, likely due to precipitation biases and warm temperatures (Wrzesien et al., 2019). We speculate that these effects may be stronger in the less glacierized basins we modeled given the dominance of snowmelt in spring, with little glacier ice melt contribution in spring or summer.
During winter months, modeled discharge in the less-glacierized basins is near-zero, in contrast to observations that show sporadic discharge. However, modeled precipitation volumes in fall and early winter exceed station observations. A possible explanation for the winter month discharge discrepancy is that because our modeled temperatures are lower than observed during winter months, precipitation events arrive as snow instead of rain, thus adding to the snowpack rather than to discharge.
In general, the lack of reliable mountain precipitation data continues to be a challenge for both climatological and hydrological models. The MERRA-2 product for example relies partly on assimilated station data and partly on model physics to produce precipitation fields for latitudes up to 62.5° (Bosilovich et al., 2015), and station data are scarce in this region, particularly at elevation. We underscore the critical need for continuous high-elevation stations in the mountainous regions of Alaska for this purpose. Additionally, there are limitations to downscaling coarse-scale meteorological forcing over complex mountain terrain. For example, MicroMet does not account for orographic effects (i.e. decreased precipitation on leeward slopes), relying instead on a simple elevation dependent precipitation adjustment factor. Moreover, the precipitation adjustment factor applies to solid precipitation only, and not rain, which may also help explain the discrepancy between modeled and observed discharge in the lesser glacierized basins.
Altogether, there is much room for improvement in the characterization of precipitation and particularly snow in complex mountain terrain with sparse observation networks. In the meantime, our model's limited ability to reproduce discharge in less-glacierized basins may lead to increased uncertainty in the magnitudes of spring and winter runoff in those basins in particular. Given our principal goal of examining changes for a 44% glacier covered domain, with an emphasis on glacier changes, we accept this cost.

Agreement With GRACE Highlights Reproduction of Large-Scale Climate Processes
The robust agreement between the model and GRACE (Figure 8), in terms of both long-term trends and time series correlation, emphasizes the model's ability to reproduce meso-and synoptic scale climatic processes that drive sub-and interannual water balance changes over glacierized terrain. We note that the mass balance rate we derive for the larger GRACE domain (−0.51 [−0.18, +0.13] m w. e. a −1 ) is less negative than that for only the Juneau Icefield for the same time period (−0.71 m w. e. a −1 ). We attribute this to inclusion in the GRACE domain of many smaller, higher-elevation glaciers with less negative mass balance rates even at their termini (∼−2 m w. e. a −1 ) relative to the large, low-elevation valley glaciers that dominate the icefield (∼−8 m w. e. a −1 ).
Our finding that modeled seasonal amplitudes for the full land + ice domain are a closer match to those from GRACE than those from the ice-only terrain is consistent with findings for the Gulf of Alaska (Beamer et al., 2016) and the Canadian Arctic Archipelago (Lenaerts et al., 2013). In both studies, seasonal amplitudes from GRACE solutions could only be reproduced by summing together model-generated mass changes over both glacierized and ice-free regions of their modeling domains.
In terms of long-term trends for the full water balance, our model results show a less negative trend than GRACE. This discrepancy is also evident in the Beamer et al. (2016) SnowModel study using MERRA-1. However, using their best-performing climate product (Climate Forecast System Reanalysis), those authors found favorable agreement between trends, suggesting that the long-term ice loss trend within GRACE has been correctly attributed (i.e., that none or little of the trend is attributable to TWS). This is also consistent with Reager et al. (2016), who found little in the way of a TWS trend along the Gulf of Alaska. These two regional studies suggest that the increasing trend we see over ice-free land in our model results is likely incorrect, particularly because the model does not account for true storage-enhancing processes (e.g., aquifer recharge, uptake into vegetation in newly deglaciated terrain) that would counteract the expected decreasing water balance from glacier ice loss. One explanation for the increase may be due to biases within our MicroMet-interpolated MERRA-2 input data, which may overproduce precipitation that is, not contributing to runoff. In particular, the model is likely generating excess, perennial snow over non-glacier high eleva-tion land cells, resulting in a positive water balance. The overproduction of snow can be linked to both (a) the overall positive precipitation biases, and (b) the cold biases we observe in air temperature fields versus those at the NOAA weather stations in Juneau and Skagway (see Section 3.2.2). This finding highlights the challenge of reproducing precipitation in mountain topography, particularly in high latitude maritime regions where snow can occur at all months at elevation.

Model Limitations
There are several sources of uncertainty within our model results. SnowModel-HydroFlow focuses largely on internal processes within the snowpack, but neglects several elements that may be important to glacier mass balance. In terms of processes that may contribute to additional ice melt, these include geothermal fluxes at the glacier ice/bed interface, as well as dynamical processes such as frictional melting from viscous heating (internal deformation of the ice) or sliding at the glacier bed (Mernild et al., 2014). Including these processes would require incorporating geothermal flux and ice dynamics components into the model, which is beyond the scope of this study on surface processes.
To prevent endless snow accumulation at high elevations from year to year, at the end of each summer SnowModel resets the snowpack to zero. While this prevents the formation of firn, which acts as a porous medium allowing for meltwater percolation and refreezing, snow also serves that role within SnowModel. The temporary loss of a porous medium is likely short-lived at these high elevation areas. However, this does not solve the issue of possible percolation beyond the summer surface into the firn (i.e., internal accumulation). Previous observational studies on the Juneau Icefield have identified the presence of internal accumulation within firn (Miller & Pelto, 1999;Pelto et al., 2013). However, the process is believed to occur only before mid-June (i.e., before isothermal conditions reign) (Miller & Pelto, 1999), and given the recent migration of transient snow lines to higher elevations (Pelto, 2019), the firn area on the icefield is overall being reduced. Our simulations are also tuned to reproduce what is thought to be the most accurate mean annual mass balance rate, suggesting that the correct amount of water equivalency still becomes available as runoff.
SnowModel also does not account for changes in glacier geometry resulting from climate forcing, either in terms of reduced area with glacier retreat, or lowered surface elevations with ice thinning. Rather, our simulations use a reference glacier surface representing conditions in the early 2010s, during which the highest-quality imagery was collected and incorporated into the National Elevation Data set (our DEM), and used to delineate the most accurate glacier outlines to date (Pfeffer et al., 2014). However, as this time period lies toward the end of our model period, it is likely that our icefield geometry is too low in elevation and too small in extent for the initial years of our simulation. The former would likely cause an overproduction of glacier ice melt and runoff due to higher temperatures at lower elevations, while the latter would cause an underproduction due to insufficient glacial extent. Quantifying each of these would require accurate DEMs and glacier outlines for our full model domain from the 1980s, which unfortunately do not exist. This presented a similar barrier to quantifying uncertainty in the Beamer et al. (2016) SnowModel study on Gulf of Alaska region runoff. In terms of glacier area, in their geodetic study, Berthier et al. (2018) employed a similar approach as ours and assumed a fixed area between 2000 and 2016, given a lack of additional outlines during that time period. The authors of that study instead included area uncertainty as one element in their overall uncertainty calculations. Because we leverage the Berthier et al. (2018) uncertainties in deriving our own model uncertainties, we believe our results in effect include this area uncertainty (e.g., in the 95% confidence intervals included in Table 4).

Glacier Change Present and Future
In a study by Ziemen et al. (2016) that modeled future mass loss scenarios for the Juneau Icefield to 2100, the authors found the only possible scenario leading to mass balance stabilization would be maintaining climate at the 1971 to 2010 level. Literature on current and future climate, however, suggests that such a climate scenario is highly unlikely. Several studies on Alaska glaciers have for example linked increasing glacier mass loss rates primarily to increases in summer air temperatures (Arendt et al., 2009;Criscitiello et al., 2010;O'Neel et al., 2014;Young et al., 2018), and indeed summer air temperatures are expected to increase as much as 5°C over northern high latitudes by 2100 (Koenigk et al., 2013). Maritime glaciers are also particularly sensitive to precipitation variations, and especially to decreasing amounts of snow deecting solar radiation (e.g., De Woul and Hock (2005)). Climate studies have found decreasing recent trends of annual snow precipitation volumes along the southeastern Alaska coast from 1979 to 2010 (Liston & Hiemstra, 2011), as well as into the future (to 2100), despite increasing precipitation overall (Shanley et al., 2015). Analysis of a downscaled gridded climate product also shows that Alaska is experiencing shifts in the rainsnow fraction towards rain (McAfee et al., 2014), a phenomenon to which coastal glaciers have been found to be especially sensitive (Moore et al., 2009).
Taken together, we see little evidence that a constant-climate scenario will occur in this region, given current and future trends in increasing air temperature and decreasing snow. As such, there is little indication that glacier mass loss acceleration in the western Juneau Icefield area will decrease or reverse. In fact, our 1980 to 2016 mass loss rate, being more negative than Ziemen et al. (2016) to begin with, may point to even stronger accelerations to 2100 than their anticipated five-fold mass loss rate increase. This could result in an even greater reduction in size than their estimated 63% volume loss and 62% area loss by 2100, an outcome that would substantially alter downstream hydrology.

Glaciological Linkage to Total Runoff
We find that mean annual total runoff from our coastal watershed domain is 20.0 km 3 a −1 for 1980 to 2016. On a seasonal basis, total runoff ranges from a minimum of 0.004 km 3 in February to a maximum of 5.0 km 3 in July (Figure 7). We observe a single peak in runoff in summer associated with glacier contributions and no secondary peak associated with spring snowmelt. This is consistent with Hill et al. (2015) who observed in a modeling study of 1960-2010 freshwater discharge a single peak in the hydrograph of the southern Gulf of Alaska region versus a dual peak in the north. Of the total runoff, 55% is sourced from glacier surfaces, a higher value than previous regional estimates for the Gulf of Alaska at 38%-47% (Beamer et al., 2016;Neal et al., 2010). The contribution of glacier volume loss to total runoff in our coastal domain is 12% for 1980 to 2016, as compared to regional Gulf of Alaska estimates of 7%-10% (Beamer et al., 2016;Hill et al., 2015;Neal et al., 2010). The larger glacier contributions here are likely due to the greater extent of ice cover in our domain (44%) relative to the larger Gulf of Alaska domain (∼17%).
Our results indicate that total annual runoff over the 36 year period of study is not correlated with annual glacier mass balance values. This shows that, in coastal environments, even large glaciers or icefields experiencing mass loss may not exert a strong control on total runoff given an overwhelming precipitation signal. This emphasizes the importance of not using annual mass balance values as a proxy for streamflow, and is supported by similar findings for another maritime Alaska glacier basin in O'Neel et al. (2014).
We also find that glacier runoff volumes are more strongly correlated with total runoff (r = 0.95, p < 0.001) than with glacier ice melt (r = 0.82, p < 0.001), suggesting that glacier runoff is more strongly controlled by overall precipitation events than by glacier ice melt. This decoupling between glacier ice melt and runoff is likely to be further enhanced in the future, given the projected change in rain/snow fraction toward rain (McAfee et al., 2014;Shanley et al., 2015), which is likely to contribute proportionally more to glacier runoff than to glacier ice melt. These findings highlight the dominant control that precipitation plays in this coastal glacierized environment.

Glacier Ice Melt and Glacier Runoff Trends Present and Future
Examining the annual volume of glacier ice melt over our study period, our results suggest a strongly increasing trend of nearly 10% per decade. Further evidence of increasing glacier ice melt rates is seen in the increasing amplitudes in Figure 11b in recent decades, as well as in the increasing anomalies toward the end of the study period in Figure 12. This finding indicates that in this high latitude maritime glacierized domain, the annual volume of glacier ice melt has not yet reached its maximum and will continue to increase to a yet unknown peak before it begins to decrease. This increasing signal is more difficult to detect (both in terms of magnitude as well as statistical metrics) in annual volumes of glacier runoff (+3% increase) and in total runoff (+1.4% increase). We expect this given proportionally increasing contributions in those two variables from precipitation, which is prone to high variability in this area (Bieniek et al., 2014). Nonetheless our findings of an increase in total runoff are consistent with an analysis of stream gauge records from the Wolverine Glacier, another maritime glacier watershed in Alaska that experienced a 23% increase in summer streamflow (i.e., a measure of total runoff) between 1966 and 2011 (O'Neel et al., 2014). While that study was based on gauge measurements and therefore lacked the ability to partition hydrological components, our modeling approach allows us to identify that glacier ice melt is most responsible for the increase in total runoff in our coastal glacierized domain.
Our findings, along with projections for strong and continued warming at high latitudes (Koenigk et al., 2013), lead us to expect that glacier runoff in the western Juneau Icefield will continue to increase before such time as the glaciers lose enough volume to reverse this trend. Although accurately predicting when this will occur would require coupling a hydrological routing model to glacier mass balance modeling projections such as those in Ziemen et al. (2016), which is beyond the scope of this hindcasting study, we speculate that given regional projections for the Gulf of Alaska of a peak water period near 2060 to 2070 (Huss & Hock, 2018), it will be several decades before the phenomenon occurs in our domain.

A Changing Hydrological Regime
Even with a strong increasing trend in annual glacier ice melt volumes, total runoff in this coastal glacierized area shows evidence of only a slightly increasing trend. Our findings instead reveal that the most prominent signs of hydrological regime change in this region are with respect to the timing and biogeochemical characteristics of the water being delivered downstream.
One indicator of these water quality changes is an increase in the magnitude of the maximum daily volume of glacier ice melt at a rate of 10% per decade. This increase has the potential, on those maximum flow days, to substantially modify freshwater conditions downstream as the proportion of glacier ice melt input grows relative to other freshwater sources. Additionally, although we do not detect robust trends in the onset, end, or subsequent length of the glacier ice melt season, our results suggest a marked increase in glacier ice melt delivery during the spring months, which in essence serves to shift periods of high glacier ice melt earlier into the year (Table 4, Figure 11). This earlier arrival signals a shift toward a hydrograph more closely resembling that of snowmelt-dominated basins. This finding is supported by regional analyses of temperature records in western North America over the past 50 years that show an asymmetry in warming of spring versus fall (Abatzoglou & Redmond, 2007). However, in projections to 2,100, Koenigk et al. (2013) found the most pronounced increases in air temperature in Alaska are likely to occur in winter and fall. We suggest, therefore, that there is potential for future increases in glacier ice melt and glacier runoff volumes in the fall season as well.
Several downstream impacts have occurred since the 1980s with a 16% increase per decade in springtime glacier ice melt and a corresponding 7% increase in glacier runoff. Given the tight relationship between stream temperature and glacier cover in this area (Fellman et al., 2014), our results suggest that stream temperatures during the spring months have likely become lower on account of the higher proportion of glacier ice melt input. In addition, we speculate there has been an increase in turbidity stemming from the influx of glacially eroded sediment along with increased glacier melt (Milner et al., 2017). Minerals and limiting nutrients contained therein are in turn likely delivered earlier and at larger magnitudes, including phosphorous, nitrogen, iron, and bioavailable organic carbon to riverine and estuarine food webs .
In rivers and streams, both temperature and water clarity are key variables for Pacific salmon spawning ground habitat selection (Lorenz & Filer, 1989), particularly given the sharp thermal limits of these species (Richter & Kolmes, 2005;Welch et al., 1998). Indeed, evidence is already mounting that populations among several Pacific salmon species are migrating to freshwater up to 0.5 days earlier per year than they did historically (Kovach et al., 2015). Although the mechanisms for the earlier timing remain complex, freshwater conditions in the riverine environment may contribute, given freshwater conditions that may support migration earlier in the year. For other populations, however, there is some concern that eventual decreased summer flows may lead to higher water temperatures and in turn lead to reduced salmonid function (Richter & Kolmes, 2005) as well as a reduction in spawning habitat (Wobus et al., 2015). These latter concerns may come to fruition after the period of peak water has passed in this domain.
In addition to altering stream conditions, the biogeophysical signature of glacier runoff also extends kilometers into Gulf of Alaska fjords, by setting up a stratified water column with fresh, cold, turbid, and generally nutrient-rich water at the ocean surface (Arimitsu et al., 2016). Therefore, changes in the timing of arrival of large volumes of glacier runoff will influence both estuary and stream conditions. In the estuary, glacially influenced environmental gradients explain much of the distribution and abundance of phytoplankton, which in turn drives higher trophic level food web structure for copepods, fish, and sea birds (Arimitsu et al., 2016).
Altogether, under continued warming and a decrease in precipitation as snow, projections continue to call for substantial and varied change to these and other hydroecological variables into the future (Shanley et al., 2015).

Conclusions
This study applied the coupled glacio-hydrological model SnowModel-HydroFlow to estimate daily freshwater runoff from 1980 to 2016 for the coastal watershed draining the western Juneau Icefield in Southeast Alaska, an area of 6,405 km 2 with 44% glacier cover. We find a strongly increasing trend in annual glacier ice melt production (9.6% per decade), with especially pronounced increases during spring months (16.5% per decade). This increase can also be detected in both glacier runoff (3.0% for annual volumes, 6.8% for spring volumes) and total runoff (1.4%, 2.7%). Together, these results suggest that this particular region has not yet passed the period of peak water associated with a persistent negative mass balance, likely on account of the extensive glacier coverage.
Unlike studies based on stream gauge data, our model results afford the opportunity to identify that glacier ice melt is the likely hydrological driver behind increases in total runoff seen over the past several decades. Moreover, our study contributes new and affirmative knowledge toward the question of whether glacier runoff trends can be detected in maritime climates with extreme precipitation amounts and high precipitation variability.
Overall in this domain, glacier runoff contributes 55% of total runoff, including 12% from non-renewable glacier volume loss. Total runoff in the domain is found not to be correlated to annual glacier mass balance, supporting the paradigm that advises against using annual balances as a proxy for glacier runoff volumes. Given projection studies that predict increasing glacier volume loss for the Juneau Icefield through 2,100, we anticipate ongoing glacier ice melt increases decades into the future, until such point as peak water is passed and the contributions of glacier ice melt and glacier runoff to the domain begin to reverse.
We find that changes in runoff timing and biogeochemical properties are the aspects of the hydrological regime undergoing the greatest changes in this coastal glacierized environment, with substantial impacts for downstream ecosystems. In particular, the earlier arrival of large volumes of glacier ice melt in spring is likely exerting an influence on stream temperature and clarity, a point of concern for downstream species such as salmon that have evolved to survive in particular freshwater conditions.
Ultimately, our results emphasize that even in maritime climates with high precipitation variability, high latitude glacierized watersheds are experiencing perceptible and ongoing hydrological regime change given persistent glacier volume loss.

Data Availability Statement
All previously published data used in this study are available from the sources mentioned within the Data section of this article. Model code for SnowModel-HydroFlow can be found at ftp://ftp.cira.colostate.edu/ ftp/Liston/JCYoung_WRR_2020/model_code/, and for the SoilBal routine at https://doi.org/10.4211/ hs.8e12debf926c4299acc782f9407512f5.