Estimation of Freshwater Discharge from the Gulf of Alaska Drainage Basins

Abstract

The freshwater discharge from catchments along the Gulf of Alaska, termed Alaska discharge, is characterized by significant quantity and variability. Owing to subarctic climate and mountainous topography, the Alaska discharge variations may deliver possible impacts beyond the local hydrology. While short-term and local discharge estimation has been frequently realized, a longer time span and a discussion on cascading impacts remain unexplored in this area. In this study, the Alaska discharge during 1982–2022 is estimated using the Soil and Water Assessment Tool (SWAT). The adequate balance between the model complexity and the functional efficiency of SWAT suits the objective well, and discharge simulation is successfully conducted after customization in melting calculations and careful calibrations. During 1982−2022, the Alaska discharge is estimated to be 14,396 ± 819 m³⋅s⁻¹⋅yr⁻¹, with meltwater contributing approximately 53%. Regarding variation in the Alaska discharge, the interannual change is found to be negatively correlated with sea surface salinity anomalies in the Alaska Stream, while the decadal change positively correlates with the North Pacific Gyre Oscillation, with reasonable time lags in both cases. These new findings provide insights into the relationship between local hydrology and regional climate in this area. More importantly, we provide rare evidence that variation in freshwater discharge may affect properties beyond the local hydrology.

1. Introduction

River discharge from continents to oceans accounts for a global average of 36,055 km³⋅yr⁻¹ [1]. It plays a non-negligible role in the Earth’s hydrological and biogeochemical cycles by delivering freshwater, essential sediments, and nutrients to the ocean [2,3,4]. The quantity of and variation in freshwater discharge are the most crucial parameters among many others [5,6]. Specifically, freshwater discharge can freshen the seawater and deliver impacts on large-scale ocean processes such as thermohaline circulation [7]. For example, a thermohaline circulation existing in the North Pacific Ocean has been shown to significantly influence the local and global climate and marine ecosystems via the redistribution of heat and the transportation of organic carbon and nutrients [8,9,10,11]. The strength of this circulation process is primarily controlled by the sea surface salinity. According to Uehara et al. [12], anomalies in the sea surface salinity may have impacts on the overturning in the northwestern Sea of Okhotsk. These anomalies can propagate along the current pathway from the Alaskan Stream, central Bering Sea, western subarctic gyre, East Kamchatka Current region, and eastern Okhotsk Sea (Figure S1). Following on from this finding, a significant negative correlation was revealed between the interannual variation in sea surface salinity in the northwestern Sea of Okhotsk and the variation in annual river discharge from the Kamchatka Peninsula, which is located upstream of the overturning along the current pathways [13]. Similarly, the role of freshwater discharge from the catchments along the Gulf of Alaska (GOA) (hereafter referred to as the Alaska discharge) was hinted at but has not yet been adequately explored. Therefore, the present study aims to investigate the correlation between the Alaska discharge and sea surface salinity anomalies in the Alaskan Stream, which is an upstream component of the salt pathway. The significant Alaska discharge into the neighboring oceans is a result of rainfall, snow and ice melt, and little human activity in the Alaskan region. Furthermore, variations in the Alaska discharge exhibit sensitivity to atmospheric phenomena over the North Pacific Ocean and to global climate change [14,15,16]. However, the long-term variations in the Alaska discharge and their connection with corresponding atmospheric variability remain little known. One potential factor influencing the hydrological processes in the catchments along the GOA is the North Pacific Oscillation (NPO), a large-scale atmospheric variability observed over the North Pacific Ocean. The NPO is characterized by a seesaw north–south dipole of sea level pressure anomalies between the Aleutian Low and North Pacific High [17,18]. Therefore, a comprehensive estimation of the Alaska discharge and its variation can contribute to a better understanding of the connections between regional discharge, large-scale atmospheric variations, and ocean processes in the North Pacific sector.

Various methods can be used to estimate discharge, including the rating curve method, the velocity-area method, hydrological modeling, and remote sensing techniques [19,20,21,22]. Among these options, hydrological models, which include empirical, conceptual, and process-based models, offer an adaptable approach to simulate the processes of water movement under various land cover and soil conditions [23,24,25]. Several studies on the Alaska discharge have primarily focused on estimating the total volume, analyzing the contributions of hydrological components, and examining seasonal variations by using various hydrological models. For instance, Wang et al. [15] simulated the Alaska discharge and its seasonal and interannual variations by developing a digital elevation model (DEM) that simplifies the main physical processes within the catchment, including rainfall-runoff, snow and ice storage, and melting processes. Neal et al. [26] estimated the annual Alaska discharge and its contribution from glaciers and icefields using a combination of measured discharge and enhanced precipitation models. Similarly, Hill et al. [27] developed regression equations to estimate the monthly Alaska discharge at high resolution based on various physical characteristics of the basin, including the area, mean elevation, and land cover, as well as meteorological characteristics such as temperature and precipitation. Moreover, Beamer et al. [28] estimated the Alaska discharge and examined its source components with good accuracy by using a suite of distributed, process-based models, with consideration of important hydrological processes such as rainfall-runoff, evapotranspiration (ET), infiltration, baseflow runoff, and snow and ice melting.

While these studies confirm that process-based hydrological models provide a comprehensive and efficient representation of real hydrological processes, there are still unresolved questions, particularly regarding the spatial variability of the model parameters. These parameters, which account for different catchment characteristics and climate forcing, are often underrepresented at the basin scale. Additionally, although these studies provide plausible historical discharge estimates (700–900 km³ during 1961–2014), an extensive exploration of recent estimations is necessary to better understand the long-term variations in the Alaska discharge, given the climatic variation. Furthermore, the influence of the Alaska discharge variations on coastal oceanic systems, such as sea surface salinity in the surrounding ocean sector, remains underexplored.

Hence, the present study builds on previous work by focusing on the following: (i) estimating the Alaska discharge during 1982–2022, using a process-based hydrological model, namely, the Soil and Water Assessment Tool (SWAT, [29,30]) with modified melting processes, (ii) investigating spatiotemporal variations in the Alaksa discharge, (iii) exploring the correlations between variations in the Alaska discharge and atmospheric phenomena, as well as sea surface salinity, to provide a more comprehensive understanding of the atmosphere–land–ocean interactions. First, the details of the study area, the SWAT application and its inputs, and the SWAT calibration and validation, as well as the statistical approach for analyzing the discharge, are introduced in Section 2. The model calibration, validation, and simulation results are then presented in Section 3, along with the connection between the Alaska discharge and atmospheric variability. The relationship between the Alaska discharge and the North Pacific Ocean is discussed in Section 4. Finally, the key conclusions are summarized in Section 5.

2. Materials and Methods

2.1. Study Area

The catchments along the GOA (hereafter referred to as the Alaska catchment) span the southern region of Alaska and a portion of Canadian territory (54–64° N, 125–157° W), bordering the GOA in the Northeastern Pacific Ocean (Figure 1a). The present study focuses on the Alaska catchment which covers an area of approximately 320,000 km², including five major basins, namely the following: the Susitna (44,391 km²), Copper (60,451 km²), Alsek (28,298 km²), Taku (11,013 km²), and Stikine (47,734 km²) basins, in addition to numerous ungauged basins (127,902 km²). The Alaska catchment has a large range of elevation, with an average height of around 1000 m and some peaks over 3000 m (Figure 1b). The region is covered by a number of mountains that contribute to the complex topography. The major mountain ranges that define the basins include the Coast Mountains in the southeast, the Wrangell-St. Elias Mountains in the east, the Chugach Mountains along the central coast, the Alaska Range stretching from east to west in the northern area, and the Kenai Mountains on the Kenai Peninsula. Particularly, the Wrangell-St. Elias Mountains, above 5000 m, are covered by extensive glaciers. Meanwhile, the Chugach Mountains along the central coast block moist Pacific air, leading to high precipitation on the coastal side and drier conditions inland. The Alaska Range, a long (~1000 km) west–east mountain range, includes the highest peak in North America and is characterized by a tough climate and significant glaciation. Due to the rugged topography, the following 5 slope classes were established across all subbasins in the Alaska catchment: 0–3%, 3–15%, 15–30%, 30–60%, and >60%, based on DEM [31] using ArcGIS 10.7.1 software. Most areas (78%) fall within the 3–60% slope range. Additionally, according to North America Land Cover 2015 [32], the study area consists of various land cover types (Table S1), dominated by forest (32.3%), particularly subpolar taiga, polar or subpolar shrubland–lichen–moss (25.5%), barren land (34.1%), glaciers (17.8%), and wetlands (3.48%).

The climate in the study area is diverse, due to the high latitude, complex topography, and proximity to the ocean [33]. The higher-latitude region is dominated by cold climates that result in long, cold winters and brief, cool summers, with the Alaska Range experiencing a polar climate [34,35]. Meanwhile, the lower-latitude region is characterized by boreal and warm temperate climates, thus leading to humid conditions associated with relatively mild winters and cool summers [33,34,36]. Lader et al. [37] used reanalysis data from 1979 to 2009 to analyze the near-surface air temperature and precipitation for the Alaska region. Their findings indicated that the average monthly temperature during winter is approximately −12 °C, with an average precipitation of around 6 cm, while the temperatures in the Alaska Range drop below −15 °C and the precipitation exceeds 15 cm. Except for parts of the Alaska Range, summer temperatures in the higher-latitude regions are approximately 10 °C, with precipitation ranging between 10–20 cm. Moreover, the southeastern region is the warmest part of the study area during winter, with average temperatures close to freezing, and a monthly average precipitation of above 20 cm. Although no uniform trend in the annual precipitation is observed across the Alaska catchment, a significant increase in the annual temperature is observed between 1957 and 2021, with a rate of roughly 0.4 °C⋅decade⁻¹ [38]. These variable climate conditions potentially influence the hydrological processes across the Alaska catchment.

These important characteristics such as topography and climate directly influence the Alaska discharge quantity and variability. The discharge is predominantly generated in the mountainous areas that receive substantial snow and rain, with the highest discharge rates observed in low-elevation regions covered by large coastal glaciers [28]. Additionally, the phase of precipitation input depends on the air temperature variability, leading to a strong seasonal cycle. The diverse climate of the Alaska catchment, varying from polar in the Alaska Range to temperate in the lower-latitude coastal regions, primarily controls the timing and magnitude of meltwater.

2.2. Model Setup and Configuration

Herein, the SWAT model is used to estimate the discharge during the study period of 1982–2022, as well as the warm-up period during 1979–1981 (methodological flowchart in Figure 2). The SWAT model is a semi-distributed process-based hydrological model that is capable of simulating the water transport in river basins. Specifically, the ArcSWAT 2012 extension (hereafter referred to as SWAT) of the ArcGIS 10.7.1-based graphical interface is used in this study to process the discharge simulation. The SWAT model delineates a basin and divides it into a series of subbasins depending on topography. These subbasins are further divided into the smallest unit, namely the hydrological response unit (HRU), based on the homogeneous type of land cover, soil, and slope. The SWAT model then calculates the water balance at the HRU level, aggregates it at the subbasin level, and finally routes it towards the channels and outlet of the basin [39]. The water balance is calculated using Equation (1) [30]:

$$S{W}_t=S{W}_0+ {{\displaystyle \sum }}_{i=1}^t\left(R_{day}-Q_{surf}-E_a-W_{seep}-Q_{gw}\right)$$

(1)

where SW_t is the final soil water content (mm), SW₀ is the initial soil water content (mm) on day 0, t is time (days), R_day is the amount of precipitation (mm) on day i, Q_surf is the amount of surface runoff (mm) on day i, E_a is the amount of evapotranspiration (mm) in day i, W_seep is the amount of water (mm) entering the vadose zone from the soil profile on day i, and Q_gw is the amount of groundwater discharge (mm) on day i.

The main hydrological processes performed in SWAT can be divided into two major phases. The first of these is the land phase of the hydrological cycle, in which the amount of water flowing into the main channel is estimated for each subbasin. The runoff, evapotranspiration, melting processes, and soil–water processes are simulated in this phase. The second division is the channel phase of the hydrological cycle which can be defined as the movement of water through the channel network to the basin outlet [30]. Various methods are used to simulate these processes numerically and to estimate the key components of the catchment water balance. The soil conservation service curve number is an essential method, and is widely used for empirically estimating surface rainfall-runoff [40]. This estimation is primarily based on the land cover and soil type of the catchment. The SWAT model is capable of estimating the runoff from frozen soil if the temperature in the first soil layer is less than 0 °C [30]. The constant land cover and soil types were used in the present study due to their minimal variation over time (see Section 2.3 for more details). Meanwhile, the Penman–Monteith method [41,42] is satisfactory for estimating potential evapotranspiration by considering various climate parameters such as temperature, humidity, wind speed, and solar radiation. The kinematic storage method [43] is used to estimate lateral flow in soil layers, accounting for the variations in conductivity, slope, and soil water content of the catchment. Additionally, the variable storage routing method [44] is used to estimate the movement of water through the channel network by considering channel characteristics such as width, length, and slope within the catchment. Therefore, the selection of SWAT as the foundational model is underpinned by its notable attributes: its complicated preprocessing is made relatively accessible through ArcGIS 10.7.1 software add-ins and its computation is relatively efficient owing to the semi-distributed structure with fewer repetitive routines. Moreover, the existing literature has verified that SWAT performs well in reproducing the hydrological processes in mountainous and glacierized basins [39,45].

The snowfall and melting processes are important components of the hydrological cycle in the present study, since the Alaska catchment is partially covered by snow and ice. The SWAT model distinguishes between rainfall and snowfall based on a comparison of the near-surface air temperature and the threshold snowfall temperature (SFTMP). Snowfall occurs when the mean daily air temperature is lower than the SFTMP, and snow can accumulate until melting occurs. Melting begins when the air temperature is higher than the threshold snowmelt temperature (SMTMP). The SWAT model allows for spatial variations in precipitation and temperature by accounting for changes in elevation at the subbasin level. A maximum of ten elevation bands can be defined in SWAT. In each elevation band, the precipitation and temperature are modified based on two lapse rates, namely the temperature lapse (TLAPS) and precipitation lapse (PLAPS). Given that the elevation of glacierized areas in the Alaska catchment is mostly above 2000 m, ten elevation bands were established in the present study.

In the default setting of SWAT, melting calculations are based on a simplistic temperature-index method that does not account for ice melting. This is potentially problematic, given the considerable contribution of meltwater from glacierized areas. Therefore, the melting process in glacierized areas is additionally modified by separate calculations based on the equilibrium-line altitude (ELA). The ELA is a theoretical line where the annual mass balance is zero and represents the boundary between the accumulation and ablation zones of a glacier. In this study, the ELA was initially estimated by using the accumulation area ratio (AAR), which is the ratio between a glacier’s accumulation zone and its total area [46]. More specifically, based on the available AAR observational data [47] for glaciers in the Alaska catchment, the average AAR value was estimated to be 0.62. The glaciers in each basin were treated as a group, rather than as individual glaciers. Next, the elevation data for the glaciers were used, along with their outlines from the Randolph Glacier Inventory (RGI) 6.0 dataset [48], to conduct the ELA calculations of the glaciers within the five basins. As the glacial volume loss makes a relatively low contribution (6–10%) to the discharge [26,27,28], the glacierized areas are treated as constant during the estimation. The ELAs were therefore estimated as 1378, 1602, 1287, 1331, and 1542 m for glaciers in the Susitna, Copper, Alsek, Taku, and Stikine basins, respectively. The average ELA for the five basins was 1428 m, which is similar to the available observed mean ELA for Alaska glaciers (1363 ± 103 m) during 1982–2020 [49]. Finally, the estimated ELA in each basin was applied to divide the elevation bands into accumulation and ablation zones within the glacierized subbasins. The lowest elevation band of the ablation zone corresponds to the location where 90% of the glaciers’ coverage is reached.

Separate melt factors were established for the ablation and accumulation zones to determine the spatially distinct melting processes, based on the methods described by Omani et al. [45]. The parameters used to estimate the snowfall, snow accumulation, and melting processes at the subbasin level were given fixed values and were well constrained by previous research and physical interpretations [45,50,51]. Thus, for the non-glacier zone, the SFTMP and SMTMP were each set as 0 °C, the maximum snowmelt factor (SMFMX) on 21 June was set as 5 mm⋅H₂O⋅°C⁻¹⋅day⁻¹, the minimum snowmelt factor (SMFMN) on December 21 was set as 3 mm⋅H₂O⋅°C⁻¹⋅day⁻¹, and the snowpack temperature lag factor (TIMP) was set as 1. For the ablation zone, the SFTMP, SMTMP, SMFMX, SMFMN, and TIMP were 1 °C, 0.5 °C, 6 mm⋅H₂O⋅°C⁻¹⋅day⁻¹, 3 mm⋅H₂O⋅°C⁻¹⋅day⁻¹, and 0.05, respectively. For the accumulation zone, the corresponding values were 2 °C, 2 °C, 7 mm⋅H₂O⋅°C⁻¹⋅day⁻¹, 4 mm⋅H₂O⋅°C⁻¹⋅day⁻¹, and 0.01, respectively. Seasonal snow melt and summer ice melt mainly occur in the ablation zone, while permanent snow and ice undergo accumulation and melting in the accumulation zone. The initial glacier depth was set at 15,300 mm by averaging the thickness of the observed glaciers (RGI 6.0) within the Alaska catchment. Overall, while this approach does not fully capture the complexities of melting processes, such as the glacial dynamics and behavior of individual glaciers, it effectively accounts for the major accumulation and melting processes driven by temperature, considering spatial variability and using reasonable melting-related parameters based on the existing literature [39,45,50]. Despite its simplicity, the SWAT model in combination with the temperature-index approach is suitable for melting calculations.

2.3. Model Inputs

The model input data are listed in

Table 1. The input data used in the SWAT model.

Data	Description	Spatiotemporal Resolution
DEM	ASTER Global Digital Elevation Model version 3	30 m grid
LC	North America Land Cover 2015	30 m grid
Soil	FAO Digital Soil Map of the World	1:5 000 000
Climate	Climate Forecast System Reanalysis (CFSR, 1979–2013); Climate Forecast System version 2 operational data (CFSv2, 2011–2022)	0.31° × 0.31° grid/daily; 0.5° × 0.5° grid/daily

. The 30 m resolution DEM data, obtained from ASTER Global Digital Elevation Model version 3, were used in SWAT to delineate basins and subbasins, and to compute the flow directions of the river system. Land cover data with a resolution of 30 m were obtained from the North America Land Cover 2015. The Alaska catchment area was divided into 9 land cover types. The Food and Agriculture Organization (FAO) Digital Soil Map of the World [52], with a scale of 1:5,000,000, was used to classify the 18 soil types that could be directly used in the SWAT soil database.

The land cover and soil inputs for SWAT were regarded as constant throughout the simulation, which could lead to uncertainties. However, no evidence that significantly contradicts the consideration of constant land cover and soil inputs was found in previous studies. More specifically, regarding the land cover distribution, most forested areas and other land cover types such as pasture, rangeland, built-up land, and cropland, did not undergo significant changes in northwestern North America, including in the Alaska catchment, during 1960–2020 [53,54,55]. Regarding the change in soil distribution, the process is more commonly long-term, while extreme and significant change can be brought about by short-term processes primarily associated with human activities. However, such short-term changes were found to be less evident in the Alaska catchment [56].

The meteorological data including daily precipitation (mm), maximum and minimum temperature (°C), solar radiation (MJ⋅m⁻²), relative humidity (%), and wind speed (m⋅s⁻¹), were used to force the hydrological processes in the SWAT simulation. Notably, the meteorological data were obtained from two separate datasets, namely the climate forecast system reanalysis (CFSR, 1979–2013, [57]) and the climate forecast system version 2 operational data (CFSv2, 2011–2022, [58]). Nevertheless, the CFSv2 can be considered as an extension of the CFSR, which has been successfully applied to simulate the hydrological processes in North America [28,57,58]. Therefore, the dual datasets were used to force the discharge simulation in this study.

2.4. Model Calibration and Validation

In this study, the SWAT Calibration and Uncertainty Program (SWAT-CUP) was used to calibrate the model parameters [59]. SWAT-CUP is an interface that applies the Sequential Uncertainty Fitting (SUFI-2) algorithm to determine the parameter values and their uncertainties via sequential and fitting processes [59]. Although many parameters influence the discharge estimation in SWAT, the availability of adequate parameter values based on observation is limited. To avoid excessive tuning, fewer parameters should be considered, with more attention being paid to the features that they represent [60]. Therefore, the present study mainly focuses on the most influential parameters based on previous studies that used SWAT to estimate the discharge in snow and glacierized basins [39,45,61,62]. These parameters account for the primary hydrological processes within the study area, including the surface runoff, evapotranspiration, and groundwater flow processes, thereby ensuring the efficiency and accuracy of the model calibration. The calibrated parameters included the curve number (CN2), which is a key parameter for estimating surface runoff, as it directly influences the quantity of precipitation that becomes runoff. The available water capacity of the soil layer (SOL_AWC, mm⋅H₂O⋅mm soil⁻¹) and the soil evaporation compensation (ESCO) are crucial to soil moisture and evaporation processes. The minimum snow water content that corresponds to 100% snow cover (SNOCOVMX, mm⋅H₂O) influences the areal depletion of snowpack coverage, which is essential for representing snowmelt timing and runoff volume. The threshold depth of water in the shallow aquifer required for return flow to occur (GWQMN, mm⋅H₂O), the groundwater delay time (GW_DELAY, days), and the baseflow alpha factor (ALPHA_BF) play an important role in routing meltwater through the subsurface, affecting the hydrological response of the watershed. The TLAPS (°C⋅km⁻¹,) and the PLAPS (mm⋅H₂O⋅km⁻¹) control the spatiotemporal distribution of precipitation, determining whether it falls as rain or snow. The model performance was evaluated by statistical metrics such as the Nash–Sutcliffe efficiency (NSE), the coefficient of determination (r²), and the percentage bias (PBIAS). The model result was considered acceptable when the critical threshold for each index was satisfied, i.e., NSE > 0.50, r² > 0.70, and |PBIAS| < 25% [63]. Positive values of PBIAS indicate an underestimation bias in the model, while negative values indicate an overestimation bias. More detailed information on SWAT can be found in the official document [64].

In the present study, the model calibration and validation were performed using limited monthly observed streamflow data from the outlets of the five basins. The observed data were obtained from the Global Runoff Data Centre (https://portal.grdc.bafg.de/applications/public.html?publicuser=PublicUser#dataDownload/Home). The streamflow gauging positions are shown in Figure 1a, and the specific details are presented in Table S2. For the Susitna basin, the streamflow observations for 1982–1988 were used for calibration, and those for 1989–1992 were used for validation. For the Copper basin, the calibration period was 1982–1985 and the validation period was 1986–1990.9. For the Alsek basin, the calibration period was 1991.7–1997 and the validation period was 1998–2012.11. For the Taku basin, the calibration period was 1987.8–1997 and the validation period was 1998–2008. For the Stikine basin, the calibration period was 1982–1997 and the validation period was 1998–2005. The model configurations including the setup of elevation bands and the calibrated parameters in each of the five basins were subsequently projected to their adjacent ungauged basins for further discharge estimation, according to the principle of spatial proximity and hydrological similarity [65,66]. Specifically, the setup of the elevation bands and the calibrated parameters in the Susitna basin were projected to its neighboring basins, namely, ungauged basins No. 1, 2, and 3 (Figure 1a). Similarly, the setup from the Copper basin was applied to ungauged basin No. 4, those from the Alsek basin to ungauged basins No. 5 and 6, and those from the Stikine basin to ungauged basins No. 7 and 8.

The results of model calibration and validation are presented in Section 3.1. The estimated Alaska discharge is then discussed as a whole in Section 3.2, and for each of the five basins and the ungauged basins in Section 3.3. Here, the discharge characteristics for each of the five individual basins are discussed to provide a better understanding of the regional heterogeneity. Additionally, the integrated discharge from the ungauged basins is discussed and compared with that from the five basins.

2.5. Statistical Analysis

Herein, the trend in the Alaska discharge was assessed by using the non-parametric Mann–Kendall Test and Sen’s slope estimator. The Mann–Kendall Test is widely used to detect significant trends in time series datasets [67,68]. A trend exists when the statistical values satisfy p < 0.05 and |Zs| > 1.96. A positive Zs value indicates an increasing trend, while a negative Zs indicates a decreasing trend. The Sen’s slope estimator determines the slope of the trend, thereby indicating the rate of change [69]. Furthermore, the relationship between the long-term variations (5-year running means) in the Alaska discharge and the North Pacific Gyre Oscillation (NPGO) index [70] (i.e., the time series of the NPGO) is determined. The NPGO is suggested to be driven by sea surface wind variability associated with the NPO, which is an atmospheric variability characterized by a seesaw north–south dipole between the Aleutian Low and the North Pacific High. Therefore, the NPGO index is referred to as the oceanic expression of the NPO [70], and is used to reveal the connection between variations in the Alaska discharge and the NPO. Moreover, the annual Alaska discharge is compared with the sea surface salinity anomalies for 1984–2008 in the Alaskan Stream, as reported by Uehara et al. [12], to investigate the effects of the Alaska discharge on the oceanic characteristics in the North Pacific sector.

3. Results

3.1. Model Calibration and Validation

The default values for the selected parameters, the methods of adjustment used, and the final calibrated values over each basin are presented in

Table 2. The SWAT parameters and their values used for calibration.

Parameter a	Range	Default Value	Calibrated Value
			Susitna	Copper	Alsek	Taku	Stikine
r_CN2	35–98	36–92 b	–0.19	–0.22	–0.29	–0.24	–0.18
r_SOL_AWC	0–1	0–0.175 b	–0.56	–0.53	–0.46	–0.46	–0.07
v_ESCO	0.01–1	0.95	0.06	0.48	0.77	0.69	0.37
v_SNOCOVMX	0–500	1	125.64	424.85	354.08	437.4	299.7
v_GWQMN	0–5000	1000	255.90	1450.01	1310.40	969.01	1257
v_GW_DELAY	0–500	31	26.60	34.87	49.23	34.94	25.86
v_ALPHA_BF	0–1	0.048	0.03	0.048	0.097	0.08	0.04
v_TLAPS	–10–10	–6.5	–4.73	–7.34	–5.82	–6.02	–7.64
v_PLAPS	–1000–1000	200	500	105.09	500	215.52	392.5

Notes: ^a Here: r represents a relative test where the parameter value is multiplied by 1+ the given value and v represents a replacement test where the parameter value is replaced by the given value. ^b Values determined by soil and land cover types in each HRU.

. The parameters CN2 and SOL_AWC were adjusted by relative tests wherein the parameter value was multiplied by 1 + the given value in the Table, as indicated by r in the parameter designation. The other parameters were modified by replacement tests wherein the parameter value was replaced by the given value, as indicated by v in the parameter designation. The calibrated parameters are within reasonable ranges. Among the calibrated parameters, ALPHA_BF is the most sensitive across the five basins (Figure S2). The model’s performance varies significantly when ALPHA_BF, baseflow alpha factor, is adjusted within the range of −0.1–0.03, whereas the model tends to stabilize and deliver satisfactory performance when its value exceeds 0.03 (Figure S3). The calibration and validation performance indices for each basin are given in

Table 3. The statistical results of the model performance.

Basin	Calibration			Validation
	NSE	r2	PBIAS (%)	NSE	r2	PBIAS (%)
Susitna	0.855	0.905	17.13	0.846	0.901	16.80
Copper	0.812	0.816	7.20	0.805	0.808	5.60
Alsek	0.787	0.841	21.40	0.845	0.877	16.10
Taku	0.827	0.869	13.46	0.872	0.879	5.45
Stikine	0.825	0.866	−10.76	0.848	0.927	−14.69

, which indicates that NSE values are greater than 0.5, r² values are greater than 0.7, and PBIAS values are within ± 25%. These results statistically demonstrate an acceptable model performance. Indeed, the simulated discharge hydrograph in Figure 3a,b agrees well with observations during both the calibration and validation periods, especially in the Alsek, Taku, and Stikine basins, which have relatively longer observational datasets than the Susitna and Copper basins. However, the results indicate that the SWAT model underestimates the discharge from the Susitna, Copper, Alsek, and Taku basins, while it overestimates the Stikine basin, especially in summer (Figure 3 and Table 3). This can be primarily attributed to uncertainties in the meltwater estimation related to the catchment characteristics, particularly the topography. Compared to the other basins, the elevation bands of the Stikine basin are predominantly distributed at higher altitudes. Although the use of elevation bands in SWAT is acceptable for estimating snow-related processes and improving hydrological performance, it tends to slightly overestimate the snow accumulation and subsequent snowmelt in high-elevation regions and underestimate them in low-elevation regions, due to uncertainties arising from the linear variation in the temperature and precipitation lapse rates [39]. Moreover, the Susitna and Copper basins largely encompass areas with slopes of less than 15% (56% and 52% of the basin area, respectively). This relatively gentle slope may introduce uncertainties in the hydrological processes by reducing the surface runoff and increasing the soil water storage and evaporation during the simulated summer season. Despite several discrepancies in the hydrograph, the calibrated model is satisfactory for representing the discharge processes, and is therefore considered to be capable of estimating the discharge from the five basins.

3.2. Alaska Discharge

The present subsection examines the combined Alaska discharge across the five basins plus the ungauged basins. The results in Figure 4 indicate that the annual Alaska discharge is approximately 14,396 ± 819 m³⋅s⁻¹ during the period of 1982–2022, with a substantial contribution from meltwater (53%; Figure 5). In particular, the finding that melt runoff makes a significantly larger contribution to the total discharge than other contributors, such as direct rainfall, is consistent with the results of Neal et al. [26] and Beamer et al. [28]. However, the estimated annual discharge is slightly less than the 19,000–26,920 m³⋅s⁻¹ reported in existing studies [15,26,27,28]. This is because the present study uses a smaller area to avoid uncertainties in the discharge estimation from flat areas near coastlines. SWAT is considered to be difficult to accurately delineate streams and catchments in such areas due to the presence of large but shallow depressions and subtle variations [71]. Additionally, the present study accounts for complicated hydrological processes in both the land phase and the channel phase, which may lead to a relatively high evapotranspiration and transmission loss.

The monthly discharge exhibits a strong seasonal cycle in response to temperature variations (Figure 4). Specifically, the monthly discharge experiences a relatively constant average baseflow of 4252 ± 906 m³⋅s⁻¹ from January to March compared to higher and more variable values in other seasons. This is primarily due to the low temperatures during January to March, which result in low meltwater (Figure 5). The Alaska discharge starts to increase rapidly in April as the temperature rises above 0 °C. This increase is closely related to the change in meltwater, which reaches its peak in May, thus leading to a peak discharge of approximately 30,000 m³⋅s⁻¹ in June (Figure 4). Thereafter, the monthly discharge declines as the amount of meltwater decreases. This decline continues in the autumn, in contrast to the relatively high rainfall in late summer to early autumn, thereby indicating a less dominant impact of rainfall on the seasonal variability of discharge compared to the effect of temperature. During the melting season (i.e., from April to August), meltwater accounts for a significant proportion (77%) of the discharge (Figure 5). Additionally, during this period, the meltwater experiences an increasing contribution (from 40% in April to 60% in August) from glacierized areas. Notably, this growing contribution against the overall decrease in meltwater from June to August suggests a possible increase in the amount of ice melt in glacierized areas, in agreement with the results of Beamer et al. [28]. The decrease in discharge continues at a drastic rate (4135 m³⋅s⁻¹⋅month⁻¹) from September to November, with meltwater contributing only 5% to the discharge in November. Finally, the monthly discharge in December returns to its baseflow level as the temperature drops below 0 °C.

While no statistically significant trend is observed in the annual discharge during 1982–2022, the results in Figure 6a indicate a decadal pattern. Moreover, a statistically significant correlation is found between the five-year running mean of the NPGO index and the Alaska discharge when a two-year time lag between them is considered (p < 0.05, r² = 0.202; Figure 6b). This implies that the decadal variation in the Alaska discharge is influenced by long-term variations in the NPO. As noted above, the surface wind variability associated with the NPO is a crucial driving force for the NPGO, thereby suggesting that the NPGO is regarded as the oceanic expression of the NPO [70]. The NPO consists of two typical phases, namely the Aleutian below phase and the Aleutian above phase [17,18], and can therefore be represented by the NPGO index. Specifically, positive phases of the NPGO index correspond to the Aleutian below phase of the NPO, and are characterized by an enhanced Aleutian low pressure over a relatively wide region, including coastal and interior North America, along with an enhanced North Pacific High across the subtropical North Pacific that develops and shifts poleward. This atmospheric variability leads to intensified westerlies that subsequently blow over the central Pacific, thereby resulting in increased storminess along with higher precipitation and milder winter temperatures along the west coast of North America [18,70]. Therefore, the hydrological processes may be subsequently influenced. Furthermore, the present study indicates that the NPO impacts the Alaska discharge with a lag of approximately two years. This time lag is deemed reasonable, as sea level pressure anomalies in the winter season primarily dominate precipitation in the form of snow, which may then influence melt runoff in the following seasons.

In addition, the monthly Alaska discharge exhibits significant interannual variations and significant increases of 65.52 and 63.71 m³⋅s⁻¹⋅yr⁻¹ in summer (June–August) and autumn (September–November), respectively, during 1982–2022 (p < 0.05; Figure 7). These results can be attributed to the increases in winter snowfall and summer temperatures across the Alaska catchment. Specifically, it has been reported that the winter precipitation increased by roughly 4.9%⋅decade⁻¹, and the summer temperature increased by around 0.42 °C⋅decade⁻¹, during 1957–2021 [38].

3.3. The Discharge of the Five Basins and the Ungauged Basins

Figure 8a,b show the total discharge from the five basins and the specific meltwater contribution, respectively, while the corresponding results for the ungauged basins are shown in Figure 8c,d. Thus, during 1982–2022, the total annual discharge from the five basins was approximately 5848 ± 533 m³⋅s⁻¹ (0.03 m³⋅s⁻¹⋅km⁻²), including a substantial contribution from meltwater (58%). Specifically, the freezing temperatures from January to March resulted in minimal meltwater (138 m³⋅s⁻¹ on average) and low rainfall (34 mm on average), thus leading to an average baseflow of 701 m³⋅s⁻¹. The discharge experienced a drastic increase from 1495 m³⋅s⁻¹ in April to 9216 m³⋅s⁻¹ in May, and reached a peak discharge of approximately 15,000 m³⋅s⁻¹ in June. From June to August, although the rainfall increased, the discharge of the five basins declined by approximately 1749 m³⋅s⁻¹⋅month⁻¹ in response to a rapid decline in meltwater. During the melting season (from April to August), the monthly discharge was approximately 9970 m³⋅s⁻¹, with a considerable contribution of approximately 78% from meltwater, 36% of which came from glacierized areas. From September to December, the discharge decreased to a baseflow level at a rate of 2309 m³⋅s⁻¹⋅month⁻¹ as the temperature declined to below 0 °C. While no statistically significant trend is detected in the total annual discharge from the five basins between 1982 and 2022, the total monthly discharge in January, August, September, and November experienced a significant increase (Table S3).

Moreover, the ungauged basins (Figure 8c,d) deliver a larger contribution of 59% to the Alaska discharge, compared to 41% from the five basins (Figure 8a,b), despite having a smaller catchment area. This finding is in agreement with the results reported by Wang et al. [15]. The annual discharge of the ungauged basins is approximately 8548 ± 481 m³⋅s⁻¹ (0.07 m³⋅s⁻¹⋅km⁻²) during 1982–2022, with meltwater contributing approximately 50%. Specifically, the discharge in the ungauged basins from January to March is five times that of the five major basins. This is primarily because the temperature in the ungauged basins fluctuates more closely around 0 °C compared to that in the five basins, thereby resulting in more meltwater (1426 m³⋅s⁻¹ on average) and rainfall (150 mm on average). Similarly to that from the five basins, the discharge from the ungauged basins begins to increase substantially in April as the temperature rises above 0 °C. However, the rate of increase from April to May is slower than that observed in the five basins, primarily due to the less drastic increase in temperature in the ungauged basins, which results in a smaller variation in meltwater. The discharge of the ungauged basins also reaches its peak in June, but at a lower rate than that in the five basins, due to a slight decrease in snowmelt during this period. From June to August, despite a slightly larger rainfall, the decrease in discharge continues at a rate of 810 m³⋅s⁻¹⋅month⁻¹. During the melting season (from April to August), the monthly discharge of the ungauged basins is approximately 11,396 m³⋅s⁻¹, and again includes a larger contribution (76%) from meltwater than that received by the five basins. Of this contribution, 39% comes from glacierized areas. From August to October, the monthly discharge of the ungauged basins is approximately 11,627 m³⋅s⁻¹, which is higher than that of the five basins due to the correspondingly higher rainfall and meltwater. The temperature drops below 0 °C from October to December, but remains higher than that in the five basins, thereby leading to relatively higher amounts of meltwater and, hence, higher baseflow. In addition, during this period, the higher rainfall contributes significantly to the higher baseflow in the ungauged basins relative to the five basins. While no statistically significant trend is detected in the annual discharge from 1982 to 2022, the monthly discharge from the ungauged basins exhibits a significant decreasing trend in May, and an increasing trend in August (Table S3).

Finally, the discharges from each of the five individual basins between 1982 and 2022 are shown in Figure 9, where differences in the quantity and variation are observed. The annual specific discharges from the Stikine, Taku, Susitna, Alsek, and Copper basins are 0.036, 0.031, 0.03, 0.029, and 0.027 m³⋅s⁻¹⋅km⁻², respectively. Thus, the specific discharge in the southern portion (Stikine and Taku basins) of the Alaska catchment is generally higher than that in the northern part (Susitna, Alsek, and Copper basins). This is predominantly due to a higher annual precipitation in the southern part. Additionally, the Stikine and Taku basins receive a higher contribution from groundwater (approximately 600 mm⋅yr⁻¹) compared to the northern part. The lowest specific discharge is observed in the Copper basin, which may be due to the limited precipitation (approximately 1800 mm⋅yr⁻¹) and low average temperature (below −2 °C).

In terms of the trends, only the Susitna basin exhibits an increase in the annual discharge during the period of 1982 to 2022. Further, the Copper, Alsek, and Taku basins each exhibit a decrease in the monthly discharge in April, while all basins except the Stikine basin exhibit an increasing trend during August (Table S4).

4. Discussion

The freshwater discharge from the Alaska catchment may deliver cascading impacts on the neighboring ocean sector. After entering the GOA, the Alaska discharge primarily joins the Alaska Coastal Current (ACC). The ACC begins in the vicinity of the Washington–British Columbia border and travels at a speed of about 0.1 m⋅s⁻¹ along the coast, carrying precipitation over the ocean and continental discharge into the Bering Sea [72,73]. The ACC flows westward and southward along the Alaska coast and partially enters the Alaskan Stream before entering the Bering Sea through Unimak Pass [74,75]. As shown in Figure 10, the interannual variation of sea surface salinity in the Alaskan Stream during the period of 1984–2008 exhibits a statistical correlation with the variation in the annual Alaska discharge when considering a time lag, thereby suggesting that the Alaska discharge may impact the Alaskan Stream with a delay of two or three months (p < 0.05, r² = 0.2). Although these time lags are derived from statistical results, they are reasonable considering that the average distance covered by the ACC in going from the coastal regions of the GOA to the location of the Alaskan Stream is more than 500 km. In addition to the freshening effect of the continental discharge, other factors such as upwelling and current advection may also affect the seawater salinity of the Alaskan Stream [76], thereby rendering the correlation less prominent. This significant correlation indicates that the Alaska discharge may influence the current circulation of the North Pacific Ocean through its impacts on the characteristics of the surrounding ocean.

5. Conclusions

The Alaska discharge is characterized by large quantity and considerable variability, due to the complex interplay between climatic variations and mountainous topography, delivering cascading impacts on the coastal circulation beyond the regional scale. The present study estimated the quantity and variability of the continental discharge from the Alaska catchments under the continual fluctuations in temperature and precipitation in recent decades. The Alaska discharge was effectively estimated by using a process-based hydrological model known as the Soil and Water Assessment Tool (SWAT) with modified melting processes. During the estimation, SWAT was able to account for the complex spatial variability of land cover and soil features at the basin scale. The primary advances of this study are as follows: (i) an update on the estimated discharge during recent years, (ii) an analysis of the temporospatial variabilities in discharge, and (iii) a first discussion on the relationship between continental discharge and large-scale oceanic circulation in this region.

The results indicate an annual Alaska discharge of 14,396 ± 819 m³⋅s⁻¹, with a strong seasonal variability. In detail, the discharge maintains a baseflow level from January to March, when the temperature is well below freezing, and begins to increase drastically in April in response to rising temperatures. A peak flow of approximately 30,000 m³⋅s⁻¹ occurs in June, after which the discharge decreases swiftly despite the relatively high rainfall during late summer to early autumn, and resumes its baseflow level in December. Meltwater plays an important role in the discharge accumulation process, contributing approximately 53% to the annual discharge and exerting a more dominant control than rainfall on the seasonal discharge variation. Notably, the ungauged basins provide a larger contribution of roughly 59% to the Alaska discharge, despite having a smaller catchment area, whereas the five basins provide the remaining 41%. The interannual variation in the Alaska discharge experiences an increasing trend for the summer and autumn, while the annual discharge exhibits no significant trend. The interannual variation in the Alaska discharge from 1982 to 2022 exhibits a decadal trend. Additionally, a significant correlation is found between the 5-year running mean of the Alaska discharge and the NPGO index, thereby suggesting that the NPO, a large-scale sea level pressure anomaly variation over the North Pacific, may affect the decadal variations in the Alaska discharge.

This study revealed an unprecedented significant correlation between the interannual variation in the Alaska discharge and sea surface salinity in the Alaskan Stream. The salinity anomalies potentially impact the overturning circulation through the current pathway, since it is situated in the upper branch of the pathway. The correlation becomes significant when a time lag of two or three months is assumed, which is reasonable considering the distance between the Alaska coastal outlets and the Alaskan Stream, along with the traveling speed of sea currents. Therefore, it is concluded that the regional discharge from the Alaska catchment could impact the variations in seawater anomalies in the North Pacific sector. Overall, the present study investigated the connection between large-scale climatic variations, Alaska discharge variations, and oceanic processes, which is expected to provide a broader understanding of the atmosphere–land–ocean relationship, especially in recent years.

Nevertheless, the present study has two major limitations. Firstly, while the temperature-index approach used to calculate meltwater is efficient, it is fairly simple. This method can be regarded as adequate for the monthly simulation presented herein, but a more complex method with high spatial resolution should be applied for the glacierized area. Secondly, parameterization might be improved based on more observation records in the calibration process for further improvement.