Estimating Stage-Frequency Curves for Engineering Design in Small Ungauged Arctic Watersheds

: The design of hydraulic structures in the Arctic is complicated by shallow relief, which cause unique runoff processes that promote snow-damming and refreeze of runoff. We discuss the challenges encountered in modeling snowmelt runoff into two coastal freshwater lagoons in Utqia˙gvik, Alaska. Stage-frequency curves with quantified uncertainty were required to design two new discharge gates that would allow snowmelt runoff flows through a proposed coastal revet-ment. To estimate runoff hydrographs arriving at the lagoons, we modeled snowpack accumulation and ablation using SnowModel which in turn was used to force a physically-based hydraulic runoff model (HEC-RAS). Our results demonstrate the successful development of stage-frequency curves by incorporating a Monte Carlo simulation approach that quantifies the variability in runoff timing and volume. Our process highlights the complexities of Arctic hydrology by incorporating significant delays in runoff onset due to localized snow accumulation and melting processes. This methodology not only addresses the uncertainty in snow-damming and refreeze processes which affect the arrival time of snowmelt inflow peaks, but is also adaptable for application in other challenging environments where secondary runoff processes are predominant.


Introduction
Snowmelt runoff can be delayed by small-scale processes like snow-damming, meltwater storage in snowpacks, and refreezes.Although these processes are common in many regions that experience snowmelt, they are amplified in low gradient, extremely cold Arctic regions [1][2][3].Common hydrologic modeling approaches used elsewhere typically ignore these processes with limited consequence.However, in the Arctic, snow-damming resulting from runoff collecting in snow-filled channels can delay runoff onset on the order of days and weeks [1].
Schramm et al. [1] applied TopoFlow to the Imnavait Creek watershed in Alaska and found that the onset of simulated discharge occurred 7 days earlier than measured discharge.They attributed this in part due to snow damming, which was not incorporated into the model due to a poor understanding of the process.Previously, Hinzman and Kane [4] also struggled with the complexities of snow dams when applying the Swedish Hydrologiska Byråns Vattenbalansavdelning (HBV) model to the same watershed.They applied an empirical calibration parameter (MAXBAS) between 1 and 5 days to account for the delay.Pohl et al. [2] used WATCLASS, a water and energy model [5], to simulate snowmelt runoff in Trail Valley Creek, located in a small Arctic watershed in Canada.They also noted the model predicted early runoff inception, which they attributed to meltwater storage and percolation through snow and snow damming in the streams.A review by Biu et al. [6] discusses the key hydrologic models used for surface Arctic hydrology, namely Topoflow, HBV, Soil and Water Assessment tool (SWAT), Ecological Model for Applied Geophysics (ECOMAG), and the Cold Regions Hydrological Model (CRHM), but interestingly, they do not mention limitations due to delay processes, perhaps because the impact of snow damming likely tends to decrease for larger watersheds [4].
The formation of snow dams is a very local process that is currently not implemented in common snowmelt runoff models and would be difficult to reasonably parameterize.Xia and Woo [7] attempted to numerically simulate snow dams, but their model was at the scale of a single snow dam, and not suitable for application at the watershed scale.This reality makes it extremely challenging to accurately calibrate and validate runoff models when modeled runoff peaks may be days or even weeks earlier than observations.The challenge of accurately modeling refreeze and snow damming delays at the watershed scale represents a significant capability gap that remains unsolved.
Typically engineering studies require stage-frequency curves to guide design decisions as a quantitative tool used to balance project cost and resilience.In the case of an outlet structure on a reservoir, a standard modeling approach would involve calibrating a runoff model to gauge data, then the routing of historic inflows to create a time-series of reservoir elevations.Often, the annual series, or the peak stage for each year, would then be converted statistically into a stage-frequency curve.The challenges related to accurately modeling runoff delays presents challenges for the calibration of the runoff model and the approach breaks down.
Recognizing these limitations, our study proposes a novel approach to creating the stage-frequency curves in an Arctic waterbody without directly modeling the complex delay processes.We bypass the traditional calibration of models to simulate lagoon stages, which has proven nearly impossible due to the unpredictable nature of snow dams and refreezing.Instead, we calibrate our model to total runoff volumes, which are not sensitive to the delay processes, allowing us to compute a range of possible lagoon stages as functions of the variability in the length and intensity of melt and runoff period.This method, which uses a Monte Carlo (MC) framework, integrates different runoff delay scenarios to derive stage-frequency curves with uncertainty for waterbodies receiving runoff through the snowmelt period.
We utilize SnowModel [8] to simulate the accumulation and ablation of the snowpack.SnowModel does not simulate runoff routing and must be paired with a routing model to link the simulation of snowpack processes to runoff fluxes across the landscape.Liston et al. use HydroFlow to route meltwater at the bottom of the snowpack simulated by SnowModel in one example [9].The objective of our methodology is to fill a critical gap in the hydrological engineering approaches for small Arctic watersheds where delay processes dominate.Our proposed method is intended to be agnostic to both the snowpack and the routing model employed.This flexibility broadens the utility of our approach, making it applicable to a wider range of Arctic hydrologic engineering scenarios.Therefore, the emphasis of this paper is focused on describing the novel method of developing stagefrequency curves and not the details of the snowpack evolution or runoff routing models.

Study Area
The US Army Corps of Engineers (USACE) Alaska District is planning to construct a coastal revetment in the city of Utqia ġvik, AK (Figure 1) to mitigate coastal erosion risk.Utqia ġvik is the northernmost incorporated community in the United States, situated on the coast of the Arctic Ocean well above the Arctic Circle on Alaska's North Slope.The rock revetment is designed to prevent erosion from wave action that is threatening the city.Besides protecting the city from the ocean, the design of the revetment must include the consideration of interior drainage capacity to accommodate seasonal snowmelt runoff.
A key concern is that the revetment could impede runoff discharge and cause interior flooding.To mitigate this risk, new discharge structures will be constructed to effectively pass snowmelt through to the ocean during the runoff period.There are two primary watersheds in Utqiaġvik (Figure 1), a 19.8 km 2 basin draining to the Middle Salt Lagoon and a 12.8 km 2 basin draining into a complex of small dams and lagoons above the Tasigarook Lagoon.Both lower lagoons are directly adjacent to the ocean and the alignment of the proposed revetement.In late May or early June of each year, the annual snowpack melts in an approximately two-week ablation and runoff period [10,11].The melt makes its way to the lagoons where it is retained until sand berms on the beach constructed to prevent seawater from entering the lagoons are breached by heavy equipment and the meltwater is released to the ocean.
The average annual precipitation in the region is approximately 170 mm and the annual runoff is roughly 110-120 mm [12,13].Snowmelt makes up the majority of the annual runoff, though the fraction of available snow water equivalent (SWE) that is converted to runoff ranges significantly from estimates as high as 90-100% [10,14] to measurements following periods of drought as low as 33% [15] and 29% [13].
Although there is a long history of scientific instrumentation in the region, there are no long-term runoff discharge measurements or lagoon stage records available.A formerly operated U.S. Geological Survey (USGS) stream gauge on Nunavak Creek (NC) (No. 15798700 operated from 1972 to 2004, Figure 1) directly adjacent to the study area provided the only long-term source of runoff data for model calibration.

Materials and Methods
A workflow was developed to produce stage-frequency curves for the lagoons located in basins which are subject to runoff delay processes from snow damming and meltwater refreeze.The workflow does not require parameterization to directly model the delay processes.The framework for this workflow consists of a snow accumulation and ablation model which calculates runoff or Liquid Water at Soil Surface (LWASS) on a 2D grid for past years and is forced with observed and reanalysis meteorological data.The There are two primary watersheds in Utqia ġvik (Figure 1), a 19.8 km 2 basin draining to the Middle Salt Lagoon and a 12.8 km 2 basin draining into a complex of small dams and lagoons above the Tasigarook Lagoon.Both lower lagoons are directly adjacent to the ocean and the alignment of the proposed revetement.In late May or early June of each year, the annual snowpack melts in an approximately two-week ablation and runoff period [10,11].The melt makes its way to the lagoons where it is retained until sand berms on the beach constructed to prevent seawater from entering the lagoons are breached by heavy equipment and the meltwater is released to the ocean.
The average annual precipitation in the region is approximately 170 mm and the annual runoff is roughly 110-120 mm [12,13].Snowmelt makes up the majority of the annual runoff, though the fraction of available snow water equivalent (SWE) that is converted to runoff ranges significantly from estimates as high as 90-100% [10,14] to measurements following periods of drought as low as 33% [15] and 29% [13].
Although there is a long history of scientific instrumentation in the region, there are no long-term runoff discharge measurements or lagoon stage records available.A formerly operated U.S. Geological Survey (USGS) stream gauge on Nunavak Creek (NC) (No. 15798700 operated from 1972 to 2004, Figure 1) directly adjacent to the study area provided the only long-term source of runoff data for model calibration.

Materials and Methods
A workflow was developed to produce stage-frequency curves for the lagoons located in basins which are subject to runoff delay processes from snow damming and meltwater refreeze.The workflow does not require parameterization to directly model the delay processes.The framework for this workflow consists of a snow accumulation and ablation model which calculates runoff or Liquid Water at Soil Surface (LWASS) on a 2D grid for past years and is forced with observed and reanalysis meteorological data.The runoff is routed through the basin terrain with a 2D hydraulic model and both runoff volumes and characteristic runoff hydrograph shapes are recorded at the lagoons.The runoff shapes represent the expected variability in runoff timing and intensity resulting from meteorological conditions in the forcing data, such as an abrupt warmup resulting in a melt out spanning a few days, or a cooler spring with a runoff period that lasts for several weeks.The runoff hydrograph shapes are scaled by volume and used to create a family of stagevolume curves for the lagoons, which vary depending on the timing of the runoff processes.Separately, the runoff volumes calculated from the forcing data are statistically fitted to produce runoff-volume frequency curves for each basin.The family of stage-volume curves (representing variability in the duration and intensity of the runoff timing) is combined with runoff volume-frequency curves (representing variability in total runoff volume for a given year) with a MC simulation to produce stage-frequency curves which represent the likely range of stages in the lagoon.The overall process for developing the stage-frequency estimates is shown in Figure 2.
Water 2024, 16, x FOR PEER REVIEW 4 of 20 runoff is routed through the basin terrain with a 2D hydraulic model and both runoff volumes and characteristic runoff hydrograph shapes are recorded at the lagoons.The runoff shapes represent the expected variability in runoff timing and intensity resulting from meteorological conditions in the forcing data, such as an abrupt warmup resulting in a melt out spanning a few days, or a cooler spring with a runoff period that lasts for several weeks.The runoff hydrograph shapes are scaled by volume and used to create a family of stage-volume curves for the lagoons, which vary depending on the timing of the runoff processes.Separately, the runoff volumes calculated from the forcing data are statistically fitted to produce runoff-volume frequency curves for each basin.The family of stage-volume curves (representing variability in the duration and intensity of the runoff timing) is combined with runoff volume-frequency curves (representing variability in total runoff volume for a given year) with a MC simulation to produce stage-frequency curves which represent the likely range of stages in the lagoon.The overall process for developing the stage-frequency estimates is shown in Figure 2.

Snow Accumulation and Ablation
SnowModel [8] was used to simulate snow accumulation and ablation in our study area (Figure 1).This model includes complex snow processes such as blowing-snow redistribution and sublimation, snowpack ripening and snowmelt precipitation, snow metamorphism, and snowpack runoff, to name a few.SnowModel has been used to accurately model snowpack in other complex snow environments that include Arctic watersheds [16].Metrological inputs to the model include air temperature, relative humidity, wind speed and direction, and precipitation data.SnowModel is designed to run on grid increments of 1 to 1000 m and temporal increments of 10 min to 1 day.
We performed snow modeling simulations from 1 September 1981 through 31 August 2021 using hourly time steps.Our modeling domain of about 179 km 2 (12,450 m by 14,400 m; 30 m grid cells) included the Middle Salt Lagoon, Tasigrook Lagoon, and NC watersheds (Figure 1).We used the 2015 North American Land Change Monitoring System (NALCMS) landcover, and a 0.46 m horizontal resolution digital elevation model derived from lidar in NAVD88 [17].
Meteorological observations from the National Weather Service station located at the Utqiaġvik Wiley Post-Will Rogers Airport station (STN 700260, Weather Bureau Army-

Snow Accumulation and Ablation
SnowModel [8] was used to simulate snow accumulation and ablation in our study area (Figure 1).This model includes complex snow processes such as blowing-snow redistribution and sublimation, snowpack ripening and snowmelt precipitation, snow metamorphism, and snowpack runoff, to name a few.SnowModel has been used to accurately model snowpack in other complex snow environments that include Arctic watersheds [16].Metrological inputs to the model include air temperature, relative humidity, wind speed and direction, and precipitation data.SnowModel is designed to run on grid increments of 1 to 1000 m and temporal increments of 10 min to 1 day.
We performed snow modeling simulations from 1 September 1981 through 31 August 2021 using hourly time steps.Our modeling domain of about 179 km 2 (12,450 m by 14,400 m; 30 m grid cells) included the Middle Salt Lagoon, Tasigrook Lagoon, and NC watersheds (Figure 1).We used the 2015 North American Land Change Monitoring System (NALCMS) landcover, and a 0.46 m horizontal resolution digital elevation model derived from lidar in NAVD88 [17].
Meteorological observations from the National Weather Service station located at the Utqia ġvik Wiley Post-Will Rogers Airport station (STN 700260, Weather Bureau Army-Navy [WBAN] 27502, Figure 1) and reanalysis data from the fifth-generation European Centre for Medium-Range Weather Forecast (ERA5) [18] were used to force the model.Although meteorological observations matching SnowModel's input precipitation variable are available within the study area, Liljedahl et al. [13] document significant underestimation of total winter precipitation measurements as observed at WBAN 27502.They calculate factors between 2.8 to 4.3 necessary to correct observed precipitation to match observed snow depths and water equivalents with significant variations year to year.To address the observational precipitation under-catch, we used hourly precipitation from the atmospheric model reanalysis product ERA5-Land [19], which is a downscaled version of ERA5 [18].ERA5 has been identified as the best reanalysis product for estimating precipitation in the Arctic [20].The hourly ERA5-Land data were downloaded with Google Earth Engine for the available period WY 1982-2021.
SnowModel uses an energy balance method to calculate surface energy exchanges and snowmelt.Air temperature, relative humidity, and wind speed and direction from the Utqia ġvik Airport are used to force the simulation of snow ablation.The density and water content of the snowpack are simulated using a single-layer snowpack evolution model with excess melt assumed to reach the base of the snowpack [8].Ultimately, SnowModel produces a gridded output of runoff, i.e., liquid water at the soil surface (LWASS) available to run downgradient.
Due to the extreme environment of our study site, there are limited physical measurements of SWE and snowmelt runoff.Consequently, to account for this deficiency of available data, we used several methods to validate the snow modeling effort.As a part of this validation, we used snow depth measurements from the Circumpolar Active Layer Monitoring (CALM) program and SWE measurements from the SnowNet field campaigns [21], and we compared these measurements to simulated SWE.We also compared our simulated snowmelt dates to observed snowmelt dates derived from Terra's Moderate Resolution Imaging Spectroradiometer (MODIS) [22] as well as observed snowmelt dates from time-lapse camera imagery from the Phenocam project [23,24].

Snowmelt Runoff
The USACE-HEC's River Analysis System (HEC-RAS) version 5.0.7 [25] was used to route the snowmelt runoff calculated by SnowModel.The model was selected to address a number of specific modeling constraints.First, the region's low gradient terrain includes the presence of many vegetated, drained thaw lake basins [13] and exhibits a "fill and spill" behavior of runoff progression [15].Shallow depressions must fill with runoff to a level that exceeds some local grade control before entering a successive depression, creating complex flow paths.In some cases, portions of watersheds can be "deranged", having no clearly defined drainage network [26].Second, the gridded runoff timeseries from SnowModel includes spatial variability resulting from wind redistribution of snow and topographic effects in the energy balance calculations.The runoff routing model must include both fine spatial topographic detail and the ability to incorporate spatially varied melt inputs to capture the complexity of the region's hydrology.
We used the HEC-RAS 2D Diffusive Wave Approximation to the Shallow Water Equations to simulate runoff in a mesh model domain.The relatively high-resolution (0.46 m horizontal) terrain data available allowed us to capture the complex topography including dry lakes and shallow depressions that will control flow patterns during runoff.HEC-RAS allows users to define a spatially variable precipitation boundary condition which is applied to the mesh.We used this feature to define the snowmelt or LWASS timeseries as the model's hydrologic boundary condition.NetCDF files of the "roff" variable produced by SnowModel were converted to the HEC Data Storage System (DSS) format using a Python script and assigned as a meteorological "Precipitation" boundary condition.Bottom roughness (Table 1) was derived from the 2016 National Landcover Database (NLCD) [27] land use classification and values in the range suggested by HEC [28].The basin includes several hydraulic structures including road culverts and small earthen dams.These structures along with bathymetric information were used to create the geometry needed for the HEC-RAS model.Bathymetric data for the lakes were limited and were estimated using the best available information.Data for Middle Salt Lagoon were derived from some incidental depth measurements made during a geotechnical drilling project [29].Depths and storage assumptions for the lagoons above the Tasigarook Lagoon were based on information published in a dam safety Periodic Inspection Report [30].
The snow ablation model captures the effects of water storage as melt moves vertically through the snowpack, but because it is decoupled from the runoff routing model, it does not directly handle snowmelt storage of meltwater after it reaches the soil surface.Other models attempt to capture the effects of runoff storage in the snowpack, but because of the dominant impacts of snow-damming, these models still miss the timing of runoff onset by several days [1].Instead of viewing our model as a precise simulation of the runoff process, an unattainable goal with currently available methods, we consider it a tool to capture the effects of spatially varied melt and terrain effects, such as storage in dry lakes, albeit with runoff unimpeded by intermediate snow storage, refreezes, and snow-damming.In the absence of these delaying processes, simulated runoff will tend to begin at concentrated locations lower in the basin and earlier than observations.Simulated peaks may be larger or smaller than observations.However, as the lack of attenuation due to snow dams during a rapid melt period may cause an earlier high peak, and storage behind snow dams can slowly accumulate during a milder melt period, a large peak from a sudden breach can still occur.
Simulated runoff from NC was used to calibrate both the snow accumulation and ablation model, and the runoff model to measured discharge at the NC gage.The watershed area above the gage is approximately 7.5 km 2 and its runoff is partially regulated by Lake Emaiksoun (area ~1.8 km 2 ).Flows recorded at the NC gage were made manually by USGS staff for the ascending limb, then a new datum was set at the stilling well recorder for the rest of the season (Personal communication, Matt Schellekens, 3 February 2022).Our calibration effort focused on runoff volumes rather than the usual peak timing and magnitude approaches.The calibrated models were then applied to the Middle Salt Lagoon and Tasigarook Lagoon watersheds.
The SnowModel snow accumulation and ablation model and HEC-RAS runoff models were both calibrated to the volumes of snowmelt discharge at the NC gage.We manually partitioned the gauge flow by extrapolating the descending limb of the runoff hydrograph.The cumulative precipitation volumes measured at the airport at the onset of runoff were significantly less than the cumulative NC snowmelt runoff volumes (Figure 3).This trend agrees with the findings of Liljedahl et al. [13] in a nearby basin.In contrast, cumulative precipitation from ERA5-Land tended to be higher than the snowmelt volumes at NC, which is expected because there will be losses due to terrain storage, infiltration, and evaporation.
A correction factor of 1.4 was applied to the ERA5-land precipitation primarily to calibrate the snow depths in SnowModel to field measurements.With this correction, all but two (1988 and 1999) of the twenty-three years ) of ERA5-Land precipitation volumes are higher than the runoff volume at NC. Errors from the NC measurements, particularly from the descending limb of the hydrograph for which discharge was not directly measured, or errors in our manual partitioning of snowmelt, may contribute to the apparent mismatch in volumes for those two years.
Water 2024, 16, x FOR PEER REVIEW 7 of 20 but two (1988 and 1999) of the twenty-three years ) of ERA5-Land precipitation volumes are higher than the runoff volume at NC. Errors from the NC measurements, particularly from the descending limb of the hydrograph for which discharge was not directly measured, or errors in our manual partitioning of snowmelt, may contribute to the apparent mismatch in volumes for those two years.

Stage Frequency and Uncertainty
The end products of this study were stage-frequency curves for two lagoons adjacent to the proposed coastal revetment.The Tasigarook and Middle Salt Lagoons were both included in the HEC-RAS model, and simulations of runoff included routing through them and their existing outlet structures.Tasigarook is a small lagoon with limited storage, and thus requires discharge through the outlet structure to prevent overtopping of roadways even for relatively small runoff volumes.In contrast, Middle Salt Lagoon has a storage capacity sufficient to store some runoff without breaching the sand berm on the beach and allowing discharge.Stage-frequency curves were developed for both lagoons assuming the existing outlet structures were fully open allowing free discharge, and for the additional case where the Middle Salt Lagoon outlet remained closed, and runoff was stored without release.
A key challenge of this study was including the impact of snow dam processes on the stage-frequency curves for the subject lagoons.Although we could extract peak stages from the lagoons from our simulated runoff results from each year, snow dams and refreezing were not directly modeled and the stages did not necessarily represent the peak stage for that specific year.To address this, we developed a unique method to obtain stage-frequency curves by separating our analysis into three parts.First, we developed volume-frequency curves which were insensitive to the delay processes.Second, we developed stage-volume curves for each watershed using a generalization of runoff

Stage Frequency and Uncertainty
The end products of this study were stage-frequency curves for two lagoons adjacent to the proposed coastal revetment.The Tasigarook and Middle Salt Lagoons were both included in the HEC-RAS model, and simulations of runoff included routing through them and their existing outlet structures.Tasigarook is a small lagoon with limited storage, and thus requires discharge through the outlet structure to prevent overtopping of roadways even for relatively small runoff volumes.In contrast, Middle Salt Lagoon has a storage capacity sufficient to store some runoff without breaching the sand berm on the beach and allowing discharge.Stage-frequency curves were developed for both lagoons assuming the existing outlet structures were fully open allowing free discharge, and for the additional case where the Middle Salt Lagoon outlet remained closed, and runoff was stored without release.
A key challenge of this study was including the impact of snow dam processes on the stage-frequency curves for the subject lagoons.Although we could extract peak stages from the lagoons from our simulated runoff results from each year, snow dams and refreezing were not directly modeled and the stages did not necessarily represent the peak stage for that specific year.To address this, we developed a unique method to obtain stage-frequency curves by separating our analysis into three parts.First, we developed volume-frequency curves which were insensitive to the delay processes.Second, we developed stage-volume curves for each watershed using a generalization of runoff responses as a function melt intensity observed in the simulated results.Last, we combined the curves to obtain peak stage-frequency curves for each lagoon using a MC simulation to capture uncertainty in both the volume-frequency curve and the stage-volume curves.We compiled the annual series of melt runoff volumes from the simulated results, based on a sliding 35-day window to capture the melt period.The selection of a 35-day window ensured that most of the runoff was captured for any given year, although the melt period was typically only 2-3 weeks in duration.The runoff volume annual series was fit with a Log-Pearson Type-III (LP3) distribution using the Bulletin 17C [31] methodologies in the USACE-HEC's Statistical Software Package (HEC-SSP, version 2.3) [32] to obtain a volume-frequency curve for each studied watershed.
Next, we investigated the effect of snow dams on runoff response to develop a generalized approach to accounting for the process of snow dam filling and breaching.Xia and Woo [7] describe in detail how snow dams fail due to erosion from seepage or overflow.The breach progression depends on the energy inputs and the rate of melt.For example, if melt is gradual, reservoirs may fill slowly and erode the dam thermally resulting in a gradual release of meltwater in the basin.On the other hand, if meteorological conditions cause a rapid melt, the dams may overtop and breach from the erosive forces of the flowing water over their crests.Overtopping dams tend to cause a rapid cascading failure progression as upstream breaches lead to failures of dams lower in the basin [7].From this we surmise that while snow damming and breaching processes may delay the onset of runoff, the basin response is still strongly tied to the rate of melt; gradual melt will result in a moderated release of runoff over a longer time span, while rapid melt will produce a shorter event with a distinct period of peak flow.
To incorporate snow dam failing modes, we identified distinctive patterns of runoff hydrograph "shapes" from the population of simulation results that represent a range of different runoff responses driven by the meteorological conditions.The patterns generally represented years with (1) sudden rapid melt with the majority of runoff occurring in just a few days, (2) a moderate period of melt ramping up to a single peak over the course of nearly a week, (3) a two-peaked response where early melt from solar radiation was followed by a second peak a few days later driven by sensible heat as the air temperatures warmed above freezing, and (4) a low and slow response with small diurnal peaks spanning over two weeks.Four different simulation years with distinct runoff hydrograph shapes were identified and normalized by volume.These normalized hydrographs encode the runoff response of each basin as a function of the intensity of the melt period without the need to consider the timing of onset which is not important for developing stage-volume curves.
To develop peak stage-volume curves for each lagoon, the normalized hydrographs were scaled to a range of volumes coinciding with the range of the runoff volume-frequency curves previously developed with Bulletin 17C.Bulletin 17C provides guidelines for estimating flood frequency and magnitude and recommends using the LP3 distribution to extrapolate flood values.Goodness of fit tests using the annual maximum dataset at Tasigarook and Middle Salt Lagoon, shown in Table 2, suggest LP3 is a good fit compared to other distributions (Generalized Pareto is generally used for partial duration data).Given its recommendation in Bulletin 17C and goodness of fit comparison, LP3 was selected to perform the frequency analysis.
The scaled hydrographs were hydraulically routed through the lagoons with the HEC-RAS model, and peak stage-volume curves were computed for each of the four runoff patterns.For the additional case, of Middle Salt Lagoon with the outlet closed, the peak stage-volume curve is insensitive to the runoff timing (or the runoff hydrograph "shape") and is therefore a single curve regardless of the runoff delay processes.
We combined the volume-frequency curves with the peak stage-volume curves, using the latter as a transfer function to produce a peak stage-frequency curve for each lagoon.The curves were combined using an MC framework similar to that used by the HEC Flood Damage Assessment (FDA) software (version 1.4.1).For each iteration of the MC simulation, a bootstrap method was used to alter the annual series of runoff volumes and create a variant of the volume-frequency curve fit by resampling the analytical distribution parameters determined from the HEC-SSP 17C analysis.The parameters for the LP3 distribution were randomly sampled to produce each realization.This curve was Water 2024, 16, 1321 9 of 18 combined with a random selection of one of the four equally weighted peak stage-volume curves representing each runoff pattern to create a unique stage-frequency curve for that iteration (Figure 4).The scaled hydrographs were hydraulically routed through the lagoons with the HEC-RAS model, and peak stage-volume curves were computed for each of the four runoff patterns.For the additional case, of Middle Salt Lagoon with the outlet closed, the peak stage-volume curve is insensitive to the runoff timing (or the runoff hydrograph "shape") and is therefore a single curve regardless of the runoff delay processes.
We combined the volume-frequency curves with the peak stage-volume curves, using the latter as a transfer function to produce a peak stage-frequency curve for each lagoon.The curves were combined using an MC framework similar to that used by the HEC Flood Damage Assessment (FDA) software (version 1.4.1).For each iteration of the MC simulation, a bootstrap method was used to alter the annual series of runoff volumes and create a variant of the volume-frequency curve fit by resampling the analytical distribution parameters determined from the HEC-SSP 17C analysis.The parameters for the LP3 distribution were randomly sampled to produce each realization.This curve was combined with a random selection of one of the four equally weighted peak stage-volume curves representing each runoff pattern to create a unique stage-frequency curve for that iteration (Figure 4).This process was repeated numerous times until the expected 1% annual exceedance probability (AEP) converged to an error tolerance of 2%.This process is described in EM This process was repeated numerous times until the expected 1% annual exceedance probability (AEP) converged to an error tolerance of 2%.This process is described in EM 1110-2-1619 [33] and implemented in the HEC Flood Damage Assessment (FDA) software [34].
Convergence testing is a way of determining whether a MC simulation has enough realizations to produce accurate results.The following approach was used to check for convergence: where is the edge of a confidence interval around the expected value for a given α, S is the standard deviation of the metric, X is the mean of the metric, and N is the sample size.α was set at 10% (setting confidence intervals at 90%) and the error tolerance was set at 2%.
The MC simulation is run until the convergence criteria is met and is therefore the termination of the simulation is independent of the number of iterations.

Snow Accumulation and Ablation
Snow model validation was carried out by comparing simulated SWE to observed SWE (converted from snow depth measurements) at the CALM site and SWE measurements from field campaigns [21].Additionally, we compared simulated SWE depletion dates with MODIS SCA from the CALM and camera sites.We also evaluated simulated snowmelt dates and compared these dates to the imagery at the time-lapse camera site.More details about this comparison are available in Wagner et al. [35].Figure 5 shows the simulated SWE at the CALM site for WY 1982-2021.The maximum simulated SWE varies from 0.077 to 0.18 m.In general, the snowmelt occurs in June.

Snowmelt Runoff
The runoff hydrographs computed at the NC gage location were generally not well matched with the observed hydrographs (Figure 6).The Kling-Gupta efficiency (KGE) and Pearson-R statistics for each year's hydrographs (Table 3) demonstrate the poor performance of the runoff model for the majority of years.There are some years, such as 2001 and 2004, where the statistic metrics, and the visual alignment of the hydrographs indicated the model performed well.Modifications to the roughness parameters in the hydraulic model changed the runoff peak timing on the order of hours while the difference between the observed and modeled peaks were typically on the order of days or weeks, limiting the effectiveness of traditional hydraulic model calibration methods.The maximum 35-day volumes of modeled runoff tended to be higher than the observed volumes (negative percent bias, or p-bias, in Table 3).This indicates that the volume of input melt minus losses to storage on the modeled terrain is larger than the observed runoff.The consistent bias toward model overestimation of volume is likely in part because of the conservative limitation of losses during the runoff process.During the CALM field campaigns only measurements of snow depths were performed.Snow depth (h s ) is related to SWE: where ρ b is the bulk density and the ρ w is the density of water.Because snow density was not measured during the CALM field campaigns, an estimation of snow density was required to derive SWE values.In 2008 and 2009, Liljedahl et al. [13] measured snow densities of 320 kg/m 3 at nearby sites.During the late season, CRREL field campaigns snow densities varied between 288 and 326 kg/m 3 .A commonly used snow density when converting to SWE is 300 kg/m 3 .Therefore, when converting the CALM snow depths to SWE using Equation ( 2), we used a snow density of 300 kg/m 3 .
During most years, SWE aligns well with the measured SWE at the CALM site.It is likely that the difference in simulated and measured SWE is because of the simplified conversion of snow depth to SWE using a constant snow density.The simulated SWE is within the range of the CRREL field campaign measurements [21] collected at a nearby site (Figure 4).

Snowmelt Runoff
The runoff hydrographs computed at the NC gage location were generally not well matched with the observed hydrographs (Figure 6).The Kling-Gupta efficiency (KGE) and Pearson-R statistics for each year's hydrographs (Table 3) demonstrate the poor performance of the runoff model for the majority of years.There are some years, such as 2001 and 2004, where the statistic metrics, and the visual alignment of the hydrographs indicated the model performed well.Modifications to the roughness parameters in the hydraulic model changed the runoff peak timing on the order of hours while the difference between the observed and modeled peaks were typically on the order of days or weeks, limiting the effectiveness of traditional hydraulic model calibration methods.The maximum 35-day volumes of modeled runoff tended to be higher than the observed volumes (negative percent bias, or p-bias, in Table 3).This indicates that the volume of input melt minus losses to storage on the modeled terrain is larger than the observed runoff.The consistent bias toward model overestimation of volume is likely in part because of the conservative limitation of losses during the runoff process.We reviewed measurements of terrain albedo made at the National Oceanic and Atmospheric Administration (NOAA) Global Monitoring Laboratory (GML) station (https://gml.noaa.gov/aftp/data/radiation/baseline/brw/(accessed on 1 March 2022)), just outside the study area and found that for the period 1998-2004 runoff at NC typically did not begin until the albedo had decayed to around 0.3 and the air temperature was above  We reviewed measurements of terrain albedo made at the National Oceanic and Atmospheric Administration (NOAA) Global Monitoring Laboratory (GML) station (https: //gml.noaa.gov/aftp/data/radiation/baseline/brw/(accessed on 1 March 2022)), just outside the study area and found that for the period 1998-2004 runoff at NC typically did not begin until the albedo had decayed to around 0.3 and the air temperature was above freezing.This suggests that the arrival of longer periods of sunlight increases the incoming shortwave solar radiation triggering melt, but that meltwater is held locally on the landscape until the exposure of low albedo terrain and availability of sensible heat from the air provide the energy input necessary to allowed it to mobilize.

Stage Frequency and Uncertainty
The snow accumulation, ablation and runoff models were applied to the basins draining to the Tasigarook and Middle Salt Lagoons.Hydrographs and the volume of runoff during a 35-day window starting at runoff onset at the lagoons were compiled for 40 years (WY 1982(WY -2021)).Basic volume statistics are summarized in Table 4. Simulated runoff hydrographs from the years 1983, 1985, 1989 and 2014 were selected to represent the range of runoff shapes that the basin could experience based on meteorological conditions.The hydrographs were normalized to the 35-day volume for each year then scaled to a range of runoff volumes spanning the volume-frequency curve (Figure 7).The individual years were selected based on the visual variation in each runoff pattern.We investigated the variation in meteorological conditions that drove the change in each pattern and found them to be primarily related to the air temperature (Figure 8).We see, in 1983, a rapid warming from below 0 °C to over 10 °C in the course of just 3 days, with minimum temperatures remaining above freezing for several days.This forcing produces the concentrated spike seen in the hydrograph shape.In 1985, air temperatures cycled through above freezing temperatures just above 5 °C, while falling back to around −5 °C at night.This forcing produced two distinct spikes in the hydrograph, about 10 days apart, when the air temperatures rose above 5 °C.In 1989, the air temperatures reached around 5 °C falling to around 0 °C or just below at night.This caused the hydrograph to have a distinct peak, like in 1983, but spread over a longer time as the air temperatures were not as high.In 2014, the air temperatures straddled the freezing point, not reaching 5 °C until around 3 weeks after the start of melt.This represents a very slow thaw and a muted runoff hydrograph.The shape of the simulated runoff hydrographs was most sensitive to the range of air temperatures relative to the freezing point, and in this particular basin, where the temperatures fall within the 5 °C to −5 °C band seems to be important.The individual years were selected based on the visual variation in each runoff pattern.We investigated the variation in meteorological conditions that drove the change in each pattern and found them to be primarily related to the air temperature (Figure 8).We see, in 1983, a rapid warming from below 0 • C to over 10 • C in the course of just 3 days, with minimum temperatures remaining above freezing for several days.This forcing produces the concentrated spike seen in the hydrograph shape.In 1985, air temperatures cycled through above freezing temperatures just above 5 • C, while falling back to around −5 • C at night.This forcing produced two distinct spikes in the hydrograph, about 10 days apart, when the air temperatures rose above 5 • C. In 1989, the air temperatures reached around 5 • C falling to around 0 • C or just below at night.This caused the hydrograph to have a distinct peak, like in 1983, but spread over a longer time as the air temperatures were not as high.In 2014, the air temperatures straddled the freezing point, not reaching 5 • C until around 3 weeks after the start of melt.This represents a very slow thaw and a muted runoff hydrograph.The shape of the simulated runoff hydrographs was most sensitive to the range of air temperatures relative to the freezing point, and in this particular basin, where the temperatures fall within the 5 • C to −5 • C band seems to be important.
The scaled hydrographs were then routed through the HEC-RAS model of the lagoons, which are controlled by existing outlet structures.Peak stages from each simulation were extracted from the lagoon stage hydrographs and used to construct stage-volume curves for each hydrograph shape (Figure 9).The single-peaked shape from 1983, representing a rapid melt, produced the highest stages as inflows greatly exceeded the outlet capacity, while the gradual melt-out of the 2014 shape produced the lowest stages.The scaled hydrographs were then routed through the HEC-RAS model of the lagoons, which are controlled by existing outlet structures.Peak stages from each simulation were extracted from the lagoon stage hydrographs and used to construct stage-volume curves for each hydrograph shape (Figure 9).The single-peaked shape from 1983, representing a rapid melt, produced the highest stages as inflows greatly exceeded the outlet capacity, while the gradual melt-out of the 2014 shape produced the lowest stages.The scaled hydrographs were then routed through the HEC-RAS model of goons, which are controlled by existing outlet structures.Peak stages from each s tion were extracted from the lagoon stage hydrographs and used to construct stag ume curves for each hydrograph shape (Figure 9).The single-peaked shape from representing a rapid melt, produced the highest stages as inflows greatly exceed outlet capacity, while the gradual melt-out of the 2014 shape produced the lowest s  The bootstrapped LP3 fit was obtained for each realization of the MC simulation by randomly perturbing the input annual series of runoff volumes for each watershed.This produced a unique volume-frequency curve each time, capturing the uncertainty in the volume annual series due to the length of the series.Most of the variation was in the high and low probability events.
The stage-frequency curves computed from the 35-day volume curve are shown in Figure 10.The metric to test convergence was the 1% AEP.The 1% AEP convergence is near 14,000 realizations for Middle Salt Lagoon and 31,500 realizations for Tasigarook Lagoon for the stages generated from the 35-day volume curve.The convergence is not depended on the number of realizations.The distribution of stages created using the stage-volume curves as a transfer function are shown in Figure 10.Based on our results, a pool stage in Middle Salt Lagoon resulting overtopping of an adjacent roadway (Cakeater Road) with an elevation of 2.3 m has approximately a 5% annual chance of exceedance.
volume annual series due to the length of the series.Most of the variation was in the high and low probability events.
The stage-frequency curves computed from the 35-day volume curve are shown in Figure 10.The metric to test convergence was the 1% AEP.The 1% AEP convergence is near 14,000 realizations for Middle Salt Lagoon and 31,500 realizations for Tasigarook Lagoon for the stages generated from the 35-day volume curve.The convergence is not depended on the number of realizations.The distribution of stages created using the stagevolume curves as a transfer function are shown in Figure 10.Based on our results, a pool stage in Middle Salt Lagoon resulting overtopping of an adjacent roadway (Cakeater Road) with an elevation of 2.3 m has approximately a 5% annual chance of exceedance.

Discussion
Snow dams, meltwater refreezes, and shallow terrain add significant complexity to modeling hydrologic processes in the Arctic.Although we have identified likely causes for runoff delay, a complete physically-based modeling framework suitable for engineering applications that captures these phenomena at a basin scale has not been developed to date.Simply ignoring the process by applying commonly used runoff modeling methods yields results that do not capture the uncertainty the delay processes introduce and is an unattractive approach and would likely produce inaccurate results.
The method presented which focuses on calibration of total volumes and treatment of uncertainty with a MC analysis demonstrates an approach that avoids the thorny issue of accounting for snow damming and runoff refreezes.This represents a significantly different approach than what other authors have pursued.Schramm et al. [1] and Pohl et al. [2] simply note the delay causes issues with timing, though their objective is mainly documenting model performance, not using the results for further application.Hinzman and Kane [4] do attempt to parameterize the delay, but with a somewhat basic empirical calibration factor.This approach is not particularly useful in studies where calibration data are not available for all years and is therefore not viable when trying to reconstruct historical stages for statistical analyses.
The limited efforts to directly model snow damming physical processes, such as those by Xia and Woo [7] are worth continuing to pursue.It would be preferable to have a more physically based model that could be calibrated to observations, and where the delay processes could be parameterized by physical conditions in the basin, and meteorological forcing.Improving these complex models may have received limited attention [6] because in larger watersheds snow damming and localized runoff delay processes may not be dominant [4].

Discussion
Snow dams, meltwater refreezes, and shallow terrain add significant complexity to modeling hydrologic processes in the Arctic.Although we have identified likely causes for runoff delay, a complete physically-based modeling framework suitable for engineering applications that captures these phenomena at a basin scale has not been developed to date.Simply ignoring the process by applying commonly used runoff modeling methods yields results that do not capture the uncertainty the delay processes introduce and is an unattractive approach and would likely produce inaccurate results.
The method presented which focuses on calibration of total volumes and treatment of uncertainty with a MC analysis demonstrates an approach that avoids the thorny issue of accounting for snow damming and runoff refreezes.This represents a significantly different approach than what other authors have pursued.Schramm et al. [1] and Pohl et al. [2] simply note the delay causes issues with timing, though their objective is mainly documenting model performance, not using the results for further application.Hinzman and Kane [4] do attempt to parameterize the delay, but with a somewhat basic empirical calibration factor.This approach is not particularly useful in studies where calibration data are not available for all years and is therefore not viable when trying to reconstruct historical stages for statistical analyses.
The limited efforts to directly model snow damming physical processes, such as those by Xia and Woo [7] are worth continuing to pursue.It would be preferable to have a more physically based model that could be calibrated to observations, and where the delay processes could be parameterized by physical conditions in the basin, and meteorological forcing.Improving these complex models may have received limited attention [6] because in larger watersheds snow damming and localized runoff delay processes may not be dominant [4].
Hydraulic modeling in small Arctic watersheds is challenging in part due to very sparse datasets available for model calibration and validation [6].Without any available lagoon stage data in the study area, we were unable to assemble an annual series of stages to compare these results to, and if those were available, a direct calculation of the stage frequency estimates could be made.
The method we describe is driven by the annual series of snowmelt runoff volumes that are not significantly impacted by runoff delay processes.In its current form, the approach does require the practitioner to select and apply probabilistic weighting to a collection of characteristic hydrographs assumed to be representative of the variability in runoff possibilities.We derived these from simulated runoff hydrographs and attempt to link the runoff patterns to meteorological forcing, namely air temperature patterns.Although this approach does have limitations due to the bias introduced, there does not appear to be an alternative method available in the literature to calculate stage-frequency curves in basins with dominant delay processes.Future improvements could include a more robust linkage between modeled snow conditions, meteorological forcing, and the selection of the characteristic runoff hydrograph.However, even in its current form, this framework could be adapted to other similar engineering applications where runoff delay processes are known to be dominant, but quantitative estimates of flow or stage are necessary for design.

Conclusions
This study sought to develop a methodology to calculate stage-frequency curves, in an environment where standard methods failed due to the complexities of snow damming and runoff refreeze.We applied a novel approach using a MC simulation which captured variability in the runoff behavior and uncertainty in the probability of total runoff volumes.Our results demonstrate that even without explicitly modeling the runoff delay processes, we can effectively capture the overall shape of the stage-frequency curve.This is driven by runoff volume, along with the variability of the stage, which is driven by the runoff intensity, a function of meteorological forcing.This is a significant improvement over attempting to fit a statistical distribution to an annual series of modeled historical stages which would likely be incorrect due to the inability to model the runoff delay processes.Our presented method provides a practical framework for engineering practitioners in the Arctic who have a need to develop stage-frequency curves which are critical for the design of engineered hydraulic structures.The inclusion of confidence intervals allows these curves to be used in designs that incorporate risk-based decision-making.The method can also be easily transferred to other small Arctic basins that have dominant runoff-delay processes.
Despite the benefits of this method, there are some limitations.This method does not provide insights on the physical details of the delay processes or attempt to identify the factors that cause them to be more dominant in certain years.In addition, the assignment of probability of occurrence of each runoff hydrograph shape relies on the judgement of the modeler, which likely introduces some bias to the variability in the stage-frequency curves.Future research efforts should seek to identify the parameters which could be used to refine the assignment of probability of each shape.Synthetic meteorological forcings could also be used to develop a larger suite of runoff hydrograph shapes, instead of relying on representative shapes selected from historical modeling.Increased field collection of lake stage data in the Arctic is desirable, though time require to assemble the multi-year records necessary to statistically validate these methods would preclude these efforts in the near term.We demonstrate that the annual series of runoff volumes follows an LP3 distribution and assume that the resampled data in the MC realization follow the parent distribution.We recognize there is some distribution uncertainty in the resampled data that may warrant more detailed investigation outside the scope of this study.Future studies may consider layered MC methods incorporating distribution uncertainty.
In conclusion, our research fills a critical need for the practicing design engineer who must calculate stage-frequency curves in basins where calibration and/or validation of standard hydrologic models are difficult due to snow damming and runoff refreeze, which to date, cannot be accurately modeled.Funding: This study was funded by the US Army Corps of Engineers Alaska District.

Figure 1 .
Figure 1.Utqiaġvik (latitude: 71.290556, longitude: −156.788611)location map.Middle Salt Lagoon and Tasigarook Lagoon watersheds are labeled with the names of the lagoons at their downstream ends, adjacent to the ocean.

Figure 1 .
Figure 1.Utqia ġvik (latitude: 71.290556, longitude: −156.788611)location map.Middle Salt Lagoon and Tasigarook Lagoon watersheds are labeled with the names of the lagoons at their downstream ends, adjacent to the ocean.

Figure 2 .
Figure 2. Schematic of method to develop stage-frequency curves for lagoons.Yellow box is input data; red boxes represent processes; blue boxes represent outputs/inputs to the processes.

Figure 2 .
Figure 2. Schematic of method to develop stage-frequency curves for lagoons.Yellow box is input data; red boxes represent processes; blue boxes represent outputs/inputs to the processes.

Figure 3 .
Figure 3.Comparison of cumulative precipitation from the NWS weather station at the Utqiaġvik Wiley Post-Will Rogers Airport (STN 700260) and ERA5-Land hourly product with snowmelt runoff at the NC gage for the period 1982-2004.Cumulative precipitation is from the start of the water year ending when runoff begins at the NC gage.

Figure 3 .
Figure 3.Comparison of cumulative precipitation from the NWS weather station at the Utqia ġvik Wiley Post-Will Rogers Airport (STN 700260) and ERA5-Land hourly product with snowmelt runoff at the NC gage for the period 1982-2004.Cumulative precipitation is from the start of the water year ending when runoff begins at the NC gage.

Figure 4 .
Figure 4. Example iteration of MC simulation combination of bootstrapped fit volume-frequency curve combined with a randomly selected stage-volume curve transfer function to obtain a single updated stage-volume curve.

Figure 4 .
Figure 4. Example iteration of MC simulation combination of bootstrapped fit volume-frequency curve combined with a randomly selected stage-volume curve transfer function to obtain a single updated stage-volume curve.

Water 2024 , 20 Figure 5 .
Figure 5. SnowModel SWE (blue line) at the CALM location.CALM site measurements are in orange diamonds while the CRREL field campaign data are represented by red triangles where minimum and maximum SWE are illustrated with red error bars.

Figure 5 .
Figure 5. SnowModel SWE (blue line) at the CALM location.CALM site measurements are in orange diamonds while the CRREL field campaign data are represented by red triangles where minimum and maximum SWE are illustrated with red error bars.

20 Figure 6 .
Figure 6.Runoff hydrographs (in cubic meter per second, cms) modeled results (solid blue line) and observed (dashed black line) at NC for the period WY 1982-2004.

Figure 6 .
Figure 6.Runoff hydrographs (in cubic meter per second, cms) modeled results (solid blue line) and observed (dashed black line) at NC for the period WY 1982-2004.

Figure 7 .
Figure 7. Scaled snowmelt runoff discharge curve shapes for the (a) Middle Salt Lagoon and (b) Tasigarook Lagoon watersheds.Each colored line represents a different total volume of runoff.

Figure 7 .
Figure 7. Scaled snowmelt runoff discharge curve shapes for the (a) Middle Salt Lagoon and (b) Tasigarook Lagoon watersheds.Each colored line represents a different total volume of runoff.

Water 2024 , 20 Figure 8 .
Figure 8. Minimum and maximum air temperatures (blue and red line, respectively), and precipitation (light blue bars) driving the simulated runoff hydrographs for the WY 1983, 1985, 1989, and 2014.

Figure 9 .
Figure 9. Stage-volume curves for (a) Middle Salt Lagoon and (b) Tasigarook Lagoon for different inflow hydrograph shapes (WY 1983, 1985, 1989, 2014, and outlet closed for Middle Salt Lagoon).Each color (black, red, green, blue) represents a different hydrograph shape.The bootstrapped LP3 fit was obtained for each realization of the MC simulation by randomly perturbing the input annual series of runoff volumes for each watershed.This produced a unique volume-frequency curve each time, capturing the uncertainty in the

Figure 8 .
Figure 8. Minimum and maximum air temperatures (blue and red line, respectively), and precipitation (light blue bars) driving the simulated runoff hydrographs for the WY 1983, 1985, 1989, and 2014.

Figure 8 .
Figure 8. Minimum and maximum air temperatures (blue and red line, respectively), and p tation (light blue bars) driving the simulated runoff hydrographs for the WY 1983, 1985, 19 2014.

Figure 9 .
Figure 9. Stage-volume curves for (a) Middle Salt Lagoon and (b) Tasigarook Lagoon for d inflow hydrograph shapes (WY 1983, 1985, 1989, 2014, and outlet closed for Middle Salt La Each color (black, red, green, blue) represents a different hydrograph shape.The bootstrapped LP3 fit was obtained for each realization of the MC simulat randomly perturbing the input annual series of runoff volumes for each watershed produced a unique volume-frequency curve each time, capturing the uncertainty

Figure 10 .
Figure 10.Example stage-frequency curve for the (a) Middle Salt, and (b) Tasigarook Lagoons including uncertainty.

Figure 10 .
Figure 10.Example stage-frequency curve for the (a) Middle Salt, and (b) Tasigarook Lagoons including uncertainty.
funding acquisition, A.W. and E.D. All authors have read and agreed to the published version of the manuscript.

Table 1 .
Manning's Roughness values used for HEC-RAS model.

Table 2 .
Goodness of fit tests comparing distributions against annual maximum flows at Tasigarook and Middle Salt Lagoon.Lower values indicate a better fit to the dataset.

Table 3 .
Statistics comparing the modeled and observed hydrographs in Nunavak Creek, shown graphically in Figure6.

Table 4 .
Simulated basin runoff volumes in million cubic meters (MCM).