Add PTF for air entry pressure (hb) computation in setup_soilmaps

docs/changelog.rst

@@ -10,6 +10,7 @@ Unreleased
 ==========
 New
 ---
+- **setup_soilmaps**: add workflow to derive air entry pressure based on PTF from Rawls & Brakensiek or Clapp & Hornberger
 - Before writing the forcing to disk, throw an error for any missing values (#748)
 
 Fixed

examples/data/wflow_upgrade/sbm/wflow_sbm_v0x.toml

@@ -58,6 +58,7 @@ cf_soil = "cf_soil"
 cfmax = "Cfmax"
 e_r = "EoverR"
 f = "f"
+hb = "hb"
 infiltcappath = "InfiltCapPath"
 infiltcapsoil = "InfiltCapSoil"
 kext = "Kext"

examples/linux64/wflow_piave_clip/wflow_sbm.toml

@@ -83,6 +83,7 @@ soil__thickness = "soil_thickness"
 soil_layer_water__brooks_corey_exponent = "soil_brooks_corey_c"
 soil_surface_water__vertical_saturated_hydraulic_conductivity = "soil_ksat_vertical"
 soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter = "soil_f"
+soil_water__air_entry_pressure_head = "soil_hb"
 
 [input.static.subsurface_water__horizontal_to_vertical_saturated_hydraulic_conductivity_ratio]
 value = 100

examples/linux64/wflow_piave_subbasin/wflow_sbm.toml

@@ -99,6 +99,7 @@ soil__thickness = "soil_thickness"
 soil_layer_water__brooks_corey_exponent = "soil_brooks_corey_c"
 soil_surface_water__vertical_saturated_hydraulic_conductivity = "soil_ksat_vertical"
 soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter = "soil_f"
+soil_water__air_entry_pressure_head = "soil_hb"
 
 [input.static.subsurface_water__horizontal_to_vertical_saturated_hydraulic_conductivity_ratio]
 value = 100

examples/linux64/wflow_piave_subbasin/wflow_sbm_results.toml

@@ -103,6 +103,7 @@ soil__thickness = "soil_thickness"
 soil_layer_water__brooks_corey_exponent = "soil_brooks_corey_c"
 soil_surface_water__vertical_saturated_hydraulic_conductivity = "soil_ksat_vertical"
 soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter = "soil_f"
+soil_water__air_entry_pressure_head = "soil_hb"
 subsurface_water__horizontal_to_vertical_saturated_hydraulic_conductivity_ratio.value = 100
 snowpack__degree_day_coefficient.value = 3.75653
 soil_surface_water__infiltration_reduction_parameter.value = 0.038

examples/linux64/wflow_piave_subbasin/wflow_sbm_results2.toml

@@ -103,6 +103,7 @@ soil__thickness = "soil_thickness"
 soil_layer_water__brooks_corey_exponent = "soil_brooks_corey_c"
 soil_surface_water__vertical_saturated_hydraulic_conductivity = "soil_ksat_vertical"
 soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter = "soil_f"
+soil_water__air_entry_pressure_head = "soil_hb"
 subsurface_water__horizontal_to_vertical_saturated_hydraulic_conductivity_ratio.value = 10
 snowpack__degree_day_coefficient.value = 3.75653
 soil_surface_water__infiltration_reduction_parameter.value = 0.038

examples/wflow_build.yml

@@ -73,6 +73,7 @@ steps:
   - setup_soilmaps:
       soil_fn: soilgrids     # source for soilmaps: {soilgrids}
       ptf_ksatver: brakensiek    # pedotransfer function to calculate hydraulic conductivity: {brakensiek, cosby}
+      ptf_hb: brakensiek    # pedotranfer function to calculate air entry pressure: {brakensiek, clapp}
 
   - setup_outlets:
       river_only: True

examples/wflow_build_basin.yml

@@ -71,6 +71,7 @@ steps:
   - setup_soilmaps:
       soil_fn: soilgrids     # source for soilmaps: {soilgrids}
       ptf_ksatver: brakensiek    # pedotransfer function to calculate hydraulic conductivity: {brakensiek, cosby}
+      ptf_hb: brakensiek    # pedotranfer function to calculate air entry pressure: {brakensiek, clapp}
 
   - setup_outlets:
       river_only: True

examples/wflow_piave_clip/wflow_sbm.toml

@@ -83,6 +83,7 @@ soil__thickness = "soil_thickness"
 soil_layer_water__brooks_corey_exponent = "soil_brooks_corey_c"
 soil_surface_water__vertical_saturated_hydraulic_conductivity = "soil_ksat_vertical"
 soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter = "soil_f"
+soil_water__air_entry_pressure_head = "soil_hb"
 
 [input.static.subsurface_water__horizontal_to_vertical_saturated_hydraulic_conductivity_ratio]
 value = 100

examples/wflow_piave_subbasin/wflow_sbm.toml

@@ -99,6 +99,7 @@ soil__thickness = "soil_thickness"
 soil_layer_water__brooks_corey_exponent = "soil_brooks_corey_c"
 soil_surface_water__vertical_saturated_hydraulic_conductivity = "soil_ksat_vertical"
 soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter = "soil_f"
+soil_water__air_entry_pressure_head = "soil_hb"
 
 [input.static.subsurface_water__horizontal_to_vertical_saturated_hydraulic_conductivity_ratio]
 value = 100

examples/wflow_piave_subbasin/wflow_sbm_results.toml

@@ -103,6 +103,7 @@ soil__thickness = "soil_thickness"
 soil_layer_water__brooks_corey_exponent = "soil_brooks_corey_c"
 soil_surface_water__vertical_saturated_hydraulic_conductivity = "soil_ksat_vertical"
 soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter = "soil_f"
+soil_water__air_entry_pressure_head = "soil_hb"
 subsurface_water__horizontal_to_vertical_saturated_hydraulic_conductivity_ratio.value = 100
 snowpack__degree_day_coefficient.value = 3.75653
 soil_surface_water__infiltration_reduction_parameter.value = 0.038

examples/wflow_piave_subbasin/wflow_sbm_results2.toml

@@ -103,6 +103,7 @@ soil__thickness = "soil_thickness"
 soil_layer_water__brooks_corey_exponent = "soil_brooks_corey_c"
 soil_surface_water__vertical_saturated_hydraulic_conductivity = "soil_ksat_vertical"
 soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter = "soil_f"
+soil_water__air_entry_pressure_head = "soil_hb"
 subsurface_water__horizontal_to_vertical_saturated_hydraulic_conductivity_ratio.value = 10
 snowpack__degree_day_coefficient.value = 3.75653
 soil_surface_water__infiltration_reduction_parameter.value = 0.038

hydromt_wflow/naming.py

@@ -253,6 +253,10 @@
         "wflow_v0": "vertical.rootdistpar",
         "wflow_v1": "soil_wet_root__sigmoid_function_shape_parameter",
     },
+    "soil_hb": {
+        "wflow_v0": "vertical.hb",
+        "wflow_v1": "soil_water__air_entry_pressure_head",
+    },
     "soil_thickness": {
         "wflow_v0": "vertical.soilthickness",
         "wflow_v1": "soil__thickness",

hydromt_wflow/wflow_sbm.py

@@ -1128,7 +1128,6 @@ def setup_glaciers(
         * **glacier_initial_leq_depth** map: storage (volume) of glacier per cell [mm]
 
         Required setup methods:
-
         * :py:meth:`~WflowSbmModel.setup_basemaps`
 
         Parameters
@@ -1924,6 +1923,7 @@ def setup_soilmaps(
         self,
         soil_fn: str = "soilgrids",
         ptf_ksatver: str = "brakensiek",
+        ptf_hb: str = "brakensiek",
         wflow_thicknesslayers: list[int] = [100, 300, 800],
         output_names: dict = {
             "soil_water__saturated_volume_fraction": "soil_theta_s",
@@ -1932,6 +1932,7 @@ def setup_soilmaps(
             "soil__thickness": "soil_thickness",
             "soil_water__vertical_saturated_hydraulic_conductivity_scale_parameter": "soil_f",  # noqa: E501
             "soil_layer_water__brooks_corey_exponent": "soil_brooks_corey_c",
+            "soil_water__air_entry_pressure_head": "soil_hb",
         },
     ):
         """
@@ -1943,8 +1944,9 @@ def setup_soilmaps(
 
         Currently, supported ``soil_fn`` is "soilgrids" and "soilgrids_2020".
         ``ptf_ksatver`` (PTF for the vertical hydraulic conductivity) options are
-        "brakensiek" and "cosby". "soilgrids" provides data at 7 specific depths,
-        while "soilgrids_2020" provides data averaged over 6 depth intervals.
+        "brakensiek" and "cosby". ``ptf_hb`` (PTF for the air entry pressure)
+        options are "brakensiek" and "clapp". "soilgrids" provides data at 7 specific
+        depths, while "soilgrids_2020" provides data averaged over 6 depth intervals.
         This leads to small changes in the workflow:
         (1) M parameter uses midpoint depths in soilgrids_2020 versus \
 specific depths in soilgrids,
@@ -1987,6 +1989,9 @@ def setup_soilmaps(
         * **meta_soil_texture** map: soil texture based on USDA soil texture triangle \
 (mapping: [1:Clay, 2:Silty Clay, 3:Silty Clay-Loam, 4:Sandy Clay, 5:Sandy Clay-Loam, \
 6:Clay-Loam, 7:Silt, 8:Silt-Loam, 9:Loam, 10:Sand, 11: Loamy Sand, 12:Sandy Loam])
+        * **soil_water__air_entry_pressure_head** map:
+            The air entry pressure [cm] at which air begins to enter the largest \
+            pores of a saturated soil during drying.
 
 
         Required setup methods:
@@ -2007,6 +2012,10 @@ def setup_soilmaps(
             Pedotransfer function (PTF) to use for calculation of ksat vertical
             (vertical saturated hydraulic conductivity [mm/day]).
             By default 'brakensiek'.
+        ptf_hb : {'brakensiek', 'clapp'}
+            Pedotransfer function (PTF) method to use for calculation of hb
+            (air entry pressure [cm]).
+            By default 'brakensiek'.
         wflow_thicknesslayers : list of int, optional
             Thickness of soil layers [mm] for wflow_sbm soil model.
             By default [100, 300, 800] for layers at depths 100, 400, 1200 and >1200 mm.
@@ -2025,6 +2034,7 @@ def setup_soilmaps(
             ds=dsin,
             ds_like=self.staticmaps.data,
             ptfKsatVer=ptf_ksatver,
+            ptf_hb=ptf_hb,
             soil_fn=soil_fn,
             wflow_layers=wflow_thicknesslayers,
         ).reset_coords(drop=True)

hydromt_wflow/workflows/ptf.py

@@ -471,3 +471,141 @@ def mean_diameter_soil(clay, silt):
     d50 = d50 / 1000  # [mm]
 
     return d50
+
+
+def hb_brakensiek(clay, sand, porosity):
+    """
+    Determine air entry pressure for soils with 5-60% clay and 5-70% sand.
+
+    The air entry pressure for soils that do not meet these requirements are
+    computed using the PTF from Clapp & Hornberger (1987).
+    Note that in Rawls & Brakensiek (1985) the air entry pressure
+    is referred to as the Brooks-Corey bubbling pressure.
+
+    Based on:
+      Rawls, W.J., and D.L. Brakensiek. 1985. Prediction of soil water
+      properties for hydrologic modeling. p. 293-299. In E.B. Jones and
+      T.J.Ward (ed.) Proc. Symp.Watershed Management in the Eighties,
+      Denver, CO. 30 Apr.-1 May 1985. Am. Soc. Civil Eng., New York mapping
+
+    Parameters
+    ----------
+    clay: float
+        clay percentage [%].
+    sand: float
+        sand percentage [%].
+    porosity: float
+        porosity of the soil [-].
+
+    Returns
+    -------
+    air_entry_pressure : float
+        based on equation from Rawls & Brakensiek (1985).
+    """
+    air_entry_pressure = np.where(
+        np.logical_and(
+            np.logical_and(clay > 5.0, clay < 60), np.logical_and(sand > 5, sand < 70)
+        ),
+        -np.exp(
+            5.3396738
+            + 0.1845038 * clay
+            - 2.48394546 * porosity
+            - 0.00213853 * (clay**2)
+            - 0.04356349 * sand * porosity
+            - 0.61745089 * clay * porosity
+            + 0.00143598 * (sand**2) * (porosity**2)
+            - 0.00855375 * (clay**2) * (porosity**2)
+            - 0.00001282 * (sand**2) * clay
+            + 0.00895359 * (clay**2) * porosity
+            - 0.00072472 * (sand**2) * porosity
+            + 0.0000054 * (clay**2) * sand
+            + 0.50028060 * (porosity**2) * clay
+        ),
+        hb_clapp(clay, sand),
+    )
+
+    return air_entry_pressure
+
+
+def hb_clapp(clay, sand):
+    """
+    Determine air entry pressure from Clapp & Hornberger (1978), Table 2.
+
+    Note that in Clapp & Hornberger (1978) the air entry pressure is referred
+    to as ψ_s.
+
+    Based on:
+        Clapp, R. B., and G. M. Hornberger (1978), Empirical equations for
+        some soil hydraulic properties, Water Resour. Res., 14(4),
+        601–604, doi:10.1029/WR014i004p00601.
+
+    Parameters
+    ----------
+    clay: float
+        clay percentage [%].
+    sand: float
+        sand percentage [%].
+
+    Returns
+    -------
+    air_entry_pressure : float
+        based on equation from Clapp & Hornberger (1978).
+    """
+    silt = 100 - (clay + sand)
+
+    air_entry_pressure = np.where(
+        np.logical_and(clay >= 40.0, sand >= 20.0, sand <= 45),
+        -40.5,  # clay
+        np.where(
+            np.logical_and(clay >= 27.0, sand >= 20.0, sand <= 45),
+            -63.0,  # clay loam
+            np.where(
+                np.logical_and(silt <= 40.0, sand <= 20.0),
+                -40.5,  # clay
+                np.where(
+                    np.logical_and(silt > 40.0, clay >= 40.0),
+                    -49.0,  # silty clay
+                    np.where(
+                        np.logical_and(clay >= 35.0, sand >= 45.0),
+                        -15.3,  # sandy clay
+                        np.where(
+                            np.logical_and(clay >= 27.0, sand < 20.0),
+                            -35.6,  # silty clay loam
+                            np.where(
+                                np.logical_and(clay <= 10.0, silt >= 80.0),
+                                -78.6,  # silt (value not included in paper,
+                                # so use value for silt loam)
+                                np.where(
+                                    (silt >= 50.0),
+                                    -78.6,  # silt loam
+                                    np.where(
+                                        np.logical_and(
+                                            clay >= 7.0, sand <= 52.0, silt >= 28.0
+                                        ),
+                                        -47.8,  # loam
+                                        np.where(
+                                            (clay >= 20.0),
+                                            -29.9,  # sandy clay loam
+                                            np.where(
+                                                (clay >= (sand - 70)),
+                                                -21.8,  # sandy loam
+                                                np.where(
+                                                    (clay >= (2 * sand - 170.0)),
+                                                    -9.0,  # loamy sand
+                                                    np.where(
+                                                        np.isnan(clay), np.nan, -12.1
+                                                    ),  # sand
+                                                ),
+                                            ),
+                                        ),
+                                    ),
+                                ),
+                            ),
+                        ),
+                    ),
+                ),
+            ),
+        ),
+    )
+
+    return air_entry_pressure