Skip to content

Overview of Numerical Methods

1. Vertical flow

Unsaturated flow is one of the central processes in MIKE SHE and in most model applications. The unsaturated zone is usually heterogeneous and characterized by cyclic fluctuations in the soil moisture as water is replenished by rainfall and removed by evapotranspiration and recharge to the groundwater table. Unsaturated flow is primarily vertical since gravity plays the major role during infiltration. Therefore, unsaturated flow in MIKE SHE is calculated only vertically in one-dimension, which is sufficient for most applications. However, this may limit the validity of the flow description in some situations, such as on steep hill slopes with contrasting soil properties in the soil profile. MIKE SHE includes an iterative coupling procedure between the unsaturated zone and the saturated zone to compute the correct soil moisture and the water table dynamics in the lower part of the soil profile

There are three options in MIKE SHE for calculating vertical flow in the unsaturated zone:

Richards Equation

The full Richards equation requires a tabular or functional relationship for both the moisture-retention curve and the effective conductivity.

The full Richards equation is the most computationally intensive, but also the most accurate when the unsaturated flow is dynamic.

Gravity Flow

The simplified gravity flow procedure assumes a uniform vertical gradient and ignores capillary forces.

The simplified gravity flow procedure provides a suitable solution when you are primarily interested in the time varying recharge to the groundwater table based on actual precipitation and evapotranspiration and not the dynamics in the unsaturated zone.

Two Layer Water Balance

The simple two-layer water balance method divides the unsaturated zone into two zones: the root zone and the zone between the roots and the water table.

The simple two-layer water balance method is suitable when the water table is shallow, and groundwater recharge is primarily influenced by evapotranspiration in the root zone.

2. UZ-SZ coupling

Specific Yield of upper SZ layer

MIKE SHE forces the specific yield of the top SZ layer to be equal to the “specific yield” of the UZ zone as defined by the difference between the specified moisture contents at saturation, \(\theta_s\), and field capacity, \(\theta_{fc}\). This correction is calculated from the UZ values in the UZ cell in which the initial SZ water table is located.

Ensuring the correct UZ thickness

The following procedure could be used to ensure that the unsaturated zone does not drop below the bottom of the first calculation layer of the saturated zone:

  1. After a simulation, create a map of grid statistics of the potential head in the first calculation layer of the saturated zone
  2. Subtract the map of the minimum potential head from the map of the bottom level of the first calculation layer of the saturated zone.
  3. View the difference map. If the difference is very small in some areas of the map (e.g. > 0.5 m), it is strongly advised to move the bottom level of the first calculation layer of the saturated zone downwards.
  4. Repeat this procedure until there are no small differences.

Evaluation of the UZ-SZ Coupling

The WM_Print log file generated by MIKE SHE should be reviewed after each simulation to evaluate the performance of the UZ module. If the user specified maximum UZ iterations is exceeded an excessive number of times and there are no problems with the soil data used in the UZ module, the UZ and SZ time step should be evaluated. Sometimes it is possible to reduce the number of times the maximum UZ iterations is exceeded by making the UZ and SZ time steps more similar. Typically the SZ to UZ time step ratio should be no larger than four.

It is also useful to save the value of Ecum as a grid series output, or as a detailed time series output at critical locations. These plots can be used to determine if there are locations or periods of time during the simulation where the Ecum term exceeds Emax. This can occur if:

  • the water table drops below the first SZ calculation layer (positive value),
  • the water table rises above the top of the first SZ calculation layer (negative value),
  • the vertical hydraulic conductivity in the upper SZ calculation layer is much greater than the saturated hydraulic conductivity used in the UZ, or if:
  • the drainage time constant is too high.

In the first two cases above, the epsilon term can exceed Emax because the UZ module cannot get rid of epsilon because there is no available storage for the error term. In the third case, the UZ and SZ hydraulic properties should be consistent, or it will be difficult for MIKE SHE to simulate consistent vertical flow rates. In the last case, the drainage time constant should be reduced to prevent excessive and unrealistic drainage outflows from the SZ module.

Limitations of the UZ - SZ coupling

The coupling between UZ and SZ is limited to the top calculation layer of the saturated zone. This implies that:

  • As a rule of thumb, the UZ soil profiles should extend to just below the bottom of the top SZ layer.
  • However, if you have a very thick top SZ layer, then the UZ profiles must extend at least to below the deepest depth of the water table.
  • If the top layer of the SZ model dries out, then the UZ model usually assumes a lower pressure head boundary equal to the bottom of the uppermost SZ layer.
  • However, if the top layer of the SZ model dries out, and you are using the Richards Equation method, then you should ensure that there is one UZ node below the bottom of the top SZ layer. Otherwise, an error may be generated if there is an upwards capillary pressure gradient.
  • All outflow from the UZ column is always added to the top node of the SZ model.
  • UZ nodes below the water table and the bottom of the top SZ layer are ignored.

More detail on interaction between the lower UZ boundary and the SZ Layer 1 is given in the Section: Lower Boundary.