Skip to content

Commit

Permalink
update crop transpiration
Browse files Browse the repository at this point in the history
  • Loading branch information
@ftomei
ftomei committed Feb 22, 2024
1 parent 9bf26cd commit a4d0441
Show file tree
Hide file tree
Showing 4 changed files with 70 additions and 90 deletions.
12 changes: 7 additions & 5 deletions agrolib/crop/crop.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<bool> isLayerStressed;
isLayerStressed.resize(nrLayers);
for (unsigned int i = 0; i < nrLayers; i++)
{
isLayerStressed[i] = false;
Expand Down Expand Up @@ -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;
}


Expand Down
35 changes: 17 additions & 18 deletions bin/CRITERIA3D/criteria3DProject.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand All @@ -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]
}
}
}
Expand All @@ -266,7 +266,7 @@ void Crit3DProject::assignETreal()
void Crit3DProject::assignPrecipitation()
{
// initialize
totalPrecipitation = 0;
totalPrecipitation = 0; // [m3]

double area = DEM.header->cellSize * DEM.header->cellSize;

Expand All @@ -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]
}
}
}
}
Expand Down
107 changes: 42 additions & 65 deletions bin/CRITERIA3D/shared/project3D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Expand Down Expand Up @@ -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))
Expand Down Expand Up @@ -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;

Expand All @@ -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;
}
Expand Down
6 changes: 4 additions & 2 deletions bin/CRITERIA3D/shared/project3D.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -120,8 +120,10 @@
std::vector <double> layerThickness; // [m]

// sink/source
std::vector <double> waterSinkSource; // [m^3/sec]
double totalPrecipitation, totalEvaporation, totalTranspiration;
std::vector <double> waterSinkSource; // [m3 s-1]
double totalPrecipitation; // [m3 h-1]
double totalEvaporation; // [m3 h-1]
double totalTranspiration; // [m3 h-1]

Project3D();

Expand Down

0 comments on commit a4d0441

Please sign in to comment.