Below is a step-by-step explanation of the infiltration implementation in infil.c
(along with references to infil.h
), from EPA SWMM5. This file handles how SWMM calculates infiltration rates from a subcatchment’s surface into the soil, given different infiltration models (Horton, Modified Horton, Green-Ampt, Modified Green-Ampt, and Curve Number).
1. Overview and Purpose
The infiltration code in SWMM 5 can handle multiple infiltration models:
- Horton (classical decay to a minimum rate),
- Modified Horton (adds a new approach to capacity regeneration),
- Green-Ampt (physically based infiltration theory),
- Modified Green-Ampt (slight variation in infiltration capacity),
- Curve Number (SCS-based approach).
At the start, SWMM determines which infiltration model each subcatchment uses. The infiltration object (one per subcatchment) is stored in a union TInfil
so that each subcatchment's infiltration parameters can be set appropriately.
Key: A subcatchment’s infiltration is computed at each runoff time step, taking into account the chosen infiltration model, the surface water depth, and infiltration states from previous timesteps (like cumulative infiltration).
2. Data Structures
2.1 THorton
typedef struct
{
double f0; // initial infiltration rate (ft/s)
double fmin; // minimum infiltration rate (ft/s)
double Fmax; // maximum total infiltration capacity (ft)
double decay; // decay coefficient of infiltration rate (1/s)
double regen; // regeneration coefficient of infiltration rate (1/s)
//--------------
double tp; // current time on infiltration curve (s)
double Fe; // cumulative infiltration (ft)
} THorton;
What:
f0
: initial infiltration rate, e.g. .fmin
: limiting/minimum infiltration rate.Fmax
: optional overall maximum infiltration capacity (ft).decay
: exponential decay parameter .regen
: exponential “drying” or “regeneration” coefficient.
State:
tp
: how long infiltration has been happening since the subcatchment became wet.Fe
: total infiltration so far in the event (ft).
2.2 TGrnAmpt
typedef struct
{
double S; // capillary suction head (ft)
double Ks; // saturated conductivity (ft/s)
double IMDmax; // max. initial moisture deficit
//--------------
double IMD; // current initial moisture deficit
double F; // total infiltration during storm (ft)
double Fu; // infiltration into upper soil zone (ft)
double Lu; // depth of upper soil zone (ft)
double T; // time until next event (s)
char Sat; // saturation flag
} TGrnAmpt;
What:
S
: capillary suction head .Ks
: saturated hydraulic conductivity.IMDmax
: maximum possible deficit in soil moisture.
State:
IMD
: current deficit.F
: total infiltration volume in current event.Fu
: infiltration into “upper zone.”Lu
: thickness of an effective upper zone of the soil.T
: used in “drying” or “recovery” time for infiltration.Sat
: whether upper zone is saturated.
2.3 TCurveNum
typedef struct
{
double Smax; // maximum infiltration capacity (ft)
double regen; // infiltration regeneration constant (1/s)
double Tmax; // maximum inter-event time (s)
//--------------
double S; // current infiltration capacity (ft)
double F; // current cumulative infiltration (ft)
double P; // current cumulative precipitation (ft)
double T; // current inter-event time (sec)
double Se; // current event infiltration capacity (ft)
double f; // infiltration rate in previous step (ft/s)
} TCurveNum;
What:
Smax
: infiltration capacity derived from the SCS curve number formula.regen
: capacity “regeneration” constant.Tmax
: maximum time after which infiltration is fully regenerated.
State:
S
: current infiltration capacity.F
: cumulative infiltration in current event.P
: cumulative precipitation since event start.T
: time since end of last event.Se
: infiltration capacity for the event.f
: infiltration rate from previous step.
2.4 TInfil
Union
typedef union TInfil {
THorton horton;
TGrnAmpt grnAmpt;
TCurveNum curveNum;
} TInfil;
- Each subcatchment gets a union that can store one of the infiltration model states.
3. Creating and Deleting Infiltration Objects
3.1 infil_create(int n)
void infil_create(int n)
{
Infil = (TInfil *) calloc(n, sizeof(TInfil));
if (Infil == NULL) ErrorCode = ERR_MEMORY;
InfilFactor = 1.0;
return;
}
- Allocates an array of
TInfil
of length n (one for each subcatchment). - Initializes a global factor
InfilFactor = 1.0
(used for monthly/time pattern adjustments to hydraulic conductivity).
3.2 infil_delete()
void infil_delete()
{
FREE(Infil);
}
- Frees the infiltration array.
4. Reading and Initializing Infiltration Parameters
4.1 infil_readParams()
int infil_readParams(int m, char* tok[], int ntoks)
- Reads infiltration parameters from one input line.
m
is the default infiltration model. The code checks the last token to see if it’s an overriding infiltration model.- Depending on the final infiltration method (Horton, Modified Horton, Green-Ampt, Modified Green-Ampt, or Curve Number), a different number of parameters is read.
- Then calls a helper function:
horton_setParams()
grnampt_setParams()
curvenum_setParams()
- On success, assigns the infiltration model to
Subcatch[j].infilModel
.
4.2 Helper Functions
-
horton_setParams(THorton* infil, double p[])
- Reads
f0
,fmin
,decay
,regen
, optionalFmax
. - Converts rates to ft/sec, times to 1/sec, etc.
- Reads
-
grnampt_setParams(TGrnAmpt* infil, double p[])
- Reads
S
,Ks
,IMDmax
. - Also calculates
Lu
from Mein’s equation (empirical).
- Reads
-
curvenum_setParams(TCurveNum* infil, double p[])
- Reads initial curve number, dryness time, etc.
- Converts CN to
Smax
.
4.3 infil_initState(int j)
- Based on the infiltration model, calls:
horton_initState()
,grnampt_initState()
, orcurvenum_initState()
.
- E.g., for Horton, sets
tp = 0.0
,Fe = 0.0
.
5. Managing Infiltration State
5.1 Get/Set State
infil_getState(int j, double x[])
: returns infiltration state variables for subcatchment j.infil_setState(int j, double x[])
: sets infiltration state variables from an arrayx[]
.
Inside, the code delegates to one of:
horton_getState()
/horton_setState()
grnampt_getState()
/grnampt_setState()
curvenum_getState()
/curvenum_setState()
This is used for Hot Start file storage.
6. Infiltration Factor Adjustment
6.1 infil_setInfilFactor(int j)
void infil_setInfilFactor(int j)
{
// Start with global hyd. cond. factor
InfilFactor = Adjust.hydconFactor;
// If subcatch j has a monthly pattern, override with pattern factor
if (j >= 0) {
p = Subcatch[j].infilPattern;
if (p >= 0 && Pattern[p].type == MONTHLY_PATTERN) {
m = datetime_monthOfYear(getDateTime(OldRunoffTime)) - 1;
InfilFactor = Pattern[p].factor[m];
}
}
}
- Applies the global hydraulic conductivity factor from climate adjustments, or possibly a subcatchment-specific infiltration pattern.
- This factor is used in infiltration equations to scale
Ks
,f0
, etc.
7. Computing the Infiltration Rate
7.1 infil_getInfil()
double infil_getInfil(int j, double tstep, double rainfall,
double runon, double depth)
- Master function that calls the correct infiltration model based on
Subcatch[j].infilModel
. - Summation: infiltration rate depends on
rainfall + runon
as the effective “surface inflow” to the infiltration process. - Some models also treat ponded
depth
as an available water volume.
8. Model-Specific Routines
The code then breaks down the infiltration into each model. They share the same function signature:
double X_getInfil(X* infil, double tstep, double irate, double depth);
Here, irate
is net infiltration “supply rate” = rainfall + runon
(ft/s), optionally minus evaporation (for some models). In SWMM’s code, some infiltration models subtract surface evaporation from infiltration supply, but the details can vary.
8.1 Horton & Modified Horton
8.1.1 horton_getInfil()
- The infiltration rate decays exponentially from an initial rate
f0
tofmin
with exponentdecay
. - A
tp
tracks how long infiltration has been going. - If there’s no water available, infiltration capacity “regenerates” at a rate
regen
. - There’s an optional overall limit
Fmax
. - The code does some Newton-Raphson iteration for partial timesteps.
8.1.2 modHorton_getInfil()
- Simplifies the concept of infiltration capacity as
f(t) = f0 - kd * Fe
, not the classic exponential but a linear decrease with cumulative infiltration. - Also includes a regeneration step when no water is available.
8.2 Green-Ampt & Modified Green-Ampt
8.2.1 grnampt_getInfil()
Distinguishes if the upper zone is already saturated (Sat=TRUE
) or unsaturated. Then calls:
grnampt_getSatInfil()
if saturated, orgrnampt_getUnsatInfil()
if unsaturated.
grnampt_getUnsatInfil()
:
- If no water, tries to “recover” the upper zone.
- If water supply < Ks, infiltration = supply.
- If supply > Ks, eventually saturates the upper zone, does partial infiltration unsaturated, partial infiltration saturated.
grnampt_getSatInfil()
:
- Uses an integral solution of the Green-Ampt equation (the “F2” function) with Newton iteration to find infiltration volume.
- The infiltration cannot exceed the actual water supply.
grnampt_getF2()
solves the integrated G-A infiltration equation.
8.3 Curve Number (SCS)
8.3.1 curvenum_getInfil()
- If it starts raining, infiltration capacity is
CN-based
. - The infiltration is
F1 = P * (1 - P/(P+Se))
, a derived expression. - If no rainfall, infiltration capacity regenerates with
regen
. - If there's ponded water but no rainfall, it uses the infiltration rate from the last step.
9. Summary of Steps
- Parse infiltration model type from input (Horton / G-A / CN, etc.).
- Set infiltration object parameters with conversion to consistent units (ft/s, etc.).
- On each timestep:
- Possibly update infiltration factor (monthly, subcatch pattern, etc.).
- Compute infiltration supply =
(rainfall + runon) + ponded water / tstep
. - Call the infiltration model’s function (Horton, G-A, or CN).
- The function returns infiltration rate in
ft/s
. - Update internal infiltration states (cumulative infiltration, dryness, etc.).
- At simulation end or after use, free infiltration objects.
Hence, the infil.c
file ensures that for each subcatchment, infiltration is computed using the selected model with consistent states and time-step increments, so that SWMM can properly route water from surface runoff into infiltration loss.
No comments:
Post a Comment