A One-Dimensional Temperature Model for a Snow-Cover

SNTHERM is a process driven, one-dimensional energy and mass balanced model, written by Rachel Jordan at the US Army Corps of Engineers, Cold Regions Research & Engineering Laboratory. This model is actively maintained and updated, and the source code is available via the internet. SNTHERM source code is available via FTP from Login as "anonymous" using any password with at least three characters, cd to SNTHERM89.


SNTHERM is a process driven, one-dimensional energy and mass balance model. SNTHERM can take initial snowpack conditions and observed meteorological conditions over a given time period, and using mathmatical equations based on known physical processes, calculate melt and other snowpack fluxes. SNTHERM is essentially an accounting program that distributes energy, mass, and momentum to and through the snowpack and underlying soil as a function of metereological driving variables. SNTHERM utilizes a variety of micro-processes that all affect the energy balance of the snowpack. SNTHERM calculates the surface energy flux of the snowpack for each time step. Energy, mass, and momentum are distributed through the snowpack based on net fluxes of those same variables. SNTHERM grows and melts "layers or nodes" in the snowpack, grows snow grains, increases or decreases density as a function of both compaction and metamorphism, decreases or increases the depth of the snowpack.

It is important to remember that SNTHERM is a one dimensional model. This means that all fluxes generated are at a point. Distributing SNTHERM's output over an area is another matter.

There are 14 parameters in SNTHERM. Parameters within SNTHERM are constant through the period of interest, and include such things as elevation, terrain position (slope, aspect), surface roughness, etc.

There are 20 state variables in SNTHERM. These are things that do change through time, and are usually what we are trying to model, such as snow density, grain size, snow depth, snow temperature, etc. SNTHERM initializes these state variables at the start of the run. That means that you must provide information on these state variables to SNTHERM at the start of each run. SNTHERM then changes these state variables as a function of the driving variables.

Driving variables in SNTHERM are a function of atmospheric forcing, and include air temperature, relative humidity, wind speed (these are for calculating the turbulent fluxes), incoming shortwave radiation, outgoing shortwave radiation, and incoming longwave radiation. SNTHERM can also model each of these radiation variables, but it is better to provide SNTHERM with the measured radiation values whenever possible.

SNTHERM than runs on the same time step as the driving meteorological variables, usually on an hourly basis.

Layers and nodes. "Layers" refers to different media, such as snow and the underlying soil. SNTHERM is normally parameterized with two layers, a snow layer and a soil layer. However, I've also had luck with running SNTHERM on a snow/glacier system. Subdivisions of each layer are called "nodes". "Nodes" correspond to what we would normally call layers in the snowpack. Each node has four components, entered on the same row, in this order:

Snow should have a minimum of ten nodes. The initial node for snow should be about 1 cm in thickness. The second node for snow should be about 2 cm in thickness. The reason for the thin nodes at the snow surface is that this is where the majority of energy, mass, and momentum transfers are taking place. The model tends to crash with thicker nodes because of a lack of convergence. In the example below, note that the first three nodes for snow all have the same temperature, density, and grain size. From the third node on, a thickness of 5-10 cm should usually be selected.

Soil should have a minimum of 5 nodes. Grain size for soils can be set to zero.

Another counter-intuitive aspect of SNTHERM is that the snow and soil properties are read from the bottom up. In the example below, The bottom line is the surface of the snowpack. The top line is the bottom of the soil.

FILENAME contains a list of standard input files that SNTHERM uses and a list of the output files it generates. The file contains the general parameters and initial nodal values of the system. Modeled snowpack fluxes are driven by meteorlogical variables contained in the file You can use any file names that you like. That is, "" can be named "" or "newdata" or anything else that you want. However, each file in FILENAME must have the appropriate format for that type of file. Sntherm reads the file names in FILENAME in order: the file in line one must have general parameters and initial nodal values, the file in line two must have the meteorlogical data, etc.

SNTHERM generates two output files. The file snow.out contains generated snowpack conditions at each time step. The file flux.out contains all modeled snowpack fluxes.

Explanation of old version line 1

(General Parameters)





Number of layers 2 snow and sand
Print-out interval for sntherm.out 1 use same interval as input data
Print-out of heat fluxes, optional 1 yes (what is this file name?)
Barometric pressure 990 use default of 990 mb
Estimate solar radiation 0 no, use values in
--- 1 yes, then must add line 6 to
Estimate incident longwave radiation 0 no, use values in
--- 1 yes, then must add line 6 to
Sloped terrain 0 no
Snow compacted by tanks 0 no tanks here
Near IR extinction coefficient 400 use this default
Measured temp data, optional 0 no (Snow temps)
Water infiltration estimates 0 NO-never turn on
Basic time step in seconds 3600 must be same as in
Estimate standard MET data 0 NO, use
Snow albedo 0.78 default value
Irreducible water saturation for snow 0.04 use default of 4% #2

(Measurement height of meteo rological instruments above GROUND in meters)
1.85,3.9,1.85, height[1..3]




Air temp 1.85 ---
Wind speed 3.9 ---
Relative humidity 1.85 --- #3

Specification of characteristics for each layer type. #4

Convergence related input #5

(Optional)Define your own material codes #6

(Optional)Required data if solar insolation is to be calculated or terrain is sloped. #7

(Optional)Required data for tank tracks.

Explanation of new version line 1

(General Parameters)

2,6,1,990.,1,1,0,0,400.,0,0,3600.,0,0.78,0.04, 999.,100.,1,1,0.9d6




new snow density algorithm 999 A value > 950 turns on
dmlimit: compaction rate 100 100 kg/m3 maximum limit
istboff: stability correction for turbulent fluxes 0 Turns on stability correction
- - - 1 Turns off; same as original version
iqturb: turbulent flux algorithm 1 turns on qturb1: original algorithm
- - - 2 turns on qturb2: new algorithm
vicosity coefficient 090d6 kg-s/m2

Variable albedo . A variable albedo algorithm can be triggered by inputting a snow albedo greater than 1. The routine is an adaptation of that described in the dissertation of Danny Marks (1988), which in turn is based on the work of Marshall and Warren (1987). The main driver in the routine is grain size. The routine does not handle dirty or shallow snow (<10cm) and is limited to clear sky conditions, but does make an adjustment for slope and solar angle.

new snow density algorithm. SNTHERM can estimate the new density based on air temperature, using the Alta function described in Anderson (1976) or in LaChapelle (1961). This option is triggered by inputting a density greater than 950. kg/m3. There are of course variables other than temperature affecting the new snow density (wind speed, crystal type, etc.), so that estimates may be off by as much as 50%

Compaction. Recent testing of the compaction algorithm, taken from Anderson (1976), indicated that adjustments were in order. The settling or destructive metamorphism component was found to significantly overpredict the compaction rate for new snow and the overburden component to compact the snow too slowly. More study is needed of the compaction rate of new snow, but as a 'Band-Aid' it is suggested that the upper snow density limit [dmlimit] for the settling process be lowered from the previously fixed value of 150 kg/m3 to around 100 kg/m3. Dmlimit is now entered in the LAYER.IN file, and the model takes as a cut-off the smaller of this value and the density of the snow at the start of the storm compacted by 15%. Limited testing of the overburden function suggests a lower value of 0.9d6 kg-s/m2 for the viscosity coeffcient, slightly less than the value of 1.0d6 kg-s/m2 used by Brun et al. (1989). The viscosity coefficient [eta0] is now an input. In Jordan (1991), the units on eta0 were incorrectly stated as N-s/m2 and the units on overburden as N/m2 instead of kg/m2. In accordance with Brun et al., I have also changed the exponential scaling factor on density in eq (29) from 0.021 to 0.023 m3/kg. The model had adjusted for mass losses due to melt by decreasing the thickness of nodes (rather than by decreasing snow density). In the heavier stages of melt this procedure contributed to over-densification of the pack, and with the exception of the top node has now been limited to snow less dense than 250 kg/m3. The compaction of wet snow and snow undergoing melt needs further study. In the later stages of the melt period, after the snow has undergone several freeze-thaw cycles, the model over- densifies the snowcover due to repeated freezing of the retained residual water.

istboff stability correction. In the original version, the stability correction for stable atmospheric conditions was turned off, as it appeared to degrade the model's ability to predict surface temperatures. Subsequent testing of the model under melt conditions, at CRREL and elsewhere, have indicated that these default parameters cause the model to over-predict melt. Further study is needed of the turbulent transfer functions, and how they might transition between cold and wet snow conditions. Meanwhile, for melt conditions, it is recommended that the stability correction be turned on. This is done through a new input switch [istboff], where the stability correction is turned off when istboff = 1, and is on when istboff = 0.

iqturb: turbulent transfer module. An alternative turbulent transfer module (qturb2.f) is now availbale. The functions in this module are more rigorous than those in qturb.f, and will significantly affect the results when either the instruments are at different levels or the three roughness lengths are not equal. An input switch [iqturb] has been added, where qturb.f is implemented for iqturb = 1 and qturb2.f is implemented for iqturb = 2. The one caveat on this module is that it has not been thoroughly tested and there may be numerical problems with the unstable case. If the instrument levels and roughness lengths are the same, it is therefore recommended that you stay with the original module.

Don Cline rewrote the output to make it more user friendly. Following is the format for the output files of this version of SNTHERM that you will use:

SNTHERM.89 rev.4 Output Modifications: 1/9/96, by Don Cline

In this version of SNTHERM.89, the output has been completely modified to eliminate all headers, collapse the profile data to weighted mean data, and to generate a few new variables. The output routine that was modified was write.f, however converge.f, subdvide.f, and combo.f also were modified to prevent the printing of convergence and node combining/subdivision messages. Flux.f was also modified to remove headers from the flux file.

The new snow.out and flux files now are straightforward, tabular output files, which may be read more conveniently into graphing software, statistical software, etc. Both files are space-delimited.

The column headers of the new snow.out file are:

1. Decimal time
2. Basestep Number (sequential numbering, 1..N) [ibasestep]
3. Snow depth (m)
4. SWE (m)
5. Mass gain/loss from atmospheric latent heat flux [dlatmass] (kg)
6. Temperature of the top snow node [t(n)] (K)
7. Temperature of the snow/soil interface [tsg] (K)
8. Air temperature at reference height [tkair] (K)
9. Phase of the top snow node [Phase(n)], (1=Thawed, 2=Frozen, 3=Melting)
10. Weighted mean temperature of snow profile [tsnowbar] (K)
11. Weighted mean density of snow profile [bwsnowbar] (kg/m^3)
12. Total depth of liquid water in the snow profile [blsnow] (Units are .01 mm)
13. Weighted mean specific heat of the snowpack [ctsnowbar] (J/kg-K)
14. Weighted mean thermal conductivity of the snowpack [thksnowbar] (W/m-K)
15. Weighted mean snow grain size [dsnowbar] (mm)
16. Grain size of the top snow node [dtopsnow] (mm)
17. Grain size of bottom snow node [dbasesnow] (mm)
18. Weighted mean temperature of soil profile [tsoilbar] (K)
19. Total depth of liquid water in soil profile [blsoil] (m)

The column headers of the new flux file are:

1. Decimal time
2. Year (last two digits)
3. Day of year
4. Hour
5. Change in heat content of the snowpack [dqest](W/m^2)
6. Net solar radiative flux [solar] (W/m^2)
7. Net long wave radiative flux [dlong] (W/m^2)
8. Ground heat flux (across snow/soil interface) [hg] (W/m^2)
9. Sensible heat flux [sensheat] (W/m^2)
10. Latent heat flux [dlatheat] (W/m^2)
11. Precipitation value [prcp] (m/hour)
12. Precipitation type
13. Air temperature at reference height (C)
14. Wind speed at reference height [wsp] (m/s)
15. Relative Humidity at reference height [rh] (%)