Author

Corey T. White

Modified

April 2, 2026

[WIP] Clay Center, KS - Soil informed rainfall excess overland flow simulation

This notebook demonstrates how to import SSURGO soil data into GRASS using the r.in.ssurgo extension. The r.in.ssurgo tool can import both local SSURGO data and data from the Soil Data Access (SDA) web service. The imported data includes soil areas, hydrologic groups, ksat maps, and mukey maps, which can be used for hydrological modeling and analysis in GRASS.

Objectives

  1. Import SSURGO soil data into GRASS using r.in.ssurgo.
  2. Visualize the imported soil maps, including soil areas, hydrologic groups, ksat maps, and mukey maps.
  3. Calculate curve number using the imported hydrologic group map and a land cover map.
  4. Calculate runoff depth and peak discharge using the curve number map, flow direction, and time of concentration.
  5. Run a basic SIMWE simulation using the imported soil data and visualize the results.

We will examine the Clay Center, KS (39.4003765, -97.0630433) site for this demonstration, which has a mapset called ssurgo_demo that we will use to import the SSURGO data and run the SIMWE model.

Set the computational region to the elevation raster and create a relief raster for hillshading.

Code 
region_text = tools.g_region(raster="elevation", flags="pac").text
print(region_text)
projection:      1 (UTM)
zone:            14
datum:           wgs84
ellipsoid:       wgs84
north:           4363578.81553571
south:           4362417.81553571
west:            666092.80807312
east:            667469.80807312
nsres:           3
ewres:           3
rows:            387
cols:            459
cells:           177633
center easting:  666781.308073
center northing: 4362998.315536

Let’s also create a relief raster for hillshading or soil maps.

Code 
tools.r_relief(input="elevation", output="relief")

Import SSURGO data using r.in.ssurgo. Note that this will import the soil areas map, hydrologic group map, and ksat maps for low, regular, and high values. The mukey column is also imported as a raster with a categorical color scheme applied. The hydrologic group map also has a categorical color scheme applied.

SDA Import

Code 
tools.r_in_ssurgo(
    soils="soil_areas_sda",
    hydgrp="hydgrp_sda",
    ksat_l="ksat_l_sda",
    ksat_r="ksat_r_sda",
    ksat_h="ksat_h_sda",
    mukey="mukey_sda",
    hzdept_r=0,
    hzdepb_r=100,
    desgnmaster="A",
    nprocs=30
)

Visualize Imported Data

The tool imports the soil polygons as a vector map (soil_areas_sda), which contains the soil polygons with attributes from the SSURGO database.

::: {#cell-SDA soil attributes .cell execution_count=6}

cat mukey mukey_int cokey comppct_r hydgrp ksat_l ksat_r ksat_h hsg
0 1 1454645 1454645 26334409 85 B 15.228 32.399997 50.795997 2
1 2 1454645 1454645 26334409 85 B 15.228 32.399997 50.795997 2
2 3 1454645 1454645 26334409 85 B 15.228 32.399997 50.795997 2
3 4 1454645 1454645 26334409 85 B 15.228 32.399997 50.795997 2
4 5 1454645 1454645 26334409 85 B 15.228 32.399997 50.795997 2

:::

MUKEY Map

Ksat Maps

Low

Low

Regular

Regular

High

High

Ksat maps for low, regular, and high values

Hydrologic Group

We can visualize the hydrologic group maps to compare the local and SDA imports. The hydrologic group map is a categorical raster with values of A, B, C, D, and A/D. The color scheme is applied using r.colors with the hydgrp color table.

Hydrologic groups

Hydrologic groups

Calculate Curve Number

Install the r.curvenumber extension if it is not already installed. This extension is available in the GRASS Addons repository and can be installed using g.extension.

To calculate the curve number, we need both a land cover map and a hydrologic soil group map. The land cover map is from the National Land Cover Database (NLCD) 2024 dataset, which was imported as a COG raster. The landcover_source parameter specifies that the land cover map is from the NLCD dataset, which allows the tool to use the appropriate lookup table for curve number values.

Code 
tools.r_import(
    output="nlcd_2024",
    input="/home/coreywhite/Downloads/Annual_NLCD_LndCov_2024_CU_C1V1.cog.tif",
    extent="region",
    resolution="region",
)

NLCD 2024 land cover map:

Code 
with gs.RegionManager(raster="nlcd_2024", w="w-1600", flags="a"):
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="nlcd_2024")
    m.d_legend(
        raster="nlcd_2024",
        title="NLCD 2024",
        flags="",
        at="5,80,2,45",
        fontsize=12,
    )
    m.show()

Calculate the hydrolic curve number using the r.curvenumber extension, which requires both a land cover map and a hydrologic soil group map. The land cover map is from the National Land Cover Database (NLCD) 2024 dataset, which was imported as a COG raster. The landcover_source parameter specifies that the land cover map is from the NLCD dataset, which allows the tool to use the appropriate lookup table for curve number values.

Code 
tools.r_curvenumber(
    landcover="nlcd_2024",
    soil="hydgrp_sda",
    landcover_source="nlcd",
    output="curvenumber",
)
Code 
with gs.RegionManager(raster="curvenumber", w="w-1600", flags="a"):
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="curvenumber")
    m.d_legend(raster="curvenumber", title="Curve Number", flags="s", at="5, 30, 2, 10", fontsize=12)
    m.show()

Calculate Runoff and time of concentration

Let’s also run r.runoff to calculate runoff depth using the curve number method. This tool requires a precipitation raster, which we can create using r.mapcalc with a constant value for precipitation.

We need the flow direction and stream maps to run r.runoff with time of concentration, which can be calculated using the r.watershed tool.

Code 
tools.r_watershed(
    elevation="elevation",
    drainage="flow_direction",
    accumulation="accumulation",
    stream="stream",
    basin="basins",
    threshold=10,
)
Code 
with gs.RegionManager(raster="accumulation", w="w-1600", flags="a"):
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="accumulation")
    m.d_legend(
        raster="accumulation",
        title="Accumulation",
        flags="s",
        at="5, 30, 2, 10",
        fontsize=12,
    )
    m.show()

    m2 = gj.Map(use_region=True)
    m2.d_shade(shade="relief", color="basins")
    m2.d_legend(
        raster="basins",
        title="Basins 50K threshold",
        flags="s",
        at="5, 30, 2, 10",
        fontsize=12,
    )
    m2.show()

Let calculate the time of concentration using the r.timeofconcentration tool, which requires an elevation raster, flow direction raster, and stream raster as inputs. The output is a raster of time of concentration values in hours.

Code 
tools.r_timeofconcentration(
    elevation="elevation",
    direction="flow_direction",
    stream="stream",
    length_min=1,  # minimum length of flow path
    tc="time_concentration",
)
Code 
with gs.RegionManager(raster="elevation", s="s-200", flags="a"):
    tools.r_colors(map="time_concentration", color="byg", flags="e")
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="time_concentration")

    m.d_legend(
        raster="time_concentration",
        at="8, 12, 20, 80",
        title="Time of Concentration (hours)",
        font="Fira Sans Condensed Light",
        fontsize=14,
        flags="s",
    )
    m.d_barscale(
        at=(20, 28),
        font="Fira Sans Condensed Light",
        fontsize=10,
        length=3,
        width_scale=1,
        units="kilometers",
        bgcolor="none",
        style="both_ticks",
        color="#f7f7f7",
        flags="n",
    )
    m.show()

Let’s check the univarite statistics for the time of concentration raster to see the range of values.

Table 1: Time of Concentration — Summary Statistics
Statistic Value
n 112,024.000
null_cells 65,609.000
cells 177,633.000
min 0.002
max 0.901
range 0.899
mean 0.020
mean_of_abs 0.020
stddev 0.059
variance 0.003
coeff_var 289.757
sum 2,265.908
first_quartile 0.005
median 0.008
third_quartile 0.014
percentile_percentile 90.000
percentile_value 0.032

NOAA Atlas Precipitation Frequency

{'request': {'lat': 39.4003765, 'lon': -97.0630433, 'data': 'intensity', 'units': 'english', 'series': 'pds', 'url': 'https://hdsc.nws.noaa.gov/cgi-bin/new/fe_text.csv?lat=39.4003765&lon=-97.0630433&data=intensity&units=english&series=pds'}, 'metadata': {'data_type': 'Precipitation intensity', 'time_series_type': 'Partial duration', 'project_area': 'Midwestern States', 'location_name_(esri_maps)': 'None', 'station_name': 'None', 'latitude': '39.4003765 Degree', 'longitude': '-97.0630433 Degree', 'elevation_(usgs)': 'None None'}, 'table': {'return_periods_years': [1, 2, 5, 10, 25, 50, 100, 200, 500, 1000], 'rows': [{'duration': '5-min', 'values': {'1': 4.84, '2': 5.72, '5': 7.2, '10': 8.46, '25': 10.3, '50': 11.7, '100': 13.1, '200': 14.7, '500': 16.7, '1000': 18.3}}, {'duration': '10-min', 'values': {'1': 3.54, '2': 4.19, '5': 5.27, '10': 6.2, '25': 7.51, '50': 8.55, '100': 9.62, '200': 10.7, '500': 12.3, '1000': 13.4}}, {'duration': '15-min', 'values': {'1': 2.88, '2': 3.4, '5': 4.28, '10': 5.04, '25': 6.1, '50': 6.95, '100': 7.82, '200': 8.73, '500': 9.96, '1000': 10.9}}, {'duration': '30-min', 'values': {'1': 2.05, '2': 2.42, '5': 3.05, '10': 3.59, '25': 4.36, '50': 4.98, '100': 5.61, '200': 6.27, '500': 7.17, '1000': 7.87}}, {'duration': '60-min', 'values': {'1': 1.33, '2': 1.57, '5': 1.97, '10': 2.33, '25': 2.85, '50': 3.28, '100': 3.73, '200': 4.21, '500': 4.88, '1000': 5.4}}, {'duration': '2-hr', 'values': {'1': 0.818, '2': 0.96, '5': 1.21, '10': 1.43, '25': 1.76, '50': 2.04, '100': 2.33, '200': 2.64, '500': 3.08, '1000': 3.44}}, {'duration': '3-hr', 'values': {'1': 0.605, '2': 0.708, '5': 0.892, '10': 1.06, '25': 1.31, '50': 1.52, '100': 1.75, '200': 1.99, '500': 2.34, '1000': 2.62}}, {'duration': '6-hr', 'values': {'1': 0.355, '2': 0.415, '5': 0.523, '10': 0.622, '25': 0.771, '50': 0.896, '100': 1.03, '200': 1.17, '500': 1.38, '1000': 1.55}}, {'duration': '12-hr', 'values': {'1': 0.202, '2': 0.237, '5': 0.299, '10': 0.355, '25': 0.436, '50': 0.504, '100': 0.575, '200': 0.651, '500': 0.757, '1000': 0.843}}, {'duration': '24-hr', 'values': {'1': 0.115, '2': 0.135, '5': 0.169, '10': 0.199, '25': 0.242, '50': 0.278, '100': 0.314, '200': 0.353, '500': 0.406, '1000': 0.449}}, {'duration': '2-day', 'values': {'1': 0.065, '2': 0.076, '5': 0.094, '10': 0.109, '25': 0.132, '50': 0.15, '100': 0.17, '200': 0.19, '500': 0.218, '1000': 0.24}}, {'duration': '3-day', 'values': {'1': 0.047, '2': 0.054, '5': 0.067, '10': 0.078, '25': 0.094, '50': 0.107, '100': 0.12, '200': 0.134, '500': 0.154, '1000': 0.169}}, {'duration': '4-day', 'values': {'1': 0.038, '2': 0.043, '5': 0.053, '10': 0.062, '25': 0.074, '50': 0.084, '100': 0.095, '200': 0.106, '500': 0.12, '1000': 0.132}}, {'duration': '7-day', 'values': {'1': 0.025, '2': 0.029, '5': 0.035, '10': 0.04, '25': 0.048, '50': 0.054, '100': 0.06, '200': 0.067, '500': 0.076, '1000': 0.083}}, {'duration': '10-day', 'values': {'1': 0.02, '2': 0.023, '5': 0.028, '10': 0.032, '25': 0.038, '50': 0.042, '100': 0.047, '200': 0.052, '500': 0.059, '1000': 0.064}}, {'duration': '20-day', 'values': {'1': 0.014, '2': 0.016, '5': 0.019, '10': 0.021, '25': 0.025, '50': 0.028, '100': 0.031, '200': 0.034, '500': 0.038, '1000': 0.041}}, {'duration': '30-day', 'values': {'1': 0.011, '2': 0.013, '5': 0.015, '10': 0.017, '25': 0.02, '50': 0.022, '100': 0.024, '200': 0.027, '500': 0.029, '1000': 0.032}}, {'duration': '45-day', 'values': {'1': 0.009, '2': 0.01, '5': 0.012, '10': 0.014, '25': 0.016, '50': 0.018, '100': 0.019, '200': 0.021, '500': 0.023, '1000': 0.024}}, {'duration': '60-day', 'values': {'1': 0.008, '2': 0.009, '5': 0.011, '10': 0.012, '25': 0.014, '50': 0.015, '100': 0.016, '200': 0.017, '500': 0.019, '1000': 0.02}}]}, 'bound': 'expected'}
duration,1,2,5,10,25,50,100,200,500,1000
5-min,4.84,5.72,7.2,8.46,10.3,11.7,13.1,14.7,16.7,18.3
10-min,3.54,4.19,5.27,6.2,7.51,8.55,9.62,10.7,12.3,13.4
15-min,2.88,3.4,4.28,5.04,6.1,6.95,7.82,8.73,9.96,10.9
30-min,2.05,2.42,3.05,3.59,4.36,4.98,5.61,6.27,7.17,7.87
60-min,1.33,1.57,1.97,2.33,2.85,3.28,3.73,4.21,4.88,5.4
2-hr,0.818,0.96,1.21,1.43,1.76,2.04,2.33,2.64,3.08,3.44
3-hr,0.605,0.708,0.892,1.06,1.31,1.52,1.75,1.99,2.34,2.62
6-hr,0.355,0.415,0.523,0.622,0.771,0.896,1.03,1.17,1.38,1.55
12-hr,0.202,0.237,0.299,0.355,0.436,0.504,0.575,0.651,0.757,0.843
24-hr,0.115,0.135,0.169,0.199,0.242,0.278,0.314,0.353,0.406,0.449
2-day,0.065,0.076,0.094,0.109,0.132,0.15,0.17,0.19,0.218,0.24
3-day,0.047,0.054,0.067,0.078,0.094,0.107,0.12,0.134,0.154,0.169
4-day,0.038,0.043,0.053,0.062,0.074,0.084,0.095,0.106,0.12,0.132
7-day,0.025,0.029,0.035,0.04,0.048,0.054,0.06,0.067,0.076,0.083
10-day,0.02,0.023,0.028,0.032,0.038,0.042,0.047,0.052,0.059,0.064
20-day,0.014,0.016,0.019,0.021,0.025,0.028,0.031,0.034,0.038,0.041
30-day,0.011,0.013,0.015,0.017,0.02,0.022,0.024,0.027,0.029,0.032
45-day,0.009,0.01,0.012,0.014,0.016,0.018,0.019,0.021,0.023,0.024
60-day,0.008,0.009,0.011,0.012,0.014,0.015,0.016,0.017,0.019,0.02

We can use the the time_concentration raster as an input to the r.runoff tool, which calculates runoff depth and peak discharge using the curve number method. The r.runoff tool also requires a precipitation raster, which we can create using r.mapcalc with a constant value for precipitation (\(50\) mm) over a 24 hour period.

Code 
df_tc_min = df_tc[["min", "mean", "max"]].apply(
    lambda x: round(x * 60, 2)
)  # convert to minutes
max_tc = df_tc_min["max"].max()
print(max_tc)

df_rf_intensity = pd.read_csv(io.StringIO(pfds_csv))
display(
    df_rf_intensity.style.hide(axis="index")
    .set_caption("Rainfall Intensity — NOAA Atlas 14, Coweeta, NC")
    .format(
        {
            "Value": lambda v: f"{v:,.3f}"
            if isinstance(v, (int, float, np.number))
            else v
        }
    )
)

# | label: rainfall_intensity_figure
# | echo: false
# | fig-cap: "NOAA Atlas 14 rainfall intensity by storm duration for Coweeta, NC, with the maximum time of concentration marked."

rf = df_rf_intensity.copy()

duration_col = next(
    (c for c in rf.columns if "duration" in c.lower()),
    rf.columns[0],
)

def _duration_to_minutes(value):
    text = str(value).strip().lower()
    match = re.search(r"(\d+(?:\.\d+)?)", text)
    if not match:
        return np.nan
    number = float(match.group(1))
    if "day" in text:
        return number * 24 * 60
    if "hr" in text or "hour" in text:
        return number * 60
    return number

import re

rf["duration_min"] = rf[duration_col].apply(_duration_to_minutes)

value_cols = [c for c in rf.columns if c not in {duration_col, "duration_min"}]
rf_long = rf.melt(
    id_vars=["duration_min"],
    value_vars=value_cols,
    var_name="return_period",
    value_name="intensity_mm_hr",
)
rf_long["intensity_mm_hr"] = pd.to_numeric(rf_long["intensity_mm_hr"], errors="coerce")
rf_long = rf_long.dropna(subset=["duration_min", "intensity_mm_hr"]).sort_values(
    ["return_period", "duration_min"]
)

sns.set_theme(style="whitegrid", context="talk")
fig, ax = plt.subplots(figsize=(10, 6))

sns.lineplot(
    data=rf_long,
    x="duration_min",
    y="intensity_mm_hr",
    hue="return_period",
    marker="o",
    linewidth=2.5,
    ax=ax,
)

ax.axvline(
    max_tc,
    color="black",
    linestyle="--",
    linewidth=1.5,
    label=f"Max Tc = {max_tc:.1f} min",
)

ax.annotate(
    f"Max Tc\n{max_tc:.1f} min",
    xy=(max_tc, rf_long["intensity_mm_hr"].max()),
    xytext=(8, -10),
    textcoords="offset points",
    ha="left",
    va="top",
    fontsize=11,
    bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="0.6", alpha=0.9),
)

ax.set_xscale("log")
ax.set_xlabel("Storm duration (minutes)")
ax.set_ylabel("Rainfall intensity (mm/hr)")
ax.set_title("Rainfall intensity–duration curves")
ax.legend(title="Return period", frameon=True, ncol=2)
sns.despine()
plt.tight_layout()
plt.show()


return_periods = rf_long["return_period"].dropna().unique()
rp_100 = next(
    (
        rp
        for rp in return_periods
        if re.search(r"100", str(rp)) and re.search(r"(yr|year)", str(rp).lower())
    ),
    next((rp for rp in return_periods if re.search(r"100", str(rp))), None),
)

if rp_100 is None:
    raise ValueError("Could not identify the 100-year storm column in rainfall intensity data.")

rf_100 = (
    rf_long.loc[rf_long["return_period"] == rp_100, ["duration_min", "intensity_mm_hr"]]
    .dropna()
    .sort_values("duration_min")
)

storm_intensity_100yr = float(
    np.interp(
        np.log10(max_tc),
        np.log10(rf_100["duration_min"].to_numpy()),
        rf_100["intensity_mm_hr"].to_numpy(),
    )
)

rain_inches_per_hour = storm_intensity_100yr

display(
    pd.DataFrame(
        {
            "return_period": [rp_100],
            "max_tc_min": [round(max_tc, 2)],
            "intensity_in_hr": [round(storm_intensity_100yr, 2)],
        }
    ).style.hide(axis="index").set_caption(
        "Interpolated rainfall intensity at maximum time of concentration"
    )
)

print(f"100-year storm intensity at max Tc ({max_tc:.2f} min): {storm_intensity_100yr:.2f} in/hr")
54.05
Table 2: Rainfall Intensity — NOAA Atlas 14, Coweeta, NC
duration 1 2 5 10 25 50 100 200 500 1000
5-min 4.840000 5.720000 7.200000 8.460000 10.300000 11.700000 13.100000 14.700000 16.700000 18.300000
10-min 3.540000 4.190000 5.270000 6.200000 7.510000 8.550000 9.620000 10.700000 12.300000 13.400000
15-min 2.880000 3.400000 4.280000 5.040000 6.100000 6.950000 7.820000 8.730000 9.960000 10.900000
30-min 2.050000 2.420000 3.050000 3.590000 4.360000 4.980000 5.610000 6.270000 7.170000 7.870000
60-min 1.330000 1.570000 1.970000 2.330000 2.850000 3.280000 3.730000 4.210000 4.880000 5.400000
2-hr 0.818000 0.960000 1.210000 1.430000 1.760000 2.040000 2.330000 2.640000 3.080000 3.440000
3-hr 0.605000 0.708000 0.892000 1.060000 1.310000 1.520000 1.750000 1.990000 2.340000 2.620000
6-hr 0.355000 0.415000 0.523000 0.622000 0.771000 0.896000 1.030000 1.170000 1.380000 1.550000
12-hr 0.202000 0.237000 0.299000 0.355000 0.436000 0.504000 0.575000 0.651000 0.757000 0.843000
24-hr 0.115000 0.135000 0.169000 0.199000 0.242000 0.278000 0.314000 0.353000 0.406000 0.449000
2-day 0.065000 0.076000 0.094000 0.109000 0.132000 0.150000 0.170000 0.190000 0.218000 0.240000
3-day 0.047000 0.054000 0.067000 0.078000 0.094000 0.107000 0.120000 0.134000 0.154000 0.169000
4-day 0.038000 0.043000 0.053000 0.062000 0.074000 0.084000 0.095000 0.106000 0.120000 0.132000
7-day 0.025000 0.029000 0.035000 0.040000 0.048000 0.054000 0.060000 0.067000 0.076000 0.083000
10-day 0.020000 0.023000 0.028000 0.032000 0.038000 0.042000 0.047000 0.052000 0.059000 0.064000
20-day 0.014000 0.016000 0.019000 0.021000 0.025000 0.028000 0.031000 0.034000 0.038000 0.041000
30-day 0.011000 0.013000 0.015000 0.017000 0.020000 0.022000 0.024000 0.027000 0.029000 0.032000
45-day 0.009000 0.010000 0.012000 0.014000 0.016000 0.018000 0.019000 0.021000 0.023000 0.024000
60-day 0.008000 0.009000 0.011000 0.012000 0.014000 0.015000 0.016000 0.017000 0.019000 0.020000

Table 3: Interpolated rainfall intensity at maximum time of concentration
return_period max_tc_min intensity_in_hr
100 54.050000 4.010000 
100-year storm intensity at max Tc (54.05 min): 4.01 in/hr
Code 
duration_min = round(max_tc, 2)
intensity_in_hr =rain_inches_per_hour
intensity_mm_hr = intensity_in_hr * 25.4

duration_hr = duration_min / 60.0
rain_mm = intensity_mm_hr * duration_hr

print(f"Storm duration: {duration_min:.4f} min")
print(f"Storm duration: {duration_hr:.4f} hr")
print(f"Rainfall intensity: {intensity_in_hr:.3f} in/hr")
print(f"Rainfall intensity: {intensity_mm_hr:.3f} mm/hr")
print(f"Rain depth: {rain_mm:.3f} mm")
tools.r_mapcalc(
    expression=f"rain = {rain_mm}",  # mm of precipitation
)
Storm duration: 54.0500 min
Storm duration: 0.9008 hr
Rainfall intensity: 4.013 in/hr
Rainfall intensity: 101.937 mm/hr
Rain depth: 91.828 mm
Code 
# tools.g_extension(extension="r.runoff", url="../../../cwhite911/grass-addons/src/raster/r.runoff")
# Make sure you have the following GRASS addons installed
# - r.accumulate
# - r.runoff

tools.r_runoff(
    rainfall="rain",
    direction="flow_direction",
    duration=duration_hr,  # hours
    curvenumber="curvenumber",
    time_concentration="time_concentration",
    runoff_depth="runoff_depth",
    runoff_volume="runoff_volume",
    upstream_area="upstream_area",
    upstream_runoff_depth="upstream_runoff_depth",
    upstream_runoff_volume="upstream_runoff_volume",
    time_to_peak="time_to_peak",
    peak_discharge="peak_discharge",
)

Let’s visualize the runoff depth

Code 
with gs.RegionManager(raster="runoff_depth", s="s-200", flags="a"):
    tools.r_colors(map="runoff_depth", color="water")
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="runoff_depth")
    m.d_legend(
        raster="runoff_depth", at="5, 30, 2, 10", title="Runoff Depth (mm)", flags="b"
    )
    m.show()

Peak discharge and time to peak discharge maps

Peak Discharge (m³/s)

Peak Discharge (m³/s)

Time to Peak (hours)

Time to Peak (hours)

Discharge

Manning’s n Map

We will calculate the Manning’s n map using the r.manning extension, which uses the National Land Cover Dataset (2024) map to assign Manning’s n values based on the lookup table from Kalyanapu et al. (2009). The output is a raster of Manning’s n values that can be used as an input to the SIMWE model.

Now we can calculate the roughness map using the r.manning tool.

Code 
tools.r_manning(
    input="nlcd_2024",
    output="mannings_n",
    landcover="nlcd",
    source="kalyanapu",
    seed=42,
)


with gs.RegionManager(raster="mannings_n", s="s-200", flags="a"):
    tools.r_colors(map="mannings_n", color="plasma", flags="")
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="mannings_n")

    m.d_legend(
        raster="mannings_n",
        at="8, 12, 20, 80",
        title="Manning's n",
        font="Fira Sans Condensed Light",
        fontsize=14,
        flags="",
    )
    m.d_barscale(
        at=(20, 28),
        font="Fira Sans Condensed Light",
        fontsize=10,
        length=3,
        width_scale=1,
        units="kilometers",
        bgcolor="none",
        style="both_ticks",
        color="#f7f7f7",
        flags="n",
    )
    m.show()

Run SIMWE model

Calcuate the slope and aspect from the elevation raster, which are required inputs for the r.sim.water tool.

Precipitation and Rainfall Excess

We will use a constant rainfall value for the simulation, which is based on an 100 year storm from NOAA Atlas 14 for Clay Center, KS. The rainfall excess is calculated as the runoff depth (mm), dervied using the SCS-CN method, divided by the 24 hours to get a rainfall excess in mm/hr.

Location Information Value
Name Clay Center, Kansas, USA*
Station name
Site ID
Latitude 39.4003765°
Longitude -97.0630433°
Elevation ft**

Rainfall excess map

Code 
# Convert runoff_depth (mm) to rainfall excess. during a 54 min 100 year storm. So the runoff depth (m) is divided by 0.9 hours (54 min) to get the rainfall excess in mm/hr. The 0.9 hours is derived from the duration of the storm, which is 54 minutes, converted to hours (54/60 = 0.9). This gives us the average rainfall excess rate in mm/hr over the duration of the storm.
# Runoff depth per cell [mm] is calculated by the SCS-CN method.


tools.r_mapcalc(
    expression=f"rainfall_excess = runoff_depth / 0.9",
)


with gs.RegionManager(raster="rainfall_excess", s="s-200", flags="a"):
    tools.r_colors(map="rainfall_excess", color="blue", flags="e")
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="rainfall_excess")

    m.d_legend(
        raster="rainfall_excess",
        title="Rainfall Excess (mm/hr)",
        at="8, 12, 20, 80",
        font="Fira Sans Condensed Light",
        fontsize=14,
        flags="",
    )
    m.d_barscale(
        at=(20, 28),
        font="Fira Sans Condensed Light",
        fontsize=10,
        length=3,
        width_scale=1,
        units="kilometers",
        bgcolor="none",
        style="both_ticks",
        color="#f7f7f7",
        flags="n",
    )
    m.show()

Calcuate the overland flow for a 47 min 100 year storm using the r.sim.water tool, which simulates overland flow using a particle-based approach. The tool requires the elevation raster, rainfall excess raster, and Manning’s n raster as inputs, along with the number of walkers to simulate and the number of iterations to run. The output is a series of depth and discharge rasters that can be visualized as maps or animations.

Code 
# 100yr storm running for 54.05 min with soil informed
# Rainfall excess
duration_min = 54.05
tools.r_slope_aspect(elevation="elevation", dx="dx", dy="dy")
tools.r_sim_water(
    elevation="elevation",
    depth="rx_100yr_depth",
    discharge="rx_100yr_discharge",
    walkers_output="rx_100yr_walkers",
    rain="rainfall_excess",  # Infiltration / Veg Intercept
    man="mannings_n",  # Mixed forest Kalyanapu et al. (2009)
    nwalkers=1000000,
    # infil="ksat_r_sda", # Infiltration of flowing water
    # niterations=1440 / 2,
    duration=duration_min, # Duration of the simulated water flow [minutes]
    output_step=2, # Time interval for creating output maps [minutes]
    nprocs=30,
    flags="t",
)

# Register the output maps into a space time dataset
gs.run_command(
    "t.create",
    output="rx_100yr_depth_sum",
    type="strds",
    temporaltype="absolute",
    title="Runoff Depth",
    description="Runoff Depth in [m]",
    overwrite=True,
)

# Get the list of depth maps
depth_list = gs.read_command(
    "g.list", type="raster", pattern="rx_100yr_depth.*", separator="comma"
).strip()

# Register the maps
gs.run_command(
    "t.register",
    input="rx_100yr_depth_sum",
    type="raster",
    start="2024-01-01",
    increment="2 minutes",
    maps=depth_list,
    flags="i",
    overwrite=True,
)
Code 
tools.t_info(input="rx_100yr_depth_sum")
ToolResult(returncode=0, stdout=' +-------------------- Space Time Raster Dataset -----------------------------+
 |                                                                            |
 +-------------------- Basic information -------------------------------------+
 | Id: ........................ rx_100yr_depth_sum@ssurgo_demo
 | Name: ...................... rx_100yr_depth_sum
 | Mapset: .................... ssurgo_demo
 | Creator: ................... coreywhite
 | Temporal type: ............. absolute
 | Creation time: ............. 2026-04-02 15:30:01
 | Modification time:.......... 2026-04-02 15:30:05
 | Semantic type:.............. mean
 +-------------------- Absolute time -----------------------------------------+
 | Start time:................. 2024-01-01 00:00:00
 | End time:................... 2024-01-01 00:54:00
 | Granularity:................ 2 minutes
 | Temporal type of maps:...... interval
 +-------------------- Spatial extent ----------------------------------------+
 | North:...................... 4363578.815536
 | South:...................... 4362417.815536
 | East:.. .................... 667469.808073
 | West:....................... 666092.808073
 | Top:........................ 0.0
 | Bottom:..................... 0.0
 +-------------------- Metadata information ----------------------------------+
 | Raster register table:...... raster_map_register_17f8bf055ce04af0b677eba9da5f19e6
 | North-South resolution min:. 3.0
 | North-South resolution max:. 3.0
 | East-west resolution min:... 3.0
 | East-west resolution max:... 3.0
 | Minimum value min:.......... 0.000401
 | Minimum value max:.......... 0.000512
 | Maximum value min:.......... 0.227084
 | Maximum value max:.......... 0.959598
 | Aggregation type:........... None
 | Number of semantic labels:.. 0
 | Semantic labels:............ None
 | Number of registered maps:.. 27
 |
 | Title:
 | Runoff Depth
 | Description:
 | Runoff Depth in [m]
 | Command history:
 | # 2026-04-02 15:30:01 
 | t.create --o output="rx_100yr_depth_sum" type="strds"
 |     temporaltype="absolute" title="Runoff Depth"
 |     description="Runoff Depth in [m]"
 | # 2026-04-02 15:30:05 
 | t.register --o -i input="rx_100yr_depth_sum"
 |     type="raster" start="2024-01-01" increment="2 minutes"
 |     maps="rx_100yr_depth.02,rx_100yr_depth.04,rx_100yr_depth.06,rx_100yr_depth.08,rx_100yr_depth.10,rx_100yr_depth.12,rx_100yr_depth.14,rx_100yr_depth.16,rx_100yr_depth.18,rx_100yr_depth.20,rx_100yr_depth.22,rx_100yr_depth.24,rx_100yr_depth.26,rx_100yr_depth.27,rx_100yr_depth.29,rx_100yr_depth.31,rx_100yr_depth.33,rx_100yr_depth.35,rx_100yr_depth.37,rx_100yr_depth.39,rx_100yr_depth.41,rx_100yr_depth.43,rx_100yr_depth.45,rx_100yr_depth.47,rx_100yr_depth.49,rx_100yr_depth.51,rx_100yr_depth.53"
 | 
 +----------------------------------------------------------------------------+
', stderr='')
Code 
depth_list = gs.read_command(
    "g.list", type="raster", pattern="rx_100yr_depth.*", separator="comma"
).strip()

max_depth = depth_list.split(",")[-1]
stats = tools.r_info(map=max_depth, format="json").json
with gs.RegionManager(raster="accumulation", s="s-200", flags="a"):
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="naip_2021_rgb@naip")
    m.d_rast(map=max_depth, values=f"0.005-{stats['max']}")
    m.d_legend(
        raster=max_depth,
        title="Deqpth [m]",
        at="8, 12, 15, 85",
        font="Fira Sans Condensed Light",
        range=[0.005, stats["max"]],
        fontsize=14,
        flags="",
    )
    m.d_barscale(
        at=(20, 28),
        font="Fira Sans Condensed Light",
        fontsize=10,
        length=3,
        width_scale=1,
        units="kilometers",
        bgcolor="none",
        style="both_ticks",
        color="#f7f7f7",
        flags="n",
    )
    m.show()

Animation of depth maps

Code 
max_depth = depth_list.split(",")[-1]
stats = tools.r_info(map=max_depth, format="json").json
m = gj.TimeSeriesMap(use_region=True, width=800)
m.d_rast(map="naip_2021_rgb@naip")
m.add_raster_series("rx_100yr_depth_sum", values=f"0.005-{stats['max']}")
m.d_legend(
    raster=max_depth,
    title="Depth [m]",
    flags="bt",
    border_color="none",
    range=[0.005, stats["max"]],
    at=(1, 50, 0, 5),
)

m.save("output/clay_center_rx_100yr_simwe_depth_animation.gif")
m.show()

Overland Flow Simulation with Soil Data

Overland Flow Simulation with Soil Data

No soils

Code 
intensity_mm_hr = 101.937
tools.r_sim_water(
    elevation="elevation",
    depth="basic_depth",
    discharge="basic_disch",
    walkers_output="basic_particles",
    rain_value=intensity_mm_hr,  # Rainfall intensity: 95.584 mm/hr
    man="mannings_n",  # Mixed forest Kalyanapu et al. (2009)
    nwalkers=1000000,
    # infil="ksat_r_sda", # Infiltration of flowing water
    niterations=duration_min,
    output_step=2,
    nprocs=30,
    flags="t",
)

gs.run_command(
    "t.create",
    output="basic_depth_sum",
    type="strds",
    temporaltype="absolute",
    title="Runoff Depth",
    description="Runoff Depth in [m]",
    overwrite=True,
)

# Get the list of depth maps
depth_list = gs.read_command(
    "g.list", type="raster", pattern="basic_depth.*", separator="comma"
).strip()

# Register the maps
gs.run_command(
    "t.register",
    input="basic_depth_sum",
    type="raster",
    start="2024-01-01",
    increment="2 minutes",
    maps=depth_list,
    flags="i",
    overwrite=True,
)
Code 
tools.t_info(input="basic_depth_sum")
ToolResult(returncode=0, stdout=' +-------------------- Space Time Raster Dataset -----------------------------+
 |                                                                            |
 +-------------------- Basic information -------------------------------------+
 | Id: ........................ basic_depth_sum@ssurgo_demo
 | Name: ...................... basic_depth_sum
 | Mapset: .................... ssurgo_demo
 | Creator: ................... coreywhite
 | Temporal type: ............. absolute
 | Creation time: ............. 2026-04-02 16:01:15
 | Modification time:.......... 2026-04-02 16:01:18
 | Semantic type:.............. mean
 +-------------------- Absolute time -----------------------------------------+
 | Start time:................. 2024-01-01 00:00:00
 | End time:................... 2024-01-01 00:54:00
 | Granularity:................ 2 minutes
 | Temporal type of maps:...... interval
 +-------------------- Spatial extent ----------------------------------------+
 | North:...................... 4363578.815536
 | South:...................... 4362417.815536
 | East:.. .................... 667469.808073
 | West:....................... 666092.808073
 | Top:........................ 0.0
 | Bottom:..................... 0.0
 +-------------------- Metadata information ----------------------------------+
 | Raster register table:...... raster_map_register_b0cb9d6dc0084c97abe2c2afa447039b
 | North-South resolution min:. 3.0
 | North-South resolution max:. 3.0
 | East-west resolution min:... 3.0
 | East-west resolution max:... 3.0
 | Minimum value min:.......... 0.000975
 | Minimum value max:.......... 0.000975
 | Maximum value min:.......... 0.252924
 | Maximum value max:.......... 1.342173
 | Aggregation type:........... None
 | Number of semantic labels:.. 0
 | Semantic labels:............ None
 | Number of registered maps:.. 27
 |
 | Title:
 | Runoff Depth
 | Description:
 | Runoff Depth in [m]
 | Command history:
 | # 2026-04-02 16:01:15 
 | t.create --o output="basic_depth_sum" type="strds"
 |     temporaltype="absolute" title="Runoff Depth"
 |     description="Runoff Depth in [m]"
 | # 2026-04-02 16:01:18 
 | t.register --o -i input="basic_depth_sum"
 |     type="raster" start="2024-01-01" increment="2 minutes"
 |     maps="basic_depth.02,basic_depth.04,basic_depth.06,basic_depth.08,basic_depth.10,basic_depth.12,basic_depth.14,basic_depth.16,basic_depth.18,basic_depth.20,basic_depth.22,basic_depth.24,basic_depth.26,basic_depth.28,basic_depth.30,basic_depth.32,basic_depth.34,basic_depth.36,basic_depth.38,basic_depth.40,basic_depth.42,basic_depth.44,basic_depth.46,basic_depth.48,basic_depth.50,basic_depth.52,basic_depth.54"
 | 
 +----------------------------------------------------------------------------+
', stderr='')
Code 
depth_list = gs.read_command(
    "g.list", type="raster", pattern="basic_depth.*", separator="comma"
).strip()

basic_max_depth = depth_list.split(",")[-1]
stats = tools.r_info(map=basic_max_depth, format="json").json
with gs.RegionManager(raster="accumulation", s="s-200", flags="a"):
    m = gj.Map(use_region=True)
    m.d_shade(shade="relief", color="naip_2021_rgb@naip")
    m.d_rast(map=basic_max_depth, values=f"0.005-{stats['max']}")
    m.d_legend(
        raster=basic_max_depth,
        title="Depth [m]",
        flags="btl",
        border_color="none",
        range=[0.005, stats["max"]],
        at=(5, 50, 2, 5),
    )
    m.show()

no soil animation

Code 
basic_max_depth = depth_list.split(",")[-1]
stats = tools.r_info(map=basic_max_depth, format="json").json
m = gj.TimeSeriesMap(use_region=True, width=800)
m.d_rast(map="naip_2021_rgb@naip")
m.add_raster_series("basic_depth_sum", values=f"0.005-{stats['max']}")
m.d_legend(
    raster=basic_max_depth,
    title="Depth [m]",
    flags="bt",
    border_color="none",
    range=[0.005, stats["max"]],
    at=(1, 50, 0, 5),
)
m.save("output/clay_center_basic_depth_simwe_depth_animation.gif")
m.show()

Overland Flow Simulation without Soil Data

Overland Flow Simulation without Soil Data