From a4d04415b62659e550be029e1451efcfb8dc4e3e Mon Sep 17 00:00:00 2001 From: ftomei Date: Thu, 22 Feb 2024 14:55:58 +0100 Subject: [PATCH] update crop transpiration --- agrolib/crop/crop.cpp | 12 +-- bin/CRITERIA3D/criteria3DProject.cpp | 35 +++++---- bin/CRITERIA3D/shared/project3D.cpp | 107 +++++++++++---------------- bin/CRITERIA3D/shared/project3D.h | 6 +- 4 files changed, 70 insertions(+), 90 deletions(-) diff --git a/agrolib/crop/crop.cpp b/agrolib/crop/crop.cpp index 48789f8ba..79b0b735e 100644 --- a/agrolib/crop/crop.cpp +++ b/agrolib/crop/crop.cpp @@ -724,7 +724,10 @@ double Crit3DCrop::computeTranspiration(double maxTranspiration, const std::vect // initialize unsigned int nrLayers = unsigned(soilLayers.size()); - bool* isLayerStressed = new bool[nrLayers]; + + // initialize vectors + std::vector isLayerStressed; + isLayerStressed.resize(nrLayers); for (unsigned int i = 0; i < nrLayers; i++) { isLayerStressed[i] = false; @@ -821,14 +824,13 @@ double Crit3DCrop::computeTranspiration(double maxTranspiration, const std::vect waterStress = 1 - (TRs / maxTranspiration); - double totalTranspiration = 0; + double actualTranspiration = 0; for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++) { - totalTranspiration += layerTranspiration[unsigned(i)]; + actualTranspiration += layerTranspiration[unsigned(i)]; } - delete[] isLayerStressed; - return totalTranspiration; + return actualTranspiration; } diff --git a/bin/CRITERIA3D/criteria3DProject.cpp b/bin/CRITERIA3D/criteria3DProject.cpp index 05eeaa234..605db661e 100644 --- a/bin/CRITERIA3D/criteria3DProject.cpp +++ b/bin/CRITERIA3D/criteria3DProject.cpp @@ -222,8 +222,8 @@ void Crit3DProject::dailyUpdateCrop() */ void Crit3DProject::assignETreal() { - totalEvaporation = 0; - totalTranspiration = 0; + totalEvaporation = 0; // [m3 h-1] + totalTranspiration = 0; // [m3 h-1] double area = DEM.header->cellSize * DEM.header->cellSize; @@ -244,18 +244,18 @@ void Crit3DProject::assignETreal() lai = 0; } - // assign real evaporation - double realEvap = assignEvaporation(row, col, lai, soilIndex); // [mm] - double evapFlow = area * (realEvap / 1000.); // [m3 h-1] - totalEvaporation += evapFlow; + // assigns actual evaporation + double actualEvap = assignEvaporation(row, col, lai, soilIndex); // [mm h-1] + double evapFlow = area * (actualEvap / 1000.); // [m3 h-1] + totalEvaporation += evapFlow; // [m3 h-1] - // assign real transpiration + // assigns actual transpiration if (lai > 0) { float degreeDays = degreeDaysMap.value[row][col]; - double realTransp = assignTranspiration(row, col, lai, degreeDays); // [mm] - double traspFlow = area * (realTransp / 1000.); // [m3 h-1] - totalTranspiration += traspFlow; + double actualTransp = assignTranspiration(row, col, lai, degreeDays); // [mm h-1] + double traspFlow = area * (actualTransp / 1000.); // [m3 h-1] + totalTranspiration += traspFlow; // [m3 h-1] } } } @@ -266,7 +266,7 @@ void Crit3DProject::assignETreal() void Crit3DProject::assignPrecipitation() { // initialize - totalPrecipitation = 0; + totalPrecipitation = 0; // [m3] double area = DEM.header->cellSize * DEM.header->cellSize; @@ -278,16 +278,15 @@ void Crit3DProject::assignPrecipitation() int surfaceIndex = indexMap.at(0).value[row][col]; if (surfaceIndex != indexMap.at(0).header->flag) { - double waterSource = 0; float prec = hourlyMeteoMaps->mapHourlyPrec->value[row][col]; if (! isEqual(prec, hourlyMeteoMaps->mapHourlyPrec->header->flag)) - waterSource += prec; - - if (waterSource > 0) { - double flow = area * (waterSource / 1000.); // [m3 h-1] - waterSinkSource[unsigned(surfaceIndex)] += flow / 3600.; // [m3 s-1] - totalPrecipitation += flow; + if (prec > 0) + { + double flow = area * (prec / 1000.); // [m3 h-1] + waterSinkSource[surfaceIndex] += flow / 3600.; // [m3 s-1] + totalPrecipitation += flow; // [m3] + } } } } diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp index 35b81208c..6ff9b4d73 100644 --- a/bin/CRITERIA3D/shared/project3D.cpp +++ b/bin/CRITERIA3D/shared/project3D.cpp @@ -766,10 +766,10 @@ void Project3D::computeWaterBalance3D(double totalTimeStep) { double previousWaterContent = soilFluxes3D::getTotalWaterContent(); - logInfo("total water [m^3]: " + QString::number(previousWaterContent)); - logInfo("precipitation [m^3]: " + QString::number(totalPrecipitation)); - logInfo("evaporation [m^3]: " + QString::number(-totalEvaporation)); - logInfo("transpiration [m^3]: " + QString::number(-totalTranspiration)); + logInfo("total water [m3]: " + QString::number(previousWaterContent)); + logInfo("precipitation [m3]: " + QString::number(totalPrecipitation)); + logInfo("evaporation [m3]: " + QString::number(-totalEvaporation)); + logInfo("transpiration [m3]: " + QString::number(-totalTranspiration)); logInfo("Compute water flow..."); soilFluxes3D::initializeBalance(); @@ -1385,15 +1385,15 @@ double Project3D::assignEvaporation(int row, int col, double lai, int soilIndex) /*! \brief assignTranspiration - * Assigns the actual crop transpiration (waterSinkSource) from the soil root zone - * \param row, col + * it computes the actual crop transpiration from the soil root zone + * and assigns to waterSinkSource * \param currentLai: leaf area index [m2 m-2] * \param currentDegreeDays: degree days sum [°C] - * \return crop transpiration sum [mm] + * \return actual transpiration [mm] */ double Project3D::assignTranspiration(int row, int col, double currentLai, double currentDegreeDays) { - double actualTranspiration = 0; + double actualTranspiration = 0; // [mm] // check lai and degree days if (currentLai < EPSILON || isEqual(currentDegreeDays, NODATA)) @@ -1458,14 +1458,14 @@ double Project3D::assignTranspiration(int row, int col, double currentLai, doubl layerTranspiration[i] = 0; } - // set water surplus stress threshold - // 0 = saturation 1 = field capacity - double waterSurplusStressFraction = 0.5; // [-] + // set water surplus stress threshold (0: saturation 1: field capacity) + double waterSurplusStressFraction = 0.5; // [-] if (currentCrop.isWaterSurplusResistant()) { waterSurplusStressFraction = 0.0; } + double rootDensityWithoutStress = 0.0; // [-] int firstRootLayer = currentCrop.roots.firstRootLayer; int lastRootLayer = currentCrop.roots.lastRootLayer; @@ -1481,84 +1481,61 @@ double Project3D::assignTranspiration(int row, int col, double currentLai, doubl double volWaterSurplusThreshold = thetaSat - waterSurplusStressFraction * (thetaSat - horizon.waterContentFC); double volWaterScarcityThreshold = horizon.waterContentFC - currentCrop.fRAW * (horizon.waterContentFC - horizon.waterContentWP); - if ((volWaterContent - volWaterSurplusThreshold) > EPSILON) + double ratio; + if (volWaterContent <= horizon.waterContentWP) { - // WATER SURPLUS - double ratio = (thetaSat - volWaterContent) / (thetaSat - volWaterSurplusThreshold); - layerTranspiration[layer] = maxTranspiration * currentCrop.roots.rootDensity[layer] * ratio; + // NO AVAILABLE WATER + ratio = 0; isLayerStressed[layer] = true; } - else if (volWaterContent <= horizon.waterContentWP) + else if (volWaterContent < volWaterScarcityThreshold) { - layerTranspiration[layer] = 0; + // WATER SCARSITY + ratio = (volWaterContent - horizon.waterContentWP) / (volWaterScarcityThreshold - horizon.waterContentWP); isLayerStressed[layer] = true; } - else if (volWaterContent < volWaterScarcityThreshold) + else if ((volWaterContent - volWaterSurplusThreshold) > EPSILON) { - // WATER SCARSITY - double ratio = (volWaterContent - horizon.waterContentWP) / (volWaterScarcityThreshold - horizon.waterContentWP); - layerTranspiration[layer] = maxTranspiration * currentCrop.roots.rootDensity[layer] * ratio; + // WATER SURPLUS + ratio = (thetaSat - volWaterContent) / (thetaSat - volWaterSurplusThreshold); isLayerStressed[layer] = true; } else { // NORMAL CONDITION + ratio = 1; + isLayerStressed[layer] = false; + rootDensityWithoutStress += currentCrop.roots.rootDensity[layer]; } - } - - /* - else - { - // normal conditions - layerTranspiration[i] = maxTranspiration * roots.rootDensity[i]; - - TRs += layerTranspiration[i]; - TRe += layerTranspiration[i]; - if ((soilLayers[i].waterContent - layerTranspiration[i]) > waterScarcityThreshold) - { - isLayerStressed[i] = false; - totRootDensityWithoutStress += roots.rootDensity[i]; - } - else - { - isLayerStressed[i] = true; - } - } + layerTranspiration[layer] = maxTranspiration * currentCrop.roots.rootDensity[layer] * ratio; + actualTranspiration += layerTranspiration[layer]; } - double totRootDensityWithoutStress = 0.0; // [-] - double redistribution = 0.0; // [mm] - // WATER STRESS [-] - double firstWaterStress = 1 - (TRs / maxTranspiration); + double waterStress = 1 - (actualTranspiration / maxTranspiration); - // Hydraulic redistribution - // the movement of water from moist to dry soil through plant roots - // TODO add numerical process - if (firstWaterStress > EPSILON && totRootDensityWithoutStress > EPSILON) + // Hydraulic redistribution: movement of water from moist to dry soil through plant roots + if (waterStress > EPSILON && rootDensityWithoutStress > EPSILON) { - // redistribution acts on not stressed roots - redistribution = MINVALUE(firstWaterStress, totRootDensityWithoutStress) * maxTranspiration; + // redistribution acts only on not stressed roots + double redistribution = maxTranspiration * std::min(waterStress, rootDensityWithoutStress); // [mm] - for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++) - { - if (! isLayerStressed[i]) - { - double addLayerTransp = redistribution * (roots.rootDensity[unsigned(i)] / totRootDensityWithoutStress); - layerTranspiration[unsigned(i)] += addLayerTransp; - TRs += addLayerTransp; - } - } + for (int layer = firstRootLayer; layer <= lastRootLayer; layer++) + if (! isLayerStressed[layer]) + layerTranspiration[layer] += redistribution * (currentCrop.roots.rootDensity[layer] / rootDensityWithoutStress); } - waterStress = 1 - (TRs / maxTranspiration); - - for (int i = roots.firstRootLayer; i <= roots.lastRootLayer; i++) + // assigns transpiration to water sink source + double area = DEM.header->cellSize * DEM.header->cellSize; // [m2] + actualTranspiration = 0; + for (int layer = firstRootLayer; layer <= lastRootLayer; layer++) { - totalTranspiration += layerTranspiration[unsigned(i)]; + double flow = area * (layerTranspiration[layer] / 1000); // [m3 h-1] + long nodeIndex = long(indexMap.at(layer).value[row][col]); + waterSinkSource.at(nodeIndex) -= (flow / 3600); // [m3 s-1] + actualTranspiration += layerTranspiration[layer]; // [mm] } - */ return actualTranspiration; } diff --git a/bin/CRITERIA3D/shared/project3D.h b/bin/CRITERIA3D/shared/project3D.h index 464e344bd..05ff4b8e2 100644 --- a/bin/CRITERIA3D/shared/project3D.h +++ b/bin/CRITERIA3D/shared/project3D.h @@ -120,8 +120,10 @@ std::vector layerThickness; // [m] // sink/source - std::vector waterSinkSource; // [m^3/sec] - double totalPrecipitation, totalEvaporation, totalTranspiration; + std::vector waterSinkSource; // [m3 s-1] + double totalPrecipitation; // [m3 h-1] + double totalEvaporation; // [m3 h-1] + double totalTranspiration; // [m3 h-1] Project3D();