A dry snow pack model

cold regions science and technology ELSEVIER

Cold Regions Science and Technology, 22 ( 1994 ) 135-148

A dry snow pack model J . M . N . T . Gray, L.W. M o r l a n d School of Mathematics, University of East Anglia, Norwich NR 4 7TJ, UK (Received 1 March 1993; accepted after revision 29 May 1993 )

Abstract A dry snow pack is viewed as an interacting mixture of rigid ice grains and air occupying the pore space in the matrix. The three-dimensional mass, m o m e n t u m and energy balances are formulated within the framework of interacting continua theory. The assumption of a c o m m o n ice and air temperature is made to reflect the very slow relative motion between the air and ice. Focusing on a one-dimensional vertical snow pack, five constitutive postulates and two external prescriptions are made to obtain a closed system of partial differential equations. We consider the response o f the snow pack to surface forcing by annual variations of heat flux and pressure. Nondimensional analysis o f the differential equations and surface conditions determines the major balances and leads to a reduced model. The equations are solved numerically for the temperature, air velocity, density and pressure through the snow pack resulting from both thermal forcing and pressure forcing at the surface.

I. Introduction

Dry snow packs are abundant on the Antarctic continent where sub-zero temperatures are maintained throughout the year, preventing the formation of any significant water content. There have been a number of experiments (e.g. MacDowall, 1960) that have measured the temperature distribution with depth over a complete annual cycle, and large numbers of measurements are made regularly of the ten meter snow temperature used as a gauge of the mean annual air temperature. The wealth of experimental field data, and the simplified mathematical structure implied by the reduced number of constituents, mean that dry snow packs provide a useful test bed in which to develop a rational model. Morland et al. (1990) laid down a rigorous theoretical framework for a four constituent phase-changing snow pack. The conservation

equations for the ice, water, water-vapour and air were derived from the principles of mixture theory (Morland, 1972, 1978) and the reduced three-dimensional theory, in which each constituent has the same common temperature T, consists of a coupled set of seventeen partial differential equations written in terms of fifty-three independent variables. Closure of the system requires a variety of physical assertions about the properties of the constituents, and about their interactions, and further simplifying assumptions are needed to develop a tractable model. There has been no further work to date on the full three-dimensional system. However, a simplified one-dimensional approach has been used by both Jordan ( 1991 ) and Bader and Weilenmann (1992), in which vertical gradients exist but the snow pack is laterally uniform with zero horizontal velocity component. This is a very useful reduction which allows a numbers of

O165-232X/94/$07.00 © 1994 Elsevier Science B.V. All rights reserved SSDIO165-232X(93)EOO19-F

136

J.M.N.T. Gray, L. W.. Morland / Cold Regions Science and Technology 22 (1994) 135-148

problems to be investigated. These models represent the state of the art in wet snow pack modelling, but the presence of the water-vapour and air is not treated particularly well in either case. Bader and Weilenmann ( 1992 ) assume that both gases have zero velocity and that there is no sensible or latent heat contribution to the energy balance. Jordan ( 1991 ) again assumes stagnant gases, but has a latent heat contribution due to sublimation entering the energy balance. In this paper we attempt to improve the description of the gaseous components within the snowpack, and determine their importance in the energy balance. It is an approximate treatment in the sense that solar penetration and compaction are ignored and the emphasis is on the seasonal time scales so that shorter term effects are not described. An order of magnitude analysis highlights the dominant balances, which are the basis of a reduced mathematical model in dimensionless variables. For illustration we determine the temperature, air velocity, density and pressure fields when the surface is subjected to periodic pressure and heat flux variations.

2. Mixture theory A natural dry snow pack consists of a porous ice matrix with the space between the crystals filled with water-vapour and dry air. The ratio of water-vapour to dry air in a unit volume of the composite gas is determined by the vapour pressure (Colbeck, 1982), and is typically of the order of 10- 2 for sub-zero temperatures. The composite gas appears, therefore, to be dominated by the dry air, and in this paper we shall consider a simplified two constituent theory, consisting of ice and dry air, in order to gain a deeper understanding of the gas contribution to the sensible heat balance. It should be noted that there is no phase change between the ice matrix and the dry air, as they are chemically inert, so the rate of mass exchange between the two constituents is identically zero. It follows that the porosity can only change by compaction of the matrix, and that latent heat contributions to the energy balance are absent in this model. With this under-

standing it should then be feasible to consider a three constituent pack including water vapour. A three-dimensional two constituent theory analogous to the reduced theory of Morland et al. (1990) is easily formulated in the absence of phase change. In each element of the snowpack there is a given proportion of each of the two constituents, defined by the volume fraction of constituent v per unit mixture volume, ~". Departing slightly from the notational conventions established by Morland et al., the constituent letters v = i, a will be used to refer to the ice and air respectively. By definition the volume fractions lie between zero and one and their sum over all constituents equals unity:

0~0V~l,

¢i'q-oa= 1

(2.1)

Partial variables are defined per unit mixture volume or area,while intrinsic variables are defined per unit constituent volume or area. Partial density and stress are related to their intrinsic counterparts by linear volume fraction scaling

p~=~p~,

~r~=c~r ~

(2.2)

where the superscript constituents letters are lower case for partial and upper case for intrinsic variables. A detailed formulation of three-dimensional mixture theory with such interpretations is presented by Morland ( 1992 ). Both constituents satisfy a mass balance equation of the form

OP~+div(p'v~) 0t

=0

( v = i , a)

(2.3)

The constituent velocity fields v" are associated with the mixture in the sense that p"v v determines the mass flux per unit mixture cross-section, and a consistent theory requires that v" defines the mean intrinsic velocity of constituent v, so here there is no distinction between partial and intrinsic variables. The momentum balance for constituent v = i, a is div (~r ~) +p"g+p~ " = 0

(2.4)

where a " is the partial stress tensor of constituent v which determines a partial pressure p" = - ~tr (a" ). The gravitational body force is denoted by g, and p ~ ~ is the interaction drag ex-

J.M.N.T. Gray, L. W. Morland / Cold Regions Science and Technology 22 (1994) 135-148

erted on constituent v by the other constituent per unit volume of mixture. The two constituents are in a state of thermal equilibrium provided that the forcing time scale is much longer than the time scale for thermal transfers between constituents to take place, which is assumed here. Thus, both constituents have a common temperature T, and there is a single bulk energy equation

, L i p Cp (ff[+r .VT)-pr -p --Di-]

which balances the local and convected heating, radiation, pressure working and diffusion. The work done by pressure is commonly neglected, and we demonstrate later that it is indeed small for the applications considered. However, it is noted that the pressure working can be significant in other situations. For each constituent v, is the specific heat at constant pressure, r ~ is the rate of external energy supply per unit volume of mixture, and K X is the thermal conductivity. The convected derivative defined as

Cp~

DJDt

D, Dt

-

0 0t

Fv".V

(2.6)

is a measure of the change following a particle of constituent v.

ever, there are essentially thirteen unknown variables in the problem: two volume fractions, two densities, two vertical stresses, two velocities, two interaction drags, two external energy supplies and one temperature. To complete the system we must make five constitutive postulates and two external prescriptions. The first two conditions are the equations of state that determine vertical stresses in the ice and the air. We assume that the ice grains are rigid and have constant intrinsic density p l = 918 kg m - 3, so the vertical s t r e s s 033 i = 0"i is arbitrary, determined by the momentum balance. The air is assumed to obey the perfect gas law

pA=pARAT

pM~ = - p ~

(3.2)

Morland (1978) determined the form of the interaction drag by Darcy's law for a gas with viscosity/z and matrix permeability k. In one dimension the horizontal components are identically zero and the vertical component is A0~a

We focus here on a one-dimensional vertical snow pack that has no lateral gradients and whose horizontal velocity components are identically zero. Let z denote a vertical co-ordinate measured in the upward sense, with vertical components of vectors denoted by a subscript 3. In this system the horizontal momentum balances are trivially satisfied. Thus, the conservation principles yield two mass balances, two vertical momentum balances and one bulk energy balance, a total of five equations, and there is in addition one volume fraction relation (Eq. 2.1 ). How-

(3.1)

where R A is the gas constant for air, and to be non-viscous, supporting no shear stress, so the stress 0"33A= _ p A . The third condition requires the interaction drag experienced by the ice due to the motion of the gas, p~i, to be equal to but opposite in sign to the interaction drag experienced by the gas while passing through the ice, p ~ i; thus

PN~=P Oz 3. Constitutive properties

137

].~ a 2

k (v~-v~)

(3.3)

which is the fourth condition. In the presence of gravity the ice matrix collapses under its own overburden pressure, leading to densification of the snowpack in the absence of phase change. Yoshida et al. (1956), Bader (1962) and Kojima ( 1964 ) have tried to interpret this behaviour by means of a linearly viscous stress relation with a "compactive" viscosity t/c, dependent on porosity and temperature. Under normal conditions 101°< qc < 1017 kg m -i s -l

(3.4)

However, Mellor (1975) concluded that the properties were of such bewildering complexity that it was sensible to adopt a greatly simplified

138

J.M.N.T. Gray, L. W. Morland/ ColdRegions Science and Technology22 (1994) 135-148

description. Therefore, instead of prescribing a snow rheology which relates volume fraction changes to stress, we adopt the simplification that the ice matrix remains stationary: v~=O

(3.5)

which is the fifth constitutive condition. This assumption precludes the withdrawl of gases due to matrix compaction, but allows us to investigate the effects of surface temperature and pressure forcing in isolation. In the same spirit, the final two conditions eliminate solar radiation effects in the near surface layers, i.e. the external energy supplies r ~= 0, ( v = i, a).

4. Model equations The conservation laws Eqs. (2.3)-(2.5) in conjunction with the volume fraction relation Eq. (2.1), the relations in Eq. (2.2) between partial and intrinsic variables, and the five constitutive properties and two external energy supply prescriptions in Section 3, define a rational closed system of equations. These are now reduced to a system of four equations in four unknown variables. The ice mass balance equation (2.3) in the absence of phase change, together with the neglect of matrix compaction (3.5), reduces to OP~-o at

(4.1)

The partial and intrinsic densities are related by the volume fraction scaling pi=f~ipi, and so, for constant intrinsic ice density pi, Eq. (4.1) implies that the ice volume fraction ~i is independent of time. The volume fraction relation (2.1) then implies that the porosity q~a,equal to 1 -Oi, is also independent of time. Thus, an initially prescribed porosity distribution Oa=Oa(Z) is fixed for all time. Spatial porosity gradients still contribute to the air mass balance (Eq. 2.3), which becomes O a ~ - q- ~ZZ(¢ pAV~) = 0

aOPA

(4.2)

The ice velocity field is identically zero by property (3.5) and the ice momentum balance (Eq. 2.4) is an equation for the partial ice stress a33i required to maintain a stationary ice matrix. This is not a physical quantity that influences the air response, so Eq. (2.4) can be ignored. Substituting the Darcy drag relation (3.3), with v3i=0, into the vertical air momentum balance, and combining the partial pressure gradient term and volume fraction gradient term, 0p A

Oz

.p Ag--k~ au~ = 0

(4.3)

The convected derivatives D~fb~/Dt (v=i, a), defined in Eq. (2.6), are greatly simplified in the absence of lateral gradients, time independent volume fractions and a stationary matrix, being simply

D~

a0Oa

~ - - v 3 ~-z'

Di~i-0 Dt

(4.4)

Thus only the air pressure contributes to the pressure working terms in the energy balance (Eq. 2.5 ), which becomes

(l__Oa)pIClpOT+oaoA A loOT aOr~ Ot ~ r Cp ~ ' [ - U 3 ~ z ) a - ~-Z[( (1-Oa)KI-I-(baKA)~I=O --p AV3a ~~Z (4.5) recalling that external energy supply to both of the constituents has been neglected. It should be noted that the volume weighted average conductivity ( 1 - 0 a ) K I + ( b a K a is not consistent with measured conductivities of the snow (Yen, 1962; Mellor, 1964). The conservation equations (4.2), (4.3) and (4.5), together with the equation of state for the air, Eq. (3.1) (the ideal gas law), now represent a closed system of four equations for the four variables pA,pA, T and va3.

5. Non-dimensional variables Natural dry snow packs are forced by the prescribed atmospheric conditions prevailing at

139

J.M.N. T Gray, L. W. Morland/ ColdRegions Science and Technology22 (I 994) 135-148

their surface. Our interest lies in predicting the temperature profile in the snow for a variety of surface pressure and temperature scenarios. On the high Antarctic plateau, surface temperatures range between 233 and 273 K, with annual mean Tm = 253 K and amplitude of fluctuation T* -~ 20 K, suggesting a temperature scaling of the form T=Tm(1 +El T)

(5.1)

where the ratio of fluctuation amplitude to mean temperatures defines e~=T*/Tm. The tilde is used to indicate a non-dimensional variable. Surface pressures lie in the range 7.8N1048.2X104 Pa (780-810 mbar) with mean p m - 8 X 10 4 Pa and fluctuation p*-~2×103 Pa. The air density responds via the ideal gas relation ( 3.1 ) and has mean density

magnitude over the depth z*. Suppose that the pressure gradient OpA/Oz in the air momentum balance (Eq. 4.3) is of order p* / z* ~ 103 kg m -2 s-2, then since the gravitational force Pmg~ 10 kg m -z s -2, the greater pressure gradient must be balanced by the Darcy drag ItOv3a/k. For air viscosity/~ = 1.45 X 10- s kg m - ~ s- ~ ( Batchelor, 1967) and ice matrix permeability k = 10 -9 m e (Shimizu, 1970 ), this balance implies an air velocity of magnitude v* = p ' k ~ (z*/t) __ 10- ~m s- 1. Observed velocities are much lower, so this is not the correct balance and the pressure gradient with depth must be smaller. This motivates a pressure scaling that separates the surface pressure variation/~h=ph(~'), which is independent of z, from a depth varying part/~(2, T). Thus pA=pm [ 1 + e2/~h(h + %/5]

Pm

Pm z- R ATm

(5.2)

which is approximately 1 kg m - 3 for R A = 2 . 8 7 × 1 0 2 J k g -I K -1 (Gill, 1982). In the high polar latitudes where perpetual daylight reigns in summer and perpetual darkness comes in winter, surface temperatures change by order T* on a seasonal time-scale t*---107 s and we introduce a non-dimensional time ~by t=t*~

t*= 107s

where z* = 1{ ~KoI t~* ]j 1/2 = 3 m

where e2=P*/Pm = 0.025 and the surface value /~(0J) =0. This demonstrates that the surface value /3h(~') is transmitted instantaneously through the snow pack. The pressure gradient is now of magnitude e3 Pm/Z*,and can be no larger than the gravitational force Pmg. Allowing this maximum magnitude, and using Eq. (5.2), shows that ~ . 3 = g z * / ( R A T m ) , ~ 4 x I O -4. The ideal gas relation (3.1) then implies a density scaling and separation of the form

(5.3)

A balance between the conduction/dO 2T/Oz 2 and the local heating plCpIOT/Ot terms in Eq. (4.5) determines the length scale z* over which surface temperature changes penetrate in time t*. The vertical co-ordinate is therefore scaled as z=z*~,

(5.5)

pA ^ {l+e2/~h'~.t_^ e = =Pmk i ~ _ ~ ] //m 3l/

(5.6)

where the non-dimensional density fi is

fi-(l+e~T)

(5.7)

(5.4)

for Kt=2.2 J m -~ s - 1 K - i (Hobbs, 1974; Glen, 1974) and CpI=2.093X 103 J kg -1 K - ' . The temperature penetration length scale z* is typically shorter than the snow depth h at which pore space ceases to be interconnected, at a porosity c a _ 0.1 when gas is trapped within the ice. Surface pressure is assumed to change by order p* in time t*, but this does not necessarily imply that the pressure will change by the same

with the surface value ~(0,t')=0. The three dimensionless parameters governing the balances are *

el=T*'007, Tm "

p*

e2=--=0.025 Pm

gz~* ~4X10-4<<1 1~3=RATm

and (5.8)

J.M.N.T. Gray, L. W. Morland / Cold Regions Science and Technology 22 (1994) 135-148

140

of which only e3 is very small; the magnitudes El, e2 are not negligible. To leading order in the small parameter e3, the air density gradient is

0PA

~lPm (l÷E2ffh) 0T

0z -

z* ( 1 + ~ 1 ~ ) 2 0-~

(5.9)

from which it follows that decreases in density with height are associated with 0T/0~>0, and vice-versa. In the latter case OpA/Oz> 0, the air column is unstably stratified, which, if the motion was not constrained to the vertical, could convectively overturn. This effect is excluded by the present one-dimensional treatment. The air velocity magnitude v* has still to be determined. We assume a scaling of the form

v~---v*~

(5.10)

and substitute this, together with the density scaling (Eq. 5.6 ), into the air mass balance equation (4.2). To lead order in the small parameter ¢3, Eq. (4.2) becomes

whole of the air column. However, the estimate of v* =pmgk/~ assumes that the Darcy drag contributes to the leading order balance, which is now contradicted. Thus this estimate of v* is not valid. The velocity magnitude v* is not then set by the momentum balance at all, which must be a pressure gradient-gravity balance, but by a balance between the two terms in the leading order mass balance equation (5.11 ). We suppose that the motion is driven by the surface conditions, and not dominated by the initial porosity gradients. Then, if surface pressure forcing dominates, ~2

~-

~-

(l+~2/~h)~

and since the surface pressure forcing is independent of 2, the surface velocity is given by g

1

(5.11) It has already been established that Pmg is the largest term in the air momentum balance equation (4.3). We first consider a balance between this and the Darcy drag/t0 av3a/k would imply an air velocity of magnitude v*=pmgk/gt ~ - 10 - 3 m s-1, and the parameter v*t*/z*,,, 10 3 ~ ~3--1 For this v*, Eq. (5.11 ) becomes to leading order in

~3

0 ( ~1+~2~h-~ O--~_O ]~v)=o

(5.12)

dp

.

(5.15) Therefore, since for a snow depth h=z*h the scaled surface velocity ~ is order unity by construction, surface pressure forcing induces a velocity magnitude v*=e2h~-2x 10-Sms -1 t*

(1+,l~) ~ oF - \ ~ - / ( l + , ,

~a 1 +ezffh

( v*t*~

(5.13)

In addition we require that at some depth, 2= 0 say by choice of origin, the pack becomes impermeable and ~=0. Applying this boundary condition we find that f(i') ---0, which in turn implies that to leading order ~= 0 throughout the

(5.16)

for a snow depth of 8 m. Alternatively, if thermal effects dominate, the major contributions to the mass balance equation ( 5.11 ) are

with solution that is a function of time only, namely

l+~ ~ ~=f(?)

(5.14)

1

~): oe

Of

Z~] 1..l_El ~ 0 = 0

(5.17)

The third term is of magnitude v*t*/z*, which is El-1 larger than the second term, therefore the leading order balance lies between the third and first terms implying a velocity magnitude

J.M.N.T. Gray, L. W. Morland / Cold Regions Science and Technology 22 (1994) 135-148

v*= ~

Z~

t*

= 2 X I 0 - 8 m s -1

(5.18)

Both pressure and temperature induced velocity magnitudes (Eqs. 5.16 and 5.18) are therefore the same, and either one can be used in the velocity scaling (Eq. 5.10). Note, however, that pressure changes here have been assumed to take place on a time scale t*= 107 s, so if pressure changes of order p* occur on a shorter time scale t +, due to the passage of a large weather system or storm, then the velocity must be normalised using v*= e2h/t +. In summary, the relevant physical magnitudes are T * = 2 0 K, p * = 2 × 103 Pa, t*= 107 S, Z'm3 m and v* = 10- 8 m s- 1. The velocity scaling procedure, described above, holds provided that the Darcy drag term It~a/33a is less than the gravitational force pag in the m o m e n t u m equation (4.3). Very rapid variations of the surface pressure, caused by strong winds blowing over surface features like dunes or ripples, can result in the Darcy drag balancing the pressure gradient. If pressure variations occur on a time scale t +, the condition for a dominant Darcy drag term is that the velocity V+~2 h -

Pmg

__ 10_3 m S_l

(5.19)

t + >1 ( i t ~ k )

and the present scaling procedure breaks down. Colbeck (1989) has investigated this process, known as wind pumping, in the case of constant temperature and constant porosity. The present model and illustrations apply to a different set of applications, of equal importance and complementary to those of Colbeck. 6. D i m e n s i o n l e s s e q u a t i o n s

The scaling relations (5.1), (5.3), (5.4), (5.5), (5.6) and (5.10) are now substituted to obtain the model equations in terms of dimensionless variables. The ideal gas relation (3.1) becomes /3=/~(1 + E, T) and the air mass balance (Eq. 4.2) is

(6.1)

~a 0 [-l+e2Ph7. a0/0

141

0 F o,al'l-~2ffh~7 0

a

+,3 {0 ~ ~ + , 1~-~(0 /~tT)} = 0

(6.2)

The air momentum balance (Eq. 4.3 ) reduces to

0/~ ( 1 "1-'2/3h ) a~= - - ~ - - k i ~ l ~ ~-e3/O --e4~Z) 0

(6.3)

where the non-dimensional parameter

E4 = . -(#/k)v* --

(6.4)

Pmg

is small provided v* <
0T

(X-t-ez/~h ~ ) ( 0 7

~0T)

~ ~0¢ a

--~6[ 1 "[-6.2ph q-e3p]v ~Z

o?

- - ~ [( ( 1 -- ¢ a) "~-6 7 ¢ a ) ~ ] = 0

(6.5)

where three more non-dimensional parameters are introduced. Firstly pmC A

65 __ p~C-pTpi,~ 5× 10 -4

(6.6)

with CAp= 1.012× 10 3 J kg -1 K -1. The importance of the pressure working term is determined by Pm e6 --PI Cip

( v*T*~ T* \ ~--]

(6.7)

which is of order 2 × 10 -4 for the velocities induced by temperature variations or slow pressure change, which are investigated here. However, for faster pressure forcing on a time scale t +, when Eq. (5.19) applies, •6 exceeds order

142

J.M.N. T. Gray, L. W. Morland / Cold Regions Science and Technology 22 (I 994) 135-148

unity and Colbeck's (1989) assumption of constant temperature does not hold when the snow porosity is non-uniform. To properly investigate the process of wind pumping and assess the contribution of the pressure working term in the energy balance, a full non-dimensional analysis, similar to that presented in this paper, is still required. Finally, the ratio of conductivities defines the third non-dimensional parameter KA '7 = ~ i - ' ~ 10 -2 (6.8) w i t h K a = 2 . 2 4 × 10 -z J m -1 s -1 K -1 The non-dimensional equations have four small parameters '3, ,4, ,5 and '6, that can be exploited to uncouple the four non-dimensional equations (6.1), (6.2), (6.3) and (6.5). Each variable is approximated by the first two terms of an asymptotic expansion in the small parameter 6,

7= 7o+67,

#=#o +6#,

/~=/~o+ a/~,

} (6.9)

which, given To, is an equation for z7o. The nondimensional pressure in the momentum balance (Eq. 6.3 ) is eliminated in favour of~o by substitution of Eq. (6.1), giving to leading order

_o ( o(l +,, 7o)) =_(!\ I 02

]

(6.13)

+,~ To/

which determines Po. The leading order pressure can then be calculated from Eq. (6.1) once Po and To are determined: Po =rio( 1 + ' , To)

(6.14)

Taken in order, the left-hand sides of Eqs. ( 6.11 ) - ( 6.14 ) are all written in terms of known or previously calculated variables, allowing To, bo, Po and/~0 to be determined in succession. The first order equations follow the same pattern but are considerably more complex. To first order the energy balance is

(1--~ a) - - ~ -- ~

(1--oa)'t-'7¢")

where

3-----max(E3,~4,,5,E6) ~ 5 X 10 -4

(6.10)

The parameters ,~, E2, and '7 are less than unity, but are retained as possible order unity quantities. Substituting the approximations (6.9) into the non-dimensional equations and matching zero andfirst order powers in ~ yields a leading and first order system of equations for the eight variables. To leading order the energy balance (Eq. 6.5 ) uncouples temperature from the other variables, giving (1--

~a)070

0

--I(

(1--q~a)+'vq~

a)070/

0 --~-J---(6.11)

so that To can be determined once suitable boundary conditions are prescribed. The leading order mass balance (Eq. 6.2) becomes

O /[t~a 1 +'2/0h ~

1

0a

= ('-~-)[ 1"l" ,2/~h] ~00~aOZ ['5"~

a [ l +'2Ph'~ { OTo

-C7) +

(6.15) the mass balance becomes

2[0.1+,2 h- 1 1

aOfiO

0

a

"3t-oa~t--[(}'3L'2,~P~hx271]l..~_,llo) 0F.a

°,~_L o (6.12 )

_ 0 7o'~

o, +"'o-Cs-z)

l+'2Ph ---. ] (1 + , , 7o)2IlVoJ

the air momentum balance reduces to

(6.16)

J.M.N. T. Gray. L. W. Morland / Cold Regions Science and Technology 22 (1994) 135-148

a 0T

1 + E2Ph .,p O ( p , ( 1 +~, To) ) -~-1~1( 1 -]-E 1 To) 2-1

(2= ( ( 1 - ~ a ) + e T ~ ) f f ~ = 4 ×

--¢,~0(~oP]P1)- (~)fio--( 8)0¢4a~vo (6.17, and the pressure p, =~lfi01F, +p, (1 +e, iPo)

(6.18)

which allow T~, 17~,~ a n d / ~ to be determined successively. The leading order solutions are accurate to order 8, which is sufficient for the purposes of this paper, and all illustrations are therefore based on a solution of Eqs. (6.11)(6.14).

7. Boundary conditions

~=0:

~o=0

(7.1)

which allows the coefficient of integration in the leading order mass balance (Eq. 6.12) to be determined. The base of the snow pack also experiences a prescribed geothermal heat flux (Paterson, 1969) of the form

Q= ( ( 1-- q~a)Kl+ o a K A ) ~ - 4 × 10 -2 J m -2 s - I

(7.2)

(7.4)

which is slightly larger than 6, but still much smaller than unity, so for the purposes of this paper the leading order basal boundary condition for the energy balance (Eq. 6.11 ) is assumed to be

e=o:

°T°--o 0~

(7.5)

The surface of the snow pack is subjected to atmospheric pressure fluctuations, determined by the prescribed function ,Ohin the pressure scaling (Eq. 5.5). In addition, sincepA=pm[ 1 + e ~ h ] at the surface, the third boundary condition requires that 13= 0 at ~ =/~. To leading order 2=~':

The non-dimensional equations (6.1), (6.2), (6.3) and (6.5) require two initial conditions and four boundary conditions in order to compute a solution. Recall that the surface of the snow pack lies at z=h, and by choice of origin the base of the snowpack lies at z = 0, where the pore space ceases to be interconnected and #3 = 0. The first boundary condition to leading order in 3 is

10 -3

143

/~0 = 0 ~ rio = 0

(7.6)

by Eq. (6.14). The final boundary condition necessary to solve the energy balance (Eq. 6.5) is a heat flux condition applied at the surface of the snow, and is chosen because it is physically more realistic than simply prescribing the surface temperature. The prescribed slow surface pressure variation /~h and the surface heat flux, both of which change with the prevailing atmospheric conditions, are the dominant forcing terms that drive the thermal and mass transport regimes within the interior of the snowpack. Two sets of surface boundary conditions are therefore considered. The first is designed to investigate the thermal regime over an annual cycle when the surface pressure is held constant, ,Oh- 0. This is achieved by prescribing an idealised sinusoidal heat flux at the snow surface 2=t7:

~= acos~- )

(7.7)

and is scaled as

Q=Q*O_

where Q * = ( - - ~ ) =

1 0 J m -2 s -~ (7.3)

Thus, the non-dimensional heat flux is

with period T= 4 equivalent to one year in physical variables. The order unity factor a is included to allow the applied heat flux to produce a maximum temperature equal to unity. To leading order the surface boundary condition for the thermally driven problem is

J.M.N. 7". Gray, L. W. Morland / Cold Regions Science and Technology 22 (1994) 135-148

144

a

~'=E:

~-=acos~,~- )

/3h--0

(7.8)

The second problem is designed to investigate the effect of slow pressure changes that occur over an annual cycle. The non-dimensional temperature ]ro is initialised to a constant throughout the whole of the snowpack and no heat flux is input at the surface, and so the temperature remains at its initial value. A slow sinusoidal surface pressure variation is then used to drive the system. To leading order in J

B(5) = ~ . ~ (cosfl(f/+~) sinhfl(~'-5) + cos,8(fi- 5) sinhfl(tT+ 5) - s i n f l ( ~ + 5 ) coshfl(~'-5) - sinfl(h-- 5) coshfl(tT+ ~) )

(8.4)

and the constants are fl=~c

s=l(cosh2fl~-cos2fl~)

and

(8.5) 0~---0

/~h = sin ~-

(7.9)

which again has period ~'= 4 (one year).

8. Analytic solution

An analytic temperature solution can be obtained in the case of constant porosity 0a(2,i,) = ¢, when the energy balance (Eq. 6.11 ) reduces to

0To 02 o Or" --c 0~2

This, for example, can be derived as the long time asymptotic solution from a formal Laplace transform solution of a corresponding initial value problem. A complete cycle of this periodic solution (period 4 time units) for snow depth h = 4 is illustrated in Fig. 1, shown between time units 8 and 12, in the form of a contour plot in space and time of the non-dimensional temperature To. The constant a has been chosen to make the maximum temperature equal to unity. The m a x i m u m heat flux is supplied at the surface at the beginning of the cycle r= 8, and it takes time to heat the snow, so that the maximum temperature does not occur until approximately 8.5 time

i

(8.1)

9 4

t

I

t

I

10

I

t

11

I

12 ,I

\

4

where the constant c = 1 + e 7 ¢ / ( 1 - ¢ ) ~ 1. The thermally driven problem, in which zero heat flux is applied at the base (Eq. 7.5 ) and a sinusoidal heat flux (Eq. 7.8) is applied at the snow surface, has a long time oscillatory solution of the form

~?

nT

To =A (ff)sin(~-) + B ( ~ ) c o s ( ~ - )

(8.2)

where the spatially dependent coefficients are A(5) = ~

a

1

\

\,,,

\

(cosfl(h'+5) sinhfl(tT-~)

+ cosfl( ~'- ~) sinhfl(h'+ 5)

I

I 9

I

I

\

\

-1 \

I

I 10

i ~ l

I 11

~'~

0 12

+ sinfl(//+ 5) coshfl(~'- 5) + sinfl(/~- ~) coshfl(//+ 2) )

( 8.3 )

Fig. 1. Non-dimensional temperature contours for periodic solution.

J.M.N.T. Gray, L. W. Morland / Cold Regions Science and Technology 22 (1994) 135-148

units. The surface temperature lags behind the heat flux by half a time unit. Heat is conducted down through the snowpack losing intensity and being phase shifted forward in time with increasing depth. At the base of the snowpack the temperature changes only slightly from zero. This analytic periodic solution acts as a good check on the numerical solution, which should approach this solution asymptotically as the transients decay.

9. Numerical illustrations A numerical model has been developed to solve the leading order system (Eqs. 6.11-6.14) subject to the boundary conditions (7.1), (7.5) and (7.6), with either thermal (Eq. 7.8), pressure (Eq. 7.9), or mixed thermal and pressure forcing applied at the snow surface. The energy balance (Eq. 6.11 ) is solved for the temperature To by a standard implicit technique. The right-hand sides of Eqs. (6.12 ) and (6.13 ) are now written in terms of known variables, and the velocity g0 and density Po can be computed by quadrature. Finally, the pressure Po is determined from the algebraic relation (6.14). To initialise the leading order problem it is assumed that the temperature is constant and that the velocity is identically zero throughout the whole of the snow pack. Given that the prescribed surface pressure is constant, the energy balance (Eq. 6.11 ) and the mass balance (Eq. 6.12 ) are trivially satisfied provided that To(5,0)-Tc

and

go(5,0)=0

(9.1)

The momentum balance (Eq. 6.13 ) together with the boundary condition (7.6) imply that the nondimensional density is 1 + •2Ph /~o - ( i ~-_~ ~-) ~ (/~- 5)

(9.2)

and then the non-dimensional pressure is

1 + ~2/~h (g--5) /~o - 1 +el ~c

(9.3)

by Eq. (6.14). Note that the constant initial sur-

145

face pressure and temperature are both assumed to be at their mean values, i.e. /~h(0)=0

and

Tc=To(5,0)=0

(9.4)

in all calculations presented here. The model is integrated forward in time from this initial steady-state, the transients soon decay, and the model quickly settles into a long time oscillatory solution. After eight non-dimensional time units the solution has reached a periodic cycle with a period of four non-dimensional time units (one year). Results are presented for the range 8-12 time units and snow depth h = 4 corresponding to Fig. 1. Although the numerical model is able to handle non-uniform porosities, for the purposes of illustration a uniform porosity ¢a=0.5 is considered. The thermally forced problem, driven by the surface boundary condition (7.8), has a temperature solution that asymptotically approaches the long time oscillatory solution (Eq. 8.2), illustrated in Fig. 1, for the case of constant porosity. This is in general agreement with Mellor's (1964) illustrations of Schytt's (1958) Antarctic measurements. To lead order in the small parameter el, though larger than 3, the mass balance equation (6.12 ) implies that

0~o 0To 02To 05 ~ 0 7 - ~05

(9.5)

and integrating with respect to 5 with the boundary condition (7.1) gives a velocity 0To Vo ~ 05

(9.6)

Thus, at the surface the velocity maxima, minima and zeros correspond, with an error of order 1, with the applied heat flux 0~Po/05= cos (nU2). The velocity cycle, illustrated in Fig. 2, has the same features as the temperature cycle, but is phase shifted backwards in time by half a time unit. The non-dimensional density, illustrated in Fig. 3, is perturbed only slightly from its linear distribution in 5 given in Eq. (9.2). The lower densities experienced between 8 and 10 time units are a direct result of the integral

146

J.M.N.T. Gray, L. W. Morland / Cold Regions Science and Technology 22 (1994) 135-148

i 9 t

I

I

fi

10

I

I

11 I

I

I

I

J

1

12 |

I

I

I

0

~ d~ 1 -I-61To

(9.7)

which is smaller when the total heat content is high. The pressure field differs only slightly from this solution and is not illustrated. The case of prescribed surface pressure forcing, Eq. (7.9), is much simpler than the thermal problem and can be solved analytically. The temperature is determined by the energy balance (Eq. 6.1 1 ) subject to the heat flux conditions (7.5), (7.9) and the initial distribution (Eq. 9.4), and is identically zero for all space and time: ~Po(e,7) = 0 I

I

I

I

I

r

t

i

9

I

i

i

I

I

10

I

i

I

i

0

I

11

12

i Fig. 2. N o n - d i m e n s i o n a l velocity contours during third cycle for thermal forcing.

to lead order in 5. The effects of surface pressure variation are transmitted instantaneously through the whole gas column, causing uniform expansion or contraction, and result in a flow of air. The air mass balance (Eq. 6.12 ), subject to the zero velocity at depth condition (7.1), can be solved analytically in the case of constant uniform temperature and porosity. Thus, the velocity is Vo=

8

9

10

11 I

I

I

I

I

I

I

I

-2

-1

i

i

i

I 9

i

i

i

i

I

i

10

11

(9.9)

4

-3

i

E2/~l dPh~ l+e2Ph dt "z

12

=-2

8

(9.8)

12

i Fig. 3. Non-dimensional density contours during third cycle for thermal forcing.

[ 1 + e2sin(~zt'/2) ] c°s

ff

for the prescribed surface forcing (Eq. 7.9 ). The air flow has a linear spatial distribution reaching its maximum modulus at the surface when d~h/dt is at a maximum or minimum, T even in this case, and the velocity is identically zero when d~h/dt=O, at/'odd. This is shown in Fig. 4. The maximum velocity magnitude is larger than in the thermal problem, which is due in part to the large snow depth at constant porosity. Real snow packs have decreasing porosity with depth, and as the matrix structure becomes more tortuous, the permeability decreases impeding the flow of air. The non-dimensional pressure and density are identical by Eq. (6.14). Their distribution, given by integrating the momentum balance (Eq. 6.13 ) with the boundary condition (7.6), is once again linear in ~ with the form

J.M.N. T. Gray, L. W. Morland / Cold Regions Science and Technology 22 (1994) 135-148

Po =/~o = (l + ~2Ph)(fi-~) 10

8 4

12

14

.L

t

147

(9.10)

illustrated in Fig. 5.

10. Conclusions

/ /

j/

0 8

9

10

11

12

Fig. 4. Non-dimensional velocity contours during third cycle for pressure forcing.

8 4

I

J

I

i

9 I

t

i

t

t

3-

lO I L

11 I

I

I

I

12 I

I

I

I

4

3

52}

z

2

]

9

i

i

i

i

I

lO

i

i

i

i

I

11

i

i

i

o

i

12

Fig. 5. Non-dimensional density contours during third cycle for pressure forcing.

A model for a two-constituent dry snow pack has been formulated from mixture theory and supplemented with additional physics to obtain a rational closed system of partial differential equations. The effects of surface temperature and pressure forcing have been assessed using nondimensional analysis for slow, time scales of order one year, variations. A reduced model has been obtained by a regular perturbation, and is described by Eqs. (6.11 ) - (6.14), allowing the temperature, velocity, density and pressure to be calculated, accurate to order 5 -~ 5 × 10 -4. The air conductivity enters the energy balance, but to leading order there are no convective terms and no working terms involving the air velocity, and so the energy balance is the same as for a stagnant gas, corresponding to the original direct assumption of Jordan (1991) and Bader and Weilenmann (1992). However, it is apparent that for faster time scale pressure variations the velocity magnitude is increased and the work done by a gas travelling through a non-uniformly porous snow pack can make a significant contribution to the energy balance. Over an annual cycle the density gradient Eq. (5.9) is unstable for half the time. One-dimensional models are not able to describe convection processes, which would overturn the gas, and further work is necessary to determine whether the velocities induced by convection could transport enough heat to affect the temperature distribution. Water vapour has not been included in this study as at these temperatures it has a low partial pressure and would occupy a very small part of the gas volume fraction. The contributions to the specific heat and conductivity are undoubtedly small, but latent heat contribution through phase change with the ice matrix could well be significant.

148

J.M.N.T. Gray, L. W.. Morland / Cold Regions Science and Technology 22 (1994) 135-148

Acknowledgements This research has been s u p p o r t e d by the Natu r a l E n v i r o n m e n t R e s e a r c h C o u n c i l , a n d h as benefitted from the continuous interaction with D r . E . M . M o r r i s at t h e B r i t i s h A n t a r c t i c S u r v e y , Cambridge.

References Bader, H., 1962. Theory of densification of dry snow on high polar glaciers, II. CRREL Research Report 108. Bader, H.P. and Weilenmann P., 1992. Modeling temperature distribution, energy and mass flow in a (phasechanging) snowpack. I. Model and case studies. Cold Reg. Sci. Technol., 20:157-181. Batchelor, G.K., 1967. An Introduction to Fluid Dynamics. Cambridge University Press, Cambridge, 594 pp. Colbeck, S.C., 1982. An overview of seasonal snow metamorphism. Rev. Geophys. Space Phys., 20( 1 ): 45-61. Colbeck, S.C., 1989. Air movement in snow due to windpumping. J. Glaciol., 35( 120): 209-213. Gill, A.E., 1982. Atmosphere-Ocean Dynamics. International Geophysics Series, Vol. 30. Academic Press, New York, 40 pp. Glen, J.W., 1974. The physics of ice. CRREL Monograph IIC2a. Hobbs, P.V., 1974. Ice Physics. Oxford University Press, Oxford. Jordan, R., 1991. A one-dimensional temperature model for a snow cover: technical documentation for SNTHERM.89. CRREL Special Report, 91-16, pp. 1-49.

Kojima, K., 1964. Densification of snow in Antarctica. In: M. Mellor (Editor), Antarctic Snow and Ice Studies. Antarctic Research Series, 2. American Geophysical Union. MacDowall, J., 1960. Some observations at Halley Bay in seismology, glaciology and meteorology. Proc. R. Soc. A, 256: 147-197. Mellor, M., 1964. Properties of snow. CRREL Monograph III-A1. Mellor, M., 1975. A review of basic snow mechanics. In: Snow Mechanics Symposium (Proceedings of the Grindelwald Symposium, April 1974). IAHS-AISH, 114:251-291. Morland, L.W., 1972. A simple constitutive theory for a fluid saturated porous solid. J. Geophys. Res., 77: 890-900. Morland, L.W., 1978. A theory of slow fluid flow through a porous thermoelastic matrix. Geophys. J. R. Astron. Soc., 55: 393-410. Morland, L.W., 1992. Flow of viscous fluids through a porous deformable matrix. Surv. Geophys., 13: 209-268. Morland, L.W., Kelly, R.J. and Morris, E.M., 1990. A mixture theory for a phase-changing snowpack. Cold Reg. Sci. Technol., 17: 271-285. Paterson, W.S.B., 1969. The Physics of Glaciers. Pergamon Press, Oxford, 211 pp. Schytt, V., 1958. Snow and ice temperatures in Dronning Maud Land. Glaciology II, Scientific results of the Norwegian-British-Swedish Expedition, IV, Norsk Polarinstitutt, Oslo. Shimizu, H., 1970. Air permeability of deposited snow. Contributions from the Institute of Low Temperature Science, No. 22. Hokkaido University, Sapporo. Yen, Y.C., 1962. Effective thermal conductivity of ventilated snow. J. Geophys. Res., 67(3): 1091-1098. Yoshida, Z. et al., 1956. Physical studies on deposited snow. Contributions from the Institute of Low Temperature Science, No. 9. Hokkaido University, Sapporo.