Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Computer modeling of temperature profiles in freezing ground Webb, Fern Marisa 2004

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2004-0682.pdf [ 11.99MB ]
Metadata
JSON: 831-1.0052394.json
JSON-LD: 831-1.0052394-ld.json
RDF/XML (Pretty): 831-1.0052394-rdf.xml
RDF/JSON: 831-1.0052394-rdf.json
Turtle: 831-1.0052394-turtle.txt
N-Triples: 831-1.0052394-rdf-ntriples.txt
Original Record: 831-1.0052394-source.json
Full Text
831-1.0052394-fulltext.txt
Citation
831-1.0052394.ris

Full Text

Computer Modeling of Temperature Profiles in Freezing Ground by Fern Marisa Webb B.Sc., The University of British Columbia, 1997  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES EARTH AND OCEAN SCIENCES We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA September, 2004 © Fern Marisa Webb, 2004  /  THE UNIVERSITY OF BRITISH COLUMBIA  FACULTY OF GRADUATE STUDIES  Library Authorization  In p r e s e n t i n g this thesis in partial fulfillment of the requirements for an a d v a n c e d d e g r e e at the University of British C o l u m b i a , I a g r e e that t h e Library shall m a k e it freely available for reference a n d study. I further a g r e e that p e r m i s s i o n for e x t e n s i v e c o p y i n g of this thesis for scholarly p u r p o s e s m a y b e g r a n t e d by t h e head of m y d e p a r t m e n t or by his or her representatives. It is u n d e r s t o o d that c o p y i n g or publication of this thesis for financial g a i n shall not be allowed without m y written p e r m i s s i o n .  T h e University of British C o l u m b i a V a n c o u v e r , BC  Canada  grad.ubc.ca/forms/?formlD=THS  page 1 of 1  last updated: 20-Jul-04  ABSTRACT  Greater than 50% of the area of Canada is underlain by permafrost, a thermal condition defined by mean ground temperatures remaining below 0 for a minimum of two consecutive years. The condition of frozen ground bears on many aspects of northern ecology, climate, engineering and society. The temperature based definition of permafrost highlights that understanding the condition of permafrost requires understanding the temperature distribution and energy balance of the ground. Physically based numerical modeling of earth systems is a tool for understanding how past geoclimate conditions have produced current features, and how prospective changes in forcing might manifest future changes in landscape or climate. 0  I have developed a numerical model to solve for a one-dimensional temperature distribution responding to time-dependent boundary conditions. Novel features of the model are a coordinate transformation which allows for a spatially mobile upper domain boundary, and a constituent mixture approach to define temperature dependent thermophysical soil properties. The model development is guided by a desire to minimize the stringency of input data requirements due to the sparse availability of quantitative information on soil properties and surface conditions. A variety of model applications are demonstrated using synthetic simulations and real world data.  n  Table Of Contents  ABSTRACT  ii  Table of Contents  iii  List of Tables  v  List of Figures  vi  Chapters  1  1 INTRODUCTION  1  1.1 Context and Goals 1.2 Terminology  1 2  2 PHYSICS OF GROUND TEMPERATURES  2.1 Energy Balance and Heat Flow 2.1.1 Energy Processes in the Ground 2.1.2 Equation Development for Energy Balance 2.2 Boundary Conditions on the Ground Temperature Profile 2.2.1 Qualitative Description and Models 2.2.2 Boundary Conditions Applied to the Governing Equation 2.3 The Thermodynamic Material Properties of the Ground 2.3.1 The Layered Soil Column 2.3.2 The Constituent Mixture of Materials 2.3.3 The Phase Change of Soil Water 2.3.4 Thermal Properties of Snow 3 NUMERICAL DEVELOPMENT  5  5 5 6 10 10 15 17 17 18 24 27 31  3.1 The Discrete Staggered Grid 3.2 The Governing Equations 3.2.1 Algebraic Manipulation 3.2.2 Finite Difference Discretization 3.3 The Boundary Conditions  31 31 31 33 34 iii  3.3.1 The Upper Boundary Condition 3.3.2 The Lower Boundary Condition 3.4 The Phase Change 3.5 The Discrete Property Fields From the Layered Column 3.6 Testing 3.6.1 Testing the Diffusion Code 3.6.2 Testing the Phase Change 3.6.3 Testing the Mixing Equations 3.6.4 Testing the Mass/Energy Conservation 4 APPLICATIONS 4.1 Synthetic Applications. 4.1.1 Sensitivity to Soil Column Parameters -Goodrich, 1982 4.1.2 Solving for Multiple Possible Soil Columns 4.2 Real World Applications 4.2.1 Mackenzie Delta 4.2.2 Mackenzie Valley 4.2.3 Southern Yukon Territory 4.3 Summary of Model Applications  34 35 35 37 39 39 40 42 44 47 47 47 54 60 62 69 74 82  5 CONCLUSIONS  83  References  84  Appendix  87  iv  List of Tables 2.1 Properties of common soil constituents 2.2 Typical bulk densities and porosities of common soil types 3.1 Comparison of measured and calculated values of thermal conductivity 4.1 Thermal properties of soil materials from Goodrich, 1982 4.2 Properties of variable layers in synthetic soil column ensemble 1 4.3 Properties of variable layers in synthetic soil column ensemble 2  v  18 21 44 48 56 58  List of Figures  1.1 1.2 1.3 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11  Relationship of frozen ground to other fields Characteristics of an air and ground temperature profile Permafrost distribution terminology Energy processes that contribute to the ground temperature Deep borehole temperature profiles in Yukon Territory Snow accumulation model Temperature profile of ground and atmosphere near the surface Effect of the spatial variable transformation on the layered domain A layered soil column : Simple three component sample soil volume Composite sample soil volume Effect of geometry on thermal conductivity Latent heat, heat capacity and thermal conductivity of water Summary of thermal conductivity relationships for dry snow Summary of heat capacity relationships for dry snow Staggered grid Effect of Shape Parameter B Bulk thermal properties with phase change Numerical solution versus analytical solution of test problems The Stefan Problem Numerical solution versus analytical solution of the Stefan Problem Effect of the spatial variable transformation on the discrete layered domain Global mass and global energy conservation Instantaneous and envelope temperature profiles Effect of soil type on temperature profiles Effect of winter snow thickness on temperature profiles Effect of length of winter season on temperature profiles Effect of surface organic layer on temperature profiles 1 Effect of surface organic layer on temperature profiles 2 Hypothetical ground section and 32 soil profiles Solution temperature profiles and histogram of active layer depth Solution temperature profiles and histogram of winter frost penetration Solution temperature profiles and histogram of active layer depth Solution temperature profiles and histogram of active layer depth  vi  1 3 4 6 11 13 14 16 17 19 20 23 26 29 30 31 33 38 40 41 42 43 46 49 50 51 52 53 54 55 57 58 59 59  4.12 4.13 4.14 4.15 4.15 4.17 4.18 4.19 4.20 4.21 4.22 4.23 4.24 4.25 4.26 4.27 4.28 4.29 4.30 4.31 4.32 4.33 4.34 4.35 4.36  Boreholes and weather stations in the Yukon Territory and Mackenzie Valley Permafrost occurence in Yukon Territory and Mackenzie Valley Boreholes and weather stations in the Mackenzie Delta Inuvik weather station records Tuktoyaktuk weather station records Soil depth in the Mackenzie Delta Soil texture in the Mackenzie Delta Model soil columns and boundary conditions for the Mackenzie Delta Solution profiles for the Mackenzie Delta regional model Comparison of borehole and solution profiles for the Mackenzie Delta region Mackenzie Valley region overview Norman Wells weather station records Soil depth in the Mackenzie Valley Soil texture in the Mackenzie Valley Model soil columns and boundary conditions for the Mackenzie Valley Solution profiles for the Mackenzie Valley model Measured borehole and modeled profiesl for the Mackenzie Valley Southern Yukon Territory overview Mayo weather station records Whitehorse weather station records Soil depth in the southern Yukon Soil texture in the southern Yukon Model soil columns and boundary conditions for the southern Yukon Solution profiles for the southern Yukon Measured borehole and model solution profiles for the southern Yukon  vii  60 61 62 63 64 65 65 66 67 68 69 70 71 71 72 73 74 75 76 77 78 78 79 80 81  1 INTRODUCTION 1.1 Context and Goals This thesis concerns applications of numerical modeling of energy balance to problems of ground temperatures and permafrost. Permafrost is a thermal ground condition ubiquitous in polar latitudes and high elevations, underlying about one half of Canada (Brown, 1970). The thickness and extent of permafrost is a factor in vegetation distribution, surface and groundwater hydrology, atmospheric carbon exchange, civil and resource engineering, and moisture exchange between the earth and atmosphere (French and Heginbottom, 1983). Additionally frozen ground is an archive of information about paleoclimate and ecology (Lachenbruch and Marshall, 1986). Figure 1.1 identifies several fields of science and technology that are affected by the condition of frozen ground.  / Frozen Ground  \ Figure 1.1  Drainage  Boreal/Tundra Ecology  Global Carbon Balance  Ground Moisture  Atmospheric Moisture  Arctic Climate  Engineering  Settlements  Resource Extraction  Relationship of frozen ground to other fields  It has been recognized that arctic climate responds sensitively to global changes in climate. Deep boreholes have recorded a warming of 2 to 4°C since the Little Ice Age, and parts of the Mackenzie Basin of Northwest Territories have warmed 1.5°C in the last century (Burn, 1998). In understanding the effects of climate change as a problem of physics the key variable of interest is temperature, both of the air and the ground. Understanding the condition of permafrost under changing climate regimes is at its core a problem of understanding energy balance and temperature fields. Modeling of ground temperatures is not a new area of study. Analytical solutions to problems of energy balance in various geometries and media were developed in the 19th century and continue to be a useful approach to exploring simple problems of freezing and thawing (Carslaw and Jaeger, 1959; Lunardini, 1996). Numerical strategies to solve more complex problems have advanced in tandem with increasing computational power.  1  Existing models of ground temperatures and permafrost generally fall into three categories: calibrated physical models of temperature diffusion in one or more dimensions (e.g. Nakano and Brown, 1972), areal distribution models based on climate statistics (Nelson and Outcalt, 1987), and semi empirical models which relate climate statistics to ground temperatures (Smith and Riseborough, 1996). Burn (1998) compares these three approaches and outlines some strengths and weaknesses of each style. The work of this thesis falls into the first category. However, it is the statistical and empirical approaches which have gained popularity in studies of climate change because they are more easily applied to the continental and hemispheric scale required by general circulation models (Nelson and Outcalt, 1987; Hinzman et al. 1998). The weakness of distributive models is that they do not provide any information on temperatures in the ground, nor are they able to illuminate the effects of transient changes in surface climate. It is only physically based numerical models which are able to capture the real world variation in boundary conditions and the changing temperatures within the permafrost at depth. The goals which motivate this thesis are to explore numerical modeling strategies that: • Have a basis in physics (as opposed to empiricism). • Have a small or manageable set of input data requirements that is not application specific. • Provide full solution of the ground temperature field in both space and time under transient surface forcing. The following quotation establishes the spirit which guides this work: "Reliable predictions of permafrost change are handicapped, at present, by a lack of accurate information on the thermal regime and physical characteristics of permafrost over significant areas. This cannot be totally compensated by good thermal models of permafrost change, but the margin of speculation can be reduced." (Lunardini, 1996) However the following caveat should not be ignored: "// would be dangerous, though attractive, to assume that computer modeling may be substituted for field investigations." (Burn, 1998) 1.2  Terminology  As with any branch of research the field of permafrost and ground temperatures has a rich but occasionally non-intuitive vocabulary. It is useful to define some terms. Figure 1.2 illustrates several characteristics of an air and ground temperature profile and explains the common acronyms. For the purpose of weather data collection, air temperature is usually measured at a height of one to two meters above the ground surface. Mean annual air temperature ( M A A T ) is the annual mean of temperature measured at this level. In moisture bearing ground that experiences seasonal freezing and thawing of an active layer, the M A G T profile will exhibit a curvature in the near surface due to the contrast in thermal properties of water in its liquid and ice states. The thermal offset is the difference in temperature between the M A S T at the top of the profile and the temperature at the inflection of the profile. 2  Mean Annual Air Temperature, (MAAT)  r i  \>\,  *— Mean Annual Surface Temperature (MAST)  7  y  Seasonal Surface Temperature Amplitude (SSTA) Thermal Offset  Mean Annual Ground Temperature (MAGT)  Depth  Figure 1.2  Characteristics of an air and ground temperature profile.  Figure 1.3  Permafrost distribution terminology.  3  Permafrost is defined as ground which remains below zero degrees Celsius for a minimum of two consecutive years. The active layer is the maximum thaw depth over permafrost during the thaw season (Anderson and Anderlund, 1973). The inflection of the M A G T profile due to thermal offset in the upper moisture bearing soil occurs at the base of the active layer. Because the coldest value of M A G T occurs at the base of the active layer it is possible to have stable permafrost even when the M A S T is above freezing (Smith and Riseborough, 1996). The terms sporadic, discontinuous and continuous are used to describe the areal extent of frozen ground. Where permafrost underlies less than 35% of the ground surface the term sporadic is applied; continuous permafrost requires greater than 90% coverage, while discontinuous permafrost is the intermediate condition (Dyke and Brooks, 2000). The transition from discontinuous to continuous permafrost extent corresponds roughly to the location of the —5 °C isotherm of M A A T (Burn, 1998). Figure 1.3 is a cartoon illustrating how frozen ground might be distributed in a transect from the continuous to sporadic permafrost zones. The frost front is the location of the zero degree Celsius isotherm in the temperature profile. In permafrost with an active layer there will be a frost front at the bottom of the permafrost that migrates exceedingly slowly and an upper frost front that migrates with seasonal temperature variations.  4  2 PHYSICS OF GROUND  TEMPERATURES  In this chapter I develop the analytical equations that will be used to solve for the ground temperature profile. The model is developed conceptually before it is expressed in discrete form for numeric coding. 2.1 E n e r g y Balance and H e a t Flow  This section discusses the equations that describe energy flow appropriate to the ground temperature problem. The energy balance equations are developed in full detail from global to local form. This development has been a common source of error in past thermal modeling literature. 2.1.1  E n e r g y Processes in the G r o u n d  The temperature profile in the soil column is determined by the physics of energy flow in response to changing boundary conditions at the surface and depth. Generally the scales of lateral and vertical heterogeneity differ enough that the subsurface can be treated as a horizontally layered system and treated analytically in only one (vertical) dimension. This one dimensional approach would be inappropriate in locations where there is significant lateral energy flow (e.g. aquifers, borders of lakes/water bodies, buried utilities, locations of extreme contrast in surface conditions over short distance). In the case of frozen ground the layers of the soil column are not only distinguished by structural and physical contrasts, but by temperature contrasts also. In particular it is possible to identify layers as frozen or unfrozen, with the border between the two corresponding roughly with the frost point isotherm. Figure 2.1 illustrates the typical layers in the soil column that can be identified on thermal and structural bases. In addition the chart indicates the thermodynamic processes active in each layer. Note that the soil column is static and does not include any biogenic or radiogenic heat sources. Conductive energy flow operates throughout the soil column, and convective heat transfer is possible anywhere there is a fluid component to the column (groundwater or air). At the margins between frozen and unfrozen layers the latent energy of the phase change of pore water contributes to the overall energy balance. Radiation and vapour processes are limited to the upper boundary of the column, and the contribution of convection (mass transfer) to heat transfer in the interior of the soil column is minimal compared to conduction. In coarse grain or organic soils non-conductive effects can occur but in practical long term measurements a purely conductive description of energy flow is adequate. In cold permafrost, heat transfer is exclusively by conduction (Lachenbruch and Marshall, 1986; Burn 1998).  5  L o w e r a t m o s p h e r e Contact layer U p p e r soil layer 1 Active layer  U p p e r soil layer 2 Frozen fue utih r tnmhililr  Frozen  -nofreewafer  D e p t h of a n n u a l temp, oscillations  Ftoztn  •• -nofreewater  Permafrost thickness  Frozen frit u alt r iittnhible  Unfrozen lithosphere  Figure 2.1. Energy processes that contribute to the ground temperature.  2.1.2 Equation Development for Energy Balance The global form of the energy balance is dU  ir  dK  +  ^  = Q  ^ +  p  '  „  '  ft  _  (">  where U is the total internal energy, K is kinetic energy, Q is the total rate of heating of the system, and P is rate of working. For the system of interest, kinetics K (e.g. strain heat), internal heat sources (e.g. radiogenic or biogenic) and applied forces (work done) P do not contribute to the energy balance. Thus the global balance equation simplifies to,  6  (2.2) Define the system as a continuum body of volume V bounded by a surface S. The heating of the system due to energy transferred through the boundaries is expressed in integral form as  Q = -J^q n d r  = -J^^d\,  2  k  k  where n is the surface normal vector and q is heat flux with units of W m conductive heat flux q , Fourier's law applies: k  k  (2.3) - 2  . For  k  dT  » =- * * : •  (2  -  4)  Therefore,  J  dx  v  \  k  dx /  (2.5)  k  The total internal energy U is an extensive property of the material volume V that can be expressed as a volume integral of the intensive time-dependent field property (pu), where p is density with units of k g m ~ and u is specific internal energy which has units ofJkg" , 3  1  U = I pud r. 3  Jv  Specific heat (at constant pressure) is defined as  Integrating Equation (2.7) between two arbitrary reference temperatures gives,  7  (2.6)  A convenient choice for the lower reference temperature is absolute zero (To = 0 K ) where the specific internal energy u is also zero. Equation (2.8) simplifies, 0  u (T) = I c(T')dT', Jo  (2.9)  where the integration limit T is a temperature field, T = T(xm,t). By the Reynold's transport theorem the time rate of change of the extensive property U(t) can be expressed as,  dU dt  'v(t)  Applying the chain rule and collecting terms,  Recognize that the second term of (2.11) vanishes by the restriction of mass conservation  dU  *=/  f  (du  du \ ,  '(« "&;K'+  v(1)  (2  -  . 12)  A n approximation of the partial time derivative of u is  du ~di  u(x ,t  + dt) — u(x ,t) Jt "  m  m  (2.13)  Expressing Equation (2.9) at the times t + dt and t,  T(x ,t+dt)  r  u(x ,t m  + dt)= /  Jo  u(x ,t)= m  /  Jo  The lower integration limits are the same so,  m  c{T')dT\  c(T')dT'.  (2.14)  (2.15)  u(x ,t m  + dt)-u(x ,t)  = /  m  c(T')dT.  (2.16)  JT(x ,t) m  The integration interval is very small, (order of (dT / dt)dt) and the change in the value of the integrand over that interval is even smaller (order of (dc/dT)(dT/dt)dt). The value of c(T) is essentially constant over the integration interval.  /  c(T')dT' * c(T(x ,t))  /  m  JT(x,„.,t)  dT',  (2.17)  JT(x ,t) m  t + dt) - u(x ,t) m  , , T(x ,t = c(T(x ,t))  + dt) - T(x ,t)  m  m  m  _ .  /oio\  (2.18)  The partial time derivative of u can be expressed in terms of the partial time derivative of the temperature field T,  ^=c(T(W))^,  (2.19)  and by a similar development of the spatial partial derivatives of u,  £  Substituting  (2.19)  and  (2.20)  lf=[ d t  =  into Equation  Jv(t)  (2.20)  « ( n w ) j £ .  Pc(T)?£  (2.12)  ot  + v pc(T)£-a*r.  (2.21)  dx  k  k  The local form of the energy balance can be derived by substituting  (2.5)  arid  (2.21)  into  (2.2).  / pc (T)— +v c (T)—- d\ = / — Jv(t) dt dx J ox \ kP  J  k  v  k  (K—-)d r dx J 3  k  (2.22)  9  For the above equation to apply to any arbitrary volume V(t) the integrands must be equal:  The product of p and c(T) is the volumetric heat capacity C(T) with units of J m ~ K . Additionally for the layered earth the governing equation can be simplified to one vertical dimension 3  - 1  Equation (2.24) is the form of the governing equation most easily solved in discrete form. 2.2 Boundary Conditions on the Ground Temperature  Profile  2.2.1 Qualitative Description and Models Lower Boundary Condition Geothermal Heating In principle, the one dimensional domain called the soil column can extend to an arbitrary depth. However it is convenient to choose the bottom boundary of the domain as a depth where the vertical temperature gradient is constant. The geothermal gradient is a consequence of energy processes deep within the earth and is relatively constant in time and space, except in areas of volcanism or hydrothermal activity. A constant heat flux or geothermal gradient is a standard choice for the lower boundary condition in ground temperature modeling (e.g. Zhang et al., 1996; Osterkamp and Romanovsky, 1999). Deep borehole temperature profiles provide a direct measure of the geothermal gradient. Figure 2.2 shows deep borehole temperature profiles for a number of locations in Yukon Territory (Taylor et al., 1982). The temperature gradient at depth for these holes falls generally within the range of 0.02 — 0.035 K m , as indicated by the dashed lines. - 1  10  Temperature, °C  Figure 2.2.  Deep borehole temperature profiles in Yukon Territory  Steeper gradients of over 0.055 K m have also been observed in deep boreholes in the Mackenzie Valley near Norman Wells. - 1  Upper Boundary Condition Four microclimatic variables control the presence of permafrost in the discontinuous zone, for example south central Yukon. These are the depth of snow cover, the thickness of the organic soil layer, the soil moisture content and the vegetation canopy. In this region the presence of permafrost is more dependent on these physical variables than on climatic conditions (Burn, 1998). However in the continuous permafrost zone, for example Northern Alaska, where frozen ground is ubiquitous it is air temperature, rather than snow or soil moisture that controls the ground temperatures (Zhang and Stamnes, 1998). For subarctic permafrost, the cycle of seasonal snow is important in two ways. First, snow is a good insulator and thus reduces heat loss from the ground during winter and delays the penetration of springtime warmth. Second the thickness of the snow distances the ground from the atmosphere, acting as a filter on the phase and amplitude of the propagation of seasonal air temperatures into the ground (Smith, 1975). Along a transect of measured ground temperature profiles in central Alaska, systematic warming of the permafrost is not correlated with a systematic change in mean annual air temperature but rather thicker accumulations of snow during the observation interval 11  (Osterkamp and Romanovsky, 1999). Thus there are two processes controlling the energy exchange at the top surface of the soil column: energy transfer by conduction and/or radiation through the atmosphere-ground interface, and the effect of accretion or ablation of material at the upper surface. The location of the upper boundary of the soil column is tricky to define. The surface condition changes seasonally with the accumulation and ablation of snow. Longer term variations may also occur with changing vegetation and surface hydrology. Both the material properties and the thicknesses of the surface snow and vegetation layers can change. Taking the lower boundary as a constant reference depth, a change in thickness of a surface layer corresponds to a change in total height of the soil column. The highest frequency changes in surface variability occur in the seasonal snow layer. Constructing an upper boundary condition for the one-dimensional model of ground temperature first requires a model of the accumulation and ablation of this seasonal layer. Seasonal Snow Generally the snow depth builds slowly from autumn through mid-winter. There can be a period of little change or even slight deflation due to densification through mid to late winter before a rapid spring melt. Four dates define this pattern: the first day of permanent seasonal snow ( P i ) , the date ( P ) that maximum snow depth ( i T ) is reached, the date of the onset of spring melting (P3) and the last day of permanent seasonal snow (P4), within the annual cycle ( P ) . Goodrich (1982) used a simple linear increase in snow depth to a midwinter maximum, thereafter remaining constant until spring melt. Based on records from Barrow, A K , Zhang et al. (1996) suggested a description of snowdepth (Hs(t)) assembled piecewise from simple functions. Both models can be expressed in the following form, 2  H (0 <t s  <Pi) = 0  d  H {Pi  m a x  (2.25a)  < U < P ) = AH  s  2  (2.25b)  max  H {P2 <td< Pz) = H S  (2.25c)  max  Hs(P <U<P ) 3  #s(P  = (l-B)H  4  4  (2.25d)  max  < U < P) = 0,  (2.25e)  where,  (2.26a)  (2.26b) For a linear snowpack buildup and melt the exponents are ni = n = 1. The empirical model determined by Zhang uses values of ni = 1/2, n = 3/2 and P = P 3 . The snow depth model is illustrated in Figure 2.3. 2  2  12  2  0  50  f00  150 200 250 Days after July 1st  300  350  Figure 2.3. Snow accumulation model. The solid line is the Zhang model with n i = 1/2, ri2 — 3/2; the dashed line is linear with n\ = rti = 1  In the Mackenzie Valley a similar pattern of snow buildup is observed with the switch from rain to snow occurring by November and an increase in snowfall until the mid-winter arctic high strengthens to block incoming humid systems. While snowfall is diminished in mid-winter, snowdepth nonetheless increases as thaw events are rare and maximum snow depth is reached in/around March. The breakdown of the arctic high results in rapid spring ablation (Dyke, 2000). The same pattern of snow buildup is also observed in Norway (Seppala, 1982).  13  Temperature  Mean Annual Air Temperature, (MAAT)  Mean Annual Surface Temperature (MAST)  Seasonal Surface Temperature Amplitude (SSTA)  Depth of Zero Annual Amplitude  Mean Annual Ground Temperature (MAGT)  Depth  T  Figure 2.4. Temperature profile of ground and atmosphere near the surface (After Lunardini, 1981)  The temperature measured at the soil surface in field studies is usually greater than the temperature in the atmosphere only a few centimeters above the ground surface. The relationship between the two temperature measurements is complex. It is atmospheric temperature that is most commonly measured in the field and that is predicted by climate modeling. What is needed is some way to relate the atmospheric temperature to the ground surface temperature T 5 . Surface energy balance modeling attempts to book-keep all the energy contributions to a contact layer between the atmosphere and the ground. A substantial number of variables and parameters are involved in tracking the conductive, convective and radiative processes that operate in this layer (Lunardini, 1981; Hinzman et al., 1998). Such a complex surface model is a challenge unto itself and not a feasible option for determining the upper boundary condition of a more general ground temperature model. While sound in applying principles of thermal physics, surface energy balance models cannot be extended beyond the site scale because of the limited database describing microclimate vegetation and terrain.  14  Zhang et al. (1996) measured ground surface and air temperatures at Barrow, A K and determined the following empirical relationship for that location:  T s  ~ { 3.4 + 0.89T  #5 = 0.  A  ( 2 - 2 7 )  A tabulation of differences between mean annual air temperatures and mean annual ground temperatures from weather stations in northern Canada shows ground temperatures from 1 to 6 °C warmer than air temperatures. There is no correlation with thickness of permafrost or M A A T (Lunardini, 1981). 2.2.2 Boundary Conditions Applied to the Governing Equation The boundary conditions on the vertical soil column which extends from an elevation of z — 0 at its base to z = H(t) at its top can be expressed as  T  =T (t) s  (2.29)  The values of the geothermal gradient G, the surface temperature Ts and the time rate of change of the surface layer thickness dHs(t)/dt and are known for all time, as are the thermal properties of the soil. The boundary conditions are applied at the levels ZB — 0 and zs = H(t). As was described, the ZB = 0 level is chosen as a depth where the geothermal temperature gradient is constant. However if the upper domain boundary is defined as the interface between the earth and the atmosphere, it moves up and down with respect to the fixed ZB = 0 depth as surface layers accumulate and erode, thus zs = H{t). Such a moving boundary condition is tricky to solve numerically. It is possible to express the governing equations in terms of a new spatial variable 7 that adjusts in time such that the lower and upper boundaries correspond to constant values of 7. Marshall (1996) describes this grid transformation process in detail, with references to similar applications in glaciological and atmospheric studies. Let z = f(j, H(t)), where H(t) is the height of the vertical domain of the soil column. The form of the function / is arbitrary (as long as it does not compromise the accuracy of the numerical scheme). The simplest form for / is a linear scaling, however it is also possible to choose a form that when discretized will concentrate grid cells where the solution 15  changes most rapidly. In the case of the ground temperature problem temperature varies most near the surface boundary. Expressing the governing equation (2.24) in terms of the new variable 7 requires new partial derivatives which are obtained by applying the chain rule:  d  d_  d'j  dz  dz dj dz d dt ~ dt  (2.31) (2.32)  Substituting (2.31) and (2.32) into (2.24) produces,  dt  dt dz dj  dz d^ \  dz  (2.33)  dj  The tilde overscripts serve as a reminder that the governing equation is now being expressed in terms of a new space variable 7. The boundary conditions on this form of the governing equation are now applied at constant values of 7 and the problem of moving energy within the system involves an internal velocity field.  z = H(0)  z = 0 Time (years)  Figure 2.5.  time (years)  Effect of the spatial variable transformation on the layered domain  Figure 2.5 illustrates how the transformation of the vertical coordinate results in a domain of constant size but with variable property fields.  16  2.3 The Thermodynamic Material Properties of the Ground In this section I further develop the equations for the material properties. 2.3.1 The Layered Soil Column The soil column is a layered system. The simplest model is a seasonal snow layer above a homogeneous halfspace. Adding layers improves the depiction of reality but increases model complexity. The limit of extending the model to include an infinite number of different layers would be the uniformly varying subsurface. Limiting the model description to five distinct layers or less seems to be an adequate compromise. The soil column can be developed in situ by the mechanical and chemical weathering of bedrock into mineral soil. Alternatively sedimentary layers can be deposited and built up. Vegetation will add organic material to the upper horizons of the soil column. Figure 2.6 illustrates a hypothetical four layer soil column.  H  n  Snow  Hs(t)  Peat or Mud  H  Sediment or Weathered Rock  H  H(t)  Deep Sediment or Bedrock  Figure 2.6.  A layered soil column  17  a  b  2.3.2 The Constituent Mixture of Materials  Soil can be considered as a constituent mixture with solid, gas and liquid components. The physical and thermodynamic properties of a sample volume of soil will depend on the individual properties of the constituents, as well as the geometry and proportions of their combination. The local continuum equation of energy balance (2.24) assumes spatially continuous material property fields, thus it is necessary to determine the bulk values for a representative soil volume made of various constituents. The Components of a Soil Mixture  Constituent  Density (kgnf ) 3  H e a t  C a p a c i t y  <MWni K-') 3  T h e r m a l  Conductivity  (Wni'K- ) 1  Quartz  2660  Other typical minerals  2650  2.0  2.9  Organic soil matter  1300  2.5  0.25  Water (liquid)  1000  4.2  0.57  Ice  920  1.9  2.2  Air  1.25  0.00125  0.025  2.0  8.8  "  Table 2.1. Properties of common soil constituents (Hillel, 1980; Kay and Goit, 1975)  Table 2.1 lists the values of the material properties of the most common soil constituents. Several contrasts and similarities in value are worth noting. It is only necessary to differentiate between types of mineral grains if the soil has a high quartz content, as this will result in a higher value of thermal conductivity than typical of other minerals. Organic material has a low density and low thermal conductivity compared to minerals. The thermal properties of water have very different values in the liquid and solid phases. The low thermal conductivities of water and air are also important. Porosity and Saturation  The total volume dV of a soil sample comprises solid and void components, and 8V^ \ Porosity d> is the ratio of void volume to total volume in a sample. Figure 2.7 shows a simple three part mixture decomposed into its fractional constituent volumes. V  18  Figure 2.7.  Simple three component sample soil volume  dV  = SV  + SV  {S)  (2.34)  (V)  Voids (pores) can contain liquids or be 'empty' (gas filled). Saturation, f is the ratio of soil water filled volume SV^ to total void volume. The unsaturated void volume 8V^ ^ is generally assumed to be gas filled A  (V)  6V  =  y(F)  8  +  y{A)  8  (  2  _  3  6  )  The preceding simple mixture example assumes that the sample volume dV is large enough to represent the bulk properties of the soil, yet small enough to apply to a local 'volume element' in the continuum sense. A n additional level of complication is introduced when the sample volume dV is large enough to encompass the macroscopic heterogeneity of the structured soil column. This development leads to the discrete formulation and numerics, recognizing that the sample volume dV is analogous to the discretization interval or the numerical cell size. Let the sample dV be a composite of two differing constituent mixtures, denoted a and b such that dV = SV + SV . a  19  0  Figure 2.8.  Composite sample soil volume  Each subvolume has a distinct porosity and saturation indicated by the subscripts a and b. Bulk values for the porosity and saturation of the whole sample volume dV can be calculated  v  dV __ 6V  _ SV^  ~ 6VW  ~  (F)  f  dV  dV  + SV[  F)  6V™  _ f <p 6V a  ~  a  a  '  K  + f <j> 6V b  4>dV  b  b  •  -  [  b  )  Mixing Equations for Density, Heat Capacity, and Thermal Conductivity In the energy balance equation the material property fields required are those of density p(z,t), heat capacity C(z,t), and thermal conductivity K(z,i), each of which depends on the porosity and saturation of the mixture. Density The bulk density p^ is the mass of a unit volume of the soil mixture, including all constituents. For the simple three constituent mixture as depicted in Figure 2.7 the expression is p( ) = B  -4>)  +  pWftj, + ^ ) ( 1 - /)</>.  (2.40)  For a volume that contains several different types of solid grains, or one that spans a compositional discontinuity, the expression is slightly more complex, however some simplifications can be assumed. Having several different types of fluids within the total volume is a very unlikely possibility, thus it is reasonable to assume in general that p( ^ and p^ ^ F  20  A  are single constant values for water and air, and it is also reasonable to assume that the mass of air contained in unsaturated pore volume is negligible. Hence,  " = £^U (i-«+,. <B,  S )  Soil  Bulk p  IB1  ( F )  /^  (2.41)  Density (kgrn' ) 3  Peat  1000  Mud  1000 - 1300  Clay, silt  1400 - 2000  Sand, gravel  1700 - 2300  Till  2000 - 2400  Table 2.2 Typical bulk densities and porosities of common soil types (Andersland and Anderson, 1978)  Heat  Capacity  The volumetric heat capacity, C ^ - * is the change in heat content per unit bulk volume. Since soil is a mixture, the bulk value of for a sample volume (dV) depends on the quantities and types of constituents. The bulk heat capacity is calculated through a volume weighted summation of constituent contributions. The heat capacity of air is three orders of magnitude smaller than that of other common soil constituents and can be ignored. The heat capacity of the fluid fraction is temperature dependent, both sensibly and through phase change. Most soil minerals have similar heat capacities so it is not necessary to discriminate between them, however organic solids in soil can have an impressive effect on the thermal properties of the soil (deVries, 1975; Hillel, 1980; Lunardini, 1981). 3  For a simple three-part mixture the expression is  (2.42) If the solid fraction is composed of a variety of grain types, identified by different values of the subscript /?, which each occupy a volume 8Vp of the total sample dV, the expression becomes more complicated: 21  ^  (  B  )  = { i : ^  5  ^ } + / ^  )  (  i  ?  )  '  ^  where  Y,W  { S) p  (2.44)  = {l-<p)dV.  P  If the total volume dV is composed of different subvolumes identified by the subscript a each with a different three-part constituent mixture, the heat capacity of the composite mixture is  C  ( B )  =  £  _  ^  )  C  ( S )  +  ^  F  J  j  .  (  2  .  4  5  )  Conflating the two expressions (2.43) and (2.45) gives complete generality in calculating a bulk value of heat capacity for a complicated mixture volume. Thermal  Conductivity  The thermal conductivity of a mixture is dependent on the geometry of the constituents. The two end member cases are series conduction and parallel conduction and the relations for calculating the conductivity in these cases are simple. The weighted geometric mean is a relation for the intermediate case that permits the most generality about the geometry of the mixture and can be used with good accuracy for complicated systems with any number of constituents (Lunardini, 1981). Thus for calculating the conductivity of a simple mixture the geometric mean applies. For the compound case - that is when the sample volume includes an interface between two different units, there is a known structure which in the case of the layered soil column with one dimensional heat flow means that the series combination applies to the bulk case. Figure 2.9 illustrates the effect of geometry on the calculated value of conductivity for a saturated soil.  22  pure rock  1 j  1 s  —arithmetic \  \  1  1 \ .  harmonic X.  1  \ .  i  I" —  - - | 1 L  _ .  pure water 0.6  1 ...1  —•—  Porosity  direction of flux  direction of flux geometric  Figure 2.9.  Effect of geometry on thermal conductivity  The expressions of bulk conductivity for parallel and series conduction are given below: Kallel  = ^ (1 -  B) _ series (  (S)  / a ^  4>)  + K^ft  + K^ \l A  +), j ± _ , (1 - m  Ki^)  1  TJ-/  c \ "T  Ki ) A  -  f)4>  (2.46) (2.47)  The geometric mean conductivity is  (2.48) Even in a saturated ice-rich soil there will in reality be small air bubbles and discontinuities in the ice that will conspire to reduce the bulk conductivity of the soil. The real world limit of thermal conductivity of a soil as <f> —> 1 is less than that of pure ice. Equation (2.48) above can therefore be expected to overestimate the conductivity by a small amount (Slusarchuk and Watson, 1975). For an arbitrary number of sub-layers within the sample volume, denoted by the subscript a, the bulk conductivity is 23  #< >  (2.49)  B  =  d  V  where  dV  (2.50) a  2.3.3 The Phase Change of Soil Water Qualitative Description The importance of pore water in soil cannot be understated. Moisture content is a controlling factor in the mechanical and thermal properties of soil. Phase change of soil water at the freezing point has a dramatic effect on soil behaviour as it introduces a step change in thermophysical characteristics with only slight change in temperature. The thermodynamics of pure water are complicated, but the soil environment rarely holds pure water in pore space. A parsimonious treatment of soil water thermodynamics would be excessively complex and ultimately futile as full details of the microscale variability of soil geometry and chemistry at any given location are impossible to know (Anderslund and Anderson, 1978; Hillel, 1980). To capture the phase dependence of the material property fields used in the governing equation (2.24), it is necessary to develop temperature dependent functions for the volumetric heat capacity and thermal conductivity of the pore water constituents in the mixing equation. The heat capacities and thermal conductivities of the ice/water constituents (previously C^ and K^) will be denoted C*(T) and K^T). When the latent heat of phase change is included in a temperature dependent heat capacity this quantity is often referred to as the apparent heat capacity. Building an apparent heat capacity for the soil provides a strategy for dealing with the energetics of phase change without having to directly track a frost front as in Goodrich (1979) or Zhang et al. (1996). The apparent heat capacity can be built in a variety of ways, drawing on empirical or theoretical insight to find a formulation that captures the important temperature dependent behaviour of soil heat capacity. Hinzman et al. (1998) and Anderson and Anderslund (1978) present two different approaches, both of which require knowledge of the soil type and water content. I present a different approach which requires knowledge of the soil porosity, saturation, and dry solid properties. F)  Thermal conductivity and density of water also have different values in the liquid and solid phases. For these properties the temperature dependence manifests simply as a step change occurring at the frost temperature.  24  Equation Development  Equation (2.9) relates the specific internal energy (more accurately the enthalpy) of a volume to its temperature T and specific heat c. The change in internal energy with changing temperature is described as the sensible heat process  u(T)= f  Jo  c(T')dT'.  (2.51)  However, changing the phase of a material will result in a change in internal energy accompanied by no change in temperature - a process labelled insensible. The specific heat of a material is not defined at the temperature of phase change Tf (fusion temperature or frost point). The change in internal energy at the fusion temperature is called the latent heat. A liquid will have more internal energy than a solid, thus the freezing process liberates latent heat while the thawing process absorbs it.  If there is a phase change at temperature Tf < T then Equation (2.51) can be modified:  u(T) = f ' c(T')dT' + L+ j J  Jo  c(T')dT',  (2.52)  Tf  where L is the latent heat of fusion with units of J k g . This equation can be re-expressed in volumetric terms: - 1  pu(T)= [ ' C(T')dT'+pL+ Jo  [  C(T')dT'.  (2.53)  JT  }  One approach to building an apparent heat capacity for water that includes the insensible energy contribution of phase change is to consider the temperature dependence in C* as having two contributions: a step change between the solid and liquid properties and a latent heat spike occurring at a single defined fusion temperature Tf.  25  350 300  %  B250  ^3.5  u CQ  g.200  f I  4.0  §3.0  150  £2.5  100  -20  a)  -10  10  0  -20  Temperature ( °C)  j  Temperature ( °C)  c  0 io 10 Temperature ( °C)  -ro  Figure 2.10. Latent heat, heat capacity and thermal conductivity of water as functions of temperature  These two contributions can be represented mathematically by the Heaviside function H(x — XQ) and the Dirac delta function 6(x — XQ) respectively. Recognize that the temperature dependence in volumetric heat capacity involves contrasts in both the specific heats and the densities of water and ice. Assuming that the phase change occurs with negligible volume change, the volumetric heat capacity of the ice-water system can be expressed  C*{T) = { %  c e  T  l> water  C*(T) = C (T)H(T  1  ice  f  l  [  (2.54)  T  >  1  f  —T) + pL8(T - T » + C (T)H(T water  - T) f  (2.55)  where the • subscript indicates a state between ice and water, L is the gravimetric latent heat with units J k g and p is an intermediate density, for example, - 1  P = -J{P  ice +  P  water)-  (2-56)  Defining C*(T) this way the following equation is equivalent to Equation (2.53)  pu(T)=  I C*(T')dT'. Jo  (2.57)  A temperature dependent thermal conductivity K(T) can be constructed in a similar manner:  26  K  ^ = { K1„ III',  <- ' 2  K  K+(T) = K  ice  + H(T-T )(K -K ). f  water  48  (2.58)  ice  In finer detail this strict step-change switch in pore water from liquid to solid phase is not observed. Above Tf the thermal conductivity of a soil is insensitive to temperature changes. Below Tf all pore water does not freeze instantaneously and the fourfold contrast between the conductivities of water and ice leads to a more gradual change in bulk conductivity. Especially for finer grained soils there can be significant amounts of supercooled unfrozen water, with up to 10% remaining liquid even down to 22K below the frost point (Penner, 1970). The consequence is that Equation (2.58) will tend to overestimate the thermal conductivity of the fluid constituent that would be observed in a real soil. Referring back to the mixing equations of Section 2.3.2 the fluid constituent can now be replaced with the temperature dependent functions  C  ( B )  (T)  (2.59) -l  6V  n  K^-^KATY^K^) / « < £ «  Ts(A)( ~f '>'1 1  a  >a  (2.60)  2.3.4 Thermal Properties of Snow The Importance of Snow Compared with most earth materials snow is an insulator. The effect of snow on the ground thermal regime is to limit the release of ground heat accumulated in summer, thus inhibiting the penetration of winter cold. The result is to raise the mean annual ground temperature ( M A G T ) above what it would be if there were no snowcover. Additionally the high surface albedo of snow affects the difference in mean air and ground temperatures. This relationship between snow cover and ground temperature has been confirmed by field studies (Burn, 1998) and by modeling (Goodrich 1982). In the discontinuous permafrost zone snow cover affects permafrost distribution. In the continuous zone, where permafrost exists everywhere, snow cover affects the temperature of the permafrost. In permafrost environments, there can be an additional insulative effect in that subarctic snow consists primarily of depth hoar. Depth hoar is a metamorphosed snow texture made up of large (up to 3 cm) prismatic grains that form in cold temperatures under steep temperature gradients. The crystals bond poorly, restricting conduction through the solid matrix. The thermal conductivity of depth hoar is less than that of unmetamorphosed snow such that a given snowpack thickness will result in a warmer ground surface temperature if the depth hoar fraction is higher (Zhang et al., 1996). Typically an arctic snowpack is thin, cold, dry, and low in density. 27  Thermal Conductivity There have been studies of the relationship between snow density and thermal conductivity published for over a century. Figure 2.11 shows the curves for functional relationships derived from fitting observations of various studies as presented in two summary papers (Mellor, 1964; Sturm et al., 1997). For contrast, the thermal conductivity of pure ice is 2.22 W m ^ K " (Andersland and Anderson, 1978). In Figure 2.11 the light dashed lines are listed in Mellor (1964), the results of 8 studies published between 1894 and 1958 from observations around the world. The functions themselves range from linear to 4th order polynomial fits. The heavy lines are the results of statistical analysis of a single set of 488 measurements taken over 12 years in high latitude studies (Sturm et al., 1997). The heavy dashed line is a logarithmic fit that gives best results at low densities: 1  K(p  < 0.6) =  IQ/ ' 2  6 5 0 P  ' " '" ,  I O  - : L  '  6 5 2  (2.61)  ).  The solid line is a piecewise continuous concatenation of a linear extrapolation and a quadratic fit which gives best results if densities vary in the full range from zero (air) to solid ice (0.917g c m ) : - 3  0.023 + Q.2Mp 0.138 - lM + 3.233/9 2  Psnow <  anow  K  snow  Psnow  1  i/ 1  t  •!/  1 1  j  1  /  !  " " " [ " "  1 1 1 1  t  Sturm etal., 1997 Abels, 1894; Goodrich. 1982 Mellor, 1964  -/ /•- i ii I/  ]  1 0.5\ O 0.4  psnow  i ;  !;> 1 / / ! //  u I  & o.e| *  /  • •'/  0.156  0.156 <  '  I/ :/// /0-M'\ —  ! // ' 1  /  I// /  'Iff  ////\/,  /  I f  f  3  ! !  1 1  1 1 1  — f--//c&P-\ i 0  Figure 2.11.  100  I  300 Density kg/m'  Summary of thermal conductivity relationships for dry snow.  28  (2.62)  Heat Capacity The snowpack can be thought of as a porous constitutive mixture in the same way that soil has been. In this case the solid matrix consists of ice grains, and for dry snow the pore space contains no water. For a solid fraction n a simple estimation of the bulk volumetric heat capacity is  C = nC  ice  + (1 - n)C .  (2.63)  air  The solid fraction is n = (1 — 4>) where 4> is porosity, and for dry snow can be estimated as  71 — P snow  I Pice*  This relationship is plotted as a solid line in Figure 2.12, along with a number of other published estimates of volumetric heat capacity for snow. Figure 2.12 demonstrates that the heat capacity of snow is a well constrained linear function of snow density. 1.2  0  1  I 100  200  300  Density,  400  500  kg/m 3  Figure 2.12. Summary of heat capacity relationships for dry snow. The light grey boxes indicate typical values listed in Andersland and Anderson (1978). The dark grey swath on the main plot is derived from the functional relationships given by Mellor (1964) for ice of different purities, the dashed line is a typical value of specific heat suggested by this source. The large dots are effective volumetric heat capacity of a snowpack containing various proportions of depth hoar following Zhang et al. (1996).  29  Together Equations (2.62) and (2.63) give a convenient way of determining the thermal properties of snow based on the single variable of snow density in the numerical model. However an important consideration is that in the very cold, dry environments typical of polar permafrost significant snow metamorphism occurs due to the strong temperature gradients of thin cold snowpacks. In these conditions the thermal properties of the snow layer will depend more on the structure of crystals and on the cohesion and strength of depth hoar layers rather than on the bulk density of the snow layer (Sturm and Johnson, 1992; Zhang et al., 1996)  30  3 NUMERICAL DEVELOPMENT 3.1 The Discrete Staggered G r i d As introduced in previous sections the model domain is a layered column extending from an elevation of ZB — 0 at depth in the earth to an elevation of zs = H(t) at the surface of the earth/bottom of the atmosphere. The surface elevation zs varies in time with the evolution of the surface boundary condition. A new spatial coordinate variable 7 was introduced with the properties that the upper and lower domain limits of 7 are constant values, i.e. 7 is scaled to z through the variable domain height H(t). W i t h the upper and lower values of 7 denned the whole domain can be sectioned into n discrete cells of equal 7-volume. The locations in 7 space of the centres of these regular cells correspond to migrating locations in z space. The temperature field is solved as a suite of temperatures applied at the cell centres. Fluxes occur at the cell faces which are located regularly midway between the cell centres in 7 space. This offset of half of a cell width between the field and the flux is reason for describing such a discretization as staggered.  Figure 3.1.  Staggered grid  3.2 The Governing Equations 3.2.1 Algebraic Manipulation The governing equation (2.33) to be solved numerically is  dt  dt dz dy  dz dy \ 31  dz dy  J  The tilde overscripts serve as a reminder of the special nature of the space variable, from this point however they will be assumed and dropped from the notation in the interest of reducing visual clutter. A transformation function z = f(y,H(t)^ relates the space variables. A simple form for the transformation function is an exponential scaling. Let H(t) be the height of the vertical domain and B is an arbitrary shape parameter:  H{t) U17 B B  dy d~ dz  (3.2) (3.3)  H(t) 7 InydH ~B~dJ  z  dt  (3.4)  The partial derivatives of the transformation function can be substituted into the governing equation, Equation (3.1)  dT  Cylny dHdT  (  m iiw^  c  dT ~dt  ' B \ \  B  V  d ('  = {w)  +  d (  dT\  ,  Q  (3  dT\  c  ,  5)  ylnydHdT  H) cd^\ d^)-^r^tW Kl  -  (3  6)  Thus the time derivative of the temperature has two contributions, a conduction term and an advection term,  ~di D T  (A)  ~dt  Cdy\ y\ny dH dT '~H~~dtdy  32  7  dy  (3.7a) (3.7b)  s* I 0.9 B = 1 5  i i  B=8 1  B=5  ~?~ ' r  / N  0.6  ts .N  "5  g o c  0.5  /  B=3  /  1/  1  I /  7  ^  /  /  V  -A /  /  / / i  /  /  /B=0.S/'  1/  I  I  / B=0.1  -i-  I 0.4  Figure 3.2.  0.5  0.6  normalized  J  0.7  0.8  0.9  Effect of the shape parameter B  The shape parameter B determines the magnitude of exponential stretching in the numerical domain. A value of B closer to zero results in a nearly linear non-dimensional scaling. High positive values of B result in more discrete cells concentrated in the upper part of the domain, however the resulting under-representation of the lower domain requires care to ensure the lower boundary condition is still applied properly. Negative values of B will concentrate cells in the lower part of the domain. 3.2.2 Finite Difference Discretization The finite difference approximations for differential equations derive from Taylor expansion of the partial derivatives in space. Details can be found in Lunardini (1981). Using centred differencing, the conductive and advective terms for the time derivative of the temperature field are  33  (3.8a) 0T< > A  Tti  7 l  ln  7 i  dH j T  - Tj. \  i+1  —H-lu\  =  x  2A  J"  7  ( 3  -  8 b )  Equations (3.8a) and (3.8b) are in a form suitable for solution in M A T L A B . 3.3 The Boundary Conditions 3.3.1 Lower Boundary Condition The lower boundary condition is a constant temperature gradient, dT/dz — G. The value of this constant G is the geothermal gradient which may vary somewhat laterally, but not with time:  B dT  dj dT G  dz-^  =  ^ )  =  =  HW  ^  7  —•  (3-10)  The governing equation expressed in discrete form for the lower domain boundary becomes  d  T  (BV  iO  1 f(^7)x |  7I  dT^  \n dHJ dTx  7l  7l  _  +  T -Tj  f  2  (  K  . ( d T s  (3.11a)  K  Substituting (3.10) into (3.11):  3T  {K  ( A )  dt!  B—'Y  {H*J) c;\ ^ -  =  ¥a  7 l  ln  7 l  Tl)  dH \ GHAj  -wK -^\-B^ -  =  7  iT2  .(3.12a)  \  - \-  + T2  B  34  TL  -  (3  12B)  3.3.2 Upper Boundary Condition The upper boundary condition on the temperature field is a known temperature. Tying this known temperature to an actual physical model of energy balance at the ground air interface is non-trivial. The details of such a model are too much for this work and will be ignored. Thus the upper boundary condition is chosen to be  T =  +\  T n  s  Rearranging (3.13) for T  dT^ dt „  (  n + 1  R  C ^ ^( K  .  (3.13)  and substituting into (3.8) where i = n :  B \\ ,  \HAjJ  + T n  s  - " ) - tt ^  2Ts  r  n  K  + ( * 7 ) „ - j ) r + (Jf7)«-lT»-i j ,  s  n  (3.14a) din  H  Tn ) 2A  dt I  T —i  (3.14b)  n  7  3.4 The Phase Change The contrast in the thermal properties of water between the liquid and solid phases has been previously described with equations involving the Heaviside and Dirac delta functions. Both these functions are inconvenient in a numerical scheme, thus I introduce simple functional approximations. The singular behaviour embodied by the Dirac function in Equation (2.55) is an unpleasant obstacle. W i t h the goal of making the equations more computationally tractable without excessive dilution of the physical reality, the Dirac function can be replaced with simple bell like function with unit area  (0  F(T) = I G(T)  T <Tf - \AT T - | A T / <T <T f  f  r > r  f  /  + f AT/  (3.15)  + | A r . /  The bell function G(T) should be simply evaluated, and well represented by a small number of points when discretized (broad peak). A temperature range A T Q , centred about the frost point T / , can be tuned to give the best results.  35  Defining  T__ = 2/ - 2 A T  (3.19a)  G  (3.19b)  T_ =  \AT  T+ = i> +  (3.19c)  G  (3.19d) 16 AT 2 Ai = Alb' A  0  =  (3.19e)  3  G  (3.19f)  the bell function G(T) is constructed as: 0 A (r-r__) Aj - ^ ( r - r / )  r < r__  2  0  G(T)  2  A (T-T++)  2  0  l 0  r__ <T<TT- <T <T+ T+ < T < 1+4.  (3.20)  The Heaviside function -H^T) is replaced by a smooth step function S(T) over the range AT . H  Taking  T-  =  T --AT  T  =  T + -AT  H  H+  f  (3.21a)  H  (3.21b)  1  f  H  9  (3.21c)  AT,H 2'  the step function 5(T) is constructed as: 0  T < THT - <T <T Tf <T < T H  i- Bo(r-r  1  J  H +  )  2  H+  T +<T H  36  f  .  (3.22)  The new functions S(T) and G(T) can be substituted into the appropriate equation.  C*(T) = C  + S(T)(C  K*(T) = K  + S(T)(K  ice  ice  water  - C. j + pLG(T)  (3.23)  - K^).  (3.24)  ce  water  It could be argued that the step change function S(T) would be better constructed to extend a width ATH to the cold side of Tf. This is because while there can often be supercooled liquid water present within soil pores at temperatures several degrees below Tf there is never 'superwarmed' ice within pores at temperatures above Tf. The step function S(T) as outlined in the equations above is anti-symmetric about Tf, but need not be. The limiting temperatures Tjj- and Tfj+ and the parabolic coefficient Bo could be modified to adjust the step curve to more closely match curves observed in real soils. As was shown in Chapter 2 the bulk thermal properties of a simple three constituent soil mixture with phase change are calculated as follows:  C \T)  = (1 - </>)C + f<f>C*(T),  (3.25)  K^\T)  = K ^ ^ K ^ T ) ' * ^  (3.26)  (B  (5)  1  - ™ .  Figure 3.3 shows a hypothetical temperature profile with two frost fronts and the corresponding temperature dependent bulk heat capacity and thermal conductivity calculated from Equations (3.25) and (3.26). The width of ATQ is 2 K as indicated by the shaded zone. It can be seen that where the temperature gradient is steepest the frost front is most sharply defined in the profile of volumetric heat capacity. 3.5 The Discrete Property Fields From The Layered Column The discrete equations (3.8), (3.12), (3.14) require a value of heat capacity appropriate to each discrete volume referenced by the index i and a value of thermal conductivity appropriate to each discrete volume referenced by the index i + | . Chapter 2 introduced the layered column of the domain and provided strategies for calculating bulk material properties of the layered soil mixture with a phase dependent water fraction. The numerical description of the layered column can be distilled into a suite of 6 parameters for each non-snow layer: layer thickness H , thermal conductivity, heat capacity, and density of the solid fraction C^f\ pi^\ porosity (j> and saturation f . The thermophysical properties of the soil solids are tied to a particular soil lithology, and thus are not independent. The density, thermal conductivity, and heat capacity of the snow layer are prescribed directly as bulk values: p , K ,C . a  a  3  S  a  S  Each discrete cell with volume AVi within the numerical domain has bulk properties determined by its location within the layered column. The mixing equations (2.49) and (2.50) developed in Chapter 2 are expressed in discrete form 37  Figure  3.3.  Bulk  thermal properties  with  phase  change.  (Note the  heat  c a p a c i t y is p l o t t e d o n a l o g a r i t h m i c scale)  C ( B )  ^  I =  ^(E*  V  «{( -^)  ^ ' ( ^ ^ V , ( £  w h e r e 6V  a  are s u b v o l u m e s w i t h i n  1  C  ...^  r  ^  +  ^« *(^)}^  *^^  (3-27)  C  ) ( 1  . . /  ) #  ..)  ',  (3.28)  AV;  Y,W  a  = AVi.  38  (3.29)  3.6 Testing 3.6.1 Testing the Diffusion Code  For uniform material constants and simple boundary conditions the diffusion equation has analytical solutions. I compare the results of the numerical model to the exact solution for two cases: a periodic steady-state surface temperature, and an instantaneous surface cooling of a uniform halfspace. The governing equation for thermal diffusion simplifies to  W  =  d*>  K  ( 3  -  3 0 )  where K = K/C is the thermal diffusivity of the material. The initial and boundary conditions for the two test cases, and the exact solutions are as summarized below. Periodic Surface Condition  The boundary condition applied at one end of a semi-infinite domain is  T {z = 0, t) = T + ATcos(ut). s  0  (3.31)  The exact solution is  T(z,t) =T + 0  ATexp  where AT is the amplitude of the periodic surface temperature, u> is the angular frequency and To is the longterm mean value of Ts(z = 0,t). Instantaneous Cooling  The initial condition of the semi-infinite domain is  T(z > 0,t = 0) = T .  (3.33)  0  The boundary conditions for t > 0 are  T{z = 0,t > 0) =  T, s  T(z -> oo,t > 0) = T . 0  39  (3.34a) (3.34b)  The exact solution is  T(z,t) = (T  -T )erfc  s  0  2\fhti  + T , 0  (3.35)  where erfc is the complementary error function.  1.5  2  2.5  Time, (years)  Figure 3.4. Numerical solution (symbols) versus analytical solution (solid lines) for the Instantaneous Surface Cooling problem (top) and the Periodic Steady State problem (bottom)  Figure 3.4 shows the comparison of numerical and analytical solutions. The solid lines are the exact solutions as determined by Equations (3.32) and (3.35). The symbols are the numerical solutions at discrete times and depths. 3.6.2 Testing the Phase Change A n established method for testing the accuracy of the numerical phase change technique is by comparison to the known solution to the Stefan Problem (Goodrich, 1978). The solidification of an initially liquid semi-infinite region is a non-linear problem with an exact solution known as the Neumann Solution to the Stefan Problem (Carslaw and 40  Jaeger, 1959). The surface of separation between the frozen and unfrozen domains is mobile and therefore a function of time. As indicated in Figure 3.5 the position of this surface is denoted x(tf).  '  \  \  \  \  solid  'dt  i i  moving frost front  \  t  o  %(t)  ]  Figure 3.5.  liquid  The Stefan Problem  The quantity is the thermal diffusion distance which is the characteristic lengthscale of the problem.The depth of the phase change front %(<) is proportional to the thermal diffusion distance, %(r.) oc y/Ki t. The constant of proportionality is expressed as 27V, thus ce  = 2NVr7~t.  (t)  x  (3.36)  The value of N is determined from the following transcendental equation which involves the thermal properties of water and ice, and the initial and boundary conditions of the desired problem. Note that in this case the assumption is that density does not vary over the phase change. Solving the Stefan problem with the density contrast requires only a slight increase in analytical complexity. (Carslaw and Jaeger, 1959; Turcotte and Schubert, 1982),  (T -T ) . . erf(iV) f  c  / A T  K \lZ:. J f T  _2  water  N  e  H  r  t  °) _ /pEjv)»  K  ice  e  rf ( c  LN^  u  e  —  v  '•'water  =  ;  (3.37)  c  /_£iitt_jv) ^ A/  "water  ice  '  where the boundary and initial conditions are  T(x,*) = 2>  (3.38a) (3.38b)  T(o,t) = T -> oo,t) = T  c  0  r(aj,o) = r . 0  41  (3.38c) (3.38d)  The latent energy-loaded heat capacity method includes the energetics of phase change without tracking the exact location of the frost front itself, this being of secondary interest. However, given the same numerical boundary conditions the temperature solution obtained using a successful apparent heat capacity should be approximately equal to that obtained from the Neumann solution after interpolating the frost front location out of the numerical solution temperature profile.  Time, years  F i g u r e 3.6. Compare the numerical solution using the apparent heat capacity with the Neumann solution for the problem of freeze-down in a semi-infinite domain. The heavy grey line is the exact Neumann solution. The thin lines are numerical solutions using the apparent heat capacity method. The plot on the right shows the discrete cell spacing versus depth for each of the numerical solutions.  Figure 3.6 shows a comparison of numerical solutions using the apparent heat capacity method and the Neumann solution. The thick grey fine is the Neumann solution, while the thin lines are numerical solutions of various coarseness. The yellow, black, and purple lines correspond to discrete staggered grids with fine resolution nearest the top of the domain, and these solutions closely track the exact Neumann solution. The blue and green curves of the progression of the frost front are 'choppier' due to the coarser spatial resolution of these discretizations. 3.6.3 T e s t i n g the M i x i n g E q u a t i o n s The mixing equations must satisfy two requirements. First the discrete layered column must maintain the correct structure as material is added and removed from the system.  42  Second the constituent mixing must generate realistic values for the bulk properties of the soil mixture. Figure 3.7 illustrates the results of applying Equations (3.27) and (3.28) to a single phase domain with three distinct soil layers and a seasonal snow layer (compare this to Figure 2.5 which illustrated conceptually the effect of the grid transformation on the layered domain). Additionally Figure 3.7 shows a comparison of a linear and exponential distribution of cells in the transformed domain. Each case has the same solution size and the same layered structure in real (z) space, however in the exponential case the discrete cells are concentrated in the upper part of the domain.  Figure 3.7 domain  Effect of the spatial variable transformation on the discrete layered  43  While there are a number of catalogs of the bulk thermal properties of common soil types, there are fewer studies which provide both the bulk properties and the corresponding constituent properties. Penner (1970) reports thermal conductivity measurement on two frozen saturated soils, comparing measured and calculated values. Table 3.1 summarizes the comparison of the measured properties of the two soils with the thermal conductivities calculated using Equation (3.26).  Soil  type  Leda Clay Sudbury Silty Clay  K  (S>  measured K (fully frozen)  Porosity  calculated K (fully frozen)  (B>  measured (partly  frozen)  calculated K ^ (fully frozen)  1.72  0.7  1.84  1.98  1.34  1.44  2.92  0.45  2.76  2.52  1.97  2.06 (nil Rvalues in W/mK)  Table 3.1 Comparison of measured and calculated values of thermal conductivity  3.6.4 Testing the Mass/Energy Conservation of the Code The  mass of the system is  M(*)= / Jv where x  m  p(x ,t)d r  (3.39)  3  m  is location in space. In the discrete domain this integral is approximated as  n  M (*)  =  X  A  ( - ) 3  40  i=l where pi(t) is a known discrete density field that is computed at each solution time. The change in mass with time can be approximated as  dAfi  _  dt  ~  -  M (t x  A l  -  At)  (3.41)  However as snow accumulates or ablates at the upper boundary the total mass of the system changes at a known rate 44  dM  dH  2  - ! T  s  =  p  s  - d T '  (  3  -  4  2  )  where ps is the density of the snow layer with thickness Hs(t). Thus it should also be possible to track the total mass at each time by book-keeping the contribution of the changing snow layer thickness over each discrete time step:  M (t) = M (t-At) 2  + p j (H (t-^))At.  2  S t  (3.43)  s  Equations (3.40) and (3.43), and Equations (3.41) and (3.42) should give the same results if the numerical scheme conserves mass. The total internal energy of the system can similarly be expressed as the integral  U(t) = / CTd r, Jv  (3.44)  3  which is approximated from the discrete solutions as n  U (t) = Y C T (t)AV (t). 1  J  i  i  (3.45)  i  i=i  The change in internal energy can be approximated as dUi dt  U(t) - U(t - At) At  (3.46)  If there is no change of phase in the domain then the internal energy of the system changes only due to conductive heat fluxes and convective energy fluxes through the domain boundaries  ^  =  C s T s ^  -  K,mj  +  (3.47)  K.Q,  dt dt V oz / s where Cs is the heat capacity of the snow layer of thickness Hs(t) at temperature Ts(t); Ks is the thermal conductivity of the snow layer; K B is the thermal conductivity of the material at the bottom of the domain and G is the geothermal temperature gradient. Again it should be possible to track the internal energy of the system by adding the contributions of the boundary fluxes over each discrete time step:  U {t) = U {t - At) + AtCsTsit)^ 2  2  45  _ K {^) S  At + K GAt. B  (3.48)  - 1 9  a)  I  i i i i  i i i !  1  i i i  -21  Time, t, (years)  Time, t, (years)  x  10  j  I  | 1 •s iw .5  i i  I I  it-  1.006  -4  ! I  I i  1.......  ... _  ...  I  I !  x I  e)  !  10  f)  i  l  I  j.  5*  S §I 1.01 fi & 5 1  I  I  "1  I  K L I  1  | - ~] < K  j  3.  J ...._ 1  /  tafinpsnsL  _ _  ...  i i  I I I  It  .. _. L .  L -  I  <znmu - -  Time, t, years  Figure 3.8.  I I I  '  1 Time, t, years  Time, t, years  1.02  <J  (  I I  I  1.002  - - -1 —  I, —  I I I  1 1 I  i „ _ I  1 1 1  Time, t, years  G l o b a l mass a n d global energy  Figures 3.8a and 3.8b show the surface boundary conditions for testing the mass and energy conservation of the numerical model. The global mass Mx{i) as calculated from Equation (3.40) is shown as a solid line i n Figure 3.8c; M ( i ) as calculated from Equation (3.43) is shown as red dots. The mass change rates dM\/dt and dM^jdt calculated from Equations (3.41) and (3.42) are shown respectively as a solid line and red dots in Figure 3.8d. Figure 3.8e and 3.8f are a similar presentation of the global energy and change in global energy; and dU\/dt appear as solid lines, U~2{t) and dU2/dt appear as red dots. 2  46  4 APPLICATIONS In chapters two and three a model was developed to solve the one dimensional vertical temperature profile for a layered soil column. The boundary conditions on the model are a fixed lower domain boundary where the temperature gradient G is known, and a mobile upper boundary where both the location H{t) and temperature Ts(t) are known but not constant in time. In addition to the foregoing boundary conditions the material properties and structure of the layered column are specified by a suite of six properties for each layer: (S) (S) (S) H , K , Ca , Pa , 4>ai faThe a subscript indicates a particular layer. The soil solid properties of thermal conductivity, heat capacity, and density are collectively dependent on the lithology of a particular layer. In this chapter I apply the model to a number of synthetic applications following the seminal Goodrich 1982 paper, I choose to parallel this paper because it rather simply and clearly illustrates the effect of different model variables on the solution to a I D temperature profile. The comparison is more one of convenience than validation. The next development is to modify the model implementation to address the question of permafrost occurrence and extent rather than the question of vertical temperature variation, again using synthetic soil profiles inspired by the Goodrich work. The final part of the chapter applies the model to some real world data from northwest Canada. a  a  4.1 Synthetic Applications 4.1.1 Sensitivity to Soil Column Parameters - Goodrich,  1982  Goodrich (1982) presents a numerical study of the effects of snow cover on ground temperatures which provides a convenient suite of comparison tests. The idealized soil column and surface boundary condition are not intended to represent a real world application and the results are essentially qualitative. The study is intended simply to illustrate the effect of varying parameters of the model. The soil column is a homogeneous profile of either fine or coarse grain soil, possibly overlain by a thin organic surface layer and the effect of a seasonal snow pack is compared to a bare ground case. Setup and Inputs Two surface temperature regimes are considered. The seasonal frost case has a M A S T of 5°C and a SSTA of 15°C. The permafrost case has a M A S T of - 1 0 ° C and SSTA of 20° C. The initial condition for both cases is a linear profile from M A S T at the surface increasing with depth following a geothermal gradient of —0.056 K m " , Each model is run to steady-state and 15 years beyond with a time step of 15 days. The vertical domain is 50m, discretized into 100 points following the concepts of Figure 3.1 and Equation (3.2); the shape parameter B has a value of 5. 1  Table 4.1 summarizes the thermal properties of the four material types considered in the study. The numbers within the grey box are those provided by Goodrich. The coarse soil corresponds to a sandy soil, the fine soil to a silty clay, and the organic soil to a peat (Andersland and Anderson, 1978).  47  Type  A-;'"  Snow  0 18  ,„  w  , HI)  A*  f  L  —  B R  0 52  Organic  1.20  "•' '-0 40  Pine Grain  221  1 13  1 92  Coarse Grain  •soi  " ,\2.19  • .-'".1.86  '  — i i l l l i P iippi P 170 ~ 18*11111 2i"'9 I l l l p l l 0.35 " j 152.0  &~  3 89  >2 40  1  Bedrock  <s>  e  250\  MM  >2 30  ,:  Sill  c  (S>  K  P  <S>  0  —  0.85  0.038  2.5  1000*  0.45  2.1  1.85  1550  0.26  3.4  1.81  2014  0.05  2.9  2.0  2650  (see  text  v  " ' " ~ '  for explanation  and  units)  Table 4.1 Thermal properties of soil material types in Goodrich (1982). A l l conductivity (K) values are in W m K , all heat capacity (C) values are in M W m ~ K , densities are in kg m ~ , soil latent heat L has units of M J m ~ . The subscript / indicates the frozen state while the subscript t indicates the thawed state. _ 1  3  _ 1  - 1  3  3  y  The material properties used by Goodrich are provided as bulk values for the thawed and frozen states however Equations (3.27) and (3.28) require only the properties of the solid constituent. The product of saturation and porosity is known as the volume wetness (0) which can be determined from the numbers provided by Goodrich in two ways:  p  (B)  %W)  (- )  W  4  la  U' ' Here W is the mass based water content, p( ) is the density of water, Ly is the volumetric latent heat of liquid water and the other properties are as described in Table 4.1. W i t h the tentative assumption that all the soil types are fully saturated ( / = 1), the volume wetness fraction 8 is equivalent to the porosity <f> and Equations (2.40), (2.42) and (2.48) can be re-arranged to solve for C > and if< >: ( 4  l b )  w  <s  s  '  c ( s }  I  S  )  = W  ( 4  CJ^H  =  48  -  2 )  { i 3 )  (  (4.4)  Equations (4.2), (4.3) and (4.4) can be solved using the frozen and thawed bulk values provided by Goodrich. The average values of the properties of the soil constituents are given in Table 4.1. However solving Equation (4.2) for the organic soil gives a curious result. The high value of soil latent heat for the organic soil suggests a high porosity and water content. However the bulk density of the organic soil provided by Goodrich is 1 7 0 k g m , so low that Equation (4.2) gives a negative value for the density of the soil solids. As listed in Table 2.2 a typical density for organic soil, for example an organic mud, is 1300 kg m , not much greater than that of water. The density given by Goodrich would only be possible if the organic soil was completely dry, in which case the given values of mass based water content and soil latent heat are not possible. Due to this inconsistency a density of 1000 kg m is chosen for the solid constituents of the organic soil. The properties of the bedrock material are also taken from Table 2.2. - 3  - 3  - 3  Presentation of Results The results of each run are presented as a triplet of curves. These curves are the envelope and mean of temperature variations versus depth. Under conditions of seasonally varying surface temperature, the instantaneous temperature profile in the ground will show the penetration of the cyclical surface anomaly, diminishing asymptotically in a whiplash shape to the geothermal gradient, Figure 4.1. Over a seasonal period the extremal excursions of the instantaneous temperature profiles trace out trumpet shaped maximum and minimum envelope curves. The M A G T profile is taken as the mean of the envelope curves. Characteristics of the ground temperature, for example the top and bottom limits of the permafrost and the depth of zero annual amplitude, are most easily identified from the envelope and mean curves. For each model run the mean and envelope curves are generated from 15 model years of integration after the solution reaches a steady state.  -10  Temperature, °C  -5  0  f|  5  10  Permafrost  Bottom of Permafrost  Figure 4.1. Relationship of instantaneous temperature profile to the envelope and mean profiles (Andersland and Anderson, 1978)  49  C o m p a r i s o n of Soil T y p e s  For both the seasonal frost and permafrost climates, the coarse grain soil exhibits a wider envelope of seasonal variation in temperature as both winter frost and summer warming penetrate deeper into the coarse soil. Figure 4.2 presents the results of the model for both the permafrost and seasonal frost cases on a common graph. The dashed curves correspond to the coarse soil, and the solid curves correspond to the fine soil. Maximum and minimum envelope curves are shown in red and blue respectively; the black curve is the mean annual ground temperature profile ( M A G T ) . The bulk conductivity of the coarse soil is almost twice that of the fine soil for the thawed state, with the result that summer heat penetrates the active layer of the coarse soil more easily than the fine soil. As a result the M A G T of the coarse soil for the permafrost case is about 1 K warmer than for the fine soil. For the seasonal frost case the difference in the M A G T profiles of the fine and coarse soils is only in the active layer where the higher porosity fine grain soil shows a larger thermal offset due to the larger ratio of frozen to thawed conductivity.  F i g u r e 4.2. Comparison of effect of soil type on temperature profile for seasonal frost and permafrost conditions. The red and blue curves are the maximum and minimum envelope curves of the seasonal temperature variation. The black curves are the M A G T profiles  50  Effect of the snow season parameters The effect of a winter snow cover on the underlying ground temperature is investigated using the snow accumulation model described by Equations (2.25) and (2.26). In all cases the parameters ni and n have respective values of 1/2 and 3/2. The length of the winter season is the number of days between the first date of permanent seasonal snow Pi, and the last day of permanent seasonal snow P 4 . Figure 2.3 illustrates these parameters. The soil column is a homogenous profile of fine grain soil and the surface boundary condition is the permafrost climate described above. 2  Thickness of Snow Figure 4.3 illustrates three model runs. The dotted line is a reference case showing the M A G T profile and envelope of seasonal variation for a fine grain soil column in a steady permafrost climate with no winter snow cover. The dashed and solid lines show the M A G T and seasonal variation profiles for the same soil column and climate with winter snow cover reaching a maximum of l m and 0.5 m respectively. In this comparison the winter season is 155 days in each case.  0 E  i..  'I  !' X '> V \. '• 1 1  1  1m snow  !  1 1 1 1 1  .  0.5m snow no snow  -30  _ _  I  i .. .. I. . - 2 5 - 2 0 J  -  '1  '1V ' f - \ — u1 1 - 1*1] 1 1.11 f~ 1 M; I ilt ! 1 1 i I 5 -  '1 10  -  5  r  1 1" ~ 1 1 1 0  j !  !  "--[—•]-—•• 5  10  15  20  Temperature, "C  Figure 4.3. Effect of thickness of winter snowfall in the permafrost climate. The dotted line is a reference profile with no winter snowfall. The dashed and solid curves are profiles with maximum winter snow cover of 0.5 m and 1 m respectively  The effect of a snowpack is to insulate the ground and limit the release of heat during the winter months. The M A G T profile is warmer by several degrees for the snow covered cases as compared to the bare ground case. The presence or absence of the winter snow layer has a more profound impact than the thickness of snow, as the doubling of the maximum snow depth results in a warming of the M A G T of less than I K .  51  Duration of Snow Cover Figure 4.4 shows the effect of the length of the snow cover season on the mean ground temperatures and in the winter heat loss. As in Figure 4.3 the dotted line is a reference profile for a fine grained soil profile in a permafrost climate with no winter snow cover. The dashed and solid curves correspond to winter seasons of 155 days and 170 days respectively. The snow cover reaches a maximum depth of 0.5 m on the same date in mid February in both cases, however for the short winter (155 days) the onset of permanent seasonal snow cover occurs later, and the onset of spring melt occurs earlier than in the long winter (170 days) case.  8  i^yrXA  - 5 Snowcover  <W-10 §:  •  o  -15  U CO  Q.-20  a  -25 -30  Figure 4.4.  -  Duration  |  '\  15S days  J  V(  170 day "\ no snow ,  ,  !  -25  -  2  j ^  \\  !  0  -  {  hj-jj  j _,  !  \  1  1  5  -  1  0  -  i  i  !  ]  1 1  11  1 1  1 1  5  0  i •  I 1 • 1 1 1  -i---i  I  1 1 1  \  I  1 1  1  10  Temperature, °C  1 1  15  20  Effect of length of winter season in the permafrost climate  The profiles for the longer winter season show the amplitude of seasonal temperature variations attenuates more rapidly with depth. The winter minimum temperature (solid blue curve) is warmer below the active layer due to the snow cover acting as an insulator and preventing escape of accumulated summer heat. Both Figures 4.3 and 4.4 illustrate that in a permafrost climate, winter snow cover is associated with warmer ground temperatures.  52  Effect of T h i n Organic Layer A column of fine grain soil is overlain with a 20 cm layer of organic material. This organic layer has a high water content and low thermal conductivity. The surface boundary conditions is a seasonal frost climate and both the bare ground, and winter snow cover cases are shown. The winter season is 155 days long with a maximum snow depth of 0.5 m. Organic Layer with No Winter Snow Cover Figure 4.5 presents a comparison between the temperature profiles of a homogenous fine grain soil column (solid lines) and a similar soil column with a shallow surface layer of organic material (dashed lines).  Figure 4.5. Effect of a thin organic surface layer, with no winter snowfall, in the seasonal frost climate. The solid lines are a reference case of fine grained soil with no surface layer  The temperature profile of the homogeneous column has a smaller thermal offset and a narrower amplitude of seasonal variation in the temperature profile compared to the column featuring the surface layer. With the organic layer the mean annual ground temperature profile is 1 to 2 degrees cooler, mainly due to extra winter cooling of the non-frozen portion of the profile. This is because the high porosity of the organic layer holds a very high water content and the contrast in thermal conductivity of frozen and unfrozen water controls the seasonal profile. During winter when the active layer (including the organic layer) is frozen, heat is more easily lost from the ground, cooling the ground overall more than in summer when the low conductivity of the liquid pore water inhibits penetration of summer heat.  53  Organic Layer with Winter Snow Cover  -5  0  5 10 Temperature, ° C  15  20  Figure 4.6. Effect of a thin organic surface layer, with winter snowfall, in the seasonal frost climate. The solid lines are a reference case of fine grained soil with no surface layer  The differences between the profiles for the homogeneous soil column (solid) and the soil column with surface layer (dashed) are very small compared to those displayed in Figure 4 . 5 . This comparison shows that the presence of a winter snowpack may have a greater effect on the ground temperatures than the presence of an organic surface layer. 4.1.2 Solving for Multiple Possible Soil Columns The case by case comparison of the previous section is a tedious approach to analyzing the ground temperature solutions, especially if the question of interest involves not so much the details of temperature at every depth in the ground, but rather the consequence of the temperature profiles at a number of location on the extent and occurence of frozen ground. The next step is to modify the implementation of the ground temperature model for application to questions of broader spatial scope. This modification is referred to as ensemble modelling as it is concerned not with a single soil column and temperature solutions, but rather a collection of several possible soil columns and temperature solutions. Conceptual development It is possible to imagine some area of interest over which there is a range of soil types from fine to coarse and discontinuous surface layers of organic material and drifted snow. In such an area each of the soil profiles and temperature solutions investigated in the Goodrich study might be found at a discrete location, however if a borehole were dug it is unlikely that the measured temperature profile would exactly match any of the solutions. Figure 4 . 7 illustrates a hypothetical inhomogeneous region where bedrock (green) is overlain by layers of sediment (greys), organic material (brown) and seasonal drifted snow (blue) of varying thicknesses and extent. Each column has a vertical extent of Ho from the solid ground surface to the lower boundary. 54  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19 20  21  22  23  24  25  26  27  28  29  30  31  32  Figure 4.7. Hypothetical ground section and 32 possible soil profiles. The green, grey, white, and orange layers are soils of different lithologies. The discontinuous blue surface patches represent drifted snow  Note that the vertical scale of this example is greatly exaggerated compared to the horizontal scale. The lower cartoon of Figure 4.7 is a depiction of 32 soil columns that each could be found somewhere within the the ground section. Each of these columns can be constructed as a four layer soil profile according to the details of Table 4.2. Model application and solutions Table 4.2 identifies the maximum and minimum value for each of: maximum snow thickness, organic layer thickness (a), upper soil layer thickness (b), upper soil layer lithology (b), and lower soil layer lithology (c). The model can be solved for all possible combinations of these variables to yield a suite of solution temperature profiles. W i t h 5 variables, each allowed two possible conditions, the size of the solution set is 2 = 32, as depicted in Figure 4.7. Note that the thickness of the lower soil layer (c) also has two possible values, however these are dependent on the thicknesses of the overlying layers. 5  Figures 4.8 and 4.9 are the solutions for the set of soil profiles described in Table 4.2 under the permafrost and seasonal frost climates previously described (section 4.1.1). The upper chart of each figure shows the seasonal temperature variation envelope curves for  55  Layer  Property  Snow  Maximum Thickness  Snow  Thermal Properties  a a b b c c  Variable  Value (see Table 4.2) [0, 0.5], (m)  K  snow  [snow]  ^ snow  Thickness  [0,0.2], (m) jJS)  Lithology, Porosity, Saturation  „(S)  t>o  fa  Thickness  [organic] [10,30], (m)  Lithology, Porosity, Saturation Kb  Lb  Pb  <t>6  lb  Thickness  [fine grain, coarse grain] IHo-(H  a  Lithology, Porosity, Saturation K! c  C  C  pc  <t>c  fc  + H )], (m) b  [coarse grain, bedrock]  Table 4.2. Properties of variable layers in Figure 4.6  each of the 32 possible soil profiles. The lower chart of each figure is a histogram of the the depth of seasonal freeze/thaw from the surface of the column (active layer depth). In each climate case the histogram of active layer depth shows a bimodal distribution which can be attributed to the discontinuous winter snow cover. For both the permafrost and seasonal frost cases the majority of solutions place the active layer depth between 1.5 and 2.5 meters. The histogram of active layer depth for the seasonal frost climate has almost half the total population of solutions with no winter freezing of the ground; this is due to the insulative effect of the winter snow cover. The input requirements listed in Table 4.2 are only slightly more extensive than what would be needed to apply the model to a single unique soil column. W i t h the exception of a few research locations that have been studied in detail, comprehensive information about the subsurface structure and thermal properties of the ground is not available. The ensemble approach to modelling is convenient in that it does not require fine detail of input and is flexible to use what input data is available.  56  Figure 4.8. Solution temperature profiles and histogram of active layer depth for ground section of Figure 4.7 in permafrost climate with drifted snow, 32 solutions  To demonstrate the hypothetical set of layer however the surface soil saturation. The addition  flexibility of the ensemble approach, Table 6 lists a different properties. In this case winter snow cover is not a possibility, layer can be organic or fine grain with two possible values of of another variable increases the size of the solution set to 64.  Figures 4.10 and 4.11 display the 64 temperature solutions in maximum and minimum seasonal temperature profiles and histograms of active layer depth. Due to the greater thickness of the surface organic layer the majority of solutions put the active layer depth between 1 and 1.5 meters for both the permafrost and seasonal frost cases.  57  -5  10  0  Temperature, C 0  1.5 'live Layer Thickness, m  Solution temperature profiles and histogram of winter frost penetration depth for ground section of Figure 4.7 in seasonal frost climate with winter snow, 32 solutions F i g u r e 4.9.  Layer  Property  Variable  Snow  n/a  nla  a  Thickness  a  Lithology, Porosity  xf cf „ <(>„  [organic, fine grain]  a  Saturation  fa  [0.5, 1]  b  Thickness  H  b  Lithology, Porosity, Saturation rfb'c'b'pb <t>6 fb  [fine grain, coarse grain]  c  Thickness  [Ho - (H + Hi,)], (m)  c  Lithology, Porosity, Saturation K f Cf'pe  Value (see Table 4.1) nla [0.5,3], (m)  P  [10,30], (m)  b  H  a  c  T a b l e 4 . 3 . Hypothetical layer properties  58  4>c fc  [coarse grain, bedrock]  Active Layer Thickness, m  F i g u r e 4.10. Solution temperature profiles and histogram of winter frost penetration depth for Table 4.3 in seasonal frost climate, 64 solutions  Active Layer Thickness, m  Figure 4.11. Solution temperature profiles and histogram of active layer depth for Table 4.3 in permafrost climate, 64 solutions  59  4.2 Real World Applications Along the northern coast of northwest Canada, permafrost is continuous and ubiquitous in occurrence with local thicknesses in excess of 600 m and near surface ground temperatures below —8°C. In southern Yukon permafrost is sporadic and discontinuous with maximum local thicknesses of 20 m and mean annual ground temperatures above — 1°C. The wide variability of ground conditions over a limited distance (10 degrees of latitude) coupled with the rapid response of Arctic climate to global warming make this part of the world especially interesting (Burn, 1998).  Figure 4.12. Location of select boreholes and weather stations in Yukon Territory and Mackenzie Valley region. Boreholes are shown in parentheses and labelled with their catalog number (Peterson and Vose, 1997; Taylor et al., 1982)  Figure 4.12 shows the locations of boreholes and weather stations in Yukon Territory and Mackenzie Valley region. In general deep boreholes were drilled for oil exploration purposes. The weather stations have reported monthly temperature and precipitation data for several decades until the early 1990s (Peterson and Vose, 1997). The oldest boreholes were drilled in the early 1940s in the Norman Wells area, but most of the listed holes date from the 1970s and mid-1980s (Taylor et al., 1982). The process of drilling boreholes warms the surrounding ground by several degrees, and this warming anomaly may not have equilibrated in all holes before the temperature profiles 60  were measured. Legal abandonment requirements for petroleum exploration sites have made many of these boreholes inaccessible now. Figure 4.13 is a map of permafrost occurrence over this area as determined from published maps (Cogley, 1995). The blue boxes outline three regions of different climate and permafrost occurrence. These regions are: the delta of Mackenzie River between Inuvik and Tuktoyaktuk where permafrost is ubiquitous, a section of the Mackenzie Valley centred around Norman Wells where permafrost is common, and the band of Yukon Territory from the town of Mayo south to the border with British Columbia where permafrost is sporadic and discontinuous.  F i g u r e 4.13. Permafrost occurrence in Yukon Territory and Mackenzie Valley region (Cogley, 1995)  Using the principles developed in section 4.1.2 the ground temperature model is applied to each of these focus areas. Data on the soil column structure are taken from published maps of soil depth and texture (Zobler, 1986; Webb et al., 1993). These maps are spatially coarse with a resolution of one degree of latitude by one degree of longitude. At this spatial scale it is difficult to compose a single soil column representative of the entire area. Qualitative and anecdotal descriptions of ground conditions are also available from a wide range of sources (e.g. French and Heginbottom, 1983; Dyke and Brooks, 2000; Smith, 1975). The information from these sources can be combined to construct a suite of multiple possible soil columns and surface climates based on the range of properties within 61  each region. Each simulation is run with a periodic steady state surface boundary condition. As before, the model integration proceeds for 15 years beyond the time necessary to reach a periodic steady state solution. Where relevant a seasonal snow cover is applied using the accumulation model of Equation (2.25) and (2.26). 4.2.1 Mackenzie Delta The delta of the Mackenzie River has some of the deepest permafrost in North America. The microclimate and surface conditions vary strongly over very short distances. Figure 4.14 is an overview map of the region showing the course of the Mackenzie River from the southeast corner to the delta complex on the coast of the Beaufort Sea. The region selected for ensemble modelling is between 67 and 70 degrees of north latitude, and 131 and 138 degrees of west longitude.  1 , 1  Beaufort  pcfe  ao  1^  I  JU--H W ,  r 7 U S  I  II7S)  0  _  Tuktoylkyuk  1  I I  ,(377) (179} I P. o  ^^t"* • Shingle Point Airport  W  ^  (  1076)  I I tms)\ a  ' J  %\ IY Akltvlk I 1  tnuvlk "fa^l j  132  136"W  W  Longitude  Figure 4.14. Lower Mackenzie River and Delta. Local weather station are labeled with a star, deep borehole locations are indicated by a circle and the borehole catalog number in parentheses (Peterson and Vose, 1997; Taylor et al. 1982)  The numerical labels in parentheses in Figure 4.14 are the locations and catalog numbers of deep boreholes (Taylor et al., 1982). The locations of five weather stations are also indicated. In the Mackenzie Delta the permafrost is continuous and typically extends to a depth of several hundred meters. 62  Climate description Figures 4.15 and 4.16 show the monthly air temperature and precipitation records for Inuvik and Tuktoyaktuk. The records for Shingle Point, Aklavik and Fort McPherson are similar. The precipitation record in the upper graph of each figure is shown as a bar graph of monthly total precipitation in millimeters water equivalent (mm w.e.) for each year from 1957 to 1991. The lower graph for each figure shows the monthly mean air temperatures as a line for each year of the period of record. The summer (snow free) season is assumed to begin and end when the mean air temperature rises above 0 °C as indicated by the blue lines joining the two graphs.  July February  Figure 4.15. 1997)  April  September November August October December  Inuvik weather station records (1957-1991) (Peterson and Vose,  Historical weather station records from Tuktoyaktuk, Shingle Point, Inuvik, Aklavik and Fort McPherson give an average M A AT for the area of —8.5°C with a seasonal amplitude of between 18 and 22 K . The snow free season extends from May to late September. Smith et al. (1998) report a M A A T of - 8 . 5 ° C for the years of 1993-1994 recorded at a site midway between Inuvik and Tuktoyaktuk, and estimate the soil surface temperature as between —5.5 to —1.7°C.  63  Tuktoyaktuk  Figure 4.16. Vose, 1997)  January March February April  May  January March February April  May  July  September November August October December  July  September November August October December  June  June  Tuktoyaktuk weather station records (1957-1989) (Peterson and  Soil description Smith (1974) describes five study sites along a short transect of a distributary of the Mackenzie River between Inuvik and Tuktoyaktuk. Soils are thin with mainly sand or silt composition with sporadic peat and moss surface cover. In the lacustrine and intertidal deposits around Tuktoyaktuk significant areas are covered by peat layers between 1 and 2 m thick (French and Heginbottom, 1983). Two years of snow cover data recorded at Inuvik show maximum snow depth of 40 to 70 cm occurring in February and May through September free of snow. However snow accumulation varies with prevailing wind and vegetation distribution. In areas with little or no vegetation snow is scoured and bare ground is exposed. Deepest drifting occurs where willows trap blowing snow; under forest canopy the snow cover is thinner and more variable due to trapping in the crowns of trees. Figures 4.17 and 4.18 are maps of soil depth and soil texture with a resolution of one degree latitude by one degree longitude (Webb et al., 1993; Zobler, 1986). French and Heginbottom (1983) report the common soil textures of the Mackenzie Delta to be silty loam and loam with increasing amounts of clay moving up the Mackenzie River Valley.  64  Longitude  Figure 4.17.  Soil depth (Zobler, 1986)  Figure 4.18.  Soil texture (Webb et al., 1993)  65  Regional model  From the available data 16 potential soil profiles are assembled as shown in Figure 4.19. There are 8 configurations of two or three soil layers labeled a, b and c. Layer a is an organic horizon and layer b is a sediment horizon with a minimum thickness of 10 cm and a maximum thickness of 3.5 m. The sediment layer texture varies from fine silt to coarse sand. The material properties for the soil layers are those listed in Table 4.1. The bottom layer labelled c in each profile is assigned bedrock values. Each soil layer configuration is solved for a bare winter case, and a snow-covered case, with snow depth reaching a maximum of 20 cm, and a winter season lasting from mid September to mid May. A mean surface temperature of —8.3°C and seasonal temperature amplitude of 20°C are chosen based on the weather station records.  b 1  Is  MAAT =  •8.3 "C  SSTA  20  -  "C  Maximum Snow Depth =  [0 0.2) meters  Organic Layer Thickness, a •  10 0.5] meters  Soil Texture, b =  [Fine silt Coarse  Soil Thickness, b =  )0.1 3.5) meters  Total Column Thickness =  400 meters  Sand)  Geothermal Gradient =  IL*  lr  Figure 4.19. Model soil columns and boundary conditions for the Mackenzie Delta calculations  The upper graph of Figure 4.20 shows the model solution maximum and minimum temperatures versus depth for the 16 potential soil profiles. The effect of snow cover is immediately obvious in the two major groupings of the maximum curves. When snow cover is absent the ground temperature is warmer and the depth of zero annual amplitude is greater. The lower graph of Figure 4.20 is a histogram of the depth of summer thaw for the suite of 16 solutions. Active layer thicknesses greater than 3 meters occur only in the absence of winter snow cover. Where snow cover is present the active layer most commonly will reach a depth of around l m . This result is supported by Smith (1975) whose field study at Reindeer Station on the Mackenzie River north of Inuvik measured active layer depths of 110.6 cm, 103.3 cm and 99.2 cm.  66  0  2  4  6 8 10 Active Layer Thickness, m  12  14  16  Seasonal maximum and minimum ground temperature solution profiles, and histogram of active layer depth for the Mackenzie Delta regional model Figure 4.20.  In Figure 4.21 the 16 mean annual ground temperature profiles as calculated from the model (blue) are plotted with the ground temperature profiles (black) measured in deep boreholes in this region. Each of the borehole profiles is a 'snapshot' measurement, and the date of record is not available in all case. For this reason comparisons between the model solution temperatures and the borehole temperatures are most relevant below the depth where seasonal variations penetrate; the borehole temperatures at depths shallower than 10 m are not reported. Additionally, following drilling of the borehole it takes several years for the heat of drilling to dissipate and the borehole temperatures to re-equilibrate to the ambient profile. In most cases the borehole temperatures have been recorded before this re-equilibration is complete, so the profiles may be offset warmer by a few degrees. This is likely the case with borehole (177) that records virtually no permafrost despite its proximity to borehole (176) and borehole (179), which are the coldest profiles in the region.  67  Temperature, (  °C)  Figure 4.21. Comparison of borehole and solution temperature profiles for the Mackenzie Delta region. The dashed curves are temperature profiles recorded in boreholes within the focus region. Each curve is labelled with a catalog number in red, and Figure 4.14 shows the location of these boreholes. The blue curves are the solution MAGT curves to the ensemble model illustrated in Figure 4.19  All the borehole profiles show a trend toward a MAST warmer than —8°C, however this is not surprising as that choice for the surface boundary condition for the model was made based on air temperature records and as was explained in Chapter 2 the MAST is usually a few degrees warmer than the MA AT. Nevertheless the solution profiles fit well within the family of borehole profiles for the Mackenzie Delta region. There are similarities in the near surface inflections of boreholes (179), (176) and (089) which in the solution profiles are associated with the presence of a winter snowcover. The solution profiles place the depth of the base of the permafrost at 330 m, in agreement with borehole (179). The broad curvature in the upper 150 m of the solution profiles for the bare winter case is also subtly echoed in the profiles for boreholes (179), (176) and (089). In the real world 68  the winter snowcover varies from year to year with seasons of no snow and seasons of deep snow. The deep borehole temperatures represent an integration of these changing surface conditions over time, so it is interesting that the measured profiles should share characteristics of both model conditions. 4.2.2 Mackenzie Valley  The valley of the Mackenzie River connecting Great Slave Lake to the Beaufort Sea traverses the discontinuous and continuous permafrost zones. The Mackenzie and Richardson mountains rise to the southwest of the valley. The uplands at Norman Wells host hydrocarbon resources that were the motivation for drilling deep boreholes in this region in the 1970s. Figure 4.22 is an overview map of he region showing the course of the Mackenzie River through the town of Norman Wells. The elevated topography of the Mackenzie Mountains is in the lower left of the figure. The region selected for modeling extends from 63.5 to 66 degrees north latitude and 124 to 130 degrees west longitude.  124 "w  Figure 4.22. Mackenzie River valley near Norman Wells. Deep borehole and weather station locations are shown. Boreholes are shown in parentheses and labelled with their catalog number. (Peterson and Vose, 1997; Taylor et al., 1982)  69  Numerical labels in parentheses in Figure 4.22 are the location of boreholes nests along the Mackenzie River valley near Norman Wells. Each borehole nest consists of one or more borehole records; these have usually been collected from separate holes within a few kilometers of each other. The Norman Wells station provides the only record of weather data in this region. In this region permafrost is discontinuous but common, with thicknesses of 45 m or more (Brown, 1970). Climate description Only the Norman Wells weather station in this region has a long record of monthly temperature and precipitation summarized in Figure 4.23 (Peterson and Vose, 1997). However field studies have shown that mean monthly temperatures along the Mackenzie Valley bottom do not vary significantly with latitude. Precipitation is light due to the rain shadow effect of the Mackenzie Mountains on westerly low pressure systems. Winter temperature inversions are common with air temperatures at 1000 m elevation up to 10°C warmer than at ground level (Dyke, 2000). Winter precipitation amounts are low enough that snow does not accumulate to significant depths and is easily transported by wind.  Figure 4.23. Vose, 1997)  Norman Wells weather station records (1944-1991)(Peterson and  70  Soil description Figures 4.24 and 4.25 are maps of soil depth and texture. Coarser thinner soils are associated with the higher elevation mountainous regions. The fluvial deposits in the river valley and glaciolucustrine soils of the plateau between the Mackenzie River and Great Bear Lake to the north are finer and deeper.  Figure 4.24.  Soil depth (Zobler, 1986)  Dominant Soil Texture Fine Medium Medium Coarse  128°W  12S°W  Longitude  Figure 4.25.  Soil texture (Webb et al., 1993)  71  tn°w  Regional model  /IMS 7" =  1-6 -1.9] "C  SSTA =  22 °C  [0 0.5] meters Organic Layer Thickness, a = [Coarse Sand]  Soil Texture, 6 =  So/7 Thickness, b =  [1.25 2.0]meters  Total Column Thickness = 200 meters •1 •0.055Km' Geothermal Gradient =  Not to scale  Figure 4.26. Model soil columns and boundary conditions for the Mackenzie Valley/Norman Wells regional model  A model suite of 6 possible soil columns and surface boundary conditions is assembled. Soil material properties are taken from Table 4.1. Two possible surface temperature boundary conditions are used. The first is to use the M A A T and SSTA to directly calculate a ground surface temperature Ts using an annual seasonal sinusoid. The second applies Equation 2.27 to calculate a warmer ground surface temperature from the air temperature. As indicated in Figure 4.26 applying Equation 2.27 results in a mean annual surface temperature of — 1.9°C. The upper graph of Figure 4.27 shows the model solution maximum and minimum temperatures versus depth for the 6 possibilities of soil column and climate. The lower graph of Figure 4.27 is a histogram of the active layer depth for the ensemble of 6 profiles. The deeper active layers are associated with the warmer M A S T . The shallowest active layer occurs when an organic layer is present.  72  1  2  3 4 Active Layer Thickness,  5  6  m  Figure 4.27. Seasonal maximum and minimum ground temperature profiles, and histogram of active layer depth for the Mackenzie Valley/Norman Wells regional model  Figure 4.28 shows the 6 model solution mean temperature profiles (blue) plotted with measured borehole temperature profiles (black). As was described for the Mackenzie Delta region, the borehole records are each an instantaneous measurement of temperature and may be affected by the heat of drilling. This possibly explains why of the 5 well records obtained at Norman Wells, borehole (088), there are two distinct groupings of temperature profiles offset by about 4 ° C . However the solution profiles of the ensemble modelling of ground temperature based on the climate statistics provided from the Norman Wells weather station fall in between the groupings of measure temperatures.  73  Temperature, C ~  0  2  2  20  40  >.  80  100  120  140  Figure 4.28. Measured borehole and model solution temperature profiles for the Mackenzie Valley  4.2.3 Southern Yukon Territory The southern Yukon region is in the Boreal Cordillera ecozone. The region selected for ensemble modeling is 60 to 64.5 degrees north latitude by 126 to 138 degrees west latitude. Figure 4.29 shows the locations of five weather stations and seven deep boreholes distributed in this region. The southwest corner rises into the heavily glaciated St. Elias Mountains while the northeast corner along the border with North West Territories rises into the Richardson and Mackenzie ranges. Between these mountainous zones the interior topography consists of deep narrow alleys incised into a tableland of uniform elevation (Burn, 1993; Smith et al., 1998). Scattered discontinuous permafrost occurs in the region, with up to 65% of ground area in the northern and eastern sections underlain by frozen ground. In the southernmost sections, along the border with British Columbia, sporadic frozen ground underlies less than 35% of the ground area. In the sporadic zone, the mean annual ground temperature is above —5°C and permafrost is less than 20m thick. In central Yukon, around Mayo,  74  Longitude  Figure 4.29. Southern Yukon Territory with weather station and deep borehole locations (Peterson and Vose, 1997; Taylor et al., 1982)  where permafrost is widespread, thicknesses of up to 40 m and active layers from 35 to 100 cm have been measured (Leverington and Duguay, 1997). In the Takhini Valley near Whitehorse the active layer is 1.4 m thick over permafrost extending to 18.5 m depth (Burn, 1998); near the southern town of Teslin permafrost less than 2 m thick has been uncovered in municipal excavations (Burn, 1994). Climate description Figures 4.30 and 4.31 present the historical weather station records for_Mayo and Whitehorse (Peterson and Vose, 1997). These records together with those of Carmacks, Watson Lake and Aishihik indicate a M A A T of around freezing to several degrees below freezing for the region. The rain shadow effect of the St. Elias Mountains and the development of a stable arctic high in the winter result in little winter precipitation in the interior of the territory. What snow does fall is dry and mobile, leading to discontinous drifted deposits whose thicknesses are dependent on microterrain features and vegetation.  75  Mayo  January  March  February  Figure 4.30. 1997)  May April  July June  September August  November  October  December  Mayo weather station records (1925-1990)(Peterson and Vose,  A particular characteristic of the winter climate of central Yukon is the development of a strong topographically enhanced temperature inversion. Dense cold air drainage coupled with minimal solar heating results in temperatures as much as 10° C cooler at valley bottoms than at 1000 m above. The weather stations are generally located in valley bottoms so it is worth noting that ground surface temperatures on hillslopes may be higher than those reported. Additionally orographic effects produce greater precipitation amounts at higher elevations than at the valley bottoms where measurements have been recorded. Permafrost is maintained in valley bottoms by cold winter temperatures and minimal snow accumulation; at higher elevations the more persistent winter snowcover insulates the ground from a nevertheless short and cool summer season, resulting in the same effect (Burn, 1994). Temperatures at the soil surface are expected to be a few degrees higher than the air temperature, however the precise relationship between soil surface temperature and air temperature is complex. The soil surface temperatures at Takhini and Mayo have been estimated at 0.65 °C and 0.75 °C respectively.  76  Whitehorse  January  March  February  Figure 4.31. Vose, 1997)  May April  July June  September August  November  October  December  Whitehorse weather station records (1942-2001 )(Peterson and  Soil description Although this is the largest area to which I have chosen to apply the model, the maps of soil texture and depth show less variability than the Mackenzie regions (Figures 4.27 and 4.28). This seeming homogeneity of the lithologic structure belies the true complexity of the geoclimate in this area. The areas of shallow soil shown in Figure 4.32 are associated with the St. Elias mountains in the southwestern corner, and the Ogilvie and Selwyn mountains in the southcentral region.  77  Longitude  Figure 4.32.  Soil depth (Zobler, 1986)  Longitude  Figure 4.33. Soil texture (Webb et al., 1993)  Regional M o d e l Figure 4.34 illustrates a simple set of eight possible models and boundary conditions which apply to the Southern Yukon. The soil column structures are inspired by the maps of soil depth and texture, Figure 4.32 and Figure 4.34. The parameters for the surface boundary condition (MAST and MATA) are chosen to illustrate the difference between applying a measured air temperature condition versus a soil surface temperature.  78  Warmer  a a a a a a a b b b b a  b c  b  c  c  c  c  MAST = l-5 I] " C MATA=[30  20] "C  Soil Texture, a = b  [Coarse  Sediment Layer, b - Coarse Sand, 10 meters Geothermal Gradient =  c  c  Medium]  Soil Thickness, a = [0.1 135] meters  -0.025Km'  1  c  Figure 4.34. Model soil columns and boundary conditions for the Southwest Yukon calculations  This set is not comprehensive in capturing all the real world variability of this region of interest, but is limited by the availability of quantitative field data. One conscious oversight is the lack of possibility for a seasonal snow layer despite the fact that snow is a critical influence on winter ground temperature in the discontinous region. The importance of snow has been previously illustrated in the synthetic simulations of section 4.1 and the regional ensemble model of the Mackenzie Delta (section 4.2.2). However there is an absence of published ground temperature profiles for areas where consistent winter snowcover is known to control ground temperatures, making it difficult to compare real world and model results for the purpose of model validation. Thus the results of this section are not intended to provide a complete study of the thermal condition of the southern Yukon but are limited to comparison with the available borehole records from Taylor et. al (1982).  79  The upper graph of Figure 4.35 shows the model solution maximum and minimum temperatures versus depth for the 8 possibilities of soil column and surface climate. The lower graph is a histogram of the active layer depth. The temperature curves are in two major groupings, corresponding to the two surface climate conditions. Only the profiles associated with the colder surface boundary condition display perennial sub-freezing ground temperatures with an over lying active layer. Those profiles associated with the MAST of 1 °C show only winter penetration of frost over unfrozen ground and no permafrost.  5  0l  ,  ,  -30  - 2 5  0  1  -20  1  -  2  1  I  5  -  1  1 0 - 5 0 Temperature? C  3 4 Active Layer Thickness, m  5  5  I  ,  ,  10  15  20  6  7  Figure 4.35. Seasonal maximum and minimum ground temperature solution profiles, and histogram of active layer depth for the Southern Yukon regional model.  80  The 8 M A G T profiles calculated from the numerical model are shown in blue in Figure 4.36. The dashed black curves are recorded borehole temperature profiles from the locations labelled in Figure 4.29. None of the available deep borehole records for the focus region show permafrost, nevertheless from field studies, river and roadcuts, and engineering projects, permafrost is known to occur. The profiles calculated from the colder surface climate condition are not matched by any measured borehole record, however the agreement with the warmer climate model is very good.  -6  -4  Temperature, f C) -2 0 2  4  6  Figure 4.36. Measured borehole and model solution temperature profiles for Southern Yukon Territory  81  4.3 Summary of Model Applications The model developed through Chapters 2 and 3 was used to determine the ground temperature fields for a variety of synthetic and real world scenarios. While none of these applications addresses a current question in permafrost research the goal has been to illustrate the types of problems that might be tackled. The ensemble-modeling strategy introduced in Section 4.1.2 is in particular a new tool. A key point is that the input requirements of the ensemble-modelling are adaptable and not strict. In fact the ensemble approach has been developed to suit the type of data that is usually available. The intention is not to solve a single trustworthy temperature profile at a well prescribed location, nor to assume minimal lateral variation in soil structure or surface conditions - but to apply the model over a broader region which does host variations in surface and sub surface conditions and tease out the importance of those variations on the range of possible temperature profile solutions. However application to the large focus region of the southwest Yukon demonstrated the limitations of operating the ensemble model on a very restricted set of input data.  82  5 CONCLUSIONS The goal of this thesis has been to develop a numerical model for solving the ground temperature field of the layered earth in response to a temporally varying surface boundary condition. By applying continuum transport theory the model is able to cope with an inhomogeneous domain whose vertical extent varies in time. Changing the phase of soil water results in an energy contribution that is captured by prescribing the thermophysical properties of water as functions of temperature. The soil material is modeled as a porous constituent mixture of mineral and organic solids and water. Bulk soil properties are calculated from the temperature dependent water properties and a small set of variables which described the layered lithology of the earth. A seasonal snow layer can accumulate and ablate at the surface boundary. The conceptual and qualitative underpinnings of the model are described first in Chapter 2. The local form governing equation is developed in full detail from the global integral form of the energy balance equation. A spatial variable transformation is introduced to allow for a mobile surface boundary. The functional forms for the temperature dependent material properties and the mixing equations for the constituent soil mixture are introduced analytically. In Chapter 3 the model equations are discretized in finite difference and coded into M A T L A B for solution by computer. Comparison with simple analytically tractable thermal diffusion problems confirms the performance of the numerical energy diffusion-advection scheme, the mixing of material properties and the phase change strategyThe flexibility of the model is demonstrated by comparison to a published study (Goodrich, 1982), and by bulk solution of potential temperature profiles of a hypothetical piece of earth with varying lithologic structure. Finally, using field data from northwest Canada, model temperature profiles are compared to measured borehole temperature data. Agreement between model solutions and measured borehole temperatures, for three focus areas is good for the continuous and widespread discontinuous permafrost zones. While the modeling efforts of Chapter 4 are all equilibrium scenarios integrated to steady state under seasonal variation, the results of Sections 3.6.1 and 3.6.2 indicate that the model will perform under conditions of transient forcing as well, for example climate change, or response to forest fires. The priority of this work has been in the conceptual and computational development and verification of modeling tools, however the application of the model to real world data demonstrates the potential for addressing research questions.  83  REFERENCES Andersland, O.B., and D . M . Anderson, 1978. Geotechnical Engineering for Cold Regions, McGraw-Hill Book Company, 566pp. Brown, R . J . E . , 1970. Permafrost in Canada - Its Influence on Northern University of Toronto Press, 234pp.  Development,  Burn, C R . 1998. Field investigations of permafrost and climatic change in northwest North America, Proceedings of the Seventh International Conference, Collection Nordicana 55, 107-119. Carslaw, H.S., and J . C . Jaeger, 1959. Conduction of Heat in Solids, Oxford University Press, 510pp. Cogley, J . G . , 1995. GGHYDRO - Global Hydrographic Data, Release 2.1., Department of Geography, Trent University, Peterborough, Ontario, Canada. Dyke, L . D . , and G . R . Brooks, 2000. The physical environment of the Mackenzie Valley, Northwest Territories: a base line for the assessment of environmental change, Geological Survey of Canada, Bulletin 547. French, H . M . , and J . A . Heginbottom, 1983. .Guidebook to Permafrost and Related Features of the Northern Yukon Territory and Mackenzie Delta, Canada, University of Alaska, 186pp. Goodrich, L . E . , 1978. Efficient numerical technique for one-dimensional thermal problems with phase change, International Journal of Heat and Mass Transfer, 21, 615-621. Goodrich, L . E . , 1982. The influence of snow cover on the ground thermal regime, Canadian Geotechnical Journal, 19, 421-43. Hillel, D . 1980. Fundamentals of Soil Physics, Academic Press, New York, 413pp. Hinzman, L . D . , D . J . Goering, and D . L . Kane, 1998. A distributed thermal model for calculating soil temperature profiles and depth of thaw in permafrost regions, Journal of Geophysical Research - Atmosphere., 103(D22), 28975-28991. Kay, B . D . , and J.B. Goit, 1975. Temperature-dependent specific heats of dry soil materials, Canadian Geotechnical Journal, 12, 209-212.  84  Lachenbruch, A . H . , and B . V . Marshall, 1986. Changing climate: geothermal evidence from permafrost in the Alaskan arctic, Science, 234, 689-696. Leverington, D . W . , and C R . Duguay, 1997. A neural network method to determine the presence or absence of permafrost near Mayo, Yukon Territory, Canada, Permafrost and Periglacial Processes, 8, 205-215. Lunardini, V . J . , 1981. Heat Transfer in Cold Climates, Van Nostrand Reinhold Company, 731pp. Lunardini, V . J . , 1996. Climatic warming and the degradation of warm permafrost, Permafrost and Periglacial Processes, 7, 311-320. Marshall, S.J., 1996. Modelling Laurentide Ice Stream Thermomechanics, PhD thesis, Department of Earth and Ocean Sciences, University of British Columbia. Mellor, M . , 1964. Properties of Snow, Cold Regions Science and Engineering, I I I ( A l ) , 65-72. Nakano, Y . , and J.Brown, 1972. Mathematical modeling and validation of the thermal regimes in tundra soils, Barrow, Alaska, Arctic and Alpine Research, 4(1), 19-38. Nelson, F . E . , and S.I. Outcalt, 1987. A computational method for prediction and regionalization of permafrost, Arctic and Alpine Research, 19, 279-288. Penner, E . 1970. Thermal conductivity of frozen soils, Canadian Journal of Earth Science, 7(3), 982-987. Peterson, T . C . , and R.S. Vose, 1997. A n overview of the Global Historical Climatology Network temperature database, Bulletin of the American Meteorological Society, 78, 28372849. Seppala, M . , 1982. A n experimental study of the formation of palsas, Proceedings of the Jth Canadian Permafrost Conference, 36-42. Slusarchuk, W . A . , and G . H . Watson, 1975. Thermal conductivity of some ice rich permafrost soils, Canadian Geotechnical Journal, 12, 413-424. Smith, M . W . 1975, Microclimatic influences on ground temperatures and permafrost distribution, Mackenzie Delta, Northwest Territories, Canadian Journal of Earth Science, 12, 1421-1438.  85  Smith, M . W . and D . W . Riseborough, 1996. Permafrost monitoring and detection of climate change, Permafrost and Periglacial Processes, 7, 301-309. Smith, C.A.S., C R . Burn, C. Tarnocai, and B . Sproule, 1998. A i r and soil temperature relations along an ecological transect through the permafrost zones of northwestern Canada, Proceedings of the Seventh International Conference on Permafrost, Collection Nordicana 55, 1009-1015. Sturm, M . and J . B . Johnson, 1992. Thermal conductivity measurements of depth hoar. Journal of Geophysical Research, 97(B2), 2129-2139. Taylor, A . E . , M . M . Burgess, A.S. Judge, and V.S. Allen, 1982. Canadian Geothermal Data Collection Northern Wells: Permafrost Temperatures and the Thickness of Permafrost, Earth Physics Branch, Energy, Mines and Resources Canada, Geothermal Series 13, 153pp. Turcotte, D . L . and, G . Schubert, 1982. Geodynamics applications of continuum physics to geological problems, John Wiley and Sons. 450pp. Osterkamp, T . E . , and V . E . Romanovsky, 1999. Evidence for warming and thawing of discontinuous permafrost in Alaska, Permafrost and Periglacial Processes, 10, 17-37. Webb, R.S., C E . Rosenzweig, and E . R . Levine, 1993. Specifying land surface characteristics in general circulation models: soil profile data set and derived water-holding capacities, Global Biogeochemical Cycles, 7, 97-108. Zhang, T., T . E . Osterkamp, and K . Stamnes, 1996. Influence of the depth hoar layer of the seasonal snow cover on the ground thermal regime, Water Resources Research, 32(7), 2075-2086. Zhang, T., and K . Stamnes, 1998. Impact of climatic factors on the active layer and permafrost at Barrow, Alaska, Permafrost and Periglacial Processes, 9, 229-246. Zobler, L . 1986. A world soil file for global climate modeling. NASA Tech. Memo. 87802, N A S A , 33pp.  86  A. Code Appendix  The model developed in this thesis was implemented using M A T L A B , enabled with O D E solver toolboxes. This appendix is not an exhaustive collection of all the code used in the application described in Chapters 3 and 4. Rather the goal of this appendix is to provide enough examples to illuminate the algorithmic approaches. The code is annotated to refer back to relevant equations and chapter sections. The code is divided into two categories: standard and tailored. The standard functions and subroutines are reusable and application independent. The tailored functions and subroutines are specifically adjusted to suit particular applications but make use of the standard functions where appropriate. List of standard functions and subroutines: • set.grid.m • set_time.m • waterice_air.m • Cstar.m • Kstar.m • rstar.m • mix_41ayers.m • topboundary.m • snow_fall.m List of tailored application routines: Testing the Phase Change: • Section3_6_2.m • Section3_6_2_Solver.m Reai World AppHcation - Mackenzie Delta: • RW_MacDelta.m • solve_RW_MacDelta.m • define_models.m  87  % Standard  function:  %set_grid: function  [gc, gf, zc, z f , Dg,  % Refer to F i g u r e  Dz,  dz]  =  set_grid(n,B,H);  (3.1)  % i n p u t argument: % n -> number o f d i s c r e t e c e l l s i n n u m e r i c a l model domain % B -> shape parameter f o r e x p o n e n t i a l g r i d t r a n s f o r m a t i o n % H -> h e i g h t o f model domain i n meters % output % % % % % % % %  exponential  gf_bot Dg  arguments: gc -> a r r a y o f c e l l c e n t r e l o c a t i o n s i n t r a n s f o r m e d space gf -> a r r a y o f c e l l f a c e l o c a t i o n s i n transformed space zc -> a r r a y o f c e l l c e n t r e l o c a t i o n s i n r e a l space z f -> a r r a y o f c e l l f a c e l o c a t i o n s i n r e a l space Dg -> c e l l width i n t r a n s f o r m e d space Dz -> f a c e to f a c e width i n r e a l space dz -> c e n t r e to c e n t r e w i t h i n r e a l space  =  = 1; g f _ t o p =  exp(B);  (gf_top - gf_bot)/n;  g f ( l ) = gf_bot; f o r j=2:n; gf (j) = gf ( j - l ) gf(n+l) = gf_top; gc(l:n) = zf = zc =  + Dg;  end  (gf(2:n+1)+gf(1:n))12 ; .  (H/B)*log(gf) ; (H/B)*log(gc);  Dz(l:n)  = zf(2:end)-zf(l:end-l) ;  dz(2:n) =  (zc(2:end)  - zc(1:end-1));  d z ( l ) = D z ( l ) / 2 ; dz(n+l) = Dz(n)/2;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % %  Standard  subroutine:  set_time.m The format f o r t h i s s u b r o u t i n e i s s t a n d a r d f o r a l l a p p l i c a t i o n s , however the v a l u e s f o r the s t a r t and end time o f the i n t e g r a t i o n i n seconds can be changed to s u i t the problem.  g l o b a l s2y s2d t i n t o u t dt nt % choose a p p r o p r i a t e v a l u e s tin_y =0; tout_y =10; dt = 2e6; % convert s2y = s2d = tin = tout =  y e a r s to seconds 60*60*24*365; 60*60*24; tin_y*s2y; tout_y*s2y;  tspan  f o r a p p l i c a t i o n here: % s t a r t time i n y e a r s % end time i n y e a r s % time increment i n seconds  % s t a r t time i n seconds % end time i n seconds  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Standard subroutine: % % water_ice_air.m % % T h e r m o p h y s i c a l p r o p e r t i e s of l i q u i d water, i c e and a i r . g l o b a l K_W  K_I  K_A  88  g l o b a l C_W C_I C_A g l o b a l r_W r _ I r_A global L Tf % Phase Change Parameters. L = 334.7e3; Tf = 273; % Ice K_I C_I r_I  l a t e n t heat melt temperature  Properties = 2.22; = 1.93e6; = 900.0;  % Water K_W c_W C_W r_W  Properties = 0.602; = 4.186e3; = 4.18e6; = 1000.0;  % A i r Properties K_A = 0.02 5; C_A = 1000; C_A = 1.25e3; r A = 0;  J/kg K  % thermal c o n d u c t i v i t y % volume heat c a p a c i t y % density  W/(m K) J/(m"3 K) kg/m~3  % % % %  thermal c o n d u c t i v i t y s p e c i f i c heat volume heat c a p a c i t y density  W/(m K) J / ( k g K) J/(m"3 K) kg/m^  % % % %  thermal c o n d u c t i v i t y s p e c i f i c heat volume heat c a p a c i t y density  W/(m K) 31(kg K) J/nT3 K) kg/m~3  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Standard function % % % %  %  function:  [heat] =  Cstar(Temp)  F u n c t i o n t o c a l c u l a t e the l a t e n t heat l o a d e d v o l u m e t r i c heat c a p a c i t y . Step change from i c e t o water phase i s a p a r a b o l i c a l l y smoothed H e a v i s i d e . L a t e n t heat i s added as a b e l l f u n c t i o n of u n i t a r e a and base w i d t h Del_TG  % i n p u t argument: Temp - p a s s e d as T(k) i n l o o p from d r i v e r % o u t p u t argument: h e a t - p a s s e d t o C_IW(k) i n d r i v e r global L Tf g l o b a l C_I C_W g l o b a l r _ l r_W  % g r a v i m e t r i c l a t e n t heat, f u s i o n temperature % v o l u m e t r i c heat c a p a c i t i e s of i c e and water % d e n s i t i e s of i c e and water  % L a t e n t Heat - B e l l F u n c t i o n G(T) % % E q u a t i o n s (3.19a - 3.19f) Del_TG Tmm = Tm = Tp = Tpp = AO Al  = Tf Tf Tf Tf  + +  % % % % %  Del_TG/2 Del_TG/4 Del_TG/4 Del_TG/2  Base w i d t h o f b e l l E q u a t i o n (3.19a) E q u a t i o n (3.19b) E q u a t i o n (3.19c) E q u a t i o n (3.19d)  % Equation % Equation  = 16/(Del_TG.~3); = 2/Del TG;  (3.19e) (3.19f)  % % Step change from i c e to water phase - Step f u n c t i o n Del_TH = THm = T f THp = T f B0 =2/  1; - Del_TH/2; + Del_TH/2; (Del_TH' 2) ;  % Equation % Equation % Equation  s  % C o n s t r u c t the b e l l % E q u a t i o n (3.20) i f (TempTm) rho = (r_I+r_W) 12; e l s e i f (Temp>Tp)  function  (3.21a) (3.21b) (3.21c)  G(T)  % Equation  (2.56)  89  function  S(T)  G(T)  rho = r_W; end c l e a r Del__T Tm Tp Temp %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Standard f u n c t i o n : function  [Hs, dH_dt, T s ] = topboundary(t)  % i n p u t argument: % t -> time i n seconds % o u t p u t arguments: % Hs -> snow t h i c k n e s s % dH_dt -> change i n snow t h i c k n e s s % Ts -> temperature a t top domain % required ancillary % snow_fall(t) % required global % s2y -> % s2d -> % Tf ->  boundary  function:  parameters: seconds p e r y e a r seconds p e r day f u s i o n temperature o f water  % r e q u i r e d model parameters: % MAAT -> annual mean temperature % SSTA -> a m p l i t u d e o f s e a s o n a l temperature % NOTES: The snow model chooses J u l y 1 s t as t h e s t a r t o f annual c y c l e . % The s i m p l e s t c h o i c e f o r t h e s u r f a c e temperature s e a s o n a l % v a r i a t i o n i s a s i n u s o i d . A c o s i n e i s used t o match t h e maximum % s u r f a c e temperature w i t h t h e mid-summer s t a r t o f t h e snow % c y c l e , an o f f s e t v a r i a b l e moves the summer maximum. g l o b a l s2y s2d d t g l o b a l MAAT SSTA Hmax % T h i c k n e s s o f snow l a y e r t p = max(0,t+dt/2); tm = max(0,t-dt/2) ; t k = max(0,t-dt) ; Hp = s n o w _ f a l l ( t p ) ; Hm = s n o w _ f a l l ( t m ) ; dH_dt = (Hp-Hm)/(dt); Hs = s n o w _ f a l l ( t ) ; % S u r f a c e temperature o f f = 0; T a i r = MAAT + S S T A * c o s ( 2 * p i * t / s 2 y - o f f ) ; i f (Hs > 0) Ts = T a i r ; else Ts = (3.4 + 0 . 8 9 * ( T a i r - T f ) ) + T f ; end  % E q u a t i o n (2.27)  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Standard f u n c t i o n function[Hs] = snow_fall(t) % T h i s f u n c t i o n c a l c u l a t e s t h e t h i c k n e s s o f t h e s e a s o n a l snow l a y e r % a t time t .  90  % i n p u t argument: t % output arguments:  -> t i m e i n s e c o n d s Hs -> snow d e p t h a t  time  t  % required global parameters: % % s2y : seconds i n 1 year % s2d : seconds i n 1 day % P : days i n 1 y e a r % %  r e q u i r e d m o d e l p a r a m e t e r s ->  previously  assigned  % % % % %  Hmax PI P2 P3 P4  maximum w i n t e r snow d e p t h i n m e t e r s f i r s t d a y o f snow i n d a y s a f t e r J u l y 1 d a t e o f maximum snow d e p t h f i r s t d a y o f s p r i n g snow m e l t l a s t d a y o f snow i n s p r i n g , d a y s a f t e r J u l y  % %  nl n2  shape shape  1  9;  % % %  parameter for parameter for  snow snow  buildup melt  procedure: -> c o n v e r t t h e i n p u t t i m e t o t h e a p p r o p r i a t e J u l i a n -> c a l c u l a t e t h e Hs a t t h a t t i m e s  global global global  day  s2y s2d P n l n2 P i P2 P3 P4 Hmax  % % Convert tj  =  time  t  (seconds)  to  Julian  day i n annual  cycle  rem(t,s2y)/s2d;  % % C a l c u l a t e snow t h i c k n e s s % Equation (2.25)  on day t j  if ( ( t j >= 0) & ( P i >= t j ) ) Hs = 0; elseif ( ( t j > P i ) & (P2 >= t j ) ) Hs = Hmax*( ( t j - P l ) / ( P 2 - P 1 ) )"nl; elseif ( ( t j > P2) & (P3 >= t j ) ) Hs = Hmax; elseif ( ( t j > P3) & (P4 >= t j ) ) Hs = Hmax*( 1 - ( ( t j - P 3 ) / ( P 4 - P 3 ) elseif ( ( t j > P4) & ( s 2 d >= t j ) ) Hs = 0; else disp('problem') end  according  to  %if  t  is  in  late  summer  %if  t  is  in  early  %if  t  is  i n mid winter  %if t )"n2 ) ; %if t  is  in  spring  is  in  e a r l y summer  winter  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Standard function  function:  [C,  Kf,  rp,  pp]  % input arguments: % z _ * -> l o c a t i o n s < <  % % % z  H c  0  % output % C % Kf % rp % pp  of > <  z arguments: -> a r r a y o f -> a r r a y o f -> a r r a y o f -> a r r a y o f  c  = mix_41ayers(z_a,  layers H b  in  real  space  ->|<H -z  z_b,  b  z_c,  -  see  H  a  z_s);  diagram H_S  z  a  b u l k h e a t c a p a c i t y on c e n t r e d g r i d b u l k t h e r m a l c o n d u c t i v i t y on s t a g g e r e d bulk density of centred g r i d p o r o s i t y on c e n t r e d g r i d  91  grid  ->I >l  % r e q u i r e d parameters and v a r i a b l e s : % % *_I/W/A -> t h e r m o p h y s i c a l p r o p e r t i e s o f i c e , water and a i r % % *_a -> t h e r m o p h y s i c a l p r o p e r t i e s o f l a y e r a s o l i d s % *_b -> t h e r m o p h y s i c a l p r o p e r t i e s o f l a y e r b s o l i d s % * _ c -> t h e r m o p h y s i c a l p r o p e r t i e s o f l a y e r c s o l i d s % % C_IW -> temperature dependent h e a t c a p a c i t y o f s o i l water % K_IW -> temperature dependent thermal c o n d u c t i v i t y o f s o i l water % C_IW -> temperature dependent d e n s i t y o f s o i l water  % % % % %  Important n o t e : The s u r f a c e l a y e r ' s ' can ONLY BE SNOW. The p r o p e r t i e s o f C_s,K_s,r_s a r e a s s i g n e d as BULK VALUES and a r e n o t mixed f u r t h e r w i t h t h e temperature dependent water f r a c t i o n s . p_s = 0 ALWAYS.  g l o b a l n Dg gc g f zc z f Dz dz H g l o b a l K_IW C_IW r_IW global global global global  K_W C_W c_W r_W  K_I K_A C_I C_A c _ l c_A r _ I r_A  global global global global  K_a K_b K_c K_s  r _ a C_a r _ b C_b r _ c C_c r _ s C_s  global  p_a f _ a p_b f_b p_c f _ c p_s f _ s  H_s H_a H_b H_c  for j=l:n i f ( (zf (j) < z_a) i_a = j ; end i f ( (zf (j) < z_b) i_b = j ; end i f ( (zf (j) < z_c) i_c = j ; end end  &  (z._a <= z f (j+1))  &  (z._b  &  (z._c <= z f ( j+1) )  < = z f ( j+1) )  for i = l : n - l i f ( z c ( i ) < z_a) & (z_a <= z c ( i + l ) ) j _ a = i+1; end i f ( z c ( i ) < z_b) & (z_b <= z c ( i + l ) ) j_b = i+1; end i f ( z c ( i ) < z_c) & (z_c <= z c ( i + l ) ) j _ c = i+1; end end if if if if if if  (z_a < (z_a > (z_b < (z_b > (z_c < (z_c >  zc(l)) zc(n)) zc(l)) zc(n)) zc(l)) zc(n))  j_a= j_a= j_b = j_b = j_c = j_c =  1; end n+1; end 1; end n+1; end 1; end n+1; end  % h e a t c a p a c i t y combines a r i t h m e t i c a l l y by volume f r a c t i o n s % E q u a t i o n (2.59) % from z=0 t o i f ( i _ c > 1) C(l:i_c-1) rp(l:i_c-l) pp(l:i_C-l) end  z=z_c = (l-p_c)*C_c + f_c*p_c*C_IW(l:i_c-l); = (l-p_c)*r_c + f_c*p_c*r_IW(l:i_c-l); = p_C;  92  % at  z=z_c  dz_c i f  =  z_c  (i_b i f  -  ==  (i_a  zf(i_c);  i_c); ==  dz_b  dz_b  i_c);  =  =  zf(i_c+l)  H_b;  dz_a  =  dz_a  H_a;  =  -  z_c;  zf(i_c+l)  dz_s  =  dz_a -  zf(i_c+l)  =  0;  dz_s  =  0;  z_c; -  z_b;  end  end; C(i_c)  =  (dz_a/Dz(i_c))*(  (l-p_a)*C_a  i f  +  . . .  (dz_c/Dz(i_c))*(  (l-p_c)*C_c  +  f_c*p_c*C_IW(i_c)  )  +  . . .  );  (dz_a/Dz(i_c))*(  (l-p_a)*r_a  +  f_a*p_a*r_IW(i_c)  )  +  . . .  (dz_b/Dz(i_c))*(  (l-p_b)*r_b  +  f_b*p_b*r_IW(i_c)  )  +  . . .  (dz_c/Dz(i_c))*(  (l-p_c)*r_c  +  f_c*p_c*r_IW(i_c)  )  +  . . .  +  . . .  (dz_b/Dz(i_c))*p_b  +  . . .  (dz_c/Dz(i_c))*p_c  ;  z=z_c >  to  z=z_b  i_c)  % at  z=z_b  dz_b  =  C(x_b)  -  == =  = = =  (l-p_b)*C_b (l-p_b)*r_b p_b;  (i_b  z_b  (i_a  >  i_b);  pp(i_b)  =  =  + +  f_b*p_b*C_IW(i_c+l:i_b-l); f_b*p_b*r_IW(i_c+l:i_b-l);  i_c)  zf(i_b);  dz_a  dz_a  =  =  H_a;  zf(i_b+l) dz_s  =  -  z_b;  zf(i_b+l)  dz_s -  =  z_a;  0; end  (dz_a/Dz(i_b))*(  (l-p_a)*C_a  +  f_a*p_a*C_lW(i_b)  (dz_b/Dz(i_b))*(  (l-p_b)*C_b  +  f_b*P_b*C_IW(i_b)  (dz_s/Dz(i_b))*( rp(i_b)  );  (dz_a/Dz(i_c))*p_a  C(i_c+l:i_b-D rp(i_c+l:i_b-l) pp(i_c+l:i_b-l)  i f  . . .  )  =  (i_b  +  f_b*p_b*C_IW(i_c)  (dz_s/Dz(i_c))*(r_s  % from  )  +  =  pp(i_c)  f_a*p_a*c_iw(i_c)  (l-p_b)*C_b  (dz_s/Dz(i_c))*(C_s rp(i_c)  +  (dz_b/Dz(i_c))*(  C_s  );  (dz_a/Dz(i_b))*(  (l-p_a)*r_a  +  f_a*p_a*r_IW(i_b)  (dz_b/Dz(i_b))*(  (l-p_b)*r_b  +  f_b*p_b*r_lW(i_b)  (dz_s/Dz(i_b))*(  r_s  (dz_a/Dz(i_b))*(  p_a  (dz_b/Dz(i_b))*(  p_b  ) + ... ) + .. .  ) + .. . ) + .. .  ); )  +  . . .  );  end % from i f  z=z_b  (i_a  >  to  z=z_a  i_b)  C(i_b+l:i_a-l)  =  (l-p_a)*C_a  +  f_a*p_a*C_IW(i_b+l:i_a-l);  rp(i_b+l:i_a-l)  =  (l-p_a)*r_a  +  f_a*p_a*r_IW(i_b+l:i_a-l);  pp(i_b+l:i_a-l)  =  p_a;  % at  z=z_a  dz_a  =  C(i_a)  (i_a  z_a =  -  >  i_b)  zf(i_a);  dz_s  (dz_a/Dz(i_a))*( (dz_s/Dz(i_a))*(  rp(i_a)  =  (dz_a/Dz(i_a))*( (dz_s/Dz(i_a))*(  pp(i_a)  =  =  zf(i_a+l)-z_a;  (l-p_a)*C_a C_s  (l-p_a)*r_a r_s  +  f_a*p_a*C_lW(i_a)  ) + ...  +  f_a*p_a*r_iw(i_a)  ) + ...  ); );  (dz_a/Dz(i_b))*p_a;  end % from i f  (i_a  z=z_a <  to  z=H  n)  C(i_a+l:n)  % there  exists  % there  i s  a  snowpack.  = C_S;  rp(i_a+l:n)  =  r_s;  pp(i_a+l:n)  =  0;  end i f  (i_a C(n)  == =  n)  no  snowpack  (l-p_a).*C_a  +  f_a*p_a*C_IW(n);  +  f_a*p_a*r_IW(n);  rp(n)  =  (l-p_a).*r_a  pp(n)  =  p_a;  end  % thermal  conductivity  % harmonically % Equation  by  combines  geometrically by  layers.  (2.60)  93  volume  fractions  and  % f r o m z=0 if (j_c >  to  z=z_c  1)  Kf(1:j_c-l) end % at  (K_C~(l-p_c) ) .*(K_IW(1:j_C-l)  =  . ~ ( p _ C * f _ c ) ) . * (K_A  / S  (p_C*(l-f_C) ) ) ;  z=z_c  if  (j_c  elseif if  ==  n+1)  (j_c  ==  (j_b if  ==  (j_a  H_b;  dz_a =  H_a;  1)  dz_c =  dz_c =  H_c;  dz_b =  zc(j_c)  -  z_c;  dz_a  0;  j_c)  dz_b =  H_b;  dz_a =  zc(j_c)  -  z_b;  ==  j _ c ) ,-  z_c  -  z c ( j _ c - l ) ;  dz_a =  H_a;  dz_b =  dz_s =  zc(j_c)  -  z_a;  =  dz_s =  dz_s =  H_s;  0;  end  end else  %  <  j_b  z_c  -  (1  dz_c  =  i f  (j_b  if  <  ==  (j_a  n+1)  z c ( j _ c - l ) ; j_c)  ==  dz_b =  dz_b =  j_c)  H_b;  dz_a =  H_a;  zc(j_c) dz_a =  -  z_c;  zc(j_c)  dz_s =  zc(j_c)  dz_a = -  0;  dz_s =  0;  z_b; -  z_a;  end  end end; Kf(j_c)  =  . . .  dz(j_c)/(dz_a/((K_a~(l-p_a)).*(K_IW(j_c)  (p_a*f_a)).*(K_A"(p_a*(l-f_a))))  dz_b/((K_b~(l-p_b)).*(K_IW(j_c)(p_b*f_b)).*(K_A~(p_b*(l-f_b))))+... dz_c/((K_C~(l-p_C)).*(K_IW(j_c).~(p_C*f_C)).*(K_A"(p_C*(l-f_C))))+... dz_s/K_S % from if  z=z_c  (j_b  >  );  to  z=z_b  j_c)  Kf(j_c+l:j_b-l)  =  (K_b  A  (l-p_b)).*(K_IW(j_c+l:j_b-l).~(p_b*f_b)).*...  (K.A"(p_b*(l-f_b))); % at i f  z=z_b  (j_b  else  %  dz_b if  ==  n+1);  dz_b =  (j_c  <  j__b <  z_b  -  z c ( j _ b - l ) ;  =  (j_a  ==  j_b);  z_b  -  z c ( j _ b - l ) ;  dz_a =  H_a;  dz_s = H_s;  z_b;  dz_s =  n+1)  dz_a =  dz_a = H_a;  zc(j_b)  dz_s =  -  zc(j_b)  -  z_a;  0;  end  end Kf(j_b)  =  dz(j_b)/(dz_a/((K_a Ml-p_a))•*(K_IW(j_b)."(P_a*f_a)).*... (K_A~(p_a*(l-f_a))))+... /  dz_b/((K_b~(l-p_b))•*(K_lW(j_b)."(p_b*f_b)).*... (K_A"(p_b*(l-f_b))))+... dz_s/K_s  );  end % from if  z=z_b  (j_a  >  to  z=z_a  j_b)  Kf(j_b+l:j_a-l)  =  (K_a"(l-p_a)).*(K_IW(j_b+l:j_a-l)."(p_a*f_a)).*...  (K_A~(p_a*(l-f_a))); % at if  z=z_a  (j_a  else dz_a  %  ==  n+1);  dz_a =  (j_b  <  j_a  z_a  -  z c ( j _ a - l ) ;  =  <  z_a  -  z c ( j _ a - l ) ;  dz_s =  H_s;  n+1) dz_s =  zc(j_a)  -  z_a;  end Kf(j_a)  =  dz(j_a)/(dz_a/((K_a"(l-p_a)).*(K_IW(j_a)."(p_a*f_a)).*... (K_A  / S  (p_a* (l-f_a) dz_s/K_s  )))+... );  end % from  z=z_a  (n+1  if  >  to  z=H  j_a)  Kf(j_a+l:n+l)  = K_S;  end  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  94  + . . .  % F e r n Webb  % % % % % % % % %  T h i s program i s t o compare the s o l u t i o n to a freezedown i n a pure water domain g i v e n by the n u m e r i c a l a p p r o a c h u s i n g the l a t e n t heat l o a d e d apparent heat c a p a c i t y v e r s u s the Neumann s o l u t i o n which g i v e s the a n a l y t i c a l s o l u t i o n d i r e c t l y . A d d i t i o n a l l y t h e r e i s a comparison between a r e g u l a r d i s c r e t i z a t i o n of the n u m e r i c a l domain v e r s u s an e x p o n e n t i a l l y d i s t r i b u t e d i r r e g u l a r d i s t r i b u t i o n t h a t c o n c e n t r a t e s c e l l s a t t h e upper ( a c t i v e ) boundary. The e f f e c t o f c o u r s e n i n g t h e r e g u l a r g r i d i s i n v e s t i g a t e d . The lower domain boundary i s a f l a t g r a d i e n t .  % % % % % %  O b t a i n i n g the Neumann s o l u t i o n i s not automatic, as i t r e q u i r e s an i t e r a t i v e p r o c e s s t o f i n d the r o o t o f a t r a n c e n d e n t a l e q u a t i o n dependent on the i n i t i a l and boundary c o n d i t i o n s . T h i s i t e r a t i o n i s coded i n two s t a g e s . The f i r s t i s a c o a r s e r o o t s e a r c h u s i n g g r a p h i c a l methods (not shown). The second i s a r e f i n i n g 'guess and t e s t ' s t a g e which seeks to minimize the d i f f e r e n c e between the l e f t and r i g h t s i d e s o f E q u a t i o n (3.37)  %  global global global global global global global global  s2y s2d n Dg gc Kf C r p K_W K_I C_W C_I c_W c _ I r_W r _ I L Tf  t i n tout dt nt tspan g f z c z f Dz dz HO B pp K_IW C_IW r_IW K_A C_A c_A r_A  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set_time; tspan = [ t i n : d t : t o u t ] ; water_ice_air; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Boundary C o n d i t i o n s G = T_o T_c H =  0; = Tf+2; = Tf-10; 27;  % % % %  Geothermal G r a d i e n t I n i t i a l homogenous temperature S u r f a c e temperature a t t>0 v e r t i c a l domain h e i g h t (m)  Equation Equation Equation  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%*%%%%%  % case  1;  n = 2 0; B = 3; [gc, gf, zc, z f , Dg, T 0 ( l : n ) = T_o;  Dz,  dz]  =  set_grid(n,B,H)  [t,T]=ode45('Section3_6_2_Solver', nt = l e n g t h ( t ) ;  tspan,  TO);  [z_fr] = find_frost_front(T); z _ f r l = H-z_fr; t l = t; z c l = zc; T l = T; % case  2;  n = 50; B = 3; [gc, g f , zc, z f , Dg, T 0 ( l : n ) = T_o;  Dz,  dz]  =  set_grid(n,B,H)  [t,T]=ode45('Section3_6_2_Solver', nt = l e n g t h ( t ) ;  tspan,  TO);  [z_fr] = find_frost_front(T); z_fr2 = H-z_fr; t2 = t ; T2 = T; zc2 = zc; %%%%*%%%%%%%%%%%%%%%%%%%%%%*%%%%%%%% % c a s e 3: n = 5 0; B = 0.1; [gc, gf, zc, z f , Dg, T 0 ( l : n ) = T_o;  Dz,  dz]  =  set_grid(n,B,H)  [t,T] =ode45 (' S e c t i o n 3 _ 6 _ 2 _ S o l v e r ' ,  tspan, TO) ,-  95  (3.3 8c) (3.38d) (3.38b)  nt = l e n g t h ( t ) ; [z_fr] = f i n d _ f r o s t _ f r o n t (T); z_fr3 = H - z _ f r ; t3 = t; zc3 = z c ; T3 = T; %%%%%%%%%%%%%%*%%%%%%*%%%%%%%%%%%%%% % case 4: n = 100; B = 0.1; [gc, gf, z c , z f , T 0 ( l : n ) = T_0;  Dg, Dz, dz]  = set_grid(n,B,H)  [t,T]=ode45('Section3_6_2_Solver', nt = l e n g t h ( t ) ;  tspan, TO);  [z_fr] = f i n d _ f r o s t _ f r o n t ( T ) ; z_fr4 = H - z _ f r ; t4 = t; zc4 = Z C ; T4 = T ; % % % % * % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 5 % case 5; n = 100; B = 3; [gc, gf, z c , z f , T 0 ( l : n ) = T_o;  Dg, Dz, dz]  = set_grid(n,B,H)  [t,T]=ode45('Section3_6_2_Solver', nt = l e n g t h ! t ) ;  tspan,  TO) ;  [z_fr] = f i n d _ f r o s t _ f r o n t ( T ) ; z_fr5 = H - z _ f r ; t5 = t; zc5 = z c ; T5 = T; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Neumann S o l u t i o n % Solve  f o r the  l o c a t i o n of  % x=0  freezing  x=X(t) - solid T_s(x,t) -  % I % % % % %  the m o b i l e  T_c T_s(x,t) T_l(x,t)  x ->  / /  -  front.  / /  liquid / / T_l(x,t) /  oo  / /  T_f  T_o  - T_c = (T_f - T_c) e r f ( x / 2 s q r t ( k _ s t) ) / e r f (N) - T_o = (T_f - T_o) e r f c ( x / 2 s q r t ( k _ l t) ) / e r f c (aN)  N = X/2sqrt(k_s  t)  E q u a t i o n (3.36)  % thermal d i f f u s i v i t y k_I = K _ I / C _ I ; k_W = K_W/C_W; % s p e c i f i c heat c_I = C _ I / r _ I ; c_W = C_W/r_W;  capacity  a = sqrt(k_I)/sqrt(k_W); % iterative  search for transcendental  root N  % 1) choose a v a l u e f o r N N = 17280.795331/100000; % 2) c a l c u l a t e l e f t and r i g h t s i d e s o f E q u a t i o n (3.37) rightside = L.*N.*sqrt(pi)./c_I; leftside = (Tf-T_c).*exp(-(N."2))./erf(N) + (K_W/K_I)*a*(Tf-T_o).*exp(-((a.*N)."2))./erfc(a.*N); % 3) m i n i m i z e the d i f f e r e n c e difn = rightside - leftside  between the  left  and r i g h t  % d e f i n e a s e t o f s o l u t i o n times t_X = l i n s p a c e ( t i n , t o u t , 1 0 0 ) ;  96  sides  % c a l c u l a t e depth of f r o s t X = 2.*N.*sqrt(k_I.*t_X);  f r o n t a t s o l u t i o n times % E q u a t i o n (3.36)  %%%**%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  figure % F i g u r e (3.6) h o l d on,plot(t_X,X,'r-') plot(tl,z_frl,t2,z_fr2,t3,z_fr3,t4,z_fr4,t5,z_fr5) t i t l e ( ' F r o s t Front progression') xlabel('time') y l a b e l ( ' F r o s t F r o n t d e p t h , X, (meters)') % % % % % m % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 4 * % % % % % % % % % % % % % %%%%%%%%%%%%%%%%%%%% FUNCTIONS *%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%m%%%%%%%%%%%%%%%%m%%%%%%%%%%%%%%% % set_grid function  [gc,gf,zc,zf,Dg,Dz,dz]  = set_grid(n,B,H)  g_bot = 1; g_top = exp(B); Dg = (g_top - g _ b o t ) / n ; gf(1) = g_bot; f o r j=2:n; g f ( j ) = g f ( j - l ) + Dg; end; g f ( n + l ) = g_top; gc(l:n)  = (gf (2:n+l)+gf ( l : n ) ) / 2 ;  zf=(H/B)*log(gf); Dz ( l : n ) = z f ( 2 : n + l ) - z f ( l : n ) ,zc=(H/B)*log(gc) ; dz(l)=Dz(1)12;  dz(2:n)=(zc(2:n)-zc(l:n-l));dz(n+l)=Dz(n)12;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % find_frost_front f u n c t i o n [z_fr] = f i n d _ f r o s t _ f r o n t ( T ) z_fr(l:nt)=H; for p=l:nt for i = l : n - l if ( T(p,i) < Tf ) i f ( T ( p , i+1) > T f ) z_fr(p) = z c ( i ) + ( z c ( i + l ) - z c ( i ) ) M T f end end if ( T(p,i) > Tf ) i f ( T(p,i+1) < T f ) z_fr(p) = zc(i) + ( z c ( i + l ) - z c ( i ) ) * ( T f end end end end  - T ( p , i ) ) / (T(p, i+1) - T ( p , i ) ) ;  - T(p,i))/(T(p,i+1)  %%%%%%%%*%*%%%%%*%%*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% f u n c t i o n dT_dt = s o l v e r _ S t e f a n ( t , T ) global global global global  n Dg gc g f zc z f B H L Tf T_c G  g l o b a l K f C K_IW C_IW r_IW T = T' ; % update upper boundary c o n d i t i o n s Ts = T _ C ; % update thermal p r o p e r t i e s : C h a p t e r 3.4 for i = l : n C_IW(i) end K_IW(1) = for i=2:n T_int = K_IW(i) end  = Cstar(T(i) ) ; Kstar(T(l)); (T(i-1)+T(i))12; = Kstar(T_int);  97  - T(p,i));  K_IW(n+l)  = Kstar(Ts);  C = C_IW; Kf = K_IW; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Equation (3.12): dT_dt(l) = (gc(l)/C(1))*(B~2/(H*Dg)"2)*( Kf(2)*gf(2)*(T(2) - T ( l ) ) Kf(1)*G*H*Dg/B ) - ( g c ( 1 ) * l o g ( g c ( 1 ) ) * d H _ d t * G / B * g f ( 1 )  );  % E q u a t i o n (3.8) : dT_dt(2:n-l) = (gc(2:n-l)./C(2:n-l)).*(B"2/(H*Dg)"2).*( Kf(3:n).*... gf (3:n) . * T ( 3 : n ) - (Kf (3:n) . *gf (3:n) + Kf (2 : n - l ) . * . . . gf(2:n-l)).*T(2:n-l) + (Kf(2:n-l).*gf(2:n-l)).*T(1:n-2) )... - (gc(2:n-l).*log(gc(2:n-l)).*dH_dt/(2*H*Dg) ).*(T(3:n) - T ( l : n - 2 ) ) ; % Equation (3.14): dT_dt(n) = (gc ( n ) / C (n) ) * (B' 2/(H*Dg) "2) * ( Kf (n+1) *gf (n+1) * (2 *Ts - T ( n ) ) (Kf(n+1)*gf(n+1) + K f ( n ) * g f ( n ) ) * T ( n ) + ( K f ( n ) * g f ( n ) ) * T ( n - l ) ) ( g c ( n ) * l o g ( g c ( n ) ) * d H _ d t / ( 2 * H * D g ) )*((2*Ts - T(n)) - T ( n - l ) ) ; ,  d T _ d t = [dT_dt] ' ;  98  % Tailored  code  % properties % t i t l e d  which  of  of  {poss_comb}  global  s2y  global  n  s2d  tin  h2y  gc  gf  zc  C  rp  pp  G H  global  Kf  global  K_IW  global  K_w  K_I  K_A  global  C_w  c_l  C_A  C_IW  is  tout  zf  the  models  possible  listed  s o i l  constructed  dt  in  columns in  Figure are  4.18.  held  subroutine  nt  Dz  dz  c_W  c_I  c_A  global  r_W  r_I  r_A  global  L  global  PI  global  K_a  global global  MAAT P3  n2  Hmax  B  SSTA P4  P  nl  r_a  C_a  p_a  K_b  r_b  C_b  p_b  f_b  K_c  r_c  C_c  p_c  f_c  p_s  f_s  global  K_s  r_s  C_s  global  H_s  H_a  H_b  global  poss_comb  f_a  H_c  np  nb  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set_time; tspan  =  [tin:dt:tout];  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% water_ice_air  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n  =  B  =  200; 5;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% define_models;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % i n i t i a l i z e  storage  T_Store(l:n)  =  array  Tf;  TX_Store(l:n)  =  Tf;  TN_Store(1:n)  =  Tf;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for  k=l:np  HO  400;  =  Hmax  =  %  =  poss_comb(k,2)  r_b  =  poss_comb(k,3)  p_b  =  poss_comb(k,4)  K_b  =  poss_comb(k,5)  C_b  =  poss_comb(k,6)  H_b  =  poss_comb(k,7)  H_c  =  H  =  %%%%% for  HO -  -  (H_a  build  the  +  grid  --%%%%%% % i n i t i a l  gf, -  height  H_b  HO;  [gc,  column  poss_comb(k,1  H_a  %%%%%%  solid  zc,  zf,  i n i t i a l  Dg,  Dz,  conditions  dz]  =  on  the  domain  =  MAAT  temperature  (H-zc(i)  end T  =  (no  snow)  set_grid(n,B,H);  i=l:n  T_0(i)  height  T_0;  99  f i e l d  —  in  The an  array  {define_models)  r_IW  global  P2  the  which  Dg  Tf  solves  each  %%%%  %%%%% for for  material  i=l:n; i=l:n;  property  C_IW(i) r_IW(i)  fields  - %%%%%%  = Cstar(T_0(i)); = rstar(T_0(i));  end end  K_IW(1) = K s t a r ( T _ 0 ( 1 ) ) ; for i=2:n; T_int = (T_0(i-1)+T_0(i))/2; K_IW(n+l) = K s t a r ( T _ 0 ( n ) ) ; % locations z_c = H_c; z_b = H_c + z_a = H_c + z_S = H_C +  H_b; H_b + H _ a ; H_b + H _ a + H _ S ;  [C,  pp]  Kf,  of  rp,  internal  K_IW(i)  = Kstar(T_int);  layers  = mix_4layers(z_a,  z_b,  z_c,  z_s);  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ODE S o l v e r [t,T]=ode45('solve_RWl',  tspan,  T) ;  Tsol = T; t s o l = t; nt = length(t) ; [TMX,Ix] [TMN,In]  = =  max(Tsol); min(Tsol);  T X _ S t o r e = [TMX; T X _ S t o r e ] ; T N _ S t o r e = [TMN; T N _ S t o r e ] ; T_Store = [T(end,:); T_Store]; clear  T TMX TMN t s o l  Tsol  end save  solution_output;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Tailored subroutine: % % define_models % % A p p l i c a t i o n to the Mackenzie d e l t a r e g i o n , % Layer thicknesses Zobler, 1986 % S o i l t e x t u r e s f r o m Webb e t a l . , 1993 % M a t e r i a l p r o p e r t i e s s u m m a r i z e d i n T a b l e 4. global global global global global global global global global global global % Snow r_s P_s f_s K_s C_s H s % a)  K_w K _ I K _ A C_w c _ i C _ A c_W c _ I c _ A r_W r _ I r _ A K_a r _ a C_a p_a K_b r _ b C_b p_b K_c r _ c C_c p_c K_s r _ s C_s p_s H_s H_a H_b H_c P i P2 P3 P4 P n l p o s s _ c o m b n p nb  Figure  4.18.  f_a f_b f_c f_s n2  = [175]; = [0]; = [1]; = 10"(2.650*r_s/1000 - 1.652); = (r_s/r_I)*C_I + (l-r_s/r_I)*C_A =[00.2]; % bl  Organic Tundra  see  Layer:  100  % Equation % Equation  (2.61) (2.62)  end  = [1000] ;  P_a f_a K_a C_a H a  = = = = =  [0.85]; [1]; [0.038]; [2.5e6]; [0 0 . 5 ] ;  % b) V a r i a b l e r b P_b f b K b C b H b % c)  = = = = =  Sediment  Layer  [(fine  [1300 1 7 5 0 ] ; [0.45 0.26] ; [1]; [2.1 3.4] ; [1.85e6 1 . 8 1 e 6 ] ; [0.1 3 . 5 ] ;  =  Bottom  r_c P_c f_c K_c C_c  %  = = = = =  b2 silty  %  b3  % % %  b3 b3 b4  % b3  layer  [2300]; [0.10]; [1]; [2.9]; [2.0e6];  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% f o r b l = 0:1 f o r b2 = 0:1 f o r b3 = 0:1 f o r b4 = 0:1  bb = b l + b2*2 + b3*2~2 + b4*2"3; poss_comb(bb+1,1) poss_comb(bb+1,2) p o s s _ c o m b ( b b + 1 , 3) p o s s _ c o m b ( b b + 1 , 4) poss_comb(bb+1,5) poss_comb(bb+1,6) poss_comb(bb+1,7)  H s(bl+1) H a(b2+l) r b(b3+l) p b(b3+l) K b(b3+l) C b(b3+l) H_b(b4+1)  end end end end [np  nb]  =  size(poss_comb)  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %  I n i t i a l  G  =  Conditions  - 0 . 0 2 5;  MAAT  =  SSTA  =20;  %  and Boundary  Tf  -  % Geothermal % Mean  8.3;  % Seasonal  Parameters  P = 365; PI = 77; P2 = P i + P4 = 33 0; P3 = P4 n l = 0.5; n2 = 1 . 5 ;  for  snow  model  % mid  150;  %  -  gradient,  Annual  Surface  Chapter  [K/m]  A i r Temperature, Temperature  [K] Amplitude,  [K]  2.2.1  September  march  % late may 35;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function  dT_dt  = s o l v e _ R W l ( t , T)  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global global  s2y s2d h2y t i n t o u t d t n Dg g c g f z c z f Dz d z  nt  101  global global  K f C r p pp G H B K_IW C_IW r _ I W  global global global global global  K_W K_I K_A C_W C _ I C_A c _ W c _ I c_A r_W r _ I r_A L T f MAAT S S T A  global global global global  K_a K_b K_c K_s  global  H_s H_a H_b H_c  global  P i P2 P3 P4 P n l n2 Hmax  r _ a C_a r _ b C_b r _ c C_c r _ s C_s  p_a f _ a p_b f _ b p_c f _ c p_s f _ s  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T  = T' ;  %  boundary c o n d i t i o n s :  [H_s, %  dH_dt,  grid  Ts] =  topboundary(t);  transformation:  H = H _ C + H_b + H_a + H _ S ; [ g c , g f , z c , z f , Dg, D z , d z ] = %  material property  for for  i=l:n; i=l:n;  C_IW(i) r_IW(i)  set_grid(n,B,H);  fields = C s t a r ( T ( i ) ) ; end = r s t a r ( T ( i ) ) ; end  K_IW(1) = K s t a r ( T ( l ) ) ,f o r i=2:n; T _ i n t = ( T ( i - 1 ) + T ( i ) ) / 2 ; K_lW(n+l) = K s t a r ( T s ) ; % locations of internal z _ c = H_c; z_b = H_c + H_b; z _ a = H_c + H_b + H_a; z _ s = H; [C,  K_IW(i)  = Kstar(T_int);end  layers  K f , r p , pp] = m i x _ 4 l a y e r s ( z _ a ,  z_b, z _ c , z _ s ) ;  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Equation dT_dt(l) =  (3.12): ( g c ( 1 ) / C ( l ) ) * ( B " 2 / ( H * D g ) " 2 ) * ( K f ( 2 ) * g f ( 2 ) * ( T ( 2 ) - T ( l ) ) - ... Kf(1)*G*H*Dg/B ) - ( g c ( 1 ) * l o g ( g c ( 1 ) ) * d H _ d t / ( 2 * H * D g ) ) * ... (G*H*Dg/(B*gf(1)) + T(2) - T ( l ) ) ;  % Equation (3.8): dT_dt(2:n-l) = (gc(2:n-l)./C(2:n-l)).*(B~2/(H*Dg)~2).*( Kf(3:n).*... gf(3:n).*T(3:n) - (Kf(3:n).*gf(3:n) + K f ( 2 : n - l ) . * . . . gf(2:n-l)).*T(2:n-l) + (Kf(2:n-l).*gf(2:n-l)).*T(1:n-2) )... - (gc(2:n-l).*log(gc(2:n-l)).*dH_dt/(2*H*Dg) ).*(T(3:n) -T(l:n-2)); % Equation dT_dt(n) =  dT_dt  =  (3.14) : ( g c ( n ) / C ( n ) ) * ( B " 2 / ( H * D g ) " 2 ) * ( K f ( n + 1 ) * g f ( n + 1 ) * ( 2 * T s - T ( n ) ) - ... (Kf(n+1)*gf(n+1) + K f ( n ) * g f ( n ) ) * T ( n ) + ( K f ( n ) * g f ( n ) ) * T ( n - l ) ) (gc(n)*log(gc(n))*dH_dt/(2*H*Dg) )*((2*Ts - T(n)) - T ( n - l ) ) ;  [dT_dt]';  102  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052394/manifest

Comment

Related Items