Introduction

Polluted snow generally contains a great many impurities which may arise from natural sources, deposition of locally created contaminants, or from long-range sources of pollutants transported within the upper atmosphere. During the winter any contaminants contained within the precipitation, or deposited later onto the snow surface, are accumulated over a period of time within the snow-fields and transmitted to water channels only in the spring when the snow melts. The impact of collected contaminants on the local environment is intensified by this melting process since pollutants held initially within a snow-pack are not released at a uniform concentration corresponding to their bulk concentration, but are released with significantly higher concentrations within the first melt water. A typical melt curve obtained by Reference Wooding and Morel-SeytouxJohannessen and others (1977) is reproduced as Figure 1 and shows the concentrations of a number of chemical impurities expressed relative to their bulk concentrations. The interaction of the snow-melt process with any pollutant concentration initially held within the snow-pack is therefore critically important in evaluating the impact of chemical contaminants accumulated over the winter period on the environment.

Fig. 1. Relative concentrations experimentally measured in the run-off from a polluted snow (from Reference CrankJohannessen and others, 1975). Vertical bars denote the extent of pollutant concentrations in the first fractions of melt water.

General studies of the movement of impurities within a liquid flow through a porous medium have been extensively orientated towards investigating the redistribution of fertilizers or transport of a tracer within a soil (Reference Colbeck and DavidsonKirda and others, 1974). Reviews of the theories of water infiltration through a soil have been given by Reference PhilipPhilip (1970), Wooding and Reference Wooding and Morel-SeytouxMorel-Seytoux (1976), and Reference ParlangeParlange (1980), whilst a review of the spreading encountered by a tracer within a flow through a porous medium (hydrodynamic dispersion) was given by Reference BearBear (1972). Percolation of melt water through a snow-pack at the start of the melting season can be described by using the theory established by Reference Colbeck and DavidsonColbeck (1972 Reference Colbeck1974) and Reference Colbeck and DavidsonColbeck and Davidson (1973). A dominant physical feature associated with the first influx of percolating water into an unsaturated porous medium is the appearance of a moving saturation wave-front, across which the level of water saturation changes abruptly. Application of the theory for tracer movement has been made by Reference ColbeckColbeck (1977) to obtain an elementary solution by neglecting any hydrodynamic dispersion and consequently, all the tracer is converted within the interstitial pore water movement. Reference BreslerBressler (1973) included dispersion effects in computed solutions for the displacement of a chemical solute from a soil by percolating water, but this used a numerical technique which only approximates the motion of the saturation wavefront formed.

A theoretical description involving all the diverse physical aspects which could affect the release of a chemical pollutant within the run-off from a melting snow-pack is a formidable task. Consequently, to enable a more modest and tractable calculation for the initial pollutant concentrations to be made, a model is presented which is based only on the dominant physical features affecting the redistribution of pollutant. The simplest, yet physically realizable, snow-melt situation arising from a uniform surface melting is evaluated in detail so that theoretical predictions can be readily obtained and compared with known experimental results. This theory is readily extendible to more complex flow situations arising from a differential heating or diurnal melting which results in the formation of more than one saturation wave-front, but such solutions would require more extensive computation.

Snow-Melt Model

Impurities within a snow-pack may either be associated with the formation of snow-flakes in the atmosphere or else from dry deposition onto the snow surface. A recent report by Reference SeipSeip (1978) suggests that over most parts of Norway the dry deposition of pollutant is relatively small and that the larger part of impurities trapped in seasonal snow-packs is associated with formation of the ice crystals. The type of ice formed for snow under normal conditions does not readily include impurities bound within the ice crystals and any such concentration may subsequently be reduced in the process of melt-metamorphism. Thus pollutants held in a seasonal snow-pack at the start of the melting season are most likely to be contained with the interstitial water content existing at the grain boundaries and to have an initial concentration which is approximately uniform with depth.

As noted earlier, the principal flow feature which emerges from the unsaturated flow of melt water through a snow-pack is the development of a saturation wave-front which either leads the percolating water, as in the case of the first influx of melt water, or else arises from within subsequent diurnal fluxes of melt water. Initially, any pollutant remains held within the pores of the snow-pack until released by the arrival of the first wavefront, and afterwards the pollutant becomes mixed within the melt water. Movement of the wavefront is well known to be governed by the condition of continuity of liquid mass across the front. A similar argument of continuity of solute across the front must also hold, which must then relate the levels of pollutant concentrations on either side of the front. Once contained within the percolating liquid any pollutant is convected with the melt water but is also subjected to two major dispersive effects modifying the convected pollutant concentrations. Chemical properties of the pollutant give rise to molecular diffusion whilst the presence of the pore structure gives rise to a mechanical dispersion of the pollutant. Analysis of these phenomena are given by Reference SlatterySlattery (1972) and Reference BreslerBear (1972) who show that the total dispersive effect on the pollutant concentration can be readily described by a modified form of Fick′s law. Owing to the mechanical dispersion of pollutant originating from the dynamic interaction of the percolating liquid with the ice matrix, the coefficient of dispersion will depend strongly upon the interstitial pore velocity. At high melt rates the effect of the mechanical dispersion may be several orders of magnitude greater than molecular diffusion, in which case the elution of pollutant is expected to depend only slightly on the chemical species, and this is in part observed in the typical melt curves displayed in Figure 1.

At the start of the melting season any application of heat to the surface of the snow-pack causes melting confined to a surface layer of snow, producing melt water which then percolates into the remainder of the snow-pack. This surface layer requires a treatment which must include temperature effects such as given by Reference Obled and RosseObled and Rosse (1977) to determine details of melt water released. Below this surface layer an isothermal percolation theory is applicable if we assume that the bulk of the snow-pack is mature.

Governing Equations

A one-dimensional theory for water percolation through a homogeneous snow-pack has been reviewed by Reference ColbeckColbeck (1978). This theory uses Darcian flow concepts applied to a gravity-dominated flow together with an experimentally determined relationship that links the flux of percolating water to the local value of water saturation held within the ice pores. Under typical fair-weather melting conditions the volume flux of melt water can be described by a relationship of the form

(1)

where S is the effective water saturation as introduced by Reference Colbeck and DavidsonColbeck and Davidson (1973) who also give typical values for tne constant k and the exponent n. Combining the experimental result given by Equation((1) with the usual equation for mass conservation yields the equation

(2)

to determine S, where z ^⋆ measures the depth into the snow-pack, t^⋆ the elapsed time, Φ the porosity, and S_i the irreducible liquid content of the snow-pack. Continuous solutions of Equation (2) are readily obtainable, but are severely limited in use due to the formation of saturation wave-fronts forming within the flow which must be given a separate treatment. If at an initial instant an abrupt volume flux of liquid is applied to the top of an otherwise undisturbed snow-pack then a wave-front will develop immediately to separate the imposed saturation level from the initial level. The speed V^⋆ of a wave-front and the percolation velocity u^⋆ at the wave-front are determined from the requirement that across the moving discontinuity in saturation levels the mass of migrating melt water must be conserved. Applying this criterion gives that

(3)

where β = S _i/(1-S _i) and S ₊, S ₋ are the saturation values directly behind and preceding the wavefront.

In general the addition of a second liquid phase within the calculation considerably increases the difficulty in describing the flow processes. Simplification is gained however by noting that most liquid pollutants present within a snow-pack consist of a weak solute held within the melt water. Defining the concentration of a pollutant species as c (volume of pollutant per unit volume of melt liquid) the volume flux of solute J comprises a flux Qc, arising from convection with the solvent, and modified by the action of hydrodynamic dispersion. In the analogous problem of saturated flow through soil, the flux due to hydrodynamic dispersion is well established (see e.g. Reference BearBear, 1972) and leads to an equation for the solute flux of the form

(4)

where D^⋆ is the dispersion coefficient. For an unsaturated flow, Equation (4) can be modified by inclusion of the local saturation level whilst the flux of melt watery Q is given by Equation (1). Expression of a solute mass balance then yields the equation governing solute concentrations as

(5)

If the flow is taken as saturated then Equation (5) reduces to the familiar dispersion equation, whilst if the saturation is held constant behind the wave-front then Equation (5) is equivalent to the equation used by equation Reference BreslerBressler (1973). Solution of the solute concentration equation is coupled to the liquid conservation Equation (2), which therefore must be solved either beforehand or in conjunction with Equation(5).

Across any saturation wave-front a solute mass balance must be achieved simultaneously with that of the liquid. Denoting the pollutant concentration directly behind the front as c ₊ and the concentration preceding the front as c ₋, the balance condition which must be satisfied is given by

(6)

Within Equation (5) the quantities D₊ ^⋆ and u ₊ ^⋆ refer to the values of the dispersion coefficient and interstitial liquid velocity compatible with a saturation level S₊ and D₋ ^⋆ and u₋ ^⋆ are similarly defined from a saturation level S ₋. Development of saturation wave-fronts within the percolating liquid flow separates regions of continuous changes in liquid saturation. These regions also become the regions of continuous variations in solute concentrations as governed by the partial differential equation (5). The solute balance condition (6) thus provides an internal boundary condition to link regions of continuous solute variation across a wave-front.

Noting that the effective saturation S and the solute concentration c are both non-dimensional quantities, it is convenient for the purposes of a general study to introduce further non-dimensional variables. Denoting the total depth of the isothermal core of the snow-pack as L ^⋆, then a non-dimensional depth and time are given by z and t defined by

(7)

where ρ_ℓ and are μ_ℓ the density and viscosity of the melt water. Other corresponding non-dimensional quantities must then he given by

(8)

where V, u, and D denote the non-dimensional wave-front speed, percolation velocity, and dispersion coefficient. Such a formulation allows a general study of the resulting pollutant redistribution without recourse to the intricate details of a specific experiment which may only be found by inference.

Writing the governing equation in terms of the non-dimensional quantities, the liquid continuity equation becomes

(9)

and the solute continuity equation is then

(10)

The region occupied by the isothermal snow-pack is taken to have an upper surface specified by z = 0 and the base of the snow-pack is therefore given as z = 1. In non-dimensional form the wave-front speed is given by

(11)

whilst the solute balance condition which must be applied across any wave-front becomes

(12)

Uniform Melt-Water Percolation

The simplest experimental situation is obtained by subjecting the surface of an otherwise undisturbed snow-pack to a constant heat flux. Outside a surface layer of snow, the effect of the applied heating is to produce a uniform percolation of melt water headed by a saturation wave-front. All the pollutant is taken to be held initially within the irreducible water content at a concentration c _i which will be signifi­cantly higher than the more easily measured bulk pollutant concentration c₀ (volume of pollutant per unit volume of snow-pack). The relation between these two quantities is straightforwardly found as

(13)

where ρ_s is the bulk density of the snow-pack and ρ_ℓ the melt-water density. This relative concentration factor would be expected in the run-off if the melt water were produced from any relaxation of the capillary retention force or if the effect of percolating melt water is to simply displace the irreducible water content. The exceptional case with no interstitial mixing corresponds to the convection theory of Reference ColbeckColbeck (1977) which, although it has the advantage of analytic results, has only a limited range of applicability.

Taking the interior of the initially undisturbed snow-pack as possessing a constant porosity and irreducible saturation level, then the initial effective saturation level is zero and the pollutant is held at a uniform concentration c _i within the initial water content. Correspondingly these give constant values to flow quantities preceding the first saturation wave-front which penetrates into the snow-pack, and, due to uniform melt-water release, the effective saturation level S behind the wave-front will remain constant. Figure 2 shows schematically the liquid flow problem and its solution. The saturation wave-front moves with a constant speed V = Sⁿ⁻¹ and separates the undisturbed region of the snow-pack from the region containing melt water percolating at a uniform rate u = S ⁿ /(S+β). No melt water will be released at the base of the snow-pack until the time of arrival of the wave-front given by t = V⁻¹. After this time, melt water will be released from the snow-pack at a uniform rate. At the base of the snow-pack, accumulation of melt water may form a layer of saturated snow which will drain by a lateral flow along the snow-soil interface, or else percolate through the underlying soil. In this study a freely draining base is assumed so as to model the recorded experimental studies.

Fig. 2. Non-dimensional flow description corresponding to an imposed uniform saturation S_ℓ at the snow-pack surface.

Once the basic liquid flow has been established, evaluation of the resulting pollutant concentrations can be considered independently. With conditions of uniform melting imposed, the resulting constant saturation level behind the moving wave-front ensures a constant value for the dispersion and percolation velocity. Redistribution of the pollutant behind the moving saturation front is given by Equation (10) which must be solved within a varying spatial region given by 0<z<Vt. For times before the saturation front reaches the base of the snow-pack, pollutant is transferred from a captive state, held within the irreducible water content, to within a higher saturation level and convected towards the base of the snow-pack. Across the saturation front the solute concentration behind the front is related to the initial solute concentration by the balance of solute condition

(14)

on evaluating the general expression (12). After the saturation front has reached the base of the snow-pack all the pollutant has been released by the percolating melt water and the solute balance condition (14) must be replaced by a condition specified in accordance with the physical melt-water run-off geometry. The subsequent pollutant redistribution is given by the solution of the solute conservation equation (10) within a fixed region of calculation. Assuming the base of the snow-pack to be freely draining results in all the solute carried by the percolating liquid near the base being simply convected within the liquid phase. Consequently the condition which must be imposed on the solute equation at the base of the snowpack corresponds to a continuity with solutions in which D = 0 in z > 1, i.e. that

(15)

At the surface of the snow-pack we assume that no pollutant is introduced within the volume flux of liquid resulting from melting. It is inconsistent with the model to impose a condition of zero pollutant concentration at the top surface since any pure melt water introduced into the snow-pack would immediately become mixed with the initial solute concentration. Hence the condition which must be applied at the uppermost surface is that no flux of solute is introduced, i.e.

(16)

The complete theoretical formulation governing the variation of pollutant concentration is summarized by Figure 3.

Fig.3. Space-time diagram showing the mathematical formulation to determine the solute concentration c(z, t).

Calculation of Pollutant Concentrations

Even for the simplest problem of uniform melting, analytic solutions describing the resultant pollutant concentrations are not available and so computed solutions are presented. Numerical results are discussed in the next section, but certain aspects of the flow can be best understood from an analytic study of the problem in which only slight surface heating is applied. At very low melt rates the effect of mechanical dispersion is small and provided molecular diffusion effects are also small then any hydrodynamic dispersion might be negligible. Under these circumstances the solute concentration equation reduces to the first-order equation

(17)

which can be readily solved using the method of characteristics. Characteristic directions for the solute concentrations are given by

(18)

Behind the wave-front the percolation velocity is constant and so the paths of characteristics within a space-time diagram are straight lines. Since the parameter β is positive, all the slopes of the characteristics are less than the slope representing the path of the wave-front. Figure 4 shows the pattern of characteristics behind the wave-front which consist of two distinct sets, one set of characteristics arises from the surface of the show-pack whilst the other set originates from the wave-front. Following from the characteristic relation given by Equation (18), along any characteristic the value of the solute concentration is constant and equal to its value, either at the snow-pack surface, or immediately behind the wave-front. The dividing characteristic is given by z = ut and corresponds to the extent of in­filtration of melt water into the draining snow-pack which arises from the snow-pack surface. Without the inclusion of dispersion effects, the surface condition (16) reduces to simply c = 0 which determines the value of the solute concentration associated with all characteristics originating from the snow-pack surface. At the wave-front the solute balance condition (14) yields that c = c _i everywhere directly behind the front, and subsequently predicts that c(z,t) = c _i everywhere within the region of the snow-pack given by ut<z⩽Vt. Interpreted physically, the effect of the imposed flux of pure water is to displace all the initial solute and water mixture, which is initially held by capillary retention, ahead of the pure water. At the base of the snow-pack all the solute would be released with a time t = u ⁻¹−v⁻¹ = βS ⁻ⁿ and at a relative concentration factor c _i/c ₀ given explicitly by Equation (13). For a typical snow-melt situation this corresponds to all the solute being released within approximately the first 5–10% of snow melted and at relative concentration factors of 10–20, which is clearly in disagreement with the experimental findings. The asymptotic behaviour for the special case of large values of the dispersion is given by

Integration gives immediately that

and hence

(19)

where the two parameters α(t) and γ(t) are independent of space variation and are determined from the two boundary conditions specified by Equations (14) and (16). Evaluating these variables gives that

which in the limit of large values of D gives a prediction of the maximum decrease in relative pollutant concentration attainable from dispersion effects as

(20)

Detailed predictions for the development of the sol­ute concentrations within a melting snow-pack when dispersion effects cannot be neglected are only available from a numerical study. Numerical tech­niques for solving the solute concentration equation are well established (see, e.g. Reference Colbeck and DavidsonCrank, 1975) however, in the snow-melt situation a complication arises in that at first a varying region of integration must be used. A solute mass-balance relation provides a boundary condition for the solute equation which must be applied on the moving boundary of the region of integration. As soon as the wave-front reaches the base of the snow-pack the region of integration re­mains fixed and the freely draining bottom condition given by Equation (15) must be applied. Initially the region of calculation behind the wave-front has zero width and thereafter expands to cover the whole snow-pack. During the initial stages of development there is difficulty in applying a numerical integration routine and so initially a local series expansion was used.

Fig. 4. Solution for the concentration c(z, t) on using the method of characteristics in the case of negligible dispersion (convection theory).

Results

The experimental results shown in fiqure 1 of Reference Johannessen, Johannessen, Dale, Gjessing, Henriksen and WrightJohannessen and others (1977), and the experimental results of other authors, do not give sufficient details of the bulk physical properties associated with the snow samples used to allow a detailed evalu­ation of the theory presented. Theoretical results are therefore displayed covering a range of values for the principal snow properties and flow velocities which could be considered as typically measurable during the early spring melting period. Reference LerameláLemmelá (1973) reports that the irreducible residual water content of a mature snow-pack is typically within the range 2–3% per unit volume, whilst the melt-water content of a seasonal snow-pack reaches a maximum of approxi­mately 5% per unit volume. In laboratory experiments these values may be significantly higher, due to the possibility of a more Intense and sustained surface heating than can be applied naturally during the early spring melting.

The dependence of the coefficient of hydrodynamic dispersion upon the interstitial melt-water velocity is complicated by the dispersion phenomena arising from two different physical mechanisms. Reference ColbeckColbeck (1977) reports some unpublished findings to suggest that for unsaturated flow through a snow-pack, D ^⋆ = u ^⋆ d ^⋆ where d ^⋆ is a typical pore size. In terms of the non-dimensional formulation this gives a value

(21)

indicating that the coefficient of dispersion will in general depend upon the level of melt-water saturation and on the dimensions of the snow-pack.

Detailed computer solutions are presented in Figure 5 to determine the development of the pollutant concentration profiles held within the percolating melt liquid behind the advancing saturation wave-front. The calculated pollutant concentrations are shown relative to the initial interstitial pollutant concen­tration c _i but can be readily compared to the much more easily determined bulk pollutant concentra­tions on using Equation (13). In the limiting case of no dispersion we obtain a convection theory in which no mixing of the imposed flux of melt water and the irreducible water and pollutant content occurs. Consequently the applied melt water, has a well-defined leading edge and all the interstitial mixture must move ahead of this melt water. No dilu­tion of the polluted water content occurs and so the initial relative concentration remains unaltered by the percolation process and is collected between the wave-front and the leading edge of the imposed melt-water flux. With a small coefficient of dispersion specified, mixing takes place mainly at the mutual interface as formed by the convection, but the concentration gradients behind the wave-front are not sufficiently large to reduce the solute concen­tration levels immediately behind the wave-front. As the dispersion parameter is increased, denoting a stronger mixing process, the dispersion mechanism becomes efficient enough to produce a concentration gradient behind the wave-front and consequently, as a result of the solute balance condition (14), this must influence the value of the pollutant concen­tration found immediately behind the wave-front. In this case more of the pollutant is held by the porous structure of the snow-pack, and diluted by the pure melt water rather than simply converted withthe average interstitial pore velocity. Further increase in the value of the dispersion parameter produces a decrease in the pollutant concentrations found im­mediately behind the wave-front with correspondingly more pollutant mixing within the pore structure. Figure 5 shows the profile development with the satur­ation parameters taken as S = 0.05 and β = 0.05, which approximately corresponds to a snow-pack in which the irreducible saturation level is 1.5% volume, and with an imposed melt-water content of approximately 3% per unit volume.

Fig.5. Theoretical development of concentration pro­files behind the saturation wave-front for different values of the dispersion parameter.

D = 0; b = 10⁻⁵; D = 10⁻⁴; D = 10⁻³.

Saturation parameters are taken as S = 0. 05 and β = 0. 05.

The value of the pollutant concentration directly behind the wave-front as it reaches the base of the snow-pack is of particular importance, since this determines the level of pollutant concentration with the first fraction of run-off. Figure 6 records the variation of this first concentration of pollutant with changes in the coefficient of hydrodynamic dis­persion. Results are displayed for a variety of im­posed liquid saturation levels between S = 0.01 and S = 0.1, with relative concentration factors shown against values of the dispersion parameter D plotted on a logarithmic scale. With a value of D below 10⁻⁵, in both sets of results shown, the relative pollutant concentrations remain unaltered from their interstitial values, whereas for values of D above 0.1, the concentration levels as predicted by the asymptotic result (20) are reached. On increasing the value of the saturation, and therefore the flow rate, a larger dispersive mechanism is required to provide any marked change in the relative concentra­tion, but ultimately a stronger dispersive mechanism produces a much reduced relative concentration. In the case of low flow rates, with S < 0.01, a very small amount of dispersion will affect the first fraction of pollutant but does not change the concen­tration significantly for any value of the dispersion; in the limit of decreasingly small flow rates the convection theory result of no initial pollutant concentration reduction from its interstitial value is recovered.

Fig. 6. Variation of the relative pollutant concen­tration c/c_i within the first melt liquid from a snow-pack on changing the dispersion parameter D. The initial irreducible saturation is given by β = 0.05.

In Figure 7 the calculations are extended to pro­vide a prediction for the distribution of pollutant concentrations present in the run-off. This is accomplished by equating the mass of liquid calcula­ted as flowing through the base of the snow-pack with the amount of the snow-pack which would have been re­quired to melt in the physical situation. This simu­lation is obviously not accurate once significant melting has taken place since the snow-pack character­istics will have changed radically, but for reasons outlined earlier such a procedure is likely to be fairly accurate near the start of the melting period. Superimposed on the melt curve are the vertical bars shown in Figure 1 which denote the extent of the measured pollutant concentration distribution in the melt run-off. Provided that sufficiently high dispersion is produced by the melt water, Figure 7 clearly shows that the experimentally observed efflux of pollutant from a snow-pack can be modelled.

Fig.7. Comparison of theoretically calculated pollut­ant concentrations with the experimental results of Figure 1. Snow parameters are taken as S = 0.05, β = 0.05 and c_i/c₀ = 10.18.

Results of the snow-melt model appear to be in agreement with the available experimental studies, however two aspects require further comment. First the percolation of melt water is not a homogeneous process; melt water often follows preferred channels through the snow-pack rather than percolating down­wards as a uniform front of water. Such inhomogenei­ties will cause pollutant and melt water to penetrate the snow-pack more quickly within the preferred channels, resulting in an earlier and less abrupt start to the melt-water run-off. The second aspect concerns the value for the coefficient of hydro­dynamic dispersion. Theories for hydrodynamic dis­persion are well proven for the saturated flow of liquids through porous media, but many results for un­saturated flow must be obtained by direct inference from the corresponding results for the saturated flow. Correspondingly the prediction for the dispersion co­efficient given by Equation (21) may need modifica­tion when used for the highly unsaturated flow con­ditions typical within a melting snow-pack.