Restore observed soil moisture forcing (GH-3)

CHANGELOG.md

@@ -157,6 +157,18 @@
   - Same functionality, lighter dependencies
   - Consolidates all parallel processing to one library
 
+### 1 Dec 2025
+- [feature] Added site-level soil observation configuration (GH-3 improvement)
+  - New `soil_observation` block in site properties for cleaner YAML configuration
+  - Soil observation metadata now correctly modelled as a site-level property (not per-surface)
+  - Maintains backwards compatibility with legacy per-surface configuration
+  - Documentation updated with preferred YAML approach and legacy fallback
+
+### 16 Nov 2025
+- [maintenance] Reorganised soil observation conversion logic from top-level `_soil_obs.py` to `util/_forcing.py`
+  - Follows SuPy's established pattern where specialized utilities live in the `util/` subdirectory
+  - No functional changes; purely organisational refactoring for better code structure
+
 ### 14 Nov 2025
 - [feature] Added `SUEWSSimulation.from_sample_data()` factory method and comprehensive OOP enhancements (#779)
   - New factory method for cleaner OOP workflow: `sim = SUEWSSimulation.from_sample_data()`
@@ -171,6 +183,11 @@
   - Four validation layers: (1) basic range [1, 366], (2) consistency (both set or both None), (3) leap year (DOY 366 only in leap years), (4) hemisphere pattern check (NH/SH typical ranges)
   - First three layers raise ERROR; hemisphere check adds INFO to report "NO ACTION NEEDED" section
   - Useful when Phase C runs standalone or via `SUEWSConfig.from_yaml()` (Phase B auto-corrects values in full pipeline)
+- [bugfix] Enabled observed soil moisture forcing (GH-3)
+  - Added soil observation metadata fields (`obs_sm_depth`, `obs_sm_cap`, `obs_soil_not_rocks`, `soildensity`) to the land-cover definition
+  - SuPy now converts observed volumetric/gravimetric `xsmd` data to soil moisture deficits before calling the SUEWS kernel
+  - Simplified approach: metadata only needed on surface 0 (Paved); other surfaces ignored
+  - Documentation updated to reflect the required inputs when `SMDMethod` = 1 or 2
 
 ### 12 Nov 2025
 - [feature] Added irrigation year-wrapping pattern detection (#843)

docs/source/inputs/tables/SUEWS_SiteInfo/Input_Options.rst

@@ -2673,49 +2673,6 @@ Input Options
 			:widths: 44 18 38
 
 
-.. option:: OBS_SMCap
-
-	:Description:
-
-		The maximum observed soil moisture. [|m^3| |m^-3| or kg |kg^-1|]
-
-	:Configuration:
-		.. csv-table::
-			:class: longtable
-			:file: csv-table/OBS_SMCap.csv
-			:header-rows: 1
-			:widths: 44 18 38
-
-
-.. option:: OBS_SMDepth
-
-	:Description:
-
-		The depth of soil moisture measurements. [mm]
-
-
-	:Configuration:
-		.. csv-table::
-			:class: longtable
-			:file: csv-table/OBS_SMDepth.csv
-			:header-rows: 1
-			:widths: 44 18 38
-
-
-.. option:: OBS_SoilNotRocks
-
-	:Description:
-
-		Fraction of soil without rocks. [-]
-
-	:Configuration:
-		.. csv-table::
-			:class: longtable
-			:file: csv-table/OBS_SoilNotRocks.csv
-			:header-rows: 1
-			:widths: 44 18 38
-
-
 .. option:: OHMCode_SummerDry
 
 	:Description:

docs/source/inputs/tables/SUEWS_SiteInfo/SUEWS_Soil.rst

@@ -12,12 +12,8 @@ Each of the non-water surface types need to link to soil characteristics specifi
 If the soil characteristics are assumed to be the same for all surface types, use a single code value to link the characteristics here with the SoilTypeCode columns in `SUEWS_NonVeg.txt` and `SUEWS_Veg.txt`.
 
 Soil moisture can either be provided using observational data in the met
-forcing file (the `xsmd` column when `SMDMethod` = 1 or 2 in `RunControl.nml`) and providing some soil properties here, or modelled by SUEWS (`SMDMethod` = 0 in `RunControl.nml`).
-
-
-.. .. caution::
-..   The option to use observational data is not operational in the current release!
-
+forcing file (the `xsmd` column when `SMDMethod` = 1 or 2 in `RunControl.nml`), or modelled by SUEWS (`SMDMethod` = 0 in `RunControl.nml`).
+When using observed soil moisture, see the YAML configuration documentation for required site-level metadata.
 
 .. DON'T manually modify the csv file below
 .. as it is always automatically regenrated by each build:

docs/source/inputs/tables/SUEWS_SiteInfo/csv-table/SUEWS_Soil.csv

@@ -5,6 +5,3 @@ No.,Column Name,Use,Description
 4,`SatHydraulicCond`,`MD`,Hydraulic conductivity for saturated soil [mm |s^-1|]
 5,`SoilDensity`,`MD`,Soil density [kg |m^-3|]
 6,`InfiltrationRate`,`O`,Infiltration rate.
-7,`OBS_SMDepth`,`O`,The depth of soil moisture measurements. [mm]
-8,`OBS_SMCap`,`O`,The maximum observed soil moisture. [|m^3| |m^-3| or kg |kg^-1|]
-9,`OBS_SoilNotRocks`,`O`,Fraction of soil without rocks. [-]

docs/source/inputs/tables/SUEWS_SiteInfo/sample-table/SUEWS_Soil.txt

@@ -1,8 +1,8 @@
-1    	 2         	 3            	 4                	 5           	 6                	 7           	 8         	 9
-Code 	 SoilDepth 	 SoilStoreCap 	 SatHydraulicCond 	 SoilDensity 	 InfiltrationRate 	 OBS_SMDepth 	 OBS_SMCap 	 OBS_SoilNotRocks
-551  	 350       	 150          	 5.00E-04         	 -999        	 -999             	 -999        	 -999      	 -999             	 ! 	 Swindon   (below   Paved)  	 Ward   et   al.   (2013)
-552  	 350       	 150          	 5.00E-04         	 -999        	 -999             	 -999        	 -999      	 -999             	 ! 	 Swindon   (below   Built)  	 Ward   et   al.   (2013)
-553  	 350       	 150          	 5.00E-04         	 -999        	 -999             	 -999        	 -999      	 -999             	 ! 	 Swindon   (below   others) 	 Ward   et   al.   (2013)
-661  	 350       	 150          	 5.00E-04         	 -999        	 -999             	 -999        	 -999      	 -999             	 ! 	 London
+1    	 2         	 3            	 4                	 5           	 6
+Code 	 SoilDepth 	 SoilStoreCap 	 SatHydraulicCond 	 SoilDensity 	 InfiltrationRate
+551  	 350       	 150          	 5.00E-04         	 -999        	 -999             	 ! 	 Swindon   (below   Paved)  	 Ward   et   al.   (2013)
+552  	 350       	 150          	 5.00E-04         	 -999        	 -999             	 ! 	 Swindon   (below   Built)  	 Ward   et   al.   (2013)
+553  	 350       	 150          	 5.00E-04         	 -999        	 -999             	 ! 	 Swindon   (below   others) 	 Ward   et   al.   (2013)
+661  	 350       	 150          	 5.00E-04         	 -999        	 -999             	 ! 	 London
 -9
 -9

docs/source/inputs/tables/SUEWS_SiteInfo/typical-general.csv

@@ -82,7 +82,4 @@ SoilStoreCap,Soil,150,Capacity of sub-surface soil store [mm],
 ,,,,
 SatHydraulicCond,Soil,0.0005,Hydraulic conductivity for saturated soil [mm s-1],
 SoilDensity,Soil,1.16,Soil density [kg m-3],
-InfiltrationRate,Soil,,Infiltration rate [mm h-1],
-OBS_SMDepth,Soil,,Depth of soil moisture measurements [mm],
-OBS_SMCap,Soil,,Maxiumum observed soil moisture [m3 m-3 or kg kg-1],
-OBS_SoilNotRocks,Soil,,Fraction of soil without rocks [-],
+InfiltrationRate,Soil,,Infiltration rate [mm h-1],

src/supy/_run.py

@@ -41,6 +41,7 @@
 from ._version import __version__ as sp_version
 
 from ._env import logger_supy
+from .util._forcing import convert_observed_soil_moisture
 
 from .util._debug import save_zip_debug
 
@@ -392,6 +393,9 @@ def run_supy_ser(
     ]
     df_forcing = df_forcing.loc[:, list_var_forcing]
 
+    # Convert observed soil moisture to deficits (if required)
+    df_forcing = convert_observed_soil_moisture(df_forcing, df_init)
+
     # grid list determined by initial states
     list_grid = df_init.index
 

src/supy/data_model/core/site.py

@@ -2328,6 +2328,84 @@ def from_df_state(cls, df: pd.DataFrame, grid_id: int) -> "SPARTACUSParams":
         return cls(**params)
 
 
+class SoilObservationConfig(BaseModel):
+    """Configuration for observed soil moisture measurements.
+
+    Required when SMDMethod = 1 (volumetric) or 2 (gravimetric).
+    These parameters describe the sensor installation and measurement setup,
+    not the physical soil properties used by the SUEWS water balance model.
+
+    Note:
+        - `depth`: Where the sensor is installed, not SUEWS's modelled soil depth
+        - `smcap`: Maximum observable moisture, not SUEWS's soil storage capacity
+        - `bulk_density`: Only used for gravimetric (SMDMethod=2) conversion
+    """
+
+    model_config = ConfigDict(title="Soil Observation Configuration")
+
+    depth: FlexibleRefValue(float) = Field(
+        gt=0,
+        description="Depth of the soil moisture sensor installation",
+        json_schema_extra={"unit": "mm", "display_name": "Observation Depth"},
+    )
+    smcap: FlexibleRefValue(float) = Field(
+        gt=0,
+        description="Maximum observable soil moisture (saturated θ or w at sensor location)",
+        json_schema_extra={"unit": "fraction", "display_name": "Saturation Capacity"},
+    )
+    soil_not_rocks: FlexibleRefValue(float) = Field(
+        gt=0,
+        le=1,
+        description="Fraction of the measured soil volume that is soil (not rocks)",
+        json_schema_extra={"unit": "dimensionless", "display_name": "Soil Fraction"},
+    )
+    bulk_density: FlexibleRefValue(float) = Field(
+        gt=0,
+        description="Soil bulk density at the measurement point (used for gravimetric conversion)",
+        json_schema_extra={"unit": "g cm^-3", "display_name": "Bulk Density"},
+    )
+
+    ref: Optional[Reference] = None
+
+    def to_df_state(self, grid_id: int) -> pd.DataFrame:
+        """Convert soil observation config to DataFrame state format."""
+        df_state = init_df_state(grid_id)
+
+        for attr in ["depth", "smcap", "soil_not_rocks", "bulk_density"]:
+            field_val = getattr(self, attr)
+            val = field_val.value if isinstance(field_val, RefValue) else field_val
+            # Use obs_sm_ prefix for clarity in df_state
+            col_name = (
+                f"obs_sm_{attr}" if attr != "bulk_density" else "obs_sm_bulk_density"
+            )
+            df_state[(col_name, "0")] = val
+
+        return df_state
+
+    @classmethod
+    def from_df_state(cls, df: pd.DataFrame, grid_id: int) -> "SoilObservationConfig":
+        """Create SoilObservationConfig from DataFrame state format."""
+        params = {}
+        field_mapping = {
+            "depth": "obs_sm_depth",
+            "smcap": "obs_sm_smcap",
+            "soil_not_rocks": "obs_sm_soil_not_rocks",
+            "bulk_density": "obs_sm_bulk_density",
+        }
+
+        for field_name, col_name in field_mapping.items():
+            if (col_name, "0") in df.columns:
+                val = df.loc[grid_id, (col_name, "0")]
+                if val is not None and not pd.isna(val) and val > -998:
+                    params[field_name] = RefValue(val)
+
+        # Return None if no valid parameters found
+        if not params or len(params) < 4:
+            return None
+
+        return cls(**params)
+
+
 class LUMPSParams(BaseModel):
     """LUMPS model parameters for surface moisture."""
 
@@ -2527,6 +2605,10 @@ class SiteProperties(BaseModel):
     snow: SnowParams = Field(
         default_factory=SnowParams, description="Parameters for snow modelling"
     )
+    soil_observation: Optional[SoilObservationConfig] = Field(
+        default=None,
+        description="Soil moisture observation metadata (required if SMDMethod > 0)",
+    )
     land_cover: LandCover = Field(
         default_factory=LandCover,
         description="Parameters for land cover characteristics",
@@ -2630,22 +2712,26 @@ def to_df_state(self, grid_id: int) -> pd.DataFrame:
         df_stebbs = self.stebbs.to_df_state(grid_id)
         df_building_archetype = self.building_archetype.to_df_state(grid_id)
 
-        df_state = pd.concat(
-            [
-                df_state,
-                df_lumps,
-                df_spartacus,
-                df_conductance,
-                df_irrigation,
-                df_anthropogenic_emissions,
-                df_snow,
-                df_land_cover,
-                df_vertical_layers,
-                df_stebbs,
-                df_building_archetype,
-            ],
-            axis=1,
-        )
+        dfs_to_concat = [
+            df_state,
+            df_lumps,
+            df_spartacus,
+            df_conductance,
+            df_irrigation,
+            df_anthropogenic_emissions,
+            df_snow,
+            df_land_cover,
+            df_vertical_layers,
+            df_stebbs,
+            df_building_archetype,
+        ]
+
+        # Optional soil observation config (only when SMDMethod > 0)
+        if self.soil_observation is not None:
+            df_soil_obs = self.soil_observation.to_df_state(grid_id)
+            dfs_to_concat.append(df_soil_obs)
+
+        df_state = pd.concat(dfs_to_concat, axis=1)
         return df_state
 
     @classmethod
@@ -2705,6 +2791,9 @@ def from_df_state(cls, df: pd.DataFrame, grid_id: int) -> "SiteProperties":
         params["stebbs"] = StebbsProperties.from_df_state(df, grid_id)
         params["building_archetype"] = ArchetypeProperties.from_df_state(df, grid_id)
 
+        # Optional soil observation config (returns None if not present)
+        params["soil_observation"] = SoilObservationConfig.from_df_state(df, grid_id)
+
         return cls(**params)
 
 

src/supy/data_model/core/surface.py

@@ -254,6 +254,14 @@ class SurfaceProperties(BaseModel):
             "display_name": "Saturated Hydraulic Conductivity",
         },
     )
+    soildensity: Optional[FlexibleRefValue(float)] = Field(
+        default=None,
+        description="Bulk soil density",
+        json_schema_extra={
+            "unit": "g cm^-3",
+            "display_name": "Soil Density",
+        },
+    )
     waterdist: Optional[WaterDistribution] = Field(
         default=None,  # TODO: Can this be None?
         description="Water distribution parameters",
@@ -346,6 +354,7 @@ def set_df_value(col_name: str, value: float):
             "statelimit",
             "wetthresh",
             "sathydraulicconduct",
+            "soildensity",
             "waterdist",
             "storedrainprm",
             "snowpacklimit",
@@ -413,6 +422,7 @@ def set_df_value(col_name: str, value: float):
                     defaults = {
                         "soildepth": 150.0,
                         "sathydraulicconduct": 0.0001,
+                        "soildensity": -999.0,
                     }
                     value = defaults.get(property, 0.0)
                 set_df_value(property, value)
@@ -464,6 +474,7 @@ def from_df_state(
             "statelimit",
             "wetthresh",
             "sathydraulicconduct",
+            "soildensity",
             "waterdist",
             "storedrainprm",
             "snowpacklimit",
@@ -521,8 +532,17 @@ def from_df_state(
                 value = df.loc[grid_id, ("kkanohm", f"({surf_idx},)")]
                 property_values["k_anohm"] = RefValue(value)
             else:
-                value = df.loc[grid_id, (property, f"({surf_idx},)")]
-                property_values[property] = RefValue(value)
+                # Check if column exists (for backwards compatibility with old tables)
+                col_key = (property, f"({surf_idx},)")
+                if col_key in df.columns:
+                    value = df.loc[grid_id, col_key]
+                    property_values[property] = RefValue(value)
+                else:
+                    # Column doesn't exist - skip if optional, raise if required
+                    field_info = cls.model_fields.get(property)
+                    if field_info and not field_info.is_required():
+                        property_values[property] = None
+                    # If required, let Pydantic validation catch it later
 
         return cls(**property_values)
 

src/supy/meson.build

@@ -36,6 +36,7 @@ py.install_sources(
             'util/_atm.py',
             'util/_debug.py',
             'util/_era5.py',
+            'util/_forcing.py',
             'util/code_manager.py',
             'util/conversion_engine.py',
             'util/_gap_filler.py',