Modelling hydrologic impacts of light absorbing aerosol deposition on snow at the catchment scale

Light absorbing impurities in snow and ice (LAISI) originating from atmospheric deposition enhance snowmelt by increasing the absorption of shortwave radiation. The consequences are a shortening of the snow duration due to increased snowmelt and, at the catchment scale, a temporal shift in the discharge generation during the spring melt season.

In this study, we present a newly developed snow algorithm for application in hydrological models that allows for an additional class of input variable: the deposition mass flux of various species of light absorbing aerosols. To show the sensitivity of different model parameters, we first use the model as a 1-D point model forced with representative synthetic data and investigate the impact of parameters and variables specific to the algorithm determining the effect of LAISI. We then demonstrate the significance of the radiative forcing by simulating the effect of black carbon (BC) deposited on snow of a remote southern Norwegian catchment over a 6-year period, from September 2006 to August 2012. Our simulations suggest a significant impact of BC in snow on the hydrological cycle. Results show an average increase in discharge of 2.5, 9.9, and 21.4 %, depending on the applied model scenario, over a 2-month period during the spring melt season compared to simulations where radiative forcing from LAISI is not considered. The increase in discharge is followed by a decrease in discharge due to a faster decrease in the catchment's snow-covered fraction and a trend towards earlier melt in the scenarios where radiative forcing from LAISI is applied. Using a reasonable estimate of critical model parameters, the model simulates realistic BC mixing ratios in surface snow with a strong annual cycle, showing increasing surface BC mixing ratios during spring melt as a consequence of melt amplification. However, we further identify large uncertainties in the representation of the surface BC mixing ratio during snowmelt and the subsequent consequences for the snowpack evolution.

1 Introduction

The representation of the seasonal snowpack is of outstanding importance in hydrological models aiming for application in cold or mountainous environments. In many mountain regions, the seasonal snowpack constitutes a major portion of the water budget, contributing up to 50 %, and even more, to the annual discharge (e.g. Junghans et al., 2011). Snowmelt plays a key role in the dynamic of the hydrology of catchments of various high mountain areas such as the Himalayas (Jeelani et al., 2012), the Alps (Junghans et al., 2011), and the Norwegian mountains (Engelhardt et al., 2014), and is an equally important contributor to streamflow generation as rain in these areas. Furthermore, timing and magnitude of the snowmelt are major predictors for flood (Berghuijs et al., 2016) and landslide (Kawagoe et al., 2009) forecasts, and important factors in water resource management and operational hydropower forecasting. Lastly, the extent and the temporal evolution of the snow cover is a controlling factor in the processes determining the growing season of plants (Jonas et al., 2008). For all these reasons, a good representation of the seasonal snowpack in hydrological models is paramount. However, there are large uncertainties in many variables specifying the temporal evolution of the snowpack, and the snow albedo is one of the most important among those due to the direct effect on the energy input to the snowpack from solar radiation (Anderson, 1976). Fresh snow reflects most of the incoming solar radiation in the near-UV and visible spectrum (Warren and Wiscombe, 1980). However, as snow ages and snow grain size increases, the snow albedo will drop as a result of the altered scattering properties of the larger snow grains (Flanner and Zender, 2006). Furthermore, ambient conditions also play a large role. The ratio of diffuse and direct incoming shortwave radiation, the zenith angle of the sun, and the albedo of the underlying ground in combination with the snow thickness can have a large impact on the snow albedo (Warren and Wiscombe, 1980). Of recent significance is the role light absorbing impurities, or particles, which absorb in the range of the solar spectrum, have on albedo when present in the snowpack (e.g. Flanner et al., 2007; Painter et al., 2007; Skiles et al., 2012). These light absorbing impurities in snow and ice (LAISI) can originate from fossil fuel combustion and forest fires in the form of black carbon, BC, and organic carbon (Bond et al., 2013; AMAP, 2015), mineral dust (Painter et al., 2012), volcanic ash (Rhodes et al., 1987), organic compounds in soils (Wang et al., 2013), and biological activity (Lutz et al., 2016), and have species-specific radiative properties.

As LAISI lower the snow albedo, the effect on the snowmelt has the potential to alter the hydrological characteristics of catchments where snowmelt significantly contributes to the water budget. Recent research investigates the impact of LAISI on discharge generation in mountain regions at different scales. Qian et al. (2011) used a global climate model to simulate the effect black carbon and dust in snow have on the hydrological cycle of the Tibetan Plateau. They found a significant impact on the hydrology, with runoff increasing during late winter/early spring and decreasing during late spring/early summer due to a trend to earlier melt dates. Oaida et al. (2015) implemented radiative transfer calculations to determine snow albedo in the Simple Simplified Biosphere (SSiB) land surface model of the Weather Research and Forecasting (WRF) regional climate model. They showed that physically based snow albedo representation can be significantly improved by considering the deposition of light absorbing aerosols on snow. Qian et al. (2009) simulated hydrological impacts due to BC deposition in the western United States using WRF coupled with chemistry (WRF-Chem). They found a decrease in net snow accumulation and spring snowmelt due to BC-in-snow induced increase in surface air temperature.

Only a few studies have developed model approaches to resolve the impact of LAISI on snowmelt discharge generation at the catchment scale. Painter et al. (2010) showed that dust, transported from remote places to the Colorado River basin, can have severe implications for the hydrological regime due to disturbances to the discharge generation from snowmelt during the spring time, shifting the peak runoff by several weeks and leading to earlier snow-free catchments and a decrease in annual runoff. Kaspari et al. (2015) simulated the impact of BC and dust in snow on glacier melt on Mount Olympus, USA, by using measured concentrations in summer horizons and determining the radiative forcing via a radiative transfer model. Results indicate enhanced melt during a year of heavy nearby forest fires, coinciding with an increase in observed discharge from the catchment.

Despite these efforts, the direct integration of deposition mass fluxes of light absorbing aerosols in a catchment model is still lacking. To date, there is no rainfall–runoff model with a focus on runoff forecast at the catchment scale that is able to consider aerosol deposition mass fluxes alongside snowfall. On the other hand, there is evidence that including the radiative forcing of LAISI has the potential to improve the quality of hydrological predictions: Bryant et al. (2013) showed that during the melt period errors in the operational streamflow prediction of the National Weather Service Colorado Basin River Forecast Center are linearly related to dust radiative forcing in snow. They concluded that implementing the effect of LAISI on the snow reflectivity could improve hydrological predictions in regions prone to deposition of light absorbing aerosols on snow, which emphasizes the need for development of a suitable model approach. Furthermore, we continuously move towards hydrological models with an increasingly complex representation of the physical processes involved in the evolution of the seasonal snowpack. Heretofore there has been little focus on the factors related to LAISI, such as the impact of aerosol deposition on snow albedo, that may alter the timing and character of discharge generation at the catchment scale.

In this study, we address this deficiency by introducing a rainfall–runoff model with a newly developed snow algorithm that allows for a new class of model input variable: the deposition mass flux of different species of light absorbing aerosols. The model integrates snowpack dynamics forced by LAISI and allows for analysis at the catchment scale. The algorithm uses a radiative transfer model for snow to account dynamically for the impact of LAISI on the snow albedo and the subsequent impacts on the snowmelt and discharge generation. Aside from enabling the user to optionally apply deposition mass fluxes as model input, the algorithm depends on standard atmospheric input variables (precipitation, temperature, shortwave radiation, wind speed, and relative humidity). To enable a critical evaluation of the newly developed snowpack algorithm, we conduct two independent analyses: (i) a 1-D sensitivity study of critical model parameters, and (ii) a catchment-scale analysis of the impact of LAISI. In both analyses we use BC in snow from wet and dry deposition as a proxy for the impact of LAISI.

We first present an overview of the hydrological model used in this study and our algorithm to treat LAISI in Sect. 2. A description of the catchment used for our study and the input data sets is given in Sect. 3. Section 4 describes the 1-D model experiments and the model settings in the case study. Lastly, our results are presented in Sect. 5 and discussed in Sect. 6.

2 Modelling framework and snowpack algorithm

In the following section we provide descriptions of the hydrologic model (Sect. 2.1) and the formulation of a novel snowpack module used for the analyses (Sect. 2.2).

2.1 Hydrologic model framework

For the analysis, we use Statkraft's hydrologic forecasting toolbox (Shyft; https://github.com/statkraft/shyft), a model framework developed for hydropower forecasting (Burkhart et al., 2016; Ghimirey, 2016; Westergren, 2016). Shyft provides the implementation of many well-known hydrological routines (conceptual parameter models, and more physically based approaches), and allows for distributed hydrological modelling. Standard model input variables are temperature, precipitation, wind speed, relative humidity, and shortwave radiation.

The methods used to simulate hydrological processes are (i) a single-equation implementation to determine the potential evapotranspiration, (ii) a newly developed snowpack algorithm using an online radiative transfer solution for snow to account for the effect of LAISI on the snow albedo, and (iii) a first-order nonlinear differential equation to calculate the catchment response to precipitation, snowmelt, and evapotranspiration. (i) and (iii) are described in more detail herein, while (ii) is described in detail in Sect. 2.2.

To determine the potential evapotranspiration, E_pot, we use the method according to Priestley and Taylor (1972):

$$\begin{array}{}\text{(1)}& E_{\mathrm{pot}}={\displaystyle \frac{a}{\mathit{\lambda }}}\cdot {\displaystyle \frac{s\left(T_{\mathrm{a}}\right)}{s\left(T_{\mathrm{a}}\right)+\mathit{\gamma }}}\cdot K_{\mathrm{n}},\end{array}$$

with a=1.26 being a dimensionless empirical multiplier, γ the psychrometric constant, s(T_a) the slope of the relationship between the saturation vapour pressure and the temperature T_a, λ the latent heat of vaporization, and K_n the net radiation.

The catchment response to precipitation and snowmelt is determined using the approach of Kirchner (2009), who describes catchment discharge from a simple first-order nonlinear differential equation. Following Kirchner (2009), we solve the log-transformed formulation

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\text{d}\left(\ln \right(Q\left)\right)}{\text{d}t}}=g\left(Q\right)({\displaystyle \frac{P-E}{Q}}-\mathrm{1})\end{array}$$

due to numerical instabilities of the original formulation. In Eq. (2), Q is the catchment discharge, E the evapotranspiration, and P the precipitation.

We assume that the sensitivity function, g(Q), has the same form as described in Kirchner (2009):

$$\begin{array}{}\text{(3)}& \ln \left(g\right(Q\left)\right)\approx c_{\mathrm{1}}+c_{\mathrm{2}}\ln \left(Q\right)+c_{\mathrm{3}}\left(\ln \right(Q){)}^{\mathrm{2}},\end{array}$$

with c₁, c₂, and c₃ being the only catchment-specific parameters, which we estimate by standard model calibration of simulated discharge against observed discharge. In contrast to Kirchner (2009)'s approach, we use the liquid water response from the snow routine instead of precipitation P in Eq. (2) (Kirchner, 2009, used snow-free catchments). The response from the snow routine can be liquid precipitation, meltwater, or a combination of both.

2.2 A new snowpack module for LAISI

To account for snow in the model, we developed a snow algorithm to solve the energy balance:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathit{\delta }F}{\mathit{\delta }t}}=K_{\text{in}}(\mathrm{1}-\mathit{\alpha })+L_{\text{in}}+L_{\text{out}}+H_{\mathrm{s}}+H_{\mathrm{l}}+R,\end{array}$$

with the incoming shortwave radiation flux K_in, the incoming and outgoing longwave radiation fluxes L_in and L_out, the sensible and latent heat fluxes H_s and H_l, and the heat contribution from rain R. $\frac{\mathit{\delta }F}{\mathit{\delta }t}$ is the net energy flux into or out of the snowpack. Fluxes are considered to be positive when directed into the snowpack and as such an energy source.

L_in and L_out are calculated using the Stephan–Boltzmann law, with L_in depending on the air temperature T_a and L_out on the snow surface temperature T_ss, calculated as $T_{\mathrm{ss}}=\mathrm{1.16}\cdot T_{\mathrm{a}}-\mathrm{2.09}$ (Hegdahl et al., 2016). The latent and sensible heat fluxes are calculated using a bulk-transfer approach that depends on wind speed, temperature, and relative humidity (Hegdahl et al., 2016).

The main addition provided in the algorithm described herein is the implementation of a radiative transfer solution for the dynamical calculation of snow albedo, α. This implementation allows a new class of model input variable, wet and dry deposition rates of light absorbing aerosols. From this, the model is able to simulate the impact of dust, black carbon, volcanic ash, or other aerosol deposition on snow albedo, snowmelt, and runoff. To account for the mass balance of LAISI, while maintaining a representation of sub-grid snow variability and snow cover fraction (SCF), a tiling approach is applied, where a grid cell's snowfall is apportioned to sub-grid units. Energy balance calculations are then conducted within each tile. Currently, a gamma distribution is used to distribute snowfall to the tiles.

Below, we introduce the radiative transfer calculations required to represent LAISI (Sect. 2.2.1), and provide further details of the sub-grid-scale tiling approach to represent snowpack spatial variability (Sect. 2.2.2).

2.2.1 Aerosols in the snowpack

Wiscombe and Warren (1980) and Warren and Wiscombe (1980) developed a robust and elegant model for snow albedo that remains today as a standard. Critical to their approach was the ability to account for (i) wide variability in ice absorption with wavelength, (ii) the forward scattering of snow grains, and (iii) both diffuse and direct beam radiation at the surface. Furthermore, and of particular importance to the success of the approach, the model relies on observable parameters.

Both the albedo of clean snow and the effect of LAISI on the snow albedo strongly depend on the snow optical grain radius r (Warren and Wiscombe, 1980), which alters as snow ages. r can be related to the snow-specific surface area A_s via

$$\begin{array}{}\text{(5)}& r={\displaystyle \frac{\mathrm{3}}{{\mathit{\rho }}_{\mathrm{ice}}\cdot A_{\mathrm{s}}}},\end{array}$$

with ρ_ice the density of ice. A_s represents the ratio of surface area per unit mass of snow grain (Roy et al., 2013).

In our model, we compute the evolution of A_s in dry snow following Taillandier et al. (2007) as

$$\begin{array}{ll}{\displaystyle }{A}_{\mathrm{s}}\left(t\right)& {\displaystyle }=[\mathrm{0.629}\cdot A_{\mathrm{s},\mathrm{0}}-\mathrm{15.0}\cdot (T_{\mathrm{s}}-\mathrm{11.2}\left)\right]-[\mathrm{0.076}\cdot A_{\mathrm{s},\mathrm{0}}\\ {\displaystyle }& {\displaystyle }-\mathrm{1.76}\cdot (T_{\mathrm{s}}-\mathrm{2.96})]\\ \text{(6)}& {\displaystyle }& {\displaystyle }\ln \left\{t+\exp \left({\displaystyle \frac{-\mathrm{0.371}\cdot A_{\mathrm{s},\mathrm{0}}-\mathrm{15.0}\cdot (T_{\mathrm{s}}-\mathrm{11.2})}{\mathrm{0.076}\cdot A_{\mathrm{s},\mathrm{0}}-\mathrm{1.76}\cdot (T_{\mathrm{s}}-\mathrm{2.96})}}\right)\right\},\end{array}$$

where t is the age of the snow layer (hours), A_s,0 is A_s at t=0 (cm² g⁻¹), and T_s is the snow temperature (^∘C). The evolution of A_s in wet snow is calculated according to Eq. (5) and Brun (1989) as

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\text{d}r}{\text{d}t}}={\displaystyle \frac{C_{\mathrm{1}}+C_{\mathrm{2}}\cdot {\mathrm{\Theta }}^{\mathrm{3}}}{r^{\mathrm{2}}\cdot \mathrm{4}\mathit{\pi }}},\end{array}$$

where $C_{\mathrm{1}}=\mathrm{1.1}\times {\mathrm{10}}^{-\mathrm{3}}$ mm³ d⁻¹ and $C_{\mathrm{2}}=\mathrm{3.7}\times {\mathrm{10}}^{-\mathrm{5}}$ mm³ d⁻¹ are empirical coefficients. Θ is the liquid water content of snow in mass percentage. A_s,0 is set to 73.0 m² kg⁻¹ (Domine et al., 2007) and we set the minimum snowfall required to reset A_s to 5 mm snow water equivalent (SWE).

To solve for the effect of light absorption from LAISI on the snow albedo, we have integrated a two-layer adaption of the Snow, Ice, and Aerosol Radiative (SNICAR) model (Flanner et al., 2007, 2009) into the energy and mass budget calculations. By providing the solar zenith angle of the sun, the snow optical grain radius r, mixing ratios of LAISI in the snow layers, and the SWE of each layer, SNICAR calculates the snow albedo for a number of spectral bands. To achieve this, SNICAR utilizes the theory from Wiscombe and Warren (1980) and the two-stream, multilayer radiative approximation of Toon et al. (1989). Following Flanner et al. (2007), our implementation of SNICAR uses five spectral bands (0.3–0.7, 0.7–1.0, 1.0–1.2, 1.2–1.5, and 1.5–5.0 µm) in order to maintain computational efficiency. Flanner et al. (2007) compared results from the 5-band scheme to the default 470-band scheme in SNICAR and concluded that relative errors are less than 0.5 %. The incident fluxes were simulated offline assuming mid-latitude winter clear- and cloudy-sky conditions.

The absorbing effect of LAISI is most efficient when the LAISI reside at or close to the snow surface (Warren and Wiscombe, 1980). As snow melts, LAISI can remain near the surface due to inefficient melt scavenging, which leads to an increase in the near-surface concentration of LAISI and thus a further decrease in the snow albedo – the so-called melt amplification (e.g. Xu et al., 2012; Doherty et al., 2013, 2016; Sterle et al., 2013). Field observations suggest that the magnitude of this effect is determined by the particle size and the hydrophobicity of the respective LAISI (Doherty et al., 2013). Conway et al. (1996) observed vertical redistribution and the effect on the snow albedo by adding volcanic ash and hydrophilic and hydrophobic BC to the snow surface of a natural snowpack. Flanner et al. (2007) used the results from Conway et al. (1996) to determine the scavenging ratios, specifying the ratio of LAISI contained in the melting snow that is flushed out with meltwater, of both hydrophilic and hydrophobic BC. They found the scavenging ratio for hydrophobic BC, k_phob, to be 0.03, and, for hydrophilic BC, k_phil, 0.2. Doherty et al. (2013) found similar results by observing BC mixing ratios close to the surface of melting snow. However, more recent studies report efficient removal of BC with meltwater (Lazarcik et al., 2017), revealing large gaps in the understanding of the process.

Figure 1 (a) Elevation versus coefficients of variation (CVs) of sub-grid snow distribution from Gisnås et al. (2016) of forest-free areas in the Atnsjoen catchment (dots) and the relationship between the CVs and the elevation resulting from simple linear regression analysis (black line). (b) Solid precipitation multiplication factors for the sub-grid snow tiles for different CVs.

To represent the evolution of LAISI mixing ratios near the snow surface, we treat LAISI in two layers in our model. The surface layer has a time-invariant maximum thickness (further called the maximum surface layer thickness). The mixing ratio of each LAISI species in this layer is calculated from a uniform mixing of the layer's snow with either falling snow with a certain mixing ratio of aerosol (wet deposition) or aerosol from atmospheric dry deposition. The second layer (bottom layer) represents the snow exceeding the maximum thickness of the surface layer. Following Krinner et al. (2006), we apply a maximum surface layer thickness of 8 mm SWE. Krinner et al. (2006) suggest this value is based on observations of 1 cm thick dirty layers in alpine firn cores used to identify summer horizons. Due to potential accumulation of LAISI in surface snow via dry deposition and melt amplification, we expect the simulated surface mixing ratios of LAISI to be sensitive to the maximum surface layer thickness of our model. For this reason, we use a factor of 2 for the maximal surface layer thickness to account for the uncertainty of this model parameter.

To allow for melt amplification in the model, we include LAISI mass fluxes between the two layers during snow accumulation and snowmelt. Generalizing Jacobson (2004)'s representation of LAISI mass loss due to meltwater scavenging for multiple snow layers, we characterize the magnitude of melt scavenging using the scavenging ratio k and calculate the temporal change in LAISI mass m_s in the surface layer as

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\text{d}{m}_{\text{s}}}{\text{d}t}}=-k{q}_{\mathrm{s}}{c}_{\mathrm{s}}+D,\end{array}$$

and the change in LAISI mass m_b in the bottom layer as

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\text{d}{m}_{\text{b}}}{\text{d}t}}=k(q_{\mathrm{s}}{c}_{\mathrm{s}}-q_{\mathrm{b}}{c}_{\mathrm{b}}).\end{array}$$

Herein, q_s and q_b are the mass fluxes of meltwater from the surface to the bottom layer and out of the bottom layer, respectively, and c_s and c_b are the mass mixing ratios of LAISI in the respective layer. D is the atmospheric deposition mass flux. A value for k of < 1 is equal to a scavenging efficiency of less than 100 % and hence allows for accumulation of LAISI in the surface layer during melt. In our analysis, we account for hydrophobic and hydrophilic BC. By following Flanner et al. (2007), we set k_phob to 0.03 and k_phil to 0.2, and account for the large uncertainty in those estimates by using an order of magnitude variation on k_phob and k_phil. Like Flanner et al. (2007), we treat aged, hydrophilic BC as sulfate coated to account for the net increase in the mass absorption cross section (MAC) by 1.5 at λ=550 nm compared to hydrophobic BC caused by the ageing of BC (reducing effect on MAC) and particle coating from condensation of weakly absorbing compounds (enhancing effect on MAC) suggested by Bond et al. (2006). As a consequence, hydrophilic BC absorbs more strongly than hydrophobic BC under the same conditions. On the other hand, hydrophilic BC undergoes a more efficient melt scavenging. The competing mechanisms are subjects of the 1-D sensitivity study in Sect. 5.1.3.

2.2.2 Sub-grid variability in snow depth and snow cover

In order to allow for explicit treatment of snow layers while representing sub-grid snow variability, we follow Aas et al. (2017), and assume that the sub-grid spatial distribution of each single event of solid precipitation follows a certain probability distribution function. From this distribution we calculate multiplication factors, which then are used to assign the snowfall of a model grid cell to a number of sub-grid computational elements, the so-called tiles (Aas et al., 2017). The snow algorithm described herein is executed for each of the tiles separately, providing a mechanism to account for snow spatial distribution while preserving conservation of mass. Therefore, variables related to the snow state, such as SWE, liquid water content, LAISI mixing ratios, and snow albedo differ among the tiles. To calculate the multiplication factors, we assume that the sub-grid redistributed snow follows a gamma distribution (see e.g. Kolberg and Gottschalk, 2010; Gisnås et al., 2016), determined by the coefficient of variation (CV) of SWE at snow maximum. Gisnås et al. (2016) used Winstral and Marks (2002)'s terrain-based parametrization to model snow redistribution in Norway by accounting for wind effects during the snow accumulation period over a digital elevation model with 10 m resolution. In the case study presented in Sect. 5.2, we use the CV values from Gisnås et al. (2016) to derive a linear relationship between the model grid cell's elevation and the corresponding CV value by simple linear regression (see Fig. 1a), which results in a R² value of 0.71 and a p value of smaller than 2.0 × 10⁻⁵ for the study area. The linear relationship is only applied to grid cells with an areal forest cover fraction of lower than or equal to 0.5. For grid cells with a forest cover fraction of higher than 0.5, a constant snow CV value of 0.17 is used, following the findings of Liston (2004) for high-latitude, mountainous forest. Examples of multiplication factors for forested grid cells and forest-free grid cells for different CV values are shown in Fig. 1b.

3 Site description, meteorologic model input, and atmospheric deposition data

We selected the unregulated upper Atna catchment for our analysis. The catchment is located in a high-elevation region of southern Norway (Fig. 2). The watershed covers an area of 463 km² and ranges in elevation from 700 m a.s.l. at the outlet at lake Atnsjoen to over 2000 m a.s.l. in the Rondane mountains in the western part of the watershed, with approximately 90 % of the area above the forest limit. The average annual precipitation in the watershed during the study period is approximately 655 mm. The mean annual discharge is approximately 11 m³ s⁻¹, with low flows of 1–3 m³ s⁻¹ during the winter months and peak flows of over 130 m³ s⁻¹ during the spring melt season.

For the meteorological model input of precipitation, temperature, relative humidity, and wind speed we use daily observations from the Norwegian Water Resources and Energy Directorate (NVE) and the Norwegian Meteorological Institute (MET). Four meteorological stations are located in the watershed at elevations between 701 and 780 m a.s.l. along the Atna river, two of these measuring precipitation and two measuring temperature. Observations of relative humidity and wind speed originate from two stations at locations close by the catchment (not shown in Fig. 2). Further information about the stations are given in Table 1. Due to poor availability of continuous solar radiation observations in Norway, we use gridded global radiation data from the Water and Global Change (WATCH) Forcing Data methodology applied to ERA-Interim reanalysis data (WFDEI; Weedon et al., 2014) with a resolution of 0.5^∘. Discharge observations are from a station located at the outlet of the catchment at lake Atnsjoen and are used for model calibration and validation. In the following, we present the development of atmospheric deposition rates of BC, which we use as a proxy for LAISI, due to a lack of available deposition rates for other species. For the 1-D sensitivity study of Sect. 5.1 we developed representative model input based on the meteorological conditions in this catchment.

Figure 2 Location of the Atnsjoen catchment in Norway (black box in the left map) and overview map of the Atnsjoen catchment (right).

Table 1 Information about observational stations.

Atmospheric deposition of black carbon from the REMO-HAM model

The wet and dry deposition rates of BC for the study area are generated using the REMO-HAM regional aerosol-climate model (Pietikäinen et al., 2012). The core of the model is a hydrostatic, 3-D atmosphere model developed at the Max Planck Institute for Meteorology in Hamburg. With the aerosol configuration, the model incorporates the HAM (Hamburg Aerosol Module) by Stier et al. (2005) and Zhang et al. (2012). HAM calculates the aerosol distributions using seven log-normal modes and includes all the main aerosol processes.

For the simulations, we follow the approach of Hienola et al. (2013), but with changes to the emission inventory: Hienola et al. (2013) used emissions based on the AeroCom emission inventory for the year 2000 (see Dentener et al., 2006). In the REMO-HAM simulations conducted herein, emissions are made by the International Institute for Applied Systems Analysis (IIASA) and are based on the Evaluating the Climate and Air Quality Impacts of Short-Lived Pollutants (ECLIPSE) V5a inventory for the years 2005, 2010, and 2015 (years in between were linearly interpolated) (Klimont et al., 2016, 2017). We also updated other emissions modules (wildfire, aviation, and shipping) following the approaches presented in Pietikäinen et al. (2015). The only difference to Pietikäinen et al. (2015) in this work is that we used the Global Fire Emissions Database (GFED) version 4 based on an updated version of van der Werf et al. (2010).

REMO-HAM was used for the same European domain as in Pietikäinen et al. (2012) using a 0.44^∘ spatial resolution (50 km), 27 vertical levels, and a 3 min time step. The ERA-Interim re-analysis data were utilized at the lateral boundaries for meteorological forcing (Dee et al., 2011) and, for the lateral aerosol forcing, data from the ECHAM-HAMMOZ global aerosol-climate model (version echam6.1.0-ham2.2) were used. ECHAM-HAMMOZ was simulated in a nudging mode, i.e. the model's meteorology was forced to follow ERA-Interim data, and the ECLIPSE emissions were used (plus the other updated emission modules shown in Pietikäinen et al., 2015). The boundaries of REMO-HAM were updated every 6 h for both meteorological and aerosol related variables. Simulations with REMO-HAM were conducted for the time period of 1 July 2004–31 August 2012 and the time period used in our analysis is from 1 September 2006 onwards. The initial state for the model was taken from the boundary data, except for the soil parameters, which were taken from a previous long-term simulation for the same domain (a so-called warm start). The output frequency of REMO-HAM was 3 h and the total BC deposition flux was calculated from the accumulated dry and wet deposition and sedimentation fluxes, and resampled to daily time resolution. Herein, dry deposition refers to the sum of REMO-HAM dry deposition and sedimentation.

Table 2 Model parameters used in the sensitivity and case study. Parameters optimized during calibration are marked with ^a. Further parameters were pre-set and not included in parameter estimation during calibration. Parameters with different values in the minimum (min), central (mid), and maximum (max) BC radiative forcing estimates are marked with ^b.

4 Modelling experiments and calibration

Our analysis is conducted in two parts. First, in a 1-D sensitivity study, we investigate the sensitivity of parameters and variables specific to the LAISI algorithm presented in Sect. 2.2. We then demonstrate the impact of BC at the catchment scale in a case study by simulating the impact of wet and dry deposition of BC on snowmelt and discharge generation in a remote southern Norwegian catchment (Sect. 5.2).

We assume uncertainties of the LAISI radiative forcing in snow to originate mainly from the model representation of surface layer thickness, melt scavenging of BC, and uncertainties in the deposition input data. To account for the uncertainties, we declare minimum (min), central (mid), and maximum (max) effect estimates for each of the critical parameters, outlined together with further model parameters in Table 2. The min, mid, and max estimates are both subjects of analysis in the sensitivity study (further described in Sect. 4.1) and used in the case study to give an uncertainty estimate of the LAISI effect on the hydrologic variables (further described in Sect. 4.2). We investigate the impact of BC impurities on the response variables by comparing the results from aerosol radiative forcing model experiments (ARF scenarios) to simulations in which all BC deposition rates are set to zero (no-ARF scenario).

4.1 One-dimensional sensitivity study experiments

The results of the 1-D sensitivity study are presented in Sect. 5.1; herein we describe the configurations to conduct our analysis. The purpose of this study is to isolate the impact of different model parameters: (i) maximum surface layer thickness (parameter max_surface_layer; see Table 2), (ii) scavenging ratio, and (iii) the impact of the scavenging ratio with respect to the BC species (parameters k_phob and k_phil).

Our approach evaluates these parameters and the evolution of the snowpack under constant melting conditions. We run the 1-D simulations with model parameters as outlined in Table 2 and forcing data based on synthetic input data. The synthetic forcing data set is based on the average meteorological conditions during the melt season from mid March until mid July of the Atnsjoen catchment. In our sensitivity experiments, all snowpacks have 250 mm SWE of snow with a mixing ratio of 35 ng g⁻¹ in both surface and bottom layer at melt onset. These values are representative of the upper 50 % of tiles at winter snow maximum in the Atnsjoen catchment during the study period of the case study. During the melt period, we exclude fresh snowfall and dry deposition, in order to isolate the effect of the tested model parameters on the snowpack evolution under melt conditions. This may lead to an underestimation of total BC mass in the snow column.

To investigate the impact of the maximum surface layer thickness (parameter max_surface_layer) of the model, we run simulations with synthetic forcing and use maximal surface layer thicknesses of 4.0 mm SWE (max estimate; see Table 2), 8.0 mm SWE (mid estimate), and 16.0 mm SWE (min estimate). Additionally, we include a single-layer model with a vertically uniform distribution of BC in the analysis and, for comparison, a simulation with clean snow.

To explore the sensitivity to scavenging ratio, we apply different BC scavenging ratios in the range of the uncertainty of hydrophilic BC, which covers a wide range from very efficient to inefficient scavenging. The scavenging ratios applied are based on the analysis conducted by Flanner et al. (2007) using data from Conway et al. (1996). The mid estimate for the hydrophilic BC scavenging ratio (k_phil=0.2) also compares well to field observations from Doherty et al. (2013). We further include in the analysis Flanner et al. (2007)'s upper bound uncertainty estimate for hydrophilic BC (2.0; efficient scavenging), the lower bound estimate (0.02; inefficient scavenging), and for comparison a scenario in which BC does not undergo any scavenging (0.0).

Hydrophilic BC absorbs more strongly than hydrophobic BC under the same conditions due to an increased MAC for hydrophilic BC resulting from ageing of the aerosol during atmospheric transport (Bond et al., 2006). On the other hand, hydrophilic BC undergoes more efficient melt scavenging (Flanner et al., 2007), which impacts the snowpack evolution significantly. To explore this competing interplay, we apply the mid estimate of the scavenging ratio of hydrophobic BC (k_phob=0.03) to both the hydrophobic BC and the hydrophilic BC species. In this manner we explore the isolated effect of the different absorption properties of the two species. We further apply the mid estimate for hydrophilic BC scavenging ratio (k_phil=0.2) to hydrophilic BC to quantify the gross effect. As in other cases, we include the no-ARF scenario to highlight the overall effect on the albedo and melt of the different scenarios.

4.2 Case study model set-up and calibration

We investigate the impact of BC aerosol deposition on the catchment hydrology of a Norwegian catchment over a study period of 6 years, from September 2006 to August 2012. The station-based input data described above (Sect. 3) are interpolated to the simulation grid cells (1 × 1 km² and accordingly smaller cells at the catchment boarders; see Fig. 2) using Shyft's interpolation algorithms. For temperature Bayesian Kriging (Diggle and Ribeiro, 2007) is used. For precipitation, BC deposition rates, wind speed, and relative humidity interpolation to the model grid cells are via inverse distance weighting. A 5 % increase in precipitation for every 100 m increase in altitude is used for the precipitation interpolation (Førland, 1979).

To calibrate the model against observed discharge, we first run a split-sample calibration (Klemes, 1986) using the first 3 years (1 September 2006 to 31 August 2009) of the study period as the calibration period and the following 3 years (1 September 2009 to 31 August 2012) for model validation. For parameter estimation, we use the BOBYQA algorithm for bound-constrained optimization (Powell, 2009). To assess the predictive efficiency of the model, we use the Nash–Sutcliffe model efficiency

$$\begin{array}{}\text{(10)}& E_{\text{NS}}=\mathrm{1}-{\displaystyle \frac{{\sum }_{t=\mathrm{0}}^T(Q_{\mathrm{o}}^t-Q_{\mathrm{s}}^t{)}^{\mathrm{2}}}{{\sum }_{t=\mathrm{0}}^T(Q_{\mathrm{o}}^t-\stackrel{\mathrm{‾}}{Q_{\mathrm{o}}}{)}^{\mathrm{2}}}},\end{array}$$

where $Q_{\mathrm{o}}^t$ and $Q_{\mathrm{s}}^t$ are the observed and simulated discharge at time t, respectively, and $\stackrel{\mathrm{‾}}{Q_{\mathrm{o}}}$ is the mean observed discharge over the assessed period.

Model calibration is run with mid estimates for all model parameters impacting the handling and effect of LAISI and aerosol depositions as simulated from REMO-HAM during model calibration. Those parameters and further model parameters, including the parameters estimated during calibration, are listed in the left column of Table 2. We investigate the uncertainty in the effect of BC on snowmelt by using the min and max effect parameter estimates from Table 2, while holding constant all other model parameters as estimated during calibration. To assess the gross effect of LAISI, we compare the simulations to equivalent simulations in which ARF is not included.

Figure 3 Snow albedo (top row of graphs; solid lines) and melt rate (top row of graphs; dashed lines), BC mixing ratio in the surface layer and factor increase in the mixing ratio during melt compared to the pre-melt BC mixing ratio (central row of graphs), and snowpack SWE (bottom row of graphs) for simulations forced with synthetic data based on average meteorological conditions during the melt season from mid March until mid July of the Atnsjoen catchment and different model configurations: (a) different values for maximum surface layer thickness; (b) scavenging ratio; and (c) BC species with different melt scavenging ratios applied (phob and phil in the legend stand for hydrophobic and hydrophilic BC, respectively). The black lines in all graphs show simulation results of model runs without ARF applied (no-ARF).

5 Results

5.1 One-dimensional sensitivity studies

5.1.1 Sensitivity to surface layer thickness

Figure 3a shows the effect of the different maximum surface layer thicknesses (parameter max_surface_layer) on the melting snowpack with other parameters set according to Table 2. The maximum surface layer thickness strongly determines the surface BC mixing ratio over the melt season. During snowmelt, surface BC increases up to a factor of circa 10, 20, and 30 for maximum surface layer thicknesses of 16.0, 8.0, and 4.0 mm SWE, compared to the pre-melt season BC mixing ratio (35 ng g⁻¹). For those three two-layer scenarios (purple, red, and green curves in Fig. 3a), the resulting differences in albedo and melt rate are small, even though the increase in the surface layer mixing ratio during the melt season differs strongly among the scenarios. Using the single layer model, the surface BC mixing ratio increases more slowly and stays comparably low in contrast to the two-layer models until shortly before meltout. This leads to a less pronounced decrease in albedo compared to the two-layer models and thus to a shorter meltout shift compared to a clean snowpack of about 5 days (yellow curves in Fig. 3a), whereas the two-layer scenarios show earlier meltouts of about 7 days.

5.1.2 Sensitivity to scavenging ratio of BC

In the range of investigated scavenging ratios, we find the sensitivity of the surface BC mixing ratio, the albedo, and the subsequent snowmelt to this parameter (Fig. 3b). When applying a melt scavenging factor typical for the lower bound of hydrophilic BC (0.02, purple lines) there is little effect compared to the scenario without melt scavenging (green lines). Both show circa a factor 30 increase in surface BC mixing ratio to the end of the melt season and only little differences in the development of albedo and snowmelt. Similar results are achieved when using the mid estimate scavenging factor for hydrophobic BC (0.03, not shown). A distinction exists when using the mid estimate scavenging factor for hydrophilic BC (0.2, red line). In contrast to no scavenging and the lower bound hydrophilic scavenging, surface BC does not increase as rapidly during the melt period and is completely flushed when applying a melt scavenging factor typical for the upper bound of hydrophilic BC (yellow line; the surface concentration drops continuously during the melt period).

The changes in the scavenging ratio lead to a considerable effect on albedo and snowmelt. Meltout is delayed by circa 0.5 (purple lines), 3 (red lines), and 8 days (yellow lines) for scavenging ratios of 0.02, 0.2, and 2.0, respectively, compared to no scavenging (green lines). Compared to the no-ARF experiment (black lines), our simulations show that the presence of BC causes an earlier meltout of about 9.5, 7, and 2 days for scavenging ratios of 0.02, 0.2, and 2.0, respectively.

5.1.3 Sensitivity to BC species

The column of graphs in Fig. 3c illustrate the net effect of the competing processes of more efficient absorption resulting from a larger MAC with more efficient wash out. A mid estimate of the scavenging ratio of hydrophobic BC (0.03) is applied and shown for the hydrophobic BC (green curve) and the hydrophilic BC (purple curves) species. These curves show the isolated effect of the different absorption properties of the two species. Further, the mid estimate scavenging ratio for hydrophilic BC (0.2) is also shown using radiative properties of hydrophilic BC to quantify the gross effect (red curves). The no-ARF scenario (black curves) highlights the overall impacts.

The isolated effect of the stronger absorption of hydrophilic BC leads to an earlier meltout by circa 2 days compared to hydrophobic BC (purple and green curves in graphs of Fig. 3c). However, when applying the mid estimate of the scavenging ratio for hydrophilic BC (0.2), the combined effects leads to a masking of the isolated effect of stronger absorption by hydrophilic BC (and vice versa). During the melt period, snow albedo, melt rate and the snowpack SWE barely differ between the scenarios with the mid estimate scavenging for hydrophobic and hydrophilic BC applied (green and red curves). This reveals that both scenarios, hydrophobic BC with low scavenging efficiency and hydrophilic BC with high scavenging efficiency, lead to an earlier meltout by roughly 7 days compared to the no-ARF scenario.

5.2 Case study: Impact of BC deposition on the hydrology of a south Norwegian catchment

5.2.1 Performance of the model

In the split-sample test, the model performance is acceptable during both calibration and validation, with Nash–Sutcliffe model efficiencies of 0.86 during the calibration period (green line in Fig. 4a) and 0.82 during the validation period (red line in Fig. 4a). However, in the winter season (November until March) the model generally underestimates the discharge and peaks in the beginning of the melt season are slightly underestimated. The scatter plot in Fig. 5 confirms the underestimation of low-flow situations. For the different scenarios explored within the case study, all LAISI-relevant parameters are fixed to mid estimates and model parameters optimized for the full period (1 September 2006 to 31 August 2012; Fig. 4b) resulting in a Nash–Sutcliffe model efficiency of 0.84. The optimized parameters are listed in Table 2. Note that switching ARF off entirely (no BC deposition) leads to a slight decrease in the model quality (Nash–Sutcliffe model efficiency of 0.83 over the whole period; not shown).

Figure 4 Simulated (green and red curves) and observed (black curve) daily discharge from the Atnsjoen catchment. (a) is showing the simulation results for 3 years of calibration (green) and 3 years of validation (red). (b) is showing the results for the 6-year calibration period. Parameters estimated in the latter are used in the case study. Parameters not included in the optimization are set to mid estimate values during the calibration process (see Table 2).

Figure 5 Comparison of observed and simulated daily discharge Q of the Atnsjoen catchment. The dashed black line demonstrates perfect agreement between simulation and observation.

5.2.2 Surface BC and albedo

For the min and mid estimate, the model simulates an average annual surface BC mixing ratio of about 18 and 71 ng g⁻¹, respectively. Our max estimate yields 198 ng g⁻¹. The evolution of surface albedo driven by BC deposition is distinct in the accumulation period vs. the melt period. During the snow accumulation period (until end of March), only slight differences in albedo are noticeable. The average annual snow albedo from 1 January until 22 March is 0.871 for the no-ARF scenario (Fig. 6a), while during the same time period, min, mid, and max estimates show relative albedo reductions of 0.003, 0.010, and 0.014, respectively from the no-ARF case. At the beginning of the melt period, surface layer concentrations of min, mid, and max estimate average to 12, 49, and 98 ng g⁻¹ (Fig. 6b).

Figure 6 (a) Simulated mean catchment snow albedo (solid lines) and snow-covered fraction (SCF; dashed lines) for the mid (red lines), min, and max (shaded) estimates and for the scenario without ARF (no-ARF; black lines) averaged over the 6-year period. (b) Mixing ratio of BC in the model surface layer for the mid (solid line), min (lower bound of the shaded area), and max (upper bound of the shaded area) estimates.

With the start of the melt season, the difference in albedo between model experiments becomes increasingly larger over time. During the melt season, the mid estimate spatially averaged surface BC mixing ratio increases from 49 to about 250 ng g⁻¹ (factor 5 increase) at the end of the melt season (beginning of July). For the max estimate, the increase is from roughly 100 to over 2500 ng g⁻¹ (factor 25 increase). The min estimate on the other hand leads to a decrease in the BC surface mixing ratio. The distinctly different surface BC mixing ratios at the end of the melt season and among the three scenarios cause large differences in albedo decrease relative to the no-ARF case of about 0.03, 0.1, and over 0.3 for the min, mid, and max estimates, respectively.

5.2.3 BC-induced radiative forcing

The radiative forcing in snow (RFS) induced by the presence of BC is calculated from the average radiative forcing over snow-bearing tiles only. The RFS represents the additional uptake of energy from solar radiation per area snow cover due to the presence of BC in the snow compared to clean snow with the same properties. Figure 7a shows the daily mean RFS and demonstrates the increase in RFS during snowmelt. Low RFS is observed during the snow accumulation period, then steadily increasing through spring snowmelt, reaching values of approximately 8, 18, and 57 W m⁻² for the min, mid, and max estimates, respectively (see the red solid line and shaded area in Fig. 7a). RFS in mid winter is small due to low surface BC mixing ratios and low solar irradiance.

Figure 7 Catchment snow-covered fraction (SCF; dashed lines), (a) simulated mean radiative forcing in snow (RFS), and (b) total daily energy uptake in the catchment due to BC (surface radiative forcing in Watts per square metre catchment area) for the mid (solid red lines), min (lower bound of the shaded area), and max (upper bound of the shaded area) estimates averaged over the 6-year period).

However, most relevant for discharge generation (see Sect. 5.2.4) is the catchment-wide total daily energy uptake due to BC, or surface radiative forcing, calculated as the mean radiative forcing over all grid cells. As the snow cover fraction (SCF) in the catchment drops during spring (dotted line and yellow shaded area in Figs. 6 and 7), the effect of the RFS on the melt generation is limited by the increasing area of bare ground. The net effect is shown in Fig. 7b. The catchment mean surface radiative forcing due to the presence of BC in snow shows a strong annual cycle and reaches a maximum of 1.3, 4.9, and 8.8 W m⁻² (min, mid, and max estimates, respectively) around the beginning of May.

5.2.4 BC impact on catchment discharge and snow storage

Figure 8a shows the simulated daily discharge and catchment SWE averaged over the 6-year simulation period for the mid (red lines), min, and max estimates (bounds of the shaded areas), and the no-ARF scenario (black lines). The differences in daily discharge and catchment SWE of the min, mid, and max estimates to the no-ARF scenario are shown in Fig. 8b. All simulations with ARF applied show higher daily discharge from the end of March until the end of May and lower discharge from the end of May until mid August relative to the no-ARF simulation. For the rest of the year, no effect on discharge is noticeable. The net impact of RFS results in a shift in the timing of discharge. Higher discharge early in the melt season is observed, yet offset by lower discharge following May. The cumulative annual discharge remains nearly identical.

Figure 8 (a) Simulated daily discharge (Q; solid lines) and catchment mean snow water equivalent (SWE; dashed lines) for the mid (red lines), min, and max (shaded) estimates and for the scenario without ARF (no-ARF; black lines) averaged over the 6-year period. (b) Differences in daily discharge and SWE between ARF scenarios and the scenario without ARF (no-ARF). The blue marker in (a) and (b) separates the periods where BC in snow has an enhancing effect (left of the marker) and a decreasing (right of the marker) effect on discharge.

Table 3 Average change in discharge during the early (22 March to 29 May) and late (30 May to 10 August) melt seasons of min, mid, and max estimates and average change in SWE during the melt season (22 March to 10 August) compared to the no-ARF scenario.

Min, mid, and max estimates all show the change from higher to lower discharge compared to the no-ARF scenario approximately at the same time (at the end of May; see the blue marker in Fig. 8). Therefore, we can quantify the absolute and relative effect of RFS on the discharge during the two periods: the early melt season from circa 22 March until 29 May and the late melt season from circa 30 May until 10 August (Fig. 8b, and see Table 3). This yields an average percentage increase in daily discharge of 2.5, 9.9, and 21.4 % for the min, mid, and max estimates for the early melt season and a decrease in discharge of −0.8, −3.1, and −6.7 % during the late melt season.

The differences in discharge among the scenarios can be explained by understanding the evolution of the snowpack. In all scenarios the catchment SWE (Fig. 8a) reaches a peak reduction relative to the no-ARF scenario of −4.6, −13.4, and −34.4 % in mid May. The average decreases in catchment SWE of the min, mid, and max estimates compared to the no-ARF scenario during the entire melt season are −2.1, −7.4, and −15.1 % (see Table 3). From mid May on, the differences in catchment SWE between scenarios drop continuously, which is equivalent to a higher catchment averaged snowmelt rate in the no-ARF scenario compared to the ARF scenarios.

6 Discussion

The objective of this work is to provide a mechanism to assess the impact of light absorbing aerosols on runoff at the catchment scale in a rainfall–runoff modelling context. Prior investigations into LAISI indicate potentially significant impacts on the cryosphere (Flanner et al., 2007) with potential impacts on water resources (Qian et al., 2009, 2011). Earlier studies on hydrologic impacts at the catchment scale have used altered radiative forcings to evaluate the impact on the timing of snowmelt and hydrology (Painter et al., 2010; Skiles et al., 2012). With the approach presented herein, we seek to fill a gap between land-surface model approaches (e.g. Oaida et al., 2015) and approaches that apply modified radiative forcing to provide a novel tool for hydrologic forecasting.

6.1 Parameter sensitivity

To assess the sensitivity of the newly introduced algorithm and parameters, we conducted a sequence of 1-D sensitivity studies. In this context, we are able to remove complexities that arise when conducting distributed simulations at the catchment scale.

We found the greatest sensitivity to lie in the parametrization of scavenging, as it relates to how likely the aerosol is to remain at the snow surface during melt. Field measurements indicate that only a fraction of BC is flushed out with the meltwater and BC can accumulate near the snow surface (e.g. Xu et al., 2012; Doherty et al., 2013, 2016; Sterle et al., 2013). Our model is able to simulate this process by taking the scavenging ratio of BC during meltwater movement into account (Eqs. 7 and 8). In the literature, the scavenging efficiency of BC is discussed controversially. Flanner et al. (2007)'s estimates for scavenging ratios of hydrophilic and hydrophobic BC, which are used in this study, are based on data from field experiments using artificially added soot (Conway et al., 1996). However, parameters derived from artificially added soot might not be directly transferable to the scavenging properties of naturally occurring BC. Even though field observations from Doherty et al. (2013) agree well with the estimates of Flanner et al. (2007), and further studies highlight the importance of BC retention in the snowpack (e.g. Xu et al., 2012; Sterle et al., 2013), a large uncertainty remains on the magnitude of this effect (Lazarcik et al., 2017). These uncertainties are identified in our simulations as results show large differences in BC evolution and day of meltout at the boundaries of the applied scavenging ratios (Fig. 3b). Compared to the no-ARF experiment (black lines), the presence of BC causes an earlier meltout for all scavenging ratios applied, spanning from 2 days (upper boundary hydrophilic scavenging, 2.0) to about 9.5 days (lower boundary hydrophilic scavenging, 0.02). Even when applying efficient melt scavenging, resulting in nearly all BC removed from the snow, the meltout still happens circa 2 days earlier compared to the no-ARF experiment.

Further complicating the effect is the fact that hydrophilic BC (which undergoes more efficient melt scavenging) has a larger MAC (enhanced absorption) compared to hydrophobic BC (Flanner et al., 2007). Our results suggest that distinguishing between species may play a secondary role in the determination of the overall impact of BC on snowmelt due to the compensating effect of stronger scavenging accompanied by stronger absorption and vice versa (Fig. 3c).

The 1-D model experiments further show that the definition of at least two layers in the snowpack model is important to allow for accumulation of impurities at the snow surface. This result in itself is not original: numerous prior studies have identified the importance of having multiple layers (Krinner et al., 2006; Flanner et al., 2007; Oaida et al., 2015). However, we further find that the model surface layer thickness (parameter max_surface_layer; see Table 2) has a great impact on the evolution of surface mixing ratios of BC, while at the same time the effect on albedo and snowmelt is small. This results from the fact that for all two-layer models the surface layer thickness is much thinner than the penetration depth of shortwave radiation. For example, in clean snow with an optical grain radius of 50 µm, the radiative intensity diminishes to 1∕e of its surface value (the so-called penetration depth) in 25.5 mm SWE. For snow with an optical grain radius of 1000 µm, the penetration depth increases to 117 mm SWE (both results from Flanner et al., 2007, assuming a wavelength of 550 nm and a solar zenith angle of 60^∘). Thus, BC in the surface layer absorb efficiently in all two-layer scenarios and the difference in the albedo is relatively large compared to the no-ARF scenario (solid black line in the top graph of Fig. 3a), but relatively small among the two-layer scenarios (solid purple, red, and green curves in the top graph of Fig. 3a). However, there is a critical difference when a single-layer model is used (yellow curves in Fig. 3a) due to the aerosol being distributed uniformly throughout the snowpack instead of allowing accumulation at the surface. Thus, a large fraction of the BC is located at depths where the radiative intensity is much lower than in the top few mm of the snowpack. This leads to a weaker absorption efficiency and a less pronounced decrease in albedo compared to the two-layer models and thus to a shorter meltout shift compared to a clean snowpack than in the two-layer scenarios.

Observations of BC in melting snow support the accumulation of BC near the surface (Xu et al., 2012; Doherty et al., 2013; Sterle et al., 2013; Delaney et al., 2015). In a sequence of snow pits, Sterle et al. (2013) showed that during the ablation season, BC mixing ratios increase significantly near the snow surface (sampled in the top 2 cm) relative to bulk BC concentrations. They suggest that most likely a large fraction of previously deposited BC becomes concentrated near the surface. Delaney et al. (2015) also report surface BC increase during melt, to which BC being trapped at the snow surface is likely to contribute. BC increases in surface snow of up to an order of magnitude (Sterle et al., 2013; Doherty et al., 2016) and more (Xu et al., 2012) have been observed in natural snow during melt. Over most of the melt period, our results show a factor increase between 5 and 15 for the two-layer scenarios, which aligns well with observations. Higher values are mainly predicted shortly before meltout, when the snowpack is typically very thin and effects on discharge generation due to high increase in surface BC should be small.

We argue therefore the importance of providing, at a minimum, a separate surface layer, but recognize that simulated surface mixing ratios of BC are highly sensitive to the thickness of this layer. Since evaluation of model predictions for BC in snow is commonly performed by comparing simulated with observed BC mixing ratios in surface snow (e.g. Flanner et al., 2007; Forsström et al., 2013), this is a critical result. Snow is often sampled in the top few centimetres (typically 2 to 5 cm, e.g. Doherty et al., 2010; Aamaas et al., 2011; Forsström et al., 2013). This raises an interesting challenge given that the surface layer assumed in models is not a measurable property of snow. A comparison of model simulations with observations should therefore include some quantification of the uncertainty resulting from the surface layer thickness parametrization.

6.2 Hydrologic response to BC deposition in a snowfall dominated catchment

We are interested in addressing the impact of BC deposition – and potentially other light absorbing aerosols – on the hydrology of snowfall dominated catchments. Studies have shown the potential impact LAISI may have on the timing of snowmelt (Skiles et al., 2012; Painter et al., 2012), while others have argued that the impact on climate may be significant (Flanner et al., 2007, 2009; Qian et al., 2009, 2011). Given the importance of snow for water resources for a significant portion of the population (Barnett et al., 2005; Sturm et al., 2017) and the rapid growth of BC emissions in certain regions of the world (e.g. Paliwal et al., 2016; Bond et al., 2013), our aim is to provide a mechanism to include this process in hydrologic forecasting to better address future impact studies.

Forsström et al. (2013) found BC seasonal mean snowpack concentrations from about 10 to 80 ng g⁻¹ for different locations and time periods in mainland Scandinavia. Generally our results are within those presented in Forsström et al. (2013), though our max estimate lies above. However, Flanner et al. (2007) evaluated the global impact of the radiative forcing of BC in snow using a model that was compared with globally distributed surface BC measurements. For southern Norway, Flanner et al. (2007) predicted annual mean surface BC mixing ratios between 46 and 215 ng g⁻¹ for the year 1998, placing our simulations fully within a reasonable range of prior reported values.

The impact resulting from BC deposition in our study is seen in the timing of the annual water balance. Inclusion of ARF generally increases early season melt and causes the snowpack to melt out earlier. Comparing the ARF and no-ARF scenarios, we see a general shift in the discharge, with the ARF scenario producing greater discharge early in the season, and having less discharge after June. Such a shift in seasonal water balance will potentially have impacts on soil moisture and agriculture (Blankinship et al., 2014), as well as on regional climate (Qian et al., 2011). While we recognize significant uncertainties associated with conceptual hydrologic modelling that may impact the applicability of these results (Beven and Binley, 1992; see also the uncertainty discussion in Sect. 6.3), we feel it provides a novel mechanism to address LAISI in a manner that, to date, is not available otherwise. As a reality check of the catchment-scale process representation, we evaluate the impact of the incorporation of BC deposition on albedo, radiative forcing, and snowpack storage.

6.2.1 Surface BC and albedo

Albedo is a critical parameter in any snowmelt model, with significant control over the energy balance. During the accumulation period, the average albedo of each scenario lies within the range of albedo of fresh snow with small optical grain radius combined with a high solar zenith angle (Gardner and Sharp, 2010) and is thus reasonable for a high-latitude snowpack during snow accumulation. The differences in snow albedo during the accumulation season are mostly due to differences in aerosol deposition and in the maximum surface layer thickness of the snowpack. The time series of mid estimate modelled surface BC is within the range of values for locations in mainland Scandinavia presented in Forsström et al. (2013) during the accumulation period. The min estimate predicts values at the lower bound and lies in the range of the background surface BC level found in Svalbard in the European High Arctic (5 ng g⁻¹, Aamaas et al., 2011; 30 ng g⁻¹, Clarke and Noone, 1985). Compared to Forsström et al. (2013), the surface BC level of the max estimate seems to exceed the range of values reasonable for mainland Scandinavia during snow accumulation and reflects a range of values that is rarely found in snowpacks outside Asia (Doherty et al., 2010; Forsström et al., 2013; Wang et al., 2013; AMAP, 2015).

At the end of the melt season, the evolution of surface BC yields reductions in albedo relative to the no-ARF case of about 0.03, 0.1, and over 0.3 for the min, mid, and max estimates, respectively. This has two reasons: first, with increasing grain radius during the melt season, the absorbing effect of BC gets more efficient due to deeper penetration of radiation into the snowpack, leading to a stronger effect of the BC deposition on albedo. Snow of larger grains has a larger extinction coefficient and more effective forward scattering properties (Flanner et al., 2007). Second, with the start of the melt season there is a widespread decrease in snow thickness, allowing BC to accumulate in the surface layer. This latter effect is strongly dependent on the applied scavenging ratios, as we demonstrated in the 1-D sensitivity study (Sect. 5.1). During the melt season, the mid estimate spatially averaged surface BC mixing ratio increases from 49 to about 250 ng g⁻¹ (factor 5 increase) at the end of the melt season (beginning of July). Observations from Forsström et al. (2013) indicate that surface BC mixing ratios around 250 ng g⁻¹ are well within the range of reasonable values for a melting Scandinavian snowpack. Furthermore, an increase in surface BC by a factor of 5 and higher during snowmelt is in line with observed BC trends in melting snow from different locations (Doherty et al., 2013, 2016; Xu et al., 2012). From this, we argue that our mid estimate simulation predicts a seasonal cycle in surface BC that is within reason.

For the max estimate, the increase is from roughly 100 to over 2500 ng g⁻¹ (factor 25 increase). This strong seasonal cycle in surface BC is beyond what is observed for both absolute BC values in Scandinavian snowpacks and increase relative to surface BC during snow accumulation. The min estimate, on the other hand, leads to a decrease in the BC surface mixing ratio. Even though many studies report an increase in surface BC during snowmelt (e.g. Conway et al., 1996; Doherty et al., 2013, 2016; Xu et al., 2012), there exist observations showing that a large fraction of BC can be flushed efficiently from the snowpack with the beginning of snowmelt (Lazarcik et al., 2017). This indicates that post-depositional enrichment processes and their significance in determining surface BC trends in melting snow require further exploration. We argue that the min estimate thus marks a reasonable lower bound estimate for the seasonal evolution of surface BC.

We recognize our max estimate results in a strong increase in surface BC mixing ratios mostly due to low BC scavenging with melt (note the strong increase from the end of March on in Fig. 6). This divergent evolution of surface BC mixing ratios in the min, mid, and max estimates reveals uncertainty in the representation of the fate of BC in snow during melt, which is also reflected in the literature (Doherty et al., 2013, 2016; Xu et al., 2012; Lazarcik et al., 2017).

6.2.2 BC-induced radiative forcing

The strong increase in RFS (Fig. 7a) and surface radiative forcing (Fig. 7b) during spring melt results from the combination of (i) the aforementioned decrease in snow albedo due to the increase in surface BC mixing ratios (e.g. melt amplification and the increasing optical grain radius in melting snow as discussed in Sect. 5.2.2) and (ii) the increasing daily solar irradiation due to a lower solar zenith angle and longer days.

The annual mean surface radiative forcings in this study are 0.284, 0.844, and 1.391 W m⁻² for the min, mid, and max estimates. Averaged over Scandinavia (including Finland), Hienola et al. (2016) calculated lower values around 0.145 W m⁻². However, Hienola et al. (2016)'s study includes large areas with shorter snow cover. Since the surface radiative forcing is strongly dependent on the snow cover evolution, higher values compared to Hienola et al. (2016) are expected due to the long lasting snow cover in our case study region. The mid estimate annual cycle of surface radiative forcing due to the presence of BC in the study region is of a similar magnitude to what is found over the Tibetan Plateau. Qian et al. (2011) report similar snow cover duration and maximum mean forcing during May of over 6 W m⁻² using a global climate model. Due to the generally much lower snow-covered fraction in Qian et al. (2011)'s study region, however, RFS is presumably significantly higher on the Tibetan Plateau compared to our study region, which is in agreement with very high levels of BC reported for the Tibetan Plateau (Qian et al., 2011). Using a stand-alone version of SNICAR, we estimated RFS based on surface BC mixing ratios from Forsström et al. (2013) measured during melt in the top 5 cm of Scandinavian snowpacks to 4.7 to 18.2 W m⁻² (95 % confidence interval; details described in Appendix Appendix A). These values agree well with our min and mid estimate RFS (Fig. 7), but are significantly lower than our max estimate predictions.

6.2.3 BC impact on catchment discharge and snow storage

We mention a shift in the seasonal water balance, with more melt early in the melt season resulting from enhanced RFS. However, from mid May the melt enhancement reduces and the differences in catchment SWE between the ARF and no-ARF scenarios decrease (Fig. 8b). One would expect with more incoming radiation, later in the season, the RFS effect to become further enhanced. However, this counter-intuitive result becomes clearer when one considers the impact of fractional snow-covered area and catchment-scale processes. The dynamics driven by the faster development of SCF (see Fig. 6a) is a limiting factor in the catchment-averaged snowmelt. By comparing Fig. 7a, which shows the RFS enhancement, with Fig. 7b, which shows total daily energy uptake in the catchment, we see that a threshold period is reached and total daily energy uptake decreases, while RFS is continually increasing. The SCF decrease with increased melt due to ARF counteracts the RFS effect itself, due to the reduction in area from which snow can actually melt. For discharge, this is manifested in the ARF scenarios as an enhancement during the beginning of the melt season attributed to RFS, whereas the decreased discharge later in the season is attributed to melt limitation caused by the faster growth of fractional bare ground areas.

Similar shifts in the annual water balance due to the impact from LAISI are reported for the Upper Colorado River Basin (Painter et al., 2010) and the Tibetan Plateau (Qian et al., 2011). Those regions are well-known hotspots of LAISI disturbance to snow cover (Painter et al., 2007; Qian et al., 2014). Our results suggest that the hydrologic cycle of regions that have not been the focus hitherto (such as Norway) might also be significantly affected by ARF.

Compared to observations, all simulations (ARF and no-ARF) tend to underestimate discharge during early melt season and overestimate discharge during late melt season (Fig. 8a). However, the magnitude of over- and underestimation strongly differs between the scenarios. By including ARF the volume error is reduced in both the early melt season (by increasing melt), and in late melt season (by subsequently decreasing melt generation in the catchment due to reduced SCF). Expressed as seasonal mean volume error for early and late melt season, the difference to observed discharge is largest for the no-ARF scenario and smallest for the max estimate. The max estimate reduces the volume error by −75.1 % during early melt season and −89.9 % during late melt season, relative to the no-ARF scenario (see Table 4). The min and mid estimates also reduce the volume error. Thus, on average, an improvement in simulated discharge is achieved during the melt season by accounting for BC RFS. Similar results are achieved when estimating model parameters using a no-ARF scenario (not shown). However, we acknowledge that further studies are needed in order to be able to confirm a general model improvement when accounting for ARF in snow dominated catchments. Certain mechanisms can lead to model improvements for the wrong reason when applying ARF (Kirchner, 2006). Structural deficits of the model might lead to a negligence of processes that are important for the spring melt generation. The implementation of ARF could then optimize the model towards the observations and counteract errors coming (partly) from a missing process that is not related to ARF. A further potential mechanism is related to the equifinality of conceptual models. These implications coming from model parameter uncertainty are discussed in Sect. 6.3 alongside with further sources of uncertainty.

Table 4 Season mean volume error in discharge during the early (22 March to 29 May) and late (30 May to 10 August) melt seasons of the no-ARF, min, mid, and max scenarios compared to observed discharge. The percentage change shows an increase (+) or decrease (−) in the volume error compared to the no-ARF volume error.

6.3 Uncertainties

There are numerous challenges associated with the development of an algorithm that mixes conceptual hydrologic parametrizations with physically based approaches. Both the literature and our analysis highlight aspects that warrant a deeper investigation of ARF-induced uncertainty. The intent with this work is to introduce a new algorithm; however, as indicated in Pappenberger and Beven (2006), we feel it is important to provide an initial assessment of the uncertainty introduced with the addition of ARF terms. To achieve this we have conducted a generalized likelihood uncertainty estimation (GLUE; Beven and Binley, 1992) which provides an assessment of the degree of variability in behavioural models resulting from equifinality.

With respect to the implementation of a physical albedo model, the treatment of the darkening effect of LAISI adds additional degrees of freedom to the parameter space of the model due to the introduction of new parameters (scavenging ratios, surface layer thickness, BC input scaling factor; see the bottom four parameters in Table 2). In order to investigate the abilities and limits of the model with and without ARF to reflect the observed discharge, we quantify the parameter uncertainty prior and posterior to the implementation of ARF calculations (Fig. 9; details in Appendix Appendix B). Uncertainties are generally largest during snowmelt and summer because various parameters only play an active role in calculating discharge during snowmelt. Including ARF calculation in the model leads to a shift of the uncertainty band to higher values during April and May, and lower values during June and July, due to increased melt under the impact of ARF. From mid May to mid June, the ARF-induced shift in the uncertainty band leads to observations being within or closer to the border of the uncertainty bands (shaded box in Fig. 9), which can be interpreted as an improvement to the model. This would imply that in the model without ARF, albedo decays not sufficiently enough during spring in order to generate enough snowmelt, resulting in an underestimation of discharge in April and May. However, we admit that further testing is needed to draw a more accurate conclusion, as discussed above. Perhaps more importantly, it appears that we have not increased uncertainty much by adding complexity. In general, both simulations with and without ARF lead to acceptable results. However, we enable the inclusion of a potentially important variable, particularly with respect to increasing emissions of light absorbing aerosols due to population growth.

Figure 9 The 95 % confidence interval of simulated discharge due to parameter uncertainty when allowing for ARF (red) and disregarding ARF (grey), calculated using the GLUE method and averaged over the 6-year simulation period. The shaded box marks the period of the melt season, where observations tend to lie outside the uncertainty bounds of the no-ARF simulations.

In our case study, further uncertainties result from mixing ratios of BC in the snowpack due to prescribed BC deposition, and LAISI other than BC not accounted for in the simulations:

• i.

  prescribed BC deposition

  In the approach presented here, we use prescribed BC deposition mass fluxes. Even though this is common practice (e.g. Goldenson et al., 2012; Lee et al., 2013; Jiao et al., 2014), it was shown by Doherty et al. (2014) that the decoupling of aerosol deposition from the water mass flux of falling snow can lead to an overestimation of surface mixing ratios by a factor of 1.5–2.5. However, we would like to highlight an important difference between our approach and the one Doherty et al. (2014) claim to be problematic: first, the high bias in surface snow BC mixing ratios described by Doherty et al. (2014) refers to global climate model simulations with prescribed aerosol deposition rates (wet and dry), where the input aerosol fields are interpolated in time from monthly means. Therefore, the episodic nature of aerosol deposition due to wet deposition is generally absent in the prescribed aerosol fields. The coupling of the interpolated fields with highly variable meteorology (in particular precipitation) results in the high bias (Doherty et al., 2014). In our case study, we use deposition fields originating from the REMO-HAM regional aerosol climate model, forced with ERA-Interim reanalysis data at the boundaries. REMO-HAM output is 3-hourly, which we re-sampled to daily means in order to have consistency between the deposition fields and the observed daily precipitation used as input data in the hydrological simulations. The daily time step allows us to preserve the episodic nature of aerosol deposition. Moreover, the daily BC wet deposition rates should not be biased due to major inaccuracies in precipitation, as REMO-HAM has been shown to reproduce the Scandinavian precipitation realistically (Pietikäinen et al., 2012). The high bias occurring when using interpolated monthly averages as input should therefore be minimized. Additionally, and significantly, Doherty et al. (2014) (and the critiques therein) address an objective with consideration to climate impacts. Our analysis is focused on the impact on the hydrological cycle. Our simulations suggest that BC RFS is mostly important during spring time, where surface BC mixing ratios are predominantly controlled by melt processes, and not by deposition processes (as shown in Figs. 3 and 6b).

• ii.

  LAISI other than BC

  By including only BC deposition in our simulation, we likely underestimate the additional effect of further LAISI species such as mineral dust (Di Mauro et al., 2015; Painter et al., 2010), mixing of the snow with soil from the underlying ground or local sources (Wang et al., 2013), and biological processes (Lutz et al., 2016). Neglecting additional RFS from LAISI other than BC is likely to result in an underestimation of the overall effect of LAISI on snowmelt and discharge generation. Especially the contribution from dust is critical since it has been shown that in many regions such as the Rocky Mountains (Painter et al., 2012), Utah (Doherty et al., 2016), the southern edge of the Himalayas (Gautam et al., 2013), and Svalbard (Forsström et al., 2013), dust can play a significant role in terms of RFS or even is the dominating LAISI. For Norway, however, analysis conducted by Forsström et al. (2013) indicates that dust might only play a minor role. By comparing samples from Svalbard and near Tromsø, Norway, Forsström et al. (2013) showed that there exits a distinctive difference between the Arctic Archipelago and the mainland. The BC mixing ratio from mineral-dust-rich Svalbard measured by the thermal/optical method used in Forsström et al. (2013) averaged about half the mixing ratio of insoluble light absorbing particulates (including dust) measured by an optical method (ISSW: Integrating Sphere/Integrating Sandwich; e.g. Doherty et al., 2010). Samples collected close to Tromsø, on the other hand, resulted in BC that averaged about 1.3 times the ILAP mixing ratios. Due to the fact that the ISSW method overestimates BC for samples containing dust, Forsström et al. (2013) argue that the comparison of both methods can be used to draw conclusions about the pollution regime. However, due to the small number of samples and the single-location analysis, this needs to be addressed more in future studies in order to identify the relative importance of different LAISI species.

With respect to our study, we acknowledge that including only BC is a shortcoming with respect to the overall effect of LAISI. However, by demonstrating the significant effect of BC on accelerating snowmelt and discharge generation, our study gives a conservative estimate of the effect of LAISI and urges a more detailed investigation.

7 Conclusions

Herein we presented a newly developed snow algorithm for application in hydrologic models that allows a new class of model input variable: the deposition rates of light absorbing aerosols. By coupling a radiative transfer model for snow to an energy-balance-based snowpack model, we are providing a tool that can be used to determine the effect of various species of LAISI at the catchment scale. In this analysis we have focused solely on BC and acknowledge it therefore likely represents a conservative estimate. This work presents a novel analysis of the impact of BC deposition to snow on the hydrologic cycle through 1-D sensitivity studies and catchment-scale hydrologic modelling. From a 1-D model study, presented in Sect. 5.1, we conclude that

• i.

  the implementation of at least two layers (a thin surface layer and a bottom layer) is of outstanding importance to capture the potential effect of melt amplification on the near-surface LAISI evolution. The parametrization of the surface layer has only a small effect on the snow albedo and melt rate as long as the surface layer thickness (in SWE) is sufficiently thin (e.g. thinner than the penetration depth of shortwave radiation). However, the evolution of the LAISI surface mixing ratio is highly sensitive to the surface layer thickness. For this reason, we suggest including a surface layer thickness variation in model studies when comparing simulated and observed LAISI mixing ratios sampled in the top few centimetres of snow.

• ii.

  The determination of how LAISI are washed out of the snowpack with meltwater has a great effect on the evolution of LAISI concentration near the surface, snow albedo, and melt rate. Due to rare observations of this effect under controlled conditions, the uncertainties are high and our findings show the need for more detailed understanding of the processes involved due to the high importance of the overall effect of LAISI on the snowpack evolution.

To demonstrate the significance of BC radiative forcing for the hydrologic cycle at the catchment scale, we demonstrated the effect of BC deposition and the subsequent implications for snowmelt and discharge generation on a remote Norwegian mountain catchment. The study indicates that inclusion of BC in snow is likely to have a significant impact on melt timing, and that the effect on the discharge generation leads to a shift in the annual water balance. Our simulations further suggest that melt amplification can have severe implications for both the snowpack evolution and the discharge regime of a catchment, which means that the seasonal cycle of the surface BC mixing ratio is of great importance. However, large uncertainties are connected with the representation of surface enrichment of BC. A more robust understanding of the fate of BC in melting snow is essential to fully assess impacts on the hydrologic cycle.

Including radiative forcing from BC in the simulations leads to a reduction in volume error during the early and late melt season in our simulations. We conclude from our study that hydrological modelling can potentially be improved by including the effect of LAISI, especially when the model approach implies a physically based representation of the snowpack in general and the snow albedo in particular. However, more research in the area of catchment-scale impact of LAISI is needed to support this. The approach and algorithm presented in this analysis provide a tool to target this in future applications.

Appendix A: Radiative forcing in snow estimated from Forsström et al. (2013)

In order to calculate radiative forcing in snow (RFS) from surface concentrations during melt reported in Forsström et al. (2013), several assumptions have been made. For each input variable, a certain reasonable range is estimated, suited to snow properties during melt conditions:

• snow optical grain radius: 500–1000 µm;

• snow density: 400–600 kg m⁻³;

• BC mixing ratio: 50–200 ng g⁻¹ (from Forsström et al., 2013).

Forsström et al. (2013) report six time series of BC surface concentrations sampled in the top 5 cm of the snowpack, all of which cover the snowmelt period at three locations in Scandinavia; however, only one location can be considered remote without pollution from local sources (Abisko, Sweden). The range of BC mixing ratios during melt is estimated from this location. Global radiation during spring is estimated to 210 W m⁻². The value has been calculated from the input time series of our study region, in order to receive comparable results. The daily mean solar zenith angle has been set to 60^∘ and BC mixing ratios below the top 5 cm to 0 ng g⁻¹, since no further information is available. The latter might lead to an underestimation of RFS, and results can be seen as a conservative estimate; 1000 realizations with SNICAR have been conducted using different input variable sets, with random values for each input variable according to a uniform distribution in the stated range. Resulting RFS values are presented as the 95 % confidence interval to 4.7 to 18.2 W m⁻². The mean is 11.2 W m⁻².

Appendix B: Parameter uncertainty with GLUE

We determine parameter uncertainty using the GLUE method (Beven and Binley, 1992). Lower and upper bounds of parameters used in the calculation are shown in Table B1. We use the Nash–Sutcliffe model efficiency (see Eq. 9) as a likelihood function and choose a threshold value of 0.74 (0.1 below the best calibration result) for accepting parameter sets as behavioural parameter sets. To identify the impact of ARF on model uncertainty, we run GLUE twice, first without ARF applied, and in a second round of simulations accounting for ARF. Random parameter sets are created by choosing parameters according to a uniform distribution in the range of the parameter bounds. For each of the two uncertainty estimations, a total of 10 000 model realizations was drawn, of which 1435 (no-ARF) and 1831 (ARF) parameter sets were rated as behavioural parameter sets. This accounts for about 14 and 18 % of the total samples, respectively.

Table B1 Model parameter bounds used in the uncertainty estimation with the GLUE method. Parameters used to determine ARF are marked with ^*.