SWMM5 gwater.c Summary

 Below is a step-by-step explanation of gwater.c, which manages groundwater (GW) flow in EPA SWMM5. SWMM can model saturated and unsaturated subsurface zones, compute exchange with surface infiltration, evapotranspiration, and flow to or from the sewer system at a node. This file implements all GW computations for each subcatchment’s groundwater component.

---

1. Purpose

//   gwater.c
//   Groundwater functions.
//   ...

• gwater.c is responsible for reading aquifer/GW parameters, initializing states, updating the moisture content and depth in an upper/lower zone model each time step, and computing lateral outflow to a connected node as well as deep percolation from the lower zone.

---

2. Key Data Structures

• TAquifer: defines aquifer properties (porosity, field capacity, conductivity, etc.).
• TGroundwater: defines the subcatchment’s local groundwater parameters (node ID, surface elevation, the user’s custom parameters a1,b1,a2,b2,a3, etc.).

The array Aquifer[] holds all defined aquifers, and each subcatchment can reference one via Subcatch[j].groundwater->aquifer.

---

3. Reading Aquifer and Subcatchment GW Parameters

3.1 gwater_readAquiferParams()

int gwater_readAquiferParams(int j, char* tok[], int ntoks)

• Parses one line of data for aquifer $j$.
• Data format (13 tokens, the 14th is optional):
  1. Aquifer ID,
  2. Porosity,
  3. Wilting point,
  4. Field capacity,
  5. Hydraulic conductivity,
  6. Conductivity slope,
  7. Tension slope,
  8. Upper evaporation fraction,
  9. Lower evap. depth,
  10. Groundwater recession constant,
  11. Aquifer bottom elevation,
  12. Water table elevation,
  13. Initial upper moisture content,
  14. (optional) monthly pattern for upper evaporation.
• Assigns them to Aquifer[j] with correct unit conversions (like dividing length by UCF(LENGTH)).

3.2 gwater_readGroundwaterParams()

int gwater_readGroundwaterParams(char* tok[], int ntoks)

• Reads the subcatchment's GW parameters for the line:

  subcatch  aquifer  node  surfElev  a1  b1  a2  b2  a3  fixedDepth
             nodeElev  bottomElev  waterTableElev  upperMoisture
  

• This sets up TGroundwater within Subcatch[j]:
  • aquifer: index of the aquifer used,
  • node: index of node that subcatch’s groundwater can flow into,
  • surfElev: subcatch's ground surface elevation,
  • a1, b1, a2, b2, a3, fixedDepth: user’s curve parameters for computing GW infiltration,
  • optionally: nodeElev, bottomElev, waterTableElev, upperMoisture.

---

4. Custom Flow Expressions

4.1 gwater_readFlowExpression(...)

• SWMM allows user‐supplied math expressions for lateral or deep GW flow:

  subcatch LATERAL <expr>
  subcatch DEEP    <expr>
  

• This function:
  1. Finds the subcatch index,
  2. Checks if the expression is for LATERAL or DEEP,
  3. Tokenizes and parses the expression (via mathexpr_create) into a syntax tree stored in Subcatch[j].gwLatFlowExpr or Subcatch[j].gwDeepFlowExpr.

4.2 gwater_deleteFlowExpression(...)

• Releases memory for the expression trees associated with one subcatch.

---

5. Validation and Initialization

5.1 gwater_validateAquifer(int j)

• Checks correctness of aquifer parameters (ensuring porosity > fieldCapacity > wiltingPoint, etc.).
• Prints error if invalid.

5.2 gwater_validate(int j)

• For subcatch j, if it has groundwater, ensures:
  • surfElev >= waterTableElev (otherwise error).

5.3 gwater_initState(int j)

• Sets the subcatch’s initial groundwater state:
  • theta = upperMoisture (the initial soil moisture).
  • lowerDepth = waterTableElev - bottomElev (the saturated thickness).
  • If that is higher than surfElev - bottomElev, reduce it to avoid overflow.
  • maxInfilVol is how much infiltration can fill the upper zone.

---

6. State Saving/Reading

6.1 gwater_getState(int j, double x[])

• Copies the subcatch’s theta, lowerDepth (absolute water table height), newFlow, and maxInfilVol into array x.

6.2 gwater_setState(int j, double x[])

• Assigns these state variables back from x[]. Used for hotstart or re‐initializing from saved states.

---

7. Primary GW Computation

7.1 gwater_getGroundwater(int j, double evap, double infil, double tStep)

1. Finds subcatch’s TGroundwater* GW and the TAquifer object.
2. Calculate infiltration rate $infil / area / tStep$.
3. Convert any surface evap. to a rate that might reduce available GW evaporation.
4. Determine total aquifer thickness, the node invert, water surface in node, etc.
5. Sets up initial conditions for an ODE system:
   • upper zone moisture theta
   • lower zone depth lowerDepth.
6. Defines constraints for maximum percolation, max outflow, etc.
7. Calls an ODE solver (odesolve_integrate(...)) with the function getDxDt(...), which calculates the time derivative:
   • $\frac{d(\theta)}{dt}$ and $\frac{d(\text{lowerDepth})}{dt}$.
8. After integration, calls getFluxes(...) to compute the actual flux rates:
   • UpperEvap, LowerEvap, UpperPerc, GWFlow, LowerLoss.
   • updateMassBal(...) updates the water balance.
9. Saves new states back to GW->theta and GW->lowerDepth.

---

7.2 getDxDt(t, x, dxdt)

• The ODE solver calls this function at each step.
• x[THETA]: the upper moisture content, x[LOWERDEPTH]: depth of the lower saturated zone.
• Calls getFluxes(...), which sets global flux variables (Infil, UpperEvap, etc.).
• Then: $$\frac{d(\theta)}{dt} = \frac{\text{(Infil - UpperEvap - UpperPerc)}}{\text{(upperDepth)}},$$ $$\frac{d(\text{lowerDepth})}{dt} = \frac{\text{(UpperPerc - LowerLoss - LowerEvap - GWFlow)}}{\text{(porosity - }\theta)},$$ if denominators are nonzero.

---

7.3 getFluxes(double theta, double lowerDepth)

1. upperDepth = total depth minus lowerDepth.
2. Calls getEvapRates(...) to compute UpperEvap & LowerEvap.
3. Calls getUpperPerc(...) to compute percolation from upper to lower zone.
4. If a user expression for DeepFlowExpr is defined, that sets LowerLoss. Otherwise use Aquifer[j].lowerLossCoeff * lowerDepth / totalDepth.
5. If a user expression for LatFlowExpr is defined, we add that to the base getGWFlow(...).
6. GWFlow is also limited to max possible flow based on the unsaturated or saturated zone capacity.

---

7.4 getEvapRates(double theta, double upperDepth)

• Disallows infiltration and evaporation at the same time: if infiltration is > 0, then no GW evaporation.
• Otherwise, partition potential evap (MaxEvap) between the upper zone (some fraction upperEvapFrac) and the lower zone if the lower zone is within the “lowerEvapDepth.”

---

7.5 getUpperPerc(double theta, double upperDepth)

• If the upper zone is deeper than 0 and theta > fieldCapacity, compute an exponent on unsaturated conductivity as: $$\text{hydcon} = \text{conductivity} \times e^{(theta - porosity)\times conductSlope}.$$
• Then modifies for tension slope. Returns this percolation rate in ft/sec.

---

7.6 getGWFlow(double lowerDepth)

• The “traditional” flow formula: $$GWFlow = \frac{(a1 * (lowerDepth - Hstar)^{b1} \;-\; a2 * (Hsw - Hstar)^{b2} \;+\; a3 * lowerDepth * Hsw) } { \text{UCF(GWFLOW)} }$$ when $(lowerDepth - Hstar) > 0$.
• If a user lateral flow expression is present, add that as well.

---

7.7 updateMassBal(area, tStep)

• Converts flux rates to volumes over subcatchment area times time step, calls massbal_updateGwaterTotals(...).

---

8. Helper Routines for Expressions

8.1 getVariableIndex(char* s)

• Maps a variable name ("HGW", "THETA", etc.) to an integer ID used in the math expression parser.

8.2 getVariableValue(int varIndex)

• Returns the current numerical value of that variable from the many global placeholders (like Hgw, Theta, etc.).

---

9. Summary

gwater.c implements SWMM’s two-zone groundwater model:

• Upper unsaturated zone with moisture content theta.
• Lower saturated zone with depth lowerDepth.

At each time step:

1. Convert infiltration from surface, potential evap, node water surface, etc.
2. Use an ODE solver to integrate the 2 state variables.
3. Compute fluxes: infiltration, evap from each zone, percolation, deep loss, and lateral GW flow to the node.
4. Update the mass balance.

This model can be augmented by custom expressions for lateral or deep flow, allowing flexible definitions for real-world groundwater behavior.