SSM v.0.2: Code changes to activate Groundwater coupling

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)
+
+
 <a name="1.5.0"></a>
 # [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.

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')

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

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

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

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;

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: