This code is computing flux rates for the layers of a porous pavement LID (low impact development). It takes in a vector of storage levels (x) and returns a vector of flux rates (f). The code uses various intermediate variables, LID layer properties, and moisture levels to calculate evaporation rates, infiltration rates, percolation rates, and exfiltration rates. It also checks for saturated layers and limits rates accordingly.
| Variable | Purpose |
|---|---|
| surfaceDepth | moisture level of surface layer |
| paveDepth | moisture level of pavement layer |
| soilTheta | moisture level of soil layer |
| storageDepth | moisture level of storage layer |
| pervFrac | fraction of LID that is pervious |
| storageInflow | inflow rate to storage layer |
| availVolume | available volume in soil layer |
| maxRate | maximum rate for a specific layer |
| paveVoidFrac | void fraction of pavement layer |
| paveThickness | thickness of pavement layer |
| soilThickness | thickness of soil layer |
| soilPorosity | porosity of soil layer |
| soilFieldCap | field capacity of soil layer |
| soilWiltPoint | wilt point of soil layer |
| storageThickness | thickness of storage layer |
| storageVoidFrac | void fraction of storage layer |
| SurfaceVolume | volume of surface layer |
| PaveVolume | volume of pavement layer |
| SoilVolume | volume of soil layer |
| StorageVolume | volume of storage layer |
| SurfaceInfil | infiltration rate of surface layer |
| PavePerc | percolation rate out of pavement layer |
| SoilPerc | percolation rate out of soil layer |
| StorageExfil | exfiltration rate out of storage layer |
| StorageDrain | underdrain flow rate |
Additionally, the code also calculates the following rates:
- Evaporation rates: getEvapRates() is called to calculate evaporation rates for the surface, pavement, soil, and storage layers.
- Pavement percolation rate: getPavementPermRate() is called to calculate the nominal rate of surface infiltration into the pavement layer.
- Soil percolation rate: getSoilPercRate() is called to calculate the percolation rate out of the soil layer.
- Storage exfiltration rate: getStorageExfilRate() is called to calculate the exfiltration rate out of the storage layer.
- Storage underdrain flow rate: getStorageDrainRate() is called to calculate the underdrain flow rate, if the coefficient for the drain is greater than 0.
Finally, the code checks for saturated layers and limits the rates accordingly. For example, if the pavement or soil layer is saturated, the storage evaporation rate is set to 0. Additionally, the infiltration rate into the pavement layer is limited by the available water in the pavement layer and the percolation rate out of the pavement layer is limited by the available water in the soil layer
void pavementFluxRates(double x[], double f[])
//
// Purpose: computes flux rates for the layers of a porous pavement LID.
// Input: x = vector of storage levels
// Output: f = vector of flux rates
//
{
//... Moisture level variables
double surfaceDepth;
double paveDepth;
double soilTheta;
double storageDepth;
//... Intermediate variables
double pervFrac = (1.0 - theLidProc->pavement.impervFrac);
double storageInflow; // inflow rate to storage layer (ft/s)
double availVolume;
double maxRate;
//... LID layer properties
double paveVoidFrac = theLidProc->pavement.voidFrac * pervFrac;
double paveThickness = theLidProc->pavement.thickness;
double soilThickness = theLidProc->soil.thickness;
double soilPorosity = theLidProc->soil.porosity;
double soilFieldCap = theLidProc->soil.fieldCap;
double soilWiltPoint = theLidProc->soil.wiltPoint;
double storageThickness = theLidProc->storage.thickness;
double storageVoidFrac = theLidProc->storage.voidFrac;
//... retrieve moisture levels from input vector
surfaceDepth = x[SURF];
paveDepth = x[PAVE];
soilTheta = x[SOIL];
storageDepth = x[STOR];
//... convert moisture levels to volumes
SurfaceVolume = surfaceDepth * theLidProc->surface.voidFrac;
PaveVolume = paveDepth * paveVoidFrac;
SoilVolume = soilTheta * soilThickness;
StorageVolume = storageDepth * storageVoidFrac;
//... get ET rates
availVolume = SoilVolume - soilWiltPoint * soilThickness;
getEvapRates(SurfaceVolume, PaveVolume, availVolume, StorageVolume,
pervFrac);
//... no storage evap if soil or pavement layer saturated
if ( paveDepth >= paveThickness ||
( soilThickness > 0.0 && soilTheta >= soilPorosity )
) StorageEvap = 0.0;
//... find nominal rate of surface infiltration into pavement layer
SurfaceInfil = SurfaceInflow + (SurfaceVolume / Tstep);
//... find perc rate out of pavement layer
PavePerc = getPavementPermRate() * pervFrac;
//... surface infiltration can't exceed pavement permeability
SurfaceInfil = MIN(SurfaceInfil, PavePerc);
//... limit pavement perc by available water
maxRate = PaveVolume/Tstep + SurfaceInfil - PaveEvap;
maxRate = MAX(maxRate, 0.0);
PavePerc = MIN(PavePerc, maxRate);
//... find soil layer perc rate
if ( soilThickness > 0.0 )
{
SoilPerc = getSoilPercRate(soilTheta);
availVolume = (soilTheta - soilFieldCap) * soilThickness;
maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap;
SoilPerc = MIN(SoilPerc, maxRate);
SoilPerc = MAX(SoilPerc, 0.0);
}
else SoilPerc = PavePerc;
//... exfiltration rate out of storage layer
StorageExfil = getStorageExfilRate();
//... underdrain flow rate
StorageDrain = 0.0;
if ( theLidProc->drain.coeff > 0.0 )
{
StorageDrain = getStorageDrainRate(storageDepth, soilTheta, paveDepth,
surfaceDepth);
}
//... check for adjacent saturated layers
//... no soil layer, pavement & storage layers are full
if ( soilThickness == 0.0 &&
storageDepth >= storageThickness &&
paveDepth >= paveThickness )
{
//... pavement outflow can't exceed storage outflow
maxRate = StorageEvap + StorageDrain + StorageExfil;
if ( PavePerc > maxRate ) PavePerc = maxRate;
//... storage outflow can't exceed pavement outflow
else
{
//... use up available exfiltration capacity first
StorageExfil = MIN(StorageExfil, PavePerc);
StorageDrain = PavePerc - StorageExfil;
}
//... set soil perc to pavement perc
SoilPerc = PavePerc;
//... limit surface infil. by pavement perc
SurfaceInfil = MIN(SurfaceInfil, PavePerc);
}
//... pavement, soil & storage layers are full
else if ( soilThickness > 0 &&
storageDepth >= storageThickness &&
soilTheta >= soilPorosity &&
paveDepth >= paveThickness )
{
//... find which layer has limiting flux rate
maxRate = StorageExfil + StorageDrain;
if ( SoilPerc < maxRate) maxRate = SoilPerc;
else maxRate = MIN(maxRate, PavePerc);
//... use up available storage exfiltration capacity first
if ( maxRate > StorageExfil ) StorageDrain = maxRate - StorageExfil;
else
{
StorageExfil = maxRate;
StorageDrain = 0.0;
}
SoilPerc = maxRate;
PavePerc = maxRate;
//... limit surface infil. by pavement perc
SurfaceInfil = MIN(SurfaceInfil, PavePerc);
}
//... storage & soil layers are full
else if ( soilThickness > 0.0 &&
storageDepth >= storageThickness &&
soilTheta >= soilPorosity )
{
//... soil perc can't exceed storage outflow
maxRate = StorageDrain + StorageExfil;
if ( SoilPerc > maxRate ) SoilPerc = maxRate;
//... storage outflow can't exceed soil perc
else
{
//... use up available exfiltration capacity first
StorageExfil = MIN(StorageExfil, SoilPerc);
StorageDrain = SoilPerc - StorageExfil;
}
//... limit surface infil. by available pavement volume
availVolume = (paveThickness - paveDepth) * paveVoidFrac;
maxRate = availVolume / Tstep + PavePerc + PaveEvap;
SurfaceInfil = MIN(SurfaceInfil, maxRate);
}
//... soil and pavement layers are full
else if ( soilThickness > 0.0 &&
paveDepth >= paveThickness &&
soilTheta >= soilPorosity )
{
PavePerc = MIN(PavePerc, SoilPerc);
SoilPerc = PavePerc;
SurfaceInfil = MIN(SurfaceInfil,PavePerc);
}
//... no adjoining layers are full
else
{
//... limit storage exfiltration by available storage volume
// (if no soil layer, SoilPerc is same as PavePerc)
maxRate = SoilPerc - StorageEvap + StorageVolume / Tstep;
maxRate = MAX(0.0, maxRate);
StorageExfil = MIN(StorageExfil, maxRate);
//... limit underdrain flow by volume above drain offset
if ( StorageDrain > 0.0 )
{
maxRate = -StorageExfil - StorageEvap;
if (storageDepth >= storageThickness ) maxRate += SoilPerc;
if ( theLidProc->drain.offset <= storageDepth )
{
maxRate += (storageDepth - theLidProc->drain.offset) *
storageVoidFrac/Tstep;
}
maxRate = MAX(maxRate, 0.0);
StorageDrain = MIN(StorageDrain, maxRate);
}
//... limit soil & pavement outflow by unused storage volume
availVolume = (storageThickness - storageDepth) * storageVoidFrac;
maxRate = availVolume/Tstep + StorageEvap + StorageDrain + StorageExfil;
maxRate = MAX(maxRate, 0.0);
if ( soilThickness > 0.0 )
{
SoilPerc = MIN(SoilPerc, maxRate);
maxRate = (soilPorosity - soilTheta) * soilThickness / Tstep +
SoilPerc;
}
PavePerc = MIN(PavePerc, maxRate);
//... limit surface infil. by available pavement volume
availVolume = (paveThickness - paveDepth) * paveVoidFrac;
maxRate = availVolume / Tstep + PavePerc + PaveEvap;
SurfaceInfil = MIN(SurfaceInfil, maxRate);
}
//... surface outflow
SurfaceOutflow = getSurfaceOutflowRate(surfaceDepth);
//... compute overall layer flux rates
f[SURF] = SurfaceInflow - SurfaceEvap - SurfaceInfil - SurfaceOutflow;
f[PAVE] = (SurfaceInfil - PaveEvap - PavePerc) / paveVoidFrac;
if ( theLidProc->soil.thickness > 0.0)
{
f[SOIL] = (PavePerc - SoilEvap - SoilPerc) / soilThickness;
storageInflow = SoilPerc;
}
else
{
f[SOIL] = 0.0;
storageInflow = PavePerc;
SoilPerc = 0.0;
}
f[STOR] = (storageInflow - StorageEvap - StorageExfil - StorageDrain) /
storageVoidFrac;
}
