diff --git a/build/create_binaries/download_test_data.jl b/build/create_binaries/download_test_data.jl
index 2c1fe00fa..694f128e5 100644
--- a/build/create_binaries/download_test_data.jl
+++ b/build/create_binaries/download_test_data.jl
@@ -46,6 +46,8 @@ forcing_sbm_gw_path = testdata(
 forcing_meuse_path = testdata(v"0.2.8", "forcing_meuse.nc", "forcing_meuse.nc")
 staticmaps_sbm_gw_path =
     testdata(v"0.2.2", "staticmaps-sbm-groundwater.nc", "staticmaps-sbm-groundwater.nc")
+instates_sbm_gw_path =
+    testdata(v"0.2.2", "instates-example-sbm-gwf.nc", "instates-example-sbm-gwf.nc")
 lake_sh_1_path = testdata(v"0.2.1", "lake_sh_1.csv", "lake_sh_1.csv")
 lake_sh_2_path = testdata(v"0.2.1", "lake_sh_2.csv", "lake_sh_2.csv")
 lake_hq_2_path = testdata(v"0.2.1", "lake_hq_2.csv", "lake_hq_2.csv")
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index cd78cf337..27ab98977 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -28,6 +28,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   model parameters are defined, and it is not expected that exposing these variables is
   required (e.g. for model coupling) while code changes for these variables (including
   struct fields) are required.
+- The local inertial routing (constant) river boundary condition `depth` at a ghost node
+  (downstream river outlet) was set at 0.0 in the code and can now be provided through the
+  model parameter netCDF file (`riverdepth_bc`). In addition, the boundary condition river
+  length `dl` at a ghost node should be set through the model parameter netCDF file
+  (`riverlength_bc`), providing this boundary condition through the `[model]` settings in
+  the TOML file (`model.riverlength_bc`) has been deprecated.
 
 ### Added
 - Total water storage as an export variable for `SBM` concept. This is the total water stored
@@ -39,6 +45,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Support different vertical hydraulic conductivity profiles for the `SBM` concept. See also
   the following sections: [The SBM soil water accounting scheme](@ref) and [Subsurface flow
   routing](@ref) for a short description.
+- Local inertial routing to `sbm_gwf` model type.
 
 ## v0.7.3 - 2024-01-12
 
diff --git a/docs/src/model_docs/lateral/local-inertial.md b/docs/src/model_docs/lateral/local-inertial.md
index 1595a2965..12e6ce176 100644
--- a/docs/src/model_docs/lateral/local-inertial.md
+++ b/docs/src/model_docs/lateral/local-inertial.md
@@ -37,7 +37,7 @@ river length [m] for river cell ``i`` and ``\alpha`` is a coefficient (typically
 and 0.7) to enhance the stability of the simulation.
 
 In the TOML file the following properties related to the local inertial model can be
-provided for the `sbm` model type:
+provided for the `sbm` and `sbm_gwf` model types:
 
 ```toml
 [model]
@@ -45,9 +45,18 @@ river_routing = "local-inertial"  # default is "kinematic-wave"
 inertial_flow_alpha = 0.5         # alpha coefficient for model stability (default = 0.7)
 froude_limit = true               # default is true, limit flow to subcritical-critical according to Froude number
 h_thresh = 0.1                    # water depth [m] threshold for calculating flow between cells (default = 1e-03)
-riverlength_bc = 1000.0           # river length [m] for boundary points (default = 1e04)
 floodplain_1d = true              # include 1D floodplain schematization (default = false)
 ```
+Two optional constant boundary conditions `riverlength_bc` and `riverdepth_bc` can be
+provided at a river outlet node (or multiple river outlet nodes) through the model parameter
+netCDF file, as follows:
+```toml
+[input.lateral.river]
+riverlength_bc = "riverlength_bc"   # optional river length [m], default = 1e04
+riverdepth_bc = "riverdepth_bc"     # optional river depth [m], default = 0.0
+```
+These boundary conditions are copied to a ghost node (downstream of the river outlet node)
+in the code.
 
 The optional 1D floodplain schematization is based on provided flood volumes as a function
 of flood depth (per flood depth interval) for each river cell. Wflow calculates from these
diff --git a/docs/src/model_docs/model_configurations.md b/docs/src/model_docs/model_configurations.md
index f190cb428..90a3ce8bc 100644
--- a/docs/src/model_docs/model_configurations.md
+++ b/docs/src/model_docs/model_configurations.md
@@ -45,48 +45,6 @@ lateral.river.lake => struct NaturalLake{T} # optional
 lateral.river.reservoir => struct SimpleReservoir{T} # optional
 ```
 
-### SBM + Local inertial river and floodplain
-By default the model type `sbm` uses the kinematic wave approach for river flow. There is
-also the option to use the local inertial model for river flow with an optional 1D
-floodplain schematization (routing is done separately for the river channel and floodplain),
-by providing the following in the TOML file:
-
-```toml
-[model]
-river_routing = "local-inertial"    # optional, default is "kinematic-wave"
-floodplain_1d = true                # optional, default is false
-```
-
-Only the mapping for the river component changes, as shown below. For an explanation about
-the type parameters between curly braces after the `struct` name see the section on the model
-parameters.
-
-```julia
-lateral.river => struct ShallowWaterRiver{T,R,L}
-```
-
-### SBM + Local inertial river (1D) and land (2D)
-By default the model type `sbm` uses the kinematic wave approach for river and overland
-flow. There is also the option to use the local inertial model for 1D river and 2D overland
-flow, by providing the following in the TOML file:
-
-```toml
-[model]
-river_routing = "local-inertial"
-land_routing = "local-inertial"
-```
-The mapping for the river and land component changes, as shown below. For an explanation
-about the type parameters between curly braces after the `struct` name see the section on
-the model parameters.
-
-```julia
-lateral.river => struct ShallowWaterRiver{T,R,L}
-lateral.land => struct ShallowWaterLand{T}
-```
-
-The local inertial approach is described in more detail in the section [Local inertial
-model](@ref local_inertial).
-
 ### SBM + Groundwater flow
 For river and overland flow the kinematic wave approach over a D8 network is used for this
 wflow\_sbm model. For the subsurface domain, an unconfined aquifer with groundwater flow in
@@ -123,6 +81,48 @@ lateral.river.lake => struct NaturalLake{T} # optional
 lateral.river.reservoir => struct SimpleReservoir{T} # optional
 ```
 
+### Local inertial river and floodplain + `sbm` and `sbm_gwf` model types
+By default the model types `sbm` and `sbm_gwf` uses the kinematic wave approach for river
+flow. There is also the option to use the local inertial model for river flow with an
+optional 1D floodplain schematization (routing is done separately for the river channel and
+floodplain), by providing the following in the TOML file:
+
+```toml
+[model]
+river_routing = "local-inertial"    # optional, default is "kinematic-wave"
+floodplain_1d = true                # optional, default is false
+```
+
+Only the mapping for the river component changes, as shown below. For an explanation about
+the type parameters between curly braces after the `struct` name see the section on the model
+parameters.
+
+```julia
+lateral.river => struct ShallowWaterRiver{T,R,L}
+```
+
+### Local inertial river (1D) and land (2D) + `sbm` and `sbm_gwf` model types
+By default the model types `sbm` and `sbm_gwf` uses the kinematic wave approach for river
+and overland flow. There is also the option to use the local inertial model for 1D river and
+2D overland flow, by providing the following in the TOML file:
+
+```toml
+[model]
+river_routing = "local-inertial"
+land_routing = "local-inertial"
+```
+The mapping for the river and land component changes, as shown below. For an explanation
+about the type parameters between curly braces after the `struct` name see the section on
+the model parameters.
+
+```julia
+lateral.river => struct ShallowWaterRiver{T,R,L}
+lateral.land => struct ShallowWaterLand{T}
+```
+
+The local inertial approach is described in more detail in the section [Local inertial
+model](@ref local_inertial).
+
 ## [wflow\_hbv](@id config_hbv)
 The Hydrologiska Byrans Vattenbalansavdelning (HBV) model was introduced back in 1972 by the
 Swedisch Meteological and Hydrological Institute (SMHI).  The HBV model is mainly used for
diff --git a/pixi.lock b/pixi.lock
index 2c729d6a0..c641b5934 100644
--- a/pixi.lock
+++ b/pixi.lock
@@ -7,14 +7,43 @@ environments:
       linux-64:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_libgcc_mutex-0.1-conda_forge.tar.bz2
       - conda: https://conda.anaconda.org/conda-forge/linux-64/_openmp_mutex-4.5-2_gnu.tar.bz2
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hd590300_5.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/ca-certificates-2024.2.2-hbcca054_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/juliaup-1.13.0-he8a937b_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/ld_impl_linux-64-2.40-h41732ed_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libexpat-2.6.2-h59595ed_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libffi-3.4.2-h7f98852_5.tar.bz2
       - conda: https://conda.anaconda.org/conda-forge/linux-64/libgcc-ng-13.2.0-h807b86a_5.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/libgomp-13.2.0-h807b86a_5.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libnsl-2.0.1-hd590300_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.45.2-h2797004_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libuuid-2.38.1-h0b41bf4_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libxcrypt-4.4.36-hd590300_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.2.13-hd590300_5.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/ncurses-6.4.20240210-h59595ed_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/openssl-3.2.1-hd590300_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/python-3.12.2-hab00c5b_0_cpython.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/readline-8.2-h8228510_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/tk-8.6.13-noxft_h4845f30_101.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/tzdata-2024a-h0c530f3_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/linux-64/xz-5.2.6-h166bdaf_0.tar.bz2
       win-64:
+      - conda: https://conda.anaconda.org/conda-forge/win-64/bzip2-1.0.8-hcfcfb64_5.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/ca-certificates-2024.2.2-h56e8100_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/juliaup-1.13.0-h975169c_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/libexpat-2.6.2-h63175ca_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/libffi-3.4.2-h8ffe710_5.tar.bz2
+      - conda: https://conda.anaconda.org/conda-forge/win-64/libsqlite-3.45.2-hcfcfb64_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/libzlib-1.2.13-hcfcfb64_5.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/openssl-3.2.1-hcfcfb64_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/python-3.12.2-h2628c8c_0_cpython.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/tk-8.6.13-h5226925_1.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/tzdata-2024a-h0c530f3_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/ucrt-10.0.22621.0-h57928b3_0.tar.bz2
       - conda: https://conda.anaconda.org/conda-forge/win-64/vc-14.3-hcf57466_18.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/vc14_runtime-14.38.33130-h82b7239_18.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/vs2015_runtime-14.38.33130-hcb4865c_18.conda
+      - conda: https://conda.anaconda.org/conda-forge/win-64/xz-5.2.6-h8d14728_0.tar.bz2
 packages:
 - kind: conda
   name: _libgcc_mutex
@@ -45,6 +74,60 @@ packages:
   license_family: BSD
   size: 23621
   timestamp: 1650670423406
+- kind: conda
+  name: bzip2
+  version: 1.0.8
+  build: hcfcfb64_5
+  build_number: 5
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/bzip2-1.0.8-hcfcfb64_5.conda
+  sha256: ae5f47a5c86fd6db822931255dcf017eb12f60c77f07dc782ccb477f7808aab2
+  md5: 26eb8ca6ea332b675e11704cce84a3be
+  depends:
+  - ucrt >=10.0.20348.0
+  - vc >=14.2,<15
+  - vc14_runtime >=14.29.30139
+  license: bzip2-1.0.6
+  license_family: BSD
+  size: 124580
+  timestamp: 1699280668742
+- kind: conda
+  name: bzip2
+  version: 1.0.8
+  build: hd590300_5
+  build_number: 5
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/bzip2-1.0.8-hd590300_5.conda
+  sha256: 242c0c324507ee172c0e0dd2045814e746bb303d1eb78870d182ceb0abc726a8
+  md5: 69b8b6202a07720f448be700e300ccf4
+  depends:
+  - libgcc-ng >=12
+  license: bzip2-1.0.6
+  license_family: BSD
+  size: 254228
+  timestamp: 1699279927352
+- kind: conda
+  name: ca-certificates
+  version: 2024.2.2
+  build: h56e8100_0
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/ca-certificates-2024.2.2-h56e8100_0.conda
+  sha256: 4d587088ecccd393fec3420b64f1af4ee1a0e6897a45cfd5ef38055322cea5d0
+  md5: 63da060240ab8087b60d1357051ea7d6
+  license: ISC
+  size: 155886
+  timestamp: 1706843918052
+- kind: conda
+  name: ca-certificates
+  version: 2024.2.2
+  build: hbcca054_0
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/ca-certificates-2024.2.2-hbcca054_0.conda
+  sha256: 91d81bfecdbb142c15066df70cc952590ae8991670198f92c66b62019b251aeb
+  md5: 2f4327a1cbe7f022401b236e915a5fef
+  license: ISC
+  size: 155432
+  timestamp: 1706843687645
 - kind: conda
   name: juliaup
   version: 1.13.0
@@ -75,6 +158,81 @@ packages:
   license_family: MIT
   size: 3859930
   timestamp: 1706564340550
+- kind: conda
+  name: ld_impl_linux-64
+  version: '2.40'
+  build: h41732ed_0
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/ld_impl_linux-64-2.40-h41732ed_0.conda
+  sha256: f6cc89d887555912d6c61b295d398cff9ec982a3417d38025c45d5dd9b9e79cd
+  md5: 7aca3059a1729aa76c597603f10b0dd3
+  constrains:
+  - binutils_impl_linux-64 2.40
+  license: GPL-3.0-only
+  license_family: GPL
+  size: 704696
+  timestamp: 1674833944779
+- kind: conda
+  name: libexpat
+  version: 2.6.2
+  build: h59595ed_0
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/libexpat-2.6.2-h59595ed_0.conda
+  sha256: 331bb7c7c05025343ebd79f86ae612b9e1e74d2687b8f3179faec234f986ce19
+  md5: e7ba12deb7020dd080c6c70e7b6f6a3d
+  depends:
+  - libgcc-ng >=12
+  constrains:
+  - expat 2.6.2.*
+  license: MIT
+  license_family: MIT
+  size: 73730
+  timestamp: 1710362120304
+- kind: conda
+  name: libexpat
+  version: 2.6.2
+  build: h63175ca_0
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/libexpat-2.6.2-h63175ca_0.conda
+  sha256: 79f612f75108f3e16bbdc127d4885bb74729cf66a8702fca0373dad89d40c4b7
+  md5: bc592d03f62779511d392c175dcece64
+  constrains:
+  - expat 2.6.2.*
+  license: MIT
+  license_family: MIT
+  size: 139224
+  timestamp: 1710362609641
+- kind: conda
+  name: libffi
+  version: 3.4.2
+  build: h7f98852_5
+  build_number: 5
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/libffi-3.4.2-h7f98852_5.tar.bz2
+  sha256: ab6e9856c21709b7b517e940ae7028ae0737546122f83c2aa5d692860c3b149e
+  md5: d645c6d2ac96843a2bfaccd2d62b3ac3
+  depends:
+  - libgcc-ng >=9.4.0
+  license: MIT
+  license_family: MIT
+  size: 58292
+  timestamp: 1636488182923
+- kind: conda
+  name: libffi
+  version: 3.4.2
+  build: h8ffe710_5
+  build_number: 5
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/libffi-3.4.2-h8ffe710_5.tar.bz2
+  sha256: 1951ab740f80660e9bc07d2ed3aefb874d78c107264fd810f24a1a6211d4b1a5
+  md5: 2c96d1b6915b408893f9472569dee135
+  depends:
+  - vc >=14.1,<15.0a0
+  - vs2015_runtime >=14.16.27012
+  license: MIT
+  license_family: MIT
+  size: 42063
+  timestamp: 1636489106777
 - kind: conda
   name: libgcc-ng
   version: 13.2.0
@@ -108,6 +266,281 @@ packages:
   license_family: GPL
   size: 419751
   timestamp: 1706819107383
+- kind: conda
+  name: libnsl
+  version: 2.0.1
+  build: hd590300_0
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/libnsl-2.0.1-hd590300_0.conda
+  sha256: 26d77a3bb4dceeedc2a41bd688564fe71bf2d149fdcf117049970bc02ff1add6
+  md5: 30fd6e37fe21f86f4bd26d6ee73eeec7
+  depends:
+  - libgcc-ng >=12
+  license: LGPL-2.1-only
+  license_family: GPL
+  size: 33408
+  timestamp: 1697359010159
+- kind: conda
+  name: libsqlite
+  version: 3.45.2
+  build: h2797004_0
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/libsqlite-3.45.2-h2797004_0.conda
+  sha256: 8cdbeb7902729e319510a82d7c642402981818702b58812af265ef55d1315473
+  md5: 866983a220e27a80cb75e85cb30466a1
+  depends:
+  - libgcc-ng >=12
+  - libzlib >=1.2.13,<1.3.0a0
+  license: Unlicense
+  size: 857489
+  timestamp: 1710254744982
+- kind: conda
+  name: libsqlite
+  version: 3.45.2
+  build: hcfcfb64_0
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/libsqlite-3.45.2-hcfcfb64_0.conda
+  sha256: 4bb24b986550275a6d02835150d943c4c675808d05c0efc5c2a22154d007a69f
+  md5: f95359f8dc5abf7da7776ece9ef10bc5
+  depends:
+  - ucrt >=10.0.20348.0
+  - vc >=14.2,<15
+  - vc14_runtime >=14.29.30139
+  license: Unlicense
+  size: 869606
+  timestamp: 1710255095740
+- kind: conda
+  name: libuuid
+  version: 2.38.1
+  build: h0b41bf4_0
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/libuuid-2.38.1-h0b41bf4_0.conda
+  sha256: 787eb542f055a2b3de553614b25f09eefb0a0931b0c87dbcce6efdfd92f04f18
+  md5: 40b61aab5c7ba9ff276c41cfffe6b80b
+  depends:
+  - libgcc-ng >=12
+  license: BSD-3-Clause
+  license_family: BSD
+  size: 33601
+  timestamp: 1680112270483
+- kind: conda
+  name: libxcrypt
+  version: 4.4.36
+  build: hd590300_1
+  build_number: 1
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/libxcrypt-4.4.36-hd590300_1.conda
+  sha256: 6ae68e0b86423ef188196fff6207ed0c8195dd84273cb5623b85aa08033a410c
+  md5: 5aa797f8787fe7a17d1b0821485b5adc
+  depends:
+  - libgcc-ng >=12
+  license: LGPL-2.1-or-later
+  size: 100393
+  timestamp: 1702724383534
+- kind: conda
+  name: libzlib
+  version: 1.2.13
+  build: hcfcfb64_5
+  build_number: 5
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/libzlib-1.2.13-hcfcfb64_5.conda
+  sha256: c161822ee8130b71e08b6d282b9919c1de2c5274b29921a867bca0f7d30cad26
+  md5: 5fdb9c6a113b6b6cb5e517fd972d5f41
+  depends:
+  - ucrt >=10.0.20348.0
+  - vc >=14.2,<15
+  - vc14_runtime >=14.29.30139
+  constrains:
+  - zlib 1.2.13 *_5
+  license: Zlib
+  license_family: Other
+  size: 55800
+  timestamp: 1686575452215
+- kind: conda
+  name: libzlib
+  version: 1.2.13
+  build: hd590300_5
+  build_number: 5
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/libzlib-1.2.13-hd590300_5.conda
+  sha256: 370c7c5893b737596fd6ca0d9190c9715d89d888b8c88537ae1ef168c25e82e4
+  md5: f36c115f1ee199da648e0597ec2047ad
+  depends:
+  - libgcc-ng >=12
+  constrains:
+  - zlib 1.2.13 *_5
+  license: Zlib
+  license_family: Other
+  size: 61588
+  timestamp: 1686575217516
+- kind: conda
+  name: ncurses
+  version: 6.4.20240210
+  build: h59595ed_0
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/ncurses-6.4.20240210-h59595ed_0.conda
+  sha256: aa0f005b6727aac6507317ed490f0904430584fa8ca722657e7f0fb94741de81
+  md5: 97da8860a0da5413c7c98a3b3838a645
+  depends:
+  - libgcc-ng >=12
+  license: X11 AND BSD-3-Clause
+  size: 895669
+  timestamp: 1710866638986
+- kind: conda
+  name: openssl
+  version: 3.2.1
+  build: hcfcfb64_1
+  build_number: 1
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/openssl-3.2.1-hcfcfb64_1.conda
+  sha256: 61ce4e11c3c26ed4e4d9b7e7e2483121a1741ad0f9c8db0a91a28b6e05182ce6
+  md5: 958e0418e93e50c575bff70fbcaa12d8
+  depends:
+  - ca-certificates
+  - ucrt >=10.0.20348.0
+  - vc >=14.2,<15
+  - vc14_runtime >=14.29.30139
+  constrains:
+  - pyopenssl >=22.1
+  license: Apache-2.0
+  license_family: Apache
+  size: 8230112
+  timestamp: 1710796158475
+- kind: conda
+  name: openssl
+  version: 3.2.1
+  build: hd590300_1
+  build_number: 1
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/openssl-3.2.1-hd590300_1.conda
+  sha256: 2c689444ed19a603be457284cf2115ee728a3fafb7527326e96054dee7cdc1a7
+  md5: 9d731343cff6ee2e5a25c4a091bf8e2a
+  depends:
+  - ca-certificates
+  - libgcc-ng >=12
+  constrains:
+  - pyopenssl >=22.1
+  license: Apache-2.0
+  license_family: Apache
+  size: 2865379
+  timestamp: 1710793235846
+- kind: conda
+  name: python
+  version: 3.12.2
+  build: h2628c8c_0_cpython
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/python-3.12.2-h2628c8c_0_cpython.conda
+  sha256: b8eda863b48ae4531635e23fd15e759d93212b6204c6847d591e25fa5fd67477
+  md5: be8803e9f75a477df61d4aabea3c1246
+  depends:
+  - bzip2 >=1.0.8,<2.0a0
+  - libexpat >=2.5.0,<3.0a0
+  - libffi >=3.4,<4.0a0
+  - libsqlite >=3.45.1,<4.0a0
+  - libzlib >=1.2.13,<1.3.0a0
+  - openssl >=3.2.1,<4.0a0
+  - tk >=8.6.13,<8.7.0a0
+  - tzdata
+  - ucrt >=10.0.20348.0
+  - vc >=14.2,<15
+  - vc14_runtime >=14.29.30139
+  - xz >=5.2.6,<6.0a0
+  constrains:
+  - python_abi 3.12.* *_cp312
+  license: Python-2.0
+  size: 16083296
+  timestamp: 1708116662336
+- kind: conda
+  name: python
+  version: 3.12.2
+  build: hab00c5b_0_cpython
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/python-3.12.2-hab00c5b_0_cpython.conda
+  sha256: ddb7a2d8d78046bda5d7631e6814f9468d2eb054e10f86f4648c9d1fdaa30c0f
+  md5: ad7b68400f3a6ebe72b00be093c7f301
+  depends:
+  - bzip2 >=1.0.8,<2.0a0
+  - ld_impl_linux-64 >=2.36.1
+  - libexpat >=2.5.0,<3.0a0
+  - libffi >=3.4,<4.0a0
+  - libgcc-ng >=12
+  - libnsl >=2.0.1,<2.1.0a0
+  - libsqlite >=3.45.1,<4.0a0
+  - libuuid >=2.38.1,<3.0a0
+  - libxcrypt >=4.4.36
+  - libzlib >=1.2.13,<1.3.0a0
+  - ncurses >=6.4,<7.0a0
+  - openssl >=3.2.1,<4.0a0
+  - readline >=8.2,<9.0a0
+  - tk >=8.6.13,<8.7.0a0
+  - tzdata
+  - xz >=5.2.6,<6.0a0
+  constrains:
+  - python_abi 3.12.* *_cp312
+  license: Python-2.0
+  size: 32312631
+  timestamp: 1708118077305
+- kind: conda
+  name: readline
+  version: '8.2'
+  build: h8228510_1
+  build_number: 1
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/readline-8.2-h8228510_1.conda
+  sha256: 5435cf39d039387fbdc977b0a762357ea909a7694d9528ab40f005e9208744d7
+  md5: 47d31b792659ce70f470b5c82fdfb7a4
+  depends:
+  - libgcc-ng >=12
+  - ncurses >=6.3,<7.0a0
+  license: GPL-3.0-only
+  license_family: GPL
+  size: 281456
+  timestamp: 1679532220005
+- kind: conda
+  name: tk
+  version: 8.6.13
+  build: h5226925_1
+  build_number: 1
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/tk-8.6.13-h5226925_1.conda
+  sha256: 2c4e914f521ccb2718946645108c9bd3fc3216ba69aea20c2c3cedbd8db32bb1
+  md5: fc048363eb8f03cd1737600a5d08aafe
+  depends:
+  - ucrt >=10.0.20348.0
+  - vc >=14.2,<15
+  - vc14_runtime >=14.29.30139
+  license: TCL
+  license_family: BSD
+  size: 3503410
+  timestamp: 1699202577803
+- kind: conda
+  name: tk
+  version: 8.6.13
+  build: noxft_h4845f30_101
+  build_number: 101
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/tk-8.6.13-noxft_h4845f30_101.conda
+  sha256: e0569c9caa68bf476bead1bed3d79650bb080b532c64a4af7d8ca286c08dea4e
+  md5: d453b98d9c83e71da0741bb0ff4d76bc
+  depends:
+  - libgcc-ng >=12
+  - libzlib >=1.2.13,<1.3.0a0
+  license: TCL
+  license_family: BSD
+  size: 3318875
+  timestamp: 1699202167581
+- kind: conda
+  name: tzdata
+  version: 2024a
+  build: h0c530f3_0
+  subdir: noarch
+  noarch: generic
+  url: https://conda.anaconda.org/conda-forge/noarch/tzdata-2024a-h0c530f3_0.conda
+  sha256: 7b2b69c54ec62a243eb6fba2391b5e443421608c3ae5dbff938ad33ca8db5122
+  md5: 161081fc7cec0bfda0d86d7cb595f8d8
+  license: LicenseRef-Public-Domain
+  size: 119815
+  timestamp: 1706886945727
 - kind: conda
   name: ucrt
   version: 10.0.22621.0
@@ -156,3 +589,45 @@ packages:
   license_family: Proprietary
   size: 749868
   timestamp: 1702511239004
+- kind: conda
+  name: vs2015_runtime
+  version: 14.38.33130
+  build: hcb4865c_18
+  build_number: 18
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/vs2015_runtime-14.38.33130-hcb4865c_18.conda
+  sha256: a2fec221f361d6263c117f4ea6d772b21c90a2f8edc6f3eb0eadec6bfe8843db
+  md5: 10d42885e3ed84e575b454db30f1aa93
+  depends:
+  - vc14_runtime >=14.38.33130
+  license: BSD-3-Clause
+  license_family: BSD
+  size: 16988
+  timestamp: 1702511261442
+- kind: conda
+  name: xz
+  version: 5.2.6
+  build: h166bdaf_0
+  subdir: linux-64
+  url: https://conda.anaconda.org/conda-forge/linux-64/xz-5.2.6-h166bdaf_0.tar.bz2
+  sha256: 03a6d28ded42af8a347345f82f3eebdd6807a08526d47899a42d62d319609162
+  md5: 2161070d867d1b1204ea749c8eec4ef0
+  depends:
+  - libgcc-ng >=12
+  license: LGPL-2.1 and GPL-2.0
+  size: 418368
+  timestamp: 1660346797927
+- kind: conda
+  name: xz
+  version: 5.2.6
+  build: h8d14728_0
+  subdir: win-64
+  url: https://conda.anaconda.org/conda-forge/win-64/xz-5.2.6-h8d14728_0.tar.bz2
+  sha256: 54d9778f75a02723784dc63aff4126ff6e6749ba21d11a6d03c1f4775f269fe0
+  md5: 515d77642eaa3639413c6b1bc3f94219
+  depends:
+  - vc >=14.1,<15
+  - vs2015_runtime >=14.16.27033
+  license: LGPL-2.1 and GPL-2.0
+  size: 217804
+  timestamp: 1660346976440
diff --git a/pixi.toml b/pixi.toml
index c3adf4db8..134454d18 100644
--- a/pixi.toml
+++ b/pixi.toml
@@ -44,6 +44,7 @@ instantiate-wflow-server = "julia --project=server --eval 'using Pkg; Pkg.instan
 
 [dependencies]
 juliaup = "*"
+python = ">=3.10"
 
 [activation]
 scripts = ["utils/env_setup.sh"]
diff --git a/src/flow.jl b/src/flow.jl
index a07f15991..ac1cc848e 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -556,14 +556,37 @@ function initialize_shallowwater_river(
     # Edge 3 => 2
     # Edge 4 => 9
     # ⋮ )
-    riverlength_bc = get(config.model, "riverlength_bc", 1.0e04)::Float64 # river length at boundary point (ghost point)
     alpha = get(config.model, "inertial_flow_alpha", 0.7)::Float64 # stability coefficient for model time step (0.2-0.7)
     h_thresh = get(config.model, "h_thresh", 1.0e-03)::Float64 # depth threshold for flow at link
     froude_limit = get(config.model, "froude_limit", true)::Bool # limit flow to subcritical according to Froude number
     floodplain_1d = floodplain
 
-    @info "Local inertial approach is used for river flow." alpha h_thresh froude_limit riverlength_bc floodplain_1d
-
+    @info "Local inertial approach is used for river flow." alpha h_thresh froude_limit floodplain_1d
+    @warn string(
+        "Providing the boundary condition `riverlength_bc` as part of the `[model]` setting ",
+        "in the TOML file has been deprecated as of Wflow v0.8.0.\n The boundary condition should ",
+        "be provided as part of the file `$(config.input.path_static)`.",
+    )
+    # The following boundary conditions can be set at ghost nodes, downstream of river
+    # outlets (pits): river length and river depth
+    index_pit = findall(x -> x == 5, ldd)
+    inds_pit = inds[index_pit]
+    riverlength_bc = ncread(
+        nc,
+        config,
+        "lateral.river.riverlength_bc";
+        sel = inds_pit,
+        defaults = 1.0e04,
+        type = Float,
+    )
+    riverdepth_bc = ncread(
+        nc,
+        config,
+        "lateral.river.riverdepth_bc";
+        sel = inds_pit,
+        defaults = 0.0,
+        type = Float,
+    )
     bankfull_elevation_2d = ncread(
         nc,
         config,
@@ -589,13 +612,15 @@ function initialize_shallowwater_river(
         ncread(nc, config, "lateral.river.n"; sel = inds, defaults = 0.036, type = Float)
 
     n = length(inds)
+
+    # set river depth h to zero (including reservoir and lake locations)
+    h = fill(0.0, n)
+
     # set ghost points for boundary condition (downstream river outlet): river width, bed
-    # elevation, manning n is copied from the upstream cell, river depth h is set at 0.0
-    # (fixed). river length at boundary point is by default 1.0e4 m (riverlength_bc).
-    index_pit = findall(x -> x == 5, ldd)
-    npits = length(index_pit)
+    # elevation, manning n is copied from the upstream cell.
     add_vertex_edge_graph!(graph, index_pit)
-    append!(dl, fill(riverlength_bc, npits))
+    append!(dl, riverlength_bc)
+    append!(h, riverdepth_bc)
     append!(zb, zb[index_pit])
     append!(width, width[index_pit])
     append!(n_river, n_river[index_pit])
@@ -639,10 +664,7 @@ function initialize_shallowwater_river(
         mannings_n_sq[i] = mannings_n * mannings_n
     end
 
-    # set depth h to zero (including reservoir and lake locations)
-    h = fill(0.0, n + length(index_pit))
     q_av = zeros(_ne)
-
     waterbody = !=(0).(reservoir_index .+ lake_index)
     active_index = findall(x -> x == 0, waterbody)
 
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index c8be4c8db..336c19d4b 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -3,10 +3,11 @@
 
 Initial part of the sbm_gwf model concept. The model contains:
     - the vertical SBM concept
-    - the following lateral components:
-        - 1-D kinematic wave for river flow
-        - 1-D kinematic wave for overland flow
-        - unconfined aquifer with groundwater flow in four directions (adjacent cells)
+    - unconfined aquifer with groundwater flow in four directions (adjacent cells)
+    - the following surface routing options:
+        - 1-D kinematic wave for river flow and 1-D kinematic wave for overland flow
+        - 1-D local inertial model for river flow (optional floodplain) and 1-D kinematic wave for overland flow
+        - 1-D local inertial model for river flow (optional floodplain) and 2-D local inertial model for overland flow
 
 The unconfined aquifer contains a recharge, river and a drain (optional) boundary.
 
@@ -30,6 +31,16 @@ function initialize_sbm_gwf_model(config::Config)
     kw_river_tstep = get(config.model, "kw_river_tstep", 0)
     kw_land_tstep = get(config.model, "kw_land_tstep", 0)
     kinwave_it = get(config.model, "kin_wave_iteration", false)::Bool
+    routing_options = ("kinematic-wave", "local-inertial")
+    floodplain_1d = get(config.model, "floodplain_1d", false)::Bool
+    river_routing = get_options(
+        config.model,
+        "river_routing",
+        routing_options,
+        "kinematic-wave",
+    )::String
+    land_routing =
+        get_options(config.model, "land_routing", routing_options, "kinematic-wave")::String
 
     nc = NCDataset(static_path)
 
@@ -99,26 +110,7 @@ function initialize_sbm_gwf_model(config::Config)
     dw = (xl .* yl) ./ dl
     sw = map(det_surfacewidth, dw, riverwidth, river)
 
-    olf = initialize_surfaceflow_land(
-        nc,
-        config,
-        inds;
-        sl = βₗ,
-        dl = dl,
-        width = sw,
-        iterate = kinwave_it,
-        tstep = kw_land_tstep,
-        Δt = Δt,
-    )
-
     graph = flowgraph(ldd, inds, pcr_dir)
-
-    # river flow (kinematic wave)
-    riverlength = riverlength_2d[inds_riv]
-    riverwidth = riverwidth_2d[inds_riv]
-    minimum(riverlength) > 0 || error("river length must be positive on river cells")
-    minimum(riverwidth) > 0 || error("river width must be positive on river cells")
-
     ldd_riv = ldd_2d[inds_riv]
     graph_riv = flowgraph(ldd_riv, inds_riv, pcr_dir)
 
@@ -126,20 +118,82 @@ function initialize_sbm_gwf_model(config::Config)
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
     frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, βₗ, n)
 
-    rf = initialize_surfaceflow_river(
-        nc,
-        config,
-        inds_riv;
-        dl = riverlength,
-        width = riverwidth,
-        reservoir_index = resindex,
-        reservoir = reservoirs,
-        lake_index = lakeindex,
-        lake = lakes,
-        iterate = kinwave_it,
-        tstep = kw_river_tstep,
-        Δt = Δt,
-    )
+    if land_routing == "kinematic-wave"
+        olf = initialize_surfaceflow_land(
+            nc,
+            config,
+            inds;
+            sl = βₗ,
+            dl,
+            width = map(det_surfacewidth, dw, riverwidth, river),
+            iterate = kinwave_it,
+            tstep = kw_land_tstep,
+            Δt,
+        )
+    elseif land_routing == "local-inertial"
+        index_river_nf = rev_inds_riv[inds] # not filtered (with zeros)
+        olf, indices = initialize_shallowwater_land(
+            nc,
+            config,
+            inds;
+            modelsize_2d,
+            indices_reverse = rev_inds,
+            xlength = xl,
+            ylength = yl,
+            riverwidth = riverwidth_2d[inds_riv],
+            graph_riv,
+            ldd_riv,
+            inds_riv,
+            river,
+            waterbody = !=(0).(resindex + lakeindex),
+            Δt,
+        )
+    end
+
+    # river flow (kinematic wave)
+    riverlength = riverlength_2d[inds_riv]
+    riverwidth = riverwidth_2d[inds_riv]
+    minimum(riverlength) > 0 || error("river length must be positive on river cells")
+    minimum(riverwidth) > 0 || error("river width must be positive on river cells")
+
+    if river_routing == "kinematic-wave"
+        rf = initialize_surfaceflow_river(
+            nc,
+            config,
+            inds_riv;
+            dl = riverlength,
+            width = riverwidth,
+            reservoir_index = resindex,
+            reservoir = reservoirs,
+            lake_index = lakeindex,
+            lake = lakes,
+            iterate = kinwave_it,
+            tstep = kw_river_tstep,
+            Δt = Δt,
+        )
+    elseif river_routing == "local-inertial"
+        rf, nodes_at_link = initialize_shallowwater_river(
+            nc,
+            config,
+            inds_riv;
+            graph = graph_riv,
+            ldd = ldd_riv,
+            dl = riverlength,
+            width = riverwidth,
+            reservoir_index = resindex,
+            reservoir = reservoirs,
+            lake_index = lakeindex,
+            lake = lakes,
+            Δt = Δt,
+            floodplain = floodplain_1d,
+        )
+    else
+        error(
+            """An unknown "river_routing" method is specified in the TOML file ($river_routing).
+            This should be "kinematic-wave" or "local-inertial".
+            """,
+        )
+    end
 
     # unconfined aquifer
     if do_constanthead
@@ -294,26 +348,33 @@ function initialize_sbm_gwf_model(config::Config)
     # setup subdomains for the land and river kinematic wave domain, if nthreads = 1
     # subdomain is equal to the complete domain
     toposort = topological_sort_by_dfs(graph)
-    toposort_riv = topological_sort_by_dfs(graph_riv)
-    index_pit_land = findall(x -> x == 5, ldd)
-    index_pit_river = findall(x -> x == 5, ldd_riv)
-    streamorder = stream_order(graph, toposort)
-    min_streamorder_land = get(config.model, "min_streamorder_land", 5)
-    subbas_order, indices_subbas, topo_subbas = kinwave_set_subdomains(
-        graph,
-        toposort,
-        index_pit_land,
-        streamorder,
-        min_streamorder_land,
-    )
-    min_streamorder_river = get(config.model, "min_streamorder_river", 6)
-    subriv_order, indices_subriv, topo_subriv = kinwave_set_subdomains(
-        graph_riv,
-        toposort_riv,
-        index_pit_river,
-        streamorder[index_river],
-        min_streamorder_river,
-    )
+    if land_routing == "kinematic-wave" || river_routing == "kinematic-wave"
+        streamorder = stream_order(graph, toposort)
+    end
+    if land_routing == "kinematic-wave"
+        toposort = topological_sort_by_dfs(graph)
+        index_pit_land = findall(x -> x == 5, ldd)
+        min_streamorder_land = get(config.model, "min_streamorder_land", 5)
+        subbas_order, indices_subbas, topo_subbas = kinwave_set_subdomains(
+            graph,
+            toposort,
+            index_pit_land,
+            streamorder,
+            min_streamorder_land,
+        )
+    end
+    if river_routing == "kinematic-wave"
+        min_streamorder_river = get(config.model, "min_streamorder_river", 6)
+        toposort_riv = topological_sort_by_dfs(graph_riv)
+        index_pit_river = findall(x -> x == 5, ldd_riv)
+        subriv_order, indices_subriv, topo_subriv = kinwave_set_subdomains(
+            graph_riv,
+            toposort_riv,
+            index_pit_river,
+            streamorder[index_river],
+            min_streamorder_river,
+        )
+    end
 
     modelmap =
         (vertical = sbm, lateral = (subsurface = subsurface_map, land = olf, river = rf))
@@ -348,29 +409,57 @@ function initialize_sbm_gwf_model(config::Config)
     # for the land domain the x and y length [m] of the grid cells are stored
     # for reservoirs and lakes indices information is available from the initialization
     # functions
-    land = (
-        graph = graph,
-        upstream_nodes = filter_upsteam_nodes(graph, pits[inds]),
-        subdomain_order = subbas_order,
-        topo_subdomain = topo_subbas,
-        indices_subdomain = indices_subbas,
-        order = toposort,
-        indices = inds,
-        reverse_indices = rev_inds,
-        xl = xl,
-        yl = yl,
-        altitude = altitude,
-    )
-    river = (
-        graph = graph_riv,
-        upstream_nodes = filter_upsteam_nodes(graph_riv, pits[inds_riv]),
-        subdomain_order = subriv_order,
-        topo_subdomain = topo_subriv,
-        indices_subdomain = indices_subriv,
-        order = toposort_riv,
-        indices = inds_riv,
-        reverse_indices = rev_inds_riv,
-    )
+    if land_routing == "kinematic-wave"
+        land = (
+            graph = graph,
+            upstream_nodes = filter_upsteam_nodes(graph, pits[inds]),
+            subdomain_order = subbas_order,
+            topo_subdomain = topo_subbas,
+            indices_subdomain = indices_subbas,
+            order = toposort,
+            indices = inds,
+            reverse_indices = rev_inds,
+            xl = xl,
+            yl = yl,
+            slope = βₗ,
+            altitude = altitude,
+        )
+    elseif land_routing == "local-inertial"
+        land = (
+            graph = graph,
+            order = toposort,
+            indices = inds,
+            reverse_indices = rev_inds,
+            xl = xl,
+            yl = yl,
+            slope = βₗ,
+            altitude = altitude,
+            index_river = index_river_nf,
+            staggered_indices = indices,
+        )
+    end
+    if river_routing == "kinematic-wave"
+        river = (
+            graph = graph_riv,
+            indices = inds_riv,
+            reverse_indices = rev_inds_riv,
+            # specific for kinematic_wave
+            upstream_nodes = filter_upsteam_nodes(graph_riv, pits[inds_riv]),
+            subdomain_order = subriv_order,
+            topo_subdomain = topo_subriv,
+            indices_subdomain = indices_subriv,
+            order = toposort_riv,
+        )
+    elseif river_routing == "local-inertial"
+        river = (
+            graph = graph_riv,
+            indices = inds_riv,
+            reverse_indices = rev_inds_riv,
+            # specific for local-inertial
+            nodes_at_link = nodes_at_link,
+            links_at_node = adjacent_links_at_node(graph_riv, nodes_at_link),
+        )
+    end
 
     model = Model(
         config,
@@ -409,7 +498,7 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
         lateral_snow_transport!(
             vertical.snow,
             vertical.snowwater,
-            lateral.land.sl,
+            network.land.slope,
             network.land,
         )
     end
@@ -461,38 +550,3 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
 
     return model
 end
-
-function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
-    @unpack lateral, config = model
-
-    reinit = get(config.model, "reinit", true)::Bool
-    do_lakes = get(config.model, "lakes", false)::Bool
-    # read and set states in model object if reinit=false
-    if reinit == false
-        instate_path = input_path(config, config.state.path_input)
-        @info "Set initial conditions from state file `$instate_path`."
-        set_states(instate_path, model, type = Float, dimname = :layer)
-        # update kinematic wave volume for river and land domain
-        @unpack lateral = model
-        # makes sure land cells with zero flow width are set to zero q and h
-        for i in eachindex(lateral.land.width)
-            if lateral.land.width[i] <= 0.0
-                lateral.land.q[i] = 0.0
-                lateral.land.h[i] = 0.0
-            end
-        end
-        lateral.land.volume .= lateral.land.h .* lateral.land.width .* lateral.land.dl
-        lateral.river.volume .= lateral.river.h .* lateral.river.width .* lateral.river.dl
-
-        if do_lakes
-            # storage must be re-initialized after loading the state with the current
-            # waterlevel otherwise the storage will be based on the initial water level
-            lakes = lateral.river.lake
-            lakes.storage .=
-                initialize_storage(lakes.storfunc, lakes.area, lakes.waterlevel, lakes.sh)
-        end
-    else
-        @info "Set initial conditions from default values."
-    end
-    return model
-end
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index cdac2694b..0994966c9 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -111,7 +111,8 @@ function initialize_sbm_model(config::Config)
 
     # check if lateral subsurface flow component is defined for the SBM model, when coupled
     # to another groundwater model, this component is not defined in the TOML file.
-    if haskey(config.input.lateral, "subsurface")
+    subsurface_flow = haskey(config.input.lateral, "subsurface")
+    if subsurface_flow
         khfrac = ncread(
             nc,
             config,
@@ -259,16 +260,23 @@ function initialize_sbm_model(config::Config)
     # setup subdomains for the land and river kinematic wave domain, if nthreads = 1
     # subdomain is equal to the complete domain
     toposort = topological_sort_by_dfs(graph)
-    index_pit_land = findall(x -> x == 5, ldd)
-    streamorder = stream_order(graph, toposort)
-    min_streamorder_land = get(config.model, "min_streamorder_land", 5)
-    subbas_order, indices_subbas, topo_subbas = kinwave_set_subdomains(
-        graph,
-        toposort,
-        index_pit_land,
-        streamorder,
-        min_streamorder_land,
-    )
+    if land_routing == "kinematic-wave" ||
+       river_routing == "kinematic-wave" ||
+       subsurface_flow
+        streamorder = stream_order(graph, toposort)
+    end
+    if land_routing == "kinematic-wave" || subsurface_flow
+        toposort = topological_sort_by_dfs(graph)
+        index_pit_land = findall(x -> x == 5, ldd)
+        min_streamorder_land = get(config.model, "min_streamorder_land", 5)
+        subbas_order, indices_subbas, topo_subbas = kinwave_set_subdomains(
+            graph,
+            toposort,
+            index_pit_land,
+            streamorder,
+            min_streamorder_land,
+        )
+    end
     if river_routing == "kinematic-wave"
         min_streamorder_river = get(config.model, "min_streamorder_river", 6)
         toposort_riv = topological_sort_by_dfs(graph_riv)
@@ -285,7 +293,7 @@ function initialize_sbm_model(config::Config)
     if nthreads() > 1
         if river_routing == "kinematic-wave"
             @info "Parallel execution of kinematic wave" min_streamorder_land min_streamorder_river
-        else
+        elseif land_routing == "kinematic-wave" || subsurface_flow
             @info "Parallel execution of kinematic wave" min_streamorder_land
         end
     end
@@ -481,7 +489,9 @@ function update_total_water_storage(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,
     return model
 end
 
-function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel}
+function set_states(
+    model::Model{N,L,V,R,W,T},
+) where {N,L,V,R,W,T<:Union{SbmModel,SbmGwfModel}}
     @unpack lateral, vertical, network, config = model
 
     reinit = get(config.model, "reinit", true)::Bool
@@ -496,12 +506,14 @@ function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmModel}
         nriv = length(network.river.indices)
         instate_path = input_path(config, config.state.path_input)
         @info "Set initial conditions from state file `$instate_path`."
-        @warn string(
-            "The unit of `ssf` (lateral subsurface flow) is now m3 d-1. Please update your",
-            " input state file if it was produced with a Wflow version up to v0.5.2.",
-        )
+        if T <: SbmModel
+            @warn string(
+                "The unit of `ssf` (lateral subsurface flow) is now m3 d-1. Please update your",
+                " input state file if it was produced with a Wflow version up to v0.5.2.",
+            )
+        end
         set_states(instate_path, model; type = Float, dimname = :layer)
-        # update zi for vertical sbm and kinematic wave volume for river and land domain
+        # update zi for vertical sbm
         zi =
             max.(
                 0.0,
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index 7e2194976..fed99f23e 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -398,10 +398,33 @@ model = Wflow.run_timestep(model)
     @test q[1622] ≈ 6.0094627478450016f-5
     @test q[43] ≈ 11.900372477232796f0
     @test q[501] ≈ 3.470259878228359f0
+    @test q[5808] ≈ 0.002199141435542889f0
     h = model.lateral.river.h_av
     @test h[1622] ≈ 0.0018099697988149294f0
     @test h[43] ≈ 0.4362704420867342f0
     @test h[501] ≈ 0.05610231297517167f0
+    @test h[5808] ≈ 0.005901095907129176f0
+end
+
+# set boundary condition local inertial routing from NetCDF file
+config.input.lateral.river.riverlength_bc = "riverlength_bc"
+config.input.lateral.river.riverdepth_bc = "riverdepth_bc"
+model = Wflow.initialize_sbm_model(config)
+model = Wflow.run_timestep(model)
+model = Wflow.run_timestep(model)
+
+@testset "change boundary condition for local inertial routing (including floodplain)" begin
+    q = model.lateral.river.q_av
+    @test sum(q) ≈ 3899.01831081148f0
+    @test q[1622] ≈ 6.008529317397127f-5
+    @test q[43] ≈ 11.900394994314782f0
+    @test q[501] ≈ 3.4702817507735366f0
+    @test q[5808] ≈ 0.060081151400640784f0
+    h = model.lateral.river.h_av
+    @test h[1622] ≈ 0.001809870969211623f0
+    @test h[43] ≈ 0.4362709200147433f0
+    @test h[501] ≈ 0.05610440991535934f0
+    @test h[5808] ≈ 2.0000006940603936f0
 end
 Wflow.close_files(model, delete_output = false)
 
diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl
index ada1a2f33..629e38389 100644
--- a/test/run_sbm_gwf.jl
+++ b/test/run_sbm_gwf.jl
@@ -38,12 +38,12 @@ model = Wflow.run_timestep(model)
     @test sbm.transpiration[1] ≈ 1.0122634204681036f0
 end
 
-@testset "overland flow" begin
+@testset "overland flow (kinematic wave)" begin
     q = model.lateral.land.q_av
     @test sum(q) ≈ 2.2298616f-7
 end
 
-@testset "river domain" begin
+@testset "river domain (kinematic wave)" begin
     q = model.lateral.river.q_av
     river = model.lateral.river
     @test sum(q) ≈ 0.034846861707851576f0
@@ -57,18 +57,133 @@ end
 @testset "groundwater" begin
     gw = model.lateral.subsurface
     @test gw.river.stage[1] ≈ 1.2123636929067039f0
-    @test gw.flow.aquifer.head[19] ≈ 1.7999999523162842f0
+    @test gw.flow.aquifer.head[17:21] ≈ [
+        1.288882128227083f0,
+        1.344816977641676f0,
+        1.7999999523162842f0,
+        1.803682542252168f0,
+        1.4049539807162532f0,
+    ]
     @test gw.river.flux[1] ≈ -50.46515313302901f0
     @test gw.drain.flux[1] ≈ 0.0
     @test gw.recharge.rate[19] ≈ -0.0014241196552847502f0
 end
 
-Wflow.close_files(model)
-
 @testset "no drains" begin
     config.model.drains = false
     delete!(Dict(config.output.lateral.subsurface), "drain")
     model = Wflow.initialize_sbm_gwf_model(config)
     @test collect(keys(model.lateral.subsurface)) == [:flow, :recharge, :river]
-    Wflow.close_files(model)
 end
+
+Wflow.close_files(model, delete_output = false)
+
+# test local-inertial option for river flow routing
+tomlpath = joinpath(@__DIR__, "sbm_gwf_config.toml")
+config = Wflow.Config(tomlpath)
+config.model.river_routing = "local-inertial"
+
+config.input.lateral.river.bankfull_elevation = "bankfull_elevation"
+config.input.lateral.river.bankfull_depth = "bankfull_depth"
+
+model = Wflow.initialize_sbm_gwf_model(config)
+model = Wflow.run_timestep(model)
+model = Wflow.run_timestep(model)
+
+@testset "river domain (local inertial)" begin
+    q = model.lateral.river.q_av
+    river = model.lateral.river
+    @test sum(q) ≈ 0.026752212882189867f0
+    @test q[6] ≈ 0.005962226036202787f0
+    @test river.volume[6] ≈ 7.498085951460259f0
+    @test river.inwater[6] ≈ 0.00017741362380231486f0
+    @test q[13] ≈ 0.0004586401133311585f0
+    @test q[5] ≈ 0.006312330297863098f0
+end
+Wflow.close_files(model, delete_output = false)
+
+# test local-inertial option for river and overland flow routing
+tomlpath = joinpath(@__DIR__, "sbm_gwf_config.toml")
+config = Wflow.Config(tomlpath)
+config.model.river_routing = "local-inertial"
+config.model.land_routing = "local-inertial"
+
+config.input.lateral.river.bankfull_elevation = "bankfull_elevation"
+config.input.lateral.river.bankfull_depth = "bankfull_depth"
+config.input.lateral.land.elevation = "wflow_dem"
+
+pop!(Dict(config.state.lateral.land), "q")
+config.state.lateral.land.h_av = "h_av_land"
+config.state.lateral.land.qx = "qx_land"
+config.state.lateral.land.qy = "qy_land"
+
+model = Wflow.initialize_sbm_gwf_model(config)
+model = Wflow.run_timestep(model)
+model = Wflow.run_timestep(model)
+
+@testset "river and land domain (local inertial)" begin
+    q = model.lateral.river.q_av
+    @test sum(q) ≈ 0.026759512926399318f0
+    @test q[6] ≈ 0.005964050519182144f0
+    @test q[13] ≈ 0.00045874863609566906f0
+    @test q[5] ≈ 0.006314336504100363f0
+    h = model.lateral.river.h_av
+    @test h[6] ≈ 0.08023321813680799f0
+    @test h[5] ≈ 0.07761589477563689f0
+    @test h[13] ≈ 0.08226917799750776f0
+    qx = model.lateral.land.qx
+    qy = model.lateral.land.qy
+    @test all(qx .== 0.0f0)
+    @test all(qy .== 0.0f0)
+end
+Wflow.close_files(model, delete_output = false)
+
+# test with warm start
+tomlpath = joinpath(@__DIR__, "sbm_gwf_config.toml")
+config = Wflow.Config(tomlpath)
+config.model.reinit = false
+
+model = Wflow.initialize_sbm_gwf_model(config)
+@unpack network = model
+
+model = Wflow.run_timestep(model)
+model = Wflow.run_timestep(model)
+
+@testset "second timestep warm start" begin
+    sbm = model.vertical
+    @test sbm.runoff[1] == 0.0
+    @test sbm.soilevap[1] == 0.2927279656884887
+    @test sbm.transpiration[1] ≈ 1.0122634204681036f0
+end
+
+@testset "overland flow warm start (kinematic wave)" begin
+    q = model.lateral.land.q_av
+    @test sum(q) ≈ 1.4411427003142072f-5
+end
+
+@testset "river domain warm start (kinematic wave)" begin
+    q = model.lateral.river.q_av
+    river = model.lateral.river
+    @test sum(q) ≈ 0.011317441400219936f0
+    @test q[6] ≈ 0.002255753266287542f0
+    @test river.volume[6] ≈ 2.1212499727096956f0
+    @test river.inwater[6] ≈ -6.888767545965421f-5
+    @test q[13] ≈ 8.664553314598283f-5
+    @test q[network.river.order[end]] ≈ 0.002258255913217909f0
+end
+
+@testset "groundwater warm start" begin
+    gw = model.lateral.subsurface
+    @test gw.river.stage[1] ≈ 1.2031171676781156f0
+    @test gw.flow.aquifer.head[17:21] ≈ [
+        1.2195537323067487f0,
+        1.2741313108104333f0,
+        1.7999999523162842f0,
+        1.5862473868459186f0,
+        1.202268433263572f0,
+    ]
+    @test gw.river.flux[1] ≈ -6.044849112655228f0
+    @test gw.drain.flux[1] ≈ 0.0
+    @test gw.recharge.rate[19] ≈ -0.0014241196552847502f0
+end
+Wflow.close_files(model, delete_output = false)
diff --git a/test/runtests.jl b/test/runtests.jl
index 23a6ab4a1..0ce2458ba 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -63,6 +63,8 @@ forcing_sbm_gw_path = testdata(
 forcing_meuse_path = testdata(v"0.2.8", "forcing_meuse.nc", "forcing_meuse.nc")
 staticmaps_sbm_gw_path =
     testdata(v"0.2.2", "staticmaps-sbm-groundwater.nc", "staticmaps-sbm-groundwater.nc")
+instates_sbm_gw_path =
+    testdata(v"0.2.2", "instates-example-sbm-gwf.nc", "instates-example-sbm-gwf.nc")
 lake_sh_1_path = testdata(v"0.2.1", "lake_sh_1.csv", "lake_sh_1.csv")
 lake_sh_2_path = testdata(v"0.2.1", "lake_sh_2.csv", "lake_sh_2.csv")
 lake_hq_2_path = testdata(v"0.2.1", "lake_hq_2.csv", "lake_hq_2.csv")