The input maps
All input and output maps in LISEM are in the format of the PCRaster GIS. This raster GIS is produced by PCRaster environmental software, and can be obtained through their website, where you will find a complete list of functions and examples.
LISEM needs a minimum of 24 maps depending on the input options selected in the interface (e.g. enabling the simulation of ditches requires additional maps, and some infiltration option require more maps than others). If needed all maps can be derived from 4 base maps:
- digital elevation model
- land use
- soil type
- impermeable areas such as tarred roads (this is needed because the road can be narrower than the gridcell size and LISEM needs to know what type of soil or land use exists adjacent to the road).
All input data can be derived form these four base maps, but off course other sources can be used to make the input maps (geostatistical interpolation, remote sensing etc.).
The spatial input data for LISEM can be structured as follows:
- Catchment maps
- Vegetation maps
- Soil surface maps
- Infiltration related maps
- Erosion/deposition related maps
- Channels
In this part all input maps are listed with the PCRaster commands to make them. As an example a catchment of INRA in Northern France will be used (the “Tavaux” catchment). Examples of PCRaster commands to create the maps are shown with each map. Since PCRaster is also capable of running macros and even spatial models, all commands given here are combined in a script, so that the complete database can be constructed in one go. The script can be found at the download page.
| NOTE: the map names given here are NOT compulsory but can be specified individually by the user in the interface. However, to enable compatibility with the older versions of LISEM the map names used here are the default names. Also the extension “.map” is NOT compulsory, just convenient. |
1. Catchment maps
The actual digital elevation model is not used in LISEM, only maps derived from it. This has the advantage that a number of important maps such as the surface drainage network (called Local Drain Direction or LDD in PCRaster) and Gradient are not automatically derived with GIS operations, but can be supplied by the user.
| Map name | Contents | Type | Unit | Range |
| LDD | Local drain direction | ldd | 1-9 (see below) | |
| AREA | Catchment boundaries | boolean | 1 | |
| ID | Area covered by raingauges | nominal | 1- n (= # of raingauges) | |
| GRAD | Slope gradient (sine of slope angle) | scalar | – | must be > 0 and <= 1 |
| OUTLET | Location of outlet and suboutlets | nominal | 0-3 |
LDD: Local Drain Direction, this map gives for each cell the direction of the surface runoff (the number corrspond to a keypad on your keyboard: 1=drainage to lower left, 6 is drainage to the right etc.) the number 5 is reserved for the catchment outlet. Note that this map cannot contain local depressions, the entire catchment has to drain towards the outlet. PCRaster provides commands that can automatically generate such a map, with removal of the local depressions (‘pits’) based on constraints given by the user. In this case the LDD is based on the steepest slope. The command is:
pcrcalc –lddin ldd.map=lddcreate(dem.map,1e10,1e10,1e10,1e10)
The 4 large values are thresholds to remove the pits and create a continuous LDD.
It is also ppossible to generate a valid LDD map based on the tillage direction. A PCRaster macro has been constructed by Takken et al. (2000) that combines major topographic concenration lines with tillage direction per field, including ‘headlands’ and ‘dead furrows’.
| NOTE: errors with LISEM related to the kinematic wave where MVs are reported (missing values) are often because the LDD does not have the same dimensions as the other maps. |
AREA: all maps are checked against this map for the number and location of non-missing value cells. If maps do not have the same size and number of valid cells, an error occurs.
pcrcalc area.map=catchment(ldd.map,pit(ldd.map))
GRAD: this is the gradient in the direction of the LDD. PLEASE NOTE: this must be the SINE of the slope. Because the slope can occur in the denominator of various algorithms, it cannot be 0 (zero). Use the slope command to create the slope from the DEM:
pcrcalc grad.map=slope(dem.map))
or make a gradient map with the tangent of the slope in the direction of the flow (it has to be covered with the slope derived from the DEM to account for edge cells):
pcrcalc grad.map=(dem.map-downstream(ldd.map,dem.map))/downstreamdist(ldd.map)
pcrcalc grad.map=cover(grad.map,slope(dem.map))
convert this map to a sine and make sure it is larger than 0::
pcrcalc grad.map=sin(atan(max(0.001, grad.map)))
ID: this map can be made in various ways, using Tiessen Polygons or zones based on relief. The numbers in this map must corrspond to the number of raingauges in the rainfall file. An example is to create a map according to Thiessen polygons from a map with pxiels that have the numbers of the rain gauges:
pcrcalc id.map=nominal(spreadzone(gauge.map,0,1))
Alternatively the relief can be taken more into account by taking the slope as a weighing factor:
pcrcalc id.map=nominal(spreadzone(gauge.map,0,grad.map))
or by using the mean elevation line between two gauges as a separation line (gauges have id numbers 1, 2 and 3 here, dem is a map):
pcrcalc a1 = maptotal(if(p1 eq 1,dem,0))
pcrcalc a2 = maptotal(if(p2 eq 2,dem,0))
pcrcalc a3 = maptotal(if(p3 eq 3,dem,0))
pcrcalc a4 = (a1+a2)/2
pcrcalc a5 = (a2+a3)/2
pcrcalc id.map = nominal(if(dem lt a4,1,if(dem lt a5 and dem ge a4, 2, 3)))
NOTE: that the gauges must be situated inside the ldd.map area otherwise it will not be included in the spread operation.
OUTLET: a map with values 0 (background) and 1 which is the main outlet. On this cells a hydrograph and sedigraph are stored as ASCII files. Two additional cells can be given (numbers 2 and 3 compulsory) to store hydrographs and sedigraphs for additional locations (e.g. sub catchments). The main outlet is found by:
pcrcalc outlet.map=pit(ldd.map)
2. Vegetation maps
Vegetation and/or crops are mainly used by the interception and splash erosion processes in LISEM. The maps are usually made by reclassifying the field boundary map. vegetation influences also the soil cohesion for which a separate map must be defined (see below).
Note that LISEM has never been tested for forest areas and should be used with care in these circumstances.
| Map name | Contents | Type | Unit | Range |
| LAI | Leaf area index | scalar | – | 0-12 |
| PER | Fraction of soil covered by vegetation | scalar | – | 0-1 |
| CH | Vegetation height | scalar | m | 0-30 |
LAI: leaf area index, used to calculate the storage capacity of the canopy. The value represents the LAI of the vegetated part of the gridcell. For instance if a grid cell has a single tree with a LAI of 6 in it, the LAI of that cell is 6.
PER: fraction of soil cover by canopy, but also litter or other types of soil cover (except for stones).
CH: crop height, used to calculate throughfall kinetic energy in the splash erosion part.
These maps can be made in similar operations if a table is available that has as the first column the id numbers of the landuse map-units, and in the other columns values for the LAI, cover fraction and vegetation height (the option “–matrixtable” is needed because of the table layout conventions in PCRaster):
pcrcalc –matrixtable lai.map = lookupscalar(landuse.map, 1, veg.tbl)
pcrcalc –matrixtable per.map = lookupscalar(landuse.map, 2, veg.tbl)
pcrcalc –matrixtable ch.map = lookupscalar(landuse.map, 3, veg.tbl)
The table may look like this (four landuse types with different IDs not necessarily sequential):
======begin of file======
| 1 | 0 | 0.1 | 0 |
| 2 | 1.2 | 0.3 | 0.3 |
| 4 | 1.23 | 0.35 | 0.4 |
| 16 | 6.1 | 0.95 | 1.8 |
======end of file======
NOTE: LISEM assumes no logic in this table! It will not check if a LAI of 6 is combined with a cover fraction of 0.
3. Soil surface maps
This is a mixed catagory with maps that determine several processes: infiltration, surface storage, overland flow (veloccity). In an agricultural catchment these maps are based on the landuse mostly, or on a combination of landuse and soils. Note that in a crop rotation cycle these values will change according to for instance crusting stage, soil cover and tillage operation.
It also includes different types of surface for which the processes are calculated separately.
| Map name | Contents | Type | Unit | Range |
| N | Manning’s n | scalar | – | 0.001 – 0.5 |
| RR | Random Roughness | scalar | cm | 0.05 – 20 |
| STONEFRC | Fraction covered with stones | scalar | – | 0 – 1 |
| CRUSTFRC | Fraction cover with a crust | scalar | – | 0 – 1 |
| ROADWIDT | width of impermeable roads | scalar | m | 0-cellsize |
N: Surface resistance to flow, expressed as Manning’s n.
RR: ‘random roughness’ used here as the standard deviation of the micro relief heights (do not use transformations to obtain a more ‘normal’ distribution).
STONEFRC: fraction of stones for which no splash erosion is calculated (used only in this process).
CRUSTFRC: fraction of the soil crusted. On gridcells with non-zero values the infiltration is calculated according to soil physical properties supplied by the user, depending on the infiltration option chosen.
ROADWIDT: width of imermeable roads for which no infiltration is calculated (is used in all infiltration modules).
pcrcalc –matrixtable n.map = lookupscalar(landuse.map, 1, surface.tbl)
pcrcalc –matrixtable rr.map = lookupscalar(landuse.map, 2, surface.tbl)
pcrcalc –matrixtable stonefrc.map = lookupscalar(landuse.map, 3, surface.tbl)
pcrcalc –matrixtable crustfrc.map = lookupscalar(landuse.map, 4, surface.tbl)
pcrcalc roadwidt.map = scalar(if (roads.map neq 0, *width*, 0))
where “surface.tbl” is a table such as used with the vegetation related maps, and *width* is a value for the road width in m.
4. Infiltration related maps
There are currently 6 infiltration options in LISEM, based on 3 different models. Depending on the choice of the user in the interface, a different set of input maps is needed.
- no infiltration: no maps needed
- SWATRE (including crusts, wheel tracks and grass strips)
- Holtan/Overton
- Green & Ampt: 1 and 2 layers and additional surface types
- Morell and Seytoux
4.1 No infiltration
The entire catchment is assumed impermeable. Can be used for testing.
4.2 Swatre
To operate the SWATRE infiltration module , the user has to set up a series of maps and tables with soil physical properties that describe the profile (max 20 layers) of each land unit. These can be e.g. soil type units or fields or parts of fields.
| Map name | Contents | Type | Unit | Range |
| PROFILE.INP | ASCII lookup file describing soil profiles of map units | |||
| name.TBL | ASCII tables with soil physical data | |||
| PROFILE.MAP | map with profile ID numbers | nominal | – | >= 1 |
| PROFCRST.MAP | map with ID number of crust profiles | nominal | – | >= 1 |
| PROFWHL.MAP | ID number of wheel track profile | nominal | – | >= 1 |
| PROFGRAS.MAP | ID number of grass strip profile | nominal | – | >= 1 |
| HEADOUT.MAP | locations of matric potential output | nominal | – | 0, 1-4 |
| INITHEAD.001 …INITHEAD.0nn | positive initial matric potential of each soil layer | scalar | cm | 0..100000 |
PROFILE.INP: when the SWATRE infiltration model is used LISEM needs a 3D structure of the soil. This structure is described in the text file PROFILE.INP. In this file the number of nodes are given for the solution of the Richards equation (using a finite difference approximation) and their depths. These are independent of the actual soil horizons. NOTE that the node depths are used in LISEM as boundaries between simulation layers (not centres). Below the nodes the soil profiles are given (see syntax below) with ‘pointers’ to the names of the tables describing the soil physical relations (pF curve and hydraulic conductivity curve) of a particular soil type.
NOTE that in this file ALL information is stored, including that of the units for crusts, wheel tracks and grass strips. So if these features are simulated the file PROFILE.INP has to be adjusted to reflect their presence.
The file PROFILE.INP has the following layout:
| ===================== | Begin of file ============== |
| 14 | number of soil layers |
| 2.5 | 1:depth of boundary between first and second layer (cm) |
| 5 | 2:depth of boundary between second and third layer |
| 10 | 3: etc. |
| 15 | |
| 20 | |
| 25 | |
| 30 | |
| 40 | |
| 50 | |
| 60 | |
| 70 | |
| 80 | |
| 90 | |
| 100 | 14: depth of final soil layer (cm) |
| 1 | description of map unit 1. Used in Swatre if PROFILE.MAP, PROFWLTR.MAP or PROFCRST.MAP has pixel value of 1 |
| ASCRP.TBL | the name of the soil physics table to be used |
| 15 | the depth down to which the above table is used (cm) |
| SUB_REST.TBL | use this table for depth 15-30 |
| 30 | |
| C_ALL.TBL | |
| 100 | |
| 2 | description of map unit 2. Used in Swatre if PROFILE.MAP, PROFWLTR.MAP or PROFCRST.MAP has pixel value of 2 |
| ASCRP_A.TBL | etc. |
| 15 | |
| SUB_REST.TBL | |
| 30 | |
| C_ALL.TBL | |
| 100 | |
| 3 | etc |
| ASCRP_C.TBL | |
| 15 | |
| SUB_REST.TBL | |
| 30 | |
| C_ALL.TBL | |
| 100 | |
| ==================== | End of file ================== |
name.TBL: in the soil above three horizons are defined, each of them with seperate soil phyisical characteristics, referring to the specific tables. The horizons used here are: 0-15 cm, 15-30 cm and 30-100 cm. The tables used for these layers contain 3 columns: soil moisture content (theta, fraction), pressure head (h in cm), and unsaturated conductivity (k in cm/day). An example:
0.01 -4.86E+24 1.95E-21
0.02 -1.04E+07 1.64E-11
0.03 -9.02E+05 2.12E-09
0.04 -2.15E+05 3.63E-08
…..etc.
0.31 -3.19E+01 6.74E-01
0.32 -2.33E+01 9.89E-01
0.33 -1.54E+01 1.54E+00
0.34 0.00E+00 5.69E+01
PROFILE.MAP: a map based on soils, fields or a combination, with ID numbers that are used to look up the corresponding soil physics tables in PROFILE.INP.
PROFCRST.MAP: a map with ID numbers on those places where there is a crust, to look up the soil profiles for crusted soils in the file PROFILE.INP. The areas that are crusted are defined in CRUSTFRC.MAP. There can be more crust profiles if there are for instance more crop types with a different developing surface. Simply include more crust profile descriptions in PROFILE.INP.
PROFGRAS.MAP: a map with ID numbers corresponding to the soil profile for grass strips. The area where strips are present are given in GRASSWID.MAP.
PROFWHL.MAP: a map with ID numbers corresponding to the soil profile for wheel tracks. The area where tracks are present are given in WHEELWID.MAP. Note that there are more maps to define wheeltracks.
HEADOUT.MAP: map with up to four cells numbered 1 to 4 (de rest having value 0). When the option matric output is checked in the infiltration options of the input menu, an ASCII table with matric potential for each layer at each timestep is created.
INITHEAD.001 – INITHEAD.020: initial matric potential for each layer (if for instance you have defined 6 layers in PROFILE.INP, there are 6 maps 001 to 006). The values are positive matric potentials in cm.
4.3 Holtan/Overton
The Holtan infiltration model is a one layer empirical model. The maps need no further description. Note that the extra infiltration toptions (wheel tracks, crusts and grass strips) are not implemented for the Holtan model.
| Map name | Contents | Type | Unit | Range |
| FC | infiltration rate saturation | scalar | mm/h | 0-100 |
| A | difference initial and saturation rate | scalar | mm/h | 0-100 |
| TP | soil porosity | scalar | cm3/cm3 | 0.01-0.60 |
| DF | depth infiltration control zone | scalar | mm | 0.001-300 |
| P | coefficient infiltration equation | scalar | – | 0.4-0.8 |
| ASM | initial soil moisture content | scalar | – | 0.20-1 |
| FP | soil moisture at field capacity | scalar | – | 0.40-0.90 |
4.4 Green & Ampt, 1 and 2 layers and additional surface types
G&A demands the hydraulic properties of the first layer. NOTE that the soil depth is only of concequence if the option “impermeable sub soil” is checked, meaning that the soil stvorage is limited to the product of the moisture deficit and the depth: (thetas1-thetai1)*soildep1. If this option is not checked the soil is assumed of unlimited depth and at moisture content thetai1 for its entire depth.
| Map name | Contents | Type | Unit | Range |
| KSAT1 | Saturated hydraulic conductivity | scalar | mm/hr | 0-1000 |
| THETAS1 | Saturated volumetric soil moisture content | scalar | – | 0-1 |
| THETAI1 | Initial volumetric soil moisture content | scalar | – | 0-1 |
| PSI1 | Soil water tension at the wetting front | scalar | cm | 0-1000 |
| SOILDEP1 | Soil depth | scalar | mm | 0-1000 |
G&A1: Additional surfaces
In addition separate ksat maps have to be defined for crusted soils, wheel tracks or grass strips. These maps are required when anuy (or all) of the extra options in the infiltration options menu are checked.
| Map name | Contents | Type | Unit | Range |
| KSATWT | wheeltrack Ksat, used with infiltration option “include wheel tracks” | scalar | mm/hr | 0-1000 |
| KSATGRAS | grass strip Ksat, used with infiltration option “include grass strips” | scalar | mm/hr | 0-1000 |
| KSATCRST | crust Ksat, used with infiltration option “include crusts” | scalar | mm/hr | 0-1000 |
NOTE: in addition the user has to define the maps WHEELWID.MAP, GRASSWID.MAP and/or CRUSTFRC.MAP.
These options are only used for the first layer of the Green and Ampt module. If a 2 layer Green and Ampt is used, the second layer will be identical for all option.
G&A2: Two layer Green & Ampt
Requires all maps of Green and Ampt 1 layer plus an identical set for 2nd layer. The depth of the second layer (soildep2) is only of concequence when the “impermeable sub-soil” option is checked (see remark above).
| Map name | Contents | Type | Unit | Range |
| KSAT2 | see KSAT1 | scalar | mm/hr | see KSAT1 |
| THETAS2 | see THETAS1 | scalar | – | see THETAS1 |
| THETAI2 | see THETAI1 | scalar | – | see THETAI1 |
| PSI2 | see PSI1 | scalar | cm | see PSI1 |
| SOILDEP2 | see SOILDEP1 | scalar | mm | see SOILDEP1 |
5. Erosion/deposition related maps
Erosion and deposition in LISEM is based on the deficit transport capacity (capacity surplus equals deposition and transport deficit equals erosion). The particles in flow have a sttling velocity that is estimated by the settling velocity of the median fraction (D50 in mu). To calculate more sediment classes the LISEM-MULTICLASS SEDIMENT version has to be used.
| Map name | Contents | Type | Unit | Range |
| AGGRSTAB | Aggregate stability | scalar | – | 0.00001-200; -1 |
| COH | Cohesion of bare soil | scalar | kPa | COH+COHADD >= 0.196 |
| COHADD | Additional cohesion by roots | scalar | kPa | COH+COHADD >= 0.196 |
| D50 | D50 value of the soil | scalar | µm | (advise:) 25-300 |
AGGRSTAB: aggregate stability is used for the splash erosion. It is the median number of the drops needed to decrease the size of the aggregates by half on a sieve (LOWE test).
COH: Cohesion of the soil (kPa) measured with a torvane. The cohesion value is used in an empirical formula to decrease the transport deficit with a factor between 0 (no erosion) and 1 (chesionless soil).
COHADD: Additional cohesion (kPa) to be added by the user to simulate the effect of plant roots on the soil strength. Note that this is a calibration value as the Torvane test is not suited to measure the effect of roots.
D50: median of the texture of the soil (µm) used to simulate the settling velocity. This value determines strongly the depostion.
WARNING: if water with sediment with e.g. D50=25 flows over a gridcell with a D50=150, all sediment will deposit, even if the sediment in flow has a D50 that is less. To avoid these kind of effects use a homogeneous D50 map or better use the version LISEM MULTICLASS SEDIMENT.
6. Channels
Channels in LISEM are meant to be man-made ditches with a width much smaller than the grid cell (typically a ditch of less than 1.5 m wide in a grid cell of 10×10 m). The channel network MUST be a continuous network connected to the outlet. Water that enters the channel in LISEM does not come out again (the channels do not overflow) and if the channels are not connected to the outlet they are effectively sinks. LISEM does a separate kinematic wave for the channel.
Simulations can also be done considering a concentrated flow line, such as a thalweg, to be a channel. The PCRaster operation accuflux can be used to determine a channel based on the ldd:
pcrcalc chanmask.map = scalar(if(accuflux(ldd,cellarea()) lt 10000, 1))
This command creates a mask of a channel (having value 1 as the channel and MV outside) based on the contributing area: all cells that have 10000 m2 of drainage area upstream are considered a channel. The mask itself is not used in LISEM but it can be used to create the channel input maps:
CAUTION: The channel form (hydraulic radius) and resistance (manning’s n) determine to a large extent the shape of the hydrograph. Channels are considered impermeable.
| Map name | Contents | Type | Units | Range |
| LDDCHAN | Local drain direction of channel network | ldd | 1-9 | |
| CHANGRAD | Channel gradient | scalar | – | 0.0001-10 6) |
| CHANMAN | Manning’s n for the channel | scalar | – | 0.001-0.6 |
| CHANCOH | Cohesion of the channel bed | scalar | kPa | > 0.196 |
| CHANWIDT | Width of channel | scalar | m | 0-cellwidth |
| CHANSIDE | Channel cross section shape | scalar | – | 0-10 |
LDDCHAN: the LDD map based on the channel mask alone. This map can opnly have one pit that has to be the same as the pit (outlet) of the LDD. A way to create this map is:
pcrcalc lddchan.map=lddcreate(dem*chanmask.map,1e10,1e10,1e10,1e10)
CHANGRAD: the gradient of the channel bed based on the channel mask alone. A way to create this map is:
pcrcalc changrad.map=slope(dem*chanmask.map)
CHANMAN: the Manning’s n (resistance to flow of the channel). Usually this is not known and a literature value (ditches with vegetation) is assumed. If the channel is used to simulate a thalweg, the manning’s n of the land use could be used (n.map, see above):
pcrcalc chanman.map=n.map*chanmask.map
CHANCOH: the cohesion (kPa) of the channel bed, resistance to flow erosion. Usually this is not known and a literature value (ditches with vegetation) is assumed, or of a concrete channel is simulated, a very high cohesion could be assumed. If the channel is used to simulate a thalweg, the cohesion of the land use could be used (coh.map and cohadd.map, see above):
pcrcalc chancoh.map=(coh.map+cohadd.map)*chanmask.map
CHANWIDT: the channel bed width in m. Because the channel shape is not necessarily rectangular, the top width at water level, and the hydraulic radius, is calculated from the side angle and the water height in the channel. This value should be smaller than the grid cell width (not more than 90% of the cell width). The width is ussually a field measured value (_width_):
pcrcalc chanwidt.map=if (chanmask.map ne 0, _width_)
CHANSIDE: tangent of the angle between the channel side and the vertical: 0 is a rectangular channel, 1 is a side angle of 45 degrees, 10 is a very wide channel (84 degrees). Note that the width of the top of the water level depends on this side angle, if the water rises the stream becomes wider, but it cannot become wider than the gridcell.
pcrcalc chanside.map=if (chanmask.map ne 0, 0)