diff --git a/examples/wflow_sbmgw_build.ini b/examples/wflow_sbmgw_build.ini
new file mode 100644
index 00000000..e93b803a
--- /dev/null
+++ b/examples/wflow_sbmgw_build.ini
@@ -0,0 +1,98 @@
+[global]
+data_libs       = []            # add optional paths to data yml files
+
+[setup_config]                  # options parsed to wflow ini file <section>.<option>
+starttime = 2010-01-01T00:00:00
+endtime = 2010-03-31T00:00:00
+timestepsecs = 86400
+input.path_forcing = inmaps-era5-2010.nc
+
+[setup_basemaps]
+hydrography_fn     = merit_hydro   # source hydrography data {merit_hydro, merit_hydro_1k}
+basin_index_fn  = merit_hydro_index # source of basin index corresponding to hydrography_fn
+upscale_method  = ihu           # upscaling method for flow direction data, by default 'ihu'
+
+[setup_rivers]
+hydrography_fn  = merit_hydro   # source hydrography data, should match basemaps source
+river_upa       = 30            # minimum upstream area threshold for the river map [km2]
+slope_len       = 2000          # length over which tp calculate river slope [m]
+
+[setup_reservoirs]
+reservoirs_fn   = hydro_reservoirs  # source for reservoirs based on GRAND: {hydro_reservoirs}; None to skip
+min_area        = 1.0           # minimum lake area to consider [km2]
+priority_jrc    = True          # if True then JRC data from hydroengine is used to calculate some reservoir attributes instead of the GRanD and HydroLAKES db.   
+
+[setup_lakes]
+lakes_fn        = hydro_lakes   # source for lakes based on hydroLAKES: {hydro_lakes}; None to skip
+min_area        = 10.0          # minimum reservoir area to consider [km2]
+
+[setup_glaciers]
+glaciers_fn     = rgi           # source for glaciers based on Randolph Glacier Inventory {rgi}; None to skip
+min_area        = 1.0           # minimum glacier area to consider [km2]
+
+[setup_riverwidth]
+precip_fn       = chelsa        # source for precip climatology used to estimate discharge: {chelsa}
+climate_fn      = koppen_geiger # source for climate classification used to estimate discharge: {koppen_geiger}
+predictor       = discharge     # predictor used in power-law w=a*predictor^b {'discharge'; 'uparea', other staticmaps}; a and b can also be set here.
+fill            = False         # if False all river width values are set based on predictor, if True only data gaps and lakes/reservoirs in observed width are filled (works only with MERIT hydro)
+min_wth         = 1             # global minimum width
+
+[setup_lulcmaps]
+lulc_fn         = globcover     # source for lulc maps: {globcover, vito, corine}
+
+[setup_laimaps]
+lai_fn          = modis_lai     # source for LAI: {modis_lai}
+
+[setup_soilmaps]
+soil_fn         = soilgrids     # source for soilmaps: {soilgrids}
+ptf_ksatver     = brakensiek    # pedotransfer function to calculate hydraulic conductivity: {brakensiek, cosby}
+
+[setup_surfacearea]
+
+[setup_aquiferthickness]
+aquiferthickness_fn = soilgrids_depthbedrock 	# source for aquifer thickness: {soilgrids} absolute depth to bedrock
+
+[setup_drains]
+drain_resistance = 2 			# [days]
+drain_depth = 0.2				# [m] below the surface
+
+[setup_riverbottomelev]
+r = 0.12 			# param to determine river depth based on river width
+p = 0.78				# param to determine river depth based on river width 
+
+[setup_riverconductance]
+
+[setup_gauges]
+gauges_fn       = grdc          # if not None add gaugemap. Either a path or known gauges_fn: {grdc}
+snap_to_river   = True          # if True snaps gauges from source to river
+derive_subcatch = False         # if True derive subcatch map based on gauges.
+
+[setup_precip_forcing]
+precip_fn       = era5          # source for precipitation.
+precip_clim_fn  = None          # source for high resolution climatology to correct precipitation.
+
+[setup_temp_pet_forcing]
+temp_pet_fn     = era5          # source for temperature and potential evapotranspiration.
+press_correction= True          # if True temperature is corrected with elevation lapse rate.
+temp_correction = True          # if True pressure is corrected with elevation lapse rate.
+dem_forcing_fn  = era5_orography # source of elevation grid corresponding to temp_pet_fn. Used for lapse rate correction.
+pet_method      = debruin       # method to compute PET: {debruin, makkink}
+skip_pet        = False         # if True, only temperature is prepared.
+
+[setup_constant_pars]
+KsatHorFrac=100
+Cfmax = 3.75653
+cf_soil = 0.038
+EoverR = 0.11
+InfiltCapPath = 5
+InfiltCapSoil = 600 
+MaxLeakage = 0
+rootdistpar = -500
+TT = 0
+TTI = 2
+TTM = 0
+WHC = 0.1
+G_Cfmax = 5.3
+G_SIfrac = 0.002
+G_TT = 1.3
+specific_yield = 0.15
\ No newline at end of file
diff --git a/hydromt_wflow/utils.py b/hydromt_wflow/utils.py
index cd840a45..dce4e2ab 100644
--- a/hydromt_wflow/utils.py
+++ b/hydromt_wflow/utils.py
@@ -2,12 +2,15 @@
 import xarray as xr
 from typing import Dict, Union
 from pathlib import Path
+import logging
 
 from hydromt.io import open_timeseries_from_table
 from hydromt.vector import GeoDataArray
+from hydromt.gis_utils import cellarea
 
+__all__ = ["read_csv_results", "surfacearea"]
 
-__all__ = ["read_csv_results"]
+logger = logging.getLogger(__name__)
 
 
 def read_csv_results(fn: Union[str, Path], config: Dict, maps: xr.Dataset) -> Dict:
@@ -167,3 +170,20 @@ def read_csv_results(fn: Union[str, Path], config: Dict, maps: xr.Dataset) -> Di
         csv_dict[f"{da.name}"] = da
 
     return csv_dict
+
+
+def surfacearea(ds_like, logger=logger):
+    """ """
+    ds_out = xr.Dataset(coords=ds_like.raster.coords)
+
+    ys, xs = ds_like.raster.ycoords.values, ds_like.raster.xcoords.values
+    xres = np.abs(np.mean(np.diff(xs)))
+    yres = np.abs(np.mean(np.diff(ys)))
+    area = np.ones((ys.size, xs.size), dtype=ys.dtype)
+    cell_res = cellarea(ys, xres, yres)[:, None] * area
+    attrs = dict(_FillValue=-9999, unit="m2")
+    dims = ds_like.raster.dims
+    cell_res = xr.Variable(dims, cell_res, attrs=attrs)
+    ds_out["wflow_surfacearea"] = cell_res
+
+    return ds_out
diff --git a/hydromt_wflow/wflow.py b/hydromt_wflow/wflow.py
index c9a85ab2..58831b69 100644
--- a/hydromt_wflow/wflow.py
+++ b/hydromt_wflow/wflow.py
@@ -29,6 +29,7 @@
     waterbodymaps,
     reservoirattrs,
     soilgrids,
+    aquifer_thickness,
     glaciermaps,
 )
 from .workflows import landuse, lai
@@ -86,6 +87,10 @@ class WflowModel(Model):
         "glacareas": "wflow_glacierareas",
         "glacfracs": "wflow_glacierfrac",
         "glacstore": "wflow_glacierstore",
+        "cellarea": "wflow_surfacearea",
+        "depth_bedrock": "SoilThickness_bedrock",
+        "conductivity_average": "conductivity_average_md",
+        "conductivity_surface": "conductivity_surface_md",
     }
     _FOLDERS = [
         "staticgeoms",
@@ -1154,6 +1159,223 @@ def setup_glaciers(self, glaciers_fn="rgi", min_area=1):
                 f"Skipping glacier procedures!"
             )
 
+    def setup_surfacearea(self):
+        """
+        Calculate cell area at model resolution [m2]
+        TODO: this actually better fits in setup_basemaps as it may be useful for other models? therefore remove function from utils?
+
+        Adds model layer:
+        * **wflow_surfacearea**: surface area of cells at model resolution [m2]
+
+        """
+        self.logger.info(f"Preparing wflow_surfacearea map.")
+        dsout = utils.surfacearea(
+            self.staticmaps,
+            logger=self.logger,
+        )
+        self.set_staticmaps(dsout)
+
+    def setup_aquiferthickness(self, aquiferthickness_fn="soilgrids_depthbedrock"):
+        """
+        Setup aquifer thickness map (for sbm+gw)
+
+        Adds model layer:
+
+        * **SoilThickness_bedrock**: aquifer thickness [mm]
+
+        Parameters
+        ----------
+        aquiferthickness_fn: data source (default: soilgrids absolute depth to bedrock) [mm]
+
+        """
+        self.logger.info(f"Preparing aquifer thickness parameter map.")
+
+        dsin = self.data_catalog.get_rasterdataset(
+            aquiferthickness_fn, geom=self.region, buffer=2
+        )
+
+        dsout = aquifer_thickness(
+            dsin,
+            self.staticmaps,
+            aquiferthickness_fn,
+            logger=self.logger,
+        )
+        self.set_staticmaps(dsout)
+
+    def setup_drains(self, drain_resistance=2.0, drain_depth=0.2):
+        """
+        Setup drain maps for the sbm+groundwater model.
+        The drain conductance is based on a resistance [d]
+        and the cell area (at model resolution).
+        The drain elevation is defined from the drain depth [m] below the DEM [m].
+
+        Adds model layer:
+
+        * **drain_elevation**: elevation of the drains [m]
+        * **drain_conductance**:  [m2/day]
+
+        Parameters
+        ----------
+        drain_resistance = 2 [d] (default)
+        drain_depth = 0.2 [m] (default)
+
+        """
+        self.logger.info(f"Preparing drain maps.")
+
+        if not self._MAPS["cellarea"] in self.staticmaps:
+            raise ValueError(
+                "The setup_drains method requires to run setup_surfacearea method first."
+            )
+
+        nodatafloat = -999
+
+        drain_elevation = xr.where(
+            self.staticmaps[self._MAPS["basins"]],
+            self.staticmaps[self._MAPS["elevtn"]] - drain_depth,
+            nodatafloat,
+        )
+        drain_elevation.raster.set_nodata(nodatafloat)
+        drain_elevation = drain_elevation.rename("drain_elevation")
+        self.set_staticmaps(drain_elevation)
+
+        drain_conductance = xr.where(
+            self.staticmaps[self._MAPS["basins"]],
+            self.staticmaps[self._MAPS["cellarea"]] / drain_resistance,
+            nodatafloat,
+        )
+        drain_conductance = drain_conductance.rename("drain_conductance")
+        self.set_staticmaps(drain_conductance)
+
+    def setup_riverbottomelev(self, r=0.12, p=0.78):
+        """
+        Estimate river bottom elevation for the sbm+groundwater model.
+        This is based on a very simple formula with coefficients for gravel bed rivers in Britain. (Neal et al., 2012)
+        TODO: check if better estimates are available.
+
+        Adds model layer:
+
+        * **river_bottom_elevation**: elevation of the bottom of the river [m]
+
+        Parameters
+        ----------
+        r = 0.12 coefficient [-] (default)
+        p = 0.78 exponent [-] (default)
+        """
+
+        if not self._MAPS["rivwth"] in self.staticmaps:
+            raise ValueError(
+                "The setup_riverbottomelev method requies to run setup_riverwidth method first."
+            )
+
+        nodatafloat = -999
+        river_depth = xr.where(
+            self.staticmaps[self._MAPS["rivwth"]],
+            r * self.staticmaps[self._MAPS["rivwth"]] ** p,
+            nodatafloat,
+        )
+        river_bottom_elevation = self.staticmaps[self._MAPS["elevtn"]] - river_depth
+        river_bottom_elevation = river_bottom_elevation.rename("river_bottom_elevation")
+        self.set_staticmaps(river_bottom_elevation)
+
+        # warning if depth of river below depth to bedrock
+        if self._MAPS["depth_bedrock"] in self.staticmaps:
+            diff_depthbedrock_depthriver = xr.where(
+                self.staticmaps[self._MAPS["rivwth"]],
+                self.staticmaps[self._MAPS["depth_bedrock"]] / 1000 - river_depth,
+                nodatafloat,
+            )
+            diff_depthbedrock_depthriver = diff_depthbedrock_depthriver.rename(
+                "diff_depthbedrock_depthriver"
+            )
+            # TODO: check add this map to staticmaps? probably not, warning enough.
+            self.set_staticmaps(diff_depthbedrock_depthriver)
+            if np.min(diff_depthbedrock_depthriver) < 0:
+                self.logger.warning(
+                    "River bottom elevation lower than absolute depth to bedrock"
+                )
+
+    def setup_riverconductance(self):
+        """
+        Estimate river infiltration and exfiltration conductance based on surface and average conductivity (KsatVer).
+
+        """
+
+        nodatafloat = -999
+
+        if not self._MAPS["rivwth"] in self.staticmaps:
+            raise ValueError(
+                "The setup_riverconductance method requires to run setup_riverwidth method first."
+            )
+        if not self._MAPS["cellarea"] in self.staticmaps:
+            raise ValueError(
+                "The setup_riverconductance method requires to run setup_surfacearea method first."
+            )
+        if not self._MAPS["depth_bedrock"] in self.staticmaps:
+            raise ValueError(
+                "The setup_riverconductance method requires to run setup_aquiferthickness method first."
+            )
+        if not self._MAPS["conductivity_surface"] in self.staticmaps:
+            raise ValueError(
+                "The setup_riverconductance method requires to run setup_soil method first."
+            )
+
+        self.logger.info(
+            f"Preparing river infiltration and exfiltration conductance maps."
+        )
+
+        # Compute average distance between ditches
+        wetted_area = (
+            self.staticmaps[self._MAPS["rivlen"]]
+            * self.staticmaps[self._MAPS["rivwth"]]
+        )
+        # TODO: wetted area based on high res data?
+
+        # check if cell_area larger than wetted area or emit warning
+        if np.min(wetted_area) < 0:
+            self.logger.warning("Wetted area in cell larger than cell area!")
+
+        L = (self.staticmaps[self._MAPS["cellarea"]] - wetted_area) / self.staticmaps[
+            self._MAPS["rivlen"]
+        ]
+
+        # compute "lumped" resistance for each cel (based on surface conductivity and average conductivity over soil profile)
+        # TODO: NB: conductivity in [m/d]
+        horizontal_flow_resistance_surf = L ** 2 / (
+            8
+            * self.staticmaps[self._MAPS["conductivity_surface"]]
+            * self.staticmaps[self._MAPS["depth_bedrock"]]
+            / 1000
+        )  # depth_bedrock from mm to m
+        horizontal_flow_resistance_aver = L ** 2 / (
+            8
+            * self.staticmaps[self._MAPS["conductivity_average"]]
+            * self.staticmaps[self._MAPS["depth_bedrock"]]
+            / 1000
+        )  # depth_bedrock from mm to m
+
+        # compute conductance
+        conductance_surf = (
+            self.staticmaps[self._MAPS["cellarea"]] / horizontal_flow_resistance_surf
+        )
+        conductance_aver = (
+            self.staticmaps[self._MAPS["cellarea"]] / horizontal_flow_resistance_aver
+        )
+
+        # mask except river cells
+        conductance_surf = xr.where(
+            self.staticmaps[self._MAPS["rivmsk"]], conductance_surf, nodatafloat
+        )
+        conductance_aver = xr.where(
+            self.staticmaps[self._MAPS["rivmsk"]], conductance_aver, nodatafloat
+        )
+        conductance_surf.raster.set_nodata(nodatafloat)
+        conductance_aver.raster.set_nodata(nodatafloat)
+
+        conductance_surf = conductance_surf.rename("infiltration_conductance_Ksurface")
+        conductance_aver = conductance_aver.rename("infiltration_conductance_Kaverage")
+        self.set_staticmaps(conductance_surf)
+        self.set_staticmaps(conductance_aver)
+
     def setup_constant_pars(self, **kwargs):
         """Setup constant parameter maps.
 
diff --git a/hydromt_wflow/workflows/soilgrids.py b/hydromt_wflow/workflows/soilgrids.py
index 16dd8362..4868b77a 100644
--- a/hydromt_wflow/workflows/soilgrids.py
+++ b/hydromt_wflow/workflows/soilgrids.py
@@ -15,7 +15,7 @@
 
 logger = logging.getLogger(__name__)
 
-__all__ = ["soilgrids", "soilgrids_sediment"]
+__all__ = ["soilgrids", "soilgrids_sediment", "aquifer_thickness"]
 
 # soilgrids_2017
 soildepth_cm = np.array([0.0, 5.0, 15.0, 30.0, 60.0, 100.0, 200.0])
@@ -425,8 +425,12 @@ def soilgrids(ds, ds_like, ptfKsatVer, soil_fn, logger=logger):
     logger.info("calculate and resample KsatVer")
     kv_sl_hr = kv_layers(ds, thetas_sl, ptfKsatVer)
     kv_sl = np.log(kv_sl_hr)
+    # integrated Ksat over depth for sbm+groundwater
+    kv_sl_integrated = average_soillayers(kv_sl, ds["soilthickness"])
     kv_sl = kv_sl.raster.reproject_like(ds_like, method="average")
     kv_sl = np.exp(kv_sl)
+    kv_sl_integrated = kv_sl_integrated.raster.reproject_like(ds_like, method="average")
+    kv_sl_integrated = np.exp(kv_sl_integrated)
 
     logger.info("calculate and resample pore size distribution index")
     lambda_sl_hr = pore_size_distrution_index_layers(ds, thetas_sl)
@@ -466,6 +470,12 @@ def soilgrids(ds, ds_like, ptfKsatVer, soil_fn, logger=logger):
     ds_out["c"] = ds_c
 
     ds_out["KsatVer"] = kv_sl.sel(sl=1).astype(np.float32)
+    ds_out["conductivity_average_md"] = (kv_sl_integrated / 1000).astype(
+        np.float32
+    )  # [m/d] conductivity averaged over soil profile
+    ds_out["conductivity_surface_md"] = (kv_sl.sel(sl=1) / 1000).astype(
+        np.float32
+    )  # [m/d] conductivity at soil surface -- unit for wflow_sbm+gw
 
     for i, sl in enumerate(kv_sl["sl"]):
         kv = kv_sl.sel(sl=sl)
@@ -543,6 +553,53 @@ def soilgrids(ds, ds_like, ptfKsatVer, soil_fn, logger=logger):
     return ds_out
 
 
+def aquifer_thickness(ds, ds_like, aquiferthickness_fn, logger=logger):
+    """
+    Returns an estimate of the aquifer thickness (wflow sbm + gw) [mm]
+
+    Adds model layer:
+
+        * **SoilThickness_bedrock**: aquifer thickness [mm]
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        Dataset containing soil properties.
+    ds_like : xarray.DataArray
+        Dataset at model resolution.
+    aquiferthickness_fn : str
+        soilgrids v1.0 {'soilgrids_depthbedrock'}
+
+    Returns
+    -------
+    ds_out : xarray.Dataset
+        Dataset containing gridded soil parameters for aquifer thickness.
+
+    """
+
+    ds_out = xr.Dataset(coords=ds_like.raster.coords)
+    ds = ds.raster.mask_nodata()
+
+    # absolute depth to bedrock for sbm+gw
+    soilthickness_bedrock_hr = ds.raster.interpolate_na("nearest")
+    soilthickness_bedrock = soilthickness_bedrock_hr.raster.reproject_like(
+        ds_like, method="average"
+    )
+    # wflow_sbm cannot handle (yet) zero soil thickness
+    soilthickness_bedrock = soilthickness_bedrock.where(
+        soilthickness_bedrock > 0.0, np.nan
+    )
+    soilthickness_bedrock.raster.set_nodata(np.nan)
+    soilthickness_bedrock = soilthickness_bedrock.raster.interpolate_na(
+        "nearest"
+    ).astype(np.float32)
+    ds_out[
+        "SoilThickness_bedrock"
+    ] = soilthickness_bedrock  # [mm] conversion already done in yml.
+
+    return ds_out
+
+
 def soilgrids_sediment(ds, ds_like, usleK_method, logger=logger):
 
     """