diff --git a/CHANGELOG.md b/CHANGELOG.md index fecb0560..0070f0de 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,11 @@ # Unreleased +**Changed:** + +- Added changes to support groundwater coupling via BMI in + [#221](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/221) + + # [1.5.0](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.5.0) - 3 Jan 2024 @@ -12,7 +18,7 @@ This version of STEMMUS_SCOPE is only compatible with [PyStemmusScope 0.3.0.](ht **Added:** - STEMMUS_SCOPE 'BMI'-like mode ([#208](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/208)): - - The executable can be run in an "interactive" mode. + - The executable can be run in an "interactive" mode. In this mode the model's initialization, time update, and finalization can be called upon separately. - The model now will write away BMI-required variables to a state file, which can be used for the Python BMI. - A dockerfile for the BMI-enabled STEMMUS_SCOPE model, to make the setup easier. diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE index 8d317f62..064d1485 100755 Binary files a/run_model_on_snellius/exe/STEMMUS_SCOPE and b/run_model_on_snellius/exe/STEMMUS_SCOPE differ diff --git a/src/+init/setBoundaryCondition.m b/src/+init/setBoundaryCondition.m index 2519ac34..5e83869d 100644 --- a/src/+init/setBoundaryCondition.m +++ b/src/+init/setBoundaryCondition.m @@ -58,7 +58,6 @@ % 2 --Specified flux(BCh); % 3 --Atmospheric forcing; NBCh = 3; - BCh = -20 / 3600; % this should be fixed in issue 183 if startsWith(initialLandcoverClass, 'Croplands') diff --git a/src/+io/readGroundwaterSettings.m b/src/+io/readGroundwaterSettings.m new file mode 100644 index 00000000..e61cad60 --- /dev/null +++ b/src/+io/readGroundwaterSettings.m @@ -0,0 +1,33 @@ +function GroundwaterSettings = readGroundwaterSettings() + %{ + added by Mostafa, headBotmLayer: head at bottom layer, received from + MODFLOW through BMI indexBotmLayer: index of bottom layer that contains + current headBotmLayer, received from MODFLOW through BMI. The concept of + modflow coupling can be found at: (preprint): + https://doi.org/10.5194/gmd-2022-221 + %} + + % Activate/deactivate Groundwater coupling + GroundwaterSettings.GroundwaterCoupling = 0; % (value = 0 -> deactivate coupling, or = 1 -> activate coupling); default = 0, if update value to = 1 -> through BMI + + % Initialize the head at the bottom layer (start of saturated zone) and the index of the layer that contains that head + GroundwaterSettings.headBotmLayer = 100.0; % head at bottom layer, received from MODFLOW through BMI + GroundwaterSettings.indexBotmLayer = 40; % index of bottom layer that contains current headBotmLayer, received from MODFLOW through BMI + + % Load model settings + ModelSettings = io.getModelSettings(); + NN = ModelSettings.NN; % Number of nodes; + NL = ModelSettings.NL; % Number of layers + + % Calculate soil layer thickness (cumulative layer thickness; e.g. 1, 2, 3, + % 4, 5, 10, 20 ......., last = total_soil_depth) + GroundwaterSettings.soilLayerThickness = zeros(NN, 1); % cumulative soil layer thickness + GroundwaterSettings.soilLayerThickness(1) = 0; + TDeltZ = flip(ModelSettings.DeltZ); + + for ML = 2:NL + GroundwaterSettings.soilLayerThickness(ML) = GroundwaterSettings.soilLayerThickness(ML - 1) + TDeltZ(ML - 1); + end + + GroundwaterSettings.soilLayerThickness(NN) = ModelSettings.Tot_Depth; +end diff --git a/src/+soilmoisture/calculateBoundaryConditions.m b/src/+soilmoisture/calculateBoundaryConditions.m index a8f1513c..1e59cdfd 100644 --- a/src/+soilmoisture/calculateBoundaryConditions.m +++ b/src/+soilmoisture/calculateBoundaryConditions.m @@ -1,5 +1,5 @@ function [AVAIL0, RHS, HeatMatrices, Precip] = calculateBoundaryConditions(BoundaryCondition, HeatMatrices, ForcingData, SoilVariables, InitialValues, ... - TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap) + TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap, GroundwaterSettings) %{ Determine the boundary condition for solving the soil moisture equation. %} @@ -12,21 +12,39 @@ Precip = InitialValues.Precip; Precip_msr = ForcingData.Precip_msr; - Precipp = 0; + % Apply the bottom boundary condition called for by BoundaryCondition.NBChB - if BoundaryCondition.NBChB == 1 % Specify matric head at bottom to be ---BoundaryCondition.BChB; - RHS(1) = BoundaryCondition.BChB; - C4(1, 1) = 1; - RHS(2) = RHS(2) - C4(1, 2) * RHS(1); - C4(1, 2) = 0; - C4_a(1) = 0; - elseif BoundaryCondition.NBChB == 2 % Specify flux at bottom to be ---BoundaryCondition.BChB (Positive upwards); - RHS(1) = RHS(1) + BoundaryCondition.BChB; - elseif BoundaryCondition.NBChB == 3 % BoundaryCondition.NBChB=3,Gravity drainage at bottom--specify flux= hydraulic conductivity; - RHS(1) = RHS(1) - SoilVariables.KL_h(1, 1); + if ~GroundwaterSettings.GroundwaterCoupling % Groundwater Coupling is not activated, added by Mostafa + if BoundaryCondition.NBChB == 1 % Specify matric head at bottom to be ---BoundaryCondition.BChB; + RHS(1) = BoundaryCondition.BChB; + C4(1, 1) = 1; + RHS(2) = RHS(2) - C4(1, 2) * RHS(1); + C4(1, 2) = 0; + C4_a(1) = 0; + elseif BoundaryCondition.NBChB == 2 % Specify flux at bottom to be ---BoundaryCondition.BChB (Positive upwards); + RHS(1) = RHS(1) + BoundaryCondition.BChB; + elseif BoundaryCondition.NBChB == 3 % BoundaryCondition.NBChB=3, Gravity drainage at bottom--specify flux= hydraulic conductivity; + RHS(1) = RHS(1) - SoilVariables.KL_h(1, 1); + end + else % Groundwater Coupling is activated, added by Mostafa + headBotmLayer = GroundwaterSettings.headBotmLayer; % head at bottom layer, received from MODFLOW through BMI + indexBotmLayer = GroundwaterSettings.indexBotmLayer; % index of bottom layer that contains current headBotmLayer, received from MODFLOW through BMI + soilLayerThickness = GroundwaterSettings.soilLayerThickness; + INBT = n - indexBotmLayer + 1; % soil layer thickness from bottom to top (opposite of soilLayerThickness) + BOTm = soilLayerThickness(n); % bottom level of all layers, still need to confirm with Lianyu + if BoundaryCondition.NBChB == 1 % Specify matric head at bottom to be ---BoundaryCondition.BChB; + RHS(INBT) = (headBotmLayer - BOTm + soilLayerThickness(indexBotmLayer)); + C4(INBT, 1) = 1; + RHS(INBT + 1) = RHS(INBT + 1) - C4(INBT, 2) * RHS(INBT); + C4(INBT, 2) = 0; + C4_a(INBT) = 0; + elseif BoundaryCondition.NBChB == 2 % Specify flux at bottom to be ---BoundaryCondition.BChB (Positive upwards); + RHS(INBT) = RHS(INBT) + BoundaryCondition.BChB; + elseif BoundaryCondition.NBChB == 3 % BoundaryCondition.NBChB=3, Gravity drainage at bottom--specify flux= hydraulic conductivity; + RHS(INBT) = RHS(INBT) - SoilVariables.KL_h(INBT, 1); + end end - % Apply the surface boundary condition called for by BoundaryCondition.NBCh if BoundaryCondition.NBCh == 1 % Specified matric head at surface---equal to hN; % h_SUR: Observed matric potential at surface. This variable diff --git a/src/+soilmoisture/solveSoilMoistureBalance.m b/src/+soilmoisture/solveSoilMoistureBalance.m index 9eebfae5..75c9fbb0 100644 --- a/src/+soilmoisture/solveSoilMoistureBalance.m +++ b/src/+soilmoisture/solveSoilMoistureBalance.m @@ -1,5 +1,5 @@ function [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip] = solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, TimeProperties, SoilProperties, ... - BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg) + BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, GroundwaterSettings) %{ Solve the soil moisture balance equation with the Thomas algorithm to update the soil matric potential 'hh', the finite difference @@ -26,7 +26,7 @@ end [AVAIL0, RHS, HeatMatrices, Precip] = soilmoisture.calculateBoundaryConditions(BoundaryCondition, HeatMatrices, ForcingData, SoilVariables, InitialValues, ... - TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap); + TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap, GroundwaterSettings); [CHK, hh, C4] = soilmoisture.solveTridiagonalMatrixEquations(HeatMatrices.C4, SoilVariables.hh, HeatMatrices.C4_a, RHS); % update structures diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 06680fad..b19f9a5b 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -245,6 +245,9 @@ % get soil constants SoilConstants = io.getSoilConstants(); + %% Groundwater coupling settings (added by Mostafa) + GroundwaterSettings = io.readGroundwaterSettings(); + %% The boundary condition information settings BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, SiteProperties.landcoverClass(1)); DSTOR = BoundaryCondition.DSTOR; @@ -297,6 +300,7 @@ % The state can be modified by the STEMMUS_SCOPE BMI. See PyStemmusScope. if strcmp(bmiMode, 'update') load([OutputPath, 'STEMMUS_SCOPE_state.mat']); % Load the workspace to be able to (continue) running the model + if KT + 1 >= TimeProperties.Dur_tot bmiMode = 'none'; % Ensure the model does not try to update. disp("Finished running the model. Updating won't do anything!"); @@ -308,6 +312,10 @@ disp('The calculations start now'); end +if GroundwaterSettings.GroundwaterCoupling == 1 % Groundwater coupling is enabled + BoundaryCondition.NBChB = 1; +end + % Actually run the model if strcmp(bmiMode, 'update') || strcmp(runMode, 'full') % Will do one timestep in "update mode", and run until the end if in "full run" mode. @@ -591,7 +599,7 @@ GasDispersivity = conductivity.calculateGasDispersivity(InitialValues, SoilVariables, P_gg, k_g); % Srt is both input and output - [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip] = soilmoisture.solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, TimeProperties, SoilProperties, BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg); + [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip] = soilmoisture.solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, TimeProperties, SoilProperties, BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, GroundwaterSettings); if BoundaryCondition.NBCh == 1 DSTOR = 0; diff --git a/src/STEMMUS_SCOPE_exe.m b/src/STEMMUS_SCOPE_exe.m index dae007ef..2150de1b 100644 --- a/src/STEMMUS_SCOPE_exe.m +++ b/src/STEMMUS_SCOPE_exe.m @@ -24,7 +24,9 @@ function STEMMUS_SCOPE_exe(config_file, runMode) 'KT', ... % Index of current time step 'SiteProperties', ... % Site properties (e.g. lat, lon) 'fluxes', ... % Atmospheric fluxes - 'TT' ... % Soil temperature over depth + 'TT', ... % Soil temperature over depth + 'SoilVariables', ... % Structure that includes different variables of soil moisture, added by Mostafa + 'GroundwaterSettings' ... % added by Mostafa }; %#ok % Variables for tracking the state of the model initialization: