Measurement of unsaturated meltwater percolation flux in seasonal snowpack using self-potential

Abstract

This paper presents a feasibility study of in situ field measurements of unsaturated meltwater percolation flux within the vertical profile of a snowpack, using the self-potential (SP) method. On-site snowmelt column tests calibrated the SP measurements. The SP data measured electrical field strength with an electrode spacing of 20 cm, and coincident water saturation (S_w) measurements using time domain reflectometry allowed calculation of SP-modeled vertical percolation flux (q_sp), expressed as Darcy velocity. The results reflected transient diurnal snowmelt dynamics, with peak flux lagging arrival of a saturation wetting front. Peak daily q_sp was 60 to >300 mm d⁻¹, whereas daily snowmelt was 20–50 mm w.e. Surface refreezing events appeared to cause upward flow, possibly representing water redistribution toward the freezing boundary. Calculated fluxes were comparable to actual fluxes, although average errors ranged from −15 to +46% compared to average of melt expected from surface energy-balance and ablation stake measurements. By advancing method development to measure unsaturated meltwater percolation flux in snowpacks this study creates opportunities to study fundamental snowmelt processes, may improve mathematical modeling and may supplement glacier mass-balance studies and studies of snowmelt interactions with avalanches, groundwater and surface water.

1 Introduction

Approximately one-sixth of the world's population relies on snowmelt or glacier runoff for their water supply, and in these areas the timing and magnitude of snowmelt exerts a major influence on water availability, and other potential effects of drought and climate change (Barnett and others, 2005). For seasonal snowpacks, the timing, rate and magnitude of meltwater percolation and retention in the snowpack are important factors controlling recharge of groundwater systems and surface runoff (Smith and others, 2014; Musselman, and others, 2017). Wet snow avalanches are a threat associated with the onset of meltwater percolation to depth (Wever and others, 2018), which triggers failure of weak snowpack layers. For glaciers, percolation of meltwater is an important factor affecting glacier mass balance (Jansson and others, 2003), particularly with respect to the lag between melting and percolation, the retention and storage of meltwater within the snowpack or firn and the delayed release of meltwater deeper into the glacier hydrologic system. The recognition of firn aquifers, which store substantial volumes of meltwater in portions of the Greenland ice sheet (Forster and others, 2014; Miège and others, 2016), has placed further emphasis on the importance of understanding meltwater percolation phenomena as part of overall glacier mass balance.

Liquid water movement through the snowpack leads to a lag between the initiation of melting, the arrival of meltwater at depth and the elution of meltwater from the snowpack. Meltwater percolation in snowpacks is a complex process that involves unsaturated flow, where the pore space is only partially saturated with water. The theoretical underpinnings for the analysis of unsaturated flow during meltwater percolation in snow were historically founded on previous research in soil physics and petroleum engineering (e.g. Richards, 1931; Brooks and Corey, 1964; Bear, and others, 1968, and more). Colbeck (1972) first applied then-existing unsaturated flow theory to the physics of downward propagation of wetting fronts in snow during diurnal snowmelt cycles, and later Colbeck (1975) extended the theory to layered snowpacks. To address different problems, researchers have discretized the snowmelt percolation problem at various scales ranging from bulk meltwater outflow models (Essery and others, 2015), 1-D layered snowpack formulations (Wever and others, 2014) and centimeter-scale 3-D mathematical models (Hirashima and others, 2014; Meyer and Hewitt, 2017). The pore-scale structure of the snowpack is understood to strongly influence the unsaturated hydraulic properties (Shimizu, 1970; Colbeck, 1975; Yamaguchi and others, 2010). Accordingly, the literature shows that meltwater percolation involves development of preferential flow paths and large spatial variability in meltwater flux both due to snowpack layering as well as capillary effects (Colbeck, 1979; Marsh and Woo, 1985; Williams and others, 2010; Avanzi and others, 2016; Webb and others, 2018). Mathematical models of unsaturated flow during snowmelt percolation have increased in complexity over time, incorporating and coupling processes such as thermal effects and re-freezing, capillary forces, hydraulic processes and multi-dimensional preferential flow (Illangasekare and others, 1990; Tseng and others, 1994; Hirashima and others, 2014; Wever and others, 2014; Meyer and Hewitt, 2017).

The research reported here involves in situ measurement of the Darcy velocity of snowmelt percolation at the decimeter scale in a vertical profile within the snowpack. The pertinent unsaturated flow physics provides context for the alternatives available to measure snowmelt percolation processes. At any representative elementary volume within the downward flow field of meltwater percolation, and dependent on snowpack structure and properties, there are interrelationships between: (1) the amount of liquid water that occupies the pore space, (2) the capillary and gravitational hydraulic forces on the water and (3) the resulting Darcy velocity, or flux, of water flow. These unsaturated flow parameters are distinct, and represent separate measurement opportunities. This study focuses on the measurement of Darcy velocity, although we will briefly review the background related to other measurements of unsaturated water in snow.

The field measurement of capillary pressure using porous cup tensiometers in soil is well established (Livingston, 1908; Campbell, 1988), although described to a more limited extent in snow (Colbeck, 1976; Wankiewicz and De Vries, 1978; Clayton, 2017). Researchers have reported measurement of volumetric liquid water content (LWC) in snow at various scales using several methods (Kinar and Pomeroy, 2015). Methods with a measurement scale on the order of centimeters have included dielectric methods (Denoth, 1989), time domain reflectometry (TDR) (Schneebeli and others, 1998; Diaz and others, 2017; Samimi and Marshall, 2017) and microwave permittivity (Mavrovic and others, 2020). Techniques to measure LWC at the bulk snowpack scale include upward looking ground-penetrating radar (Heilig and others, 2015) and GPS signal attenuation (Koch and others, 2019).

Although the field measurement of LWC in snow is well established, the literature does not describe the field measurement of Darcy velocity in snow during unsaturated meltwater percolation. At the lab scale, Walter and others (2013) measured pore-scale water percolation velocities in snow in a laboratory study using micrometer-sized fluorescent particles. Snowmelt lysimeters, which collect meltwater from the base of the snowpack, often over a scale of several m² or more, are effective tools for monitoring bulk meltwater production (Haupt, 1969; Kattelmann, 2000; Webb and others, 2018). Kulessa and others (2012) used the electrical self-potential (SP) geophysical method in a laboratory snowmelt column for the purpose of estimating the vertical Darcy velocity during meltwater percolation in snow. Building upon the study of Kulessa and others (2012), this research involved development and deployment of field methods using SP to determine Darcy velocity at discrete locations in a vertical profile within a glacier seasonal snowpack during snowmelt percolation. Clayton (2017) presented a discussion paper that described initial field pilot-testing of the measurement of vertical percolation flux using the SP-field method in snow.

Other uses of SP in glacial environments have been reported that evaluated lateral flow of meltwater. Kulessa and others (2003) measured SP in subglacial meltwater drainage systems to assess flow magnitude and direction. Thompson and others (2016) collected SP measurements in a 25 m × 25 m grid on a glacier surface to assess bulk glacier meltwater production and lateral flow in thin seasonal snowpack overlying glacier ice. They also measured bulk meltwater discharge at a downgradient lysimeter and found a good qualitative temporal correlation between SP signals and diurnal meltwater production.

The use of SP to determine flow in porous media is based on the measurement of a millivolt (mV)-level electrical potential induced by the movement of pore water through a medium with a surface charge (Revil and others, 1999; Linde and others, 2007; Revil and others, 2007; Revil and Jardani, 2013). The surface charge of the solid matrix becomes balanced against an accumulation of oppositely charged ions in the pore fluid, and flow of these charged ions in the pore fluid results in an electrical field of a strength that is proportional to the flow velocity and other fluid and porous media properties. The magnitude of the difference in charge between the solid surface and the fluid is referred to as zeta potential. Zeta potential is a sensitive function of pore water pH and electrical conductivity, which is important in snowmelt systems (Kulessa and others, 2012). They inferred a reversal in sign of zeta potential as a function of pH and electrical conductivity, with zeta potential values between −0.1 and +0.02 for electrical conductivities below ~10⁻⁴ S m⁻¹.

In the SP method, the electrical potential attributed to fluid flow is referred to as streaming potential. Streaming potential is distinguished from possible electrochemical or thermoelectric potentials resulting from chemical or thermal gradients in pore water. Streaming potentials are of primary importance in SP applications such as characterizing subsurface flow in groundwater hydrology (Soueid Ahmed and others, 2016), while in many applications of SP, electrochemical or thermoelectric potentials are of primary importance. For example, electrochemical potentials are important in using the SP method for mineral exploration (Biswas, 2017), and thermoelectric potentials are important in geothermal systems (Corwin and Hoover, 1979). For isothermal snowmelt systems, electrochemical and thermoelectric potentials are expected to be small (Kulessa and others, 2012).

2 Objectives

The objectives of the exploratory research reported herein were to develop field methods and implement a feasibility study using SP in the field to measure the decimeter-scale unsaturated flow velocity in a vertical profile during meltwater percolation in snow. The focus of this study was on two areas related to field application: (1) evaluate the potential quantitative accuracy of the method, and limitations thereof, and (2) evaluate the methods' capabilities to qualitatively evaluate flow transients during diurnal melting cycles and their relationships to short-term changes in weather and energy-balance drivers of snowmelt. SP was not evaluated as a tool for the collection of long-term (months to seasons) data. As discussed above, the literature reports the use of mathematical models to evaluate snowmelt percolation and field measurement of LWC (alternately expressed later as water saturation (S _w = LWC/porosity)) within the snowpack. However, LWC measurements may not reflect transient changes in meltwater percolation velocity, since they are distinct parameters. The measurement of meltwater percolation flow velocity offers the opportunity to directly observe transient flow behaviors that may not be understood well or may not be accounted for within the framework of mathematical models. The literature does not describe techniques for the field measurement of the transient unsaturated flow velocity within the snowpack, and new field methods will inform the understanding of important fundamental meltwater percolation processes.

3 Methods

3.1 Unsaturated flow formulation

This study explores using SP-field measurements to calculate meltwater percolation flux vertically downward within the snowpack under unsaturated flow conditions. The unsaturated flow formulation builds upon the study of Kulessa and others (2012), who extended SP theory to the problem of gravity-driven unsaturated flow in a vertical melting snow column. They derived an equation to describe the SP signal generated by the unsaturated flow percolating through the snowmelt column, which is shown below (Eqn (1)), with an added correction (personal communication from Kulessa, A. and Revil, A., 2020) to add the measurement length term (L) to the numerator:

(1)$$\psi _{\rm m} = \displaystyle{{\varepsilon \zeta } \over {\sigma _{\rm w}}}\displaystyle{{S_{\rm w}} \over {S_{\rm e}^n }}\displaystyle{L \over {kA}}Q$$

where ψ _m is the streaming potential between measurement and reference electrodes (V); ɛ is the meltwater dielectric permittivity (F m⁻¹) (constant 7.8 × 10⁻⁹ Farads per m (Kulessa and others, 2012)); ζ is the zeta potential (V); σ _w is the meltwater electrical conductivity (S m⁻¹); S _w is the water saturation (=LWC/porosity); S _e is the effective water saturation; n is the pore size parameter; k is the intrinsic permeability (m²); A is the column cross-sectional area; Q is the volumetric flow rate (m³ d⁻¹) and L is the distance between reference and measurement electrodes (m).

Kulessa and others (2012) verified Eqn (1) by using it to model SP signals and comparing the results to a laboratory snowmelt column experiment. Eqn (1) shows that the relationship between streaming potential (ψ _m) and the fluid flow rate is affected by a significant number of variables that will have to be measured or estimated in the field to determine flow rate from SP data. To solve this problem, the methods included the use of an in-field snowmelt column test to estimate certainty of the parameters in Eqn (1) which were then applied to calculate percolation fluxes within the snowpack. We will assume that thermoelectric and electrochemical potentials are zero, and the appropriateness of these assumptions is addressed in subsequent sections.

We can calculate the meltwater flux, or Darcy velocity, as a function of the streaming potential and other parameters, by rearranging Eqn (1), as follows:

(2)$$q_{{\rm sp}} = \displaystyle{{\psi _{\rm m}/L} \over {( \varepsilon \zeta /\sigma _{\rm w}) ( S_{\rm w}/kS_{\rm e}^n ) }}$$

where q _sp = Q/A (m d⁻¹).

The SP-modeled meltwater Darcy velocity, or flux (q _sp), has units of mm d⁻¹ w.e., although for brevity we will omit the term w.e. and report q _sp as mm d⁻¹. For this study, we measured SP signals in mV, and report the electrical field strength (E, mV m⁻¹) measured across two electrodes of spacing, L (m):

(3)$$E = \psi _{\rm m}/L$$

We measured the SP electrical field strength directly between paired electrodes, without any distant isolated reference electrode. We refer to the electrical field strength hereinafter as SP-field strength or simply SP.

Following the unsaturated flow formulation used by Kulessa and others (2012), we use the Brooks and Corey (1964) constitutive relationships, as follows:

(4)$$S_{\rm e} = \displaystyle{{( {S_{\rm w}-S_{\rm r}} ) } \over {( {1-S_{\rm r}} ) }}$$

(5)$$k_{\rm r} = S_{\rm e}^n $$

(6)$$S_{\rm e} = \left({\displaystyle{{H_{\rm o}} \over {H_{\rm c}}}} \right)^\lambda $$

where S _r is the residual water saturation, k _r is the relative permeability, H _c is the capillary pressure head, H _o is the air entry capillary head and n = (2 + 3λ)/λ.

In Eqns (4) and (5), the exponents n and λ are characteristic of the pore size distribution of the porous media that control the relationships between permeability, saturation and capillary pressure. Pore size distribution also controls snow permeability (k, m²) which we estimated from density and grain size (Colbeck and Anderson, 1982; Kulessa and others, 2012) using the following, from Shimizu (1970):

(7)$$k = 0.077d^2{\rm e}^{{-}0.0078\rho _{\rm s}}$$

where d is the mean snow grain diameter (m) and ρ _s is the snow density (kg m⁻³).

3.2 Field study location

A number of students from the Juneau Icefield Research Program (JIRP) and myself undertook this field feasibility study during two separate short summer field sessions, each ~1 week in duration, in 2018 and 2019. We conducted the field study at the Matthes-Llewellyn glacier divide (Fig. 1), ~5 km north of the border between Alaska, USA, and British Columbia, Canada. The site is located at a transition between a predominately maritime climate to the south and west and a predominately continental climate to the north and east. The Juneau Icefield is a temperate glacier system, and during the summer melt season the glacier system is isothermal, with snow and ice temperatures at 0°C throughout. The elevation at the study site is ~1850 m, and since the location is a glacier divide the surface slope is essentially flat. JIRP conducted ice-penetrating radar studies during 2018 and 2019 that indicated the glacier thickness was ~900 m at the site. Nearby glacier mass-balance study pits in 2018 and 2019 indicated that the late July depth of the seasonal snow layer ranged from ~3–4 m. Matthes Glacier flows south as part of the Taku Glacier system, which was the subject of a glacier mass-balance reanalysis study (McNeil and others, 2020) that further describes the climatological and glaciological setting of the Juneau Icefield.

Fig. 1. Juneau Icefield study site location map. Source: Esri, Maxar, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN and the GIS User Community.

3.3 Weather and energy-balance measurements

We measured weather and energy balance to characterize conditions during the fieldwork that would drive snowmelt and meltwater percolation. Energy-balance data served as a basis calculated expected melt for comparison to SP-modeled percolation fluxes, as described in Section 3.7. A temporary meteorological station was co-located with the snowmelt percolation study which measured wind speed and direction, temperature, relative humidity, as well as four components of surface energy balance including longwave (LW) and shortwave (SW) radiation inputs and outputs. We installed a visually read ablation stake to monitor actual loss of surface snow, consisting of an aluminum avalanche probe with 0.5 cm markings. Table 1 summarizes the meteorological station measurement equipment specifications. A Campbell Scientific CR-10x data logger with 12 V solar power supply, housed in a waterproof Pelican™ case logged all data acquisition at 15 s intervals with 20 samples averaged over each recorded 5-min data period.

Table 1. Sensors used in the meteorological station

The balance of energy (Q, W m⁻²) at the snow surface drives snowmelt. This includes LW and SW radiation inputs and outputs as well as energy associated with air temperature, evaporation/condensation and precipitation (Hock, 2005), as follows:

(8)$$Q_{\rm N} + Q_{\rm H} + Q_{\rm L} + Q_{\rm G} + Q_{\rm R} + Q_{\rm M} = 0$$

where Q _N is the net radiation energy flux (solar); Q _H is the sensible heat energy flux (temperature); Q _L is the latent heat energy flux (evaporation/sublimation/condensation); Q _G is the ground heat energy flux (geothermal); Q _R is the rain heat energy flux (precipitation) and Q _M is the energy consumed by melt.

We can neglect ground heat energy flux, based on the isothermal glacier characteristics. Sensible and latent heat energy fluxes are categorized as turbulent energy fluxes because they rely on near-surface atmospheric turbulence for energy transfer between the snow surface and the atmosphere. Turbulent energy fluxes are generally less than net radiation over long time periods, but are significant contributors to short-term transient variation in melting (Hock, 2005). Turbulent energy fluxes are challenging to estimate using bulk methods, because of uncertainty in using a transfer coefficient to represent turbulence profiles that are spatially and temporally variable (Hock, 2005). Nonetheless, to support a more complete estimate of expected transient meltwater percolation flux based on surface energy balance, we estimated Q _H according to the methods used in the factorial snow model (Essery, 2015), as follows:

(9)$$Q_{\rm H} = \rho _{\rm a}c_pC_{\rm H}U_{\rm a}( {T_{\rm s}-T_{\rm a}} ) $$

where ρ _a is the air density; c_p is the heat capacity of air = 1005 J K⁻¹ Kg⁻¹; C _H is the dimensionless transfer coefficient; U _a is the wind speed at measurement height z; T _s is the surface temperature and T _a is the air temperature at measurement height z; and where

(10)$$C_{\rm H} = f_{\rm H}k^2\left[{\ln \displaystyle{{z_{\rm U}} \over {z_0}}\;{\rm ln}\displaystyle{{z_{{\rm TQ}}} \over {z_{{\rm oh}}}}} \right]^{{-}1}$$

where f _H is the atmospheric stability function ( = 1.0); k is the von Karman constant ( = 0.4); z _U = z_TQ is the measurement height for wind speed and humidity; z ₀ is the surface momentum roughness length and z _oh = 0.1z ₀.

The measurement height for wind and humidity at the meteorological station was 1.0 m. Surface momentum roughness length for snow surfaces can vary over a wide range of values, from 0.004 to 70 mm (Hock, 2005). We use values of surface momentum roughness of z ₀ of 10^−4.6 m for latent heat transfer and 10⁻⁶ m for sensible heat transfer, taken from an eddy covariance study at a mountain glacier in western Canada by Fitzpatrick and others (2017).

Following Essery (2015), latent heat flux was determined as follows:

(11)$$Q_{\rm L} = L_{\rm e}S$$

where L _e is the latent heat of evaporation/condensation = 2.835 × 10⁶ J kg⁻¹ and S is the mass rate of vapor transfer (kg m⁻² S⁻¹), and where

(12)$$S = \rho _{\rm a}C_{\rm H}U_{\rm a}[ {Q_{{\rm sat}}-Q_{\rm a}} ] $$

where Q _sat is the specific humidity at water vapor saturation, at melting snow surface and Q _a is the specific humidity of air at measurement height.

We used Eqns (8) through (12) and field meteorological measurements to calculate expected melt rate (M, in m w.e. s⁻¹) (Eqn (13)) (Hock, 2005). In conventional units, Eqn (13) reduces to: Q _M (W m⁻²) × 0.26 = M (mm w.e. d⁻¹):

(13)$$M = \displaystyle{{Q_{\rm M}} \over {\rho _{\rm w}L_{\rm f}}}$$

where ρ _w is the density of water (kg m⁻¹) and L _f is the latent heat of fusion = 334 000 J kg⁻¹.

3.4 Water retention curve measurement

We measured water retention curves relating capillary pressure and water saturation as a baseline metric of unsaturated hydraulic characteristics of the seasonal snowpack, using a sample from Taku Glacier on 21 July 2018, and a sample from the Matthes-Llewellyn divide on 25 July 2018. The testing apparatus (Clayton, 2017) involved an 8 cm diameter isothermal test cell to maintain the snow core sample at 0°C, which prevented melting during the drainage experiment. Fitting the capillary pressure–saturation curves in Eqn (6) determined the Brooks and Corey (1964) pore size distribution parameters λ, n = (2 + 3λ)/λ, S _r, and H _o. Section 4.2.2 addresses uncertainty in parameters determined from the water retention curve measurements, in the context of snowmelt column test results.

3.5 Snowmelt column tests

We conducted snowmelt column tests in the field during both the 2018 and 2019 field sessions, to calibrate the SP measurements and determine values of zeta potential and n for application to Eqn (2) in evaluating in situ measurements. The snowmelt column consisted of a 0.6 m long schedule 40 PVC pipe with an inside diameter of 0.15 m, and with a beveled end to facilitate coring into the snowpack. We inserted the column vertically through the upper 0.6 m of the snowpack and subsequently excavated and removed it for testing. A sealed hot water reservoir placed at the top of the column induced controlled melting. An insulating cover and shading minimized annular column melting from ambient conditions. A tipping-bucket rain gage on the meteorological station measured the meltwater volume exiting the base of the column (Fig. 2). SP and TDR sensors in the lower half of the column measured ψ _m, expressed as SP electrical field strength (E), and S _w during meltwater percolation. We inserted 10 cm long SP electrodes through holes in the column spaced 0.2 m vertically with the reference (−) voltage electrode on top, and the (+) electrode placed 10 cm above the base of the column. A TDR waveguide was centered in between the electrodes. Section 3.6 describes the specifications of the SP and TDR sensors. We also used a Hanna Instruments model 98129 tester to measure σ _w and pH in the column effluent water. These data were not temperature compensated, and we performed daily in-field 2-point instrument calibrations.

Fig. 2. Snowmelt column test setup. The PVC column was fitted with SP and TDR sensors, and was placed over the tipping bucket rain gage on the meteorological station. A sealed hot water reservoir was placed at the top of the column to drive controlled melting. The column is shaded and wrapped in an insulating cover to minimize annular column melting from ambient conditions. Photo credit: Andrew Opila.

The column tests allowed an empirical comparison of observed meltwater effluent from the column to meltwater flux within the column modeled using Eqn (2). We applied the measured E and S _w data to Eqn (2) to develop a best fit of SP-modeled q _sp vs observed meltwater effluent from the column. This allowed the determination of an expected value of zeta potential and an evaluation of the sensitivity of the Brooks and Corey pore size parameter, n. We evaluated several values of n using Eqn (2) for each column test, in order to assess the sensitivity of this parameter to calculate transient meltwater flux. For each value of n evaluated, an expected value of zeta potential was then determined by calibrating measured flux to calculated flux, with a sum of residuals of zero over the test period.

The empirical approach adopted to evaluate the column test data involved applying a fixed value of zeta potential to solve Eqn (2) over the full duration of each column test, and subsequently to the in situ measurements. However, in reality zeta potential varies with changes in pH and σ _w (Revil and others, 1999; Kulessa and others, 2012). We measured these water quality parameters over time in meltwater eluted from the column tests. However, methods for the continuous sampling and measurement of these water quality parameters from pore water within the snowpack remain for future development. Therefore, without transient in situ measurement of pH and σ _w, the estimation of zeta potential is likely approximate. It is ultimately a limitation of this feasibility study that we are using zeta potential as a fitting parameter, and not directly measuring the electrochemical properties of the system. Therefore, the value of zeta potential determined may not be explicitly representative of the actual electrochemical properties of the system. This is an area for further research. Accordingly, we cannot resolve the potential errors associated with the estimation of zeta potential in this feasibility study, although Section 5.2 provides further discussion.

3.6 In situ snow sensors

We placed snow sensors horizontally into undisturbed snow through the vertical wall of a snow study pit for the measurement of SP and S _w. Sensors were pre-cooled to 0°C, and immediately after placement of the snow sensors, the snow study pit was backfilled with snow to maintain the 0°C isothermal system and prevent melt-out of the sensors. We inserted the SP electrodes to a distance of 1 m beyond the snow study pit sidewall. We placed TDR waveguides used for the measurement of S _w to a distance of ~0.2–0.3 m into undisturbed snow. This spatial alignment minimized the potential disruption of the flow field through undisturbed snow in the vicinity of the SP electrodes. However, since the SP and TDR data were not perfectly spatially coincident, spatial variability under meltwater percolation conditions could therefore contribute to errors when we combined the data into Eqn (2).

The SP electrodes consisted of 8-mm diameter round Ag–AgCl sintered biomedical electroencephalogram electrodes (BioMed Electrodes, model BME-8). The electrodes were embedded in epoxy, at a 45° angle, in the tip of a 9-mm diameter, 1.2 m long fiberglass wand (Figs 3, 4). We placed each SP electrode in the snow with the angled tip facing upward to maximize contact with downward percolating meltwater. In order to provide effective contact between the electrode and the snowpack, we placed a tension device at the snow study pit vertical wall, consisting of a bovine corkscrew trochar and elastic cord system (Fig. 4). This system placed mild forward pressure onto the electrode tip.

Fig. 3. 9-mm diameter SP electrodes fabricated using Ag–AgCl sintered biomedical electrodes embedded in epoxy at the tip of a fiberglass wand. White paint shows wear from snow insertion.

Fig. 4. Photo of in situ SP electrode wand, 1.2 m in length and with tip pressure device adapted from a bovine corkscrew trochar and elastic cord.

The placement of Ag–AgCl sintered electrodes directly into unsaturated porous media for SP measurements was described by Jougnot and Linde (2013). They found that electrode effects generating electrical potential would primarily be expected due to temperature variations, which are not present in the isothermal snowpack. Static testing of the electrodes in a no-flow snow meltwater bath resulted in stable measured electrical potential between the electrode pairs of <0.5 mV over a 48-h period.

In order to measure the electrical field strength (E) at specific locations, we wired the SP electrodes in vertically spaced pairs, with the high voltage (+) electrode of each pair oriented 20 cm physically below the reference (−) electrode. There was no isolated reference electrode. We used a bubble level to insert the electrodes horizontally into the snowpack to approximately maintain the desired spacing at the electrode tip. We placed a TDR waveguide centered vertically in between each electrode pair. We used the SP and TDR data to calculate the value of q _sp for each 5-min data logger time step using Eqn (2) and the parameters determined from the snowmelt column test.

Campbell Scientific Model CS655 TDR waveguides, which have two 120 mm long rods set at 32 mm spacing, provided measurements of S _w. Diaz and others (2017) previously used Campbell Scientific Model CS650 TDR waveguides to measure LWC in snow, which employs the same sensor technology and electronics. The TDR waveguide SDI-12 output signal reported the measurements of apparent dielectric permittivity (P), which we converted to S _w using an approximate function (Eqn (14)). Equation (14) is a purely empirical expression fitted to experimental data (Fig. 5) collected by placing a CS655 TDR waveguide into a 0.15 m long × 0.15 m diameter test cell. The calibration experiment involved measuring P while applying known changes in water saturation to the test cell under both drainage and imbibition conditions. We used Eqn (14) to process field TDR data collected during both the 2018 and 2019 field seasons:

(14)$$S_{\rm w} = \displaystyle{{( \;C_1^P /C_2-\;1\;) /100 + 0.0012} \over \Phi }$$

where C ₁ = 8, C ₂ = 9 and Φ is the porosity.

Fig. 5. Measurements of apparent dielectric permittivity produced using the CS655 TDR waveguides under controlled drainage and imbibition conditions of known water saturation (S _w). Empirical function Eqn (14) determined by visual fit to data.

3.7 Error analysis approach

This study involved determination of the Darcy velocity of unsaturated meltwater percolation flux, by solving Eqn (2) using measured SP electrical field strength (E) and S _w data. Thompson and others (2016) evaluated sources of error in the various parameters in Eqn (1). They noted that there is great difficulty in estimating the actual measurement uncertainty of individual variables, and therefore they considered hypothetical uncertainty ranges in these parameters within their error analysis. Due to the difficulty in determining individual measurement uncertainties, we have undertaken an empirical error analysis consisting of two parts. The percent difference (PD) was determined between the cumulative SP-modeled percolation flux over time (in mm w.e.) using Eqn 2, and (a) ablation measurements of snow loss using a snow stake, and (b) calculations of cumulative expected melt (M) using the surface energy-balance measurements. Measured precipitation was added to the above.

4 Results

These results present feasibility study data obtained during two separate ~1-week field sessions on the Juneau Icefield, in the summer of 2018 and 2019. As described above, the objectives of this study pertained to (a) development of field methods for the use of SP to measure unsaturated meltwater percolation Darcy velocity in snow, and (b) evaluation of the quantitative and qualitative potential of the SP method to measure short-term flow transients during diurnal melting cycles and their relationships to changes in weather and energy-balance drivers of snowmelt. Consistent with these objectives, the field data are short-term in nature. Field data presented from the 2018 field session illustrate transient conditions over a single diurnal cycle, and data from 2019 represent a period of approximately four diurnal cycles.

4.1 Water retention curves

Water retention curves for measured H _c vs S _w (Fig. 6a) for the two samples evaluated during 2018 resulted in air entry pressures (H _o) of 6 and 9 cm, respectively. Figure 6b shows transformed data plotted as H _c/H _o vs S _e with λ = 5 (n = 3.4) determined from a visual best fit to both datasets using the BC model (Eqn (6)) (Brooks and Corey, 1964) at S _r = 0.001. The water retention properties are a function of the snow grain size (Yamaguchi and others, 2010) as well as density, snow crystal morphology and other factors including the wettability of the ice crystals. Figure 7 shows a microphotograph typical of Matthes Glacier snow samples tested, depicting fully rounded melt forms with an average grain size of ~1.3 mm and a very narrow range of grain sizes with an absence of grains smaller than ~1.0 mm. Yamaguchi and others (2010) measured water retention curves for snow samples of various grain sizes and found λ = 6, H _o = 6 cm at d = 1.5 mm, and λ = 4, H _o = 10 cm at d = 1.1 mm. The values of H _o and λ are within the ranges found by Yamaguchi and others for the average grain size of 1.3 mm. The value of at S _r = 0.001 derived from our data is lower than the range of 0.02–0.07 previously reported for wet snow (Colbeck, 1974; Yamaguchi and others, 2010) and commonly used in models (Leroux and Pomeroy, 2019). Although the literature is sparse in this regard, it is possible that the low value of at S _r may reflect the very narrow range of the snow grain size and absence of small grain sizes in the summer snowpack on the icefield.

Fig. 6. Raw capillary head and saturation data (a) with dashed lines indicating observed air entry capillary head (Matthes Glacier: H _o = 6, Taku Glacier: H _o = 9), and transformed data (b) fitted to Eqns (4) and (5). Visual best fit of combined datasets to Eqn (6), at values of n = 3.4 (λ = 5) and S _r = 0.001.

Fig. 7. Microphotograph of fully rounded melt forms typical of the midsummer seasonal snowpack crystal morphology on the Juneau Icefield. Grid is 1 mm, and average grain diameter is estimated at 1.3 mm.

4.2 2018 season results

4.2.1 Snow pit profile

We excavated a snow study pit (Figs 8, 9) during the 2018 season to evaluate snow structure and properties and to emplace in situ sensors. As shown in Fig. 9, we placed paired SP electrodes with a 20 cm spacing, at depths of 25 and 45 cm below the snow surface. For the 2018 experiment, we selected this shallow depth to minimize effects of flow heterogeneity resulting from the presence of ice lenses and other snow structure at greater depth, and thereby provide a rapid and clean response in SP data to transient changes in surface melting. No snow structure was observable in the shallow study pit such as ice lenses or laminations consisting of variation in grain size. Snow morphology was consistent with the microphotograph of rounded melt forms shown in Fig. 7. Snow densities measured at an adjacent pit ranged from 550 to 575 kg m⁻³, compared to the average measurement of 560 kg m⁻³ at four mass-balance pits measured during 2018 within 4 km of the study site. These values are comparable to an average snow density of 540 kg m⁻³ measured in the historical JIRP mass-balance record including hundreds of samples over multiple decades (Pelto and others, 2013). Snow temperatures were 0°C throughout the profile. Although we placed four pairs of SP electrodes, as shown in Fig. 9, wiring damage occurred to three electrodes during transport. Therefore, in situ measurements presented below in Section 4.2.4 are from the single pair of functioning electrodes, E1–E2.

Fig. 8. Photo of 2018 snow pit after installation of SP and TDR instrumentation, before backfilling. Electrodes were placed at depth of 25 and 45 cm below snow surface. Inset depicts TDR waveguide with two 120 mm length rods.

Fig. 9. Vertical snow pit profile from 2018 field season showing instrumentation layout for in situ measurements. No ice lenses or other snowpack structure were observable. Electrode (E1–E8) pairs were used for SP measurements (SP1–SP4) at locations shown. TDR waveguides were placed at locations indicated. Red crosses represent three electrodes which failed due to wiring damage. Scale in cm.

4.2.2 Snowmelt column testing

Figures 10a, b show transient data collected during a snowmelt column test over a five and one half-hour time period from 12:00 to 17:30 on 25 July 2018. Figures 10c, d show results evaluated over a focused time period during the last 4-h of the test, used for analysis below. The addition of more heat to the column test at 15:30 (~1:45 duration in Figs 10c, d) illustrates the results of transient increased melting conditions.

Fig. 10. Field column test results from 2018 season. Data collected over the period of entire column test include (a) S _w and EC, and (b) E, cumulative and instantaneous (5 min average) observed effluent flux from column. Time period shaded in blue represents the duration over which data analysis was performed. Data for S _w and E over the duration of analysis (c) show the response to an additional pulse of column melting initiated at ~1.75 h duration. A comparison of observed column effluent flux and SP-modeled fluxes (q _sp) at three different values of n (d) shows the best fit of SP-modeled flux to the data, obtained at ζ = 1.73 × 10⁻⁷ and n = 1.7.

During the testing, meltwater flux observed exiting the base of the column reached a maximum rate of >1000 mm d⁻¹ and a minimum of 200 mm d⁻¹ at the end of testing (Fig. 10b). Approximately 100 mm w.e. of cumulative meltwater flux eluted from the column over the ~5-h period. Water quality data (Fig. 10a) showed that electrical conductivity (σ _w) decreased asymptotically from 3 × 10⁻⁵ to 2.5 × 10⁻⁶ S m⁻¹. pH was relatively stable, increasing slightly from 6.1 to 6.6 s.u. S _w varied between 0.07 and 0.13, increasing sharply in response to the additional pulse of melting applied during the test, reflecting wetting front propagation. SP electrical field strength (E) was a minimum of 50 mV m⁻¹ and a maximum of 650 mV m⁻¹, which for the 0.2 m electrode spacing represents raw SP signals of 10–130 mV.

Measured E and S _w data over the focused period shown in Figure 10c were applied to Eqn (2) resulting in SP-modeled meltwater percolation flux (q _sp) that can be compared to the observed values (Fig. 10d). Values of other parameters applied to Eqn (2) included σ _w measurements shown in Fig. 10a, and k determined using Eqn (7) and d = 1.3 mm. An S _r of 0.001 was used from Fig. 6. Three different values of the BC pore size distribution parameter (n = 3.4, n = 1.7 and n = 0.25) were used to evaluate sensitivity to this parameter. For each of these n values, an expected value of zeta potential (ζ = 2.36 × 10⁻⁹, 1.73 × 10⁻⁷ and 7.26 × 10⁻⁹, respectively) was then determined by calibrating measured to observed flux, with a sum of residuals of zero over the test period. Visual observation of the results (Fig. 10d) shows that the value of n = 3.4, determined from water retention curves, resulted in an overestimation of the transient response to the additional pulse of melting at a duration of 2 h in Fig. 10d. A value of n = 1.7 resulted in the closest fit to observed transient fluxes, and resulted in a value of ζ = 1.73 × 10⁻⁷ V to obtain a sum of residuals of zero. These values were applied to 2018 season in situ measurements described below in Section 4.2.4. The small, estimated value of ζ = 1.73 × 10⁻⁷ V may be consistent with the results of Kulessa and others (2012), who modeled expected values of ζ with changes in σ _w and pH during a natural snowmelt experiment, and found a reversal in sign of ζ at pH ~ 6.3 s.u. and σ _w ~ 2 × 10⁻⁶.

The empirical determination of n described above is justified, as opposed to direct application of n determined from water retention data, since the Brooks and Corey (1964) model for relative permeability (i.e. k _r = S _eⁿ), like other unsaturated flow relative permeability formulations, was developed based on empirical data from soil/mineral porous media systems and is untested for snow. Furthermore, the only non-linear term in Eqn (2), the Brooks–Corey term S _eⁿ, is very sensitive to the value of n, as illustrated in Fig. 10d.

4.2.3 Energy balance and weather data

Figures 11a, b show the time series of energy balance and weather data collected over a period from the evening of 31 July 2018 to the morning of 2 August 2018. Calculated turbulent energy fluxes Q _L and Q _H were generally small over the period, with maximum values of 39 and 24 W m⁻², respectively, although cumulatively they contributed to ~30% of total melting (Table 2). The air temperature was a minimum of 4°C and the dew point was above 0°C throughout the period, thus Q _H and Q _L were positive throughout the period. Observations of boot penetration indicated that snow surface refreezing occurred during the period from 21:00 on 31 July 2018 to 09:15 on 1 August 2018. This time period reflects a lag of ~1 h past the period during which net radiation (Q _N) and energy available for melt (Q _M) were below zero, resulting from overnight mostly clear sky conditions and diminished LW energy input to the snowpack.

Fig. 11. Results from 2018 field season, including (a) surface energy balance, (b) weather, (c) E and S _w and (d) meltwater fluxes. Q _R was a maximum of 1 W m⁻² and is not shown. The period highlighted in the yellow-shaded area represents the time period when total energy balance (Q _M) was negative, and pink represent rain. Surface refreezing and melting was observed via boot penetration observations at times indicated. Cumulative fluxes are shown that integrate all non-negative values for the period from 1 August 2018 08:00 to 2 August 2018 00:00. Calculated melt rate determined from Q _M (Eqn (13)), cumulative calculated melt and cumulative reduction in snow w.e. from ablation stake measurements are shown for comparison.

Table 2. Energy-balance summary for 24-h period of 1 August 2018

SW energy input under clear skies, with an average albedo of 0.65, was the main driver of daytime melting on 1 August 2018. After 17:00, and again at 18:20, LW energy input increased from the incremental arrival of increasingly foggy conditions. Increasing winds in this time frame also increased calculated values of turbulent energy fluxes Q _L and Q _H, which resulting in increased Q _M. Average and total surface energy fluxes for the 24-h period of 1 August 2018 (Table 2) show that 70% of the total energy available for melting was derived from net radiation (SW and LW), 16% was derived from latent heat flux and 14% from sensible heat flux. The average air temperature over this period was 4.7°C.

4.2.4 In situ measurements

In situ snowmelt measurements from 00:00 on 1 August 2018 to 06:30 on 2 August 2018 are shown in Fig. 11c, including S _w and SP electrical field strength, E (mV m⁻¹) calculated across electrode pair E1–E2. SP-modeled flux (q _sp) through the snowpack and melt rate calculated from surface energy balance (Eqn (13)) is shown in Fig. 11d. The SP-modeled q _sp values shown in Fig. 11d were determined using Eqn (2) and the following parameters: n = 1.7, S _r = 0.001, ɛ = 7.8 × 10⁻⁹ F m⁻¹, ζ = 1.73 × 10⁻⁷ V, σ _w = 2.5 × 10⁻⁶ S m⁻¹ and k = 1.65 × 10⁻⁹ m², based on measured values of d = 1.3 mm and ρ _s = 560 kg m⁻³. Also shown in Fig. 11d are cumulative fluxes over the time period following the initial observation of surface melt on 1 August 2018, including cumulative q _sp (non-negative values only), cumulative melt calculated from surface energy balance and cumulative snow loss observed at the ablation stake, all expressed in mm w.e. Section 5.1.1 presents discussion and interpretation of the data.

With regard to error analysis of the 2018 data, the cumulative SP-modeled flux of 18.5 mm w.e. was 5% less than the stake-measured ablation of 19.6 mm w.e., and was 23% less than the cumulative melt of 24.1 mm w.e. calculated from surface energy balance (Fig. 11d). The average error was −15%, compared to an average value from the ablation stake and calculated melt.

4.3 2019 season results

4.3.1 Snow pit profile

We excavated a snow study pit (Figs 12, F13) during the 2019 field session to evaluate snow structure and properties and to emplace in situ sensors. The base of the snow study pit was at a depth of 130 cm, and at this depth a water saturated layer ~10 mm in thickness was present overlying a continuous ice lens. Several smaller ice lenses were present within the snow profile, as shown in Figure 13. We measured snow density at four depths, ranging from 547 to 575 kg m⁻³, with an average value of 560 kg m⁻³. Snow grain size and crystal structure were consistent with the microphotograph of rounded melt forms shown in Fig. 7. Snow temperatures were 0°C throughout the profile. We placed SP electrodes in pairs spaced 20 cm, at depths of 40 and 60 cm below the snow surface, and at 80 and 100 cm below the surface. In situ SP measurements presented below in Section 4.3.4 are from the four pairs of electrodes, E1–E2, E3–E4, E5–E6 and E7–E8, designated as SP1, SP2, SP3 and SP4, respectively. We paired TDR sensors with SP electrodes as shown in Fig. 13.

Fig. 12. Photo of 2019 snow pit after installation of SP and TDR instrumentation, before backfilling.

Fig. 13. Snow pit profile from 2019 field season showing snow structure and instrumentation layout for in situ measurements. Profile is oriented vertically with dimensions shown in cm. Gray cross symbols represent density sample locations and results. Electrode (E1–E8) pairs were used for SP measurements (SP1–SP4) at locations shown. TDR waveguides were placed at locations indicated.

4.3.2 Snowmelt column testing

Figures 14a, b show transient data collected during the snowmelt column test over a ~2.5-h period from 10:30 to 13:00 on 25 July 2019. Figures 14c, d show results evaluated over a focused time period during ~1 h of the test, used for analysis below. The addition of more heat to the column test at 12:15 (~0:50 duration in Figs 14c, d) illustrates the response to transient increased melting conditions.

Fig. 14. Field column test results from 2019 season. Data collected over the time period of entire column test include (a) S _w and EC, and (b) E, and cumulative and instantaneous (5 min average) observed effluent flux from column. Time period shaded in blue represents the duration over which data analysis was performed. Data for S _w and E over the duration of analysis (c) show the response to an additional pulse of column melting initiated at ~50 min duration. A comparison of observed and SP-modeled fluxes (q _sp) at three different values of n (d) shows the best fit of SP-modeled flux to the data, obtained at ζ = 6.60 × 10⁻⁶ and n = 0.25.

During the testing, meltwater flux observed exiting the base of the column reached a maximum rate of 750 mm d⁻¹ and a minimum of 400 mm d⁻¹ at the end of testing (Fig. 14b). Approximately 50 mm w.e. of cumulative meltwater flux eluted from the column over the 2.5-h period. Water quality of the column effluent (Fig. 14a) shows that water electrical conductivity (σ _w) decreased asymptotically from 1.3 × 10⁻⁵ to 2 × 10⁻⁶ S m⁻¹. pH was relatively stable, increasing slightly from 6.3 to 6.5 s.u. S _w varied between 0.056 and 0.029, and did not increase in response to the additional pulse of melting applied during the test. SP electrical field strength (E) was a minimum of 55 mV m⁻¹ and a maximum of 372 mV m⁻¹.

Measured E and S _w data over the focused period shown in Fig. 14c were applied to Eqn (2) resulting in a comparison of SP-modeled meltwater percolation flux (q _sp) to observed values (Fig. 14d). Values of other parameters applied to Eqn (2) included σ _w measurements shown in Fig. 14a, and k determined using Eqn (7) and d = 1.3 mm. An S _r of 0.001 was used from Fig. 6. Three different values of the BC pore size distribution parameter (n = 3.4, n = 1.7 and n = 0.25) were used to evaluate sensitivity to this parameter. For each of these n values, zeta potential (ζ = 1.26 × 10⁻¹⁰, 4.35 × 10⁻⁸ and 6.60 × 10⁻⁶, respectively) was then determined by calibrating measured to observed flux, with a sum of residuals of zero over the test period. Visual observation of the results (Fig. 14d) shows that the value of n = 3.4, determined from water retention curves, does not provide the best fit to observed results. A value of n = 0.25 resulted in the closest fit to observed transient fluxes, which resulted in an expected value of ζ = 6.60 × 10⁻⁶ V to obtain a sum of residuals of zero. These values were subsequently applied to 2019 season in situ measurements described in Section 4.3.4.

4.3.3 Energy balance and weather data

Figures 15a, b show the time series of energy balance and weather data collected over a period from the evening of 23 July 2019 to the morning of 28 July 2019. Calculated turbulent energy fluxes Q _L and Q _H were generally small over the period, with maximum values of 13 and 20 W m⁻², respectively, although cumulatively they contributed to ~14% of total melting (Table 3). Q _R was a maximum of 6 W m⁻². A dew point below 0°C for a brief period in the early hours of 25 July 2019 resulted in a brief period of calculated evaporation, otherwise no negative Q _L occurred. Observation of boot penetration indicated snow surface refreezing during the overnight period spanning the early hours of both 24 July 2019 and 25 July 2019, reflecting net radiation (Q _N) and energy available for melt (Q _M) values below zero, resulting from clearing sky conditions and diminished LW energy input. Precipitation in the form of 1–3 mm of rain occurred daily over the period.

Fig. 15. Results from 2019 field season, including (a) surface energy balance, (b) weather, (c) S _w, (d) E and (e) meltwater fluxes. Time periods highlighted in pink represent rain, and yellow-shaded areas represent conditions where total energy balance (Q _M) was below zero, and snow surface was refrozen. Measured SP values are unadjusted from actual measurements. Dashed lines in (e) represent cumulative fluxes from q _sp and from calculated surface melt, over a 2-d period from 25 July 2019 00:00 to 27 July 2019 00:00. Calculated melt rate determined from Q _M (Eqn (13)), cumulative calculated melt + precip., and cumulative reduction in snow w.e. from ablation stake measurements are shown for comparison.

Table 3. Energy-balance summary for 72-h period of 24–26 July 2019

Figure 15a shows that SW energy input was the largest driver of daytime melting during each diurnal period. Skies were mostly variably cloudy, with an average albedo of 0.72 over the period. Hourly timescale variations in net radiation resulted from variously clear skies, and cloudy or foggy conditions. Average and total surface energy fluxes for the 72-h period of 24–26 July 2019 (Table 3) show that 86% of the total energy available for melting was from net radiation (SW and LW), 6% was from latent heat flux and 8% from sensible heat flux. A total energy input of 0.2% was from rain. The average air temperature over this period was 2.6°C.

4.3.4 In situ measurements

In situ snowmelt measurements from the evening of 23 July 2019 to the afternoon of 27 July 2019 including S _w and SP electrical field strength, E (mV m⁻¹) calculated across four electrode pairs are shown in Figs 15c, d. Figure 15e shows SP-modeled unsaturated meltwater flux (q _sp) through the snowpack, as well as calculated melt rate determined from Eqn (13). The q _sp values shown in Fig. 15e were determined using Eqn (2) and the following parameters: n = 0.25, S _r = 0.001, ɛ = 7.8 × 10⁻⁹ F m⁻¹, ζ = 6.60 × 10⁻⁶ V, σ _w = 2.0 × 10⁻⁶ S m⁻¹ and k = 1.65 × 10⁻⁹ m², based on measured values of d = 1.3 mm and ρ _s = 560 kg m⁻³. Section 5.1.2 presents a discussion and interpretation of the data.

The error analysis of the 2019 season data involved comparing the cumulative SP-modeled flux to cumulative melt calculated from surface energy balance, and cumulative snow loss observed at the ablation stake, all expressed in mm w.e (Fig. 15e). Table 4 shows the resulting total cumulative flux values and an average value for the four sensor locations over the 48-h period from 25 July 2019 00:00 to 27 July 2019 00:00. Also, in Table 4 are the PD values between the cumulative SP-modeled fluxes, the stake-measured ablation (49.2 mm w.e.), and cumulative melt calculated from surface energy balance (39.4 mm w.e.). Approximately 2 mm w.e. of precipitation occurred as rain over this period, which is included in the cumulative totals above. The average error in cumulative Darcy flux was +40% relative to the ablation stake and +52% relative to the calculated melt based on surface energy balance, with an average of +46% for the two methods.

Table 4. Error analysis of 2019 season calculated cumulative flux for 48-h period of 25–26 July 2019

PD, percent difference.

^a Calculated melt plus precipitation = 39.4 mm w.e.

^b Snow loss at ablation stake plus precipitation = 49.2 mm w.e.

5 Discussion and conclusions

5.1 Qualitative data interpretation

The data reported herein provide an opportunity to evaluate the dynamic relationships between surface energy-balance drivers for melt, and the subsequent observed transient unsaturated flow conditions within the snowpack. The primary objective of this research was method development for the measurement of unsaturated flow Darcy flux within the snowpack. However, ultimately the method may have utility in studying a variety of fundamental short-timescale transient flow processes. Although many of the transient variations in meltwater percolation data (Figs 11, 15) are seemingly random, at other times the meltwater percolation data obtained showed transient variations that appear to be qualitatively correlated with short-term changes in weather or surface energy balance. The following subsections provide a qualitative discussion of several potentially relevant transient conditions observed in the 2018 data and the 2019 data.

5.1.1 Qualitative evaluation of 2018 data

During the early hours of 1 August 2018, the mV-level SP data were negative in value (Fig. 11c), which may reflect upward flow in response to the observed overnight surface refreezing. The 2019 data discussed below and data from Clayton (2017) also included negative SP data in snow during overnight periods, under conditions of diurnal surface refreezing. Several prior authors have recognized the phenomena of upward unsaturated water flow toward a freezing front in soil in association with frost heaving (Taber, 1930; O'Neill and Miller, 1985; Iwata and others, 2010; Rempel, 2010). The magnitude of observed negative q _sp (~−50 mm d⁻¹) on 1 August 2018 was greater than the magnitude of negative calculated melt rate (~−10 mm d⁻¹), based on surface energy balance (Fig. 11d). It is unclear if this difference reflects measurement errors, unsaturated transient flow dynamics or some combination of the two.

Following development of positive surface energy balance (Q _M) at 08:00 on 1 August 2018 and subsequent observation of surface melting by boot penetration at 09:15, measured SP data ultimately transitioned into positive values (Fig. 11c). Maximum surface energy balance and calculated melt occurred at ~13:45 on 1 August 2018. The maximum S _w value of 0.09 (equivalent LWC = 0.035) at the measurement depth of 0.35 m, occurred at ~14:45, reflecting a 1-h lag behind maximum surface melting. Peak q _sp values of ~90 mm d⁻¹ occurred at 15:45, which lagged peak S _w by 1 h. Subsequent transient changes in q _sp and S _w occurred from ~17:00 to 20:00, which correlated qualitatively to variations in net radiation in that time frame, but with an ~30 min lag time to the observed flow response at depth. Note that the magnitude of fluctuations of q _sp in this time frame were larger than the magnitude of changes in calculated melt, which may represent unstable unsaturated flows, associated with preferential fingered flow during transient percolation events (Gerke and others, 2010).

The data from 1 August 2018 showed that large changes in water saturation were associated with the initial wetting front propagation, but that S _w was relatively insensitive to subsequent transient variations in meltwater flux later in the day. This is consistent with fundamental unsaturated flow physics, where changes in saturation represent changes in water storage within the porous media, while changes in flow are independent and may occur without change in storage.

5.1.2 Qualitative evaluation of 2019 data

During the overnight period into the early hours of 24 July 2019 (Fig. 15), the surface energy balance (Q _M) was negative, we observed surface refreezing, and three out of the four SP sensor data were also negative in value. As discussed above for the 2018 data, this may reflect upward flow in response to the observed overnight surface refreezing. However, during the 2019 study period, there were other times when negative SP values occurred while surface energy balance remained positive and we did not observe surface refreezing. These results may represent spatial variability in flow processes, and/or drift in SP measurements over the study period. We also cannot rule out the possibility that variations in water chemistry led to a reversal in the sign of zeta potential, as was inferred in snowmelt experiments by Kulessa and others (2012).

During the morning of 24 July 2019 (Fig. 15), the onset of positive surface energy balance and surface melting was followed by a rapid response in SP data, while trends in S _w were initially unchanged. Beginning at 05:20, 2.3 mm of rain fell over 1.5 h, at an average precipitation rate of 37 mm d⁻¹. Subsequently, S _w1 and S _w2 exhibited rapid increases, with the deeper sensor, S _w2, lagging S _w1 by 20 min, which likely reflects percolation of this precipitation into the snowpack. The less intense rain events on subsequent days resulted in more subdued responses in saturation and flux.

Measured S _w data during 2019 (Fig. 15c) reflected propagation of a wetting front each day in response to diurnal melting, with the shallower sensors responding before the deeper sensor at the corresponding location. The time differential in diurnal wetting front propagation across the 40 cm. vertical spacing ranged from 1 to 3 h indicating a wetting front propagation velocity of 0.13–0.4 m h⁻¹, which is comparable to the range of 0.2–0.5 m h⁻¹ observed by Clayton (2017) and 0.05–0.2 m h⁻¹ observed by Samimi and Marshall (2017). Similar to the 2018 data discussed above, the 2019 data also suggest that large changes in water saturation are associated with the initial wetting front propagation, but that S _w is relatively insensitive to subsequent variations in meltwater flux later in the day.

The values of S _w measured during the 2019 study period ranged from 0.008 (S _w4 on 25 July 2019) to 0.45 (S _w2 on 24 July 2019). Note that for the disparate data at these two monitoring locations, the values of S _w tend to scale with the paired SP electrical field strength values, E4 and E2. The ranges in these values are consistent with Eqn (2), where for the same Darcy flux value, a greater water saturation (S _w) results in a greater electrical field strength (E). Despite the large differences in S _w and E between these locations, the calculated cumulative meltwater percolation flux over a 2-d period at each location was similar (within ~ 28%), suggesting that the results are realistic.

5.2 Uncertainty and limitations

An important limitation of the exploratory methods described herein pertains to the multi-step process that involves in-field column testing as a means to empirically calibrate the in situ measurements. It would be ideal if we could identify accurate and physically representative values for each of the parameters in Eqn (2) at each point in time and space over a monitoring period of interest. We can measure SP electrical field strength and water saturation, but with the current state of knowledge and tools, we must estimate zeta potential, meltwater electrical conductivity, residual saturation, the Brooks–Corey parameter (n) and relative permeability in order to determine the Darcy velocity of meltwater percolation flux. This challenging reality led to the adoption of an empirical approach, where we estimated some of these values from results of a column test conducted in the field. However, the snowpack properties determined from the column tests may vary spatially and temporally in melting snowpacks, therefore the column tests may need to be repeated using multiple samples from different locations or times.

Although calibration of parameters in Eqn (2) to column tests in the field is in a sense empirically robust, it may belie the physical meaning of some of the variables. In particular, the value of zeta potential (ζ) was determined empirically from the column tests, using the value of the BC parameter n, as a qualitative fitting parameter. This simplifying assumption resulted in a reasonable approximate empirical fit of SP-modeled column test fluxes to observed data, but ultimately negated the explicit physical meaning of the derived value of ζ, which may not accurately represent the actual electrochemical properties of the system. Also, the calculation of meltwater percolation flux assumed that the value of ζ was constant over the measurement period, while previous research has shown that ζ is a sensitive function of pH and electrical conductivity of the meltwater, which will change over time (Revil and others, 1999; Kulessa and others, 2012). This problem may be less pronounced in glacier snowpacks that are thicker and have a more extended melt cycle and may therefore be closer to equilibrium chemistry by midsummer, as compared to rapidly melting thin spring snowpacks. The above identified limitations suggest that further evaluation of the temporal changes in meltwater chemistry and associated sensitivity of ζ in association with in situ measurements will be important in future research. This will require development of methods for continuous sampling and chemical analysis of pore water from within the snowpack.

Average errors in the measurement of Darcy velocity determined in this study were −15% in 2018 and +46% in 2019. In the preliminary pilot study, Clayton (2017) reported similar errors in SP-calculated Darcy flux ranging from +26 to −47% over 24-h measurement cycles, with an average error in calculated flux of +8% over a 4-d period. The results here included highly transient variability in SP signals, at times appearing capricious in nature. Theoretically, contributions to observed SP signals include not only the streaming potential from fluid flow, but also electrochemical potentials from gradients in ionic solutes (Revil and others, 2007). Since the glacier system is isothermal at 0°C, we do not expect thermoelectric potentials, but electrode effects cannot be completely ruled out (Jougnot and Linde, 2013). Spatial variability of meltwater percolation in snowpacks and development of preferential flow paths and temporally variable flow magnitudes have been well documented (Colbeck and others, 1979; Marsh and Woo, 1985; Williams and others, 2010; Avanzi and others, 2016; Webb and others, 2018). It is therefore unclear how much of the observed error relates to measurement errors such as undetermined electrochemical potentials, and how much relates to the inherent variability of point measurements within a spatially and temporally variable flow field. Furthermore, there is uncertainty in the 3-D spatial patterns of SP electrical field strength in the vicinity of discrete preferential flow paths in melting snowpacks, and whether the correlation length scales for flow and electrical field strength are equal. Measurement of SP electrical field strength at multiple scales in a 2-D gridded vertical profile could potentially shed light on the spatial variability of unsaturated flow processes in a heterogeneous snowpack, either qualitatively or quantitatively. Also, in future research the addition of a spatially distant reference electrode may help to understand whether larger scale effects from telluric currents or subglacial drainage contribute to the measured SP-field strength.

Although the methods used herein minimized disturbance of the natural snowpack structure, the emplacement of measurement sensors in the snowpack is nonetheless invasive. Avanzi and others (2014) described snow sensor meltout and development of an air gap around surface-emplaced sensors with the insertion point exposed to solar radiation. In this study, we placed the sensors laterally through the vertical wall of a snow profile pit which was then rapidly backfilling the pit with 0°C snow, which minimizes the potential for sensor meltout. Nonetheless, the integrity of sensor contact with snow is transient and may degrade after a period of days or weeks. This limits the appropriate monitoring period for this method to short-term periods on the order of days to weeks, although alternative methods may potentially extend the viable monitoring period for SP measurements in snow.

5.3 Conclusions

The exploratory research reported herein involved field measurement of the Darcy velocity, or flux, during unsaturated flow within melting snowpacks at the decimeter-scale. The results confirmed the feasibility to measure mV-level electrical SP within melting snowpacks to determine downward percolation flux in a vertical profile. The data showed a strong qualitative relationship of transient variations in flux to diurnal and shorter timescale snowmelt dynamics. With further development the method might be appropriate to more long-term snowmelt monitoring.

The calculated fluxes were comparable to actual fluxes, with average errors in cumulative fluxes ranging from −15 to +46% compared to the average of expected melt calculated from surface energy-balance and stake measurements of snow ablation. Errors may have reflected a combination of uncertainty in parameter estimation, the presence of undetermined electrochemical potentials or electrode effects and inherent variability in point measurements within a spatially and temporally variable flow field. It is unclear how much of the error observed relates to measurement error vs flow heterogeneity. There was no attempt at mathematical modeling of the results obtained, but this may be useful in future research.

Identification of the dominant sources of error and improvement in the accuracy of the method will require further research. The method determined zeta potential as an empirical fitting parameter derived from snowmelt column tests, and the values determined may not be explicitly physically meaningful. Furthermore, the simplifying assumption that meltwater chemistry and zeta potential were constant over time within the snowpack is a limitation necessitated by a lack of methods for continuous measurement of pH and electrical conductivity of meltwater within the snowpack. Research to address this gap and characterize transient changes in zeta potential will improve the technique. Water retention curves measured from snowpack samples indicated a value of the Brooks and Corey (1964) pore size distribution parameter (n = 3.4) that was very similar to that identified in previous snowmelt studies. However, since relative permeability functions such as the Brooks and Corey k _r = S _eⁿ were determined empirically from soil and mineral systems, we cannot be certain of the applicability of n derived from water retention curves to the above k _r function for snow. The column test results showed that alternative values of n were more descriptive of transient changes in measured fluxes, suggesting a need for further research in this area.

Surface energy-balance measurements and calculations over variable summer weather conditions at the Matthes-Llewellyn glacier divide on the Juneau Icefield in 2018–19 indicated that net radiation (SW and LW) accounted for 70–86% of energy available for melting. Latent heat flux accounted for 6–16% of melting and sensible heat flux accounted for 8–14% of melting. Peak diurnal meltwater flux lagged peak surface melting, as expected due to the time required for flow from the surface to the depth of measurement. We also found the diurnal peak SP-measured flux was later than arrival of the unsaturated flow wetting front. Large changes in S _w accompanied initial arrival of the diurnal wetting front, but S _w was relatively insensitive to subsequent variations in meltwater flux later in the day. Therefore, the data suggest that changes in S _w are not a reasonable surrogate for flux.

We observed upward unsaturated flow of water within the snowpack at depths up to 1 m, at times of surface refreezing in response to negative surface energy balance. Although previous literature has not reported this phenomenon in snow, an analogous process involving upward flow toward a freezing front during frost heaving in soils has been previously reported. This process may be important to fundamental meltwater percolation dynamics, mathematical modeling of unsaturated flow in snow, and considerations of snowmelt runoff and glacier mass balance.

This was a feasibility study, and the methods presented require further development. Nonetheless, the use of SP in the field to measure the unsaturated flux of percolating meltwater within snowpacks may be a relevant tool for the study of several important topics. Field measurements can provide valuable ground-truth for mathematical modeling and can elucidate fundamental processes that are poorly understood. Measurement of snowmelt percolation flux can support studies of snowmelt–groundwater and snowmelt–runoff interactions, as well as improve forecasting of wet snow avalanches. The technique also has potential to supplement glacier mass-balance studies by quantitative measurement of fluxes associated with processes such as meltwater percolation, internal accumulation and firn aquifer formation.