added weight factors for input data and linear reservoir for vadose z…

…one to delay groundwater recharge

.gitignore

@@ -123,8 +123,10 @@ startup.py
 *.xlsx
 
 # Examples
+examples/*/*/output/*
 examples/*/*/results/
 examples/*/*/figures/
+examples/*/*/*/output/*
 examples/*/*/*/results/
 examples/*/*/*/figures/
 examples/*/*/test_*

examples/hillslope_scale/oneD_distributed_tutorial/merge_output.py

@@ -23,11 +23,11 @@
     config = yaml.safe_load(file)
 
 # merge model output into single file
-path = str(base_path / f"{config['identifier']}.*.nc")
+path = str(base_path_output / f"{config['identifier']}.*.nc")
 diag_files = glob.glob(path)
-states_hm_file = base_path / f"{config['identifier']}.nc"
-if not os.path.exists(states_hm_file):
-    with h5netcdf.File(states_hm_file, "w", decode_vlen_strings=False) as f:
+output_file = base_path_output / f"{config['identifier']}.nc"
+if not os.path.exists(output_file):
+    with h5netcdf.File(output_file, "w", decode_vlen_strings=False) as f:
         f.attrs.update(
             date_created=datetime.datetime.today().isoformat(),
             title="RoGeR simulations (Tutorial)",

examples/hillslope_scale/oneD_distributed_tutorial/netcdf_to_csv.py

@@ -41,8 +41,9 @@ def main():
     # write output to csv
     vars_sim = config["OUTPUT_RATE"] + config["OUTPUT_COLLECT"]
     for var_sim in vars_sim:
-        df_var_sim = pd.DataFrame(index=date_sim, columns=[f"sim_{x}" for x in range(config["nx"])])
-        df_var_sim.iloc[:, :] = ds_sim[var_sim].isel(y=0).values.T
+        df_var_sim = pd.DataFrame(index=date_sim[1:], columns=[f"sim_{x}" for x in range(config["nx"] * config["ny"])])
+        df_var_sim.index = pd.Series(df_var_sim.index).dt.floor("D")
+        df_var_sim.iloc[:, :] = ds_sim[var_sim].values.reshape((len(date_sim), config["nx"] * config["ny"]))[1:, :]
         df_var_sim.to_csv(base_path_output / f"{var_sim}.csv", sep=";", index=True, index_label="date")
 
     return

examples/hillslope_scale/oneD_distributed_tutorial/oneD.py

@@ -203,6 +203,24 @@ def set_parameters_setup(self, state):
                 at[2:-2, 2:-2],
                 self._read_var_from_csv("kf", self._base_path, "parameters.csv").reshape(settings.nx, settings.ny),
             )
+            # weight factor of precipitation (-)
+            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),
+            )
+            # 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),
+            )
+            # 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),
+            )
 
         @roger_routine
         def set_parameters(self, state):
@@ -286,13 +304,13 @@ 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_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_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_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_distributed_tutorial/write_parameters.py

@@ -85,10 +85,14 @@ def main(nrows, ncols):
             # write parameters to dataframe
             df_params.loc[:, param] = values.flatten()
 
+    df_params.loc[:, 'prec_weight'] = 1
+    df_params.loc[:, 'ta_weight'] = 1
+    df_params.loc[:, 'pet_weight'] = 1
+
     # write parameters to csv
     df_params.columns = [
-        ["", "[mm]", "[-]", "[1/m2]", "[1/m2]", "[mm]", "[-]", "[-]", "[-]", "[mm/hour]", "[mm/hour]"],
-        ["lu_id", "z_soil", "slope", "dmph", "dmpv", "lmpv", "theta_ac", "theta_ufc", "theta_pwp", "ks", "kf"],
+        ["", "[mm]", "[-]", "[1/m2]", "[1/m2]", "[mm]", "[-]", "[-]", "[-]", "[mm/hour]", "[mm/hour]", "[-]", "[-]", "[-]"],
+        ["lu_id", "z_soil", "slope", "dmph", "dmpv", "lmpv", "theta_ac", "theta_ufc", "theta_pwp", "ks", "kf", "prec_weight", "ta_weight", "pet_weight"],
     ]
     df_params.to_csv(base_path / "parameters.csv", index=False, sep=";")
     return

examples/hillslope_scale/oneD_with_groundwater_boundary_condition_distributed_tutorial/config.yml

@@ -32,6 +32,7 @@ OUTPUT_RATE:
 - "cpr_ss" # capillary rise from groundwater into subsoil (mm/day)
 - "dS"  # change of storage volume (mm/day)
 - "q_snow"  # snow melt (mm/day)
+- "q_re"  # groundwater recharge (mm/day)
 
 # simulated hydrologic storages
 OUTPUT_COLLECT:

examples/hillslope_scale/oneD_with_groundwater_boundary_condition_distributed_tutorial/merge_output.py

@@ -4,6 +4,7 @@
 import datetime
 from pathlib import Path
 import numpy as onp
+import yaml
 import roger
 
 base_path = Path(__file__).parent
@@ -16,19 +17,24 @@
 if not os.path.exists(base_path_figs):
     os.mkdir(base_path_figs)
 
+# load configuration file
+file_config = base_path / "config.yml"
+with open(file_config, "r") as file:
+    config = yaml.safe_load(file)
+
 # merge model output into single file
-path = str(base_path / "SVAT.*.nc")
+path = str(base_path_output / f"{config['identifier']}.*.nc")
 diag_files = glob.glob(path)
-states_hm_file = base_path / "states_hm.nc"
-if not os.path.exists(states_hm_file):
-    with h5netcdf.File(states_hm_file, "w", decode_vlen_strings=False) as f:
+output_file = base_path_output / f"{config['identifier']}.nc"
+if not os.path.exists(output_file):
+    with h5netcdf.File(output_file, "w", decode_vlen_strings=False) as f:
         f.attrs.update(
             date_created=datetime.datetime.today().isoformat(),
             title="RoGeR simulations (Tutorial)",
             institution="University of Freiburg, Chair of Hydrology",
             references="",
             comment="First timestep (t=0) contains initial values. Simulations start are written from second timestep (t=1) to last timestep (t=N).",
-            model_structure="SVAT model with free drainage",
+            model_structure="1D model with groundwater boundary condition",
             roger_version=f"{roger.__version__}",
         )
         # collect dimensions

examples/hillslope_scale/oneD_with_groundwater_boundary_condition_distributed_tutorial/netcdf_to_csv.py

@@ -41,8 +41,9 @@ def main():
     # write output to csv
     vars_sim = config["OUTPUT_RATE"] + config["OUTPUT_COLLECT"]
     for var_sim in vars_sim:
-        df_var_sim = pd.DataFrame(index=date_sim, columns=[f"sim_{x}" for x in range(config["nx"])])
-        df_var_sim.iloc[:, :] = ds_sim[var_sim].isel(y=0).values.T
+        df_var_sim = pd.DataFrame(index=date_sim[1:], columns=[f"sim_{x}" for x in range(config["nx"] * config["ny"])])
+        df_var_sim.index = pd.Series(df_var_sim.index).dt.floor("D")
+        df_var_sim.iloc[:, :] = ds_sim[var_sim].values.reshape((len(date_sim), config["nx"] * config["ny"]))[1:, :]
         df_var_sim.to_csv(base_path_output / f"{var_sim}.csv", sep=";", index=True, index_label="date")
 
     return

examples/hillslope_scale/oneD_with_groundwater_boundary_condition_distributed_tutorial/oneD.py

@@ -204,6 +204,30 @@ def set_parameters_setup(self, state):
                 at[2:-2, 2:-2],
                 self._read_var_from_csv("kf", self._base_path, "parameters.csv").reshape(settings.nx, settings.ny),
             )
+            # specific yield at the soil surface (-)
+            vs.n0 = update(
+                vs.n0,
+                at[2:-2, 2:-2],
+                self._read_var_from_csv("n0", self._base_path, "parameters.csv").reshape(settings.nx, settings.ny),
+            )
+            # weight factor of precipitation (-)
+            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),
+            )
+            # 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),
+            )
+            # 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),
+            )
 
         @roger_routine
         def set_parameters(self, state):
@@ -212,9 +236,24 @@ def set_parameters(self, state):
             if (vs.month[vs.tau] != vs.month[vs.taum1]) & (vs.itt > 1):
                 vs.update(calc_parameters_surface_kernel(state))
 
-        @roger_routine
+        @roger_routine(
+            dist_safe=False,
+            local_variables=["S_vad"],
+        )
         def set_initial_conditions_setup(self, state):
-            pass
+            vs = state.variables
+
+            # storage of vadose zone (mm)
+            vs.S_vad = update(
+                vs.S_vad,
+                at[2:-2, 2:-2, :vs.taup1],
+                0.5 * vs.n0[2:-2, 2:-2, npx.newaxis] * (self._read_var_from_nc("Z_GW", self._input_dir, "forcing.nc")[0, 0, 0] - vs.z_soil[2:-2, 2:-2, npx.newaxis])
+            )
+            vs.S_vad = update(
+                vs.S_vad,
+                at[:, :, :],
+                npx.where(vs.S_vad < 0, 0, vs.S_vad),
+            )
 
         @roger_routine
         def set_initial_conditions(self, state):
@@ -284,13 +323,13 @@ 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_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_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_day, at[:, :, :], vs.PET[npx.newaxis, npx.newaxis, vs.itt_forc : vs.itt_forc + 6 * 24] * vs.pet_weight[:, :, npx.newaxis]
                 )
                 vs.z_gw = update(vs.z_gw, at[2:-2, 2:-2, vs.tau], vs.Z_GW[npx.newaxis, npx.newaxis, vs.itt_forc])
                 vs.itt_forc = vs.itt_forc + 6 * 24
@@ -391,6 +430,11 @@ def after_timestep_kernel(state):
             at[2:-2, 2:-2, vs.taum1],
             vs.S_s[2:-2, 2:-2, vs.tau],
         )
+        vs.S_vad = update(
+            vs.S_vad,
+            at[2:-2, 2:-2, vs.taum1],
+            vs.S_vad[2:-2, 2:-2, vs.tau],
+        )
         vs.S = update(
             vs.S,
             at[2:-2, 2:-2, vs.taum1],
@@ -515,6 +559,7 @@ def after_timestep_kernel(state):
             S_rz=vs.S_rz,
             S_ss=vs.S_ss,
             S_s=vs.S_s,
+            S_vad=vs.S_vad,
             S=vs.S,
             z_sat=vs.z_sat,
             z_wf=vs.z_wf,

examples/hillslope_scale/oneD_with_groundwater_boundary_condition_distributed_tutorial/param_bounds.yml

@@ -37,4 +37,8 @@ ks:
 # hydraulic conductivity of bedrock (mm/hour)
 kf:
 - 0.1
-- 10
+- 1000
+# specific yield at the soil surface (-)
+n0:
+- 0.1
+- 0.3