fix soil fraction

ARPA-SIMC · Dec 19, 2024 · 298ef31 · 298ef31
1 parent 9586d44
commit 298ef31
Show file tree

Hide file tree

Showing 5 changed files with 41 additions and 38 deletions.
diff --git a/agrolib/soil/soil.cpp b/agrolib/soil/soil.cpp
@@ -152,6 +152,7 @@ namespace soil
 
         this->fieldCapacity = NODATA;
         this->wiltingPoint = NODATA;
+        this->waterContentSAT = NODATA;
         this->waterContentFC = NODATA;
         this->waterContentWP = NODATA;
 
@@ -223,19 +224,17 @@ namespace soil
 
         horizonPtr = horizonPointer;
 
-        double hygroscopicHumidity = -2000;     // [kPa]
-        double waterContentHH = soil::thetaFromSignPsi(hygroscopicHumidity, *horizonPtr);
-
-        // [-]
-        soilFraction = horizonPtr->getSoilFraction();
-
         // [mm]
-        SAT = horizonPtr->vanGenuchten.thetaS * soilFraction * thickness * 1000;
-        FC = horizonPtr->waterContentFC * soilFraction * thickness * 1000;
-        WP = horizonPtr->waterContentWP * soilFraction * thickness * 1000;
-        HH = waterContentHH * soilFraction * thickness * 1000;
+        SAT = horizonPtr->waterContentSAT * thickness * 1000;
+        FC = horizonPtr->waterContentFC * thickness * 1000;
+        WP = horizonPtr->waterContentWP * thickness * 1000;
         critical = FC;
 
+        // hygroscopic humidity
+        double hygroscopicHumidity = -2000;     // [kPa]
+        double volWaterContentHH = soil::thetaFromSignPsi(hygroscopicHumidity, *horizonPtr);    // [m3 m-3]
+        HH = volWaterContentHH * horizonPtr->getSoilFraction() * thickness * 1000;
+
         return true;
     }
 
@@ -413,8 +412,10 @@ namespace soil
     // estimate bulk density from total porosity
     double estimateBulkDensity(const Crit3DHorizon &horizon, double totalPorosity, bool increaseWithDepth)
     {
-        if (int(totalPorosity) == int(NODATA))
+        if (isEqual(totalPorosity, NODATA))
+        {
             totalPorosity = (horizon.vanGenuchten.refThetaS);
+        }
 
         double specificDensity = estimateSpecificDensity(horizon.organicMatter);
         double refBulkDensity = (1 - totalPorosity) * specificDensity;
@@ -630,7 +631,7 @@ namespace soil
      * \brief Compute volumetric water content from signed water potential
      * \param signPsi water potential       [kPa]
      * \param horizon
-     * \return volumetric water content     [m^3 m-3]
+     * \return volumetric water content     [m3 m-3]
      */
     double thetaFromSignPsi(double signPsi, const Crit3DHorizon &horizon)
     {     
@@ -710,14 +711,14 @@ namespace soil
 
 
     /*!
-     * \brief return current volumetric water content [m3 m^3]
+     * \brief return current volumetric water content [-]
      */
     double Crit1DLayer::getVolumetricWaterContent()
     {
-        // waterContent [mm]
-        // thickness [m]
-        double theta = waterContent / (thickness * soilFraction * 1000);
-        return theta;
+        // thickness [m] -> mm
+        double soilThickness = thickness * soilFraction * 1000.;
+
+        return waterContent / soilThickness;
     }
 
 
@@ -727,7 +728,7 @@ namespace soil
     double Crit1DLayer::getDegreeOfSaturation()
     {
         double theta = getVolumetricWaterContent();
-        return (theta - horizonPtr->vanGenuchten.thetaR) / (horizonPtr->vanGenuchten.thetaS - horizonPtr->vanGenuchten.thetaR);
+        return SeFromTheta(theta, *horizonPtr);
     }
 
 
@@ -963,8 +964,10 @@ namespace soil
 
         horizon.fieldCapacity = soil::getFieldCapacity(horizon.texture.clay, soil::KPA);
         horizon.wiltingPoint = soil::getWiltingPoint(soil::KPA);
-        horizon.waterContentFC = soil::thetaFromSignPsi(horizon.fieldCapacity, horizon);
-        horizon.waterContentWP = soil::thetaFromSignPsi(horizon.wiltingPoint, horizon);
+
+        horizon.waterContentSAT = horizon.vanGenuchten.thetaS * horizon.getSoilFraction();
+        horizon.waterContentFC = soil::thetaFromSignPsi(horizon.fieldCapacity, horizon) * horizon.getSoilFraction();
+        horizon.waterContentWP = soil::thetaFromSignPsi(horizon.wiltingPoint, horizon) * horizon.getSoilFraction();
 
         return true;
     }
@@ -999,8 +1002,9 @@ namespace soil
             psiMin = std::min(psiMin, horizon.dbData.waterRetention[i].water_potential);
             thetaMax = std::max(thetaMax, horizon.dbData.waterRetention[i].water_content);
         }
-        // add theta sat if minimum observed value is greater than 3 kPa
-        bool addThetaSat = ((thetaMax < horizon.vanGenuchten.thetaS) && (psiMin > 3));
+
+        // add theta sat if minimum observed value is greater than 5 kPa
+        bool addThetaSat = ((thetaMax < horizon.vanGenuchten.thetaS) && (psiMin > 5));
 
         // set values
         unsigned int nrValues = nrObsValues;

diff --git a/agrolib/soil/soil.h b/agrolib/soil/soil.h
@@ -140,8 +140,11 @@
 
             double fieldCapacity;               /*!<  [kPa]     */
             double wiltingPoint;                /*!<  [kPa]     */
+
+            double waterContentSAT;             /*!<  [m3 m-3]*/
             double waterContentFC;              /*!<  [m3 m-3]*/
             double waterContentWP;              /*!<  [m3 m-3]*/
+
             double PH;                          /*!<  [-]       */
             double CEC;                         /*!<  [meq/100g]*/
 

diff --git a/agrolib/soilWidget/soilWidget.cpp b/agrolib/soilWidget/soilWidget.cpp
@@ -754,13 +754,13 @@ void Crit3DSoilWidget::setInfoTextural(int nHorizon)
     }
     else
     {
-        if (mySoil.horizon[unsigned(nHorizon)].vanGenuchten.thetaS == NODATA)
+        if (mySoil.horizon[unsigned(nHorizon)].waterContentSAT == NODATA)
         {
             satValue->setText(QString::number(NODATA));
         }
         else
         {
-            satValue->setText(QString::number(mySoil.horizon[unsigned(nHorizon)].vanGenuchten.thetaS, 'f', 3));
+            satValue->setText(QString::number(mySoil.horizon[unsigned(nHorizon)].waterContentSAT, 'f', 3));
         }
 
         if (mySoil.horizon[unsigned(nHorizon)].waterContentFC == NODATA)

diff --git a/bin/CRITERIA3D/criteria3DProject.cpp b/bin/CRITERIA3D/criteria3DProject.cpp
@@ -925,15 +925,12 @@ double Crit3DProject::getSoilVar(int soilIndex, int layerIndex, soil::soilVariab
         return soilList[unsigned(soilIndex)].horizon[hIndex].wiltingPoint;
     else if (myVar == soil::soilWaterPotentialFC)
         return soilList[unsigned(soilIndex)].horizon[hIndex].fieldCapacity;
+    else if (myVar == soil::soilWaterContentSat)
+        return soilList[unsigned(soilIndex)].horizon[hIndex].waterContentSAT;
     else if (myVar == soil::soilWaterContentFC)
         return soilList[unsigned(soilIndex)].horizon[hIndex].waterContentFC;
-    else if (myVar == soil::soilWaterContentSat)
-        return soilList[unsigned(soilIndex)].horizon[hIndex].vanGenuchten.thetaS;
     else if (myVar == soil::soilWaterContentWP)
-    {
-        double signPsiLeaf = -160;      //[m]
-        return soil::thetaFromSignPsi(signPsiLeaf, soilList[unsigned(soilIndex)].horizon[hIndex]);
-    }
+        return soilList[unsigned(soilIndex)].horizon[hIndex].waterContentWP;
     else
         return NODATA;
 }

diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp
@@ -1691,14 +1691,14 @@ bool Project3D::computeSurfaceWaterContent(double &wcSum, long &nrVoxels)
     {
         for (int col = 0; col < indexMap.at(0).header->nrCols; col++)
         {
-            long nodeIndex = indexMap.at(0).value[row][col];
-            if (nodeIndex != indexMap.at(0).header->flag)
+            long surfaceNodeIndex = indexMap.at(0).value[row][col];
+            if (surfaceNodeIndex != indexMap.at(0).header->flag)
             {
-                double wc = getCriteria3DVar(volumetricWaterContent, nodeIndex);    //[m]
+                double surfaceWC = getCriteria3DVar(volumetricWaterContent, surfaceNodeIndex);    // [m]
 
-                if (wc != NODATA)
+                if (! isEqual(surfaceWC, NODATA))
                 {
-                    wcSum += wc * voxelArea;        // [m3]
+                    wcSum += surfaceWC * voxelArea;        // [m3]
                     nrVoxels++;
                 }
             }
@@ -2117,8 +2117,7 @@ double Project3D::assignTranspiration(int row, int col, double currentLai, doubl
 
         // [m3 m-3]
         double volWaterContent = getCriteria3DVar(volumetricWaterContent, nodeIndex);
-        double thetaSat = horizon.vanGenuchten.thetaS;
-        double volWaterSurplusThreshold = thetaSat - waterSurplusStressFraction * (thetaSat - horizon.waterContentFC);
+        double volWaterSurplusThreshold = horizon.waterContentSAT - waterSurplusStressFraction * (horizon.waterContentSAT - horizon.waterContentFC);
         double volWaterScarcityThreshold = horizon.waterContentFC - currentCrop.fRAW * (horizon.waterContentFC - horizon.waterContentWP);
 
         double ratio;
@@ -2137,7 +2136,7 @@ double Project3D::assignTranspiration(int row, int col, double currentLai, doubl
         else if  ((volWaterContent - volWaterSurplusThreshold)  > EPSILON)
         {
             // WATER SURPLUS
-            ratio = (thetaSat - volWaterContent) / (thetaSat - volWaterSurplusThreshold);
+            ratio = (horizon.waterContentSAT - volWaterContent) / (horizon.waterContentSAT - volWaterSurplusThreshold);
             isLayerStressed[layer] = true;
         }
         else