Modelling mineral dust emissions from proglacial valleys of the St. Elias Mountains, Canada

Proglacial valleys of western Canada and Alaska demonstrate extensive historical and contemporary records of mineral dust emissions. These contributions remain unresolved by current dust emission modelling and unaccounted for in global emission estimates. We developed and evaluated a sub-km implementation of the Weather Research and Forecasting model with Chemistry (WRF-Chem) capable of simulating dust emissions from proglacial valleys of the St. Elias Mountains, Canada. Modelling these dust sources required precise treatment of surface characteristics and wind dynamics to accurately resolve surface erodibility, emission rates and aerosol dispersion within this mountainous terrain. Land-surface inputs were overhauled, with explicit treatment of glaciofluvial deposit heterogeneity and inundation conditions. Simulations covering 5–19 d periods across 2019–2022 were evaluated against in situ meteorological and dust emission measurements, camera stations and surface-based Doppler LiDAR data. A total emission rate of 1.0×10⁴ kg km⁻² d⁻¹ was estimated from erodible deposits across 47 d of simulation, providing a first-order estimate of emissions from these valleys during dusty periods. Seasonal-dependent skill in reproducing surface meteorology and in-valley vertical dispersion is demonstrated, modifying dust dispersion. Emission dynamics from a variety of glaciofluvial deposits were successfully reproduced, however the sensitivity of the emission scheme to soil texture is discussed in light of glaciofluvial deposit heterogeneity and dataset scarcity. The successful model implementation under extreme topographic conditions and arguably the most severely constrained deposits (channel width: 0.1–3 km; sidewalls up to +1.7 km) supports the potential of this approach to simulate dust emissions from currently unaccounted for proglacial valleys across northwest North America and other regions.

1 Introduction

Mineral dust constitutes a key aerosol in the climate system (Kok et al., 2023; Shao et al., 2011a). Airborne mineral dust can alter cloud properties and microphysics (Creamean et al., 2013; Karydis et al., 2017; Nenes et al., 2014; Rosenfeld et al., 2014), modify the atmospheric radiative balance (Choobari et al., 2014; Highwood and Ryder, 2014; Miller et al., 2014) and influence biogeochemical cycles (Maher et al., 2010; Ravi et al., 2011). Mineral dust deposition on snow or ice surfaces can reduce surface albedo, enhancing melting rates (Clarke and Noone, 2007; McCutcheon et al., 2021; Nagorski et al., 2019). Proglacial landscapes are key sources of contemporary high-latitude (≥50° N and ≥40° S) dust emissions (Arnalds et al., 2016; Baddock et al., 2024; Bullard et al., 2016, 2023; Meinander et al., 2022). These glacigenic sources frequently remain coupled to active fluvial systems, resulting in coinciding aeolian, glacial and fluvial processes governing dust emissions (Bullard, 2013; Crusius et al., 2011). The contribution of emissions from these sources to regional climate modification is becoming increasingly established (Barr et al., 2023; Groot Zwaaftink et al., 2016; Koffman et al., 2021), however they remain excluded from almost all dust emission modelling efforts (Crusius, 2021; Kok et al., 2021a).

The coastal mountains of southern Alaska and western Canada exhibit some of the highest rates of glacier erosion and sediment evacuation globally (Hallet et al., 1996; Sheaf et al., 2003), sourcing extensive aeolian deposits across Alaska (Bettis et al., 2003; Black, 1951; Pewe, 1975). Contemporary dust emissions, previously only evidenced in grey literature (e.g. Moffit, 1938; Trainer, 1961), have now been identified in many of the region's proglacial valleys by loess studies (Muhs et al., 2013, 2016), remote sensing approaches (Huck et al., 2023; Sayedain et al., 2023; Schroth et al., 2017) and in situ investigations (Bachelder et al., 2019; Crusius et al., 2011; Nickling, 1978). Studies indicate the ubiquitous presence of locally-sourced dust on snowpacks at lower elevation sites in the St. Elias Mountains (2.0–2.8 km a.s.l) (Zdanowicz et al., 2006) and southern Alaska (∼2.1 km a.s.l) (Koffman et al., 2022). Dust loading at higher elevation sites (3–4 km a.s.l) is primarily of Asian provenance, though locally-sourced dust is present (Koffman et al., 2022; Zdanowicz et al., 2006). The impacts of aerosol deposition are notably absent from glacier melt projections (Nagorski et al., 2019), hence the impacts of local dust deposition in enhancing regional glacier melting are currently unknown. Current high-resolution modelling efforts covering the region (Dowell et al., 2022; Milbrandt et al., 2016) are underprepared to generate representative dust emissions from these sources, lacking appropriate land-surface treatment. While recent work has included dispersion modelling from coastal outwash deposits (Barr et al., 2023), no targeted effort to simulate dust emissions from these proglacial valleys has been performed.

Multiple challenges need to be overcome to explicitly resolve dust emissions from proglacial valleys of the St. Elias Mountains. A regional-scale modelling approach is necessary to delineate the confined glaciofluvial deposits sourcing dust emissions (Bellamy et al., 2025b). Knowledge of deposit sediment texture, fundamental to simulating dust emissions (Shao, 2008), is wholly absent, continuing a prevalent insufficiency of soil characteristic datasets in aeolian studies and emission modelling efforts (Darmenova et al., 2009; Kok et al., 2014b; Laurent et al., 2006; Meinander et al., 2022; Zender et al., 2003). Sedimentary facies within proglacial braided channels are a complex mosaic of cobbles, gravels and fine deposits, dependent on time-variant hydrological conditions and sediment supply, and susceptible to frequent reworking, erosion and deposition (Maizels, 1993; Marren, 2005; Rust, 1972). Deposit configurations are further dependent on topographical-constraints and the presence of proglacial lakes modifying catchment connectivity and acting as sediment sinks with varying depositional behaviour (Bennett et al., 2002; Carrivick and Tweed, 2013; Schiefer and Gilbert, 2008). No simple model exists to ascertain channel bed lithologies in proglacial river channels. Situated in proglacial landscapes, diurnal–seasonal variability in water level and seasonal snow cover (Bullard, 2013; Marren, 2005) join common surface-limiting factors of the aeolian system, such as soil moisture, vegetation cover and aerodynamic roughness (McKenna Neuman, 1993; Wiggs et al., 2004; Wolfe and Nickling, 1993).

The mountainous terrain bounding these deposits further complicates attempts to simulate representative surface winds. Atmospheric processes over complex terrain are challenging to simulate (Chow et al., 2019; Rohanizadegan et al., 2023; Serafin et al., 2018), and local mountain wind systems (Zardi and Whiteman, 2013) have been associated with strong wind events in this region (Bellamy et al., 2025a; Crusius et al., 2024). Post-emission, the dispersion and eventual fate of suspended material is governed by exchange processes in complex terrain, determining the total consequences of dust emissions (Baddock et al., 2017; Uno et al., 2009). Model capacity to reproduce vertical exchange mechanisms and convective boundary layer (CBL) dynamics (Serafin et al., 2018) is key to accurately reproducing emission trajectories and subsequently understand aerosol impacts.

In this work, we develop and evaluate a high-resolution, i.e. valley-resolving, modelling approach capable of simulating dust emissions from proglacial valleys of the St. Elias Mountains, Canada. We detail land-surface data preparation and modifications to existing dust emission schemes to enable representation of confined glaciofluvial surfaces. We target high emission periods from late-spring to early autumn with three simulation periods to evaluate model surface meteorology, dust emission activity and transport processes against a comprehensive array of ground-based observations. We present emission estimates for the region and argue for the applicability of this approach to investigate other glacigenic dust sources across the region.

2 Methodology

2.1 Study region

This investigation focused on several proglacial valleys of the southeast St. Elias Mountains situated in Yukon, Canada (Fig. 1). The region exhibits a sub-arctic continental climate with a mean annual precipitation of 265 mm (1991–2020) (ECCC, 2023), despite close proximity (∼170 km) to the maritime climate of southern Alaska (3000–4000 mm) (Wendler et al., 2017). Steep-sided valleys range from 23–45 km long, channel widths of 0.1–3.0 km, and sidewall heights up to 1.3–1.7 km above the valley floor (Bellamy et al., 2025a). Outwash deposits within these proglacial valleys are the source of frequent dust emissions (Bellamy et al., 2025b).

Figure 1 Overview of the study region and location of the two model domains (d01; d02) (Source: Maxar, Esri, USGS | Powered by Esri).

2.2 WRF-Chem

To simulate dust emissions, we employed the Weather Research and Forecasting coupled with Chemistry (WRF-Chem v.4.6.0) model, a non-hydrostatic atmospheric model widely employed for regional-scale numerical weather prediction and research (Grell et al., 2005; Skamarock et al., 2021). With several dust emission schemes implemented, this model is frequently applied in aeolian research (Darmenova et al., 2009; Hao et al., 2024; Kang et al., 2011; LeGrand et al., 2019). Two domains, covering several proglacial valleys in the study region are illustrated in Fig. 1. The two principal challenges to simulating dust emissions from these sources were: (1) delineating erodible glaciofluvial deposits; (2) reproducing complex terrain meteorology for emissions and aerosol dispersion. A modelling approach designed to overcome these challenges is outlined in the following sections, detailing run-time settings, land-surface dataset preparation and the dust emission scheme.

2.2.1 Run-time settings and parameterisations

The model run-time settings and parameterisations used in simulations are detailed in Table 1 and 2, respectively. Simulations were performed covering three ∼ 1–3 week-long periods (s01–03) between 2019 and 2022, reflecting periods of concurrent in situ data availability and major dust activity (detailed in Sect. 2.3). Selected study periods covered principal emission seasons, from late spring–summer–early autumn conditions. Simulation 1 (s01) comprised two short periods within May 2019, where intermittent data availability and emission events intersected (see Supplement 2). Domain d02 (Fig. 1) was only employed in a separate 2022 simulation, corresponding to in situ data availability. Meteorological input was sourced from the High-Resolution Rapid Refresh–Alaska (HRRR-AK) product, a 3 km implementation of the WRF model with real-time data assimilation and hourly forecasts updated every 3 h (Dowell et al., 2022). Hourly input files were compiled from forecast timesteps of 0 to 2 h (Blaylock, 2022), with timesteps of 3–5 h used when the initial forecast file was not available. The application of HRRR-AK output to initialise simulations permitted a cost-efficient, single domain set-up with a horizontal resolution of 600 m: sufficient to resolve both complex topography and erodible deposits. To resolve issues with the turbulence grey zone at this resolution (Chow et al., 2019), we employed the Shin and Hong (2015) “scale-aware” CBL scheme, which accounts for resolved and parameterised turbulence as a function of model resolution. The effectiveness of these scale-aware schemes for parameterising turbulence within the grey zone remains under review (Doubrawa and Muñoz-Esparza, 2020; Giani et al., 2022; Honnert et al., 2020), however they constitute the best approach for studies at this resolution to date. Spin-up periods of 24 h were used, with no chemistry initialised.

Table 1 Run-time settings for WRF-Chem (v4.6.0) simulations. For namelist options see the WRF Users' Guide.

Table 2 Parameterisation for WRF-Chem (v4.6.0) simulations.

2.2.2 Preparing model input

Substantial improvements to default land-surface datasets were required for this implementation (Table 3). Modifications to land-use, soil classification and the erodible source function are further detailed in this section.

Table 3 Updated input data for simulations.

A binary erodible source factor was employed to constrain erodible glaciofluvial deposits, identifying surfaces prone to frequent seasonal inundation and exposure, inferring fluvial deposition of glacigenic sediment. Surface water mapping was performed with Landsat 5–8 retrievals following Feyisa et al. (2014), with river conditions and exposed outwash deposits constrained for each simulation date. Following Bellamy et al. (2025b), potential erodible area was identified as surfaces inundated during summer and subsequently exposed in autumn a minimum of 5 years across 40 years of available retrievals (1984–2023). The criteria identified regularly reworked surfaces only but was sufficiently lenient to capture conditions following a major hydrological reorganisation in 2016 (Shugar et al., 2017). The erodible source factor is subsequently defined for each simulation period as potential erodible area exposed at the time of the simulation, incorporating date-specific water cover conditions.

No datasets were available that provide an adequate representation of outwash deposit characteristics. To overcome this, four outwash deposit classifications were added to the land surface scheme: very coarse; coarse; fine; very fine (Table 4). We link these classifications to similar lithofacies groupings in Maizels (1993). Each group was prescribed a representative sediment texture which were used to update model soil inputs and soil descriptions in the emission scheme. Representative examples of each classification are provided in Fig. 2, while a domain overview is available in Fig. S1 in the Supplement.

Figure 2 Overview of the modified Shao (2004) scheme in WRF-Chem developed in this work. Modifications were applied to the drag partition parameterisation (${\mathit{\lambda }}_r,u_{*}\left(r\right)$), further detailed in Sect. 2.2.3.

Table 4 Additional land use classifications for outwash plains and prescribed sediment textures. Description indicates corresponding geomorphological unit classifications in Maizels (1993). D, h are the diameter and height of particles, respectively.

Outwash deposit classification was performed based on field observations (Bellamy et al., 2026) and surficial geological maps (Rampton, 1977; Rampton and Paradis, 1982). In lieu of detailed information of surface deposits in the upper Kaskawulsh river valley (60.82° N, 138.51° W to 60.69° N, 137.97° W), all outwash units are characterised as coarse. Representative soil texture classifications were derived from multiple sediment samples from each deposit grouping, collected within the A'ą̈y Chù valley (see Supplement Fig. S2). For scheme consistency, minimally and fully disturbed volume size distributions were obtained following the procedure of Shao et al. (2011b) using a MasterSizer 3000 (Malvern Instruments) and parameterised as lognormal distributions (supplementary 1). The newly implemented particle-size distributions (PSDs) for sandy loam, sand and silt soil classes increased the fully-disturbed ≤20 µm mass from 1.6 %–14.7 %, <0.1 %–9.8 %, and 8.8 %–64.3 %, respectively, compared to existing soil descriptions for the Shao (2004) scheme in WRF-Chem (see Supplement Fig. S3). Roughness lengths (z₀) were attributed to each new land use classification based on wind-tunnel investigations of aerodynamic roughness for gravel and sand surfaces (Table 4) (Dong et al., 2001, 2002).

Preliminary modelling and sensitivity testing for s01 indicated a significant overestimation of reflected shortwave radiation on the valley floor, inducing an underestimated surface energy budget with influence on dust emission estimates (Supplement 3). Accordingly, the albedo parameter for the four glaciofluvial deposits was set to 0.14 (compared to 0.38 for “barren or sparsely vegetated”) based on the average albedo at solar noon as measured at the surface meteorological station DV (Sect. 2.3) across all simulation periods. This adaptation improved outgoing radiative fluxes (Supplement 3), however replacing the single value does not improve on domain-wide albedo discrepancies due to variable surface moisture or deposit type. Remaining characteristics required by the unified Noah land-surface model (LSM) in WRF-Chem remain identical to the “barren or sparsely vegetated” classification.

2.2.3 Dust emission scheme

We employed the University of Cologne (UoC) Shao (2004) scheme as implemented in WRF-Chem v4.6.0 to simulate dust emissions. While the dust emission flux is predicated on Shao (2004), the complete scheme draws from Shao and Lu (2000), Shao et al. (2011b) and references herein. This scheme determines threshold friction velocity (u_*t) for saltation, calculates saltation flux across a range of saltating particle sizes, and subsequently evaluates the dust emission flux and passes dust concentrations in five particle size bins (0.2–2.0, 2.0–3.6, 3.6–6.0, 6.0–12.0, 12.0–20.0 µm) to the wider WRF-Chem model for transport. Employed extensively for dust emission modelling (Darmenova et al., 2009; Edwards et al., 2022; Klose et al., 2021; LeGrand et al., 2019; Park et al., 2007) this scheme is transcribed in Fig. 2. Modifications to enhance the treatment of surface erodibility in our study region were implemented and are detailed below.

Shao (2004) incorporates the variable soil plastic pressure (P) (also referred to as soil penetrometer resistance) to represent soil binding strength in determining the mass ejected by saltation bombardments (Lu and Shao, 1999). This term is currently held constant across all soil types despite being soil texture and time-variant (Zobeck, 1991), as noted by LeGrand et al. (2019). P is currently estimated by model fitting (Klose et al., 2019; Shao, 2004) as appropriate penetration resistance measurements (e.g. Rice et al., 1997) are not widely reported. Based on the estimates available in Rice et al. (1997) we approximated P for sandy soils as 3×10³ N m⁻² and silt soils as 3×10⁴ N m⁻². Constraining P further as a function of sediment texture as pertinent for saltation is currently not possible. The dimensionless tuning coefficient (c_y), used to adjust the emission scheme for specific soil conditions (Klose et al., 2019; Shao et al., 2011b), was fixed at $c_y=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ (Shao, 2004).

Previous approaches to constrain gravel cover in drag partitioning (Klose et al., 2021; Laurent et al., 2006; Leung et al., 2023) were inapplicable at this scale of study. In lieu of this, we assigned gravel characteristics during our outwash deposit characterisation (Table 4) and derived a second drag coefficient following Raupach et al. (1993). The drag partition scheme (Fig. 2) was modified to apply an inverse correction for vegetation and gravel shear stress partitioning to u_* instead of u_*t (Kok et al., 2014a). Gravel scheme coefficients (β_g; λ_g; m_g – see Fig. 2) were set following the results of Tan et al. (2019). This effectively corresponds to a fixed decrease of u_* by a factor of 2.15 and 3.56 for the coarse and very coarse outwash deposits, respectively. Separately, the soil-moisture correction to u_*t was implemented with a modified soil moisture (0.1θ_v) to account for discrepancies between model top soil layer conditions and those of surface grains subject to aeolian processes, as frequently performed in dust emission modelling (Darmenova et al., 2009; Klose et al., 2021).

2.3 In situ observations for model evaluation

To fully test model performance, we evaluated model skill in reproducing: (1) near-surface meteorology; (2) dust emissions (location, timing and fluxes); (3) vertical mixing and aerosol transport within the mountainous terrain. Data collected during field campaigns in 2019, 2021 and 2022 are valuable in evaluating these processes, and are detailed separately in the following sections. An overview of station installations in the study region is given in Fig. 3, with instrumentation summarised in Table 5.

Figure 3 An overview of camera, MET (meteorological station) and RS (remote sensing) installations in the region. The two model domains (d01; d02) are indicated. Landsat 9 imagery (7 July 2023) courtesy of the U.S. Geological Survey.

Table 5 Station details. Data availability during years indicated as yes (Y), partial (y) or no (–).

2.3.1 Surface meteorology

Surface meteorology was evaluated across all simulation periods using three meteorological stations in the A'ą̈y Chù valley, covering late spring–early autumn conditions. Near-surface wind diagnostics were re-evaluated at sensor height for each station following the surface layer scheme (Jiménez et al., 2012) (Eq. 24):

$$\begin{array}{}\text{(1)}& u\left(z\right)=u_a{\displaystyle \frac{\ln \left(\frac{z+z_{\mathrm{0}}}{z_{\mathrm{0}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z+z_{\mathrm{0}}}{L}\right)+{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{0}}}{L}\right)}{\ln \left(\frac{z_a+z_{\mathrm{0}}}{z_{\mathrm{0}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_a+z_{\mathrm{0}}}{L}\right)+{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{0}}}{L}\right)}}\end{array}$$

where u(z) is the wind speed at diagnostic height z, u_a is wind speed at the lowest model level of height z_a, z₀ is the roughness length, L is the Obukhov length, and ψ_m is the momentum stability function. Radiative forcing and turbulent fluxes were evaluated at DV, with turbulent fluxes measured with an eddy covariance system (Aubinet et al., 2012). Due to an incomplete sonic anemometer record, friction velocity (u_*) at DV was calculated using a stability-corrected law-of-the-wall approach (e.g. Dupont et al., 2018) with wind and temperature profile measurements at DV (Table 5) over 30 min averaging periods. Linear regression fits of R²<0.9 were discarded. Limited in consistent stability measurements across the three years, we relied on a simple two criteria momentum similarity function (Bonan, 2019) to address the impact of stability on surface layer profiles:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& {\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z}{L}}\right)=\\ & \left\{\begin{array}{ll}\mathrm{2}\ln \left(\frac{\mathrm{1}+\mathit{\chi }}{\mathrm{2}}\right)+\ln \left(\frac{\mathrm{1}+{\mathit{\chi }}^{\mathrm{2}}}{\mathrm{2}}\right)-\mathrm{2}{\tan }^{-\mathrm{1}}\left(\mathit{\chi }+\frac{\mathit{\pi }}{\mathrm{2}}\right)& \frac{z}{L}<\mathrm{0}\\ -\mathrm{5}\left(\frac{z}{L}\right)& \frac{z}{L}\ge \mathrm{0}\end{array}\right.\end{array}\end{array}$$

where: $\mathit{\chi }={\left(\mathrm{1}-\mathrm{16}\left(\frac{z}{L}\right)\right)}^{\mathrm{1}/\mathrm{4}}$.

2.3.2 Boundary-layer meteorology

Vertical profiles of horizontal wind speed, wind direction and backscatter were obtained from a Halo Photonics Doppler LiDAR deployed at KLRS in May 2019 (Table 5). Significant signal attenuation limits consistent retrievals to ∼ 700–1500 m, a common limitation for locations with low background aerosol levels (Wiggs et al., 2022). We derive dust mass concentrations from LiDAR backscatter measurements following Ansmann et al. (2019), deriving a novel conversion parameter for 1.548 µm Doppler LiDAR from AERONET (Giles et al., 2019; Holben et al., 1998) retrievals at the Kluane Lake site (61.027° N 138.411° W), with reference to Sayedain et al. (2023). Further details are provided in Supplement 4. LiDAR vertical profiles and WRF vertical profiles from the closest grid point to KLRS were averaged or interpolated, respectively, to 30 m levels for direct comparison. The calculation of LiDAR diurnal averages is restricted to vertical bins with a majority of QC-passed retrievals during simulation periods.

2.3.3 Dust emission flux

The gradient method (Gillette et al., 1972) has seen extensive use in measuring the vertical dust flux (F_d) within aeolian literature (Ishizuka et al., 2014; Khalfallah et al., 2020; Nickling and Gillies, 1993; Shao et al., 2011b; Sow et al., 2009). F_d ($\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) was calculated using 30 min averaged measurements from two optical particle counters situated on a tower at the DV site, at heights of 3.1 and 6.1 m during 2019 and 2021 (Table 5) such that:

$$\begin{array}{}\text{(3)}& F_{\mathrm{d}}=u_{*}k{\displaystyle \frac{C_{\mathrm{1}}-C_{\mathrm{2}}}{\ln \left(\frac{z_{\mathrm{2}}}{z_{\mathrm{1}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{2}}}{L}\right)+{\mathit{\psi }}_{\mathrm{m}}\left(\frac{z_{\mathrm{1}}}{L}\right)}}\end{array}$$

where C₁ and C₂ are dust concentration measurements at heights z₁ and z₂, ψ_m is the momentum stability function and L is the Obukhov length. ψ_m was determined following Eq. (2) (Sect. 2.3.1). We estimate and correct for the effect of gravitational settling following González-Flórez et al. (2023) to infer the surface emission flux from the gradient-method derived diffusive flux, further detailed in Bellamy et al. (2026). Size-resolved dust fluxes across 22 optical bins between 0.28–10 µm were summed to yield F_d. We converted from optical to geometric diameters (0.25–9.81 µm) following Huang et al. (2021), based on an OPC laser wavelength of λ=630 nm, 90° ± 90° and dust refractive index of 1.54–0.001i. For the direct comparison of observed vs. simulated F_d, WRF size-resolved emission fluxes were summed across the same geometric range, partitioning fluxes in size bins 1 and 4 assuming a log-normal distribution of particle volume in each bin (μ=1.1, σ=0.25; μ=9, σ=0.8, respectively). Missing data periods within the outlined simulation periods are evaluated in Supplement 2, though they are sparse (19.4 %) and predominantly limited to the early morning hours (6–10 h), the calmest period of the day.

2.3.4 Camera-derived dust observations

Camera stations permit a broader assessment of emission activity across heterogenous glaciofluvial deposits. Model capacity to simulate emissions in adjacent valleys was evaluated through camera observations and a manual binary identification of dust emission activity. A qualitative emission log was derived from regular 10 min photos obtained in 2022 by C1 and C2 to evaluate emissions simulated in the A'ą̈y Chù valley and at the Alsek confluence (Fig. 3). Camera records detail periods of visible dust suspension, examples of which are given in Fig. 4. Near contiguous periods of emissions were classified into events. Visually observed precipitation was also logged to discuss model precipitation in relation to dust emissions. Due to the qualitative nature of camera-derived observations, this assessment was limited to evaluating the reliability of simulated emission timings across the broader region. As the approach required the visual identification of suspended dust, weaker events may be absent from camera-based emission logs.

Figure 4 Examples of dust activity from deployed camera stations C1 (top) and C2 (bottom) during September 2022.

3 Results

The first subsection (Sect. 3.1) evaluates modelled surface meteorology, emission dynamics and vertical dispersion against a suite of in situ observations to analyse model performance. The second subsection (Sect. 3.2) summarises modelling output with total emission estimates and investigates dispersion characteristics. We conclude on model estimate shortcomings, in light of the evaluation results, in Sect. 4.1.

3.1 Model evaluation

3.1.1 Surface meteorology

An evaluation of meteorological variables at DV across all simulation periods is shown in Fig. 5. Simulated net radiation (Fig. 5i–l) was slightly overestimated vs. DV (mean bias of +23 W m⁻²), with daily discrepancies attributed to disparate cloud cover. The close reproduction of diurnal transitions across all simulation periods indicates good performance of topographic shading routines. When observations were available, sensible heat fluxes exhibited a minor overestimation (mean bias of +19 W m⁻²) due to persisting discrepancies from cloud-cover differences. During model spin-up (24 h – not shown), soil moisture decreased rapidly, falling 8 %–11 % at DV between simulation periods. The adjustment continued in day 1 of s01 (Fig. 5e), suggesting insufficient spin-up time for soil moisture adjustment. This may have contributed to the subdued dynamic instability during 7 May (Fig. 5m), inhibiting surface wind speeds and dust emissions. Simulated soil moisture content remained overestimated across all simulation periods. With simulations focused on major emission periods, model precipitation was not comprehensively evaluated but is discussed briefly in Sect. 3.1.3. Simulated winds (Fig. 5a–d) closely followed observations across simulation periods, with a tendency to underestimate (mean bias of −1.24 m s⁻¹), and performing best in reproducing winds at DV during s02 (June–July 2021).

Figure 5 Comparison of hourly mean observed vs. simulated meteorological variables at DV across all simulation periods for: (a–d) wind speed (4.5 m), (e–h) soil volumetric moisture content (VMC) at 10 cm (DV) vs. 0–10 cm (WRF), (i–l) net radiation and (m–p) sensible heat flux (H). Period-wise bias ($\stackrel{\mathrm{‾}}{e}$) for each variable is indicated. All results presented in local time (UTC −7).

An evaluation of diurnal wind speed distribution and bias in Fig. 6 indicated varying model skill to simulate surface winds between study periods and locations. Wind speed distributions for s02 (Fig. 6b, d, f) compared the best, with a mean bias of −0.03 and −0.19 m s⁻¹ for MV and UV, respectively. Estimates at DV for s02 are underestimated (mean bias of −1.25 m s⁻¹), a feature persisting across all simulations at this station (Fig. 6a–c). During s03 (Fig. 6c, e, g) poorer model performance is evident (mean bias of −1.58 to −1.86 m s⁻¹), underestimating surface winds despite covering a period subject to intense synoptic forcing (see Supplement 6). Notable diurnal bias at MV and UV during s02 and s03 (Fig. 6i–n) suggests that model capacity to reproduce valley CBL dynamics may have underlain wind speed discrepancies and induced seasonal variability in model ability to capture surface forcing.

Figure 6 (a–g) Observed and simulated wind speed distributions at respective sensor height (z) for each simulation periods (s01: 2019 – 12 d; s02: 2021 – 11 d; s03: 2022 – 19 d). The s02 comparison only compares records from 24 June 2021–3 July 2021 due to station data availability. (h–n) WRF wind speed bias vs. surface observations, shaded interval indicates 1σ.

3.1.2 Vertical dust fluxes on the A'ą̈y Chù delta

Friction velocity (u_*) and the vertical dust flux (F_d) (0.25–9.81 µm) were evaluated against measurements at DV during s01 and s02, as presented in Fig. 7. The underestimation of simulated surface winds (previous section) is further apparent in the u_* comparison, with a mean bias of −0.08 m s⁻¹ at DV. Despite this, of 27 d of observed dust emission events on the A'ą̈y Chù delta during this period, all were simulated in WRF. Across both simulation periods, we measured a total vertical dust flux of 0.166 kg m⁻² at DV, with a daily average of 5.9×10³ kg m⁻² d⁻¹, notwithstanding observational gaps. In contrast, the closest grid point to DV simulated a total vertical dust flux of $\mathrm{6.2}\times {\mathrm{10}}^{-\mathrm{3}}$ kg m⁻², approximately 4 % of that measured. In Fig. 7, simulated values of F_d at DV (soil class: sand) and F_d averaged over emissive grid cells of the delta (soil class: 56 % silt; 44 % sand) are illustrated separately to examine sediment texture sensitivity in the Shao (2004) scheme. Simulated F_d at DV was substantially underestimated (mean bias of −86 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) however fluxes averaged across the delta compare more favourably (mean bias of +26 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), notwithstanding u_* discrepancies. Below we highlight the causes of emission flux variability across the delta and the main factors underlying discrepancies at the point-scale comparison.

Figure 7 Comparison of measured (DV) and WRF-simulated variables: at the DV grid point (WRF.DV) and averaged over emissive grid cells of the A'ą̈y Chù delta (WRF.Delta). (a–c) Friction velocity (u_*). (d–f) Vertical dust flux(F_d) (0.25–9.81 µm).

The impact of soil moisture varied across the delta and influenced F_d substantially at low u_*–u_*t values. Sediment texture heterogeneity, combined with simulated θ_v overestimates (Fig. 5e–h) and the texture-dependent soil moisture corrections of Fécan et al. (1999), also yielded moisture-induced spatial discrepancies in u_*t. Across all simulation periods (47 d) fine outwash deposits (sand) experienced moisture-induced heightening of u_*t for 45 d, despite the modified soil moistures used in dust emission parameterisation (Fig. 2), while very fine deposits (silt) remained unaffected. For s01, with no rain and dry conditions, a soil moisture-induced u_*t modification of +0 % to +40 % was effected at DV (average +11 %) under the modified emission scheme (Fig. 2), despite dry surface conditions being observed. Under this effect, u_*t ranged from 0.21–0.29 m s⁻¹. For very low values of u_*–u_*t (<0.1 m s⁻¹) this can induce F_d reductions of 10¹–10⁴ vs. dry soil u_*t: the source of the extremely low values for the DV comparison in Fig. 7d–f. While minor precipitation events did occur during s02 on the 21, 22 and 30 June, we note an average soil moisture correction of +16 % was effected for the simulation duration.

For higher values of u_*–u_*t, sediment texture prescribed to the two-outwash classes present on the delta becomes the key control on variability, with differing sub-20 µm particle mass and variable dust production as a function of u_*–u_*t, following the aggregate disintegration scheme γ (Shao, 2001). Under ideal conditions at low u_*–u_*t, minimally-disturbed PSD definitions account for F_d discrepancies on the order of 10¹–10². For silt outwash, between minimally- and fully-disturbed PSDS, sub-20 µm particle mass only increase ∼30 % (Supplement Fig. S3), yielding limited influence on F_d as a function of u_*–u_*t. For sand outwash deposits as defined in this study (Table 4), sub-20 µm particle mass increased ∼1800 % for the same fixed u_* range as defined by γ, with 95 % of the 9.75 % ≤20 µm mass aggregated as saltation-sized particles, only released in the fully-disturbed PSD. The spatial differences in emission across the delta narrow with increasing u_*–u_*t.

Concerning the grid point comparison at DV, the identified soil moisture adjustments to u_*t are insufficient to account for persistent F_d underestimation during major dust events in high u_* conditions. During the peak emission in s02 (24 June 2021 18:00), observations provide $u_{*}=\mathrm{0.685}$ m s⁻¹, F_d=2033 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, while simulated u_* at DV was 0.301 m s⁻¹ with F_d=1.16 $\mathrm{\mu }\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. A 1D implementation of S04 suggests bringing u_* in-line with observations would increase F_d by ∼10², insufficient to match observations. The implemented κ=1, representing a relatively resistant soil to disaggregation (Shao et al., 2011b), may have been too low for these sand deposits and simulated saltation failed to release fine-grained materials as observed. However, field measurements of aggregate stability (κ) are limited and time-variant (Klose et al., 2019), challenging the implementation of soil-specific coefficients in modelling efforts. The impact of underestimated u_* on simulated fluxes at DV is exacerbated by the aggregate disintegration efficacy in this scheme. Separately, we note that in applying stress partitioning correction for u_*(r) instead of u_*t (Fig. 2), the aggregate disintegration mechanism applied over an extended range of u_*, for the very coarse and coarse outwash deposits subject to enhanced stress partitioning. Sediment textures were measured during 2022 (see Supplement 1), however with outwash deposits prone to rapid reworking by active sands and exposing finer sediment layers, fine sediment exposure may have been higher during 2019 and 2021 (Bellamy et al., 2026), and fine content may be underestimated in the average soil PSD.

3.1.3 Adjacent valley emissions

Having evaluated emission fluxes against point measurements at DV, we consider model skill in reproducing dust events from other valleys. In Fig.  8, we compare simulated emissions against a timelapse photography-derived log of dust emissions, providing a qualitative record of emission activity at the A'ą̈y Chù delta and the Alsek confluence (∼54 km away) during s03 (September 2022). This period is characterised by strong emissions early in the month (8–9 September), with a wetter mid-month period, followed by a recurrence of major emissions from 20 September onwards. Emissions within the A'ą̈y Chù were more frequent than at the Alsek confluence, with 16 vs. 9 individual events identified within the 19 d simulation period, respectively. During periods of observed dust activity at the A'ą̈y Chù delta and Alsek confluence, substantial (>100 µg m⁻³) dust concentrations were simulated for 31.3 and 14.7 h, respectively, corresponding to 47 % and 27 % of the periods. Conversely, the same dust concentration threshold was exceeded 48.5 and 24.0 h, at each respective location, without emissions being observed. This false positive simulated dust activity can be partially attributed to missing precipitation (e.g. evening of 9 September 2022). During s03, the Alsek confluence was more prone to precipitation, potentially inhibiting similar emission activity to that observed in the A'ą̈y Chù between 13 and 16 September 2022. This discrepancy between observed and simulated precipitation, and the subsequent influence on surface erodibility, was an anticipated limitation due to the complex terrain.

Figure 8 Dust observations and precipitation record derived from 10 min photography from C1 (A'ą̈y Chù) and C2 (Alsek confluence) during s03. Simulated dust concentration at 50 m was averaged for the area of the Alsek confluence covered by the viewshed of C2 (60.561–60.670° N, 137.772–137.855° W), and for the A'ą̈y Chù delta covered by C1 viewshed (60.988–61.022° N, 138.446–138.553° W). Simulated precipitation index indicates precipitation ≥0.1 mm h⁻¹ averaged over viewsheds.

3.1.4 Profile meteorology

To evaluate model skill in reproducing CBL behaviour and dust dispersion, a diurnal comparison of LiDAR retrievals and WRF near-surface profiles at KLRS is given in Fig. 9. LiDAR observations evidenced consistent south-westerly winds from the direction of the A'ą̈y Chù valley over KLRS in May 2019 (Fig. 9c), however in s01 near-surface south-easterly winds were a persistent feature (Fig. 9d). Upon inspection, simulated south-easterly winds appeared to be enhanced by a nighttime temperature contrast between the lake surface (13 °C) and immediately adjacent boreal forest with skin temperatures 4–5 °C lower. Lake temperature measurements from previous years (Crookshanks and Gilbert, 2008; McKnight et al., 2021) suggest simulated lake temperatures were overestimated for this time of year. This may pose a persistent challenge to reproducing surface meteorology in the lake vicinity and influence late-spring and early autumn dust emission periods, requiring modifications to the initialising datasets.

Figure 9 Diurnal mean profiles from LiDAR and WRF output at KLRS during s01 (May 2019) for: (a–b) horizontal wind speed; (c–d) wind direction; (e–f) vertical velocity; (g–h) dust concentration.

Dust was simulated over KLRS throughout the diurnal cycle, through nighttime to the early morning at DV (Fig. 9h). While LiDAR-derived dust concentrations indicate that dust remains concentrated near the surface (≲500 m – Fig. 9g), simulated dust was dispersed beyond 1000 m, with the largest discrepancies occurring during the night (Fig. 9h). Despite this, average w_vert values observed/simulated for the afternoon and night compare well. We note that the paucity of negative w_vert in Fig. 9e is likely associated with limited aerosol content in these stable layers. A side-by-side comparison of observed vs. simulated dust concentration profiles is given in Fig. 10. While comparison on an hourly basis emphasised substantial model deviations, the timeseries in Fig. 10 suggests broader agreement in the timing of suspended dust over KLRS. The two records indicate that the simulated vertical dispersion of emissions was regularly overestimated at KLRS during s01. Overestimated lake surface water temperatures, as previously established, may have contributed to enhanced plume dispersal. Further reference to Fig. 5d–e suggests simulated elevated dust may have originated from emissions at the glacier terminus and upper-valley, being dispersed over Vulcan Mountain or advected down-valley. Further comparisons of LiDAR and WRF retrievals are included in the Supplement 5.

Figure 10 LiDAR- and WRF-derived dust concentration profiles retrieved/simulated at KLRS during s01 (May 2019).

3.2 Model results

3.2.1 Emission maps

A summary of cumulative emissions and average dust concentrations for each simulation period is illustrated in Fig. 11. Total emissions by glaciofluvial deposit classification are summarised in Table 6. Emissions occurred across the four glaciofluvial classifications, moderated by gravel-induced shear stress partitioning for the coarser deposits (Fig. 2). Emissions from very fine deposits dominated total mass loading, accounting for 54 % of total emissions. The strongest dust emission source was the A'ą̈y Chù delta, consisting of fine and very fine deposits. Spatial discrepancies in F_d due to these deposits (i.e. sediment texture) were also evident. During simulation periods, emissions from the A'ą̈y Chù delta were predominantly carried northward by southeasterly winds in the Shakwak Trench, however some dispersion to the southeast was noted. Dust sources in the upper A'ą̈y Chù, the terminal moraine and the Kaskawulsh valley, and subsequent transport indicated by concentration plots demonstrates the potential for dust dispersal up adjoining tributaries and over the Kaskawulsh glacier, while strong emissions from fine deposits in the terminal moraine were carried over the peak of Vulcan Mountain (60.89° N, 138.46° W; ∼2.7 km a.s.l). Emissions from the terminal moraine peaked in s01, though, we are limited in our capacity to validate fluxes from these deposits. We highlight the extensive emissions simulated for coarse deposits (no other classes present) at the Alsek confluence and Dusty River valley during s03 (Fig. 11d), which due to rock-induced drag corrections applied, indicates marked erosivity. Emissions from the Alsek confluence were dispersed up-valley to the northeast, but more predominantly northward up Ferguson Creek (60.70° N, 137.89° W) or through the Dezadeash River valley (60.69° N, 137.77° W). In summary, a dust (0.2–20 µm) emission rate of 1.0×10⁴ kg km⁻² d⁻¹ was simulated for erodible surfaces across all simulation periods.

Figure 11 (top row) Summed dust (0.2–20.0 µm) emissions (kg) from 0.6 km² grid cells across simulation periods in d01 (a–c) and d02 (d). (bottom row) Mean dust concentration at 50 m during active emission periods (non-zero emissions anywhere within domain) in d01 (e–g) and d02 (h). Percentage of period with active emissions indicated in each figure. Contours kriged from model topography.

Table 6 Summary of emission loading (3 s.f.) by land-use classification for simulations in d01 (s01–03) and d02 (s03). Dust emission rates are calculated for the total area of outwash deposit surfaces classified as erodible (Sect. 2.2.2) – simulation dependent.

3.2.2 Vertical dispersion

We explored simulated vertical dispersion in the A'ą̈y Chù valley, including the delta where vertical transport from a strong emission source immediately contends with convective dynamics over a large proglacial lake. Figure 12 presents cross-sections of model output averaged during periods of non-zero simulated dust emissions at DV. While vector averaged wind conditions suggest nondescript synoptic forcing on s01–s02, neither period was bereft of major synoptic forcing (see Supplement 5); average conditions indicate the heterogeneity of aloft conditions under which emissions at DV occur. During s01, emissions from the glacier terminus and upper valley dispersed well above ridge height (>1.5 km a.s.l). Dust remained at high altitudes through to the lower valley, where emissions climbed slowly during and after traversing the lake (Fig. 12e, h, k). Infrequent up-valley emissions from the delta also occurred (e.g. afternoon 24 May). During s02, dust was well-dispersed within the valley (Fig. 12g–h), however emissions primarily failed to traverse the lake, instead carried in variable directions along the Shakwak Trench (Fig. 11f). During s03, vertical dispersion within the valley (Fig. 12j) was strongly restricted relative to previous periods, despite occurring under intense synoptic forcing. Emissions from the upper-valley failed to advect down-valley, in contrast to s01–s02. Upon exiting the valley, vertical dispersion occurred steadily and substantial dust concentrations persisted to the edge of the Ruby Range (Fig. 12k – km 28) and beyond. Variability in the vertical dispersion of emitted dust between simulation periods is apparent. With s01–s03 situated within late spring, summer and early autumn of their respective years, we highlight these “seasonal” differences in emission trajectories between periods.

Figure 12 Cross sections of average horizontal wind speed and dust concentrations for (d–f) s01, (g–i) s02 and (j–l) s03. Contours range from 10–200 µg m⁻³. Grey contouring indicates a parallel track of background valley sidewalls; (–) partial transparency indicates foreground topography (axis not to scale). Landsat 9 imagery (7 July 2023) courtesy of the U.S. Geological Survey.

4 Discussion

The simulations detailed in the previous section targeted coinciding periods of high dust activity and observation availability to obtain regional dust emissions estimates while permitting an evaluation of model performance. In the following sections we discuss: (1) model estimate reliability considering evaluation results and model construction; (2) insights gained regarding investigating proglacial landscapes with modelling; (3) the potential to target other known high-latitude dust sources with this approach.

4.1 Model performance: how accurate are simulated emissions?

4.1.1 Total dust estimates

We assess the model emission estimates presented in Sect. 3.2 following model evaluation findings (Sect. 3.1). Surface winds were underestimated during s01 and s03 in the A'ą̈y Chù valley, reducing simulated emissions. At DV, an underestimated u_* (Fig. 7) contributed to the notable underestimation of F_d from fine (sand) glaciofluvial deposits, however analysis in Sect. 3.1.2 indicates soil moisture overestimates prescribed sediment fine content at DV contribute to discrepancies. Eschewing spatial variability and comparing domain-wide averages, the average estimate of 1.0×10⁴ kg km⁻² d⁻¹ remains within an order of magnitude of the average vertical dust flux measured at DV during s01 and s02 of 5.9×10³ kg km⁻² d⁻¹. Substantial emissions were simulated from currently exposed proglacial lake deposits at the Kaskawulsh terminus (Fig. 11a–c) with dust transported up-glacier. Uncertainty remains in ascribing sediment characteristics to these lacustrine deposits, given their morphological complexity (Carrivick and Tweed, 2013). Sediment textures prescribed across the terminal moraine based on field observations likely do not capture full sediment heterogeneity in this geomorphologically-active proglacial environment, hence emission fluxes cannot be well validated. Nevertheless, we highlight the high erosivity simulated at these proximal dust sources. For emissions at the Alsek confluence we were limited to evaluating emission timing, however simulated emissions broadly aligned with available observations (Fig. 8), with discrepancies principally attributed to model precipitation and erodibility conditions.

4.1.2 Outwash deposit characterisation

The mandatory overhaul of land-surface descriptions required the prescription of representative aerodynamic roughness lengths and sediment textures to mapped outwash deposit classes (Table 4), two key variables controlling the spatial pattern and strength of emissions. For the classes defined in this work, gravel-induced roughness across classes reduced u_* by a factor of 1–3.56 during emission parameterisation, while the assigned sediment textures had minimally-disturbed <20 µm particle counts of 0.9 %–48.2 %. Both ranges can induce order-of-magnitude changes in the simulated dust emission flux; however, this captures the heterogeneity of deposits within glaciofluvial system we attempted to resolve (Maizels 1993). Within the A'ą̈y Chù valley, surface characteristics were well constrained with field observations and sediment sampling. Larger uncertainty exists in ascribing sediment characteristics within the Kaskawulsh and Alsek River systems. Surficial geological mapping, where available, is dated (Rampton and Paradis, 1982) and precedes major seasonal inundation experienced since a 2016 hydrological reorganization (Shugar et al., 2017) that enhanced river levels and likely extended deposits of erodible fine sediments (Bellamy et al., 2025b). Considerable uncertainty subsequently remains in the characterisation of these surfaces and the final emission estimates presented. For studies modelling dust emissions from these active river systems, uncertainty in these prescribed surface properties is likely to remain persistent.

4.1.3 Transport processes and dispersion

Beyond emissions, aerosol dispersion within complex terrain and the altitude suspended dust can attain is an important crux in determining the residence time and fate of emitted particles (Baddock et al., 2017; Knippertz et al., 2009). Our evaluation of vertical dispersion was limited to a single location immediately adjacent to a large lake, however simulated vertical dispersion appears to be overestimated (Sect. 3.1.4). Clear discrepancies were seen in simulating surface winds (Fig. 6) and aerosol dispersion (Fig. 12) within the valley between simulation periods, allocated as seasonal variability. During s03, limited aerosol mixing occurred despite intense and suitably orientated synoptic forcing. Considering the underestimated surface winds, it is likely that vertical dispersion was underestimated due to shortcomings in resolving CBL dynamics within the valley. Future HLD studies of emissions in complex terrain should not omit consideration of the dynamics underlying vertical dispersion post-emission.

This work's model domains, focused on reproducing meteorological forcing and dust emissions at targeted dust sources at relatively low computational cost, were not extensive enough to fully evaluate emission trajectories and deposition. Nevertheless, we highlight extensive northward dispersion simulated from the A'ą̈y Chù delta and the potential for up-glacier dispersion from exposed proglacial lake deposits. Diurnal LiDAR comparison at KLRS (Fig. 9) suggested that simulations overestimated early morning and afternoon-evening vertical dispersion. However, even at low altitudes (Fig.  11), we highlight significant dispersion at least 15 km up-glacier and over 24 km along the Shakwak Trench to the model domain edge. A dispersion model could be quickly applied, benefitting from improved WRF meteorological fields, to interrogate wider transport pathways.

4.2 Lessons learned for initialising proglacial dust emission modelling

4.2.1 Land-use

A focus of this study was the production of detailed land-surface datasets from remote sensing products. Modifications detailed in Sect. 2.2.2 were necessary to simulate any emissions from these deposits, let alone resolving accurate surface characteristics. Within this, resolution improvements to vegetation cover input (LAI and FPAR) proved influential for the dust emission scheme given the narrow channel width (∼ 2–3 grid-cells), and the influence of vegetation cover within the shear stress partitioning scheme (e.g. 10 % vegetation cover yields a 65 % reduction to u_*). Due to the binary source factor and riparian vegetation abounding channel deposits, specific care is required in setting spatial alignment and interpolation options for continuous and categorical land-surface datasets. Application of dynamic land-surface characteristics, for example dynamic retrievals of surface albedo (e.g. Seiler et al., 2021), should be considered. Given the propensity of proglacial lakes, variable water cover, or marine-adjacent positioning of deposits, one should ensure the accurate mapping of surface for appropriate treatment of atmospheric dynamics above water surfaces.

4.2.2 Surface characteristics

The principal challenge for proglacial dust modelling resides in attributing sediment textures in active glaciofluvial settings, with deposits exhibiting strong spatial heterogeneity as a function of sediment supply and hydrological processes, alongside temporal variance induced by migrating channels (Maizels, 1993; Marren, 2005). This challenge is not unique to this landscape; large spatial heterogeneity in aeolian activity as a function of landscape unit is a common feature across dust emission regions (Bullard et al., 2011; Zender et al., 2003), where despite widespread potential deflation surfaces, hotspots arguably contribute most critically to total emission loadings (Arnalds et al., 2016; Gillette, 1999; Prospero et al., 2002). The nature of braided channels means that within an area of coarse-grained river terraces (gravel bars), exposed ephemeral channels are primarily silt and much more susceptible to aeolian erosion. At this resolution, the model is unable to resolve these fine channels of high erodibility within a coarser landscape. Depending on river morphology, this issue may persist for larger proglacial sources also and may merit a novel sub-grid scale approach to erodibility in future efforts.

Performing simulations at a resolution where fluvial deposit heterogeneity can be explicitly resolved, a double drag partitioning approach was implemented to include gravel coverage on glaciofluvial deposits. This simple approach helped ensure more representative emissions from coarse deposits in the study region, however extensive mapping of glaciofluvial sediment textures would be required to develop a continuous gravel shear stress partitioning scheme. Remote sensing approaches to estimate surface roughness (Marticorena et al., 2004; Prigent et al., 2005, 2012) have been employed in global models (Darmenova et al., 2009; Klose et al., 2021; Laurent et al., 2006; Leung et al., 2023; Zhang et al., 2022), but remain inapplicable at this resolution of study. Resolving gravel-induced roughness in these valleys remains limited by observations of sufficient resolution to characterise gravel coverage (∼ cm scale) across fluvial deposits. The mechanical strength of surfaces as pertinent for aeolian processes remains largely unconstrained across modelling efforts, despite a variety of observation approaches being developed (e.g. Rice et al., 1997; Li et al., 2010). The selection of frequently reworked surfaces with S_b was implemented to counter uncertainty for surface crusting, fine sediment exhaustion and shortcomings in resolving sparse vegetation with an LAI approach. Despite this, surface crusting was not explicitly resolved in our model, and the erodibility of certain deposits was likely overestimated. The physically-explicit scheme of Shao (2004) offers versatility in simulating emissions from heterogeneous glaciofluvial deposits. However, further application remains restricted in the absence of advancements in mapping sediment texture across dynamic fluvial landscapes, and constraining soil parameter estimates (c_y; P; γ) to explicitly treat soil aggregation state. Separately, we note that the application of schemes of enhanced physical representation for specific sources demands more detail in how model sediment PSDs have been evaluated and parameterized. For glacigenic deposits especially, enhanced fine fractions in soil class vs. model-prescribed soil PSDs may exist (e.g. glaciofluvial sands vs. beach sand).

4.2.3 Soil moisture

Soil moisture discrepancies were anticipated in this study; this approach endeavoured to represent well-drained glaciofluvial deposits from surfaces previously ascribed inaccurate sediment textures or masked as water in land-surface datasets. Results in Sect. 3.1.1 suggested a minimum of 24 h is required for soil moisture adjustment in simulations to realistic values under warm conditions, to be extended in periods of reduced insolation. Soil moisture contents may still fail to reach accurate values regardless of simulation length (Fig. 5), suggesting that LSM drainage characteristics do not permit obtaining accurate levels across short term (∼ weekly) simulations. Beyond surface meteorology, the impact on certain dust-emitting sediments was pronounced (Sect. 3.2.1). Accurately resolving soil moisture conditions is more critical for proglacial dust sources as glaciofluvial deposits can persist as viable emission sources under wet maritime climates (Crusius et al., 2011; Schroth et al., 2017). Sediment drainage and the wind-induced moisture loss critical for deposit erodibility (Crusius, 2021; Schroth et al., 2017) require accurate representation to appropriately treat deposit erodibility and simulate dust emissions. Furthermore, extended investigation of proglacial emissions must contend with time-variant snow cover during transitional seasons, ensuring the land-surface scheme appropriately accounts for snow melting in soil moisture conditions. Pre-simulation modification of initial soil moisture values across erodible surfaces may benefit future “case-study” efforts of simulating proglacial dust emissions.

4.3 Extending to other regions

4.3.1 Approach transferability

Model development and testing was performed with the overarching objective of applying this approach to adjacent proglacial landscapes, becoming a tool to estimate emission loading from these unstudied valleys. This study's application of WRF-Chem has been implemented and tested under arguably the most extreme topographic conditions, for glaciofluvial deposits constrained in narrow channels (0.1–3 km wide), bordered by sheer topography (1.3–1.7 km high). Adjacent dust-producing valleys across the St. Elias Mountains (Alsek, Donjek, Chitina) are all wider, and the Alaskan coastal delta sources are further removed from immediate topography (Crusius et al., 2011; Schroth et al., 2017). While the choice of 3 km HRRR-AK output to initialise and bound the model permitted lower computational costs during this work, the approach is equally implementable with more common reanalysis products such as 0.25° ERA5 (Hersbach et al., 2020) with a series of nested domains; the core of this study's approach lies in the surface modifications required to yield representative emissions from these glaciofluvial deposits. Initialising from coarser datasets may induce sensitivity to lateral boundary conditions (terrain smoothing/discontinuity), requiring additional boundary buffer distances or encompassing the entire region of complex topography so boundaries can be situated over smooth terrain (Warner et al., 1997).

The approach to constrain erodible surfaces to fluvially active deposits (S_b) is highly scalable to assess erodibility across the region. Currently employed with Landsat 5–8 imaging, retrieval timings permit reliable constraint of surface conditions approximately biweekly. The method identifies surfaces subject to regular ephemeral inundation but does not explicitly index long-term deposition features (e.g. proglacial lakes) that may be exposed following major hydrological reorganisation events (Clague and Shugar, 2023). However, currently exposed proglacial lake surfaces were still identified as erodible in this study. The method cannot capture diurnal variability in meltwater that will reduce erodibility diurnally during summer months, however we note that larger glacier systems (such as the Icefield Ranges) have a greater component of baseflow, and at high-latitudes – extended daylight hours – dampening diurnal variability (Marren, 2005).

4.3.2 Discerning missing regional emissions

Contemporary dust emissions have been identified in many Alaskan proglacial valleys by loess studies (Muhs et al., 2013, 2016) and in situ investigations (Barr et al., 2023; Crusius et al., 2011; Schroth et al., 2017), yet these regions continue to be excluded from current global estimates. Explicit calls to resolve these valleys in dust emission modelling efforts (Crusius, 2021; Crusius et al., 2024) have not yet been met. It is unlikely these sources can be explicitly resolved in global-scale approaches due to source spatial extent and the complexity of surface controls, however regional applications such as this work can be leveraged to estimate missing emission loading. While Crusius et al. (2024) suggests the principal obstacle for simulating dust emissions from this region is the lack of high-resolution meteorological data, such data is readily available at 3 km resolution (Dowell et al., 2022) and is likely suitable for simulating emissions from the wide coastal outwash plains. Instead, due to the heterogeneity of depositional facies, the lack of detailed sediment texture (and potentially surface drainage properties) may be key to ensuring representative emissions. The upcoming operational RRFS (Rapid Refresh Forecast System) (NOAA, 2026) is planning to provide hourly forecasts at 3 km resolution across CONUS-AK with a dust emission module based on the FENGSHA scheme (e.g. Zhang et al., 2022). With operational models at a resolution likely sufficient to reproduce meteorological forcing for the larger proglacial valleys and coastal dust sources, dust emissions may remain unresolved due to unrealistic soil texture, deposit roughness, and potentially soil moisture dynamics. Regardless of improvements to model resolution, simulating representative emissions from western Canada and Alaska remains foremost limited by appropriate land-surface treatments. Large potential exists to employ this approach for recently identified glacigenic dust emission from proglacial valleys of Svalbard (Di Mauro et al., 2023; Spolaor et al., 2021; Tobo et al., 2019) and Greenland (Bullard et al., 2023).

Sources across the Canadian Arctic are frequently identified in assessments of emission potential at high latitudes (Meinander et al., 2022) and identified as proximal sources in deposition studies (Zdanowicz et al., 1998). These areas are anticipated to be less operationally complex to simulate due to the less severe topography in contrast to the study site presented in this paper. Conversely, polythermal glacier systems of the Canadian Arctic (Copland et al., 2003) likely yield active fluvial deposits of limited extents (e.g. Pissart et al., 1977; Ranjbar et al., 2021) compared to warm-based glacier systems of coastal Yukon and Alaska, hence model spatial resolution may be challenged to resolve deposits. Due to poor availability of sediment data, these Canadian Arctic “hotspots” require further study before surface erodibility can be appropriately determined. Future application of this modelling approach endeavours to establish an annual dust loading estimate for the St. Elias Mountains, integrating a time-variant source function for a continuous run. Large potential exists to employ this approach for other proglacial dust sources situated within complex terrain (Kok et al., 2021b; Crusius et al., 2024).

5 Conclusions

A dust emission modelling approach covering four proglacial valleys of the St. Elias Mountains (Canada) was developed, resolving mineral dust emissions from this region for the first time. Three simulation periods were investigated in May 2019, June–July 2021 and September 2022, targeting major emission events in the A'ą̈y Chù, Kaskawulsh, Dusty and Alsek River valleys. Model evaluation against in situ meteorological observations, aerosol monitors and Doppler LiDAR retrievals permitted a holistic assessment of model capacity to simulate emission and transport processes within mountainous terrain. It is highly feasible to simulate dust emission timing and magnitude from topographically-confined proglacial valleys (channel width ∼ 1–2 km; sidewalls up to +1.7 km with 40 % slope) with appropriate care given to surface characteristics and model parameterisations. Depending on the complexity of the emission scheme employed, descriptions of sediment texture can be key to simulating emissions. The approach developed in this work offers an option to future research efforts to resolve these emissions, however they will continue to be challenged to better resolve surface heterogeneity across glaciofluvial deposits. With seasonal discrepancies in vertical dispersion in complex terrain, further work is necessary to ensure transport processes can be appropriately treated across CBL conditions to determine the fate of emitted material.

This work provides a first estimate of mineral dust emissions for several proglacial valleys in the region, contributing significantly to a rising need to appraise emission activity from high-latitude sources under global rapid cryosphere loss. Dust emission estimates for high-latitude sources remain exceedingly rare. The approach developed in this work is highly applicable to proglacial dust sources across northwest North America and beyond, currently unresolved by existing modelling efforts. Forthcoming work may extend this approach for annual dust emission estimates from this region and a comprehensive assessment of emission pathways, as pertinent for regional climate impacts.