diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402.webp
deleted file mode 100644
index e9e07c5..0000000
Binary files a/content/tutorials/parallelization/SIMWE_images/basin_030401010402.webp and /dev/null differ
diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_both.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_both.webp
index b930c1d..257ef8e 100644
Binary files a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_both.webp and b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_both.webp differ
diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_cn.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_cn.webp
new file mode 100644
index 0000000..f9c156b
Binary files /dev/null and b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_cn.webp differ
diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_erdep.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_erdep.webp
index d3d9993..1a342a4 100644
Binary files a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_erdep.webp and b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_erdep.webp differ
diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_manning.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_manning.webp
new file mode 100644
index 0000000..6ca06c3
Binary files /dev/null and b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_manning.webp differ
diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_masked.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_masked.webp
deleted file mode 100644
index 2327880..0000000
Binary files a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_masked.webp and /dev/null differ
diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_runoff.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_runoff.webp
new file mode 100644
index 0000000..0653385
Binary files /dev/null and b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_runoff.webp differ
diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_simwe.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_simwe.webp
index 4092362..bb814c9 100644
Binary files a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_simwe.webp and b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_simwe.webp differ
diff --git a/content/tutorials/parallelization/SIMWE_images/basin_030401010402_soils.webp b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_soils.webp
new file mode 100644
index 0000000..a273ac5
Binary files /dev/null and b/content/tutorials/parallelization/SIMWE_images/basin_030401010402_soils.webp differ
diff --git a/content/tutorials/parallelization/SIMWE_images/webmap.webp b/content/tutorials/parallelization/SIMWE_images/webmap.webp
index d06455a..515397b 100644
Binary files a/content/tutorials/parallelization/SIMWE_images/webmap.webp and b/content/tutorials/parallelization/SIMWE_images/webmap.webp differ
diff --git a/content/tutorials/parallelization/SIMWE_parallelization.qmd b/content/tutorials/parallelization/SIMWE_parallelization.qmd
index 6747649..0217259 100644
--- a/content/tutorials/parallelization/SIMWE_parallelization.qmd
+++ b/content/tutorials/parallelization/SIMWE_parallelization.qmd
@@ -3,12 +3,13 @@ title: "Parallelization of overland flow and sediment transport simulation"
 author: "Anna Petrasova"
 date: 2025-09-25
 date-modified: today
-categories: [Python, advanced, parallelization, hydrology, erosion]
+categories: [Python, advanced, soils, parallelization, hydrology, erosion]
 description: >
   In this tutorial we will run overland flow
   and erosion simulation (SIMWE) split by subwatersheds
   in parallel.
 image: SIMWE_images/thumbnail.webp
+bibliography: references.bib
 format:
   ipynb: default
   html:
@@ -41,9 +42,9 @@ subwatersheds, each approximately 100 km² in size) and run simulations
 independently for each unit.
 
 This approach demonstrates how to leverage new Python API features introduced
-in GRASS 8.5—specifically, the [Tools API](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.tools.html#)
-and the [MaskManager](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.script.html#grass.script.MaskManager)
-and [RegionManager](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.script.html#grass.script.RegionManager)
+in GRASS 8.5—specifically, the [Tools API](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.tools.html#)
+and the [MaskManager](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#grass.script.MaskManager)
+and [RegionManager](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#grass.script.RegionManager)
 context managers—which simplify region and mask handling
 and are especially useful in parallel workflows.
 
@@ -136,7 +137,7 @@ gs.create_project(nlcd_project, filename=nlcd_filename)
 session = gj.init(nlcd_project)
 ```
 
-Link the NLCD raster with [r.external](https://grass.osgeo.org/grass-devel/manuals/r.external.html). This command creates a virtual raster without importing the full dataset, allowing us to work with just the small sections needed from the larger, nationwide NLCD file.
+Link the NLCD raster with [r.external](https://grass.osgeo.org/grass-stable/manuals/r.external.html). This command creates a virtual raster without importing the full dataset, allowing us to work with just the small sections needed from the larger, nationwide NLCD file.
 
 ```{python}
 tools = Tools()
@@ -219,7 +220,7 @@ tools.v_in_ogr(
 ```
 
 Next, we want to extract the subwatersheds along the river.
-If we simply overlap (with [v.select](https://grass.osgeo.org/grass-devel/manuals/v.select.html))
+If we simply overlap (with [v.select](https://grass.osgeo.org/grass-stable/manuals/v.select.html))
 the river and the subwatersheds, we will miss some of them
 because the river data don't always overlap or touch the subwatersheds.
 
@@ -239,7 +240,7 @@ basin_map.show()
 
 ![](SIMWE_images/basins.webp)
 
-So instead we will [buffer](https://grass.osgeo.org/grass-devel/manuals/v.buffer.html) the river first:
+So instead we will [buffer](https://grass.osgeo.org/grass-stable/manuals/v.buffer.html) the river first:
 
 ```{python}
 tools.v_buffer(input="river", output="river_buffer", distance=0.01)
@@ -265,12 +266,23 @@ The rest of the workflow will be done in a CRS used in North Carolina (EPSG 3358
 Since we want our project to use a different CRS more suitable for our study area (EPSG:3358 for North Carolina), we will create it now:
 
 ```{python}
-nc_project = path / "NC"
+nc_project = Path(path) / "NC"
 gs.create_project(nc_project, epsg="3358")
 session = gj.init(nc_project)
 ```
 
-We will reproject the river subwatersheds vector with [v.proj]((https://grass.osgeo.org/grass-devel/manuals/v.proj.html)):
+In preparation for parallel importing of soil data,
+we will specify the database connection to use per-map sqlite database
+instead of a single database for the entire mapset (the default).
+This will prevent concurrent writing to a single database.
+
+```{python}
+tools.db_connect(
+    driver="sqlite", database="$GISDBASE/$LOCATION_NAME/$MAPSET/vector/$MAP/sqlite.db"
+)
+```
+
+We will reproject the river subwatersheds vector with [v.proj]((https://grass.osgeo.org/grass-stable/manuals/v.proj.html)):
 
 ```{python}
 tools.v_proj(input="river_basins", project="hydro")
@@ -278,14 +290,15 @@ tools.v_proj(input="river_basins", project="hydro")
 
 ## Downloading NED elevation data
 
-We will use GRASS addon [r.in.usgs](https://grass.osgeo.org/grass-devel/manuals/addons/r.in.usgs.html) to download and import NED layer using [USGS TNM Access API](https://apps.nationalmap.gov/tnmaccess/).
+We will use GRASS addon [r.in.usgs](https://grass.osgeo.org/grass-stable/manuals/addons/r.in.usgs.html) to download and import NED layer using [USGS TNM Access API](https://apps.nationalmap.gov/tnmaccess/).
 
 ```{python}
 tools.g_extension(extension="r.in.usgs")
 ```
 
-First, we will set the computational region to match the river subwatersheds and then list available datasets
-for that region (the output at the time of creating this tutorial, it may vary later on):
+First, we will set the computational region to match the river subwatersheds
+and then list the latest available datasets for that region
+(the output at the time of creating this tutorial, it may vary later on):
 
 ```{python}
 tools.g_region(vector="river_basins")
@@ -293,6 +306,7 @@ print(
     tools.r_in_usgs(
         product="ned",
         ned_dataset="ned13sec",
+        ned_release="current",
         flags="i",
         output_directory=os.getcwd(),
     ).text
@@ -302,39 +316,27 @@ print(
 ```text
 USGS file(s) to download:
 -------------------------
-Total download size:	5.32 GB
-Tile count:	17
+Total download size:	2.90 GB
+Tile count:	6
 USGS SRS:	wgs84
 USGS tile titles:
-USGS 1/3 Arc Second n36w080 20240611
-USGS 1/3 Arc Second n36w081 20240611
-USGS 1/3 Arc Second n36w082 20220504
-USGS 1/3 Arc Second n36w082 20220512
-USGS 1/3 Arc Second n36w082 20240611
-USGS 1/3 Arc Second n37w080 20210305
-USGS 1/3 Arc Second n37w081 20210305
-USGS 1/3 Arc Second n37w081 20240611
-USGS 1/3 Arc Second n37w082 20210305
-USGS 1/3 Arc Second n37w082 20220512
-USGS 1/3 Arc Second n37w082 20240611
-USGS 13 arc-second n36w080 1 x 1 degree
-USGS 13 arc-second n36w081 1 x 1 degree
-USGS 13 arc-second n36w082 1 x 1 degree
-USGS 13 arc-second n37w080 1 x 1 degree
-USGS 13 arc-second n37w081 1 x 1 degree
-USGS 13 arc-second n37w082 1 x 1 degree
+USGS 1/3 Arc Second n36w080 20100929
+USGS 1/3 Arc Second n36w081 20100929
+USGS 1/3 Arc Second n36w082 20100929
+USGS 1/3 Arc Second n37w080 20100929
+USGS 1/3 Arc Second n37w081 20100929
+USGS 1/3 Arc Second n37w082 20100929
 -------------------------
 ```
 
-We will select the 2024 data and use the `title_filter` option to filter them.
-We will download and import the datasets.
+Now we will download and import the datasets.
 
 ```{python}
 tools.r_in_usgs(
     product="ned",
     ned_dataset="ned13sec",
+    ned_release="current",
     output_directory=os.getcwd(),
-    title_filter="20240611",
     output_name="ned",
 )
 ```
@@ -345,8 +347,8 @@ Let's now compute the erosion/deposition for one of the subwatersheds.
 
 We begin our analysis by selecting a single HUC12 subwatershed from the set of intersected basins. This subset will be used to test the erosion and deposition simulation workflow. For all outputs, we will use unique names with the subwatershed code as a prefix. That will help us with parallelization later on.
 
-Use [v.db.select](https://grass.osgeo.org/grass-devel/manuals/v.db.select.html) to retrieve the list of HUC12 codes from the attribute table that has "areasqkm" and "huc12" attributes.
-Extract the smallest HUC12 polygon using [v.extract](https://grass.osgeo.org/grass-devel/manuals/v.extract.html).
+Use [v.db.select](https://grass.osgeo.org/grass-stable/manuals/v.db.select.html) to retrieve the list of HUC12 codes from the attribute table that has "areasqkm" and "huc12" attributes.
+Extract the smallest HUC12 polygon using [v.extract](https://grass.osgeo.org/grass-stable/manuals/v.extract.html).
 
 ```{python}
 # Select for example the smallest huc12
@@ -360,152 +362,237 @@ tools.v_extract(
 )
 ```
 
-Set the computational region with [g.region](https://grass.osgeo.org/grass-devel/manuals/g.region.html) to match the bounds of the selected subwatershed.
+Set the computational region with [g.region](https://grass.osgeo.org/grass-stable/manuals/g.region.html) to match the extent of the basin and align resolution with the elevation raster.
 
 ```{python}
-tools.g_region(vector=f"basin_{huc12}")
+tools.g_region(raster="ned", vector=f"basin_{huc12}")
 ```
 
-Display the NED data for that region:
+Convert the selected vector to a raster using [v.to.rast](https://grass.osgeo.org/grass-stable/manuals/v.to.rast.html).
 
 ```{python}
-basin_map = gj.Map(current_region=True)
-basin_map.d_rast(map="ned")
-basin_map.d_vect(map=f"basin_{huc12}", fill_color="none")
-basin_map.show()
+tools.v_to_rast(input=f"basin_{huc12}", output=f"basin_{huc12}", use="val")
+```
+
+Apply a raster mask with [MaskManager](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#grass.script.MaskManager) to restrict all subsequent raster operations to this subwatershed.
+MaskManager is new in GRASS version 8.5. Alternatively, you can use [r.mask](https://grass.osgeo.org/grass-stable/manuals/r.mask.html) to set a mask.
+
+```{python}
+mask = gs.MaskManager(mask_name=f"basin_{huc12}")
+mask.activate()
 ```
 
-![](SIMWE_images/basin_030401010402.webp)
 
-Set the region to match the extent of the basin and align resolution with the elevation raster.
-Convert the selected vector to a raster using [v.to.rast](https://grass.osgeo.org/grass-devel/manuals/v.to.rast.html).
+## Mannings roughness
+
+Next, we will estimate Manning's roughness coefficient from the NLCD landcover
+using [r.manning](https://grass.osgeo.org/grass-stable/manuals/addons/r.manning.html) addon tool.
+
+First, reproject NLCD from a "nlcd" project to this project (see [NLCD legend](https://www.usgs.gov/media/images/annual-nlcd-land-cover-change-legend)):
 
 ```{python}
-tools.g_region(raster="ned", vector=f"basin_{huc12}")
-tools.v_to_rast(input=f"basin_{huc12}", output=f"basin_{huc12}", use="val")
+tools.r_proj(project="nlcd", mapset="PERMANENT", input="nlcd", output=f"nlcd_{huc12}")
+
+nlcd_map = gj.Map()
+nlcd_map.d_rast(map=f"nlcd_{huc12}")
+nlcd_map.show()
 ```
 
-Apply a raster mask with [MaskManager](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.script.html#grass.script.MaskManager) to restrict all subsequent raster operations to this subwatershed.
-MaskManager is new in GRASS version 8.5. Alternatively, you can use [r.mask](https://grass.osgeo.org/grass-devel/manuals/r.mask.html) to set a mask.
+![](SIMWE_images/basin_030401010402_nlcd.webp)
+
+To derive manning's roughness,
+we will use lookup table [@kalyanapu2009] that is more suitable for shallow flow.
 
 ```{python}
-mask = gs.MaskManager(mask_name=f"basin_{huc12}")
-mask.activate()
+tools.g_extension(extension="r.manning")
+tools.r_manning(
+    input=f"nlcd_{huc12}",
+    output=f"mannings_{huc12}",
+    landcover="nlcd",
+    source="kalyanapu",
+)
 
-basin_map = gj.Map(current_region=True)
-basin_map.d_rast(map="ned")
-basin_map.d_legend(
-    raster="ned",
+manning_map = gj.Map()
+manning_map.d_rast(map=f"mannings_{huc12}")
+manning_map.d_legend(
+    raster=f"mannings_{huc12}",
     flags="t",
     at=[10, 15, 40, 95],
-    title="Elevation [m]",
+    title="Roughness",
     fontsize=12,
 )
-basin_map.show()
+manning_map.show()
 ```
 
-![](SIMWE_images/basin_030401010402_masked.webp)
+![](SIMWE_images/basin_030401010402_manning.webp)
+
+## Rainfall excess
 
-Reproject NLCD from a "nlcd" project to this project (see [NLCD legend](https://www.usgs.gov/media/images/annual-nlcd-land-cover-change-legend)):
+Then we will estimate rainfall excess and to do that we will first process
+soil data for our area using [r.in.ssurgo](https://grass.osgeo.org/grass-stable/manuals/addons/r.in.ssurgo.html) addon tool.
+Since we need a fairly large area, we will first download the soil data itself for North Carolina
+from [Gridded Soil Survey Geographic (gSSURGO) Database](https://www.nrcs.usda.gov/resources/data-and-reports/gridded-soil-survey-geographic-gssurgo-database)
+and then have r.in.ssurgo extract Hydrologic Soil Group for our area.
+
+::: {.callout-note title="What is Hydrologic Soil Group?"}
+Hydrologic Soil Group (HSG) classifies soils based on their minimum
+infiltration rate and potential for rainfall runoff.
+For smaller areas, r.in.ssurgo can download the data directly using the
+[Soil Data Access (SDA) API](https://sdmdataaccess.nrcs.usda.gov/).
+:::
 
 ```{python}
-tools.r_proj(project="nlcd", mapset="PERMANENT", input="nlcd", output=f"nlcd_{huc12}")
+tools.g_extension(extension="r.in.ssurgo")
+tools.r_in_ssurgo(
+    ssurgo_path="gSSURGO_NC.zip",
+    soils=f"soils_{huc12}",
+    hydgrp=f"hydgrp_{huc12}",
+    nprocs=1,
+)
 
-basin_map = gj.Map()
-basin_map.d_rast(map=f"nlcd_{huc12}")
-basin_map.show()
+soil_map = gj.Map()
+soil_map.d_rast(map=f"hydgrp_{huc12}")
+soil_map.d_legend(
+    raster=f"hydgrp_{huc12}",
+    at=[0, 30, 0, 10],
+    title="HSG",
+    fontsize=12,
+    flags="cn",
+)
+soil_map.show()
 ```
 
-![](SIMWE_images/basin_030401010402_nlcd.webp)
+![](SIMWE_images/basin_030401010402_soils.webp)
 
-We will create a file `mannings.txt` to reclassify NLCD to manning's coefficients. The values are suggestions from the [r.sim.water](https://grass.osgeo.org/grass-devel/manuals/r.sim.water.html) documentation. The format is described in [r.recode](https://grass.osgeo.org/grass-devel/manuals/r.recode.html) documentation.
+With the soil and landcover data available, we will generate Curve Number,
+which is an empirical parameter to estimate surface runoff.
+This will allow us to translate rainfall and land characteristics
+into the runoff input for SIMWE.
 
-```text
-%%writefile mannings.txt
-11:11:0.001
-21:21:0.0404
-22:23:0.0678
-24:24:0.0404
-31:31:0.0113
-41:41:0.36
-42:42:0.32
-43:52:0.4
-52:52:0.4
-71:71:0.368
-81:82:0.325
-90:90:0.086
-95:95:0.1825
-```
+::: {.callout-note title="What is a Curve Number?"}
+The [SCS Curve Number (CN)](https://www.hec.usace.army.mil/confluence/rasdocs/ras1dtechref/7.0/overview-of-optional-capabilities/modeling-precipitation-and-infiltration/curve-number)
+method is an empirical approach
+developed by the USDA Soil Conservation Service (now NRCS),
+for estimating direct runoff (effective rainfall) from a rainfall event.
+
+The curve number is a dimensionless parameter (0–100) that summarizes
+a surface's runoff potential based on soil type (hydrologic soil group),
+land cover, and antecedent moisture conditions.
+Low values (~30–50) mean high infiltration/low runoff (e.g. sandy soils, forest);
+high values (~80–98) mean low infiltration/high runoff (e.g. clay, pavement, urban surfaces).
+:::
 
 ```{python}
-tools.r_recode(input=f"nlcd_{huc12}", output=f"mannings_{huc12}", rules="mannings.txt")
+tools.g_extension(extension="r.curvenumber")
+tools.r_curvenumber(
+    landcover=f"nlcd_{huc12}",
+    landcover_source="nlcd",
+    soil=f"hydgrp_{huc12}",
+    output=f"curvenumber_{huc12}",
+)
+
+cn_map = gj.Map()
+cn_map.d_rast(map=f"curvenumber_{huc12}")
+cn_map.d_legend(
+    raster=f"curvenumber_{huc12}",
+    at=[10, 15, 40, 95],
+    flags="t",
+    title="CN",
+    fontsize=12,
+)
+cn_map.show()
 ```
 
-Similarly, we will reclassify rainfall excess (rainfall minus infiltration).
-These rainfall excess values were derived using the SCS Curve Number method for 50 mm/h of precipitation, assigning typical Curve Numbers based on common hydrologic soil groups for each NLCD land cover class;
-they provide a rough estimate and should be refined with local soil and land management data when available.
+![](SIMWE_images/basin_030401010402_cn.webp)
+
+To get an approximate rainfall depth, we will model a 10-year storm event
+that lasts 1 hour, using data from
+[NOAA Atlas 14](https://hdsc.nws.noaa.gov/pfds/) Precipitation Frequency Data Server.
+
+```{python}
+tools.g_extension(extension="r.noaa.atlas14")
+duration = 60
+depths = tools.r_noaa_atlas14(
+    mode="point", statistic="depth", units="metric", format="json"
+)
+depth = [
+    row["values"]
+    for row in depths["table"]["rows"]
+    if row["duration"] == f"{duration}-min"
+][0]["10"]
+print(f"Depth: {depth} mm")
+```
 
 ```text
-%%writefile runoff.txt
-11:11:0.0
-21:21:3.3     # Developed Open Space (CN ~68, HSG B)
-22:22:7.9     # Developed Low Intensity (CN ~81, HSG C)
-23:23:19.8    # Developed Medium Intensity (CN ~90, HSG C)
-24:24:27.9    # Developed High Intensity (CN ~94, HSG D)
-31:31:6.1     # Barren (CN ~75, HSG B)
-41:41:0.0     # Forest (CN ~55, HSG A) — P < 0.2S → Q = 0
-42:42:0.0
-43:43:0.0
-52:52:3.3     # Shrub/Scrub (CN ~68, HSG B)
-71:71:2.0     # Grassland (CN ~61, HSG A)
-81:81:5.1     # Pasture (CN ~74, HSG B)
-82:82:14.7    # Cropland (CN ~82, HSG C)
-90:90:0.0     # Woody Wetlands (CN ~70, HSG D) — often saturated, still low Q here
-95:95:0.0     # Herbaceous Wetlands
+Depth: 54.0 mm
 ```
 
+We will then create an input rainfall depth raster, assuming uniform rainfall:
+
 ```{python}
-tools.r_recode(input=f"nlcd_{huc12}", output=f"runoff_{huc12}", rules="runoff.txt")
+tools.r_mapcalc(expression=f"rainfall_{huc12} = {depth}")
 ```
 
-Then, compute slope components (dx and dy) from the elevation raster to support hydrologic modeling.
+With that we estimate runoff (rainfall excess) with [r.runoff](https://grass.osgeo.org/grass-stable/manuals/r.runoff.html),
+which estimates for each cell total runoff depth in mm for the duration of the storm
+using SCS Curve Number method.
 
 ```{python}
-tools.r_slope_aspect(elevation="ned", dx=f"dx_{huc12}", dy=f"dy_{huc12}")
+tools.g_extension(extension="r.runoff")
+duration_h = duration / 60
+tools.r_runoff(
+    rainfall=f"rainfall_{huc12}",
+    duration=duration_h,
+    curve_number=f"curvenumber_{huc12}",
+    runoff_depth=f"runoff_depth_{huc12}",
+)
+tools.r_mapcalc(expression=f"runoff_{huc12} = runoff_depth_{huc12} / {duration_h}")
+
+runoff_map = gj.Map()
+runoff_map.d_rast(map=f"runoff_{huc12}")
+runoff_map.d_legend(
+    raster=f"runoff_{huc12}",
+    at=[10, 15, 40, 95],
+    flags="t",
+    title="Runoff [mm/h]",
+    fontsize=12,
+)
+runoff_map.show()
 ```
 
+![](SIMWE_images/basin_030401010402_runoff.webp)
+
+## Hydrologic and sediment modeling
+
 Run hydrologic and sediment simulations for the selected subwatershed.
-First, simulate overland flow using [r.sim.water](https://grass.osgeo.org/grass-devel/manuals/r.sim.water.html)
+First, simulate overland flow using [r.sim.water](https://grass.osgeo.org/grass-stable/manuals/r.sim.water.html)
 with inputs for topography, Manning’s coefficients, and rainfall intensity.
 Running the simulation may take a while.
 
 ```{python}
-simulation_time = 180
 tools.r_sim_water(
     elevation="ned",
-    dx=f"dx_{huc12}",
-    dy=f"dy_{huc12}",
     depth=f"depth_{huc12}",
-    niterations=simulation_time,
+    niterations=duration,
     man=f"mannings_{huc12}",
     rain=f"runoff_{huc12}",
 )
 
-basin_map = gj.Map()
-basin_map.d_rast(map=f"depth_{huc12}")
-basin_map.d_legend(
+simwe_map = gj.Map()
+simwe_map.d_rast(map=f"depth_{huc12}")
+simwe_map.d_legend(
     raster=f"depth_{huc12}",
     flags="t",
     at=[10, 15, 40, 95],
     title="Depth [m]",
     fontsize=12,
 )
-basin_map.show()
+simwe_map.show()
 ```
 
 ![](SIMWE_images/basin_030401010402_simwe.webp)
 
-Then, define parameters for sediment transport and run [r.sim.sediment](https://grass.osgeo.org/grass-devel/manuals/r.sim.sediment.html)
+Then, define parameters for sediment transport and run [r.sim.sediment](https://grass.osgeo.org/grass-stable/manuals/r.sim.sediment.html)
 to compute erosion and deposition patterns. We will trim the resulting raster's edges to avoid extreme values at the edge.
 The additional parameters for sediment erosion modeling are based on [WEPP](https://www.ars.usda.gov/midwest-area/west-lafayette-in/national-soil-erosion-research/docs/wepp/research/) model, here we use just a single-value estimate.
 
@@ -516,20 +603,18 @@ tools.r_mapcalc(expression=f"shear_stress_{huc12} = 0.01")
 tools.r_sim_sediment(
     elevation="ned",
     water_depth=f"depth_{huc12}",
-    dx=f"dx_{huc12}",
-    dy=f"dy_{huc12}",
     detachment_coeff=f"detachment_coef_{huc12}",
     transport_coeff=f"transport_capacity_{huc12}",
     shear_stress=f"shear_stress_{huc12}",
-    niterations=simulation_time,
+    niterations=duration,
     erosion_deposition=f"erdep_tmp_{huc12}",
-    nwalkers=5 * gs.region()["cells"],
+    nwalkers=10 * gs.region()["cells"],
 )
 tools.r_grow(input=f"erdep_tmp_{huc12}", output=f"erdep_{huc12}", radius=-2.01)
 
-basin_map = gj.Map()
-basin_map.d_rast(map=f"erdep_{huc12}")
-basin_map.d_legend(
+erdep_map = gj.Map()
+erdep_map.d_rast(map=f"erdep_{huc12}")
+erdep_map.d_legend(
     raster=f"erdep_{huc12}",
     flags="t",
     at=[10, 15, 40, 95],
@@ -538,20 +623,20 @@ basin_map.d_legend(
     title="Erosion/dep. [kg/m2s]",
     fontsize=12,
 )
-basin_map.show()
+erdep_map.show()
 ```
 
 ![](SIMWE_images/basin_030401010402_erdep.webp)
 
 Calculate total erosion for the subwatershed.
-Negative values from [r.sim.sediment](https://grass.osgeo.org/grass-devel/manuals/r.sim.sediment.html)
+Negative values from [r.sim.sediment](https://grass.osgeo.org/grass-stable/manuals/r.sim.sediment.html)
 output are treated as erosion and converted to sediment mass.
-Summary statistics are computed using [r.univar](https://grass.osgeo.org/grass-devel/manuals/r.univar.html).
+Summary statistics are computed using [r.univar](https://grass.osgeo.org/grass-stable/manuals/r.univar.html).
 
 ```{python}
 # Erosion: extract negative values and convert to positive mass [kg/m^2]
 tools.r_mapcalc(
-    f"erosion_{huc12} = if(erdep_{huc12} < 0, abs(erdep_{huc12}) * {simulation_time * 60}, null())"
+    expression=f"erosion_{huc12} = if(erdep_{huc12} < 0, abs(erdep_{huc12}) * {duration}, null())"
 )
 erosion = tools.r_univar(map=f"erosion_{huc12}", format="json")
 print(erosion["mean"])
@@ -585,7 +670,7 @@ Each subwatershed needs to set different computational region and mask.
 However those setting are usually global for each mapset.
 So, to use different regions and masks for each parallel process, we will use the region and mask context managers.
 
-* Computational region is handled using `RegionManager`, a [context manager for setting and managing computational region](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.script.html#grass.script.RegionManager),
+* Computational region is handled using `RegionManager`, a [context manager for setting and managing computational region](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#grass.script.RegionManager),
 making it possible to have custom region for the current process. This feature is available only since GRASS 8.5.
 
     ```python
@@ -594,7 +679,7 @@ making it possible to have custom region for the current process. This feature i
         tools.r_sim_water(...)
     ```
 
-* Masking is handled using `MaskManager`, a [context manager for setting and managing raster mask](https://grass.osgeo.org/grass-devel/manuals/libpython/grass.script.html#grass.script.MaskManager),
+* Masking is handled using `MaskManager`, a [context manager for setting and managing raster mask](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#grass.script.MaskManager),
 making it possible to have custom mask for the current process. This feature is available only since GRASS 8.5.
 
     ```python
@@ -612,15 +697,17 @@ The resulting statistics for each subwatershed are collected and stored in a JSO
 %%writefile script.py
 import os
 import json
+from pathlib import Path
 from multiprocessing import Pool, cpu_count
+from functools import partial
 from tqdm import tqdm
 
 import grass.script as gs
 from grass.tools import Tools
 
 
-def compute(huc12):
-    simulation_time = 180
+def compute(huc12, precipitation_depth, duration, ssurgo_path):
+    duration_h = duration / 60
     # set overwrite to True to rerun the workflow
     tools = Tools(overwrite=True, quiet=True)
     tools.v_extract(
@@ -643,39 +730,40 @@ def compute(huc12):
             input="nlcd",
             output=f"nlcd_{huc12}",
         )
-        tools.r_recode(
-            input=f"nlcd_{huc12}",
-            output=f"mannings_{huc12}",
-            rules="mannings.txt",
-        )
-        tools.r_recode(
-            input=f"nlcd_{huc12}",
-            output=f"runoff_{huc12}",
-            rules="runoff.txt",
-        )
         with gs.MaskManager(mask_name=f"basin_{huc12}"):
             # Run actual computation with active mask.
-            tools.r_recode(
+            tools.r_manning(
                 input=f"nlcd_{huc12}",
                 output=f"mannings_{huc12}",
-                rules="mannings.txt",
+                landcover="nlcd",
+                source="kalyanapu",
             )
-            tools.r_recode(
-                input=f"nlcd_{huc12}",
-                output=f"runoff_{huc12}",
-                rules="runoff.txt",
+            tools.r_in_ssurgo(
+                ssurgo_path=ssurgo_path,
+                soils=f"soils_{huc12}",
+                hydgrp=f"hydgrp_{huc12}",
+                nprocs=1,
             )
-            tools.r_slope_aspect(
-                elevation="ned",
-                dx=f"dx_{huc12}",
-                dy=f"dy_{huc12}",
+            tools.r_curvenumber(
+                landcover=f"nlcd_{huc12}",
+                landcover_source="nlcd",
+                soil=f"hydgrp_{huc12}",
+                output=f"curvenumber_{huc12}",
+            )
+            tools.r_mapcalc(expression=f"rainfall_{huc12} = {precipitation_depth}")
+            tools.r_runoff(
+                rainfall=f"rainfall_{huc12}",
+                duration=duration_h,
+                curve_number=f"curvenumber_{huc12}",
+                runoff_depth=f"runoff_depth_{huc12}",
+            )
+            tools.r_mapcalc(
+                expression=f"runoff_{huc12} = runoff_depth_{huc12} / {duration_h}"
             )
             tools.r_sim_water(
                 elevation="ned",
-                dx=f"dx_{huc12}",
-                dy=f"dy_{huc12}",
                 depth=f"depth_{huc12}",
-                niterations=simulation_time,
+                niterations=duration,
                 man=f"mannings_{huc12}",
                 rain=f"runoff_{huc12}",
             )
@@ -686,14 +774,12 @@ def compute(huc12):
             tools.r_sim_sediment(
                 elevation="ned",
                 water_depth=f"depth_{huc12}",
-                dx=f"dx_{huc12}",
-                dy=f"dy_{huc12}",
                 detachment_coeff=f"detachment_coef_{huc12}",
                 transport_coeff=f"transport_capacity_{huc12}",
                 shear_stress=f"shear_stress_{huc12}",
-                niterations=simulation_time,
+                niterations=duration,
                 erosion_deposition=f"erdep_tmp_{huc12}",
-                nwalkers=5 * region["cells"],
+                nwalkers=10 * region["cells"],
             )
             tools.r_grow(
                 input=f"erdep_tmp_{huc12}",
@@ -702,7 +788,7 @@ def compute(huc12):
             )
             # Erosion: extract negative values and convert to positive mass [kg/m^2]
             tools.r_mapcalc(
-                expression=f"erosion_{huc12} = if(erdep_{huc12} < 0, abs(erdep_{huc12}) * {simulation_time * 60}, 0)",
+                expression=f"erosion_{huc12} = if(erdep_{huc12} < 0, abs(erdep_{huc12}) * {duration}, 0)",
             )
             erosion = tools.r_univar(map=f"erosion_{huc12}", format="json", flags="e")
             return {
@@ -720,11 +806,29 @@ if __name__ == "__main__":
     tools.g_gisenv(set="NPROCS=1")
     basins = tools.v_db_select(format="json", map="river_basins")["records"]
     huc12s = [basin["huc12"] for basin in basins]
+    # single depth value for entire area
+    duration = 60
+    precipitation = tools.r_noaa_atlas14(
+        mode="point", statistic="depth", units="metric", format="json"
+    )
+    precipitation_depth = [
+        row["values"]
+        for row in precipitation["table"]["rows"]
+        if row["duration"] == f"{duration}-min"
+    ][0]["10"]
+    ssurgo_path = str(Path(__file__).parent / "gSSURGO_NC.zip")
+    worker = partial(
+        compute,
+        precipitation_depth=precipitation_depth,
+        duration=duration,
+        ssurgo_path=ssurgo_path,
+    )
     # set the number of processes to be used for the computation in total
     with Pool(processes=cpu_count()) as pool:
-        result = list(tqdm(pool.imap(compute, huc12s), total=len(huc12s)))
+        result = list(tqdm(pool.imap(worker, huc12s), total=len(huc12s)))
         with open("result.json", "w") as fp:
             json.dump(result, fp)
+
 ```
 
 Now execute the script, it will take some time.
diff --git a/content/tutorials/parallelization/references.bib b/content/tutorials/parallelization/references.bib
new file mode 100644
index 0000000..ff0eb20
--- /dev/null
+++ b/content/tutorials/parallelization/references.bib
@@ -0,0 +1,8 @@
+@article{kalyanapu2009,
+  title={Effect of land use-based surface roughness on hydrologic model output.},
+  author={Kalyanapu, Alfred J and Burian, Steven J and McPherson, Timothy N},
+  journal={Journal of Spatial Hydrology},
+  volume={9},
+  number={2},
+  year={2009}
+}
\ No newline at end of file