Merge commit '0af1e82186a498b4095b6b35f1273a4fdd058c73'

ARPA-SIMC · Oct 27, 2023 · 1284215 · 1284215
2 parents a5d88ae + 0af1e82
commit 1284215
Show file tree

Hide file tree

Showing 5 changed files with 29 additions and 13 deletions.
diff --git a/agrolib/criteriaModel/criteria1DCase.cpp b/agrolib/criteriaModel/criteria1DCase.cpp
@@ -34,6 +34,7 @@
 #include "criteria1DCase.h"
 #include "soilFluxes3D.h"
 #include "soil.h"
+#include "crop.h"
 
 
 Crit1DOutput::Crit1DOutput()
@@ -561,11 +562,25 @@ bool Crit1DCase::computeDailyModel(Crit3DDate &myDate, std::string &error)
         {
             double rootDensityNorm;
             if (rootDensityMax == 0)
+            {
                 rootDensityNorm = 0;
+            }
             else
+            {
                 rootDensityNorm = crop.roots.rootDensity[l] / rootDensityMax;
+            }
+
+            double rootCohesion;
+            if (isEqual(crop.roots.rootsAdditionalCohesion, NODATA))
+            {
+                rootCohesion = 0;
+            }
+            else
+            {
+                rootCohesion = crop.roots.rootsAdditionalCohesion * rootDensityNorm;
+            }
 
-            soilLayers[l].factorOfSafety = soilLayers[l].computeSlopeStability(unit.slope, rootDensityNorm);
+            soilLayers[l].factorOfSafety = soilLayers[l].computeSlopeStability(unit.slope, rootCohesion);
         }
     }
 

diff --git a/agrolib/crop/cropDbTools.cpp b/agrolib/crop/cropDbTools.cpp
@@ -67,6 +67,12 @@ bool loadCropParameters(const QSqlDatabase &dbCrop, QString idCrop, Crit3DCrop &
     myCrop.roots.rootDepthMax = query.value("root_depth_max").toDouble();
     myCrop.roots.actualRootDepthMax = myCrop.roots.rootDepthMax;
 
+    if (! getValue(query.value("roots_additional_cohesion"), &(myCrop.roots.rootsAdditionalCohesion)))
+    {
+        // default: no mechanical effect of roots
+        myCrop.roots.rootsAdditionalCohesion = 0;
+    }
+
     getValue(query.value("degree_days_root_increase"), &(myCrop.roots.degreeDaysRootGrowth));
     if (myCrop.roots.degreeDaysRootGrowth == NODATA)
     {

diff --git a/agrolib/crop/root.cpp b/agrolib/crop/root.cpp
@@ -60,6 +60,7 @@ void Crit3DRoot::clear()
     rootLength = NODATA;
     rootDepth = NODATA;
     rootDensity.clear();
+    rootsAdditionalCohesion = NODATA;
 }
 
 
@@ -219,8 +220,6 @@ namespace root
         if (currentDD <= 0) return 0.;
         if (currentDD > myRoot->degreeDaysRootGrowth) return maxRootLength;
 
-        double halfDevelopmentPoint = myRoot->degreeDaysRootGrowth * 0.5 ;
-
         if (myRoot->growth == LINEAR)
         {
             myRootLength = maxRootLength * (currentDD / myRoot->degreeDaysRootGrowth);
@@ -230,20 +229,14 @@ namespace root
             double logMax, logMin,deformationFactor;
             double iniLog = log(9.);
             double filLog = log(1 / 0.99 - 1);
-            double k,b;
-            k = -(iniLog - filLog) / (emergenceDD - myRoot->degreeDaysRootGrowth);
-            b = -(filLog + k * myRoot->degreeDaysRootGrowth);
+            double k = -(iniLog - filLog) / (emergenceDD - myRoot->degreeDaysRootGrowth);
+            double b = -(filLog + k * myRoot->degreeDaysRootGrowth);
 
             logMax = (myRoot->actualRootDepthMax) / (1 + exp(-b - k * myRoot->degreeDaysRootGrowth));
             logMin = myRoot->actualRootDepthMax / (1 + exp(-b));
             deformationFactor = (logMax - logMin) / maxRootLength ;
             myRootLength = 1.0 / deformationFactor * (myRoot->actualRootDepthMax / (1.0 + exp(-b - k * currentDD)) - logMin);
         }
-        else if (myRoot->growth == EXPONENTIAL)
-        {
-            // not used in Criteria Bdp
-            myRootLength = maxRootLength * (1.- exp(-2.*(currentDD/halfDevelopmentPoint)));
-        }
 
         return myRootLength;
     }
@@ -571,5 +564,6 @@ namespace root
 
         return true;
     }
+
 }
 
diff --git a/agrolib/crop/root.h b/agrolib/crop/root.h
@@ -33,6 +33,7 @@
         int firstRootLayer;                 /*!< [-]  */
         int lastRootLayer;                  /*!< [-]  */
         std::vector<double> rootDensity;    /*!< [-]  */
+        double rootsAdditionalCohesion;      /*!< [kPa] Cr = roots reinforcement (RR) derived from a model */
 
         /*! state variables */
         double rootDepth;                   /*!<  [m]  current root depth */

diff --git a/agrolib/soil/soil.cpp b/agrolib/soil/soil.cpp
@@ -747,7 +747,7 @@ namespace soil
      * \return factor of safety FoS [-]
      * if fos < 1 the slope is unstable
      */
-    double Crit3DLayer::computeSlopeStability(double slope, double normalizedRootDensity)
+    double Crit3DLayer::computeSlopeStability(double slope, double rootCohesion)
     {
         double suctionStress = -waterPotential * getDegreeOfSaturation();    // [kPa]
 
@@ -760,7 +760,7 @@ namespace soil
         double frictionEffect =  tanFrictionAngle / tanAngle;
 
         double unitWeight = horizonPtr->bulkDensity * GRAVITY;                   // [kN m-3]
-        double cohesionEffect = 2 * horizonPtr->effectiveCohesion / (unitWeight * depth * sin(2*slopeAngle));
+        double cohesionEffect = 2 * (horizonPtr->effectiveCohesion + rootCohesion) / (unitWeight * depth * sin(2*slopeAngle));
 
         double suctionEffect = (suctionStress * (tanAngle + 1/tanAngle) * tanFrictionAngle) / (unitWeight * depth);