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}}$.
When Infiltration Basin mode is selected, SoakSIM uses the basin geometry and soil parameters entered in the UI rather than the soakwell catalogue. Basin storage is computed from the trapezoidal cross-section defined by the base dimensions, side slopes, and maximum depth.
Let $L_b$ and $W_b$ be the basin base length and width (m), $m$ be the side slope ratio ($1:m$), and $D_b$ be the maximum basin depth (m). The plan dimensions at depth $z$ are:
$$L(z) = L_b + 2mz, \qquad W(z) = W_b + 2mz$$The plan area at depth $z$ is therefore:
$$A(z) = L(z)\,W(z) = (L_b + 2mz)(W_b + 2mz)$$The storage volume up to depth $z$ is the integral of plan area with respect to depth:
$$V(z) = \int_0^z A(\xi)\,d\xi = L_bW_bz + m(L_b+W_b)z^2 + \frac{4}{3}m^2 z^3$$The full basin capacity is $V(D_b)$. During routing, basin depth is recovered by numerically inverting $V(z)$ so that the reported depth never exceeds the user-specified maximum depth.
Basin infiltration is driven by the user-entered saturated conductivity proxy (the infiltration rate input, converted from m/day to m/s) with the Soil Moderation Factor applied internally:
$$k_{\text{adj}} = k_{\text{sat}}\,U$$To account for early-time wetting-front behaviour, SoakSIM applies a simplified Green-Ampt-style enhancement:
$$f_{GA}(t) = k_{\text{adj}}\left(1 + \frac{\psi\,\Delta\theta}{F(t)}\right)$$where $\psi$ is the capillary suction head (m), $\Delta\theta$ is the initial moisture deficit, and $F(t)$ is cumulative infiltrated depth (m) per unit active area. This keeps infiltration high at the start of the storm and decreases it as the basin wets up.
If a groundwater depth is supplied, SoakSIM applies an additional damping factor that reduces basin infiltration as the groundwater mound rises toward the basin base:
$$f_{GW}(t) = \max\!\left(0.05,\; \min\!\left(1,\; \frac{D_{gw} - h_{\text{mound}}(t)}{D_{gw}}\right)\right)$$The effective infiltration rate becomes:
$$Q_{\text{infil}}(t) = f_{GA}(t)\,f_{GW}(t)\,A_{\text{active}}(t)$$where $A_{\text{active}}(t)$ is the wetted basin area at the current storage depth, including the sloping side walls. The groundwater mound itself is tracked using a simplified Hantush-style recession term, with the specific yield input controlling how quickly a given infiltrated volume raises the mound.
When the Include clogged layer resistance option is enabled, SoakSIM applies a two-layer series-resistance model using user inputs:
The underlying Green-Ampt conductivity term is $K_u(t)$ and the representative underlying flow-path thickness is $z_u(t)$. The effective conductivity is:
$$K_{\text{eff}}(t) = \frac{z_c + z_u(t)}{\dfrac{z_c}{K_c} + \dfrac{z_u(t)}{K_u(t)}}$$The infiltration capacity used in routing is then computed from $K_{\text{eff}}(t)$ (and groundwater damping where applicable). This makes clogging physically explicit rather than using only a fixed blanket reduction factor.
For basin mode, SoakSIM executes three comparison runs in parallel with the main Green-Ampt/Hantush run. These are adaptations inspired by the methodologies outlined in the PCSump Version 6.1 User Guide and related literature (e.g., Andreyev & Wiseman, 1989), and are provided solely for side-by-side reference.
This adaptation implements a quasi-static representation of the dimensionless recovery method from Andreyev & Wiseman (1989). It tracks horizontal groundwater mounding based on specific yield and vertical Green-Ampt infiltration, dynamically bounding infiltration capacity as the separation to the mound top decreases.
For deep groundwater conditions, this adaptation uses the Green-Ampt infiltration equation applied over a dynamically varying wetted infiltration area, calculating instantaneous capacity based on the actual water level.
This adaptation applies Darcy's equation across a user-defined clogged layer. Infiltration is entirely controlled by the hydraulic conductivity of the clogged layer ($K_c$) and its thickness ($z_c$), with active surface area calculated dynamically as water rises.
In SoakSIM, these three comparison variants are computed silently and presented in a table under the performance graph. The user can toggle the displayed final graph between the models without rerunning the simulation.
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