Mathematical modeling of temporal changes in snow-firn properties in the cold season

Abstract

The paper describes a non-linear dynamic thermomechanical model of the snow–firn development in the high mountain environment during the cold season. The model allows the estimation of the thermal and mechanical development of snow in an inhomogeneous stratified snowpack and in the upper part of the firn layer, the development of the snow–firn structure and transformation of snow into firn. The non-linear mathematical model consists of a heat problem, water-vapor diffusion problem, densification of the snow–firn layer, structural equations and strength conditions of the metamorphosed snow–firn depending on its temperature, structural parameters, density and rate of sublimation. The model includes conditions of snow transformation into firn and its structural changes. These equations are non-linear and coupled. The iterative finite-difference numerical method was used for the calculations and involved a full mathematical model. Numerical experiments were implemented using the complete model as well as real. The numerical computations were made for various meteorological data, temperature gradients and in the wide range of structural parameters. The calculation results were found to be in good agreement with the results of observations.

List of Parameters

Introduction

Dry snow is a mixture composed of ice grains connected with bonds and a gas component (dry air and water vapor). Alley and others (1982) analyzed the sphericity of grains as a function of depth for coarse-grained and fine-grained firn and found that the sphericity idealization approach is better for fine grains than for coarse firn. The symmetrical packing of equal sphere-grains connected by cylindrical bonds was chosen as the idealization of the structure model, The existence of periodical and ordered disposition of ice crystals with average number of bonds per grain i _k = 4 in the density range 240–506 kg m⁻³ is shown by Samoylyuk (1992). The good agreement of the netlike model with measurements of polycrystalline ice and snow is shown by Watanabe (1993). The porous ice matrix can be approximated by spheres with average radius r connected by cylindrical bonds with lengths l _w, l _h (vertical and horizontal direction) and diameters d _w , d _h. The snow consists of layers with different structures and properties. Weather changes produce layers with different thermomechanical and structural properties. The influence of the thermal conditions and variations in air humidity is significant in the upper part of the snow–firn cover. Snow is transformed into firn through densification and recrystallization processes. The temperature gradients define the mass transfer and modification of the structure. With depth increase, firn transformation is practically isothermal. The structure difference defines transformation of the thermomechanical properties of the snow-firn at the layer boundaries. Heat and mass transfer in the snow-firn is the heat conduction on the ice matrix with evaporation/condensation and diffusion of the water vapor on the pore space. Concentration of the water vapor is a function of the temperature and curvature of the grains and bond surface. The processes of the volume, surface and grain-boundary diffusion as well as transport of mass by evaporation and condensation lead to recrystallization of the structure. The mass transfer disturbs the equilibrium state of the water vapor in porous space. This process leads to the condensation on the surface of the coarse grains and to the evaporation of the fine grains and bonds. The rate of growth of the grains is proportional to the degree of super-saturation. The vertical diameter of bonds increases through the densification of snow and plastic flow of the bonds. Anisotropic changes in the vertical and horizontal directions result from the viscous flow and evaporation of the bonds.

Model Description

A one-dimensional non-uniform snowpack is considered. The mathematical model for the non-uniform stratified snow cover was investigated by Golubev and Guseva (1987) and Guseva and Golubev (1989, 1990). The essential difference between this model and other models is the possibility of calculating the sublimation/evaporation of the different elements of the structure and estimating the changes in diameter due to recrystallization. Compared with previous models which allow estimation of the phase-changing in the major balance (e.g. Bader and Weilenmann, 1992; Gray and others, 1995), the mathematical model is considerably extended through its description of the influence of wind-pumping on structural parameters and on densification. The origin of the vertical axis x = 0 is defined at a depth H(t) from the surface. The temperature can be described with the heat balance equation:

(1)

where the coefficients of heat conductivity and capacity of the snow K _s, C _s linearly depend on snow–firn density ? _S and f(x) is distribution of absorbing sun radiation, j = j ^diff + j ^conv (j ^diff is the sublimation rate through the diffusion of the vapor, j ^conv is the convective part through wind-pumping). The speed V _a is defined by Albert (1993) as Darcy speed of the airflow.

The boundary conditions for Equation (1) are as follows: (1) At the snow surface, x = H(t) was applied as the condition for the heat flux which is defined through the turbulent heat exchange and through the heat for evaporation from the snow surface:

(2) At the depth of the snow–firn x = 0, the condition ?T/?x = 0/_{x = 0}.

(3) If the origin of the vertical axis is defined at soil level, then Stephan’s condition of heat balance is defined on the phase boundary in the soil (melted or frozen state) (Guseva and Golubev, 1989). Initial conditions for the temperature are as follows: T(x) = f ₁(t). The water-vapor pressure e _v above elements of structure is defined by the Magnus relation:

(2)

where r' is radius of the surface curvature of the elements of the ice matrix. This radius is different for the grain, vertical or horizontal bonds. The rate of the growth of the radii of the ice-matrix elements is given by the Arrhenius equation (Stephenson, 1967):

where r', r ₀ are the current and initial curvature radii. The rate of sublimation/evaporation is defined through the processes of diffusion and convective flux (wind-pumping) in the porous space:

(3)

where j = j ⁺ + j ⁻ is the sum of the positive part of the mass transfer (corresponds to condensation) and negative part of the mass transfer (corresponds to evaporation) above the common sublimation/evaporation surface S _f of the grains and bonds per one unit of volume.

The boundary conditions for Equation (3) are as follows: (1) At the surface x = H(t), the condition for the water-vapor flux (air–snow-cover)

was applied.

(2) At the depth x = 0 the condition for the water-vapor flux

was used.

(3) If the origin of the vertical axis is defined at soil level, then the condition for the water vapor is defined on the snow–soil boundary (Guseva and Golubev, 1989).

The snow-cover thickness depends on the snow density variation and the quantity of precipitation. The densification can be described as follows (Yosida, 1963; Voitkovskiy, 1977):

(4)

where ? _a is the additional mechanical wind pressure and ? _k exponentially depends on the snow density (Yosida, 1963). Using the coefficients of the structure rigidity b _W = 0.5 d _w /r, b _h = 0.5 d _h /r and the coefficients of the structure friability k _w = (2r + l _w)/2r, k _h = 0.5(2r + l _h)/r and average coefficients b _m = 0.5(b _v + b _h), k _m = 0.5(k _w + k _h), the number of bonds per grain, the number of grains per unit volume and the Poisson Coefficient for snow can be written as (Voitkovskiy and Golubev, 1987):

(5)

(6)

where C is the parameter depending on the average number of bonds per grain i _k. Using the Arrhenius relation for the ice-matrix elements growth in snow and ice, Darcy speed of the airflow and Magnus condition, it is possible to write the conditions for rate of sublimation/evaporation on the surface of the grain j _z , on the surface of the vertical bonds j _w, and surface of the horizontal bonds j _h, as follows:

(8)

where n _w + n _h, = 1, the sublimation surface of the vertical bonds can be written as S _w = pd _w l _w and sublimation surface of the horizontal bonds as S _h = pd _h l _h; the sublimation surface of the grain as S _z = 4pr ² - S _w(i _k/2)n _w – S _h(i _k/2)n _h; and the water-vapor pressure above horizontal bonds e _h, above vertical bonds e _w and above grains e _zare described by Equation (2). Here we use the hypothesis that the sublimation rate on the surface of ice-matrix elements is the result of the equalizing of the difference between average water-vapor pressure in pore space and concentration above the structure elements and sublimation/evaporation through wind-pumping.

The variation of bond diameters as a result of all processes — sublimation, condensation, viscous flow — can be written as follows:

(9a)

(9b)

The change in grain diameter is a result of the growth of coarse grains through evaporation of fine parts of the ice matrix (Magnus relation) and changes of the rate of the evaporation condensation (Equations (8)):

(9c)

For the viscous flow of vertical and horizontal bonds, the following conditions apply:

(10)

where K _p is constant, P _s ² is constant load (the relative plastic deformation of cylindric model with uniaxial compression stress and possible lateral extension under constant load; Guseva and Golubev, 1990) and P _a is an additional mechanical wind-pumping pressure on the surface of the snow pack.

The vertical relative contact section as a function of the structural parameters can be written as follows (Golubev, 1982):

(11)

where A and a are defined by the main directions of the ice matrix and the average number of bonds. For horizontal relative contact surface the equation is analogous to Equation (11).

The initial conditions are as follows:

(12)

The system of Equations (1)–(12) is a full system, and it is the basis for the calculation of the snow strength as a function of structural parameters. The equations for snow strength are as follows:

The system of Equations (1)–(12) consists of coupled partial differential equations with boundary and initial conditions. The mathematical model proposes the modeling of the temporal changing of the main snowpack parameters, such as temperature, water-vapor pressure in the snow cover, density profiles, structure and thermomechanical properties for the mountain conditions and for the upper snow–firn layers of the glacier. This model allows the calculation of the sublimation/condensation rate. The humidity of air in the pore space is defined by the temperature and snow density. The heat and mass transfer leads to the evaporation and condensation of the elements of the ice matrix. The growth of the average grains and elements of the ice matrix is proportional to the supersaturation of gas inside the snowpack.The character of this growth strongly depends on the temperature and its gradients, density and curvature of the matrix elements. In the boundaries between different layers (contact grains with different radius) and in the snow-ground or snow–firn boundary, mass transfer process is essential. The influence of the temporal structure changes on the temperature and density profiles is insignificant. But even small variations in the temperature, its gradients and density practically define the temporal changes of the structural parameters. The temperature change and its gradients stipulate variation of the mass transfer in the snow cover. Inside the snowpack a boundary between the depth hoar and the densification crust is formed. At the layer boundaries, narrow zones develop where sublimation or condensation dominate. The processes in the upper part of the snow cover are defined by the variations of meteorological parameters including wind-pumping.

Numerical Simulations

The mathematical model consists of the following problems: the heat problem (calculation of temperature); the wind-pumping problem (calculation of the vertical component of the wind speed through the snowpack); the equilibrium water-vapor pressure, which can be calculated by the Magnus relation; the water-vapor mass-transfer problem (calculation of the intensity of mass transfer and evaporation/condensation rate); the densification problem (calculation of the snow–firn density profile); the equations for structure parameter changes (calculation of the grain radius, number of grains per unit volume, vertical and horizontal bonds diameters, number of bonds per average grain in every layer, contact sections); and strength profiles equations.

The problem of Equations (1)–(12) is solved numerically. A one-dimensional, finite-difference balance iterative numerical method (Samarskiy, 1983) similar to permafrost processes has been used (Grigoryan and others, 1987). This approach makes it possible to obtain the vertical profiles of different non-linear parameters, which vary with time. The program allows for calculations with different grid size (more than 0.03 m). It can be used for different calculation periods and different meteorological and initial data. It requires as initial data the vertical profiles of temperature, density and ice-matrix structure parameters, which are as follows: the numbers of grains per unit volume, number of bonds per grain, grain radius, diameters and lengths of bonds. The calculations also require the following meteorological data: temporal temperature of the air near the surface; air humidity; solar radiation and albedo; quantity and density of precipitation and its structure parameters; and for the wind-pumping problem the value of the additional wind load, the air pressure on the surface, the wind velocity at the surface, and frequency and amplitude of wind pulsations must be defined. These are the standard meteorological data with some additional information on snow structure and wind conditions. If soil were included in the model, the soil heat parameters and density would be needed in order to perform the calculations. The synoptic parameters would have to be input into the model every 12 h.

The results of the calculations obtained by the dynamic thermomechanical model of the formation of the structure, thermal and mechanical properties of snow–firn in a non-uniform layered stratified snowpack are as follows: temperature, density, rate of sublimation/evaporation, structure parameters (grain radius, relation of the bonds diameters, grain number per unit volume, bond number per grain, structure parameters in vertical and horizontal directions, relative bond diameters, contact sections, strengths). Numerical experiments were performed using the complete model and real examples for mountain (Golubev and Guseva, 1987; Guseva and Golubev, 1989, 1990) and glacier conditions. The numerical computations were made for various meteorological data and temperature gradients as well as in a wide range of structural parameters. The calculation results agree with observations.

Conclusion

A mathematical model for the monitoring of the main properties and changes of the snow–firn has been proposed. With this model it is possible to predict the variation and development of these parameters in different weather scenarios. It is clear that more detailed testing of the model for the various initial and weather conditions is necessary. It would be very interesting to compare calculations of some physical situations with detailed observations of the different conditions.

The model of heal and mass transfer enables the estimation of the time and place of depth-hoar layer formation and calculation of the snow strength in different layers of snow cover and its metamorphism variation. The mathematical model and the program have a practical significance as a tool for studying the density, temperature, and structure changes in the snow–firn inhomogeneous layer. The results of the mathematical model can enhance understanding of the radar-scattering results and its models. They can help to explain the complex radar-altimeter waveforms which depend on profiles of the density, temperature and structure parameters. Application of satellite microwave measurements would provide spatial and temporal information on the structure and morphology of polar ice and snow, which is significant for the modeling of ice-sheet dynamics (Rott and others, 1993). The results of the numerical modeling can be helpful in studying the capabilities of detecting small-scale variations of snow and ice morphology. The temporal profiles of the structure, temperature and mechanical property changes can be obtained. These calculations can be calibrated by comparison with the results of the airborne radar altimetry.

The model can be used for prediction and glacio-climatic monitoring of the changes in the snow–firn layer for different synoptic conditions.