Saturday, January 21, 2023

SWMM 5.2.2 Code for LID Function pavementFluxRates

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.

VariablePurpose
surfaceDepthmoisture level of surface layer
paveDepthmoisture level of pavement layer
soilThetamoisture level of soil layer
storageDepthmoisture level of storage layer
pervFracfraction of LID that is pervious
storageInflowinflow rate to storage layer
availVolumeavailable volume in soil layer
maxRatemaximum rate for a specific layer
paveVoidFracvoid fraction of pavement layer
paveThicknessthickness of pavement layer
soilThicknessthickness of soil layer
soilPorosityporosity of soil layer
soilFieldCapfield capacity of soil layer
soilWiltPointwilt point of soil layer
storageThicknessthickness of storage layer
storageVoidFracvoid fraction of storage layer
SurfaceVolumevolume of surface layer
PaveVolumevolume of pavement layer
SoilVolumevolume of soil layer
StorageVolumevolume of storage layer
SurfaceInfilinfiltration rate of surface layer
PavePercpercolation rate out of pavement layer
SoilPercpercolation rate out of soil layer
StorageExfilexfiltration rate out of storage layer
StorageDrainunderdrain 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;
}

No comments:

GitHub code and Markdown (MD) files Leveraging

 To better achieve your goal of leveraging your GitHub code and Markdown (MD) files for your WordPress blog or LinkedIn articles, consider t...