From 6f1e1fd71964e07f5b5dc576bc99337a002c0fd8 Mon Sep 17 00:00:00 2001 From: ftomei Date: Wed, 14 Aug 2024 13:02:47 +0200 Subject: [PATCH] soil cracks --- agrolib/meteo/meteo.h | 2 +- agrolib/soilFluxes3D/header/soilFluxes3D.h | 1 + agrolib/soilFluxes3D/soilFluxes3D.cpp | 74 +++++++++++---- bin/CRITERIA3D/criteria3DProject.cpp | 104 +++++++++++++++------ bin/CRITERIA3D/shared/project3D.cpp | 4 + 5 files changed, 137 insertions(+), 48 deletions(-) diff --git a/agrolib/meteo/meteo.h b/agrolib/meteo/meteo.h index 312795546..5c6e4a78d 100644 --- a/agrolib/meteo/meteo.h +++ b/agrolib/meteo/meteo.h @@ -108,7 +108,7 @@ enum criteria3DVariable {volumetricWaterContent, waterTotalPotential, waterMatricPotential, availableWaterContent, degreeOfSaturation, soilTemperature, soilSurfaceMoisture, bottomDrainage, waterDeficit, waterInflow, waterOutflow, - factorOfSafety, minimumFactorOfSafety, surfacePond}; + factorOfSafety, minimumFactorOfSafety, surfacePond, maximumVolumetricWaterContent}; const std::map MapDailyMeteoVar = { diff --git a/agrolib/soilFluxes3D/header/soilFluxes3D.h b/agrolib/soilFluxes3D/header/soilFluxes3D.h index d9724bc94..cd3550b56 100644 --- a/agrolib/soilFluxes3D/header/soilFluxes3D.h +++ b/agrolib/soilFluxes3D/header/soilFluxes3D.h @@ -59,6 +59,7 @@ __EXTERN int DLL_EXPORT __STDCALL setWaterSinkSource(long index, double sinkSource); __EXTERN double DLL_EXPORT __STDCALL getWaterContent(long nodeIndex); + __EXTERN double DLL_EXPORT __STDCALL getMaximumWaterContent(long nodeIndex); __EXTERN double DLL_EXPORT __STDCALL getAvailableWaterContent(long nodeIndex); __EXTERN double DLL_EXPORT __STDCALL getWaterDeficit(long index, double fieldCapacity); __EXTERN double DLL_EXPORT __STDCALL getTotalWaterContent(); diff --git a/agrolib/soilFluxes3D/soilFluxes3D.cpp b/agrolib/soilFluxes3D/soilFluxes3D.cpp index a46aa6537..f601d69ec 100644 --- a/agrolib/soilFluxes3D/soilFluxes3D.cpp +++ b/agrolib/soilFluxes3D/soilFluxes3D.cpp @@ -640,23 +640,47 @@ int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve, /*! * \brief getWaterContent * \param nodeIndex - * \return water content at the surface: [m] surface water level; and sub-surface: [m^3 m^-3] volumetric water content + * \return water content at the surface: [m] surface water level; and sub-surface: [m3 m-3] volumetric water content */ double DLL_EXPORT __STDCALL getWaterContent(long nodeIndex) { - if (nodeListPtr == nullptr) return(MEMORY_ERROR); - if ((nodeIndex < 0) || (nodeIndex >= myStructure.nrNodes)) return(INDEX_ERROR); + if (nodeListPtr == nullptr) + return MEMORY_ERROR; + if ((nodeIndex < 0) || (nodeIndex >= myStructure.nrNodes)) + return INDEX_ERROR; if (nodeListPtr[nodeIndex].isSurface) + { /*! surface */ return (nodeListPtr[nodeIndex].H - nodeListPtr[nodeIndex].z); + } else + { /*! sub-surface */ - return (theta_from_Se(nodeIndex)); + return theta_from_Se(nodeIndex); + } } -/*! + /*! + * \brief getMaximumWaterContent + * \param nodeIndex + * \return only sub-surface: [m3 m-3] maximum volumetric water content (theta sat) + */ + double DLL_EXPORT __STDCALL getMaximumWaterContent(long nodeIndex) + { + if (nodeListPtr == nullptr) + return MEMORY_ERROR; + if ((nodeIndex < 0) || (nodeIndex >= myStructure.nrNodes)) + return INDEX_ERROR; + if (nodeListPtr[nodeIndex].isSurface) + return INDEX_ERROR; + + return (nodeListPtr[nodeIndex].Soil->Theta_s); + } + + + /*! * \brief getAvailableWaterContent * \param index * \return available water content (over wilting point) @@ -664,15 +688,21 @@ int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve, */ double DLL_EXPORT __STDCALL getAvailableWaterContent(long index) { - if (nodeListPtr == nullptr) return(MEMORY_ERROR); - if ((index < 0) || (index >= myStructure.nrNodes)) return(INDEX_ERROR); + if (nodeListPtr == nullptr) + return MEMORY_ERROR; + if ((index < 0) || (index >= myStructure.nrNodes)) + return INDEX_ERROR; - if (nodeListPtr[index].isSurface) - /*! surface */ - return (nodeListPtr[index].H - nodeListPtr[index].z); - else - /*! sub-surface */ - return MAXVALUE(0.0, theta_from_Se(index) - theta_from_sign_Psi(-160, index)); + if (nodeListPtr[index].isSurface) + { + // surface + return (nodeListPtr[index].H - nodeListPtr[index].z); + } + else + { + // sub-surface + return MAXVALUE(0.0, theta_from_Se(index) - theta_from_sign_Psi(-160, index)); + } } @@ -680,19 +710,25 @@ int DLL_EXPORT __STDCALL setHydraulicProperties(int waterRetentionCurve, * \brief getWaterDeficit * \param index * \param fieldCapacity - * \return water deficit at surface: 0; sub-surface: [m^3 m^-3] + * \return water deficit at surface: 0; sub-surface: [m3 m-3] */ double DLL_EXPORT __STDCALL getWaterDeficit(long index, double fieldCapacity) { - if (nodeListPtr == nullptr) return(MEMORY_ERROR); - if ((index < 0) || (index >= myStructure.nrNodes)) return(INDEX_ERROR); + if (nodeListPtr == nullptr) + return MEMORY_ERROR; + if ((index < 0) || (index >= myStructure.nrNodes)) + return INDEX_ERROR; if (nodeListPtr[index].isSurface) - /*! surface */ - return (0.0); + { + // surface + return 0; + } else - /*! sub-surface */ + { + // sub-surface return (theta_from_sign_Psi(-fieldCapacity, index) - theta_from_Se(index)); + } } diff --git a/bin/CRITERIA3D/criteria3DProject.cpp b/bin/CRITERIA3D/criteria3DProject.cpp index 529708f75..2e88644f9 100644 --- a/bin/CRITERIA3D/criteria3DProject.cpp +++ b/bin/CRITERIA3D/criteria3DProject.cpp @@ -388,27 +388,28 @@ void Crit3DProject::assignPrecipitation() // Water infiltration into soil cracks -// input: precipitation [mm] +// input: rainfall [mm] // returns the water remaining on the surface after infiltration into soil cracks float Crit3DProject::checkSoilCracking(int row, int col, float precipitation) { const double MAX_CRACKING_DEPTH = 0.6; // [m] - const double MIN_WATER_WATENTIAL = 800; // [kPa] + const double MIN_VOID_VOLUME = 0.2; // [m3 m-3] + const double MAX_VOID_VOLUME = 0.25; // [m3 m-3] - // 1) check soil + // check soil int soilIndex = getSoilIndex(row, col); if (soilIndex == NODATA) return precipitation; - // 2) check pond - long nodeIndex = indexMap.at(0).value[row][col]; - double currentPond = getCriteria3DVar(surfacePond, nodeIndex); - if (precipitation <= currentPond) + // check pond + long surfaceNodeIndex = indexMap.at(0).value[row][col]; + double currentPond = getCriteria3DVar(surfacePond, surfaceNodeIndex); // [mm] + double minimumPond = currentPond * 0.5; // [mm] + if (precipitation <= minimumPond) return precipitation; - // 3) check clay soil + // check clay double maxDepth = std::min(computationSoilDepth, MAX_CRACKING_DEPTH); // [m] - bool isFineFraction = true; int lastFineHorizon = NODATA; int h = 0; @@ -437,38 +438,85 @@ float Crit3DProject::checkSoilCracking(int row, int col, float precipitation) if (maxDepth < 0.2) return precipitation; - // 4) check water potential - double step = 0.05; // [m] - double depth = step; - double waterPotentialSum = 0; + // compute void volume + double stepDepth = 0.05; // [m] + double currentDepth = stepDepth; // [m] + double voidVolumeSum = 0; int nrData = 0; - while (depth <= maxDepth ) + while (currentDepth <= maxDepth ) { - int layerIndex = getSoilLayerIndex(depth); + int layerIndex = getSoilLayerIndex(currentDepth); long nodeIndex = indexMap.at(layerIndex).value[row][col]; - double waterPotential = getCriteria3DVar(waterMatricPotential, nodeIndex); - waterPotentialSum += waterPotential; - depth += step; + + double VWC = getCriteria3DVar(volumetricWaterContent, nodeIndex); // [m3 m-3] + double maxVWC = getCriteria3DVar(maximumVolumetricWaterContent, nodeIndex); // [m3 m-3] + + // TODO: coarse fragment + voidVolumeSum += (maxVWC - VWC); + + currentDepth += stepDepth; nrData++; } - double avgWaterPotential = -soil::metersTokPa(waterPotentialSum / nrData); // [kPa] - if (avgWaterPotential < MIN_WATER_WATENTIAL) + double avgVoidVolume = voidVolumeSum / nrData; // [m3 m-3] + if (avgVoidVolume <= MIN_VOID_VOLUME) return precipitation; // THERE IS A SOIL CRACK - double downWater = precipitation - currentPond; // [mm] - double surfaceWater = precipitation - downWater; // [mm] + double crackRatio = std::min(1.0, (avgVoidVolume - MIN_VOID_VOLUME) / (MAX_VOID_VOLUME - MIN_VOID_VOLUME)); + + double maxInfiltration = precipitation * crackRatio; // [mm] + double surfaceWater = precipitation - maxInfiltration; // [mm] + surfaceWater = std::max(surfaceWater, minimumPond); + double downWater = precipitation - surfaceWater; // [mm] + int lastLayer = getSoilLayerIndex(maxDepth); double area = DEM.header->cellSize * DEM.header->cellSize; // [m2] - int layerIndex = getSoilLayerIndex(maxDepth); - nodeIndex = indexMap.at(layerIndex).value[row][col]; - double flow = area * (downWater / 1000.); // [m3 h-1] - waterSinkSource[nodeIndex] += flow / 3600.; // [m3 s-1] - totalPrecipitation += flow; // [m3] + // 1) infiltration (0.1 mm of water for each soil cm) + for (int l = 1; l <= lastLayer; l++) + { + if (downWater <= 0) + break; + + double layerThick_cm = layerThickness[l] * 100; // [cm] + double layerWater = layerThick_cm * 0.1; // [mm] + layerWater = std::min(layerWater, downWater); + + long nodeIndex = indexMap.at(l).value[row][col]; + double flow = area * (layerWater / 1000.); // [m3 h-1] + waterSinkSource[nodeIndex] += flow / 3600.; // [m3 s-1] + totalPrecipitation += flow; // [m3] + + downWater -= layerWater; + } + + // 2) accumulation on the crack bottom (0.5 mm of water for each soil cm) + for (int l = lastLayer; l > 0; l--) + { + if (downWater <= 0) + break; + + double layerThick_cm = layerThickness[l] * 100; // [cm] + double layerWater = layerThick_cm * 0.5; // [mm] + layerWater = std::min(layerWater, downWater); + + long nodeIndex = indexMap.at(l).value[row][col]; + double flow = area * (layerWater / 1000.); // [m3 h-1] + waterSinkSource[nodeIndex] += flow / 3600.; // [m3 s-1] + totalPrecipitation += flow; // [m3] + + downWater -= layerWater; + } - return surfaceWater; + if (downWater > 0) + { + return surfaceWater + downWater; + } + else + { + return surfaceWater; + } } diff --git a/bin/CRITERIA3D/shared/project3D.cpp b/bin/CRITERIA3D/shared/project3D.cpp index f366e52cf..cd22ff41a 100644 --- a/bin/CRITERIA3D/shared/project3D.cpp +++ b/bin/CRITERIA3D/shared/project3D.cpp @@ -2251,6 +2251,10 @@ double getCriteria3DVar(criteria3DVariable myVar, long nodeIndex) { crit3dVar = soilFluxes3D::getWaterContent(nodeIndex); } + else if (myVar == maximumVolumetricWaterContent) + { + crit3dVar = soilFluxes3D::getMaximumWaterContent(nodeIndex); + } else if (myVar == availableWaterContent) { crit3dVar = soilFluxes3D::getAvailableWaterContent(nodeIndex);