Source code for imod.flow.cap
import os
import pathlib
import jinja2
from imod.flow.pkgbase import Package
[docs]
class MetaSwap(Package):
"""
The MetaSWAP package (CAP), provides the input to be converted to a
MetaSWAP model, which is an external model code used to simulate the
unsaturated zone.
Note that only two-dimensional DataArrays with dimensions ``("y", "x")``
should be supplied to this package. In the current implementation
time-related files are provided as external files ("extra files"). Similar
to the iMODFLOW implementation of the projectfile. For now these need to be
provided as a path.
MetaSWAP is developed by Alterra, Wageningen as part of the SIMGRO model
code. The SIMGRO framework is intended for regions with an undulating
topography and unconsolidated sediments in the (shallow) subsoil. Both
shallow and deep groundwater levels can be modelled by MetaSWAP. This
model is based on a simplification of ‘straight Richards’, meaning that no
special processes like hysteresis, preferential flow and bypass flow are
modelled. Snow is not modelled, and neither the influence of frost on the
soil water conductivity. A perched watertable can be present in the SVAT
column model, but interflow is not modelled. Processes that are typical for
steep slopes are not included. The code contains several parameterized
water management schemes, including irrigation and water level management.
References:
* Van Walsum, P. E. V., 2017a. SIMGRO V7.3.3.2, Input and Output reference
manual. Tech. Rep. Alterra-Report 913.3, Alterra, Wageningen. 98 pp.
* Van Walsum, P. E. V., 2017b. SIMGRO V7.3.3.2, Users Guide. Tech. Rep.
Alterra-Report 913.2, Alterra, Wageningen. 111 pp.
* Van Walsum, P. E. V. and P. Groenendijk, 2008. "Quasi Steady-State
Simulation on of the Unsaturated Zone in Groundwater Modeling of Lowland
Regions." Vadose Zone Journal 7: 769-778.
* Van Walsum, P. E. V., A. A. Veldhuizen and P. Groenendijk, 2016. SIMGRO
V7.2.27, Theory and model implementation. Tech. Rep. Alterra-Report 913.1,
Alterra, Wageningen. 93 pp 491.
Parameters
----------
boundary : int or xr.DataArray of ints
2D boundary, used to specify active MetaSWAP elements, similar to
ibound in the Boundary package
landuse : int or xr.DataArray of ints
Landuse codes, referred to in the lookup table file luse_mswp.inp
rootzone_thickness : float or xr.DataArray of floats
Rootzone thickness in cm (min. value is 10 centimeter).
soil_physical_unit : int or xr.DataArray of ints
Soil Physical Unit, referred to in the lookup table file fact_mswp.inp.
meteo_station_number : float or xr.DataArray of ints
Meteo station number, referred to by mete_mswp.inp.
surface_elevation : float or xr.DataArray of floats
Surface Elevation (m+MSL)
sprinkling_type : int or xr.DataArray of ints
Sprinkling type ("Artificial Recharge Type" in iMOD manual):
* 0 = no occurrence
* 1 = from groundwater
* 2 = from surface water
sprinkling_layer : int or xr.DataArray of ints
Number of modellayer from which water is extracted ("Artificial
Recharge Location" in iMOD manual)
sprinkling_capacity : float or xr.DataArray of floats
Sprinkling capacity (mm/d) sets the maximum amount extracted for
sprinkling ("Artificial Recharge Capacity" in iMOD manual)
wetted_area : float or xr.DataArray of floats
Total area (m2) occupied by surface water elements. Values will be
truncated by maximum cellsize.
urban_area : float or xr.DataArray of floats
Total area (m2) occupied by urban area. Values will be truncated by
maximum cellsize.
ponding_depth_urban : float or xr.DataArray of floats
Ponding Depth Urban Area (m), specifying the acceptable depth of the
ponding of water on the surface in the urban area before surface runoff
occurs.
ponding_depth_rural : float or xr.DataArray of floats
Ponding Depth Rural Area (m), specifying the acceptable depth of the
ponding of water on the surface in the rural area before surface runoff
occurs.
runoff_resistance_urban : float or xr.DataArray of floats
Runoff Resistance Urban Area (day), specifying the resistance surface
flow encounters in the urban area. The minimum value is equal to the
model time period.
runoff_resistance_rural : float or xr.DataArray of floats
Runoff Resistance Rural Area (day), specifying the resistance surface
flow encounters in the rural area. The minimum value is equal to the
model time period.
runon_resistance_urban : float or xr.DataArray of floats
Runon Resistance Urban Area (day), specifying the resistance surface
flow encounters to a model cell from an adjacent cell in the urban
area. The minimum value is equal to the model time period.
runon_resistance_rural : float or xr.DataArray of floats
Runon Resistance Rural Area (day), specifying the resistance surface
flow encounters to a model cell from an adjacent cell in the rural
area. The minimum value is equal to the model time period.
infiltration_capacity_urban : float or xr.DataArray of floats
the infiltration capacity (m/d) of the soil surface in the urban area.
The range is 0-1000 m/d. The NoDataValue -9999 indicates unlimited
infiltration is possible.
infiltration_capacity_rural : float or xr.DataArray of floats
the infiltration capacity (m/d) of the soil surface in the urban area.
The range is 0-1000 m/d. The NoDataValue -9999 indicates unlimited
infiltration is possible.
perched_water_table : float or xr.DataArray of floats
Depth of the perched water table level (m)
soil_moisture_factor : float
Soil Moisture Factor to adjust the soil moisture coefficient. This
factor may be used during calibration. Default value is 1.0.
conductivity_factor : float
Conductivity Factor to adjust the vertical conductivity. This factor
may be used during calibration. Default value is 1.0.
lookup_and_forcing_files : list of pathlib.Path or str
List of paths to files required by MetaSWAP. This a list of
lookup tables and meteorological information that is required by
MetaSwap. Note that MetaSwap looks for files with a specific name, so
calling "luse_svat.inp" something else will result in errors. To view
the files required, you can call: ``print(MetaSwap()._required_extra)``
"""
_pkg_id = "cap"
_variable_order = [
"boundary",
"landuse",
"rootzone_thickness",
"soil_physical_unit",
"meteo_station_number",
"surface_elevation",
"sprinkling_type",
"sprinkling_layer",
"sprinkling_capacity",
"wetted_area",
"urban_area",
"ponding_depth_urban",
"ponding_depth_rural",
"runoff_resistance_urban",
"runoff_resistance_rural",
"runon_resistance_urban",
"runon_resistance_rural",
"infiltration_capacity_urban",
"infiltration_capacity_rural",
"perched_water_table",
"soil_moisture_factor",
"conductivity_factor",
]
_template_projectfile = jinja2.Template(
"0001, ({{pkg_id}}), 1, {{name}}, {{variable_order}}\n"
'{{"{:03d}".format(variable_order|length)}}, {{"{:03d}".format(n_entry)}}\n'
"{%- for variable in variable_order%}\n" # Preserve variable order
"{%- for layer, value in package_data[variable].items()%}\n"
"{%- if value is string %}\n"
# If string then assume path:
# 1 indicates the layer is activated
# 2 indicates the second element of the final two elements should be read
# 1.000 is the multiplication factor
# 0.000 is the addition factor
# -9999 indicates there is no data, following iMOD usual practice
'1, 2, {{"{:03d}".format(layer)}}, 1.000, 0.000, -9999., {{value}}\n'
"{%- else %}\n"
# Else assume a constant value is provided
'1, 1, {{"{:03d}".format(layer)}}, 1.000, 0.000, {{value}}, ""\n'
"{%- endif %}\n"
"{%- endfor %}\n"
"{%- endfor %}\n"
# Section for EXTRA FILES comes below
'{{"{:03d}".format(lookup_and_forcing_files|length)}},extra files\n'
"{%- for file in lookup_and_forcing_files %}\n"
"{{file}}\n"
"{%- endfor %}\n"
)
# TODO: Check which of these actually are required.
_required_extra = [
"fact_svat.inp",
"luse_svat.inp",
"mete_grid.inp",
"para_sim.inp",
"tiop_sim.inp",
"init_svat.inp",
"comp_post.inp",
"sel_key_svat_per.inp",
]
[docs]
def __init__(
self,
boundary,
landuse,
rootzone_thickness,
soil_physical_unit,
meteo_station_number,
surface_elevation,
sprinkling_type,
sprinkling_layer,
sprinkling_capacity,
wetted_area,
urban_area,
ponding_depth_urban,
ponding_depth_rural,
runoff_resistance_urban,
runoff_resistance_rural,
runon_resistance_urban,
runon_resistance_rural,
infiltration_capacity_urban,
infiltration_capacity_rural,
perched_water_table,
lookup_and_forcing_files,
soil_moisture_factor=1.0,
conductivity_factor=1.0,
):
super().__init__()
self.dataset["boundary"] = boundary
self.dataset["landuse"] = landuse
self.dataset["rootzone_thickness"] = rootzone_thickness
self.dataset["soil_physical_unit"] = soil_physical_unit
self.dataset["meteo_station_number"] = meteo_station_number
self.dataset["surface_elevation"] = surface_elevation
self.dataset["sprinkling_type"] = sprinkling_type
self.dataset["sprinkling_layer"] = sprinkling_layer
self.dataset["sprinkling_capacity"] = sprinkling_capacity
self.dataset["wetted_area"] = wetted_area
self.dataset["urban_area"] = urban_area
self.dataset["ponding_depth_urban"] = ponding_depth_urban
self.dataset["ponding_depth_rural"] = ponding_depth_rural
self.dataset["runoff_resistance_urban"] = runoff_resistance_urban
self.dataset["runoff_resistance_rural"] = runoff_resistance_rural
self.dataset["runon_resistance_urban"] = runon_resistance_urban
self.dataset["runon_resistance_rural"] = runon_resistance_rural
self.dataset["infiltration_capacity_urban"] = infiltration_capacity_urban
self.dataset["infiltration_capacity_rural"] = infiltration_capacity_rural
self.dataset["perched_water_table"] = perched_water_table
self.dataset["soil_moisture_factor"] = soil_moisture_factor
self.dataset["conductivity_factor"] = conductivity_factor
self.lookup_and_forcing_files = lookup_and_forcing_files
def _force_absolute_path(self, f):
"""Force absolute path, because projectfile cannot handle relative paths"""
return str(pathlib.Path(f).resolve())
def _render_projectfile(self, **kwargs):
"""
Returns
-------
rendered : str
The rendered projfectfile part,
if part of PkgGroup: for a single boundary condition system.
"""
lookup_and_forcing_files = [
self._force_absolute_path(file) for file in self.lookup_and_forcing_files
]
kwargs["lookup_and_forcing_files"] = lookup_and_forcing_files
return self._template_projectfile.render(**kwargs)
def check_lookup_and_forcing_files(self):
"""Check for presence of required MetaSWAP input files."""
filenames = {os.path.basename(f) for f in self.lookup_and_forcing_files}
missing_files = set(self._required_extra) - filenames
if len(missing_files) > 0:
raise ValueError(f"Missing MetaSWAP input files {missing_files}")
def _pkgcheck(self, active_cells=None):
# Dataset.dims does not return a tuple, like DataArray does.
# http://xarray.pydata.org/en/stable/generated/xarray.Dataset.dims.html
# Frozen(SortedKeysDict).keys() does not preserve ordering in keys
# Therefore convert to set and check for symmetric difference,
# to see if any dimensions are missing or unexpected dimensions are included.
dims = set(self.dataset.dims.keys())
if len(dims.symmetric_difference(("y", "x"))) > 0:
raise ValueError(f'Dataset dims not ("y", "x"), instead got {dims}')
self.check_lookup_and_forcing_files()