User Tools

Site Tools


start:hype_model_description:hype_land

Land routines

This section explains the computations in the land routines of HYPE. If you want an interactive overview of how the routines simulates runoff please have a look at the HYPE Runoff Explorer.

Basic assumptions

A subbasin in HYPE is divided into classes depending on land use, soil type etc. The classes are called SLC:s which stands for Soil type Land use Combination. This division can be compared to hydrological response units. Lakes and river is also classes, but this section consider only land classes.

HYPE model soil routine contains (up to) three soil layers. The number of soil layers and their respective lower limits (soillayerdepth in meters, figure 1) are listed by class in GeoClass.txt. It is possible to have a different number of layers and depth of these different classes.

Figure 1: Soil layer with depth variables (m)

The parameters of water retention in the soil are wilting point (wcwp), field capacity (wcfc) and effective porosity (wcep). These are exclusive and are specified in units of depth (figure 2). The model allocates the water retention capacity evenly between soil layers, depending on their thickness, if only one value set. It is also possible to specify one value per parameter and layer. The model uses water holding capacity in mm for the layers. These are represented in the model code and this document (figure 3) as wp(i), fc(i) and ep(i). The model parameters for the water storage capacity wcfc, wcwp, wcep (and layer depending alternative wcwp1, wcwp2, wcwp3, wcfc1, wcfc2, wcfc3, wcep1, wcep2, wcep3) depend on soil type.

Figure 2: Water retention parameters.
Figure 3: Soil layers water content (mm) and parameters of the water storage capacity (mm).

The soil layer water content (mm) is represented in the document as soil(i), where i = layer, but the value also depends on the subbasin and class (figure 3). The initial value of soil water is set to wp + fc. Optionally the initial value may be set to saturation (wp + fc + ep). The upper soil layer can hold more water than the porevolume. Standing water is not treated as a separate pool.

Tile drainage can be placed in any soil layer (figure 4). The depth of the drainage pipe in meters (tiledepth) is specified in GeoClass.txt for each class. A depth of 0 m is interpreted as drainage pipes are missing.

Figure 4: Drainage pipe and creek bottom in the third soil layer.

Another depth (streamdepth) is specified in GeoClass.txt; this is the maximum depth of the drainage to stream (or ditch). Soil water below this level does not contribute to the local runoff. Note that all land runoff goes through the local river (and possibly local lake) and then the main river before it reaches the outlet lake. There is no direct runoff to lakes.

Overview of flow paths

Figure 5: Illustration of flowpaths in the soil in the HYPE model.

Diagnostic variables

Some additional output variables are calculated from the soil state variables.

Ground water level

A water table is calculated for each soil layer (gwat) from the proportion of water-filled pores of effective porosity. If the ground water reaches above the surface, the water is calculated with 100% porosity.

IF(soil(k)-wp(k)-fc(k)>0.0)
  gwat(k) = (soil(k)-wp(k)-fc(k)-ep(k))/ep(k) * soillayerthick(k) – 
             soillayerdepth(k-1)
IF(gwat(1) > 0) gwat(1) = (soil(1)-wp(1)-fc(1)-ep(1))/1000.

The water table measured as a negative from ground surface to bottom. A positive ground water table means that the land is under water. The lowest soil layer that is not completely filled with water is defined as the “official” ground water table layer and is the one being printed. Soil layers above this may have water in its effective porosity.

Soil moisture deficit

Soil moisture deficit (smdf) is calculated for the root zone, i.e. the upper two soil layers. It is the water (in mm) needed to fill the soil (soil) to field capacity (wp+fc).

smdf=\sum{k=1}{2}{MAX delim{lbrace}{fc_k+wp_k-soil_k,~0}{rbrace}}

Frost depth

Ground frost depth (which is only calculated if the soil temperature is less than zero) depends on soil temperature (soiltemp), but also on soil water content (soil), field capacity (fc), wilting point (wp) and two parameters frost and sfrost.

frostdepth = frost * sfrost * soiltemp * (fc+wp) / soil

There are two parameters in order to be able to choose if you want the frost depth to be land use dependent or soil dependent. The not used parameter is set to one.

Modules (file) Procedures Section
soil_processes (soil_proc.f90) initiate_soil_water_state basic assumptions
initiate_soil_water
calculate_groundwater_table ground water level
calculate_soil_moisture_deficit soil moisture deficit
calculate_frostdepth frost depth

Soil water

Groundwater runoff

Runoff depends on the water table in relation to the drainage depth. Runoff occurs when soil water reaches above field capacity of the soil layers situated above the stream level. Runoff depends on water in the effective porosity (also used to calculate the groundwater table) and a recession coefficient (rc). If the soil is not saturated, runoff from the soil layer depends only on the water of that soil layer. Runoff occurs from all three soil layers (runoff1-3) down to the drainage depth, which is the level of the ditch or creek.

IF(soil(1) > fc(1) + wp(1))
  runoff1 = rc(1) * (soil(1)-wp(1)-fc(1))
IF(soil(2) > fc(2) + wp(2))
  runoff2 = rc(2) * (soil(2)-wp(2)-fc(2))
IF(soil(3) > fc(3) + wp(3))
  runoff3 = rc(3) * (soil(3)-wp(3)-fc(3))

If for example the stream drainage level (streamdepth) is in the third layer (figure 6) and the soil layer is not saturated, the runoff depends on the water level above the stream drainage depth and the following equation replaces the one above.

deltah = (soil(3)-wp(3)-fc(3))/ep(3) * soillayerthick(3) – 
         (soillayerdepth(3) - streamdepth)
IF(deltah>0.)THEN
  runoff3 = rc(3) * deltah / soillayerthick(3) * ep(3)
ELSE
  runoff3 = 0.
ENDIF
Figure 6: Runoff from the third soil layer with a stream.

Similarly the equations change for the other layers if the drainage is in them. Soil layers that lye entirely below the creek level have no groundwater runoff.

If a soil layer is saturated, i.e. soil>=fc+wp+ep, the runoff of the soil layer depends also on the water in the soil layer(s) above. For example if the drainage is in soil layer 3 and both soil layer 2 and 3 are saturated, the groundwater table in soil layer 1 determines the runoff of soil layer 3. The runoff is limited to the water in the layer above field capacity.

deltah = (soil(3)-wp(3)-fc(3)) / ep(3) * soillayerthick(3) – 
         (soildepth(3) - streamdep)
IF(soil(3)-wp(3)-fc(3)-ep(3)>=0.)THEN      !saturated
  deltah2 = soil(2)-wp(2)-fc(2)
  IF(deltah2>0.) deltah = deltah + deltah2 / ep(2) * soillayerthick(2)
  IF(soil(2)-wp(2)-fc(2)-ep(2)>=0.)THEN
    deltah2 = soil(1)-wp(1)-fc(1)
    IF(deltah2>0.) deltah = deltah + deltah2 / ep(1) * soillayerthick(1)
  ENDIF
ENDIF
IF(deltah>0.)THEN
  runoff3 = MIN(soil(3)-wp(3)-fc(3),rc(3)*deltah/ soillayerthick(3)*ep(3))
ELSE
  runoff3 = 0.
ENDIF

The recession coefficient is calculated from two parameters, rrcs1 and rrcs2 which depend on soil type and a parameter that is general rrcs3. The recession coefficient is assumed to decrease with depth and the parameters indicate the coefficient value in the topmost layer (rrcs1) and in the bottom layer (rrcs2). If rrcs2 is not specified it is assumed to be similar to the recession in the topmost layer. The third parameter adjusts the upper layer recession due to the subbasin gradient (slope).

rrcs1 = rrcs1*(1+rrcscorr) + rrcs3*slope

rrcs2 = rrcs2*(1+rrcscorr)

The correction factor corrects rrcscorr parameters rrcs1 and rrcs2 for different parameter regions (parreg). It is defined as an increase. Note that the recession is limited to one. The recession is assumed to diminish exponentially and values of rrcs1 and rrcs2 applies to the midpoint of each layer.

rc(d)=rrcs*e^{-b*d}

The variable b is an auxiliary variable.

b = {log(rrcs1/rrcs2)}/{(soillayerdepth(3) - {soillayerthick(3)}/2.) –
                       {soillayerthick(1)}/2.}

The result is:

rc(1) = rrcs1

rc(2) = rrcs1*exp(-b*({soillayerthick(1)}/2.+{soillayerthick(2)}/2.))

rc(3) = rrcs2

Runoff through drainage pipes

Runoff in the drainage pipes occurs when the water table (the percentage of filled pores of the effective porosity) rises above the pipe's depth (figure 7). Runoff depends on the groundwater surface elevation over the pipe (deltah, m), and a recession coefficient trrcs. Recession parameter trrcs depends on soil type, while drainage pipe level depends on the class. The recession parameter is adjusted with the correction parameter rrcscorr for different parameter regions (parreg). It is defined as an increase.

trrcs = trrcs*(1+rrcscorr)

Depending on which soil layer drainage pipe is in, the runoff will be calculated for water in that soil layer. For the soil layer k (soil(k) is the water content in soil layer k) runoff is calculated as the parameter trrcs times the water found in the effective porosity of the layer and of the overlying soil layers if it is full.

deltah = (soil(k)-wp(k)-fc(k))/ep(k) * soillayerthick(k) -  
         (soillayerdepth(k) - tiledepth)
IF(soil(k)-wp(k)-fc(k)-ep(k)>=0.) deltah= deltah + (soil(k-1)-wp(k-1)-
                                fc(k-1))/ep(k-1) * soillayerthick(k-1)
IF(deltah>0.)
  runoffd = trrcs * deltah / soillayerthick(k) * ep(k)
  
Figure 7: Illustration for calculation of runoff through the drainage pipes.

Diversion of surface runoff and macropore flow

Potential infiltration (infilt) is the sum of rain and snowmelt.

infilt = rainfall + melt

If the potential infiltration is greater than the infiltration capacity, a threshold (mactrinf), and water in the upper layer is larger than another threshold (mactrsm) then macropore flow (macroflow) and surface runoff (infoverflow) occurs. This runoff is caused by an inadequate infiltration capacity of the soil. These two flows are calculated as a percentage (macrate respective srrate) of the potential infiltration and subtracted from this.

macroflow = macrate * (infilt - mactrinf)

infoverflow= srrate * (infilt - mactrinf)

infilt = infilt - macroflow – infoverflow

All the four aforementioned parameters are soil type dependent. If macrate and srrate together are greater than one, they are weighted prior to any calculation are done.

macroflow = (macrate / {macrate + srrate}) * (infilt - mactrinf)

infoverflow= (srrate / {macrate + srrate}) * (infilt - mactrinf)

Infiltration

The base model for infiltration is described above in section Diversion of surface runoff and macropore flow. An alternative model for infiltration consideres the effect of frozen soil. It is based on Zhao and Gray (1999). This model redirects all or part of the infiltration, calculated with the base model, to macropore flow and overland flow proportionally to their respective model parameters (macrate and srrate).

If the minimum daily temperature is less than 10 degrees and the infiltration is larger than 5mm/d an icelens is created. In this case, and as long as the maximum daily temperature is below zero, the icelens redirect all infiltration to surface runoff and macropore flow.

If there is no icelens, but the soil temperature of the upper soil layer (soiltemp) is below zero the infiltration is restricted. The potential infiltration (potinfilt) depends on a model parameter (bfroznsoil) that is soil type dependent. It also depend on the “opportunity time” (t0), which is an estimate of the time with possible infiltration in hours.

t0= {lbrace}{
 \matrix{2}{2}
    {1 {for\ 0.65*snow<6}
    {0.65*snow-5} {for\ 0.65*snow>6} }}

potinfilt= {lbrace}{
 \matrix{2}{2}
    {0 {for\ soil>=pw}
    {bfroznsoil*{0.99}^{2.92}*({1-{soil}/{pw}})^{1.64}*({{0-{soiltemp}}/{273.15}})^{-0.45}*{t0}^{0.44}}/{t0/24} {for\ soil<pw} }}

where snow is snow water equivalent, soil is soil water of the upper soil layer, pw is pore volume (in mm) and calculated from the model parameters for water holding capacity (wcwp, wcfc and wcep).

If the infiltration is greater than the potential infiltration it is then reduced to the potential infiltration.

If the calculated infiltration is greater than zero, it is added to the upper layer soil water. This is done regardless of whether there is space in the soil pores there. If the water exceeds the water pore volume it is assumed to lie on the ground, but it still belongs to the upper soil layer and is mixed and thus has the same concentrations.

Percolation

The flow of water downward through the soil layers is only done by water over field capacity (water in the effective porosity). A maximum percolation (mm/d) limits the flow between soil layers. For the upper soil layer it is mperc1, and for the second soil layer it is mperc2. These parameters are soil type dependent. Flow is also limited by how much water the lower layer can receive.

Drainage from soil layer 1 to soil layer 2 is

perc1x = MIN((soil(1)-wp(1)-fc(1)),mperc1)

but if there is not enough capacity in soil layer 2 the drainage is instead

perc1 = wp(2)+fc(2)+ep(2)-soil(2)+perc2

and fills the second soil layer.

Drainage from soil layer 2 to soil layer 3 can be at most

perc2x = MIN(wp(3)+fc(3)+ep(3)-soil(3),mperc2)

for that is what the soil layer 3 can receive. If soil layer 2 does not reach field capacity with perc1x added

perc2 = 0

but if soil layer 2 with perc1x added exceeds field capacity

perc2 = MIN(soil(2)+perc1-wp(2)-fc(2),perc2x)

Regional groundwater flow is created by additional percolation from soil (see Section on Regional groundwater flow).

Upwelling

Flow may enter the lowest soil layer, i.e. regional groundwater flow. Upwelling to soil layers above may occur if the soil layer is filled.

Saturated surface runoff

Surface runoff due to a high ground water table (satoverflow) occurs when the water table in the upper soil layer reaches above the surface. It depends on a parameter srrcs which is dependent on land use. The recession parameter is corrected with the correction factor rrcscorr for different parameter regions (parreg). It is defined as an increase.

srrcs = srrcs*(1+rrcscorr)

satoverflow = MAX(srrcs * (soil(1)-wp(1)-fc(1)-ep(1)),0.)

Runoff is removed from the uppermost soil layer. The total surface runoff (due to high ground water table and low infiltration capacity) is calculated and printed.

Macropore flow

Macropore flow occurs when the potential infiltration and water in the upper soil layer is large (see Diversion of surface runoff and macropore flow above). It is caused by a limited infiltration capacity of the soil. Macropore flow (macroflow) is added to the layer in which the water table is located (see Diagnostic variables above). The water is added to this layer only until it is full and the excess is trapped in the layer above.

Thus water does not flow up into the layer above when macropore flow is larger than the empty space in the soil layer with the water table, as in the case of groundwater inflow. Instead the excess flow stays in the soil layer above before reaching the soil layer of the water table. This distinction is important for the substances following the macropore flow.

References

Zhao, L., and D.M. Gray 1999. Estimating snowmelt infiltration into frozen soils, Hydrological Processes, 13:1827-1842.

Modules (file) Procedures Section
soil_processes (soil_proc.f90) initiate_soil_water groundwater runoff
calculate_soil_runoff
calculate_tile_drainage runoff thorugh drainage pipes
infiltration diversion of surface runoff and macropore flow, infiltration
percolation percolation
add_macropore_flow macropore flow
general_water_concentration (general_wc.f90) inflow_lowest_soillayer upwelling

Snow routines

Snow melt

As default, snow melting occurs when the temperature is greater than the threshold temperature. The amount of snow (in mm) that melts (melt) depends on the snowmelt parameter cmlt, threshold temperature parameter ttmp and air temperature (temp).

melt = MIN(cmlt*(temp-ttmp),snow)

The parameters cmlt and ttmp are related to land use. The parameter ttpi is a general parameter, but was in earlier versions always equal to one. The parameter ttpd is general.

Snow melt models

Alternative snowmelt models exist, but are not fully described here yet.

Model 0 (default)

Temperature index model, without snow cover scaling

melt = MIN(cmlt*(temp-ttmp),snow)

Model 1

Temperature index model, with snow cover scaling

melt = MIN(cmlt*(temp-ttmp),snow)*snowcover

Model 2

Temperature and radiation index model, with snow cover scaling

Snow cover

Normally snow is assumed to cover the whole class if present. Alternatively if parameters are given, snow cover fraction (fsc) within a class is calculated based on snow water equivalent (snow). The formulation is based on Samuelsson et al. (2006). During snow build up the snow cover increase as a function of snow water equivalent until a maximum value (general parameter fscmax) is reached.

fsc=fscmax×tanh(0.1×snow)

It is also possible to specify a minimum snow cover (general parameter fscmin). As soon as the fractional snow cover area reaches above a certain threshold (fscmax-fsclim), the snow cover area is determined by another relation that represents the redistribution of snow during winter. In this case snow cover is dependent on maximum snow pack during the winter (snowmax) and a snow cover redistribution factor that is dependent on variation in elevation (stdelev, the standard deviation of elevation within the subbasin) and land use.

fsc= snow / {snowmax×fscdist}

fscdist=fscdist0+fscdist1×stdelev

The snow distribution factor (fscdist) is determined by three land use dependent parameters; fscdist0 and fscdist1 in the linear equation and a maximum value (fscdistmax). Also in this case the snowcover is limited by the maximum and minimum value parameters. When the end of the snow season approaches (defined by general parameter fsck1) the snowmax variable is gradually decreased in order to be reset before next winter season:

snowmax=snowmax-(fsck1×snowmax-snow)×{1-e^{-fsckexp×ts}}/{fsck1}

{snow}/{snowmax}<fsck1

The equation depends on two general parameters, fsck1 and fsckexp, where fsckexp depend on time (ts is seconds per timestep of simulation).

For winters when the snow pack not reach the definition of large snow pack, the first equation is used during the whole season.

Soil temperature and snow depth

Soil layer temperature (soiltemp) is calculated as a weight of three temperatures; previous time step soil layer temperature, soil temperature at deep depth (deeptemp) and air temperature (temp). The weight of the deep soil is constant (0.001), while the weight of the air temperature (weightair) depends on snow depth (snowdepth) and parameters. The soil memory (soilmem) depends on depth and land use, with parameters surfmem and depthrel. The memory of deep soil temperature is a general parameter (deepmem).

soilmem = {lbrace}{
 \matrix{2}{2}
    {deepmem {for\ deeptemp}
    surfmem×e^{-depthrel×depth(k)} {for\ soil\ layer\ k} }}

weigth_{air}={1}/{soilmem+10×snowdepth}

deeptemp=weight_{air}*temp+(1-weight_{air})×deeptemp

soiltemp=weight_{air}*temp+(1-weight_{air}-weight_{deep} )×soiltemp+weight_{deep}×deeptemp

In the default snow depth model, snow density (snowdens) depends on the snow's age in days (snowage). Snow density for fresh snow (sdnsnew) and the increase of density with snow age (snowdensdt) are general parameters (~ 0.1 and ~0.002). The snow's age increases by one every time step, but are weighted with age (0) for any new snow.

snowage = (snowage + 1)*oldsnow/(oldsnow+snowfall)

snowdens = sdnsnew + snowdensdt * snowage

snowdepth = 0.1 * snow/snowdens

In the alternative snow depth model, snow density (snowdens) is calculated by a compacting factor. Snow density for fresh snow (sdnsnew), maximum snow density (sdnsmax), compactation rate for low temperatures (sdnsrate) and additional compactation for high temperature (sdnsradd) are all general parameters. The change in snowdensity (densdt) due to compactation each time step is calculated as:

densdt = sdnsrate*(sdnsmax-snowdens)

for cold days (temperature is below threshold temperature parameter ttmp), and

densdt = (sdnsrate+sdnsradd)*(sdnsmax-snowdens)

for warm days.

Modules (file) Procedures Section
soil_processes (soil_proc.f90) calculate_snow snow melt
calculate_snowmelt
snowalbedo_function
calculate_fractional_snowcover snow cover
calculate_snowdepth soil temperature and snow depth
calculate_soiltemp
calculate_weigthed_temperature

Glaciers

The glacial class' area ( class_{area} ) is divided into snowfield and glacier ( glac_{area} ). The division may vary with simulation time, but the glacier cannot be larger than the slc class' area. At simulation start the area is considered to be all glacier, unless specific information on which year the class area is valid is given. The initial glacier volume is calculated based on glacier area or set for each glacier. If the initial year is given the initial volume is valid for that date. In this case, year of the area and/or volume and the average annual change in massbalans (annmb in m/yr) is given as input data or model parameters. The initial glacier volume is then adjusted according to the difference in years (yeardiff) between simulations start and glacier area data.

The glacier volume - area relationship:

glac_{vol} = coef*{glac_{area}}^{exp}

coef = glacvcoef*{e}^{c}

The initial glacier volume if calculated from class area:

glac_{vol} = coef*{class_{area}}^{exp} + yeardiff*{annmb*class_{area}/{glacdens*1000}}

The glacier area is correspondingly calculated from the glacier ice volume, but is limited by the class area. The relation is:

glac_{area} = ({glac_{vol}*{1/{coef}}})^{1/exp}

The equation coefficients coef and exp can have different values for specific glaciers. The first coefficient coef is calculated as the product of EXP( c), where c is a glacier volume correction and a general parameter (glacvcoef/glacvcoef1) depending on glacier type. The second coefficient, exp is a general parameter (glacvexp/glacvexp1) depending on glacier type. Glacier density (glacdens) is a general model parameter (m3 water / m3 ice).

Glaciers are divided into four types. The default type is (a small) glacier, the alternatives are ice cap (large glacier), ice sheet and infinite glacier. Glacier type is given as input, or determined by the glacier area (a threshold (glac2arlim, a general parameter). The glacier area is used to determine the glacier type if it is not given as input and the threshold parameter is set. The glaciers will then be divided into (small) glaciers and ice caps.

The default glacier type is (a small) glacier, while the ice cap represent larger glaciers. The difference between the two are in the parameter völues for the volume-area relationship. The default values for the volume-area relationship are based on Radic and Hock (2010).

The ice sheet glacier do not simulate an snow field, and the glacier will have a constant area. The ice sheet glacier type is thus not depending on the volume-area relationship, except for initial volume if that is not given. In that case the initial volume of the ice sheet glacier will used the ice cap parameter values.

The infinite glacier type is independent on glacier volume, because the glacier melting is not limited by the glacier volume. It has a constant area (the class area) and do not simulate an snow field.

For snowfields, all soil processes are calculated as for the common land classes. Snow pack is only calculated for the snowfield, i.e. snow depth is assumed zero on glaciers. For glaciers, all precipitation is added to the glacier ice. Glacier melting is calculated as:

melt = cmlt*(temp-ttmp)

where cmlt is the general parameter glaccmlt and ttmp the general parameter glacttmp.

An alternative glacier melting model exist, which depend on temperature (temp) and radiation (swrad) and has a refreezing component (goverened by general parameter crefr). The radiation component depend on the albedo of the glacier, which in turn depends on if the glacier is snow covered (snowcov) or not, and general parameters for rate (cmrad) and glacier ice albedo (glacalb). The refreezing component is limited to the glacier melt.

melt= {lbrace}{
 \matrix{2}{2}
    {cmrad*swrad*(1-albedo_glacier)+cmlt*(temp-ttmp) {for\ temp>=ttmp}
    cmrad*swrad*(1-albedo_glacier)-crefr*cmlt*(ttmp-temp) {for\ temp<ttmp} }}

albedo_{glacier}=albedo_snow*snowcov+glacalb*(1-snowcov)

Glacier evaporation is optional and set by an model option (snowevaporation). The evaporation is determined as a fraction of the potential evaporation by the general model parameter fepotglac.

Glacier melting, snow melting on snowfield and/or rain are added to the soil. Thus the soil represents the whole glacier class area. The concentrations in the glacial melt water are zero. This means that any atmospheric deposition of nutrients is lost on the glacier.

Modules (file) Procedures
glacier_soilmodel (glacier_soilmodel.f90) soilmodel_3
initiate_glacier
initiate_glacier_state
calculate_glacier_area
get_glacier_parameters

References

Samuelsson, P., S. Gollvik, and A. Ullerstig, 2006. The land-surface scheme of the Rossby Centre regional atmospheric climate model (RCA3), SMHI Report Meteorologi Nr 122, 25 pp.

Radic, V. and Hock, R., 2010. Regional and global volumes of glaciers derived from statisti-cal upscaling of glacier inventory data, J. Geophys. Res., 115, F01010, doi:10.1029/2009JF001373.

start/hype_model_description/hype_land.txt · Last modified: 2017/11/09 07:33 by cpers