SoakSIM sizes soakwell infiltration systems to Australian standards by coupling the ILSAX hydrological model (DRAINS Manual, 2018) with ARR 2019 design inputs and a mass-balance soakwell routing engine consistent with the Stormwater Management Manual for Western Australia (Argue, 2004).
Each storm duration and AEP combination is simulated for every ARR temporal pattern. The adopted design pattern is selected using the ARR probability-neutral method (Section 9). The critical duration is the one producing the greatest peak soakwell depth.
Design rainfall depths are sourced from the Bureau of Meteorology Intensity-Frequency-Duration (IFD) Design Rainfall Data System (Revised 2016, based on ARR 2019). SoakSIM scrapes the official BOM web interface at:
http://www.bom.gov.au/water/designRainfalls/revised-ifd/
For a given site coordinate (latitude, longitude), the system retrieves a matrix of depth values (mm) across:
The depth $P$ (mm) for a given duration $D$ and AEP is the total rainfall expected to fall in that duration with that annual probability of exceedance.
ARR 2019 replaces the old single design temporal pattern with an ensemble of 10 temporal patterns per duration. These are expressed as cumulative fraction arrays $\{F_1, F_2, \ldots, F_n\}$ where $F_n = 1.0$.
Patterns are fetched from the ARR Data Hub:
https://data.arr-software.org/
Each pattern represents one historically observed way that rainfall can be distributed across the storm duration. Using all 10 patterns avoids bias toward any single distribution shape.
Each temporal pattern is combined with the BOM design depth to create a time-varying rainfall hyetograph.
The result is a depth series $\{p_1, p_2, \ldots, p_N\}$ in mm, where $N = D / \Delta t$.
SoakSIM uses the ILSAX hydrological model as described in the DRAINS User Manual (June 2018, Section 8.3.2). Three surface types are modelled independently:
| Surface | Loss Method | Description |
|---|---|---|
| Paved (DCIA) | Depression storage only | Directly-connected impervious — roofs, sealed driveways. Only a small initial abstraction (typically 1 mm). |
| Supplementary | Depression storage only | Indirectly-connected impervious areas that drain to pervious surfaces first. |
| Grassed (pervious) | Horton infiltration + depression storage | Lawns, gardens. Uses full Horton infiltration capacity curve. |
Depression storage acts as a bucket that must fill before any runoff occurs:
$$p_j^{\text{excess}} = \max\left(0,\; p_j - S_{\text{remaining}}\right)$$where $S_{\text{remaining}}$ decreases as rainfall fills the depression stores. Typical values: 1 mm (paved), 1 mm (supplementary), 5 mm (grassed).
The Horton equation describes the soil's time-varying infiltration capacity:
$$f(t) = f_c + (f_0 - f_c)\,e^{-k\,t}$$where:
| Soil Type | $f_0$ (mm/hr) | $f_c$ (mm/hr) |
|---|---|---|
| A (Sandy) | 250 | 25 |
| B (Sandy Loam) | 200 | 13 |
| C (Clay Loam) | 125 | 6 |
| D (Heavy Clay) | 75 | 3 |
The Antecedent Moisture Condition (AMC 1–4) adjusts the starting infiltration rate. Higher AMC = wetter soil = lower starting capacity.
Rather than iterating over the Horton equation, SoakSIM uses Watson's method which splits infiltration into a diminishing component and a constant component:
$$\Delta F_{\text{cap}} = \left(1 - e^{-k\,\Delta t}\right)\left(F_{\text{max}}^{\text{dim}} - F_d\right) + f_c \cdot \Delta t$$where $F_{\text{max}}^{\text{dim}} = (f_0 - f_c) / k$ is the total diminishing infiltration capacity, and $F_d$ is the accumulated diminishing infiltration so far. The actual infiltration for each timestep is:
$$f_{\text{actual}} = \min\left(p_j,\; \Delta F_{\text{cap}}\right)$$Excess rainfall (runoff) from the grassed surface for each timestep is:
$$p_j^{\text{excess}} = p_j - f_{\text{actual}}$$Each surface type has a travel time computed from the kinematic wave equation (Ragan & Duru, 1972; DRAINS Eq. 8.3):
$$t_e = 6.94 \,\frac{(L \cdot n^*)^{0.6}}{I^{0.4} \cdot S^{0.3}}$$where:
Total time of entry = $t_e$ + any user-specified additional travel time.
A linear time-area histogram represents the proportion of catchment area contributing runoff at each timestep. The runoff hydrograph is generated by convolving excess rainfall with this histogram:
$$Q_n = \frac{A}{360}\sum_{i} h_i \cdot I_{n-i}$$where:
The three surface discharges (paved, supplementary, grassed) are summed to give the total catchment hydrograph.
SoakSIM routes the runoff hydrograph through the soakwell system using a forward-stepping mass balance at each computational timestep, consistent with the principles of the Stormwater Management Manual for WA (Argue, 2004).
where $S$ = stored volume (m³), $Q_{\text{in}}$ = inflow rate (m³/s), $Q_{\text{infil}}$ = infiltration rate (m³/s).
Concrete soakwells are modelled as hollow cylinders (void ratio = 1.0):
$$V_{\text{storage}} = \frac{\pi}{4}\,d^2\,H$$where $d$ = internal diameter (m) and $H$ = internal height (m).
The base slab is not fully open — it has a circular hole whose diameter is 300 mm less than the soakwell diameter:
$$A_{\text{base}} = \frac{\pi}{4}\,(d - 0.3)^2$$Side-wall infiltration occurs through louvre slots cast into the concrete (not the full wall surface). Each louvre is 150 mm wide × 50 mm high, arranged on a 300 mm × 300 mm grid. The open-area ratio of the wall is:
$$\alpha = \frac{150 \times 50}{300 \times 300} = \frac{7500}{90000} \approx 8.33\%$$The effective side infiltration area is:
$$A_{\text{side}} = \alpha \times \pi\,d\,H = 0.0833 \times \pi\,d\,H$$Infiltration at each timestep depends on the water depth within the soakwell:
$$Q_{\text{infil}} = k_h \cdot R_f \cdot A_{\text{active}}$$where:
where $h$ = current water depth. When the soakwell is full, $h = H$, so the entire side wall is active. When nearly empty, only the base and a small side wall portion contribute.
If stored water exceeds the soakwell capacity ($S > V_{\text{storage}}$), the excess ponds on the surface with a much larger plan area ($\approx 50 \times A_{\text{base}}$). All side walls remain fully submerged, and water drains back into the soakwell as it empties.
After the storm ends ($Q_{\text{in}} = 0$), the mass balance continues until $S \approx 0$. The total time from storm end to empty is the numerical drain time.
When the user provides a Depth to Groundwater input, SoakSIM automatically selects the most appropriate routing model:
| Depth to groundwater | Model selected | Reason |
|---|---|---|
| Not provided | Standard mass-balance (Step 6) | No groundwater data |
| > 2.0 m | Standard mass-balance (Step 6) | GW coupling negligible |
| ≤ 2.0 m and soakwell base above GWT | Roldin coupled model (this section) | GW mound affects performance |
| ≤ 2.0 m and soakwell base below GWT | Roldin coupled model with submersion penalty | Soakwell partially submerged |
The selected model and the reason for selection are displayed in the results panel.
When a shallow water table is present, the infiltration rate from a soakwell is progressively reduced as the local groundwater mound rises toward the soakwell base. The model is based on the simplified soakaway approach of Roldin et al. (2013), which couples a soakwell mass balance with the Hantush (1967) strip-basin mounding equation and Van Genuchten (1980) soil moisture retention. Roldin et al. (2013) demonstrated that this approach matches a full 2D Richards-equation model to within ~5% at approximately 600× less computational cost.
SoakSIM uses soil parameters from the Carsel & Parrish (1988) Van Genuchten catalogue. The user selects a soil texture class; SoakSIM looks up the corresponding α, n, θr, θs and Sy values.
| Soil Texture | θr | θs | α (m⁻¹) | n | Sy |
|---|---|---|---|---|---|
| Sand | 0.045 | 0.430 | 14.50 | 2.68 | 0.27 |
| Loamy Sand | 0.057 | 0.410 | 12.40 | 2.28 | 0.25 |
| Sandy Loam | 0.065 | 0.410 | 7.50 | 1.89 | 0.20 |
| Loam | 0.078 | 0.430 | 3.60 | 1.56 | 0.13 |
| Silt Loam | 0.067 | 0.450 | 2.00 | 1.41 | 0.13 |
| Sandy Clay Loam | 0.100 | 0.390 | 5.90 | 1.48 | 0.13 |
| Clay Loam | 0.095 | 0.410 | 1.90 | 1.31 | 0.08 |
| Clay | 0.068 | 0.380 | 0.80 | 1.09 | 0.03 |
The Van Genuchten (1980) water content at suction head $h$ (m, positive = unsaturated) is:
$$\theta(h) = \theta_r + (\theta_s - \theta_r)\,S_e(h), \qquad S_e(h) = \left[1 + (\alpha h)^n\right]^{-m}, \quad m = 1 - \tfrac{1}{n}$$At each timestep, SoakSIM computes the current vertical separation between the soakwell base and the top of the groundwater mound:
$$\Delta h(t) = \max\!\bigl(0,\; D_{gw} - h_{\text{mound}}(t)\bigr)$$where $D_{gw} = $ depth to groundwater minus soakwell height (clearance from soakwell base to static GWT, m) and $h_{\text{mound}}(t)$ is the current mound height above the static GWT (m).
A three-point average moisture content is used to represent the soil column above the mound (following Roldin et al., 2013):
$$\bar{\theta}(\Delta h) = \frac{\theta(\Delta h) + \theta\!\left(\tfrac{\Delta h}{2}\right) + \theta_s}{3}$$The average drainable porosity (specific yield) of this column is:
$$S_{y,\text{eff}}(\Delta h) = \max\!\bigl(0.01,\; \theta_s - \bar{\theta}(\Delta h)\bigr)$$The infiltration reduction factor is the ratio of the current moisture deficit to the full specific yield:
$$f_{\text{mound}}(t) = \frac{S_{y,\text{eff}}\!\bigl(\Delta h(t)\bigr)}{S_y}$$This factor captures the physical coupling: when the soil column above the mound is nearly dry (deep mound, large $\Delta h$), $S_{y,\text{eff}} \approx S_y$ and $f_{\text{mound}} \approx 1$ — full infiltration capacity. As the mound rises toward the soakwell base, the soil saturates, $S_{y,\text{eff}} \to 0$, and infiltration is progressively suppressed.
The effective infiltration rate at each timestep is:
$$Q_{\text{infil}}(t) = k_h \cdot R_f \cdot f_{\text{mound}}(t) \cdot A_{\text{active}}(t)$$Infiltrated water does not immediately reach the water table — it must first travel through the vadose zone. SoakSIM estimates the travel time as:
$$\tau_{\text{lag}} = \frac{D_{gw}}{k_h\,/\,(\theta_s - \theta(D_{gw}))} \qquad \text{(if } D_{gw} > 0.1\text{ m)}$$The numerator $D_{gw}$ is the travel distance; the denominator is the piston velocity of the wetting front. Infiltrated volumes are held in a buffer and only credited as groundwater recharge after $\tau_{\text{lag}}$ seconds have elapsed.
Once recharge reaches the water table it forms a mound that spreads laterally according to the Hantush (1967) strip-basin equation adapted for a circular soakwell. The soakwell base is approximated as a strip with equivalent half-width:
$$B = \sqrt{A_{\text{cross-section}} / \pi}$$For a recharge event $j$ with flux $w_j$ (m/s, spread over footprint area $A_{\text{footprint}}$) starting at time $t_j$ and lasting $\Delta t_j$, the contribution to the squared water-table height at time $t$ is:
$$\Delta\!\left(h^2\right)_j = \frac{w_j}{k_h}\!\left[F(m_c,\; t - t_j,\; B) - F\!\left(m_c,\; \max(0,\, t - t_j - \Delta t_j),\; B\right)\right]$$The Hantush $F$-function and its constituent Hantush $S$-function are:
$$F(\tau) = 2\,m_c\,\tau\;\cdot S\!\left(\frac{B}{\sqrt{4\,m_c\,\tau}}\right)$$ $$S(a) = \bigl(1 - 2a^2\bigr)\,\mathrm{erf}(a) + 2a^2 - \frac{2a}{\sqrt{\pi}}\,e^{-a^2}$$The quasi-static transmissivity factor $m_c$ is updated each timestep:
$$m_c(t) = k_h \cdot S_{y,\text{eff}}\!\bigl(\Delta h(t)\bigr) \cdot \bar{h}_{gw}(t), \qquad \bar{h}_{gw}(t) = h_{gw,0} + \tfrac{1}{2}\,h_{\text{mound}}(t)$$where $h_{gw,0}$ = initial height of water table above the impermeable aquifer base (m).
The total mound height at timestep $t$ is computed by superposing all recharge events from the full history:
$$h^2_{\text{mound,new}}(t) = h^2_{gw,0} + \sum_j \Delta\!\left(h^2\right)_j(t)$$ $$h_{\text{mound}}(t) = \min\!\left(\max\!\left(0,\;\sqrt{h^2_{\text{mound,new}}} - h_{gw,0}\right),\; D_{gw}\right)$$The mound is capped at $D_{gw}$ — the physical clearance between the soakwell base and the static water table — since it cannot rise into the soakwell itself without triggering submersion effects.
The superposition sum over the full recharge history is evaluated at every timestep. Naively this is an O(n²) computation (n² calls to erf for a simulation with n timesteps). SoakSIM vectorises the entire history array using NumPy/SciPy so that each timestep requires a single batched erf call, giving a ~30–50× speedup over the scalar loop. For an extended simulation (> 7 days for sandy soils), an age-based pruning step removes events whose Hantush argument has converged to within 0.5% of its asymptote.
When the soakwell base is at or below the static water table ($D_{gw} \leq 0$), SoakSIM applies an additional area fraction penalty to account for the portion of the soakwell that is submerged in the saturated zone and cannot contribute to infiltration:
$$f_{\text{sub}} = \min\!\left(1.0,\;\frac{d_{\text{GWT}}}{H_{sw}}\right)$$where $d_{\text{GWT}}$ is the depth to groundwater (from ground surface) and $H_{sw}$ is the soakwell height. The infiltration rate becomes $Q_{\text{infil}} = k_h \cdot R_f \cdot f_{\text{mound}} \cdot f_{\text{sub}} \cdot A_{\text{active}}$.
The analytical emptying time for a circular soakwell (Argue, 2004; Stormwater Management Manual for WA) is:
where:
In soakwell sizing, a Soil Moderation Factor ($U$) is applied to field-measured hydraulic conductivity to account for the difference between point-scale test results and actual areal conductivity (Engineers Australia, 2006).
| Soil Type | Factor $U$ | Rationale |
|---|---|---|
| Sand | 0.5 | Field tests overestimate areal $k_h$ → halve the tested rate |
| Sandy Clay | 1.0 | No correction needed |
| Medium / Heavy Clay | 2.0 | Field tests underestimate areal $k_h$ → double the tested rate |
SoakSIM applies this factor as the Soil Moderation Factor input. The user-entered infiltration rate should be the raw field-measured rate; SoakSIM then applies $U$ internally.
ARR 2019 provides 10 temporal patterns per duration. To avoid arbitrarily selecting the worst or best case, SoakSIM adopts the probability-neutral approach:
This ensures the design is based on a pattern that is neither the most extreme nor the most benign — consistent with ARR guidance that no individual temporal pattern has a defined probability.
The critical duration is then the adopted pattern (across all durations) that produces the greatest peak soakwell depth.
When the solver mode is enabled, SoakSIM automatically determines the number and size of soakwells required. The procedure is:
For each candidate soakwell size, compute an initial unit count:
$$n = \left\lceil \frac{\forall}{V_{\text{sw}} / U + I_{\text{unit}} \cdot D} \right\rceil$$where:
The formula estimate is then verified by routing:
| Mode | Description |
|---|---|
| Fewest units | Minimise total soakwell count across all catalogue sizes |
| Fixed height | Only consider sizes up to a specified height |
| Fixed diameter | Only consider sizes up to a specified diameter |
Size a soakwell system for a 500 m² (0.05 ha) residential roof in Perth (−31.95°, 115.86°) for the 5% AEP (20-year ARI), 60-minute storm duration. The site has sandy soil with a field-tested hydraulic conductivity of 1.2 m/day (50 mm/hr).
From BOM IFD for Perth at 5% AEP, 60-minute duration:
$$P = 33.5 \;\text{mm}$$ARR provides 10 temporal patterns for the 60-min duration. One example pattern (cumulative fractions at 5-minute intervals) might be:
[0.03, 0.07, 0.12, 0.18, 0.28, 0.52, 0.72, 0.84, 0.92, 0.96, 0.99, 1.00]
The incremental depths (mm per 5-min interval) are these fractions multiplied by 33.5 mm, then differenced:
[1.0, 1.3, 1.7, 2.0, 3.4, 8.0, 6.7, 4.0, 2.7, 1.3, 1.0, 0.3]
Note the characteristic "burst" in the middle of the storm.
For a fully paved (roof) catchment with depression storage = 1 mm:
Total runoff volume (time-area routing):
$$V = \frac{A \times \text{excess depth}}{1000} = \frac{500 \times 32.5}{1000} = 16.25 \;\text{m}^3$$The time-area convolution distributes this over the 60-minute storm producing a hydrograph with a peak of approximately 0.0045 m³/s.
Using the Argue (2004) sizing formula with a 1200 × 1200 mm soakwell ($d = H = 1.2$ m):
Convert $k_h$: $1.2$ m/day $= 1.389 \times 10^{-5}$ m/s
Convert $\tau$: storm time-base $= 60$ min
Soil moderation factor for sand: $U = 0.5$
$$d_{\text{required}} = \sqrt{\frac{\forall}{\frac{\pi}{4}\left(H + 120 \cdot k_h \cdot \tau \cdot U\right)}}$$ $$d_{\text{required}} = \sqrt{\frac{16.25}{\frac{\pi}{4}\left(1.2 + 120 \times 1.389 \times 10^{-5} \times 60 \times 0.5\right)}}$$ $$= \sqrt{\frac{16.25}{0.785 \times (1.2 + 0.05)}} = \sqrt{\frac{16.25}{0.981}} = \sqrt{16.56} = 4.07 \;\text{m}$$This means a single soakwell of diameter 4.07 m would be needed — far larger than standard products. Hence, multiple standard soakwells are required.
SoakSIM's optimiser selects from the catalogue. For 1200 × 1200 mm soakwells:
Initial estimate (net-storage formula):
$$n = \left\lceil \frac{16.25}{1.357/0.5 + 0.226 \times 1.0} \right\rceil = \left\lceil \frac{16.25}{2.714 + 0.226} \right\rceil = \left\lceil \frac{16.25}{2.940} \right\rceil = \lceil 5.53 \rceil = 6$$SoakSIM then routes all 10 temporal patterns through 6 × 1200 × 1200 soakwells, ranks by peak depth, and adopts the 4th highest. If any adopted pattern spills, it increments to 7, then 8, etc. until no spill occurs.
For a single 1200 × 1200 soakwell, check the emptying time:
$$T = -\frac{4.6 \times 1.2}{4 \times 1.389 \times 10^{-5}} \;\log_{10}\!\left(\frac{1.2/4}{1.2 + 1.2/4}\right)$$ $$T = -\frac{5.52}{5.556 \times 10^{-5}} \times \log_{10}\!\left(\frac{0.3}{1.5}\right) = -99{,}352 \times (-0.699)$$ $$T = 69{,}447 \;\text{s} = 19.3 \;\text{hours}$$| Parameter | Value |
|---|---|
| Design storm | 5% AEP, 60-minute duration |
| Catchment area | 500 m² (0.05 ha), 100% paved |
| Design rainfall depth | 33.5 mm |
| Runoff volume | ~16.25 m³ |
| Soakwell size | 1200 × 1200 mm |
| Number required (formula) | 6 |
| Total storage | 6 × 1.357 = 8.14 m³ |
| Infiltration credit (1 hr) | 6 × 0.226 = 1.36 m³ |
| Argue drain time | ~19.3 hours (< 24 hr ✓) |
| Soil moderation factor | 0.5 (sand) |
SoakSIM v1.0 — © 2024–2026 Innealta Engineering