Hydrologic Regime Changes in a High-Latitude Glacierized Watershed under Future Climate Conditions

A calibrated conceptual glacio-hydrological monthly water balance model (MWBMglacier) was used to evaluate future changes in water partitioning in a high-latitude glacierized watershed in Southcentral Alaska under future climate conditions. The MWBMglacier was previously calibrated and evaluated against streamflow measurements, literature values of glacier mass balance change, and satellite-based observations of snow covered area, evapotranspiration, and total water storage. Output from five global climate models representing two future climate scenarios (RCP 4.5 and RCP 8.5) was used with the previously calibrated parameters to drive the MWBMglacier at 2 km spatial resolution. Relative to the historical period 1949–2009, precipitation will increase and air temperature in the mountains will be above freezing for an additional two months per year by mid-century which significantly impacts snow/rain partitioning and the generation of meltwater from snow and glaciers. Analysis of the period 1949–2099 reveals that numerous hydrologic regime shifts already occurred or are projected to occur in the study area including glacier accumulation area, snow covered area, and forest vulnerability. By the end of the century, Copper River discharge is projected to increase by 48%, driven by 21% more precipitation and 53% more glacial melt water (RCP 8.5) relative to the historical period (1949–2009).


Introduction
Cold regions, which occupy over half the land area of the Northern hemisphere, are particularly vulnerable to hydrologic regime shifts due to the sensitivity of snow and ice to warming temperatures [1].Cold regions are broadly defined as land areas with either average temperature of the coldest month below freezing, annual peak snow accumulation of at least 0.3 m, ice that prevents navigation for 100 days, or 0.3 m of frost penetration once per decade [2].A regime shift occurs when the cumulative effect of small perturbations cause an abrupt change in system behavior [3] and the point at which a regime shift occurs may not be readily identifiable until after the fact.However, leading indicators are sometimes apparent, such as an increase in variability, or a delayed return to equilibrium following disturbance [3].
As mountainous cold regions warm, potential regime shifts include the loss of snow and glaciers, permafrost thaw, and transitioning from snow-to rain-dominated systems.The depth and extent of snow cover govern where and to what extent glaciers grow, where glacier surfaces are exposed, and the magnitude of spring runoff from snow melt.As snow melt ceases, the exposed glacier ice melts, generating streamflow throughout the summer until temperatures drop in the early fall.Many high-latitude mountainous watersheds are projected to generate more streamflow in the coming decades as glacier ablation increases in a warming climate; however, as glacier mass is exhausted, this surge in glacier runoff will be reversed [4][5][6].
Glacier regime shifts are often associated with hydrologic regime shifts such as river avulsion [7], sea level rise [8], changes to the timing and quantity of streamflow [1], and increased frequency and severity of flooding and droughts [9].Other climate-related hydrologic changes that are independent of glacier changes may also lead to hydrologic regime shifts in high latitude glacierized watersheds.As permafrost thaws, the melt water contributions to streamflow are relatively small, but altered connectivity between surface water and groundwater can spontaneously drain wetlands and surface features, or increase groundwater contributions to streamflow [10,11].As wetlands drain and soils dry, ecological feedbacks with hydrological implications may occur such as forest disturbance [12], vegetation change [13], and increased or decreased transpiration [14].
Advance knowledge of hydrologic regime change is important for planning and protection of aquatic habitat, recreation, water supply, erosion control, wildfire management, and infrastructure.Where remote settings and rugged terrain impede the establishment of robust hydro-climatic monitoring networks, parsimonious hydrologic models serve an important role to understand how a watershed responds to perturbations.The current research utilizes glacio-hydrological simulations to explore how gradual changes in temperature and precipitation may trigger regime shifts in a heavily-glacierized mountainous watershed in Southcentral Alaska.Using a calibrated monthly water balance model, MWBMglacier, we simulate the water cycle during historical and future time periods to test two hypotheses: (i) that the prevalence of hot and dry conditions will increase the duration of time the forests are vulnerable to disturbance such as wildfire; and (ii) that hydrologic regime shifts associated with glacio-hydro-climatic interactions will be detectible during the 1949-2099 glacio-hydrologic simulations.

Materials and Methods
The materials and methods employed in this study are presented in the following four subsections: study area, data, glacio-hydrologic model, and regime change identification.

Study Area
The Copper River originates from hundreds of glaciers situated in the Alaska Range, Chugach Mountains, Wrangell Mountains and St. Elias Mountains (Figure 1).The river flows 470 km southward, through a break in the Chugach Mountains, to form a broad delta at the confluence with the Gulf of Alaska.Million Dollar Bridge marks the outlet of the 60,800 km 2 watershed.Elevations range from near sea level at the watershed outlet to 5832 m.The Copper River Watershed (CRW) comprises four ecoregions [15].The Chugach-St Elias ecoregion, occupying 48% of the CRW to the south and east, is characterized by steep, rugged mountains, shallow soil, discontinuous permafrost, and a temperate coastal climate with very high precipitation (up to 7000 mm/year) [15,16].The Alaska Range and Wrangell Mountain ecoregions, occupying 10% and 15%, respectively, of the CRW have a cold and dry interior climate, steep and rugged mountains, shallow soil, and discontinuous permafrost [15].The cold and dry Copper Basin ecoregion, referred to herein as the Forested River Valley, is a wetland complex occupying 27% of the CRW and is characterized by relatively flat terrain with discontinuous permafrost and associated ponds and thaw lakes [15,16].

Data
Data sources used in the current study are summarized in Table 1.SNAP data sets of historical and future climate were selected because they are the only products currently available that offer historical and future monthly climate data processed using the same methods, with the required spatial extent (Alaska and Canada), the desired fine spatial resolution (2 km grid cells), and the required period of record (1950 to 2100).SNAP future climate data sets represent four climate scenarios, Representative Concentration Pathways (RCP) 2.6, 4.5, 6.0, and 8.5 [17].This study used RCP 8.5 and RCP 4.5.RCP 8.5 represents a "business as usual" scenario of climatic change and provides a plausible view of Copper River hydrology if emissions continue at the current pace.RCP 4.5 represents a scenario in which greenhouse gases increase until mid-Century and then stabilize.
Table 1.Sources of data.

Information
Use Source

Data
Data sources used in the current study are summarized in Table 1.SNAP data sets of historical and future climate were selected because they are the only products currently available that offer historical and future monthly climate data processed using the same methods, with the required spatial extent (Alaska and Canada), the desired fine spatial resolution (2 km grid cells), and the required period of record (1950 to 2100).SNAP future climate data sets represent four climate scenarios, Representative Concentration Pathways (RCP) 2.6, 4.5, 6.0, and 8.5 [17].This study used RCP 8.5 and RCP 4.5.RCP 8.5 represents a "business as usual" scenario of climatic change and provides a plausible view of Copper River hydrology if emissions continue at the current pace.RCP 4.5 represents a scenario in which greenhouse gases increase until mid-Century and then stabilize.For each climate scenario, five models from the Climate Model Intercomparison Project Phase 5 (CMIP5) [25] are available from SNAP, along with a five-model mean.The five CMIP5 models in the SNAP ensemble were selected from a group of 22 using the methods described in Walsh et al. [26].SNAP downscaled the five CMIP5 models from variable coarse resolutions to a uniform 2 km × 2 km resolution using the delta method (the five CMIP5 models have original spatial resolutions ranging from 6760 to 38,613 km 2 at the center of the study area, 62 N-145 W using https://www.nhc.noaa.gov/gccalc.shtmlaccessed 10 June 2017).Bias correction was performed by SNAP using PRISM data (http://www.prism.oregonstate.edu/)as the reference climate.Considerable biases can be introduced by applying different downscaling methods [27,28].The SNAP future data sets span years 2006 through 2099 and a companion data set of historical temperature and precipitation covers the period 1901-2009.Regime change analysis was performed using composite simulation results for both the historical data period (1949 to 2009), and simulations using GCM data for 2010 to 2099.
Eight additional climate data sets were constructed for use in this study (Table 2).Historical mean monthly climate data sets were constructed for the period 2006 through 2099 using mean monthly temperature and precipitation from the historical period 1949 through 2009 (scenario HTP).To isolate the effect of temperature change alone, two additional data sets were constructed, using modeled future temperature from the two five-model means (RCP4.5 and 8.5) along with mean monthly precipitation from the historical period (MM45THP and MM85THP).In an effort to identify the temperature rise at which a regime change would be detectible when considering all glaciers in the CRW as a whole, a temperature sensitivity analysis was performed.As part of the temperature sensitivity analysis, historical precipitation for the period 1949-2009 was paired with historical temperature for the same period, but the temperature values were increased by 2, 4, 6, 8, and 10 • C (scenarios S2 through S10).The temperature and precipitation data sets described above provided the forcing data for 21 different model simulations (Table 2).Precipitation adjustments were not included in this sensitivity analysis because all CMIP5 models project precipitation to increase, which mitigates, to varying degrees, the amount of glacier ablation caused by temperature rise.

Glacio-Hydrological Modeling
The conceptual water balance model selected for use in this study, the MWBMglacier [29], represents the partitioning of water among snow, glaciers, runoff, soil, and the atmosphere (Figure 2).

Glacio-Hydrological Modeling
The conceptual water balance model selected for use in this study, the MWBMglacier [29], represents the partitioning of water among snow, glaciers, runoff, soil, and the atmosphere (Figure 2).The MWBMglacier performs monthly accounting of water storage and fluxes as represented in Equations ( 1) and (2).
The model accounts for water stored in snow (Snowstor), glaciers (Glstor), soil (Prestor), and Delay, which temporarily stores runoff in soil or glaciers for delayed release.The only input flux of water to the system is precipitation (P), while output fluxes include runoff (RO), and transfers to the atmosphere in the form of evapotranspiration (ET) and sublimation (S), collectively referred to here as ETS.
The following overview of the MWBMglacier calculations [29] reference parameters that are described in Table 3. Model response units (MRUs) are discrete spatial regions of any size or shape The MWBMglacier performs monthly accounting of water storage and fluxes as represented in Equations ( 1) and (2).
The model accounts for water stored in snow (Snowstor), glaciers (Glstor), soil (Prestor), and Delay, which temporarily stores runoff in soil or glaciers for delayed release.The only input flux of water to the system is precipitation (P), while output fluxes include runoff (RO), and transfers to the atmosphere in the form of evapotranspiration (ET) and sublimation (S), collectively referred to here as ETS.
The following overview of the MWBMglacier calculations [29] reference parameters that are described in Table 3. Model response units (MRUs) are discrete spatial regions of any size or shape presumed to have homogeneous hydromorphologic characteristics.In the Copper River application, MRUs are uniform 2 km × 2 km grid cells that correspond to the resolution of the SNAP climate data.Steps 1 through 9, below, summarize model calculations performed monthly for each MRU, and Step 10 is performed at 10-year intervals.
Step 1: Snow and Rain Partitioning At air temperatures below TSNOW, all P is snow, and above TRAIN, all P is rain; between TSNOW and TRAIN, linear interpolation determines the mix of snow and rain.This is represented by partitioning step t in Figure 2.
Step 2: Snow Accumulation and Ablation Fresh snowfall is added to the snowpack after accounting for sublimation.The MWBMglacier uses a single parameter, SUBFAC, to account for sublimation from blowing snow, canopy sublimation, and sublimation from the snow surface.A fraction of new snowfall is partitioned into sublimation by multiplying new snowfall by the parameter SUBFAC (partitioning step S F in Figure 2).When current air temperature exceeds the melt temperature (TMELT), potential snow melt is calculated using a temperature-index degree-day melt rate algorithm where the parameter SMELTR defines the maximum rate of snow melt in units of mm per degree-day.SMELTR is adjusted monthly by the ratio of daylight hours in the current month to daylight hours in June; monthly day length is calculated using latitude [30].Degree days are the daily positive differences between current air temperature and TMELT.Actual snowmelt is the lesser of potential snowmelt or the amount of snow in the snow pack.The melt deficit is the difference between potential snowmelt and actual snowmelt.
Step 3: Direct Runoff Direct runoff (DRO) is surface water from precipitation that is unable to infiltrate (e.g., due to impervious surfaces) and immediately becomes runoff; it reduces the amount of rain and melting snow that is able to infiltrate.In the MWBMglacier, the parameter DROFAC is multiplied by the sum of rainfall and melting snow to produce direct runoff for every MRU (partitioning step Df in Figure 2).In the Copper River application, permafrost MRUs are assigned a higher value of DROFAC than permeable soil (0.4 for permafrost and 0.01 to 0.1 for permeable soil).
Step 4: Evapotranspiration Total available soil water (ASW) is the sum of infiltrating water from precipitation plus pre-existing soil moisture that is not bound to the soil.The percent soil saturation in the prior month influences the amount of water considered to be bound to the soil (and unavailable for ET).Actual evapotranspiration (ET) is the lesser of potential evapotranspiration (PET) or ASW.PET is calculated for non-glacierized regions using the Hamon equation [2].In the Copper River application, the Hamon coefficient was calibrated to Penman Monteith PET calculated using the Valiantzas equations [3,4] with CFSR data.
Step 5: Soil Saturation Runoff If, after accounting for ET, the total amount of water in soil exceeds the soil water holding capacity (AWC), soil is considered to be saturated and any excess water becomes soil runoff (SRO).
Step 6: Glacier Accumulation and Ablation Glaciers gain mass when snow accumulates from year to year.Any snow that remains on a glacier at the end of September is transferred from snow storage (Snowstor) to glacier storage (Glstor).Glacier runoff (GRO) consists of glacier base flow (GRO B ) plus glacier surface runoff (GRO S ).GRO B is ice melt discharged from the base of the glacier that occurs every month, regardless of temperature, as specified by the parameter GBFLOW.GRO S is the sum of ice melt from the exposed glacier surface plus any rain or snowmelt on the glacier that was not diverted to direct runoff in Step 3. If, after Step 2, a melt deficit exists on a glacierized MRU, either because no snow was present or because not enough snow was present, the melt deficit is satisfied by melting exposed glacier ice.Ice melt from the glacier surface is calculated using a degree-day equation as for snow, but with a glacier-specific rate of melt, GMELTR.GMELTR is adjusted for day length as described for SMELTR.The ice within any glacierized MRU is assumed to be of unlimited thickness, and the model tracks the change in thickness each time step to determine change in glacier mass.The simulation results are not sensitive to the assumption of unlimited ice thickness because the model relies upon change in thickness, rather than absolute thickness, to evaluate changes in glacier mass balance and changes in glacier area.When calculations of glacier area change dictate the conversion of a glacierized MRU to a non-glacierized MRU (discussed in Step 10, below), the ice thickness remains greater than zero but the MRU's change in glacier mass will be zero.Sublimation from the surface of ice is not simulated, which could result in underestimation of glacier thickness change equivalent in extent to the transfer of water from glacier ice to the atmosphere.Ice sublimation will be explored in subsequent phases of model development.
Step 7: Storage and Delayed Release of Runoff from Soil and Glaciers The function of the Delay conceptual reservoir is to simulate the lag time as water travels overland and through porous media.Each month, newly generated runoff from SRO and GRO s is combined with water previously stored, and Delay is updated as shown in Equation (3).
where m and m − 1 represent the current month, and the prior month, respectively.
Step 8: Total Runoff Total runoff (totRO) in the current month is calculated using Equation (4).
Step 9: Water Balance Calculation and Model Output Each month, the model produces a breakdown of the water balance into the following components: Transfers to the atmosphere in the form of sublimation and evapotranspiration.
In addition, glacier balanced-budget equilibrium line altitude (ELA) is calculated using median glacier elevation as a proxy [31].Accumulation area ratio (AAR) is calculated annually as the ratio of glacierized area with snow present in September to the total glacierized area.
Step 10: Change in Glacier Area The volume-area scaling method [32] is used to estimate the change in glacier area associated with simulated change in glacier volume every ten years to generate smooth transitions as the glacier equilibrates [33].Change in glacier volume (dV) is calculated as area (A) times change in thickness, and the corresponding change in glacier area (dA) is calculating using the equation, with 1.375 for the exponent constant, γ, and 0.034 km (3−2γ) for the scaling parameter, c [32,34]; If an individual glacier comprises multiple MRUs, volume-area scaling calculations are performed on that glacier's total area and volume rather than on individual MRUs.The model utilizes ten parameters, monthly time series of temperature and precipitation, and five types of input data (described in Table 3) for each model response unit (MRU).The model was previously calibrated and evaluated to historical observations in the CRW [29].In this study, additional model performance evaluations were performed to ensure that GCMs adequately represent actual hydroclimatic conditions.Using the first 11 years of simulation, 2006 through 2016, simulated runoff was compared to measured streamflow at the watershed outlet, and simulated change in storage was compared to estimates derived from the Gravity Recovery and Climate Experiment (GRACE) satellite observations.Mean monthly data over the 2006-2016 period were used for evaluation because future climate data are not intended to reproduce historical weather on a daily or monthly basis, but are expected to have similar statistics.In contrast, during calibration and evaluation of the model using historical observed climate data, monthly time series were used.

Identification of Hydrologic Regime Changes
Automated analytical methods were employed to characterize the MWBMglacier output and detect when changes in time series data support the presence of an abrupt change or regime shift marked by a tipping point [35].Automated visualization of abrupt change in hydrologic signals over time permitted the identification of one or more points in time when the slope or mean of the data plotted over time changed abruptly; in this study, a minimum of 20 years was required between change points [36].Potential regime shifts detected through the visualization of signal change were further evaluating using the cumulative sum, a test that identifies when a time series mean drifts beyond pre-specified boundaries, defined here as five standard deviations above or below a target mean [37]; unless noted otherwise, the target mean and standard deviation were calculating using the first 30 years of data in the time series.The non-parametric Wilcoxon Rank Sum test was also used to test the null hypothesis that 2 samples of data from different time periods were from the same population; a finding of p < 0.05 reinforced the identification of a regime shift.In this study, all references to a "significant trend" with time are based on the nonparametric Mann Kendall test [38] yielding p < 0.05.Mann Kendall and Wilcoxon Rank Sum tests were performed only on time series data lacking a seasonal signal, such as annual totals or an annual time series of data from a single month.
As glaciers recede, the quantity of surface melt water first increases, then decreases depending on the balance between increasing rates of ice melt and decreasing glacier area [39].To identify when a glacier's melt rate is no longer sufficient to compensate for loss of glacier area, we calculated the original rate (S 1 ) of glacier thickness change over a ten-year period (H 1 ), along with the original glacier area (A 1 ), and the new area (A 2 ).To generate the same volume (V) of melt water with the new glacier area, where the new thickness change (H 2 ) required over the subsequent ten years is calculated by rearranging Equation ( 7).

Results
Results are presented in two sections: glacio-hydrological modeling of the CRW using a calibrated and evaluated water balance model during historical and future time periods [40], followed by the identification of regime changes.

Evaluation of Model Performance 2006-2016
The model was previously calibrated to historical observations over the time period 1949 to 2009 by Valentin et al. [29].In that study, simulated streamflow was calibrated and evaluated against measured streamflow at 8 USGS stream gauges identified in Figure 1, simulated changes in glacier mass balance were calibrated to data in the published literature, and simulated snow-covered area (SCA) was calibrated to MODIS satellite-based observations.Further evaluations were performed by comparing simulated ET to MODIS satellite-based observations of ET and by comparing simulated changes in monthly total water storage (∆TWS) to estimates from GRACE satellite data.Sources of data used for calibration and validation are listed in Table 1.Acceptable calibration for time series data was defined as Nash Sutcliffe Efficiency (NSE) greater than 0.7, normalized root mean square error (RSR) less than 0.7, and percent bias (PBIAS) within 25% of zero [41].Calibration and evaluation results from the prior study [29] are summarized below.

•
Streamflow calibration and evaluation was acceptable at the six stream gauges whose catchment areas included glaciers (Million Dollar, Tonsina, Klutina, Tazlina, and Gakona) with NSE ranging from 0.79 to 0.89, RSR from 0.33 to 0.45, and PBIAS from 1.4% to 29% at those gauges.Streamflow calibration was not acceptable at the two unglacierized tributaries (Gulkana and Sinona).

•
Calibration of SCA to MODIS satellite-based observations (MOD10CM) was acceptable in all four ecoregions (NSE 0.89-0.93),but when the results were further evaluated by elevation band, calibration was unacceptable in the two regions where elevations exceeded 2500 m.

•
Twenty-three comparisons were made between simulated mass balance change (∆mb) at 12 individual glaciers during specific time periods where data were available in the published literature (reference data).The mean ∆mb for all reference data was −0.64 ± 0.36 m water equivalent per year (mwe/year) and corresponding simulated results were −0.72 ± 0.35 mwe/year.

•
A comparison of monthly simulated ET versus MODIS ET resulted in NSE of 0.68 with an average annual difference of 121 mm/year, with MOD16 observations generally higher than simulated ET.

•
The annual difference between GRACE and simulated ∆TWS was 50 mm or less on a monthly basis (well within the reported error associated with the GRACE data), and, on an annual basis, GRACE ∆TWS was on average 54 mm higher than simulated ∆TWS.The mean annual ∆TWS deviation represents 4% of mean annual precipitation flux.
In the present study, CMIP5 models of future climate overlapped with 11 years of the historical record which permitted additional evaluation of model performance.CMIP5 models of the historical time period are not expected to reproduce daily or monthly conditions, but CMIP5 models should be able to reproduce general statistics, such as mean monthly values, from the observed data.The first evaluation performed in this study compared simulated streamflow to measured streamflow at the watershed outlet using mean monthly data calculated for the period 2006-2016 (Figure 3a).The second evaluation compared simulated changes in total water storage (TWS) to estimates of TWS based on the GRACE satellite-based measurements of earth's gravity field (Figure 3b).Simulations of RCP scenarios MM45 and MM85 produced almost identical mean monthly results, and both compared favorably with observed runoff with the exception that the onset of simulated runoff in the spring occurs one month too early.Evaluation of total water storage against GRACE observations were acceptable and well within the reported range of uncertainty (Figure 3b).The mean annual change in TWS is −86 mm per year based on GRACE data and −61 mm per year based on the MWBMglacier simulations.
Water 2017, 9, x FOR PEER REVIEW 10 of 23 In the present study, CMIP5 models of future climate overlapped with 11 years of the historical record which permitted additional evaluation of model performance.CMIP5 models of the historical time period are not expected to reproduce daily or monthly conditions, but CMIP5 models should be able to reproduce general statistics, such as mean monthly values, from the observed data.The first evaluation performed in this study compared simulated streamflow to measured streamflow at the watershed outlet using mean monthly data calculated for the period 2006-2016 (Figure 3a).The second evaluation compared simulated changes in total water storage (TWS) to estimates of TWS based on the GRACE satellite-based measurements of earth's gravity field (Figure 3b).Simulations of RCP scenarios MM45 and MM85 produced almost identical mean monthly results, and both compared favorably with observed runoff with the exception that the onset of simulated runoff in the spring occurs one month too early.Evaluation of total water storage against GRACE observations were acceptable and well within the reported range of uncertainty (Figure 3b).The mean annual change in TWS is −86 mm per year based on GRACE data and −61 mm per year based on the MWBMglacier simulations.1) are illustrated by gray shading in Figure 3a.Error in GRACE estimates (illustrated by gray shading in Figure 3b) accompanied the data provided by the University of Colorado GRACE Data Portal (listed in Table 1).The reported error associated with individual monthly estimates of total water storage was doubled because the change in monthly storage was calculated for illustration in Figure 3b.

Summary of Future Climate
At the end of the century, scenario RCP 8.5 indicates that temperatures will be above freezing for two additional months each year in all ecoregions and at all elevations in the study area relative to the present time.Recent work by Young et al. (2017) [12] determined environmental factors linked to wildfire occurrence; one of the factors identified for the CRW Forested River Valley was temperature in excess of 13.4 °C.In the Forested River Valley, temperatures will be above 13.4 °C for three months by the end of the century (RCP 8.5 model mean) compared to only one month during the historical period .By the end of the century, air temperatures will be consistently  1) are illustrated by gray shading in Figure 3a.Error in GRACE estimates (illustrated by gray shading in Figure 3b) accompanied the data provided by the University of Colorado GRACE Data Portal (listed in Table 1).The reported error associated with individual monthly estimates of total water storage was doubled because the change in monthly storage was calculated for illustration in Figure 3b.

Summary of Future Climate
At the end of the century, scenario RCP 8.5 indicates that temperatures will be above freezing for two additional months each year in all ecoregions and at all elevations in the study area relative to the present time.Recent work by Young et al. (2017) [12] determined environmental factors linked to wildfire occurrence; one of the factors identified for the CRW Forested River Valley was temperature in excess of 13.4 • C. In the Forested River Valley, temperatures will be above 13.4 • C for three months by the end of the century (RCP 8.5 model mean) compared to only one month during the historical period .By the end of the century, air temperatures will be consistently above freezing for two additional months per year relative to current conditions.Precipitation is projected to rise uniformly in all regions of the study area for both RCP 4.5 and RCP 8.5.

Snow
Although the CRW will experience more winter precipitation in future years, less of it will fall as snow due to warming temperatures, resulting in a statistically significant decline in annual snowfall.The snow line, simulated as the minimum elevation where snow persists year-round, rises at a significant rate in the Alaska Range (7 m/year), the Wrangell Mountains (7 m/year) and the Chugach-St Elias Mountains (4 m/year).April snow covered area will decrease by more than 50%, from 25% to 12% in the Wrangell Mountains, and from 30% to 13% in the Chugach and St. Elias Mountains, with the decline commencing by the 2020s in both regions.Presently, snow covered area in April is near zero in the Alaska Range and the forested river valley and will remain so in the future.above freezing for two additional months per year relative to current conditions.Precipitation is projected to rise uniformly in all regions of the study area for both RCP 4.5 and RCP 8.5.

Snow
Although the CRW will experience more winter precipitation in future years, less of it will fall as snow due to warming temperatures, resulting in a statistically significant decline in annual snowfall.The snow line, simulated as the minimum elevation where snow persists year-round, rises at a significant rate in the Alaska Range (7 m/year), the Wrangell Mountains (7 m/year) and the Chugach-St Elias Mountains (4 m/year).April snow covered area will decrease by more than 50%, from 25% to 12% in the Wrangell Mountains, and from 30% to 13% in the Chugach and St. Elias Mountains, with the decline commencing by the 2020s in both regions.Presently, snow covered area in April is near zero in the Alaska Range and the forested river valley and will remain so in the future.Rising median glacier elevation highlights that CRW glaciers are out of balance with their environment; a stable median glacier elevation would indicate a glacier in equilibrium with its environment [31].The simulated median glacier elevation rose 47 m (0.8 m/year) during the historical period (1949 to 2009) but it is projected to rise another 111 to 170 m (1.2 to 1.9 m/year) between 2010 and 2100 based on RCP 4.5 and RCP 8.5 five-model means, respectively.The accumulation zone (area where snow accumulates from year to year) increased from roughly 50 to 60% from 1940 through 1980, then began a period of steady decline that is projected to continue and approach 25% by the end of the century.Increased winter precipitation mitigates, to some extent, the temperature-related acceleration of glacier loss because snow insulates the surface of the glacier, protecting it from surface melt.However, under the examined scenarios of temperature and precipitation, glacier area is projected to decline by an additional 10% ± 5% (RCP 4.5) to 15% ± 5% (RCP 8.5) by the end of this century (Figure 5).Rising median glacier elevation highlights that CRW glaciers are out of balance with their environment; a stable median glacier elevation would indicate a glacier in equilibrium with its environment [31].The simulated median glacier elevation rose 47 m (0.8 m/year) during the historical period (1949 to 2009) but it is projected to rise another 111 to 170 m (1.2 to 1.9 m/year) between 2010 and 2100 based on RCP 4.5 and RCP 8.5 five-model means, respectively.The accumulation zone (area where snow accumulates from year to year) increased from roughly 50 to 60% from 1940 through 1980, then began a period of steady decline that is projected to continue and approach 25% by the end of the century.Increased winter precipitation mitigates, to some extent, the temperature-related acceleration of glacier loss because snow insulates the surface of the glacier, protecting it from surface melt.However, under the examined scenarios of temperature and precipitation, glacier area is projected to decline by an additional 10 ± 5% (RCP 4.5) to 15 ± 5% (RCP 8.5) by the end of this century (Figure 5).The difference between the solid red line and the solid black line shows the mitigating effect of high winter precipitation on glacier loss.Gray shading illustrates the standard deviation of the five CMIP5 models.

Future Runoff
The cold and dry interior portion of the CRW (the Alaska Range, Wrangell Mountains, and the Forested River Valley) covers half the land area of the CRW, but generates only 20% of the runoff.The Chugach and St Elias Mountains, having a coastal temperate rainforest climate, generate 80% of CRW runoff.In the future, the interior will experience decreasing streamflow in January and February, increasing streamflow in the spring (March, April and May) as snow melt occurs earlier in the year, no streamflow change in June and July, and increasing streamflow in August through December (Figure 6).The coastal mountains will generate increased streamflow in all months (Figure 7).The magnitude of peak annual monthly streamflow is projected to increase in the Wrangell, Chugach and St-Elias regions, but will decrease in the Alaska Range.Monthly flows generated in glacierized mountain regions peak around July, and in un-glacierized snow-dominated region, peak flows occur between April and May.In the future, variability in the timing of peak monthly flow occurs as snow melts earlier, and glacier melt increases during the summer.Based on the temperature sensitivity analysis, for each degree of temperature rise, total runoff increases 68 mm/year.The difference between the solid red line and the solid black line shows the mitigating effect of high winter precipitation on glacier loss.Gray shading illustrates the standard deviation of the five CMIP5 models.

Future Runoff
The cold and dry interior portion of the CRW (the Alaska Range, Wrangell Mountains, and the Forested River Valley) covers half the land area of the CRW, but generates only 20% of the runoff.The Chugach and St Elias Mountains, having a coastal temperate rainforest climate, generate 80% of CRW runoff.In the future, the interior will experience decreasing streamflow in January and February, increasing streamflow in the spring (March, April and May) as snow melt occurs earlier in the year, no streamflow change in June and July, and increasing streamflow in August through December (Figure 6).The coastal mountains will generate increased streamflow in all months (Figure 7).The magnitude of peak annual monthly streamflow is projected to increase in the Wrangell, Chugach and St-Elias regions, but will decrease in the Alaska Range.Monthly flows generated in glacierized mountain regions peak around July, and in un-glacierized snow-dominated region, peak flows occur between April and May.In the future, variability in the timing of peak monthly flow occurs as snow melts earlier, and glacier melt increases during the summer.Based on the temperature sensitivity analysis, for each degree of temperature rise, total runoff increases 68 mm/year.The difference between the solid red line and the solid black line shows the mitigating effect of high winter precipitation on glacier loss.Gray shading illustrates the standard deviation of the five CMIP5 models.

Future Runoff
The cold and dry interior portion of the CRW (the Alaska Range, Wrangell Mountains, and the Forested River Valley) covers half the land area of the CRW, but generates only 20% of the runoff.The Chugach and St Elias Mountains, having a coastal temperate rainforest climate, generate 80% of CRW runoff.In the future, the interior will experience decreasing streamflow in January and February, increasing streamflow in the spring (March, April and May) as snow melt occurs earlier in the year, no streamflow change in June and July, and increasing streamflow in August through December (Figure 6).The coastal mountains will generate increased streamflow in all months (Figure 7).The magnitude of peak annual monthly streamflow is projected to increase in the Wrangell, Chugach and St-Elias regions, but will decrease in the Alaska Range.Monthly flows generated in glacierized mountain regions peak around July, and in un-glacierized snow-dominated region, peak flows occur between April and May.In the future, variability in the timing of peak monthly flow occurs as snow melts earlier, and glacier melt increases during the summer.Based on the temperature sensitivity analysis, for each degree of temperature rise, total runoff increases 68 mm/year.

ET Deficit and Forest Vulnerability
Transfers to the atmosphere are projected to increase 4% to 12% (RCP 4.5 and 8.5, respectively), resulting from an increase in ET as temperatures rise and a decrease in sublimation as snowfall decreases over the coming decades (Table 4).The ET deficit (PET-ET) increases by 37% in the Forested River Valley.Forest vulnerability is identified here as conditions with both high temperature (>13.4 °C) and a positive moisture deficit [12].Historically, these hot and dry forest conditions were prevalent only during July, but, starting mid-century, most of the forested area will be hot and dry for three months (June, July and August), tripling the duration of forest stress.All water budget fluxes (precipitation, transfers to the atmosphere, and streamflow) increase at statistically significant rates, as well as contributions to streamflow from melting glacier ice (Figure 8).Increasingly negative changes in storage (illustrated for RCP 8.5 by the black dots in Figure 8) reflect the transfer of water from perennial storage in snow and ice to runoff and the atmosphere.

ET Deficit and Forest Vulnerability
Transfers to the atmosphere are projected to increase 4% to 12% (RCP 4.5 and 8.5, respectively), resulting from an increase in ET as temperatures rise and a decrease in sublimation as snowfall decreases over the coming decades (Table 4).The ET deficit (PET-ET) increases by 37% in the Forested River Valley.Forest vulnerability is identified here as conditions with both high temperature (>13.4 • C) and a positive moisture deficit [12].Historically, these hot and dry forest conditions were prevalent only during July, but, starting mid-century, most of the forested area will be hot and dry for three months (June, July and August), tripling the duration of forest stress.All water budget fluxes (precipitation, transfers to the atmosphere, and streamflow) increase at statistically significant rates, as well as contributions to streamflow from melting glacier ice (Figure 8).Increasingly negative changes in storage (illustrated for RCP 8.5 by the black dots in Figure 8) reflect the transfer of water from perennial storage in snow and ice to runoff and the atmosphere.Trend analysis of annual components of the water budget shows a steady acceleration in intensity of the water cycle and statistically significant increases in annual precipitation, runoff, transfers to the atmosphere, and melting glacier ice.Between the present decade (2010-2019) and the last decade of the century total runoff will increase by 17% (RCP 4.5) to 48% (RCP 8.5) but the fraction of runoff generated by melting glacier ice will remain around 20-21% (Table 4).Trend analysis of annual components of the water budget shows a steady acceleration in intensity of the water cycle and statistically significant increases in annual precipitation, runoff, transfers to the atmosphere, and melting glacier ice.Between the present decade (2010-2019) and the last decade of the century total runoff will increase by 17% (RCP 4.5) to 48% (RCP 8.5) but the fraction of runoff generated by melting glacier ice will remain around 20-21% (Table 4).

Regime Change
A progression of regime changes was detected, beginning with temperature then followed by snow, glacier contributions to runoff, the timing of peak runoff, and ending with forest vulnerability (Table 5).Change point visualization revealed distinct trends in temperature (Figure 9a): a cooling phase coincident with the 30-year PDO cool phase that occurred in Alaska from the mid-1940s to the mid-1970s, followed by a multi-decade PDO warm phase, and a subsequent future period of accelerated warming.The cumulative sum test (Figure 9b) indicated that mean temperature fell below "normal" (defined here as five times the 1950-2000 standard deviation) in the early 1970s, at the tail end of the PCO cool phase, and it shows that starting around 2020, temperatures will consistently exceed "normal".It is important to note that the future temperature projections for both RCP 4.5 and RCP 8.5 demonstrate less annual and multi-decadal variability than the historical period (for example, see Figure 9a).

Decade Description of Regime Change for RCP 8.5, 5-Model Mean 1970s
A 30-year period of increasingly cool temperatures associated with the Pacific Decadal Oscillation (PDO) transitions to a 30-year warm PDO warm phase.1980s Glacier AAR begins to decline.
2010s April snow covered area in the Wrangell, Chugach and St Elias Mountains peaks in the 2010s near 30% then begins a long, steady, decline ending the century near 10%.2020s In the mid-2020s, hot and dry forests conditions become prevalent in the month of June.2040s In the 2040s, glacier ice contributions to streamflow in the Tonsina River stops accelerating and begins to decline due to loss of critical glacier mass

Regime Change
A progression of regime changes was detected, beginning with temperature then followed by snow, glacier contributions to runoff, the timing of peak runoff, and ending with forest vulnerability (Table 5).Change point visualization revealed distinct trends in temperature (Figure 9a): a cooling phase coincident with the 30-year PDO cool phase that occurred in Alaska from the mid-1940s to the mid-1970s, followed by a multi-decade PDO warm phase, and a subsequent future period of accelerated warming.The cumulative sum test (Figure 9b) indicated that mean temperature fell below "normal" (defined here as five times the 1950-2000 standard deviation) in the early 1970s, at the tail end of the PCO cool phase, and it shows that starting around 2020, temperatures will consistently exceed "normal".It is important to note that the future temperature projections for both RCP 4.5 and RCP 8.5 demonstrate less annual and multi-decadal variability than the historical period (for example, see Figure 9a).The AAR represents the glacier's capacity to grow by accumulating snow from year to year.The global balanced-budget AAR (where mass gain equals mass loss) is 58% [5] and AAR for CRW glaciers examined as a whole exceeded 60% when the PDO cool phase ended and the cool phase began.Starting around 1980, the AAR began a steady decline that is projected to reach 25% by the end of the century (MM85).
In the Wrangell Mountains and the Chugach St. Elias Mountains, despite an increase in winter precipitation, April snow covered area (SCA) will decrease by more than 50%, from 25% to 12% in the Wrangell Mountains, and from 30% to 13% in the Chugach St Elias Mountains, with the decline commencing by the 2020s in both regions, indicating that these regions are crossing a critical temperature threshold that affects partitioning of precipitation into rain and snow, and the onset of snow melt (MM85).Mean April temperatures hover near the freezing point in the 2010s, explaining the onset of April SCA decline as warming advances.Presently, April SCA is near zero in the Alaska Range and the forested river valley and will remain so in the future.
The Tonsina and Tazlina Rivers are two of many glacier-fed anadromous tributaries to the Copper River flowing from the Chugach Mountains.Spring runoff is projected to increase in both rivers, attributable to increasing winter precipitation and snow melt.However, in the late summer, when flow is dominated by ice melt rather than snow melt, divergent trends in runoff from ice melt reveal that only Tonsina reaches a tipping point (a switch from increasing to decreasing ice runoff) during the simulation period due to critical loss of glacier area (Figure 10).In this example, both the Tonsina and the Tazlina catchment areas lost 4 km 2 of glacierized area after 3 decades of simulation under scenario MM85.From that point forward, glacier ice runoff decreased in the Tonsina but continued to increase in the Tazlina.The explanation is due to the relative loss of glacierized area in each catchment (10% and 1% of total glacierized area for Tonsina and Tazlina catchments, respectively).The ten-year average rates of glacier ice melt were similar for Tonsina and Tazlina catchments (−518 and −580 mm/year, respectively).To compensate for the km 2 of glacier area lost, the The AAR represents the glacier's capacity to grow by accumulating snow from year to year.The global balanced-budget AAR (where mass gain equals mass loss) is 58% [5] and AAR for CRW glaciers examined as a whole exceeded 60% when the PDO cool phase ended and the cool phase began.Starting around 1980, the AAR began a steady decline that is projected to reach 25% by the end of the century (MM85).
In the Wrangell Mountains and the Chugach St. Elias Mountains, despite an increase in winter precipitation, April snow covered area (SCA) will decrease by more than 50%, from 25% to 12% in the Wrangell Mountains, and from 30% to 13% in the Chugach St Elias Mountains, with the decline commencing by the 2020s in both regions, indicating that these regions are crossing a critical temperature threshold that affects partitioning of precipitation into rain and snow, and the onset of snow melt (MM85).Mean April temperatures hover near the freezing point in the 2010s, explaining the onset of April SCA decline as warming advances.Presently, April SCA is near zero in the Alaska Range and the forested river valley and will remain so in the future.
The Tonsina and Tazlina Rivers are two of many glacier-fed anadromous tributaries to the Copper River flowing from the Chugach Mountains.Spring runoff is projected to increase in both rivers, attributable to increasing winter precipitation and snow melt.However, in the late summer, when flow is dominated by ice melt rather than snow melt, divergent trends in runoff from ice melt reveal that only Tonsina reaches a tipping point (a switch from increasing to decreasing ice runoff) during the simulation period due to critical loss of glacier area (Figure 10).In this example, both the Tonsina and the Tazlina catchment areas lost 4 km 2 of glacierized area after 3 decades of simulation under scenario MM85.From that point forward, glacier ice runoff decreased in the Tonsina but continued to increase in the Tazlina.The explanation is due to the relative loss of glacierized area in each catchment (10% and 1% of total glacierized area for Tonsina and Tazlina catchments, respectively).The ten-year average rates of glacier ice melt were similar for Tonsina and Tazlina catchments (−518 and −580 mm/year, respectively).To compensate for the km 2 of glacier area lost, the Tonsina catchment required an additional −57 mm/year of ice melt but it produced only −30 mm/year, which explains the start of the downward trend.To compensate for km 2 of glacier area loss in the Tazlina catchment, an additional −5.8 mm/year of ice melt was required.The Tazlina melt rate actually increased by −22.8 mm/year, explaining the upward trend in runoff from ice melt.A temperature sensitivity analysis was performed for the time period 1960-2009 by uniformly increasing the historical temperature record by 2, 4, 6, 8 and 10 degrees C in order to identify the temperature rise at which a regime change is detectable in glacier ice melt trends.Simulations using the unaltered temperature record show that a statistically significant increasing trend in runoff from ice melt occurred from 1960 to 2000.With a 2 degree rise in temperature, the trend with time is essentially flat (slightly positive but not statistically significant).With temperature increases from 4 to 10 degrees, the trend with time becomes increasingly negative and statistically significant (Table A temperature sensitivity analysis was performed for the time period 1960-2009 by uniformly increasing the historical temperature record by 2, 4, 6, 8 and 10 degrees C in order to identify the temperature rise at which a regime change is detectable in glacier ice melt trends.Simulations using the unaltered temperature record show that a statistically significant increasing trend in runoff from ice melt occurred from 1960 to 2000.With a 2 degree rise in temperature, the trend with time is essentially flat (slightly positive but not statistically significant).With temperature increases from 4 to 10 degrees, the trend with time becomes increasingly negative and statistically significant (Table 6), suggesting that a regime change occurs with temperature rise approximately 2 degrees C, marking the transition from increasing to decreasing annual runoff from ice.Annual peak runoff in all regions of the CRW demonstrates significant trends over time, but the direction of the trend and the month of the peak vary by region.In the Alaska Range, peak runoff will decline, and the timing of peak runoff demonstrates significant change in mid-century.Peak runoff in the Alaska Range consistently occurs in June until mid-century, then it becomes quite variable, ranging from April to September.Peak runoff will increase in all other regions of the CRW.The timing of peak runoff in the Wrangells will consistently range from June to August with no apparent trend with time.In the Chugach and St Elias Mountains, peak runoff will continue to occur in July.In the Forested River Valley, where peak runoff is controlled by snow melt, the date of peak runoff will shift from May to April as snow melt occurs earlier in the spring.At the outlet of the watershed, peak runoff will continue to occur in July with occasional peaks occurring in August.
Total streamflow in the CRW is projected to increase sharply through mid-century for both RCP 4.5 and RCP 8.5 as temperatures continue to rise, however, the projections diverge in the second half of the century for the two scenarios; while streamflow continues to increase in both scenarios, it increases at a much faster pace in RCP 8.5 relative to RCP 4.5 (Figure 11).Annual peak runoff in all regions of the CRW demonstrates significant trends over time, but the direction of the trend and the month of the peak vary by region.In the Alaska Range, peak runoff will decline, and the timing of peak runoff demonstrates significant change in mid-century.Peak runoff in the Alaska Range consistently occurs in June until mid-century, then it becomes quite variable, ranging from April to September.Peak runoff will increase in all other regions of the CRW.The timing of peak runoff in the Wrangells will consistently range from June to August with no apparent trend with time.In the Chugach and St Elias Mountains, peak runoff will continue to occur in July.In the Forested River Valley, where peak runoff is controlled by snow melt, the date of peak runoff will shift from May to April as snow melt occurs earlier in the spring.At the outlet of the watershed, peak runoff will continue to occur in July with occasional peaks occurring in August.
Total streamflow in the CRW is projected to increase sharply through mid-century for both RCP 4.5 and RCP 8.5 as temperatures continue to rise, however, the projections diverge in the second half of the century for the two scenarios; while streamflow continues to increase in both scenarios, it increases at a much faster pace in RCP 8.5 relative to RCP 4.5 (Figure 11).Historically in July, the Forested River Valley has been predominantly hot and dry, but less so in other months.In the future, the forested river valley will be predominantly hot and dry during June, July, and August.Change point analysis indicates that a regime change for forest vulnerability may occur by 2027 for June, and by 2050 for August (Figure 12); this is supported by analysis of slope, analysis of means, and the Wilcoxon rank sum test.Prior to the 2020s, less than 20% of the forest is hot and dry on average in August, exceeding 50% only rarely.Between the 2020s and the 2050s, the Historically in July, the Forested River Valley has been predominantly hot and dry, but less so in other months.In the future, the forested river valley will be predominantly hot and dry during June, July, and August.Change point analysis indicates that a regime change for forest vulnerability may occur by 2027 for June, and by 2050 for August (Figure 12); this is supported by analysis of slope, analysis of means, and the Wilcoxon rank sum test.Prior to the 2020s, less than 20% of the forest is hot and dry on average in August, exceeding 50% only rarely.Between the 2020s and the 2050s, the vulnerable fraction jumps to 50% then appears to stabilize between 50% and 70% in the second half of the century.

Discussion
The modeling framework used in this study was previously applied in the CRW [29].In that study, calibration and evaluation of the model were performed using observations of runoff from USGS stream gauges, satellite-based observations of ET (MOD16) [24] snow extent (MOD10CM) [23], and GRACE-derived estimates of total storage change, and by comparing glacier mass balance change to values reported in the scientific literature by Larsen et al. [42], Arendt et al. and Das et al. [43].In the present study, the additional evaluation of the first 11 years of simulations to observed runoff and GRACE storage was also successful.Our simulations of AAR, which ranged from 36% to 69% during the period 1949 to 2000, resemble the mean and variability of values reported by Dyurgerov [5] for approximately 20 glaciers worldwide for which observations were available.At the projected trajectory of AAR decline (−0.37% per year), glaciers will lose the ability to replenish

Discussion
The modeling framework used in this study was previously applied in the CRW [29].In that study, calibration and evaluation of the model were performed using observations of runoff from USGS stream gauges, satellite-based observations of ET (MOD16) [24] snow extent (MOD10CM) [23], and GRACE-derived estimates of total storage change, and by comparing glacier mass balance change to values reported in the scientific literature by Larsen et al. [42], Arendt et al. and Das et al. [43].
In the present study, the additional evaluation of the first 11 years of simulations to observed runoff and GRACE storage was also successful.Our simulations of AAR, which ranged from 36% to 69% during the period 1949 to 2000, resemble the mean and variability of values reported by Dyurgerov [5] for approximately 20 glaciers worldwide for which observations were available.At the projected trajectory of AAR decline (−0.37% per year), glaciers will lose the ability to replenish mass before the end of the next century.Our projections of change in glacier ELA are also within the range of values reported by Dyurgerov [5].
If the projected two-month extension of time when temperatures are above freezing is realized, it will trigger a regime shift in the CRW with regard to snow, glaciers, runoff, and forest vulnerability.Between the current decade and the end of the century, April snow extent is projected to decrease by 60% (both RCP scenarios).Total runoff in March will experience a regime shift in mid Century when snowmelt and rainfall begin to supplement winter base flow.In the coastal mountains, runoff in all months will increase significantly, and the same is true for the interior mountains with one exception: winter base flow in January and February will decrease as glacierized area declines.A regime shift with regard to forest health is projected to occur mid-century, when most of the forested river valley will experience hot and dry conditions for three months, versus one month now.Young et al. (2017) found that wildfire occurrence in the CRW escalates when forests experience both air temperature of 13.4 • C or greater and a positive moisture deficit.Wildfire is not the only risk to forests under stress; heat stress and moisture deficit also make forests vulnerable to insect and microbial infestations.The dramatic variability in the timing of peak runoff in the Alaska Range, which could point to a regime shift, is unexplained but does not appear to be related to future variability in rainfall; with only one exception from 2006 to 2009, peak precipitation occurs in June, July and August where peak runoff is projected to be variable from April to September in the second half of the century.Glacier regime changes appears to be well underway, with AAR declining since 1980, ELA rising, glacier area decreasing, and mass balance increasingly negative.
An abrupt change in CRW glacier extent was not detected for the watershed as a whole.However, when small glaciers that feed tributaries lose mass, abrupt changes in tributary runoff may occur, as demonstrated by the comparison of the Tonsina and Tazlina Rivers.Hagg et al. [4] reported that temperature-related increases in the rate of ice melt compensate for reductions in glacier area; however, after a certain point, the rising rate of ice melt cannot compensate for loss of glacier area as demonstrated by the Tonsina River example, where an additional 57 mm of glacier thickness change would have been required to offset the 4 km 2 of glacier area loss, but the glacier melt rate increased by only 30 mm, triggering the start of a long term trend of declining ice melt runoff from that glacier.Glacierized area that is between 10% and 40% of the catchment area serves to stabilize runoff against precipitation fluctuations; during this century, the glacierized area of the CRW (MM85) will decrease from 20% to 15% so an increase in year to year variability of streamflow is not expected at the watershed outlet.However, in tributaries such as the Tonsina and Tazlina Rivers, glacierized area is only 4% and 5.3% of the catchment area (MM85), making them more vulnerable to precipitation variability.
Flow regime changes in the Copper River and its anadromous tributaries are consequential to the fisheries that are the foundation of economic, nutritional, and cultural life in Southcentral Alaska.Anadromous fish travel upstream to spawn at the place of their birth, starting the journey at predictable times of year, and it is uncertain how changes in the magnitude and timing of flow in the Copper River may impact their migration; during the period 1988-2009, most sockeye entered the Copper River in June when monthly runoff was on the order of 160 mm, but at the end of the century, June runoff may be 50% higher at 240 mm.Increased streamflow is known to significantly increase the time required for salmon to reach spawning areas in the Upper Copper River [44].With peak monthly flow projected to increase and total annual flow projected to double, identification of the critical point at which high flows prevent anadromous fish from reaching their destination, or make winter habitat unsuitable, is beyond the scope of this research, but crossing that tipping point would have significant consequences for the local economy and residents.Increased streamflow and potential channel migration is also an important consideration for flooding or erosion of roads, bridges and other infrastructure [16].
Considerable uncertainty is likely associated with our result.A primary source of uncertainty is the limited availability of meteorological stations in Alaska to form the basis for extrapolation to unmonitored areas [45].The downscaling of coarse-resolution CMIP5 global circulation models introduces bias, and some deficiencies are known to exist in CMIP5 simulations of northern latitudes, such as incomplete simulation of the effects of aerosols on clouds [1], overestimation of active layer thickness and the permafrost-carbon feedback, and generally poor representation of the pacific decadal oscillation [46].Further inaccuracies may arise from errors in measurement and interpolation of observed data used for input and calibration, processes not represented by the hydrologic model such as groundwater exchanges with surface water, and simplification of hydrologic and glacial processes by the selected conceptual modeling framework.It is unknown if parameters calibrated to historical conditions will remain appropriate for future climate conditions [47].Additional issues are equifinality, covariance of parameters, and model nonlinearity, which can be better understood using Fourier amplitude sensitivity test (FAST) or using numerous parameter ensembles with generalized likelihood uncertainty estimates (GLUE); these solutions were beyond the present project scope.

Conclusions
Using a calibrated conceptual glacio-hydrological water balance model, we evaluated future water portioning in a high-latitude glacierized watershed in Southcentral Alaska.The parsimonious model, the MWBMglacier, was previously calibrated and evaluated successfully against ground measurements of streamflow and satellite-based observations of snow covered area (MODIS), evapotranspiration (MODIS), and total water storage (GRACE).Additionally, simulated glacier changes in mass balance and area were compared favorably to data in published literature.Additional evaluations against observed data were performed for the first 11 years of future simulations.Multiple models of future climate from SNAP representing scenarios RCP 4.5 and RCP 8.5, and monthly 2 km resolution, were used to drive the model.
Regime shifts in glaciers, snow, runoff, and forest vulnerability were identified as having already occurred, or were projected to occur in the study area under plausible scenarios of future climate.Although runoff from glacier melt increases throughout the watershed under warming temperatures, we observed the point at which glacier ice runoff in one tributary begins to decline due to the critical loss of glacierized area, a fate that is likely for other glacierized catchments as glacier loss continues into the next century.The hot and dry conditions that increase forest vulnerability to disease and wildfire will occur during June, July and August in the future, as opposed to only in July presently.Streamflow will increase by almost 50% as precipitation and glacier runoff increase, so, in the coming decades, the Copper River may not be as vulnerable to the extended droughts and rising stream temperatures predicted for snow-dominated rivers in lower latitudes.However, given the importance of salmon fisheries to the cultural and economic well-being of Southcentral Alaska, the impacts that dramatically increased spring and summer discharge may have on anadromous fish migrating up the Copper River and its tributaries warrant further research.

Water 2017, 9 , 23 Figure 1 .
Figure 1.Site map showing the Copper River watershed, locations of glaciers, the river and its tributaries, four ecoregions, mountain ranges, and USGS stream gauges.The watershed outlet is at the Million Dollar Bridge.

Figure 1 .
Figure 1.Site map showing the Copper River watershed, locations of glaciers, the river and its tributaries, four ecoregions, mountain ranges, and USGS stream gauges.The watershed outlet is at the Million Dollar Bridge.

Water 2017, 9 ,
x FOR PEER REVIEW 5 of 23

Figure 2 .
Figure 2. The MWBMglacier schematic.Storage elements, indicated by boxes, account for water stored in soil, snow, glaciers, and in temporary storage for delayed runoff.Triangles indicate partitioning steps (t uses temperature to partition precipitation into snow and rain, SF partitions snow into sublimation and additions to snowpack, DF partitions the sum of rain and snow melt into direct runoff and additions to soil storage, and RF partitions available runoff into runoff during the current month, or runoff that is temporarily stored for delayed runoff).

Figure 2 .
Figure 2. The MWBMglacier schematic.Storage elements, indicated by boxes, account for water stored in soil, snow, glaciers, and in temporary storage for delayed runoff.Triangles indicate partitioning steps (t uses temperature to partition precipitation into snow and rain, S F partitions snow into sublimation and additions to snowpack, D F partitions the sum of rain and snow melt into direct runoff and additions to soil storage, and R F partitions available runoff into runoff during the current month, or runoff that is temporarily stored for delayed runoff).

Figure 3 .
Figure 3. Evaluation of model performance, using mean monthly data for the period 2006-2016, was performed by comparing simulated output from MM45 and MM85 with: (a) runoff at the watershed outlet; and (b) TWS change based on GRACE satellite data.USGS estimates of stream gauge error (15%, see Table1) are illustrated by gray shading in Figure3a.Error in GRACE estimates (illustrated by gray shading in Figure3b) accompanied the data provided by the University of Colorado GRACE Data Portal (listed in Table1).The reported error associated with individual monthly estimates of total water storage was doubled because the change in monthly storage was calculated for illustration in Figure3b.

Figure 3 .
Figure 3. Evaluation of model performance, using mean monthly data for the period 2006-2016, was performed by comparing simulated output from MM45 and MM85 with: (a) runoff at the watershed outlet; and (b) TWS change based on GRACE satellite data.USGS estimates of stream gauge error (15%, see Table1) are illustrated by gray shading in Figure3a.Error in GRACE estimates (illustrated by gray shading in Figure3b) accompanied the data provided by the University of Colorado GRACE Data Portal (listed in Table1).The reported error associated with individual monthly estimates of total water storage was doubled because the change in monthly storage was calculated for illustration in Figure3b.

3. 1
.4.Glaciers Storage of water in glaciers declines at statistically significant rates in all simulations.Glacier thickness loss (Figure 4) increases from 50 mm/year to 200 (MM4.5)and 350 (MM8.5) by 2100.To calculate the mass of water leaving glacier storage and entering the hydrologic cycle, annual glacier thickness change was multiplied by glacierized area; the cumulative volume of water leaving glaciers and entering the hydrologic cycle is 144 and 189 gigatons (Gt) from 2006 to 2099 for the RCP 4.5 and RCP 8.5 five-model means (MM45 and MM85), respectively.Based on the temperature sensitivity analysis, an additional 60 mm of glacier thickness is lost for each degree of temperature rise.Water 2017, 9, x FOR PEER REVIEW 11 of 23

3. 1
.4.Glaciers Storage of water in glaciers declines at statistically significant rates in all simulations.Glacier thickness loss (Figure 4) increases from 50 mm/year to 200 (MM4.5)and 350 (MM8.5) by 2100.To calculate the mass of water leaving glacier storage and entering the hydrologic cycle, annual glacier thickness change was multiplied by glacierized area; the cumulative volume of water leaving glaciers and entering the hydrologic cycle is 144 and 189 gigatons (Gt) from 2006 to 2099 for the RCP 4.5 and RCP 8.5 five-model means (MM45 and MM85), respectively.Based on the temperature sensitivity analysis, an additional 60 mm of glacier thickness is lost for each degree of temperature rise.

Figure 5 .
Figure 5. Projected changes in glacierized area based on 8 climate scenarios for: (a) RCP 4.5; and (b) RCP 8.5.The five dashed lines represent results from all CMIP5 models, and the solid black line is based on simulations using climate forcing data from the five-model mean.The red lines with circles represent HTP45 and HTP85 simulations.The red solid lines represent MM45THP and MM85THP simulations.The difference between the two red lines reflects the impact of temperature rise alone.The difference between the solid red line and the solid black line shows the mitigating effect of high winter precipitation on glacier loss.Gray shading illustrates the standard deviation of the five CMIP5 models.

Figure 5 .
Figure 5. Projected changes in glacierized area based on 8 climate scenarios for: (a) RCP 4.5; and (b) RCP 8.5.The five dashed lines represent results from all CMIP5 models, and the solid black line is based on simulations using climate forcing data from the five-model mean.The red lines with circles represent HTP45 and HTP85 simulations.The red solid lines represent MM45THP and MM85THP simulations.The difference between the two red lines reflects the impact of temperature rise alone.The difference between the solid red line and the solid black line shows the mitigating effect of high winter precipitation on glacier loss.Gray shading illustrates the standard deviation of the five CMIP5 models.

Water 2017, 9 , 23 Figure 5 .
Figure 5. Projected changes in glacierized area based on 8 climate scenarios for: (a) RCP 4.5; and (b) RCP 8.5.The five dashed lines represent results from all CMIP5 models, and the solid black line is based on simulations using climate forcing data from the five-model mean.The red lines with circles represent HTP45 and HTP85 simulations.The red solid lines represent MM45THP and MM85THP simulations.The difference between the two red lines reflects the impact of temperature rise alone.The difference between the solid red line and the solid black line shows the mitigating effect of high winter precipitation on glacier loss.Gray shading illustrates the standard deviation of the five CMIP5 models.

Figure 6 .
Figure 6.Trends in monthly runoff for scenarios RCP 4.5 (a) and RCP 8.5 (b) from the interior portion of the watershed, which comprises the Alaska Range, the Wrangell Mountains, and the forested river valley.Gray shaded regions display one standard deviation around the mean of five CMIP5 models, and where a solid black trendline is present, it indicates the presence of a statistically significant trend in the five-model mean simulations (MM45 or MM85).Runoff during high-flow months (identified by blue axes) is illustrated at a different scale than during low-flow months (black axes).

Figure 6 .
Figure 6.Trends in monthly runoff for scenarios RCP 4.5 (a) and RCP 8.5 (b) from the interior portion of the watershed, which comprises the Alaska Range, the Wrangell Mountains, and the forested river valley.Gray shaded regions display one standard deviation around the mean of five CMIP5 models, and where a solid black trendline is present, it indicates the presence of a statistically significant trend in the five-model mean simulations (MM45 or MM85).Runoff during high-flow months (identified by blue axes) is illustrated at a different scale than during low-flow months (black axes).

Figure 6 .
Figure 6.Trends in monthly runoff for scenarios RCP 4.5 (a) and RCP 8.5 (b) from the interior portion of the watershed, which comprises the Alaska Range, the Wrangell Mountains, and the forested river valley.Gray shaded regions display one standard deviation around the mean of five CMIP5 models, and where a solid black trendline is present, it indicates the presence of a statistically significant trend in the five-model mean simulations (MM45 or MM85).Runoff during high-flow months (identified by blue axes) is illustrated at a different scale than during low-flow months (black axes).

Figure 7 .
Figure 7. Trends in monthly runoff for scenarios RCP 4.5 (a) and RCP 8.5 (b) from the coastal portion of the watershed, which comprises the Chugach and St. Elias Mountains.Gray shaded regions represent one standard deviation around the mean of five CMIP5 models, and where a solid black trendline is present, it indicates the presence of a statistically significant trend in the five-model mean simulations (MM45 or MM85).Runoff during high-flow months (identified by blue axes) is illustrated at a different scale than during low-flow months (black axes).

Figure 7 .
Figure 7. Trends in monthly runoff for scenarios RCP 4.5 (a) and RCP 8.5 (b) from the coastal portion of the watershed, which comprises the Chugach and St. Elias Mountains.Gray shaded regions represent one standard deviation around the mean of five CMIP5 models, and where a solid black trendline is present, it indicates the presence of a statistically significant trend in the five-model mean simulations (MM45 or MM85).Runoff during high-flow months (identified by blue axes) is illustrated at a different scale than during low-flow months (black axes).

Water 2017, 9 ,
x FOR PEER REVIEW 15 of 23

Figure 8 .
Figure 8. Future water budget from the RCP 8.5 five-model mean simulation.The following components of the water budget increase at a statistically significant rate (p < 0.05): precipitation, runoff from rain and snow melt, runoff from glacier ice melt, and transfers to the atmosphere (ET and sublimation).

Figure 8 .
Figure 8. Future water budget from the RCP 8.5 five-model mean simulation.The following components of the water budget increase at a statistically significant rate (p < 0.05): precipitation, runoff from rain and snow melt, runoff from glacier ice melt, and transfers to the atmosphere (ET and sublimation).

1980s 23 Figure 9 .
Figure 9. Illustration of temperature regime changes detected in the RCP 8.5 five-model mean (MM85) using: change point visualization (a); and cumulative sum analysis (b).

Figure 9 .
Figure 9. Illustration of temperature regime changes detected in the RCP 8.5 five-model mean (MM85) using: change point visualization (a); and cumulative sum analysis (b).

Water 2017, 9 ,
x FOR PEER REVIEW 17 of 23

Figure 10 .
Figure 10.Trends in runoff from glacier ice melt (a,b), total runoff (c,d), and the ratio of ice melt to total runoff in two glacier-fed tributaries in the Chugach Mountains: the Tonsina River gage (left column) and the Tazlina River gage (right column).Results displayed are based on RCP 8.5 fivemodel mean (MM85), with the standard deviation of individual CMIP5 models shown in gray.In scenario RCP 4.5 (not shown) ice melt from both glaciers stops increasing in mid-Century.

Figure 10 .
Figure 10.Trends in runoff from glacier ice melt (a,b), total runoff (c,d), and the ratio of ice melt to total runoff in two glacier-fed tributaries in the Chugach Mountains: the Tonsina River gage (left column) and the Tazlina River gage (right column).Results displayed are based on RCP 8.5 five-model mean (MM85), with the standard deviation of individual CMIP5 models shown in gray.In scenario RCP 4.5 (not shown) ice melt from both glaciers stops increasing in mid-Century.

Figure 11 .
Figure 11.Future trends in total streamflow for RCP 4.5 and RCP 8.5 (MM45 and MM85).Total streamflow is the sum of runoff from all sources.The vertical green lines indicate break points in the slope of total runoff.

Figure 11 .
Figure 11.Future trends in total streamflow for RCP 4.5 and RCP 8.5 (MM45 and MM85).Total streamflow is the sum of runoff from all sources.The vertical green lines indicate break points in the slope of total runoff.

Water 2017, 9 , 23 Figure 12 .
Figure 12.The percentage of forests that are vulnerable (hot and dry) in June, July, and August using historical (1949-2009) and future (2010-2099) climate scenarios.Hot is defined as temperature greater than 13.4 C, and dry is defined as PET > ET.Blue lines show future results for the five-model mean simulations (MM45 and MM85) and gray shading illustrates the standard deviation of the five individual CMIP5 models.

Figure 12 .
Figure 12.The percentage of forests that are vulnerable (hot and dry) in June, July, and August using historical (1949-2009) and future (2010-2099) climate scenarios.Hot is defined as temperature greater than 13.4 C, and dry is defined as PET > ET.Blue lines show future results for the five-model mean simulations (MM45 and MM85) and gray shading illustrates the standard deviation of the five individual CMIP5 models. snap.uaf.edu/).

Table 1 .
Sources of data.

•
Storage in snow, glaciers, soil, and the Delay lagging reservoir, • Precipitation partitioning into snow and rain, • Runoff generated from rain, snowmelt, and glacier ice melt, and •

Table 4 .
Components of the water budget summarized by decade for RCP 4.5 and RCP 8.5.The values shown represent the five-model mean and the uncertainty represents the standard deviation of the five CMIP5 models.

Table 4 .
Components of the water budget summarized by decade for RCP 4.5 and RCP 8.5.The values shown represent the five-model mean and the uncertainty represents the standard deviation of the five CMIP5 models.

Table 5 .
Regime changes in the Copper River Watershed.

Table 5 .
Regime changes in the Copper River Watershed.

Table 6 .
Temperature sensitivity of annual runoff from melting glacier ice.

Table 6 .
Temperature sensitivity of annual runoff from melting glacier ice.