MohidWiki - User contributions [en]

Module PorousMedia

2012-01-17T18:32:31Z

Carina:

==Overview==
Module Porous Media is responsible for handling all water fluxes in soil. In particular the infiltration flow is obtained from the potential water column available for this process. The potential water column is given by the [[Module Basin]] (the module that handles interface between modules) as reported below:

<math>PIC=TF+IWC</math>

where:
:{|
| ''PIC'' || is the Potential Infiltration column (m)
|-
| ''TF'' || Throughfall is the percentage of the rain that reaches the soil (m)
|-
| ''IWC'' || is the Initial Water column (m)
|}

The infiltration flux is calculated by the Buckingham-Darcy equation using head gradient between surface runoff and first soil layer (Jury et al,1991).

==Main Processes==

===Porous Media Geometry===
Porous Media is a 3D domain delimited in its upper limit by topography and lower limit by soil bottom (defined by user).
In terms of soil definition it can be defined vertical horizons to link correspond to real soil horizons with different hydraulic carachteristics.

<htm><a href="http://content.screencast.com/users/jovem/folders/Jing/media/b9ef6c79-c1d1-46ea-94f7-623471c15883/MohidLandSoilProfile.png"><img src="http://content.screencast.com/users/jovem/folders/Jing/media/b9ef6c79-c1d1-46ea-94f7-623471c15883/MohidLandSoilProfile.png" width="650" height="436" border="0" /></a></htm>

===Water Flow===
Soil contains a large distribution of pore sizes and channels through which water may flow. In general the water flow determination is based on the mass conservation and momentum equation. In the case of soil it is assumed that the forces of inertia are almost zero; therefore the balance is reduced to the forces of pressure, gravity and viscous. The equation that describes the flow through soil is the Buckingham Darcy equation (Jury et al,1991).

<math>v=-K\left ( \theta \right )\left ( \frac{\partial H}{\partial z} \right )=-K\left ( \theta \right )\left ( \frac{\partial h}{\partial z} + 1 \right )\,\,\,\,(1.1)\\ </math>

where:
:{|
|''v'' || is the water velocity at the cell interface (m/s)
|-
| ''H'' || is the hydraulic head (m)
|-
| ''θ'' || is the water content (m3/m3)
|-
| ''K''|| is the hydraulic conductivity (m/s)
|}

The hydraulic head is given by the formula:

<math>H=h+p+z\,\,\,\,(1.2)</math>

where:
:{|
|-
| ''h'' || is the hydraulic head (m)
|-
| ''p'' || is hydrostatic pressure (m)
|-
| ''z''|| is the topography (m)
|-
|}

When in saturated conditions, hydraulic head is zero and hydrostatic pressure may occur (if water is at rest or decelerating). In unsaturated conditions hydrostatic pressure is zero and hydarulic head exists.

The soil is a very complex system, made up of a heterogeneous mixture of solid, liquid, and gaseous material. The liquid phase consists of soil water, which fills part or all of the open spaces between the soil particles. Therefore it is possible to divided the soil in two layers:

*Saturated soil <math>\Longrightarrow</math> The soil pores are filled by water

*Unsaturated one <math>\Longrightarrow</math> The soil pores are filled by water and air

In the first case the equation of Buckingham Darcy is simplified to the Darcy law and the parameter associated for its resolution are connected with the saturated layer. On the other hand for the resolution of the equation (1.1) a description of the characteristics of the the unsaturated layer is needed.

==Vadose Zone==

Many vadose flow and transport studies require description of unsaturated soil hydraulic proprieties over a wide range of pressure head. The hydraulic proprieties are described using the porous size distribution model of Maulem (1976) for hydraulic conductivity in combination with a water retention function introduced by Van Genuchten (1980).

===Water content===

Water content is the quantity of water contained in the soil (called '''soil moisture'''). It is given as a volumetric basis and it is defined mathematically as:

:<math>\theta = \frac{V_w}{V_T}\hspace{5cm}(1.3) </math>

:<math>V_T = V_s + V_v = V_s + V_w + V_a\hspace{1.6cm}(1.4)</math>

where:
:{|
| ''<math>\theta</math>'' || is water content (m 3/m 3)
|-
|''<math>V_w</math> '' || is the volume of water (m s3)
|-
| ''<math>V_T</math>'' || is the total volume (m 3)
|-
| ''<math>V_s</math>'' || is the soil volume (m 3)
|-
| ''<math>V_a</math>'' || is the air space (m 3)
|}

====Initial Condition====
The user may define a water content initialization or choose the model to compute in the unsaturated area heads so that soil water is at [[Field Capacity]].
The below explains the latter.

Once determined the aquifer level ('''water table''') the water content is associated at each cells by the following criteria:
:{|
*If the cell is located above the water table <math>\theta=\theta_{s}</math>
|-
*If the cell is located over the water table <math>\theta=\theta_{ns}</math>
|}

where:
:{|
| ''<math>\theta=\theta_{s}</math>'' || is the water content in the saturated soil (m3/m3)
|-
| ''<math>\theta=\theta_{ns}</math>'' || is the water content in the non saturated soil (m3/m3)
|-
|}

In order to calculate field capacity by the equation (1.7) the evaluation of the suction head is needed :

:<math> h=-DWZ\cdot 0.5\hspace{2cm}for\,\, the\,\, cells\,\, immediately\,\, above\,\, the\,\, water\,\, table \hspace{5cm} (1.5)</math>

:<math>h=-(-DZZ-h)\hspace{1.5cm}for\,\, the\,\, other\,\, cells\hspace{9.75cm} (1.6) </math>

[[Image:Figure05.jpg|thumb|center|300px|Figure 1: Suction Head Calculation]]

As shown in the picture the suction head is calculated in order to maintain the same total head (H = z + p + h) in the cells in agreement with the field capacity definition.

===Water retention===
The model use for characterizing the shape of water retention curves is the van Genuchten model:

:<math>\theta(h) = \theta_r + \frac{\theta_s - \theta_r}{\left[ 1+(\alpha |h|)^n \right]^{1-1/n}}\hspace{0.5cm}\Longrightarrow\hspace{0.5cm}h(\theta)=\left |\frac{(S_{E}^{-1/n}-1)^{1/n}}{\alpha} \right |\hspace{3cm}(1.7)</math>

:<math>S_{E}=\frac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}}\hspace{11cm}(1.8)</math>

where:
:{|
|''<math>\theta(h)</math>'' || is the the water retention curve (m 3/m 3)
|-
| ''<math>\theta_s</math>'' || is the saturated water content (m 3/m 3)
|-
| ''<math>\theta_r</math>'' || is the residual water content (m 3/m 3)
|-
| ''h''|| is the suction head (m)
|-
| ''<math>\alpha</math>'' || is related to the inverse of the air entry (m -1)
|-
| ''n''|| is a measure of the pore-size distribution n>1 (dimensionless)
|-
| ''<math>S_{E}</math>''|| is the effective saturation (dimensionless)
|}

===Saturated and Unsaturated Conductivity===

The saturated conductivity is given depending on the type of soil; instead the unsaturated conductivity is obtained from the suction head by the Maulem model:

<math>K(\theta)=K_{s}\cdot Se^{L}\cdot (1-(1-Se^{1/m})^{m})^{2}</math>

where:
:{|
|''<math>K(\theta)</math>'' || is the unsaturated conductivity (m/s)
|-
| ''<math>K_{s}</math>'' || is the saturated conductivity(m/s)
|-
| ''L'' || empirical pore-connectivity (m)
|-
| ''m''|| m=1-1/n(dimensionless)
|}

==Evapotranspiration==

Some water may be extracted from the soil because of the evaporation and transpiration processes, which become a sink in soil water profile. These two processes are currently named Evapotranspiration.

Potential Evapotranspiration may be modeled using the Penmann Monteith equation.
Also, if vegetation exists, a differentiation between Potential Transpiration and Potential Evaporation is done using LAI.
These computation are made in [[Module Basin]] since this is the module that handles water fluxes in the interface betwen modules.

However, not all of the potential water that can be evaporated or transpired will be in fact removed from the soil. The water that will really leave the soil through these processes is calculated: i) if vegetation exists, effective transpiration is computed in module Vegetation; ii) effective evaporation is computed in the Porous Media module.
In Figure below it can be seen that the actual transpiration and evaporation are then used in Porous Media module to compute the new water content.
The actual evaporation, which happens only at the soil surface, is calculated based on: i) a pressure head limit or ii) a soil conductivity limit, chosen by the user. It allows the model not to evaporate any surface water, even if it is available, when the soil head gets below the assigned value (i) or limits evaporation velocity to layer unsaturated conductivity (ii).

[[Image:evapotranspiration fluxogram.jpg|thumb|center|400px|Evapotranspiration fluxogram in Mohid Land model]]
Remind that Feddes is one option for computing effective transpiration in plants in [[module Vegetation]].

==Iteration Process==

Once obtained the water fluxes a balance on the water volume of each cell is apply in order to obtain the new water content <math>\theta</math>. The balance applied is the following:

'''Horizontal Direction X'''

<math>\theta=(\theta\cdot V_{cell}+(FluxU_{(i,j,k)}\cdot ComputeFace_{(i,j,k)}-Flux_{(i,j+1,k)}\cdot ComputeFace_{(i,j+1,k)}\cdot \Delta t) V_{cell}</math>

'''Horizontal Direction Y'''

<math>\theta=(\theta\cdot V_{cell}+(FluxV_{(i,j,k)}\cdot ComputeFace_{(i,j,k)}-FluxV_{(i+1,j,k)}\cdot ComputeFace_{(i+1,j,k)})\cdot \Delta t)\cdot V_{cell}</math>

'''Vertical Direction W'''

<math>\theta=(\theta\cdot V_{cell}+(FluxW_{(i,j,k)}\cdot ComputeFace_{(i,j,k})-FluxW_{(i,j,k+1)}\cdot ComputeFace_{(i,j,k+1}))\cdot \Delta t)\cdot V_{cell}</math>

'''Transpiration Flux'''

<math>\theta=(\theta\cdot V_{cell}-(TranspFlux_{(i,j,k)}\Delta t)\cdot V_{cell}</math>

'''Evaporation Flux'''

<math>\theta=(\theta\cdot V_{cell}-(EvapFlux_{(i,j,k)}\Delta t)\cdot V_{cell}</math>

'''Infiltration Flux'''

<math>\theta=(\theta\cdot V_{cell}-(UnsatK\cdot Area_{cell}\cdot(1-Imp) \Delta t)\cdot V_{cell}</math>

where:

:{|

|''<math>\theta</math>'' || is the water content (m 3/m 3)
|-
|''<math>V_{cell}</math>'' || is the cell volume (m 3)
|-
|''<math>Area_{cell}</math>'' || is the cell area (m 2)
|-
|''ComputeFace'' || is the computed face (dimensionless)
|-
|''<math>\Delta t</math>'' || is the time step (s)
|-
|''FluxX'' || is the flux in X direction (m 3/s)
|-
|''FluxV'' || is the flux in Y direction (m 3/s)
|-
|''FluxW'' || is the flux in W direction (m 3/s)
|-
|''TranspFlux'' || is the transpiration flux (m 3/s)
|-
|''EvapFlux'' || is the evaporation flux (m 3/s)
|-
|''UnsatK'' || is unsaturated conductivity (m /s)
|-
|''Imp'' || is the percentage of impermeable soil of the cell (dimensionless)
|}

The value of <math>\theta</math> so obtained is compared with one used in the volumes calculation and the iterative process stop when:

<math>(\theta^{'})-(\theta^{new})<\,\, Tolerance\hspace{3cm} (1.9) </math>

where:

:{|
|''<math>\theta^{'}</math>'' || is the water content of the previous iteration (m 3/m 3)
|}

If the equation (1.9) is not satisfy the temporal step is divided in half and the new value of <math>\theta</math> is used for solving the equation (1.7) for restarting the calculation process. The iteration process is stoped when the tolerance desired is reached.

[[Image:Figure06.jpg|thumb|center|300px|Figure 2: Time step reduction]]

==Other Features==
===How to Generate Info needed in Porous Media===
====SoilMap====
Model needs to know soil ID in each cell and layer to pick hydraulic properties from that type of soil. In Pedology soil includes more than one horizon, each with different soil properties. Here Soil is used has a unit of soil hidraulic properties, i.e., to define a soil with three horizons one has to create three SoilID (see below). This also means that if a watershed has at least one soil with three horizons one has to create three soil maps. Each soil map will be infact the map of each horizon of the soils. The grid cells with only one horizon will have the same SoilID in all maps, the grid cells with three horizons will have a different SoilID in each map.

Options:
* Constant value
* Soil Grid. One possible option is to associate with soil shape file. In this case can use MOHID GIS going to menu [Tools]->[Shape to Grid Data] and provide: i) the grid (model grid), ii) the soil shape file and iii) the corespondence between soil codes and soil ID defined in data file.

Soil ID must be defined in [[Module_FillMatrix|Module FillMatrix]] standards for each soil horizon defined (grid example):
<beginhorizon>
KLB : 1
KUB : 10

<beginproperty>
NAME : SoilID
DEFAULTVALUE : 1
INITIALIZATION_METHOD : ASCII_FILE
REMAIN_CONSTANT : 1
FILENAME : ..\..\GeneralData\PorousMedia\SoilID200m.dat
<endproperty>

<beginproperty>
<endproperty>
..
<endhorizon>

'''Remarks'''

All the soil ID's appearing in the soil grid(s) must be defined in the PorousMedia data file in terms of hydraulic properties:
<beginsoiltype>
ID : 1
THETA_S : 0.3859
THETA_R : 0.0476
N_FIT : 1.39
SAT_K : 3.5556e-6
ALPHA : 2.75
L_FIT : 0.50
<endsoiltype>

<beginsoiltype>
ID : 2
...
<endsoiltype>
...

====Soil Bottom====
The soil depth must be known by the model. This is computed by the model from terrain altitude (topography) and soil bottom altitude. As so, a soil bottom grid is needed.

Options:
* Grid File.
Soil depth (and soil bottom altitude, the effective grid needed) can be defined with a constant depth or estimated from slope [[HOW TO SoilBottom LINK]]. When the soil depth is estimated as a function of slope, soil depth will be smaller in ares with higher slope. In this areas only the surface layers of the soil will be considered (see [[Module_PorousMedia#Porous_Media_Geometry|Porous Media Geometry]]).

Define the grid just generated, in the porous media data file with:
BOTTOM_FILE : ..\..\GeneralData\PorousMedia\BottomLevel.dat

====Water Level====
Options:
*Grid File.
The water table altitude represents the initial altitude of the water table.
It is recommended to do a spin-up run to estabilize water level and then do a continuous simulation starting with the final water table achieved.
Use the following blocks with [[Module_FillMatrix|Module FillMatrix]] standards:
<beginwaterlevel>
NAME : waterlevel
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 0
FILENAME : ..\..\GeneralData\PorousMedia\WaterLevel0.50.dat
<endwaterlevel>

====Impermeability====
Impermeability values (0 - completely permeable, 1 - impermeable) must be provided.

Options:
* Constant Value.
* Grid File. One possible option is to associate with land use shape file. In this case can use MOHID GIS going to menu [Tools]->[Shape to Grid Data] and provide: i) the grid (model grid), ii) the land use shape file and iii) the corespondence between land use codes and Impermeability values.
Use the following blocks with [[Module_FillMatrix|Module FillMatrix]] standards:
<beginimpermeablefraction>
NAME : impermeablefraction
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 1
FILENAME : ..\..\GeneralData\PorousMedia\AreaImpermeavel.dat
<endimpermeablefraction>

==Outputs==

===Timeseries===

Theta - is water content of selected cell (vol water/ vol soil)

relative_water_content - content in selected cell. between zero and one (zero is residual water content and one is saturated water content)

VelW_[m/s] - vertical velocity in the bottom face of the selected cell

VelW_Corr_[m/s] - vertical velocity in the bottom face of the selected cell that may be corrected if oversaturation occurs. if no correction occurs is the same as previous.

InF_Vel_[m/s] - infiltration velocity (in the soil surface)

Head_[m] - Suction in selected cell

Conductivity_[m/s] - Conductivity in selected cell

level_water_table_[m] - water table altitude

water_table_depth_[m] - water table depth (from soil surface)

Hydro_Pressure_[m] - hydrostatic pressure in selected cell

Final_Head_[m] - Soil water charge in selected cell

[Check Mohid Land Heights and Levels to understand some of the outputs]

GW_flow_to_river_total_[m3/s] - Ground water flow to river if it is a river point

Surface_Evaporation_Flux_[m3/s] - evaporation flux (in the soil surface)

Transpiration_Flux_[m3/s] - transpiration flux in selected cell

[http://screencast.com/t/I8We0B4j2uJ Check Mohid Land Heights and Levels to understand the difference between level and height results]

==References ==
*Jury,W.A.,Gardner,W.R.,Gardner,W.H., 1991,Soil Physics
*Van Genuchten, M.T., A closed form equation for predicting the hydraulic conductivity of unsaturated soils
*Wu,J.,Zhang, R., Gui,S.,1999, Modelling soil water movement with water uptake by roots, Plant and soil 215: 7-17
*Marcel G.Schaap and Martinus Th. van Genuchten, A modified Maulem van Genuchten Formulation for Improved Description of Hydraulic Conductivity Near Saturation, 16 December 2005

==Data File ==

===Keywords===
Keywords read in the Data File

Keyword : Data Type Default !Comment

BOTTOM_FILE : char - !Path to Bottom Topography File
START_WITH_FIELD : logical 1 !Sets Theta initial Field Capacity
CONTINUOUS : logical 0 !Continues from previous run
STOP_ON_WRONG_DATE : logical 1 !Stops if previous run end is different from actual
!Start
OUTPUT_TIME : sec. sec. sec. - !Output Time
TIME_SERIE_LOCATION : char - !Path to File which defines Time Series
CONTINUOUS_OUTPUT_FILE : logical 1 !Writes "famous" iter.log
CONDUTIVITYFACE : integer 1 !Way to interpolate conducivity face
!1 - Average, 2 - Maximum, 3 - Minimum, 4 - Weigthed
HORIZONTAL_K_FACTOR : real 1.0 !Factor for Horizontal Conductivity = Kh / Kv
CUT_OFF_THETA_LOW : real 1e-6 !Disables calculation when Theta is near ThetaR
CUT_OFF_THETA_HIGH : real 1e-15 !Set Theta = ThetaS when Theta > ThetaS - CUT_OFF_THETA_HIGH
MIN_ITER : integer 2 !Number of iterations below which the DT is increased
MAX_ITER : integer 3 !Number of iterations above which the DT is decreased
LIMIT_ITER : integer 50 !Number of iterations of a time step (for restart)
THETA_TOLERANCE : real 0.001 !Converge Parameter
INCREASE_DT : real 1.25 !Increase of DT when iter < MIN_ITER
DECREASE_DT : real 0.70 !Decrease of DT when iter > MAX_ITER

<beginproperty>
NAME : Theta / waterlevel

see Module FillMatrix for more options

<endproperty>

===Sample===

Some keywords of the PorousMedia input file:
BOTTOM_FILE : ..\General Data\Other\PorousMedia\SoilBottom_2cells.dat
TIME_SERIE_LOCATION : ..\General Data\TimeSeries\TimeSeriesLocation3D_2m.dat

COMPUTE_SOIL_FIELD : 1

OUTPUT_TIME : 0 86400
CUT_OFF_THETA_HIGH : 1e-15
START_WITH_FIELD : 1

LIMIT_EVAP_WATER_VEL : 0
LIMIT_EVAP_HEAD : 0
HEAD_LIMIT : -100
THETA_HYDRO_COEF : 0.98

<beginsoiltype>
ID : 1
THETA_S : 0.43
THETA_R : 0.078
SAT_K : 2.888e-6
N_FIT : 1.56
ALPHA : 3.6
L_FIT : 0.50
<endsoiltype>

!----- Hydraulic Soil Properties
<beginhorizon>
KLB : 1
KUB : 10

<beginproperty>
NAME : SoilID
INITIALIZATION_METHOD : CONSTANT
DEFAULTVALUE : 1
<endproperty>

<beginproperty>
NAME : Theta
INITIALIZATION_METHOD : CONSTANT
DEFAULTVALUE : 0.30
<endproperty>
<endhorizon>

<beginwaterlevel>
NAME : waterlevel
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 1.
REMAIN_CONSTANT : 1
FILENAME : ..\General Data\Initial Conditions\InitialWaterLevel_2cells.dat
<endwaterlevel>

<beginimpermeablefraction>
NAME : impermeablefraction
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 1
FILENAME : ..\General Data\Other\PorousMedia\InitialImpermeabilization_2cells.dat
<endimpermeablefraction>

[[Category:MOHID Land]]

Module PorousMedia

2012-01-17T18:29:15Z

Carina:

==Overview==
Module Porous Media is responsible for handling all water fluxes in soil. In particular the infiltration flow is obtained from the potential water column available for this process. The potential water column is given by the [[Module Basin]] (the module that handles interface between modules) as reported below:

<math>PIC=TF+IWC</math>

where:
:{|
| ''PIC'' || is the Potential Infiltration column (m)
|-
| ''TF'' || Throughfall is the percentage of the rain that reaches the soil (m)
|-
| ''IWC'' || is the Initial Water column (m)
|}

The infiltration flux is calculated by the Buckingham-Darcy equation using head gradient between surface runoff and first soil layer (Jury et al,1991).

==Main Processes==

===Porous Media Geometry===
Porous Media is a 3D domain delimited in its upper limit by topography and lower limit by soil bottom (defined by user).
In terms of soil definition it can be defined vertical horizons to link correspond to real soil horizons with different hydraulic carachteristics.

<htm><a href="http://content.screencast.com/users/jovem/folders/Jing/media/b9ef6c79-c1d1-46ea-94f7-623471c15883/MohidLandSoilProfile.png"><img src="http://content.screencast.com/users/jovem/folders/Jing/media/b9ef6c79-c1d1-46ea-94f7-623471c15883/MohidLandSoilProfile.png" width="650" height="436" border="0" /></a></htm>

===Water Flow===
Soil contains a large distribution of pore sizes and channels through which water may flow. In general the water flow determination is based on the mass conservation and momentum equation. In the case of soil it is assumed that the forces of inertia are almost zero; therefore the balance is reduced to the forces of pressure, gravity and viscous. The equation that describes the flow through soil is the Buckingham Darcy equation (Jury et al,1991).

<math>v=-K\left ( \theta \right )\left ( \frac{\partial H}{\partial z} \right )=-K\left ( \theta \right )\left ( \frac{\partial h}{\partial z} + 1 \right )\,\,\,\,(1.1)\\ </math>

where:
:{|
|''v'' || is the water velocity at the cell interface (m/s)
|-
| ''H'' || is the hydraulic head (m)
|-
| ''θ'' || is the water content (m3/m3)
|-
| ''K''|| is the hydraulic conductivity (m/s)
|}

The hydraulic head is given by the formula:

<math>H=h+p+z\,\,\,\,(1.2)</math>

where:
:{|
|-
| ''h'' || is the hydraulic head (m)
|-
| ''p'' || is hydrostatic pressure (m)
|-
| ''z''|| is the topography (m)
|-
|}

When in saturated conditions, hydraulic head is zero and hydrostatic pressure may occur (if water is at rest or decelerating). In unsaturated conditions hydrostatic pressure is zero and hydarulic head exists.

The soil is a very complex system, made up of a heterogeneous mixture of solid, liquid, and gaseous material. The liquid phase consists of soil water, which fills part or all of the open spaces between the soil particles. Therefore it is possible to divided the soil in two layers:

*Saturated soil <math>\Longrightarrow</math> The soil pores are filled by water

*Unsaturated one <math>\Longrightarrow</math> The soil pores are filled by water and air

In the first case the equation of Buckingham Darcy is simplified to the Darcy law and the parameter associated for its resolution are connected with the saturated layer. On the other hand for the resolution of the equation (1.1) a description of the characteristics of the the unsaturated layer is needed.

==Vadose Zone==

Many vadose flow and transport studies require description of unsaturated soil hydraulic proprieties over a wide range of pressure head. The hydraulic proprieties are described using the porous size distribution model of Maulem (1976) for hydraulic conductivity in combination with a water retention function introduced by Van Genuchten (1980).

===Water content===

Water content is the quantity of water contained in the soil (called '''soil moisture'''). It is given as a volumetric basis and it is defined mathematically as:

:<math>\theta = \frac{V_w}{V_T}\hspace{5cm}(1.3) </math>

:<math>V_T = V_s + V_v = V_s + V_w + V_a\hspace{1.6cm}(1.4)</math>

where:
:{|
| ''<math>\theta</math>'' || is water content (m 3/m 3)
|-
|''<math>V_w</math> '' || is the volume of water (m s3)
|-
| ''<math>V_T</math>'' || is the total volume (m 3)
|-
| ''<math>V_s</math>'' || is the soil volume (m 3)
|-
| ''<math>V_a</math>'' || is the air space (m 3)
|}

====Initial Condition====
The user may define a water content initialization or choose the model to compute in the unsaturated area heads so that soil water is at [[Field Capacity]].
The below explains the latter.

Once determined the aquifer level ('''water table''') the water content is associated at each cells by the following criteria:
:{|
*If the cell is located above the water table <math>\theta=\theta_{s}</math>
|-
*If the cell is located over the water table <math>\theta=\theta_{ns}</math>
|}

where:
:{|
| ''<math>\theta=\theta_{s}</math>'' || is the water content in the saturated soil (m3/m3)
|-
| ''<math>\theta=\theta_{ns}</math>'' || is the water content in the non saturated soil (m3/m3)
|-
|}

In order to calculate field capacity by the equation (1.7) the evaluation of the suction head is needed :

:<math> h=-DWZ\cdot 0.5\hspace{2cm}for\,\, the\,\, cells\,\, immediately\,\, above\,\, the\,\, water\,\, table \hspace{5cm} (1.5)</math>

:<math>h=-(-DZZ-h)\hspace{1.5cm}for\,\, the\,\, other\,\, cells\hspace{9.75cm} (1.6) </math>

[[Image:Figure05.jpg|thumb|center|300px|Figure 1: Suction Head Calculation]]

As shown in the picture the suction head is calculated in order to maintain the same total head (H = z + p + h) in the cells in agreement with the field capacity definition.

===Water retention===
The model use for characterizing the shape of water retention curves is the van Genuchten model:

:<math>\theta(h) = \theta_r + \frac{\theta_s - \theta_r}{\left[ 1+(\alpha |h|)^n \right]^{1-1/n}}\hspace{0.5cm}\Longrightarrow\hspace{0.5cm}h(\theta)=\left |\frac{(S_{E}^{-1/n}-1)^{1/n}}{\alpha} \right |\hspace{3cm}(1.7)</math>

:<math>S_{E}=\frac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}}\hspace{11cm}(1.8)</math>

where:
:{|
|''<math>\theta(h)</math>'' || is the the water retention curve (m 3/m 3)
|-
| ''<math>\theta_s</math>'' || is the saturated water content (m 3/m 3)
|-
| ''<math>\theta_r</math>'' || is the residual water content (m 3/m 3)
|-
| ''h''|| is the suction head (m)
|-
| ''<math>\alpha</math>'' || is related to the inverse of the air entry (m -1)
|-
| ''n''|| is a measure of the pore-size distribution n>1 (dimensionless)
|-
| ''<math>S_{E}</math>''|| is the effective saturation (dimensionless)
|}

===Saturated and Unsaturated Conductivity===

The saturated conductivity is given depending on the type of soil; instead the unsaturated conductivity is obtained from the suction head by the Maulem model:

<math>K(\theta)=K_{s}\cdot Se^{L}\cdot (1-(1-Se^{1/m})^{m})^{2}</math>

where:
:{|
|''<math>K(\theta)</math>'' || is the unsaturated conductivity (m/s)
|-
| ''<math>K_{s}</math>'' || is the saturated conductivity(m/s)
|-
| ''L'' || empirical pore-connectivity (m)
|-
| ''m''|| m=1-1/n(dimensionless)
|}

==Evapotranspiration==

Some water may be extracted from the soil because of the evaporation and transpiration processes, which become a sink in soil water profile. These two processes are currently named Evapotranspiration.

Potential Evapotranspiration may be modeled using the Penmann Monteith equation.
Also, if vegetation exists, a differentiation between Potential Transpiration and Potential Evaporation is done using LAI.
These computation are made in [[Module Basin]] since this is the module that handles water fluxes in the interface betwen modules.

However, not all of the potential water that can be evaporated or transpired will be in fact removed from the soil. The water that will really leave the soil through these processes is calculated: i) if vegetation exists, effective transpiration is computed in module Vegetation; ii) effective evaporation is computed in the Porous Media module.
In Figure below it can be seen that the actual transpiration and evaporation are then used in Porous Media module to compute the new water content.
The actual evaporation, which happens only at the soil surface, is calculated based on: i) a pressure head limit or ii) a soil conductivity limit, chosen by the user. It allows the model not to evaporate any surface water, even if it is available, when the soil head gets below the assigned value (i) or limits evaporation velocity to layer unsaturated conductivity (ii).

[[Image:evapotranspiration fluxogram.jpg|thumb|center|400px|Evapotranspiration fluxogram in Mohid Land model]]
Remind that Feddes is one option for computing effective transpiration in plants in [[module Vegetation]].

==Iteration Process==

Once obtained the water fluxes a balance on the water volume of each cell is apply in order to obtain the new water content <math>\theta</math>. The balance applied is the following:

'''Horizontal Direction X'''

<math>\theta=(\theta\cdot V_{cell}+(FluxU_{(i,j,k)}\cdot ComputeFace_{(i,j,k)}-Flux_{(i,j+1,k)}\cdot ComputeFace_{(i,j+1,k)}\cdot \Delta t) V_{cell}</math>

'''Horizontal Direction Y'''

<math>\theta=(\theta\cdot V_{cell}+(FluxV_{(i,j,k)}\cdot ComputeFace_{(i,j,k)}-FluxV_{(i+1,j,k)}\cdot ComputeFace_{(i+1,j,k)})\cdot \Delta t)\cdot V_{cell}</math>

'''Vertical Direction W'''

<math>\theta=(\theta\cdot V_{cell}+(FluxW_{(i,j,k)}\cdot ComputeFace_{(i,j,k})-FluxW_{(i,j,k+1)}\cdot ComputeFace_{(i,j,k+1}))\cdot \Delta t)\cdot V_{cell}</math>

'''Transpiration Flux'''

<math>\theta=(\theta\cdot V_{cell}-(TranspFlux_{(i,j,k)}\Delta t)\cdot V_{cell}</math>

'''Evaporation Flux'''

<math>\theta=(\theta\cdot V_{cell}-(EvapFlux_{(i,j,k)}\Delta t)\cdot V_{cell}</math>

'''Infiltration Flux'''

<math>\theta=(\theta\cdot V_{cell}-(UnsatK\cdot Area_{cell}\cdot(1-Imp) \Delta t)\cdot V_{cell}</math>

where:

:{|

|''<math>\theta</math>'' || is the water content (m 3/m 3)
|-
|''<math>V_{cell}</math>'' || is the cell volume (m 3)
|-
|''<math>Area_{cell}</math>'' || is the cell area (m 2)
|-
|''ComputeFace'' || is the computed face (dimensionless)
|-
|''<math>\Delta t</math>'' || is the time step (s)
|-
|''FluxX'' || is the flux in X direction (m 3/s)
|-
|''FluxV'' || is the flux in Y direction (m 3/s)
|-
|''FluxW'' || is the flux in W direction (m 3/s)
|-
|''TranspFlux'' || is the transpiration flux (m 3/s)
|-
|''EvapFlux'' || is the evaporation flux (m 3/s)
|-
|''UnsatK'' || is unsaturated conductivity (m /s)
|-
|''Imp'' || is the percentage of impermeable soil of the cell (dimensionless)
|}

The value of <math>\theta</math> so obtained is compared with one used in the volumes calculation and the iterative process stop when:

<math>(\theta^{'})-(\theta^{new})<\,\, Tolerance\hspace{3cm} (1.9) </math>

where:

:{|
|''<math>\theta^{'}</math>'' || is the water content of the previous iteration (m 3/m 3)
|}

If the equation (1.9) is not satisfy the temporal step is divided in half and the new value of <math>\theta</math> is used for solving the equation (1.7) for restarting the calculation process. The iteration process is stoped when the tolerance desired is reached.

[[Image:Figure06.jpg|thumb|center|300px|Figure 2: Time step reduction]]

==Other Features==
===How to Generate Info needed in Porous Media===
====SoilMap====
Model needs to know soil ID in each cell and layer to pick hydraulic properties from that type of soil. In Pedology soil includes more than one horizon, each with different soil properties. Here Soil is used has a unit of soil hidraulic properties, i.e., to define a soil with three horizons one has to create three SoilID (see below). This also means that if a watershed has at least one soil with three horizons one has to create three soil maps. Each soil map will be infact the map of each horizon of the soils. The grid cells with only one horizon will have the same SoilID in all maps, the grid cells with three horizons will have a different SoilID in each map.

Options:
* Constant value
* Soil Grid. One possible option is to associate with soil shape file. In this case can use MOHID GIS going to menu [Tools]->[Shape to Grid Data] and provide: i) the grid (model grid), ii) the soil shape file and iii) the corespondence between soil codes and soil ID defined in data file.

Soil ID must be defined in [[Module_FillMatrix|Module FillMatrix]] standards for each soil horizon defined (grid example):
<beginhorizon>
KLB : 1
KUB : 10

<beginproperty>
NAME : SoilID
DEFAULTVALUE : 1
INITIALIZATION_METHOD : ASCII_FILE
REMAIN_CONSTANT : 1
FILENAME : ..\..\GeneralData\PorousMedia\SoilID200m.dat
<endproperty>

<beginproperty>
<endproperty>
..
<endhorizon>

'''Remarks'''

All the soil ID's appearing in the soil grid(s) must be defined in the PorousMedia data file in terms of hydraulic properties:
<beginsoiltype>
ID : 1
THETA_S : 0.3859
THETA_R : 0.0476
N_FIT : 1.39
SAT_K : 3.5556e-6
ALPHA : 2.75
L_FIT : 0.50
<endsoiltype>

<beginsoiltype>
ID : 2
...
<endsoiltype>
...

====Soil Bottom====
The soil depth must be known by the model. This is computed by the model from terrain altitude (topography) and soil bottom altitude. As so, a soil bottom grid is needed.

Options:
* Grid File.
Soil depth (and soil bottom altitude, the effective grid needed) can be defined with a constant depth or estimated from slope [[HOW TO SoilBottom LINK]]. When the soil depth is estimated as a function of slope, soil depth will be smaller in ares with higher slope. In this areas only the surface layers of the soil will be considered (see [[Module_PorousMedia#Porous_Media_Geometry]]).

Define the grid just generated, in the porous media data file with:
BOTTOM_FILE : ..\..\GeneralData\PorousMedia\BottomLevel.dat

====Water Level====
Options:
*Grid File.
The water table altitude represents the initial altitude of the water table.
It is recommended to do a spin-up run to estabilize water level and then do a continuous simulation starting with the final water table achieved.
Use the following blocks with [[Module_FillMatrix|Module FillMatrix]] standards:
<beginwaterlevel>
NAME : waterlevel
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 0
FILENAME : ..\..\GeneralData\PorousMedia\WaterLevel0.50.dat
<endwaterlevel>

====Impermeability====
Impermeability values (0 - completely permeable, 1 - impermeable) must be provided.

Options:
* Constant Value.
* Grid File. One possible option is to associate with land use shape file. In this case can use MOHID GIS going to menu [Tools]->[Shape to Grid Data] and provide: i) the grid (model grid), ii) the land use shape file and iii) the corespondence between land use codes and Impermeability values.
Use the following blocks with [[Module_FillMatrix|Module FillMatrix]] standards:
<beginimpermeablefraction>
NAME : impermeablefraction
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 1
FILENAME : ..\..\GeneralData\PorousMedia\AreaImpermeavel.dat
<endimpermeablefraction>

==Outputs==

===Timeseries===

Theta - is water content of selected cell (vol water/ vol soil)

relative_water_content - content in selected cell. between zero and one (zero is residual water content and one is saturated water content)

VelW_[m/s] - vertical velocity in the bottom face of the selected cell

VelW_Corr_[m/s] - vertical velocity in the bottom face of the selected cell that may be corrected if oversaturation occurs. if no correction occurs is the same as previous.

InF_Vel_[m/s] - infiltration velocity (in the soil surface)

Head_[m] - Suction in selected cell

Conductivity_[m/s] - Conductivity in selected cell

level_water_table_[m] - water table altitude

water_table_depth_[m] - water table depth (from soil surface)

Hydro_Pressure_[m] - hydrostatic pressure in selected cell

Final_Head_[m] - Soil water charge in selected cell

[Check Mohid Land Heights and Levels to understand some of the outputs]

GW_flow_to_river_total_[m3/s] - Ground water flow to river if it is a river point

Surface_Evaporation_Flux_[m3/s] - evaporation flux (in the soil surface)

Transpiration_Flux_[m3/s] - transpiration flux in selected cell

[http://screencast.com/t/I8We0B4j2uJ Check Mohid Land Heights and Levels to understand the difference between level and height results]

==References ==
*Jury,W.A.,Gardner,W.R.,Gardner,W.H., 1991,Soil Physics
*Van Genuchten, M.T., A closed form equation for predicting the hydraulic conductivity of unsaturated soils
*Wu,J.,Zhang, R., Gui,S.,1999, Modelling soil water movement with water uptake by roots, Plant and soil 215: 7-17
*Marcel G.Schaap and Martinus Th. van Genuchten, A modified Maulem van Genuchten Formulation for Improved Description of Hydraulic Conductivity Near Saturation, 16 December 2005

==Data File ==

===Keywords===
Keywords read in the Data File

Keyword : Data Type Default !Comment

BOTTOM_FILE : char - !Path to Bottom Topography File
START_WITH_FIELD : logical 1 !Sets Theta initial Field Capacity
CONTINUOUS : logical 0 !Continues from previous run
STOP_ON_WRONG_DATE : logical 1 !Stops if previous run end is different from actual
!Start
OUTPUT_TIME : sec. sec. sec. - !Output Time
TIME_SERIE_LOCATION : char - !Path to File which defines Time Series
CONTINUOUS_OUTPUT_FILE : logical 1 !Writes "famous" iter.log
CONDUTIVITYFACE : integer 1 !Way to interpolate conducivity face
!1 - Average, 2 - Maximum, 3 - Minimum, 4 - Weigthed
HORIZONTAL_K_FACTOR : real 1.0 !Factor for Horizontal Conductivity = Kh / Kv
CUT_OFF_THETA_LOW : real 1e-6 !Disables calculation when Theta is near ThetaR
CUT_OFF_THETA_HIGH : real 1e-15 !Set Theta = ThetaS when Theta > ThetaS - CUT_OFF_THETA_HIGH
MIN_ITER : integer 2 !Number of iterations below which the DT is increased
MAX_ITER : integer 3 !Number of iterations above which the DT is decreased
LIMIT_ITER : integer 50 !Number of iterations of a time step (for restart)
THETA_TOLERANCE : real 0.001 !Converge Parameter
INCREASE_DT : real 1.25 !Increase of DT when iter < MIN_ITER
DECREASE_DT : real 0.70 !Decrease of DT when iter > MAX_ITER

<beginproperty>
NAME : Theta / waterlevel

see Module FillMatrix for more options

<endproperty>

===Sample===

Some keywords of the PorousMedia input file:
BOTTOM_FILE : ..\General Data\Other\PorousMedia\SoilBottom_2cells.dat
TIME_SERIE_LOCATION : ..\General Data\TimeSeries\TimeSeriesLocation3D_2m.dat

COMPUTE_SOIL_FIELD : 1

OUTPUT_TIME : 0 86400
CUT_OFF_THETA_HIGH : 1e-15
START_WITH_FIELD : 1

LIMIT_EVAP_WATER_VEL : 0
LIMIT_EVAP_HEAD : 0
HEAD_LIMIT : -100
THETA_HYDRO_COEF : 0.98

<beginsoiltype>
ID : 1
THETA_S : 0.43
THETA_R : 0.078
SAT_K : 2.888e-6
N_FIT : 1.56
ALPHA : 3.6
L_FIT : 0.50
<endsoiltype>

!----- Hydraulic Soil Properties
<beginhorizon>
KLB : 1
KUB : 10

<beginproperty>
NAME : SoilID
INITIALIZATION_METHOD : CONSTANT
DEFAULTVALUE : 1
<endproperty>

<beginproperty>
NAME : Theta
INITIALIZATION_METHOD : CONSTANT
DEFAULTVALUE : 0.30
<endproperty>
<endhorizon>

<beginwaterlevel>
NAME : waterlevel
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 1.
REMAIN_CONSTANT : 1
FILENAME : ..\General Data\Initial Conditions\InitialWaterLevel_2cells.dat
<endwaterlevel>

<beginimpermeablefraction>
NAME : impermeablefraction
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 1
FILENAME : ..\General Data\Other\PorousMedia\InitialImpermeabilization_2cells.dat
<endimpermeablefraction>

[[Category:MOHID Land]]

Module PorousMedia

2012-01-17T18:27:07Z

Carina:

==Overview==
Module Porous Media is responsible for handling all water fluxes in soil. In particular the infiltration flow is obtained from the potential water column available for this process. The potential water column is given by the [[Module Basin]] (the module that handles interface between modules) as reported below:

<math>PIC=TF+IWC</math>

where:
:{|
| ''PIC'' || is the Potential Infiltration column (m)
|-
| ''TF'' || Throughfall is the percentage of the rain that reaches the soil (m)
|-
| ''IWC'' || is the Initial Water column (m)
|}

The infiltration flux is calculated by the Buckingham-Darcy equation using head gradient between surface runoff and first soil layer (Jury et al,1991).

==Main Processes==

===Porous Media Geometry===
Porous Media is a 3D domain delimited in its upper limit by topography and lower limit by soil bottom (defined by user).
In terms of soil definition it can be defined vertical horizons to link correspond to real soil horizons with different hydraulic carachteristics.

<htm><a href="http://content.screencast.com/users/jovem/folders/Jing/media/b9ef6c79-c1d1-46ea-94f7-623471c15883/MohidLandSoilProfile.png"><img src="http://content.screencast.com/users/jovem/folders/Jing/media/b9ef6c79-c1d1-46ea-94f7-623471c15883/MohidLandSoilProfile.png" width="650" height="436" border="0" /></a></htm>

===Water Flow===
Soil contains a large distribution of pore sizes and channels through which water may flow. In general the water flow determination is based on the mass conservation and momentum equation. In the case of soil it is assumed that the forces of inertia are almost zero; therefore the balance is reduced to the forces of pressure, gravity and viscous. The equation that describes the flow through soil is the Buckingham Darcy equation (Jury et al,1991).

<math>v=-K\left ( \theta \right )\left ( \frac{\partial H}{\partial z} \right )=-K\left ( \theta \right )\left ( \frac{\partial h}{\partial z} + 1 \right )\,\,\,\,(1.1)\\ </math>

where:
:{|
|''v'' || is the water velocity at the cell interface (m/s)
|-
| ''H'' || is the hydraulic head (m)
|-
| ''θ'' || is the water content (m3/m3)
|-
| ''K''|| is the hydraulic conductivity (m/s)
|}

The hydraulic head is given by the formula:

<math>H=h+p+z\,\,\,\,(1.2)</math>

where:
:{|
|-
| ''h'' || is the hydraulic head (m)
|-
| ''p'' || is hydrostatic pressure (m)
|-
| ''z''|| is the topography (m)
|-
|}

When in saturated conditions, hydraulic head is zero and hydrostatic pressure may occur (if water is at rest or decelerating). In unsaturated conditions hydrostatic pressure is zero and hydarulic head exists.

The soil is a very complex system, made up of a heterogeneous mixture of solid, liquid, and gaseous material. The liquid phase consists of soil water, which fills part or all of the open spaces between the soil particles. Therefore it is possible to divided the soil in two layers:

*Saturated soil <math>\Longrightarrow</math> The soil pores are filled by water

*Unsaturated one <math>\Longrightarrow</math> The soil pores are filled by water and air

In the first case the equation of Buckingham Darcy is simplified to the Darcy law and the parameter associated for its resolution are connected with the saturated layer. On the other hand for the resolution of the equation (1.1) a description of the characteristics of the the unsaturated layer is needed.

==Vadose Zone==

Many vadose flow and transport studies require description of unsaturated soil hydraulic proprieties over a wide range of pressure head. The hydraulic proprieties are described using the porous size distribution model of Maulem (1976) for hydraulic conductivity in combination with a water retention function introduced by Van Genuchten (1980).

===Water content===

Water content is the quantity of water contained in the soil (called '''soil moisture'''). It is given as a volumetric basis and it is defined mathematically as:

:<math>\theta = \frac{V_w}{V_T}\hspace{5cm}(1.3) </math>

:<math>V_T = V_s + V_v = V_s + V_w + V_a\hspace{1.6cm}(1.4)</math>

where:
:{|
| ''<math>\theta</math>'' || is water content (m 3/m 3)
|-
|''<math>V_w</math> '' || is the volume of water (m s3)
|-
| ''<math>V_T</math>'' || is the total volume (m 3)
|-
| ''<math>V_s</math>'' || is the soil volume (m 3)
|-
| ''<math>V_a</math>'' || is the air space (m 3)
|}

====Initial Condition====
The user may define a water content initialization or choose the model to compute in the unsaturated area heads so that soil water is at [[Field Capacity]].
The below explains the latter.

Once determined the aquifer level ('''water table''') the water content is associated at each cells by the following criteria:
:{|
*If the cell is located above the water table <math>\theta=\theta_{s}</math>
|-
*If the cell is located over the water table <math>\theta=\theta_{ns}</math>
|}

where:
:{|
| ''<math>\theta=\theta_{s}</math>'' || is the water content in the saturated soil (m3/m3)
|-
| ''<math>\theta=\theta_{ns}</math>'' || is the water content in the non saturated soil (m3/m3)
|-
|}

In order to calculate field capacity by the equation (1.7) the evaluation of the suction head is needed :

:<math> h=-DWZ\cdot 0.5\hspace{2cm}for\,\, the\,\, cells\,\, immediately\,\, above\,\, the\,\, water\,\, table \hspace{5cm} (1.5)</math>

:<math>h=-(-DZZ-h)\hspace{1.5cm}for\,\, the\,\, other\,\, cells\hspace{9.75cm} (1.6) </math>

[[Image:Figure05.jpg|thumb|center|300px|Figure 1: Suction Head Calculation]]

As shown in the picture the suction head is calculated in order to maintain the same total head (H = z + p + h) in the cells in agreement with the field capacity definition.

===Water retention===
The model use for characterizing the shape of water retention curves is the van Genuchten model:

:<math>\theta(h) = \theta_r + \frac{\theta_s - \theta_r}{\left[ 1+(\alpha |h|)^n \right]^{1-1/n}}\hspace{0.5cm}\Longrightarrow\hspace{0.5cm}h(\theta)=\left |\frac{(S_{E}^{-1/n}-1)^{1/n}}{\alpha} \right |\hspace{3cm}(1.7)</math>

:<math>S_{E}=\frac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}}\hspace{11cm}(1.8)</math>

where:
:{|
|''<math>\theta(h)</math>'' || is the the water retention curve (m 3/m 3)
|-
| ''<math>\theta_s</math>'' || is the saturated water content (m 3/m 3)
|-
| ''<math>\theta_r</math>'' || is the residual water content (m 3/m 3)
|-
| ''h''|| is the suction head (m)
|-
| ''<math>\alpha</math>'' || is related to the inverse of the air entry (m -1)
|-
| ''n''|| is a measure of the pore-size distribution n>1 (dimensionless)
|-
| ''<math>S_{E}</math>''|| is the effective saturation (dimensionless)
|}

===Saturated and Unsaturated Conductivity===

The saturated conductivity is given depending on the type of soil; instead the unsaturated conductivity is obtained from the suction head by the Maulem model:

<math>K(\theta)=K_{s}\cdot Se^{L}\cdot (1-(1-Se^{1/m})^{m})^{2}</math>

where:
:{|
|''<math>K(\theta)</math>'' || is the unsaturated conductivity (m/s)
|-
| ''<math>K_{s}</math>'' || is the saturated conductivity(m/s)
|-
| ''L'' || empirical pore-connectivity (m)
|-
| ''m''|| m=1-1/n(dimensionless)
|}

==Evapotranspiration==

Some water may be extracted from the soil because of the evaporation and transpiration processes, which become a sink in soil water profile. These two processes are currently named Evapotranspiration.

Potential Evapotranspiration may be modeled using the Penmann Monteith equation.
Also, if vegetation exists, a differentiation between Potential Transpiration and Potential Evaporation is done using LAI.
These computation are made in [[Module Basin]] since this is the module that handles water fluxes in the interface betwen modules.

However, not all of the potential water that can be evaporated or transpired will be in fact removed from the soil. The water that will really leave the soil through these processes is calculated: i) if vegetation exists, effective transpiration is computed in module Vegetation; ii) effective evaporation is computed in the Porous Media module.
In Figure below it can be seen that the actual transpiration and evaporation are then used in Porous Media module to compute the new water content.
The actual evaporation, which happens only at the soil surface, is calculated based on: i) a pressure head limit or ii) a soil conductivity limit, chosen by the user. It allows the model not to evaporate any surface water, even if it is available, when the soil head gets below the assigned value (i) or limits evaporation velocity to layer unsaturated conductivity (ii).

[[Image:evapotranspiration fluxogram.jpg|thumb|center|400px|Evapotranspiration fluxogram in Mohid Land model]]
Remind that Feddes is one option for computing effective transpiration in plants in [[module Vegetation]].

==Iteration Process==

Once obtained the water fluxes a balance on the water volume of each cell is apply in order to obtain the new water content <math>\theta</math>. The balance applied is the following:

'''Horizontal Direction X'''

<math>\theta=(\theta\cdot V_{cell}+(FluxU_{(i,j,k)}\cdot ComputeFace_{(i,j,k)}-Flux_{(i,j+1,k)}\cdot ComputeFace_{(i,j+1,k)}\cdot \Delta t) V_{cell}</math>

'''Horizontal Direction Y'''

<math>\theta=(\theta\cdot V_{cell}+(FluxV_{(i,j,k)}\cdot ComputeFace_{(i,j,k)}-FluxV_{(i+1,j,k)}\cdot ComputeFace_{(i+1,j,k)})\cdot \Delta t)\cdot V_{cell}</math>

'''Vertical Direction W'''

<math>\theta=(\theta\cdot V_{cell}+(FluxW_{(i,j,k)}\cdot ComputeFace_{(i,j,k})-FluxW_{(i,j,k+1)}\cdot ComputeFace_{(i,j,k+1}))\cdot \Delta t)\cdot V_{cell}</math>

'''Transpiration Flux'''

<math>\theta=(\theta\cdot V_{cell}-(TranspFlux_{(i,j,k)}\Delta t)\cdot V_{cell}</math>

'''Evaporation Flux'''

<math>\theta=(\theta\cdot V_{cell}-(EvapFlux_{(i,j,k)}\Delta t)\cdot V_{cell}</math>

'''Infiltration Flux'''

<math>\theta=(\theta\cdot V_{cell}-(UnsatK\cdot Area_{cell}\cdot(1-Imp) \Delta t)\cdot V_{cell}</math>

where:

:{|

|''<math>\theta</math>'' || is the water content (m 3/m 3)
|-
|''<math>V_{cell}</math>'' || is the cell volume (m 3)
|-
|''<math>Area_{cell}</math>'' || is the cell area (m 2)
|-
|''ComputeFace'' || is the computed face (dimensionless)
|-
|''<math>\Delta t</math>'' || is the time step (s)
|-
|''FluxX'' || is the flux in X direction (m 3/s)
|-
|''FluxV'' || is the flux in Y direction (m 3/s)
|-
|''FluxW'' || is the flux in W direction (m 3/s)
|-
|''TranspFlux'' || is the transpiration flux (m 3/s)
|-
|''EvapFlux'' || is the evaporation flux (m 3/s)
|-
|''UnsatK'' || is unsaturated conductivity (m /s)
|-
|''Imp'' || is the percentage of impermeable soil of the cell (dimensionless)
|}

The value of <math>\theta</math> so obtained is compared with one used in the volumes calculation and the iterative process stop when:

<math>(\theta^{'})-(\theta^{new})<\,\, Tolerance\hspace{3cm} (1.9) </math>

where:

:{|
|''<math>\theta^{'}</math>'' || is the water content of the previous iteration (m 3/m 3)
|}

If the equation (1.9) is not satisfy the temporal step is divided in half and the new value of <math>\theta</math> is used for solving the equation (1.7) for restarting the calculation process. The iteration process is stoped when the tolerance desired is reached.

[[Image:Figure06.jpg|thumb|center|300px|Figure 2: Time step reduction]]

==Other Features==
===How to Generate Info needed in Porous Media===
====SoilMap====
Model needs to know soil ID in each cell and layer to pick hydraulic properties from that type of soil. In Pedology soil includes more than one horizon, each with different soil properties. Here Soil is used has a unit of soil hidraulic properties, i.e., to define a soil with three horizons one has to create three SoilID (see below). This also means that if a watershed has at least one soil with three horizons one has to create three soil maps. Each soil map will be infact the map of each horizon of the soils. The grid cells with only one horizon will have the same SoilID in all maps, the grid cells with three horizons will have a different SoilID in each map.

Options:
* Constant value
* Soil Grid. One possible option is to associate with soil shape file. In this case can use MOHID GIS going to menu [Tools]->[Shape to Grid Data] and provide: i) the grid (model grid), ii) the soil shape file and iii) the corespondence between soil codes and soil ID defined in data file.

Soil ID must be defined in [[Module_FillMatrix|Module FillMatrix]] standards for each soil horizon defined (grid example):
<beginhorizon>
KLB : 1
KUB : 10

<beginproperty>
NAME : SoilID
DEFAULTVALUE : 1
INITIALIZATION_METHOD : ASCII_FILE
REMAIN_CONSTANT : 1
FILENAME : ..\..\GeneralData\PorousMedia\SoilID200m.dat
<endproperty>

<beginproperty>
<endproperty>
..
<endhorizon>

'''Remarks'''

All the soil ID's appearing in the soil grid(s) must be defined in the PorousMedia data file in terms of hydraulic properties:
<beginsoiltype>
ID : 1
THETA_S : 0.3859
THETA_R : 0.0476
N_FIT : 1.39
SAT_K : 3.5556e-6
ALPHA : 2.75
L_FIT : 0.50
<endsoiltype>

<beginsoiltype>
ID : 2
...
<endsoiltype>
...

====Soil Bottom====
The soil depth must be known by the model. This is computed by the model from terrain altitude (topography) and soil bottom altitude. As so, a soil bottom grid is needed.

Options:
* Grid File.
Soil depth (and soil bottom altitude, the effective grid needed) can be defined with a constant depth or estimated from slope [[HOW TO SoilBottom LINK]]. When the soil depth is estimated as a function of slope, soil depth will be smaller in ares with higher slope. In this areas only the surface layers of the soil will be considered (see [[Porous Media Geometry]]).

Define the grid just generated, in the porous media data file with:
BOTTOM_FILE : ..\..\GeneralData\PorousMedia\BottomLevel.dat

====Water Level====
Options:
*Grid File.
The water table altitude represents the initial altitude of the water table.
It is recommended to do a spin-up run to estabilize water level and then do a continuous simulation starting with the final water table achieved.
Use the following blocks with [[Module_FillMatrix|Module FillMatrix]] standards:
<beginwaterlevel>
NAME : waterlevel
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 0
FILENAME : ..\..\GeneralData\PorousMedia\WaterLevel0.50.dat
<endwaterlevel>

====Impermeability====
Impermeability values (0 - completely permeable, 1 - impermeable) must be provided.

Options:
* Constant Value.
* Grid File. One possible option is to associate with land use shape file. In this case can use MOHID GIS going to menu [Tools]->[Shape to Grid Data] and provide: i) the grid (model grid), ii) the land use shape file and iii) the corespondence between land use codes and Impermeability values.
Use the following blocks with [[Module_FillMatrix|Module FillMatrix]] standards:
<beginimpermeablefraction>
NAME : impermeablefraction
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 1
FILENAME : ..\..\GeneralData\PorousMedia\AreaImpermeavel.dat
<endimpermeablefraction>

==Outputs==

===Timeseries===

Theta - is water content of selected cell (vol water/ vol soil)

relative_water_content - content in selected cell. between zero and one (zero is residual water content and one is saturated water content)

VelW_[m/s] - vertical velocity in the bottom face of the selected cell

VelW_Corr_[m/s] - vertical velocity in the bottom face of the selected cell that may be corrected if oversaturation occurs. if no correction occurs is the same as previous.

InF_Vel_[m/s] - infiltration velocity (in the soil surface)

Head_[m] - Suction in selected cell

Conductivity_[m/s] - Conductivity in selected cell

level_water_table_[m] - water table altitude

water_table_depth_[m] - water table depth (from soil surface)

Hydro_Pressure_[m] - hydrostatic pressure in selected cell

Final_Head_[m] - Soil water charge in selected cell

[Check Mohid Land Heights and Levels to understand some of the outputs]

GW_flow_to_river_total_[m3/s] - Ground water flow to river if it is a river point

Surface_Evaporation_Flux_[m3/s] - evaporation flux (in the soil surface)

Transpiration_Flux_[m3/s] - transpiration flux in selected cell

[http://screencast.com/t/I8We0B4j2uJ Check Mohid Land Heights and Levels to understand the difference between level and height results]

==References ==
*Jury,W.A.,Gardner,W.R.,Gardner,W.H., 1991,Soil Physics
*Van Genuchten, M.T., A closed form equation for predicting the hydraulic conductivity of unsaturated soils
*Wu,J.,Zhang, R., Gui,S.,1999, Modelling soil water movement with water uptake by roots, Plant and soil 215: 7-17
*Marcel G.Schaap and Martinus Th. van Genuchten, A modified Maulem van Genuchten Formulation for Improved Description of Hydraulic Conductivity Near Saturation, 16 December 2005

==Data File ==

===Keywords===
Keywords read in the Data File

Keyword : Data Type Default !Comment

BOTTOM_FILE : char - !Path to Bottom Topography File
START_WITH_FIELD : logical 1 !Sets Theta initial Field Capacity
CONTINUOUS : logical 0 !Continues from previous run
STOP_ON_WRONG_DATE : logical 1 !Stops if previous run end is different from actual
!Start
OUTPUT_TIME : sec. sec. sec. - !Output Time
TIME_SERIE_LOCATION : char - !Path to File which defines Time Series
CONTINUOUS_OUTPUT_FILE : logical 1 !Writes "famous" iter.log
CONDUTIVITYFACE : integer 1 !Way to interpolate conducivity face
!1 - Average, 2 - Maximum, 3 - Minimum, 4 - Weigthed
HORIZONTAL_K_FACTOR : real 1.0 !Factor for Horizontal Conductivity = Kh / Kv
CUT_OFF_THETA_LOW : real 1e-6 !Disables calculation when Theta is near ThetaR
CUT_OFF_THETA_HIGH : real 1e-15 !Set Theta = ThetaS when Theta > ThetaS - CUT_OFF_THETA_HIGH
MIN_ITER : integer 2 !Number of iterations below which the DT is increased
MAX_ITER : integer 3 !Number of iterations above which the DT is decreased
LIMIT_ITER : integer 50 !Number of iterations of a time step (for restart)
THETA_TOLERANCE : real 0.001 !Converge Parameter
INCREASE_DT : real 1.25 !Increase of DT when iter < MIN_ITER
DECREASE_DT : real 0.70 !Decrease of DT when iter > MAX_ITER

<beginproperty>
NAME : Theta / waterlevel

see Module FillMatrix for more options

<endproperty>

===Sample===

Some keywords of the PorousMedia input file:
BOTTOM_FILE : ..\General Data\Other\PorousMedia\SoilBottom_2cells.dat
TIME_SERIE_LOCATION : ..\General Data\TimeSeries\TimeSeriesLocation3D_2m.dat

COMPUTE_SOIL_FIELD : 1

OUTPUT_TIME : 0 86400
CUT_OFF_THETA_HIGH : 1e-15
START_WITH_FIELD : 1

LIMIT_EVAP_WATER_VEL : 0
LIMIT_EVAP_HEAD : 0
HEAD_LIMIT : -100
THETA_HYDRO_COEF : 0.98

<beginsoiltype>
ID : 1
THETA_S : 0.43
THETA_R : 0.078
SAT_K : 2.888e-6
N_FIT : 1.56
ALPHA : 3.6
L_FIT : 0.50
<endsoiltype>

!----- Hydraulic Soil Properties
<beginhorizon>
KLB : 1
KUB : 10

<beginproperty>
NAME : SoilID
INITIALIZATION_METHOD : CONSTANT
DEFAULTVALUE : 1
<endproperty>

<beginproperty>
NAME : Theta
INITIALIZATION_METHOD : CONSTANT
DEFAULTVALUE : 0.30
<endproperty>
<endhorizon>

<beginwaterlevel>
NAME : waterlevel
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 1.
REMAIN_CONSTANT : 1
FILENAME : ..\General Data\Initial Conditions\InitialWaterLevel_2cells.dat
<endwaterlevel>

<beginimpermeablefraction>
NAME : impermeablefraction
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 1
FILENAME : ..\General Data\Other\PorousMedia\InitialImpermeabilization_2cells.dat
<endimpermeablefraction>

[[Category:MOHID Land]]

Module PorousMedia

2012-01-17T18:26:13Z

Carina:

==Overview==
Module Porous Media is responsible for handling all water fluxes in soil. In particular the infiltration flow is obtained from the potential water column available for this process. The potential water column is given by the [[Module Basin]] (the module that handles interface between modules) as reported below:

<math>PIC=TF+IWC</math>

where:
:{|
| ''PIC'' || is the Potential Infiltration column (m)
|-
| ''TF'' || Throughfall is the percentage of the rain that reaches the soil (m)
|-
| ''IWC'' || is the Initial Water column (m)
|}

The infiltration flux is calculated by the Buckingham-Darcy equation using head gradient between surface runoff and first soil layer (Jury et al,1991).

==Main Processes==

===Porous Media Geometry===
Porous Media is a 3D domain delimited in its upper limit by topography and lower limit by soil bottom (defined by user).
In terms of soil definition it can be defined vertical horizons to link correspond to real soil horizons with different hydraulic carachteristics.

<htm><a href="http://content.screencast.com/users/jovem/folders/Jing/media/b9ef6c79-c1d1-46ea-94f7-623471c15883/MohidLandSoilProfile.png"><img src="http://content.screencast.com/users/jovem/folders/Jing/media/b9ef6c79-c1d1-46ea-94f7-623471c15883/MohidLandSoilProfile.png" width="650" height="436" border="0" /></a></htm>

===Water Flow===
Soil contains a large distribution of pore sizes and channels through which water may flow. In general the water flow determination is based on the mass conservation and momentum equation. In the case of soil it is assumed that the forces of inertia are almost zero; therefore the balance is reduced to the forces of pressure, gravity and viscous. The equation that describes the flow through soil is the Buckingham Darcy equation (Jury et al,1991).

<math>v=-K\left ( \theta \right )\left ( \frac{\partial H}{\partial z} \right )=-K\left ( \theta \right )\left ( \frac{\partial h}{\partial z} + 1 \right )\,\,\,\,(1.1)\\ </math>

where:
:{|
|''v'' || is the water velocity at the cell interface (m/s)
|-
| ''H'' || is the hydraulic head (m)
|-
| ''θ'' || is the water content (m3/m3)
|-
| ''K''|| is the hydraulic conductivity (m/s)
|}

The hydraulic head is given by the formula:

<math>H=h+p+z\,\,\,\,(1.2)</math>

where:
:{|
|-
| ''h'' || is the hydraulic head (m)
|-
| ''p'' || is hydrostatic pressure (m)
|-
| ''z''|| is the topography (m)
|-
|}

When in saturated conditions, hydraulic head is zero and hydrostatic pressure may occur (if water is at rest or decelerating). In unsaturated conditions hydrostatic pressure is zero and hydarulic head exists.

The soil is a very complex system, made up of a heterogeneous mixture of solid, liquid, and gaseous material. The liquid phase consists of soil water, which fills part or all of the open spaces between the soil particles. Therefore it is possible to divided the soil in two layers:

*Saturated soil <math>\Longrightarrow</math> The soil pores are filled by water

*Unsaturated one <math>\Longrightarrow</math> The soil pores are filled by water and air

In the first case the equation of Buckingham Darcy is simplified to the Darcy law and the parameter associated for its resolution are connected with the saturated layer. On the other hand for the resolution of the equation (1.1) a description of the characteristics of the the unsaturated layer is needed.

==Vadose Zone==

Many vadose flow and transport studies require description of unsaturated soil hydraulic proprieties over a wide range of pressure head. The hydraulic proprieties are described using the porous size distribution model of Maulem (1976) for hydraulic conductivity in combination with a water retention function introduced by Van Genuchten (1980).

===Water content===

Water content is the quantity of water contained in the soil (called '''soil moisture'''). It is given as a volumetric basis and it is defined mathematically as:

:<math>\theta = \frac{V_w}{V_T}\hspace{5cm}(1.3) </math>

:<math>V_T = V_s + V_v = V_s + V_w + V_a\hspace{1.6cm}(1.4)</math>

where:
:{|
| ''<math>\theta</math>'' || is water content (m 3/m 3)
|-
|''<math>V_w</math> '' || is the volume of water (m s3)
|-
| ''<math>V_T</math>'' || is the total volume (m 3)
|-
| ''<math>V_s</math>'' || is the soil volume (m 3)
|-
| ''<math>V_a</math>'' || is the air space (m 3)
|}

====Initial Condition====
The user may define a water content initialization or choose the model to compute in the unsaturated area heads so that soil water is at [[Field Capacity]].
The below explains the latter.

Once determined the aquifer level ('''water table''') the water content is associated at each cells by the following criteria:
:{|
*If the cell is located above the water table <math>\theta=\theta_{s}</math>
|-
*If the cell is located over the water table <math>\theta=\theta_{ns}</math>
|}

where:
:{|
| ''<math>\theta=\theta_{s}</math>'' || is the water content in the saturated soil (m3/m3)
|-
| ''<math>\theta=\theta_{ns}</math>'' || is the water content in the non saturated soil (m3/m3)
|-
|}

In order to calculate field capacity by the equation (1.7) the evaluation of the suction head is needed :

:<math> h=-DWZ\cdot 0.5\hspace{2cm}for\,\, the\,\, cells\,\, immediately\,\, above\,\, the\,\, water\,\, table \hspace{5cm} (1.5)</math>

:<math>h=-(-DZZ-h)\hspace{1.5cm}for\,\, the\,\, other\,\, cells\hspace{9.75cm} (1.6) </math>

[[Image:Figure05.jpg|thumb|center|300px|Figure 1: Suction Head Calculation]]

As shown in the picture the suction head is calculated in order to maintain the same total head (H = z + p + h) in the cells in agreement with the field capacity definition.

===Water retention===
The model use for characterizing the shape of water retention curves is the van Genuchten model:

:<math>\theta(h) = \theta_r + \frac{\theta_s - \theta_r}{\left[ 1+(\alpha |h|)^n \right]^{1-1/n}}\hspace{0.5cm}\Longrightarrow\hspace{0.5cm}h(\theta)=\left |\frac{(S_{E}^{-1/n}-1)^{1/n}}{\alpha} \right |\hspace{3cm}(1.7)</math>

:<math>S_{E}=\frac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}}\hspace{11cm}(1.8)</math>

where:
:{|
|''<math>\theta(h)</math>'' || is the the water retention curve (m 3/m 3)
|-
| ''<math>\theta_s</math>'' || is the saturated water content (m 3/m 3)
|-
| ''<math>\theta_r</math>'' || is the residual water content (m 3/m 3)
|-
| ''h''|| is the suction head (m)
|-
| ''<math>\alpha</math>'' || is related to the inverse of the air entry (m -1)
|-
| ''n''|| is a measure of the pore-size distribution n>1 (dimensionless)
|-
| ''<math>S_{E}</math>''|| is the effective saturation (dimensionless)
|}

===Saturated and Unsaturated Conductivity===

The saturated conductivity is given depending on the type of soil; instead the unsaturated conductivity is obtained from the suction head by the Maulem model:

<math>K(\theta)=K_{s}\cdot Se^{L}\cdot (1-(1-Se^{1/m})^{m})^{2}</math>

where:
:{|
|''<math>K(\theta)</math>'' || is the unsaturated conductivity (m/s)
|-
| ''<math>K_{s}</math>'' || is the saturated conductivity(m/s)
|-
| ''L'' || empirical pore-connectivity (m)
|-
| ''m''|| m=1-1/n(dimensionless)
|}

==Evapotranspiration==

Some water may be extracted from the soil because of the evaporation and transpiration processes, which become a sink in soil water profile. These two processes are currently named Evapotranspiration.

Potential Evapotranspiration may be modeled using the Penmann Monteith equation.
Also, if vegetation exists, a differentiation between Potential Transpiration and Potential Evaporation is done using LAI.
These computation are made in [[Module Basin]] since this is the module that handles water fluxes in the interface betwen modules.

However, not all of the potential water that can be evaporated or transpired will be in fact removed from the soil. The water that will really leave the soil through these processes is calculated: i) if vegetation exists, effective transpiration is computed in module Vegetation; ii) effective evaporation is computed in the Porous Media module.
In Figure below it can be seen that the actual transpiration and evaporation are then used in Porous Media module to compute the new water content.
The actual evaporation, which happens only at the soil surface, is calculated based on: i) a pressure head limit or ii) a soil conductivity limit, chosen by the user. It allows the model not to evaporate any surface water, even if it is available, when the soil head gets below the assigned value (i) or limits evaporation velocity to layer unsaturated conductivity (ii).

[[Image:evapotranspiration fluxogram.jpg|thumb|center|400px|Evapotranspiration fluxogram in Mohid Land model]]
Remind that Feddes is one option for computing effective transpiration in plants in [[module Vegetation]].

==Iteration Process==

Once obtained the water fluxes a balance on the water volume of each cell is apply in order to obtain the new water content <math>\theta</math>. The balance applied is the following:

'''Horizontal Direction X'''

<math>\theta=(\theta\cdot V_{cell}+(FluxU_{(i,j,k)}\cdot ComputeFace_{(i,j,k)}-Flux_{(i,j+1,k)}\cdot ComputeFace_{(i,j+1,k)}\cdot \Delta t) V_{cell}</math>

'''Horizontal Direction Y'''

<math>\theta=(\theta\cdot V_{cell}+(FluxV_{(i,j,k)}\cdot ComputeFace_{(i,j,k)}-FluxV_{(i+1,j,k)}\cdot ComputeFace_{(i+1,j,k)})\cdot \Delta t)\cdot V_{cell}</math>

'''Vertical Direction W'''

<math>\theta=(\theta\cdot V_{cell}+(FluxW_{(i,j,k)}\cdot ComputeFace_{(i,j,k})-FluxW_{(i,j,k+1)}\cdot ComputeFace_{(i,j,k+1}))\cdot \Delta t)\cdot V_{cell}</math>

'''Transpiration Flux'''

<math>\theta=(\theta\cdot V_{cell}-(TranspFlux_{(i,j,k)}\Delta t)\cdot V_{cell}</math>

'''Evaporation Flux'''

<math>\theta=(\theta\cdot V_{cell}-(EvapFlux_{(i,j,k)}\Delta t)\cdot V_{cell}</math>

'''Infiltration Flux'''

<math>\theta=(\theta\cdot V_{cell}-(UnsatK\cdot Area_{cell}\cdot(1-Imp) \Delta t)\cdot V_{cell}</math>

where:

:{|

|''<math>\theta</math>'' || is the water content (m 3/m 3)
|-
|''<math>V_{cell}</math>'' || is the cell volume (m 3)
|-
|''<math>Area_{cell}</math>'' || is the cell area (m 2)
|-
|''ComputeFace'' || is the computed face (dimensionless)
|-
|''<math>\Delta t</math>'' || is the time step (s)
|-
|''FluxX'' || is the flux in X direction (m 3/s)
|-
|''FluxV'' || is the flux in Y direction (m 3/s)
|-
|''FluxW'' || is the flux in W direction (m 3/s)
|-
|''TranspFlux'' || is the transpiration flux (m 3/s)
|-
|''EvapFlux'' || is the evaporation flux (m 3/s)
|-
|''UnsatK'' || is unsaturated conductivity (m /s)
|-
|''Imp'' || is the percentage of impermeable soil of the cell (dimensionless)
|}

The value of <math>\theta</math> so obtained is compared with one used in the volumes calculation and the iterative process stop when:

<math>(\theta^{'})-(\theta^{new})<\,\, Tolerance\hspace{3cm} (1.9) </math>

where:

:{|
|''<math>\theta^{'}</math>'' || is the water content of the previous iteration (m 3/m 3)
|}

If the equation (1.9) is not satisfy the temporal step is divided in half and the new value of <math>\theta</math> is used for solving the equation (1.7) for restarting the calculation process. The iteration process is stoped when the tolerance desired is reached.

[[Image:Figure06.jpg|thumb|center|300px|Figure 2: Time step reduction]]

==Other Features==
===How to Generate Info needed in Porous Media===
====SoilMap====
Model needs to know soil ID in each cell and layer to pick hydraulic properties from that type of soil. In Pedology soil includes more than one horizon, each with different soil properties. Here Soil is used has a unit of soil hidraulic properties, i.e., to define a soil with three horizons one has to create three SoilID (see below). This also means that if a watershed has at least one soil with three horizons one has to create three soil maps. Each soil map will be infact the map of each horizon of the soils. The grid cells with only one horizon will have the same SoilID in all maps, the grid cells with three horizons will have a different SoilID in each map.

Options:
* Constant value
* Soil Grid. One possible option is to associate with soil shape file. In this case can use MOHID GIS going to menu [Tools]->[Shape to Grid Data] and provide: i) the grid (model grid), ii) the soil shape file and iii) the corespondence between soil codes and soil ID defined in data file.

Soil ID must be defined in [[Module_FillMatrix|Module FillMatrix]] standards for each soil horizon defined (grid example):
<beginhorizon>
KLB : 1
KUB : 10

<beginproperty>
NAME : SoilID
DEFAULTVALUE : 1
INITIALIZATION_METHOD : ASCII_FILE
REMAIN_CONSTANT : 1
FILENAME : ..\..\GeneralData\PorousMedia\SoilID200m.dat
<endproperty>

<beginproperty>
<endproperty>
..
<endhorizon>

'''Remarks'''

All the soil ID's appearing in the soil grid(s) must be defined in the PorousMedia data file in terms of hydraulic properties:
<beginsoiltype>
ID : 1
THETA_S : 0.3859
THETA_R : 0.0476
N_FIT : 1.39
SAT_K : 3.5556e-6
ALPHA : 2.75
L_FIT : 0.50
<endsoiltype>

<beginsoiltype>
ID : 2
...
<endsoiltype>
...

====Soil Bottom====
The soil depth must be known by the model. This is computed by the model from terrain altitude (topography) and soil bottom altitude. As so, a soil bottom grid is needed.

Options:
* Grid File.
Soil depth (and soil bottom altitude, the effective grid needed) can be defined with a constant depth or estimated from slope [[HOW TO SoilBottom LINK]]. When the soil depth is estimated as a function of slope, soil depth will be smaller in ares with higher slope. In this areas only the surface layers of the soil will be considered (see [Porous Media Geometry]).

Define the grid just generated, in the porous media data file with:
BOTTOM_FILE : ..\..\GeneralData\PorousMedia\BottomLevel.dat

====Water Level====
Options:
*Grid File.
The water table altitude represents the initial altitude of the water table.
It is recommended to do a spin-up run to estabilize water level and then do a continuous simulation starting with the final water table achieved.
Use the following blocks with [[Module_FillMatrix|Module FillMatrix]] standards:
<beginwaterlevel>
NAME : waterlevel
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 0
FILENAME : ..\..\GeneralData\PorousMedia\WaterLevel0.50.dat
<endwaterlevel>

====Impermeability====
Impermeability values (0 - completely permeable, 1 - impermeable) must be provided.

Options:
* Constant Value.
* Grid File. One possible option is to associate with land use shape file. In this case can use MOHID GIS going to menu [Tools]->[Shape to Grid Data] and provide: i) the grid (model grid), ii) the land use shape file and iii) the corespondence between land use codes and Impermeability values.
Use the following blocks with [[Module_FillMatrix|Module FillMatrix]] standards:
<beginimpermeablefraction>
NAME : impermeablefraction
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 1
FILENAME : ..\..\GeneralData\PorousMedia\AreaImpermeavel.dat
<endimpermeablefraction>

==Outputs==

===Timeseries===

Theta - is water content of selected cell (vol water/ vol soil)

relative_water_content - content in selected cell. between zero and one (zero is residual water content and one is saturated water content)

VelW_[m/s] - vertical velocity in the bottom face of the selected cell

VelW_Corr_[m/s] - vertical velocity in the bottom face of the selected cell that may be corrected if oversaturation occurs. if no correction occurs is the same as previous.

InF_Vel_[m/s] - infiltration velocity (in the soil surface)

Head_[m] - Suction in selected cell

Conductivity_[m/s] - Conductivity in selected cell

level_water_table_[m] - water table altitude

water_table_depth_[m] - water table depth (from soil surface)

Hydro_Pressure_[m] - hydrostatic pressure in selected cell

Final_Head_[m] - Soil water charge in selected cell

[Check Mohid Land Heights and Levels to understand some of the outputs]

GW_flow_to_river_total_[m3/s] - Ground water flow to river if it is a river point

Surface_Evaporation_Flux_[m3/s] - evaporation flux (in the soil surface)

Transpiration_Flux_[m3/s] - transpiration flux in selected cell

[http://screencast.com/t/I8We0B4j2uJ Check Mohid Land Heights and Levels to understand the difference between level and height results]

==References ==
*Jury,W.A.,Gardner,W.R.,Gardner,W.H., 1991,Soil Physics
*Van Genuchten, M.T., A closed form equation for predicting the hydraulic conductivity of unsaturated soils
*Wu,J.,Zhang, R., Gui,S.,1999, Modelling soil water movement with water uptake by roots, Plant and soil 215: 7-17
*Marcel G.Schaap and Martinus Th. van Genuchten, A modified Maulem van Genuchten Formulation for Improved Description of Hydraulic Conductivity Near Saturation, 16 December 2005

==Data File ==

===Keywords===
Keywords read in the Data File

Keyword : Data Type Default !Comment

BOTTOM_FILE : char - !Path to Bottom Topography File
START_WITH_FIELD : logical 1 !Sets Theta initial Field Capacity
CONTINUOUS : logical 0 !Continues from previous run
STOP_ON_WRONG_DATE : logical 1 !Stops if previous run end is different from actual
!Start
OUTPUT_TIME : sec. sec. sec. - !Output Time
TIME_SERIE_LOCATION : char - !Path to File which defines Time Series
CONTINUOUS_OUTPUT_FILE : logical 1 !Writes "famous" iter.log
CONDUTIVITYFACE : integer 1 !Way to interpolate conducivity face
!1 - Average, 2 - Maximum, 3 - Minimum, 4 - Weigthed
HORIZONTAL_K_FACTOR : real 1.0 !Factor for Horizontal Conductivity = Kh / Kv
CUT_OFF_THETA_LOW : real 1e-6 !Disables calculation when Theta is near ThetaR
CUT_OFF_THETA_HIGH : real 1e-15 !Set Theta = ThetaS when Theta > ThetaS - CUT_OFF_THETA_HIGH
MIN_ITER : integer 2 !Number of iterations below which the DT is increased
MAX_ITER : integer 3 !Number of iterations above which the DT is decreased
LIMIT_ITER : integer 50 !Number of iterations of a time step (for restart)
THETA_TOLERANCE : real 0.001 !Converge Parameter
INCREASE_DT : real 1.25 !Increase of DT when iter < MIN_ITER
DECREASE_DT : real 0.70 !Decrease of DT when iter > MAX_ITER

<beginproperty>
NAME : Theta / waterlevel

see Module FillMatrix for more options

<endproperty>

===Sample===

Some keywords of the PorousMedia input file:
BOTTOM_FILE : ..\General Data\Other\PorousMedia\SoilBottom_2cells.dat
TIME_SERIE_LOCATION : ..\General Data\TimeSeries\TimeSeriesLocation3D_2m.dat

COMPUTE_SOIL_FIELD : 1

OUTPUT_TIME : 0 86400
CUT_OFF_THETA_HIGH : 1e-15
START_WITH_FIELD : 1

LIMIT_EVAP_WATER_VEL : 0
LIMIT_EVAP_HEAD : 0
HEAD_LIMIT : -100
THETA_HYDRO_COEF : 0.98

<beginsoiltype>
ID : 1
THETA_S : 0.43
THETA_R : 0.078
SAT_K : 2.888e-6
N_FIT : 1.56
ALPHA : 3.6
L_FIT : 0.50
<endsoiltype>

!----- Hydraulic Soil Properties
<beginhorizon>
KLB : 1
KUB : 10

<beginproperty>
NAME : SoilID
INITIALIZATION_METHOD : CONSTANT
DEFAULTVALUE : 1
<endproperty>

<beginproperty>
NAME : Theta
INITIALIZATION_METHOD : CONSTANT
DEFAULTVALUE : 0.30
<endproperty>
<endhorizon>

<beginwaterlevel>
NAME : waterlevel
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 1.
REMAIN_CONSTANT : 1
FILENAME : ..\General Data\Initial Conditions\InitialWaterLevel_2cells.dat
<endwaterlevel>

<beginimpermeablefraction>
NAME : impermeablefraction
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 1
FILENAME : ..\General Data\Other\PorousMedia\InitialImpermeabilization_2cells.dat
<endimpermeablefraction>

[[Category:MOHID Land]]

Module PorousMedia

2012-01-17T18:25:33Z

Carina:

==Overview==
Module Porous Media is responsible for handling all water fluxes in soil. In particular the infiltration flow is obtained from the potential water column available for this process. The potential water column is given by the [[Module Basin]] (the module that handles interface between modules) as reported below:

<math>PIC=TF+IWC</math>

where:
:{|
| ''PIC'' || is the Potential Infiltration column (m)
|-
| ''TF'' || Throughfall is the percentage of the rain that reaches the soil (m)
|-
| ''IWC'' || is the Initial Water column (m)
|}

The infiltration flux is calculated by the Buckingham-Darcy equation using head gradient between surface runoff and first soil layer (Jury et al,1991).

==Main Processes==

===Porous Media Geometry===
Porous Media is a 3D domain delimited in its upper limit by topography and lower limit by soil bottom (defined by user).
In terms of soil definition it can be defined vertical horizons to link correspond to real soil horizons with different hydraulic carachteristics.

<htm><a href="http://content.screencast.com/users/jovem/folders/Jing/media/b9ef6c79-c1d1-46ea-94f7-623471c15883/MohidLandSoilProfile.png"><img src="http://content.screencast.com/users/jovem/folders/Jing/media/b9ef6c79-c1d1-46ea-94f7-623471c15883/MohidLandSoilProfile.png" width="650" height="436" border="0" /></a></htm>

===Water Flow===
Soil contains a large distribution of pore sizes and channels through which water may flow. In general the water flow determination is based on the mass conservation and momentum equation. In the case of soil it is assumed that the forces of inertia are almost zero; therefore the balance is reduced to the forces of pressure, gravity and viscous. The equation that describes the flow through soil is the Buckingham Darcy equation (Jury et al,1991).

<math>v=-K\left ( \theta \right )\left ( \frac{\partial H}{\partial z} \right )=-K\left ( \theta \right )\left ( \frac{\partial h}{\partial z} + 1 \right )\,\,\,\,(1.1)\\ </math>

where:
:{|
|''v'' || is the water velocity at the cell interface (m/s)
|-
| ''H'' || is the hydraulic head (m)
|-
| ''θ'' || is the water content (m3/m3)
|-
| ''K''|| is the hydraulic conductivity (m/s)
|}

The hydraulic head is given by the formula:

<math>H=h+p+z\,\,\,\,(1.2)</math>

where:
:{|
|-
| ''h'' || is the hydraulic head (m)
|-
| ''p'' || is hydrostatic pressure (m)
|-
| ''z''|| is the topography (m)
|-
|}

When in saturated conditions, hydraulic head is zero and hydrostatic pressure may occur (if water is at rest or decelerating). In unsaturated conditions hydrostatic pressure is zero and hydarulic head exists.

The soil is a very complex system, made up of a heterogeneous mixture of solid, liquid, and gaseous material. The liquid phase consists of soil water, which fills part or all of the open spaces between the soil particles. Therefore it is possible to divided the soil in two layers:

*Saturated soil <math>\Longrightarrow</math> The soil pores are filled by water

*Unsaturated one <math>\Longrightarrow</math> The soil pores are filled by water and air

In the first case the equation of Buckingham Darcy is simplified to the Darcy law and the parameter associated for its resolution are connected with the saturated layer. On the other hand for the resolution of the equation (1.1) a description of the characteristics of the the unsaturated layer is needed.

==Vadose Zone==

Many vadose flow and transport studies require description of unsaturated soil hydraulic proprieties over a wide range of pressure head. The hydraulic proprieties are described using the porous size distribution model of Maulem (1976) for hydraulic conductivity in combination with a water retention function introduced by Van Genuchten (1980).

===Water content===

Water content is the quantity of water contained in the soil (called '''soil moisture'''). It is given as a volumetric basis and it is defined mathematically as:

:<math>\theta = \frac{V_w}{V_T}\hspace{5cm}(1.3) </math>

:<math>V_T = V_s + V_v = V_s + V_w + V_a\hspace{1.6cm}(1.4)</math>

where:
:{|
| ''<math>\theta</math>'' || is water content (m 3/m 3)
|-
|''<math>V_w</math> '' || is the volume of water (m s3)
|-
| ''<math>V_T</math>'' || is the total volume (m 3)
|-
| ''<math>V_s</math>'' || is the soil volume (m 3)
|-
| ''<math>V_a</math>'' || is the air space (m 3)
|}

====Initial Condition====
The user may define a water content initialization or choose the model to compute in the unsaturated area heads so that soil water is at [[Field Capacity]].
The below explains the latter.

Once determined the aquifer level ('''water table''') the water content is associated at each cells by the following criteria:
:{|
*If the cell is located above the water table <math>\theta=\theta_{s}</math>
|-
*If the cell is located over the water table <math>\theta=\theta_{ns}</math>
|}

where:
:{|
| ''<math>\theta=\theta_{s}</math>'' || is the water content in the saturated soil (m3/m3)
|-
| ''<math>\theta=\theta_{ns}</math>'' || is the water content in the non saturated soil (m3/m3)
|-
|}

In order to calculate field capacity by the equation (1.7) the evaluation of the suction head is needed :

:<math> h=-DWZ\cdot 0.5\hspace{2cm}for\,\, the\,\, cells\,\, immediately\,\, above\,\, the\,\, water\,\, table \hspace{5cm} (1.5)</math>

:<math>h=-(-DZZ-h)\hspace{1.5cm}for\,\, the\,\, other\,\, cells\hspace{9.75cm} (1.6) </math>

[[Image:Figure05.jpg|thumb|center|300px|Figure 1: Suction Head Calculation]]

As shown in the picture the suction head is calculated in order to maintain the same total head (H = z + p + h) in the cells in agreement with the field capacity definition.

===Water retention===
The model use for characterizing the shape of water retention curves is the van Genuchten model:

:<math>\theta(h) = \theta_r + \frac{\theta_s - \theta_r}{\left[ 1+(\alpha |h|)^n \right]^{1-1/n}}\hspace{0.5cm}\Longrightarrow\hspace{0.5cm}h(\theta)=\left |\frac{(S_{E}^{-1/n}-1)^{1/n}}{\alpha} \right |\hspace{3cm}(1.7)</math>

:<math>S_{E}=\frac{\theta-\theta_{r}}{\theta_{s}-\theta_{r}}\hspace{11cm}(1.8)</math>

where:
:{|
|''<math>\theta(h)</math>'' || is the the water retention curve (m 3/m 3)
|-
| ''<math>\theta_s</math>'' || is the saturated water content (m 3/m 3)
|-
| ''<math>\theta_r</math>'' || is the residual water content (m 3/m 3)
|-
| ''h''|| is the suction head (m)
|-
| ''<math>\alpha</math>'' || is related to the inverse of the air entry (m -1)
|-
| ''n''|| is a measure of the pore-size distribution n>1 (dimensionless)
|-
| ''<math>S_{E}</math>''|| is the effective saturation (dimensionless)
|}

===Saturated and Unsaturated Conductivity===

The saturated conductivity is given depending on the type of soil; instead the unsaturated conductivity is obtained from the suction head by the Maulem model:

<math>K(\theta)=K_{s}\cdot Se^{L}\cdot (1-(1-Se^{1/m})^{m})^{2}</math>

where:
:{|
|''<math>K(\theta)</math>'' || is the unsaturated conductivity (m/s)
|-
| ''<math>K_{s}</math>'' || is the saturated conductivity(m/s)
|-
| ''L'' || empirical pore-connectivity (m)
|-
| ''m''|| m=1-1/n(dimensionless)
|}

==Evapotranspiration==

Some water may be extracted from the soil because of the evaporation and transpiration processes, which become a sink in soil water profile. These two processes are currently named Evapotranspiration.

Potential Evapotranspiration may be modeled using the Penmann Monteith equation.
Also, if vegetation exists, a differentiation between Potential Transpiration and Potential Evaporation is done using LAI.
These computation are made in [[Module Basin]] since this is the module that handles water fluxes in the interface betwen modules.

However, not all of the potential water that can be evaporated or transpired will be in fact removed from the soil. The water that will really leave the soil through these processes is calculated: i) if vegetation exists, effective transpiration is computed in module Vegetation; ii) effective evaporation is computed in the Porous Media module.
In Figure below it can be seen that the actual transpiration and evaporation are then used in Porous Media module to compute the new water content.
The actual evaporation, which happens only at the soil surface, is calculated based on: i) a pressure head limit or ii) a soil conductivity limit, chosen by the user. It allows the model not to evaporate any surface water, even if it is available, when the soil head gets below the assigned value (i) or limits evaporation velocity to layer unsaturated conductivity (ii).

[[Image:evapotranspiration fluxogram.jpg|thumb|center|400px|Evapotranspiration fluxogram in Mohid Land model]]
Remind that Feddes is one option for computing effective transpiration in plants in [[module Vegetation]].

==Iteration Process==

Once obtained the water fluxes a balance on the water volume of each cell is apply in order to obtain the new water content <math>\theta</math>. The balance applied is the following:

'''Horizontal Direction X'''

<math>\theta=(\theta\cdot V_{cell}+(FluxU_{(i,j,k)}\cdot ComputeFace_{(i,j,k)}-Flux_{(i,j+1,k)}\cdot ComputeFace_{(i,j+1,k)}\cdot \Delta t) V_{cell}</math>

'''Horizontal Direction Y'''

<math>\theta=(\theta\cdot V_{cell}+(FluxV_{(i,j,k)}\cdot ComputeFace_{(i,j,k)}-FluxV_{(i+1,j,k)}\cdot ComputeFace_{(i+1,j,k)})\cdot \Delta t)\cdot V_{cell}</math>

'''Vertical Direction W'''

<math>\theta=(\theta\cdot V_{cell}+(FluxW_{(i,j,k)}\cdot ComputeFace_{(i,j,k})-FluxW_{(i,j,k+1)}\cdot ComputeFace_{(i,j,k+1}))\cdot \Delta t)\cdot V_{cell}</math>

'''Transpiration Flux'''

<math>\theta=(\theta\cdot V_{cell}-(TranspFlux_{(i,j,k)}\Delta t)\cdot V_{cell}</math>

'''Evaporation Flux'''

<math>\theta=(\theta\cdot V_{cell}-(EvapFlux_{(i,j,k)}\Delta t)\cdot V_{cell}</math>

'''Infiltration Flux'''

<math>\theta=(\theta\cdot V_{cell}-(UnsatK\cdot Area_{cell}\cdot(1-Imp) \Delta t)\cdot V_{cell}</math>

where:

:{|

|''<math>\theta</math>'' || is the water content (m 3/m 3)
|-
|''<math>V_{cell}</math>'' || is the cell volume (m 3)
|-
|''<math>Area_{cell}</math>'' || is the cell area (m 2)
|-
|''ComputeFace'' || is the computed face (dimensionless)
|-
|''<math>\Delta t</math>'' || is the time step (s)
|-
|''FluxX'' || is the flux in X direction (m 3/s)
|-
|''FluxV'' || is the flux in Y direction (m 3/s)
|-
|''FluxW'' || is the flux in W direction (m 3/s)
|-
|''TranspFlux'' || is the transpiration flux (m 3/s)
|-
|''EvapFlux'' || is the evaporation flux (m 3/s)
|-
|''UnsatK'' || is unsaturated conductivity (m /s)
|-
|''Imp'' || is the percentage of impermeable soil of the cell (dimensionless)
|}

The value of <math>\theta</math> so obtained is compared with one used in the volumes calculation and the iterative process stop when:

<math>(\theta^{'})-(\theta^{new})<\,\, Tolerance\hspace{3cm} (1.9) </math>

where:

:{|
|''<math>\theta^{'}</math>'' || is the water content of the previous iteration (m 3/m 3)
|}

If the equation (1.9) is not satisfy the temporal step is divided in half and the new value of <math>\theta</math> is used for solving the equation (1.7) for restarting the calculation process. The iteration process is stoped when the tolerance desired is reached.

[[Image:Figure06.jpg|thumb|center|300px|Figure 2: Time step reduction]]

==Other Features==
===How to Generate Info needed in Porous Media===
====SoilMap====
Model needs to know soil ID in each cell and layer to pick hydraulic properties from that type of soil. In Pedology soil includes more than one horizon, each with different soil properties. Here Soil is used has a unit of soil hidraulic properties, i.e., to define a soil with three horizons one has to create three SoilID (see below). This also means that if a watershed has at least one soil with three horizons one has to create three soil maps. Each soil map will be infact the map of each horizon of the soils. The grid cells with only one horizon will have the same SoilID in all maps, the grid cells with three horizons will have a different SoilID in each map.

Options:
* Constant value
* Soil Grid. One possible option is to associate with soil shape file. In this case can use MOHID GIS going to menu [Tools]->[Shape to Grid Data] and provide: i) the grid (model grid), ii) the soil shape file and iii) the corespondence between soil codes and soil ID defined in data file.

Soil ID must be defined in [[Module_FillMatrix|Module FillMatrix]] standards for each soil horizon defined (grid example):
<beginhorizon>
KLB : 1
KUB : 10

<beginproperty>
NAME : SoilID
DEFAULTVALUE : 1
INITIALIZATION_METHOD : ASCII_FILE
REMAIN_CONSTANT : 1
FILENAME : ..\..\GeneralData\PorousMedia\SoilID200m.dat
<endproperty>

<beginproperty>
<endproperty>
..
<endhorizon>

'''Remarks'''

All the soil ID's appearing in the soil grid(s) must be defined in the PorousMedia data file in terms of hydraulic properties:
<beginsoiltype>
ID : 1
THETA_S : 0.3859
THETA_R : 0.0476
N_FIT : 1.39
SAT_K : 3.5556e-6
ALPHA : 2.75
L_FIT : 0.50
<endsoiltype>

<beginsoiltype>
ID : 2
...
<endsoiltype>
...

====Soil Bottom====
The soil depth must be known by the model. This is computed by the model from terrain altitude (topography) and soil bottom altitude. As so, a soil bottom grid is needed.

Options:
* Grid File.
Soil depth (and soil bottom altitude, the effective grid needed) can be defined with a constant depth or estimated from slope [[HOW TO SoilBottom LINK]]. When the soil depth is estimated as a function of slope, soil depth will be smaller in ares with higher slope. In this areas only the surface layers of the soil will be considered (see Porous Media Geometry).

Define the grid just generated, in the porous media data file with:
BOTTOM_FILE : ..\..\GeneralData\PorousMedia\BottomLevel.dat

====Water Level====
Options:
*Grid File.
The water table altitude represents the initial altitude of the water table.
It is recommended to do a spin-up run to estabilize water level and then do a continuous simulation starting with the final water table achieved.
Use the following blocks with [[Module_FillMatrix|Module FillMatrix]] standards:
<beginwaterlevel>
NAME : waterlevel
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 0
FILENAME : ..\..\GeneralData\PorousMedia\WaterLevel0.50.dat
<endwaterlevel>

====Impermeability====
Impermeability values (0 - completely permeable, 1 - impermeable) must be provided.

Options:
* Constant Value.
* Grid File. One possible option is to associate with land use shape file. In this case can use MOHID GIS going to menu [Tools]->[Shape to Grid Data] and provide: i) the grid (model grid), ii) the land use shape file and iii) the corespondence between land use codes and Impermeability values.
Use the following blocks with [[Module_FillMatrix|Module FillMatrix]] standards:
<beginimpermeablefraction>
NAME : impermeablefraction
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 1
FILENAME : ..\..\GeneralData\PorousMedia\AreaImpermeavel.dat
<endimpermeablefraction>

==Outputs==

===Timeseries===

Theta - is water content of selected cell (vol water/ vol soil)

relative_water_content - content in selected cell. between zero and one (zero is residual water content and one is saturated water content)

VelW_[m/s] - vertical velocity in the bottom face of the selected cell

VelW_Corr_[m/s] - vertical velocity in the bottom face of the selected cell that may be corrected if oversaturation occurs. if no correction occurs is the same as previous.

InF_Vel_[m/s] - infiltration velocity (in the soil surface)

Head_[m] - Suction in selected cell

Conductivity_[m/s] - Conductivity in selected cell

level_water_table_[m] - water table altitude

water_table_depth_[m] - water table depth (from soil surface)

Hydro_Pressure_[m] - hydrostatic pressure in selected cell

Final_Head_[m] - Soil water charge in selected cell

[Check Mohid Land Heights and Levels to understand some of the outputs]

GW_flow_to_river_total_[m3/s] - Ground water flow to river if it is a river point

Surface_Evaporation_Flux_[m3/s] - evaporation flux (in the soil surface)

Transpiration_Flux_[m3/s] - transpiration flux in selected cell

[http://screencast.com/t/I8We0B4j2uJ Check Mohid Land Heights and Levels to understand the difference between level and height results]

==References ==
*Jury,W.A.,Gardner,W.R.,Gardner,W.H., 1991,Soil Physics
*Van Genuchten, M.T., A closed form equation for predicting the hydraulic conductivity of unsaturated soils
*Wu,J.,Zhang, R., Gui,S.,1999, Modelling soil water movement with water uptake by roots, Plant and soil 215: 7-17
*Marcel G.Schaap and Martinus Th. van Genuchten, A modified Maulem van Genuchten Formulation for Improved Description of Hydraulic Conductivity Near Saturation, 16 December 2005

==Data File ==

===Keywords===
Keywords read in the Data File

Keyword : Data Type Default !Comment

BOTTOM_FILE : char - !Path to Bottom Topography File
START_WITH_FIELD : logical 1 !Sets Theta initial Field Capacity
CONTINUOUS : logical 0 !Continues from previous run
STOP_ON_WRONG_DATE : logical 1 !Stops if previous run end is different from actual
!Start
OUTPUT_TIME : sec. sec. sec. - !Output Time
TIME_SERIE_LOCATION : char - !Path to File which defines Time Series
CONTINUOUS_OUTPUT_FILE : logical 1 !Writes "famous" iter.log
CONDUTIVITYFACE : integer 1 !Way to interpolate conducivity face
!1 - Average, 2 - Maximum, 3 - Minimum, 4 - Weigthed
HORIZONTAL_K_FACTOR : real 1.0 !Factor for Horizontal Conductivity = Kh / Kv
CUT_OFF_THETA_LOW : real 1e-6 !Disables calculation when Theta is near ThetaR
CUT_OFF_THETA_HIGH : real 1e-15 !Set Theta = ThetaS when Theta > ThetaS - CUT_OFF_THETA_HIGH
MIN_ITER : integer 2 !Number of iterations below which the DT is increased
MAX_ITER : integer 3 !Number of iterations above which the DT is decreased
LIMIT_ITER : integer 50 !Number of iterations of a time step (for restart)
THETA_TOLERANCE : real 0.001 !Converge Parameter
INCREASE_DT : real 1.25 !Increase of DT when iter < MIN_ITER
DECREASE_DT : real 0.70 !Decrease of DT when iter > MAX_ITER

<beginproperty>
NAME : Theta / waterlevel

see Module FillMatrix for more options

<endproperty>

===Sample===

Some keywords of the PorousMedia input file:
BOTTOM_FILE : ..\General Data\Other\PorousMedia\SoilBottom_2cells.dat
TIME_SERIE_LOCATION : ..\General Data\TimeSeries\TimeSeriesLocation3D_2m.dat

COMPUTE_SOIL_FIELD : 1

OUTPUT_TIME : 0 86400
CUT_OFF_THETA_HIGH : 1e-15
START_WITH_FIELD : 1

LIMIT_EVAP_WATER_VEL : 0
LIMIT_EVAP_HEAD : 0
HEAD_LIMIT : -100
THETA_HYDRO_COEF : 0.98

<beginsoiltype>
ID : 1
THETA_S : 0.43
THETA_R : 0.078
SAT_K : 2.888e-6
N_FIT : 1.56
ALPHA : 3.6
L_FIT : 0.50
<endsoiltype>

!----- Hydraulic Soil Properties
<beginhorizon>
KLB : 1
KUB : 10

<beginproperty>
NAME : SoilID
INITIALIZATION_METHOD : CONSTANT
DEFAULTVALUE : 1
<endproperty>

<beginproperty>
NAME : Theta
INITIALIZATION_METHOD : CONSTANT
DEFAULTVALUE : 0.30
<endproperty>
<endhorizon>

<beginwaterlevel>
NAME : waterlevel
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 1.
REMAIN_CONSTANT : 1
FILENAME : ..\General Data\Initial Conditions\InitialWaterLevel_2cells.dat
<endwaterlevel>

<beginimpermeablefraction>
NAME : impermeablefraction
INITIALIZATION_METHOD : ASCII_FILE
DEFAULTVALUE : 0
REMAIN_CONSTANT : 1
FILENAME : ..\General Data\Other\PorousMedia\InitialImpermeabilization_2cells.dat
<endimpermeablefraction>

[[Category:MOHID Land]]

Module Geometry

2010-09-30T17:05:26Z

Carina:

== Overview ==
Module Geometry handles the vertical discretization in MOHID. It was designed to divide the water column (in MOHID Water) or the soil compartment (in MOHID Land) in different vertical coordinates: Sigma, Cartesian, Lagrangian, Fixed Spacing, Harmonic, etc. A subdivision of the vertical domain into different sub-domains using different vertical coordinate systems is also possible.

== General options ==
*Minimum depth

== Vertical coordinate system ==

=== Sigma ===
[[Image:SZZ.gif|center|425px|thumb|Vertical sigma mesh]]

=== Cartesian ===
The Cartesian coordinate can be used with or without [[shaved cells]]. Cartesian builds layers from bottom to top with fixed depth of each layer with the possibility of colapsing layers from top to bottom. This hapens if for example water level reduces. The bottom layers are the last to colapse.

=== Fixspacing ===
The [[Fixed Spacing coordinates|Fixed Spacing]] coordinate allows the user to study flows close to the bottom.

=== Lagrangian ===
The [[Lagrangian coordinate]] moves the upper and lower faces with the vertical flow velocity.

=== Harmonic ===
The Harmonic coordinate works like the Cartesian coordinate, just that the horizontal faces close to the surface expand and collapse depending on the variation of the surface elevation. This coordinate was implemented in the geometry module to simulate reservoirs.

=== Fixsediment ===

=== SigmaTop ===

=== Cartesiantop ===
Cartesiantop is equal to cartesian but builds layers from top to bottom. This type of coodinates are used for [[Mohid_Land|Mohid Land]]. The top is the topography and the bottom is the non-porousmedia (rock). This means that in lower depth soils there will be less layers than in higher depth soils.

== Distances ==
[[Image:Mohid_distances.JPG|center|250px|thumb|MOHID syntax for distances]]
===Public routines===
; ModuleHorizontalGrid: GetHorizontalGrid(HorizontalGridID, XX_IE, YY_IE, XX_Z, YY_Z,XX_U, YY_U, XX_V, YY_V, XX_Cross, YY_Cross, DXX, DYY, DZX, DZY, DUX, DUY, DVX, DVY, XX, YY, XX2D_Z, YY2D_Z, STAT)
; ModuleGeometry: GetGeometryDistances(GeometryID, SZZ, DZZ, DWZ, DUZ, DVZ, DZI, DZE,ZCellCenter, ActualTime, STAT)

== Areas ==

== Volumes ==
[[Image:VolumeDeControlo.gif|center|250px|thumb|Single T-cell control volume]]

<htm>
<a href="http://content.screencast.com/users/GRiflet/folders/Jing/media/cf704a90-e3ba-4068-8843-8cae69f38bbc/Arakawa C grid.png"><img src="http://content.screencast.com/users/GRiflet/folders/Jing/media/cf704a90-e3ba-4068-8843-8cae69f38bbc/Arakawa C grid.png" width="544" height="191" border="0" /></a>
</htm>

== Bathymetry consistency diagnostic ==

Once the vertical discretization is imposed and the bathymetry is chosen, the bottom layer can yield stability problems when using shaved cells. You can have very thin bottom cell next to a very wide bottom cell. To diagnose the existence of such problematic cells, a [[geometry diagnostic]] tool was developed.

[[Category:Modules]]
[[Category:MOHID Base 2]]

== Input data file ==
===Keywords===
IMPERMEABILITY : 0/1 - !Consider impermeable cell faces
IMPER_COEF_U : real [1] !
IMPER_COEFX_U : real [0] !
IMPER_COEF_V : real [1] !
IMPER_COEFX_V : real [0] !

<begindomain>
ID : int - !Domain ID
TYPE : char - !Type of vertical coordinate of the domain
!Multiple options: FIXSPACING, SIGMA,
!LAGRANGIAN, CARTESIAN, HARMONIC, FIXSEDIMENT, CARTESIANTOP.
LAYERS : int - !Number of layers
EQUIDISTANT : real [0] !Equidistant layers spacing in meters
LAYERTHICKNESS : real vector - !If not equidistant specifies layers thickness
!starting from bottom layer (e.g. 50. 20. 10. 5.)
TOLERANCEDEPTH : real [0.05] !Just for SIGMA,ISOPYCNIC, LAGRANGIAN coordinates
TOTALTHICKNESS : real - !Total domain thickness
!(Just for FIXSPACING, FIXSEDIMENT, SOIL_TOPLAYER)
EMPTY_TOP_LAYERS : int [0] !Number of empty layers counting from top
DOMAINDEPTH : real
<enddomain>

Module Geometry

2010-09-30T17:03:14Z

Carina:

== Overview ==
Module Geometry handles the vertical discretization in MOHID. It was designed to divide the water column (in MOHID Water) or the soil compartment (in MOHID Land) in different vertical coordinates: Sigma, Cartesian, Lagrangian, Fixed Spacing, Harmonic, etc. A subdivision of the vertical domain into different sub-domains using different vertical coordinate systems is also possible.

== General options ==
*Minimum depth

== Vertical coordinate system ==

=== Sigma ===
[[Image:SZZ.gif|center|425px|thumb|Vertical sigma mesh]]

=== Cartesian ===
The Cartesian coordinate can be used with or without [[shaved cells]]. Cartesian builds layers from bottom to topwith fixed depth of each layer with the possibility of colapsing layers from top to bottom. This hapens if for example water level reduces. The bottom layers are the last to colapse.

=== Fixspacing ===
The [[Fixed Spacing coordinates|Fixed Spacing]] coordinate allows the user to study flows close to the bottom.

=== Lagrangian ===
The [[Lagrangian coordinate]] moves the upper and lower faces with the vertical flow velocity.

=== Harmonic ===
The Harmonic coordinate works like the Cartesian coordinate, just that the horizontal faces close to the surface expand and collapse depending on the variation of the surface elevation. This coordinate was implemented in the geometry module to simulate reservoirs.

=== Fixsediment ===

=== SigmaTop ===

=== Cartesiantop ===
Cartesiantop is equal to cartesian but builds layers from top to bottom. This type of coodinates are used for [[Mohid_Land|Mohid Land]]. The top is the topography and the bottom is the non-porousmedia (rock). This means that in lower depth soils there will be less layers than in higher depth soils.

== Distances ==
[[Image:Mohid_distances.JPG|center|250px|thumb|MOHID syntax for distances]]
===Public routines===
; ModuleHorizontalGrid: GetHorizontalGrid(HorizontalGridID, XX_IE, YY_IE, XX_Z, YY_Z,XX_U, YY_U, XX_V, YY_V, XX_Cross, YY_Cross, DXX, DYY, DZX, DZY, DUX, DUY, DVX, DVY, XX, YY, XX2D_Z, YY2D_Z, STAT)
; ModuleGeometry: GetGeometryDistances(GeometryID, SZZ, DZZ, DWZ, DUZ, DVZ, DZI, DZE,ZCellCenter, ActualTime, STAT)

== Areas ==

== Volumes ==
[[Image:VolumeDeControlo.gif|center|250px|thumb|Single T-cell control volume]]

<htm>
<a href="http://content.screencast.com/users/GRiflet/folders/Jing/media/cf704a90-e3ba-4068-8843-8cae69f38bbc/Arakawa C grid.png"><img src="http://content.screencast.com/users/GRiflet/folders/Jing/media/cf704a90-e3ba-4068-8843-8cae69f38bbc/Arakawa C grid.png" width="544" height="191" border="0" /></a>
</htm>

== Bathymetry consistency diagnostic ==

Once the vertical discretization is imposed and the bathymetry is chosen, the bottom layer can yield stability problems when using shaved cells. You can have very thin bottom cell next to a very wide bottom cell. To diagnose the existence of such problematic cells, a [[geometry diagnostic]] tool was developed.

[[Category:Modules]]
[[Category:MOHID Base 2]]

== Input data file ==
===Keywords===
IMPERMEABILITY : 0/1 - !Consider impermeable cell faces
IMPER_COEF_U : real [1] !
IMPER_COEFX_U : real [0] !
IMPER_COEF_V : real [1] !
IMPER_COEFX_V : real [0] !

<begindomain>
ID : int - !Domain ID
TYPE : char - !Type of vertical coordinate of the domain
!Multiple options: FIXSPACING, SIGMA,
!LAGRANGIAN, CARTESIAN, HARMONIC, FIXSEDIMENT, CARTESIANTOP.
LAYERS : int - !Number of layers
EQUIDISTANT : real [0] !Equidistant layers spacing in meters
LAYERTHICKNESS : real vector - !If not equidistant specifies layers thickness
!starting from bottom layer (e.g. 50. 20. 10. 5.)
TOLERANCEDEPTH : real [0.05] !Just for SIGMA,ISOPYCNIC, LAGRANGIAN coordinates
TOTALTHICKNESS : real - !Total domain thickness
!(Just for FIXSPACING, FIXSEDIMENT, SOIL_TOPLAYER)
EMPTY_TOP_LAYERS : int [0] !Number of empty layers counting from top
DOMAINDEPTH : real
<enddomain>

GenerateGrid

2010-02-26T19:12:34Z

Carina:

== Overview ==
GenerateGrid is a [[Mohid Support Tools|MOHID support tool]] used to create regular grids with a non-constant grid size. It reads an input data file named ''GridGenerator.dat''.

== Input data file (GridGenerator.dat)==
Below is a sample input data file.

OUTPUT_FILE : Sample.grd !Output file name
ORIGIN_X : -9 !Coordinate X of the lower left corner
ORIGIN_Y : 39 !Coordinate Y of the lower left corner
GRID_ANGLE : 0 !Rotation angle of the grid
COORD_TIP : 4 !Coordinate type

<begin_grid_xx>
NUMBER_OF_SEGMENTS : 1 !Number of segments in the XX direction
NODES : 0. 2 !Coordinates starting in 0. of the XX nodes
RESOLUTION : 0.1 0.1 !Resolution at each of the defined XX nodes
<end_grid_xx>

<begin_grid_yy>
NUMBER_OF_SEGMENTS : 2 !Number of segments in the YY direction
NODES : 0. 1.5 3 !Coordinates starting in 0. of YY nodes
RESOLUTION : 0.1 0.03 0.15 !Resolution at each of the defined YY nodes
<end_grid_yy>

With these options a grid file is created, with origin in -9ºW and 39ºN (the approach is independent of the type of coordinates), with a constant spacing in the XX direction (resolution of 0.1º) and a non-constant spacing in the YY direction. In the YY direction, 2 segments are defined. The first starts at 0º from the origin and ends at 1.5º from the origin. The resolution of this segment starts at 0.1º and progressively decreases to 0.03º at the end of the segment. The second segment starts at the end of the first and ends at 3º from the origin. This means that the entire grid is 3º long in the YY direction. The resolution of the second segment starts at 0.03º and progressively decreases to 0.15º at the end of the segment.

==Fast start==
How to create a grid. Go to MOHID Gis. Create NewData items choosing simple grid, then adjust coordinates.

[[Category:Tools]]