Saturday, December 28, 2024

SWMM5 infil.c Summary

 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. f0f_0.
  • fmin: limiting/minimum infiltration rate.
  • Fmax: optional overall maximum infiltration capacity (ft).
  • decay: exponential decay parameter α\alpha.
  • 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 (ψ)(\psi).
  • 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

  1. horton_setParams(THorton* infil, double p[])

    • Reads f0, fmin, decay, regen, optional Fmax.
    • Converts rates to ft/sec, times to 1/sec, etc.
  2. grnampt_setParams(TGrnAmpt* infil, double p[])

    • Reads S, Ks, IMDmax.
    • Also calculates Lu from Mein’s equation (empirical).
  3. 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(), or curvenum_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 array x[].

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()

  1. The infiltration rate decays exponentially from an initial rate f0 to fmin with exponent decay.
  2. A tp tracks how long infiltration has been going.
  3. If there’s no water available, infiltration capacity “regenerates” at a rate regen.
  4. There’s an optional overall limit Fmax.
  5. The code does some Newton-Raphson iteration for partial timesteps.

8.1.2 modHorton_getInfil()

  1. Simplifies the concept of infiltration capacity as f(t) = f0 - kd * Fe, not the classic exponential but a linear decrease with cumulative infiltration.
  2. 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, or
  • grnampt_getUnsatInfil() if unsaturated.

grnampt_getUnsatInfil():

  1. If no water, tries to “recover” the upper zone.
  2. If water supply < Ks, infiltration = supply.
  3. If supply > Ks, eventually saturates the upper zone, does partial infiltration unsaturated, partial infiltration saturated.

grnampt_getSatInfil():

  1. Uses an integral solution of the Green-Ampt equation (the “F2” function) with Newton iteration to find infiltration volume.
  2. 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()

  1. If it starts raining, infiltration capacity is CN-based.
  2. The infiltration is F1 = P * (1 - P/(P+Se)), a derived expression.
  3. If no rainfall, infiltration capacity regenerates with regen.
  4. If there's ponded water but no rainfall, it uses the infiltration rate from the last step.

9. Summary of Steps

  1. Parse infiltration model type from input (Horton / G-A / CN, etc.).
  2. Set infiltration object parameters with conversion to consistent units (ft/s, etc.).
  3. 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.).
  4. 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:

A comprehensive explanation of how minimum travel distance relates to link length in InfoSewer

In hydraulic modeling of sewer networks, the minimum travel distance is a fundamental parameter that affects how accurately the model can si...