This code is a function called "biocellFluxRates" that calculates flux rates from the layers of a bio-retention cell LID. It takes in two input arrays, "x" and "f", which represent the vector of storage levels and the vector of flux rates, respectively. The function uses several intermediate variables, including "availVolume" and "maxRate". It also makes use of several properties of the LID layers, such as "soilThickness", "soilPorosity", "soilFieldCap", "soilWiltPoint", "storageThickness", and "storageVoidFrac".
The code first retrieves the moisture levels from the input vector "x" and converts them to volumes. It then calculates the ET rates, soil layer perc rate, and exfiltration rate out of the storage layer. The code also checks for underdrain flow and limits the perc rate and exfiltration rate by available water. The surface infil is limited by unused soil volume. The code also checks for special cases where the storage layer is not present, both layers are full, or either layer is not full and limits the rates accordingly.
A table of the variables used in this function and their descriptions is as follows:
| Variable Name | Description |
|---|---|
| x[SURF], x[SOIL], x[STOR] | Input vector of storage levels |
| f[SURF], f[SOIL], f[STOR] | Output vector of flux rates |
| surfaceDepth | Moisture level variable for the surface layer |
| soilTheta | Moisture level variable for the soil layer |
| storageDepth | Moisture level variable for the storage layer |
| availVolume | Intermediate variable for available volume |
| maxRate | Intermediate variable for maximum rate |
| soilThickness | Property of the soil layer representing its thickness |
| soilPorosity | Property of the soil layer representing its porosity |
| soilFieldCap | Property of the soil layer representing its field capacity |
| soilWiltPoint | Property of the soil layer representing its wilting point |
| storageThickness | Property of the storage layer representing its thickness |
| storageVoidFrac | Property of the storage layer representing its void fraction |
| SurfaceVolume | Converted surface moisture level to volume |
| SoilVolume | Converted soil moisture level to volume |
| StorageVolume | Converted storage moisture level to volume |
| SurfaceInfil | Rate of inflow to the surface layer |
| SurfaceEvap | Rate of evaporation from the surface layer |
| SoilPerc | Rate of percolation from the soil layer |
| SoilEvap | Rate of evaporation from the soil layer |
| StorageExfil | Rate of exfiltration from the storage layer |
| StorageDrain | Rate of underdrain flow from the storage layer |
| Tstep | Time step |
| theLidProc | Pointer to LID process data |
void biocellFluxRates(double x[], double f[])
//
// Purpose: computes flux rates from the layers of a bio-retention cell LID.
// Input: x = vector of storage levels
// Output: f = vector of flux rates
//
{
// Moisture level variables
double surfaceDepth;
double soilTheta;
double storageDepth;
// Intermediate variables
double availVolume;
double maxRate;
// LID layer properties
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];
soilTheta = x[SOIL];
storageDepth = x[STOR];
//... convert moisture levels to volumes
SurfaceVolume = surfaceDepth * theLidProc->surface.voidFrac;
SoilVolume = soilTheta * soilThickness;
StorageVolume = storageDepth * storageVoidFrac;
//... get ET rates
availVolume = SoilVolume - soilWiltPoint * soilThickness;
getEvapRates(SurfaceVolume, 0.0, availVolume, StorageVolume, 1.0);
if ( soilTheta >= soilPorosity ) StorageEvap = 0.0;
//... soil layer perc rate
SoilPerc = getSoilPercRate(soilTheta);
//... limit perc rate by available water
availVolume = (soilTheta - soilFieldCap) * soilThickness;
maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap;
SoilPerc = MIN(SoilPerc, maxRate);
SoilPerc = MAX(SoilPerc, 0.0);
//... exfiltration rate out of storage layer
StorageExfil = getStorageExfilRate();
//... underdrain flow rate
StorageDrain = 0.0;
if ( theLidProc->drain.coeff > 0.0 )
{
StorageDrain = getStorageDrainRate(storageDepth, soilTheta, 0.0,
surfaceDepth);
}
//... special case of no storage layer present
if ( storageThickness == 0.0 )
{
StorageEvap = 0.0;
maxRate = MIN(SoilPerc, StorageExfil);
SoilPerc = maxRate;
StorageExfil = maxRate;
//... limit surface infil. by unused soil volume
maxRate = (soilPorosity - soilTheta) * soilThickness / Tstep +
SoilPerc + SoilEvap;
SurfaceInfil = MIN(SurfaceInfil, maxRate);
}
//... storage & soil layers are full
else if ( soilTheta >= soilPorosity && storageDepth >= storageThickness )
{
//... limiting rate is smaller of soil perc and storage outflow
maxRate = StorageExfil + StorageDrain;
if ( SoilPerc < maxRate )
{
maxRate = SoilPerc;
if ( maxRate > StorageExfil ) StorageDrain = maxRate - StorageExfil;
else
{
StorageExfil = maxRate;
StorageDrain = 0.0;
}
}
else SoilPerc = maxRate;
//... apply limiting rate to surface infil.
SurfaceInfil = MIN(SurfaceInfil, maxRate);
}
//... either layer not full
else if ( storageThickness > 0.0 )
{
//... limit storage exfiltration by available storage volume
maxRate = SoilPerc - StorageEvap + storageDepth*storageVoidFrac/Tstep;
StorageExfil = MIN(StorageExfil, maxRate);
StorageExfil = MAX(StorageExfil, 0.0);
//... 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 perc by unused storage volume
maxRate = StorageExfil + StorageDrain + StorageEvap +
(storageThickness - storageDepth) *
storageVoidFrac/Tstep;
SoilPerc = MIN(SoilPerc, maxRate);
//... limit surface infil. by unused soil volume
maxRate = (soilPorosity - soilTheta) * soilThickness / Tstep +
SoilPerc + SoilEvap;
SurfaceInfil = MIN(SurfaceInfil, maxRate);
}
//... find surface layer outflow rate
SurfaceOutflow = getSurfaceOutflowRate(surfaceDepth);
//... compute overall layer flux rates
f[SURF] = (SurfaceInflow - SurfaceEvap - SurfaceInfil - SurfaceOutflow) /
theLidProc->surface.voidFrac;
f[SOIL] = (SurfaceInfil - SoilEvap - SoilPerc) /
theLidProc->soil.thickness;
if ( storageThickness == 0.0 ) f[STOR] = 0.0;
else f[STOR] = (SoilPerc - StorageEvap - StorageExfil - StorageDrain) /
theLidProc->storage.voidFrac;
}
