diff --git a/agrolib/soil/soil.cpp b/agrolib/soil/soil.cpp index 1d1e261ed..19f6abcbb 100644 --- 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 index 236fff3b0..59656ff02 100644 --- 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 index fc09667b7..216cb32f3 100644 --- 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 index 90df788e5..99653c741 100644 --- 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 index 143e6d042..c9a324d86 100644 --- 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