added Basic Model interface (BMI) and Docker Image for svat_tutorial.…

… Modified diagnostics to use output for BMI.

.gitignore

@@ -183,6 +183,7 @@ conda_build.sh
 build.sh
 meta.yaml
 
-# Docker
-/tmp
-Dockerfile
+# local development
+/bmiroger
+roger/bmimodels/
+examples/plot_scale/boadkh/

doc/equations/hydrological_model/storages/surface.rst

@@ -121,6 +121,8 @@ Short description of land use (lu_id)
 - `577`: summer grain pea
 - `578`: winter grain pea
 - `579`: winter rye
+- `580`: clover (growing only)
+- `581`: clover (continued)
 - `598`: no crop
 - `599`: bare
 - `6`: vineyard

examples/hillslope_scale/oneD_distributed_routing_tutorial/README.md

@@ -1,13 +1,16 @@
 # Tutorial for the spatially-distributed 1D-model
 
+Lateral transfer between grid cells occurs with a single flow direction. The flow direction is based on the topography
+
 Brief overview of the files:
 
 - `input`: Contains input data
 - `output`: Contains output data
 - `figures`: Contains figures
 - `config.yml`: File to set settings and output variables
+- `param_bounds.yml`: File to set settings and output variables
+- `write_parameters.py`: Generates random model parameters based on given parameter boundaries (`param_bounds.yml`) and writes model parameter file
 - `parameters.csv`: File that contains the model parameters
-- `write_parameters.py`: Generates model parameters and writes model parameter file
 - `oneD.py`: 1D-model setup
 - `merge_output.py`: Merges output into a single NetCDF file
 - `netcdf_to_csv.py`: Writes output to csv files

examples/hillslope_scale/oneD_distributed_routing_tutorial/config.yml

@@ -15,9 +15,9 @@ x_origin: 0
 y_origin: 0
 
 # --- Diagnostics --------------------------------------------------------------
+OUTPUT_FREQUENCY: 86400  # output frequency in seconds
 # simulated hydrologic fluxes
 OUTPUT_RATE:
-- "prec"
 - "aet"  # actual evapotranspiration (mm/day)
 - "transp"  # transpiration (mm/day)
 - "evap_soil"  # soil evaporation (mm/day)
@@ -27,24 +27,11 @@ OUTPUT_RATE:
 - "inf_ss"  # infiltration into subsoil (mm/day)
 - "q_rz"  # percolation from root zone (mm/day)
 - "q_ss"  # percolation from subsoil (mm/day)
-- "q_sub"  # lateral subsurface runoff (mm/day)
 - "cpr_rz" # capillary rise from subsoil into root zone (mm/day)
 - "dS"  # change of storage volume (mm/day)
 - "q_snow"  # snow melt (mm/day)
-- "q_sub_out"
-- "q_sub_in"
-- "q_sur_out"
-- "q_sur_in"
-- "q_sur"
-- "inf_mat_pot"
-- "inf_rz"
-- "inf_mp"
-- "q_sof"
-- "q_hof"
-- "q_sub"
-- "q_sub_rz"
-- "q_sub_ss"
-- "q_pot_rz"
+- "q_sur_out" # lateral surface outflow (mm/day)
+- "q_sub_out"  # lateral subsurface outflow (mm/day)
 
 # simulated hydrologic storages
 OUTPUT_COLLECT:
@@ -54,17 +41,5 @@ OUTPUT_COLLECT:
 - "theta_ss"  # soil water content of subsoil (-)
 - "theta"  # soil water content (-)
 - "S_snow"  # storage volume snow cover (mm)
-- "S_s"
-- "S"
-- "z_sat"
-- "z_wf"
-- "z0"
-- "S_fp_rz"
-- "S_fp_ss"
-- "S_lp_rz"
-- "S_lp_ss"
-- "S_ac_rz"
-- "S_ac_ss"
-- "S_ufc_rz"
-- "S_ufc_ss"
-- "S_zsat_rz"
+- "z_sat"  # saturation water level (mm)
+- "z0"  # surface water level (mm)

examples/hillslope_scale/oneD_distributed_routing_tutorial/oneD.py

@@ -58,6 +58,9 @@ def set_settings(self, state):
             settings = state.settings
             settings.identifier = self._config["identifier"]
 
+            # output frequency (in seconds)
+            settings.output_frequency = self._config["OUTPUT_FREQUENCY"]
+
             # total grid numbers in x- and y-direction
             settings.nx, settings.ny = self._config["nx"], self._config["ny"]
             # derive total number of time steps from forcing
@@ -545,14 +548,14 @@ def set_diagnostics(self, state, base_path=tmp_dir):
 
             diagnostics["rate"].output_variables = self._config["OUTPUT_RATE"]
             # values are aggregated to daily
-            diagnostics["rate"].output_frequency = 24 * 60 * 60  # in seconds
+            diagnostics["rate"].output_frequency = self._config["OUTPUT_FREQUENCY"]  # in seconds
             diagnostics["rate"].sampling_frequency = 1
             if base_path:
                 diagnostics["rate"].base_output_path = base_path
 
             diagnostics["collect"].output_variables = self._config["OUTPUT_COLLECT"]
             # values are aggregated to daily
-            diagnostics["collect"].output_frequency = 24 * 60 * 60  # in seconds
+            diagnostics["collect"].output_frequency = self._config["OUTPUT_FREQUENCY"]  # in seconds
             diagnostics["collect"].sampling_frequency = 1
             if base_path:
                 diagnostics["collect"].base_output_path = base_path
@@ -715,7 +718,7 @@ def after_timestep_kernel(state):
             at[2:-2, 2:-2, vs.taum1],
             vs.h[2:-2, 2:-2, vs.tau],
         )
-        # set to 0 for numerical errors
+        # set to 0 to avoid numerical errors
         vs.z0 = update(
             vs.z0,
             at[2:-2, 2:-2, vs.tau],
@@ -751,22 +754,25 @@ def after_timestep_kernel(state):
             at[vs.taum1],
             vs.doy[vs.tau],
         )
-        # set to 0 for numerical errors
+        # set to 0 to avoid numerical errors
         vs.S_fp_rz = update(
             vs.S_fp_rz,
             at[2:-2, 2:-2],
             npx.where((vs.S_fp_rz > -1e-6) & (vs.S_fp_rz < 0), 0, vs.S_fp_rz)[2:-2, 2:-2],
         )
+        # set to 0 to avoid numerical errors
         vs.S_lp_rz = update(
             vs.S_lp_rz,
             at[2:-2, 2:-2],
             npx.where((vs.S_lp_rz > -1e-6) & (vs.S_lp_rz < 0), 0, vs.S_lp_rz)[2:-2, 2:-2],
         )
+        # set to 0 to avoid numerical errors
         vs.S_fp_ss = update(
             vs.S_fp_ss,
             at[2:-2, 2:-2],
             npx.where((vs.S_fp_ss > -1e-6) & (vs.S_fp_ss < 0), 0, vs.S_fp_ss)[2:-2, 2:-2],
         )
+        # set to 0 to avoid numerical errors
         vs.S_lp_ss = update(
             vs.S_lp_ss,
             at[2:-2, 2:-2],

examples/hillslope_scale/oneD_distributed_tutorial/config.yml

@@ -15,6 +15,7 @@ x_origin: 0
 y_origin: 0
 
 # --- Diagnostics --------------------------------------------------------------
+OUTPUT_FREQUENCY: 86400  # output frequency in seconds
 # simulated hydrologic fluxes
 OUTPUT_RATE:
 - "prec"

examples/hillslope_scale/oneD_distributed_tutorial/oneD.py

@@ -58,6 +58,9 @@ def set_settings(self, state):
             settings = state.settings
             settings.identifier = self._config["identifier"]
 
+            # output frequency (in seconds)
+            settings.output_frequency = self._config["OUTPUT_FREQUENCY"]
+
             # total grid numbers in x- and y-direction
             settings.nx, settings.ny = self._config["nx"], self._config["ny"]
             # derive total number of time steps from forcing
@@ -207,19 +210,25 @@ def set_parameters_setup(self, state):
             vs.prec_weight = update(
                 vs.prec_weight,
                 at[2:-2, 2:-2],
-                self._read_var_from_csv("prec_weight", self._base_path, "parameters.csv").reshape(settings.nx, settings.ny),
+                self._read_var_from_csv("prec_weight", self._base_path, "parameters.csv").reshape(
+                    settings.nx, settings.ny
+                ),
             )
             # weight factor of air temperature (-)
             vs.ta_weight = update(
                 vs.ta_weight,
                 at[2:-2, 2:-2],
-                self._read_var_from_csv("ta_weight", self._base_path, "parameters.csv").reshape(settings.nx, settings.ny),
+                self._read_var_from_csv("ta_weight", self._base_path, "parameters.csv").reshape(
+                    settings.nx, settings.ny
+                ),
             )
             # weight factor of potential evapotranspiration (-)
             vs.pet_weight = update(
                 vs.pet_weight,
                 at[2:-2, 2:-2],
-                self._read_var_from_csv("pet_weight", self._base_path, "parameters.csv").reshape(settings.nx, settings.ny),
+                self._read_var_from_csv("pet_weight", self._base_path, "parameters.csv").reshape(
+                    settings.nx, settings.ny
+                ),
             )
 
         @roger_routine
@@ -304,13 +313,22 @@ def set_forcing(self, state):
                     vs.doy, at[1], self._read_var_from_nc("DOY", self._input_dir, "forcing.nc")[vs.itt_forc]
                 )
                 vs.prec_day = update(
-                    vs.prec_day, at[:, :, :], vs.PREC[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] * vs.prec_weight[:, :, npx.newaxis]
+                    vs.prec_day,
+                    at[:, :, :],
+                    vs.PREC[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24]
+                    * vs.prec_weight[:, :, npx.newaxis],
                 )
                 vs.ta_day = update(
-                    vs.ta_day, at[:, :, :], vs.TA[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] * vs.ta_weight[:, :, npx.newaxis]
+                    vs.ta_day,
+                    at[:, :, :],
+                    vs.TA[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24]
+                    * vs.ta_weight[:, :, npx.newaxis],
                 )
                 vs.pet_day = update(
-                    vs.pet_day, at[:, :, :], vs.PET[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] * vs.pet_weight[:, :, npx.newaxis]
+                    vs.pet_day,
+                    at[:, :, :],
+                    vs.PET[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24]
+                    * vs.pet_weight[:, :, npx.newaxis],
                 )
                 vs.itt_forc = vs.itt_forc + 6 * 24
 

examples/hillslope_scale/oneD_event_distributed_routing_tutorial/README.md

@@ -0,0 +1,112 @@
+# Tutorial for the spatially-distributed 1D-event-model
+
+Lateral transfer between grid cells occurs with a single flow direction. The flow direction is based on the topography.
+
+
+Brief overview of the files:
+
+- `input`: Contains input data
+- `output`: Contains output data
+- `figures`: Contains figures
+- `config.yml`: File to set settings and output variables
+- `parameters.csv`: File that contains the model parameters
+- `param_bounds.yml`: File to set settings and output variables
+- `write_parameters.py`: Generates random model parameters based on given parameter boundaries (`param_bounds.yml`) and writes model parameter file
+- `oneD_event.py`: 1D-event-model setup
+- `merge_output.py`: Merges output into a single NetCDF file
+- `netcdf_to_csv.py`: Writes output to csv files
+- `make_figures_and_tables.py`: Produces figures and tables
+
+# Date requirements
+
+The following information is required to run the model. 
+
+## Meteorological input data
+The required meteorological input data is loaded from the input folder. The input folder should contain a tab-delimited text file
+for precipitation (`PREC.txt`; 10 minutes time steps), air temperature (`TA.txt`; daily time steps) and potential evapotranspiration (`PET.txt`; daily time steps).
+
+Format of `PREC.txt` (PREC in mm/10 minutes):
+| YYYY  | MM    | DD    | hh    | mm    | PREC  |
+| ------| ------| ------| ------| ------| ------|
+| 2023  | 1     | 1     | 0     | 0     | 0     |
+| 2023  | 1     | 1     | 0     | 10    | 0.3   |        
+| 2023  | 1     | 1     | 0     | 20    | 1.0   |
+| 2023  | 1     | 1     | 0     | 30    | 0.5   |        
+| 2023  | 1     | 1     | 0     | 40    | 0.4   |
+| 2023  | 1     | 1     | 0     | 50    | 0.7   |        
+| 2023  | 1     | 1     | 1     | 0     | 0.6   |
+| ...   | ...   | ...   | ...   | ...   | ...   |
+
+Format of`TA.txt` (TA in degC):
+| YYYY  | MM    | DD    | hh    | mm    | TA    |
+| ------| ------| ------| ------| ------| ------|
+| 2023  | 1     | 1     | 0     | 0     | 2     |
+| 2023  | 1     | 1     | 0     | 0     | 3     |        
+| ...   | ...   | ...   | ...   | ...   | ...   |
+
+Format of `PET.txt` (PET in mm/day):
+| YYYY  | MM    | DD    | hh    | mm    | PET   |
+| ------| ------| ------| ------| ------| ------|
+| 2023  | 1     | 1     | 0     | 0     | 2     |
+| 2023  | 1     | 1     | 0     | 0     | 2.1   |        
+| ...   | ...   | ...   | ...   | ...   | ...   |
+
+
+where YYYY is the year, MM is the month, DD is the day, hh is the hour and mm is the minute.
+
+## Model parameters
+The model parameters (`parameters.csv`) are loaded from the same directory as `oneD.py`.
+
+Random model parameters can be generated running the following script:
+```
+python write_parameters.py
+```
+
+Format of `parameters.csv`:
+| lu_id | z_soil   | slope  | dmph | dmpv  | lmpv  | theta_pwp | theta_ufc | theta_ac | ks  | kf   |
+| ------| ---------| -------| -----| ------| ------| ----------| ----------| ---------| ----| -----|
+| 8     | 1000     | 0.05   | 25   | 25    | 200   | 0.2       | 0.11      | 0.09     | 5   | 2500 |  
+| 8     | 1000     | 0.06   | 20   | 30    | 300   | 0.18      | 0.1       | 0.08     | 6   | 2500 |
+| ...   | ...      | ...    | ...  | ...   | ...   | ...       | ...       | ...      | ... | ...  |
+
+where *lu_id* is the land cover, *z_soil* is the soil depth (mm), *slope* is the surface slope (-), *dmph* is the density of horizontal macropores (1/$m^2$), *dmpv* is the density of vertical macropores (1/$m^2$), *lmpv* is the length of vertical macropores (mm), *theta_pwp* is soil water content of the permanent wilting point (-), *theta_ufc* is soil water content of the usable field capacity (-), *theta_ac* is soil water content of the air capacity (-), *ks* is the saturated hydraulic conductivity (mm/hour) and *kf* is the hydraulic conductivity of the bedrock (mm/hour).
+
+## Model settings
+Name of model experiment and spatial discretization are defined in `config.yml`.
+
+## Model output
+The variables written to the output files are defined in `config.yml`. Available variables
+are listed [here](https://roger.readthedocs.io/en/latest/reference/variables.html#available-variables). Generally, storage variables
+should be defined for `OUTPUT_COLLECT` and flux variables for `OUTPUT_RATE`.
+
+# Workflow
+
+The following workflow briefly describes the model application:
+
+1. Prepare the meteorological input data (see Meteorological input data).
+
+2. Generates the model parameters for the simulation:
+```
+python write_parameters.py
+```
+
+3. If required data is ready, the following script runs the simulation:
+
+```
+python oneD.py
+```
+
+4. After calculation is done, the simulation results can be merged into a single NetCDF-file:
+```
+python merge_output.py
+```
+
+5. (Optional) The following script converts the model output from NetCDF to csv.
+```
+python netcdf_to_csv.py
+```
+
+6. Figures and tables can be produced with the following script:
+```
+python make_figures_and_tables.py
+``````

examples/hillslope_scale/oneD_event_distributed_routing_tutorial/config.yml

@@ -0,0 +1,40 @@
+# --- Settings --------------------------------------------------------------
+# name of the model run; is used to name the output files
+identifier: ONED-event
+# total grid numbers in x-direction
+nx: 1
+# total grid numbers in y-direction
+ny: 20
+# spatial discretization in x-direction (in meters)
+dx: 1
+# spatial discretization in y-direction (in meters)
+dy: 1
+# origin of spatial grid
+x_origin: 0
+# origin of spatial grid
+y_origin: 0
+
+# --- Diagnostics --------------------------------------------------------------
+OUTPUT_FREQUENCY: 600  # output frequency in seconds
+# simulated hydrologic fluxes
+OUTPUT_RATE:
+- "inf_mat_rz"  # soil matrix infiltration into root zone (mm/10 mins)
+- "inf_mp_rz"  # macropore infiltration into root zone (mm/10 mins)
+- "inf_sc_rz"  # shrinkage crack infiltration into root zone (mm/10 mins)
+- "inf_ss"  # infiltration into subsoil (mm/10 mins)
+- "q_rz"  # percolation from root zone (mm/10 mins)
+- "q_ss"  # percolation from subsoil (mm/10 mins)
+- "dS"  # change of storage volume (mm/10 mins)
+- "q_sur_out" # lateral surface outflow (mm/10 mins)
+- "q_sub_out"  # lateral subsurface outflow (mm/10 mins)
+
+# simulated hydrologic storages
+OUTPUT_COLLECT:
+- "S_rz"  # storage volume of root zone (mm)
+- "S_ss"  # storage volume of subsoil (mm)
+- "theta_rz"  # soil water content of root zone (-)
+- "theta_ss"  # soil water content of subsoil (-)
+- "theta"  # soil water content (-)
+- "S_snow"  # storage volume snow cover (mm)
+- "z_sat"  # saturation water level (mm)
+- "z0"  # surface water level (mm)