Attempt v.0.3 to modify STEMMUS_SCOPE for MODFLOW coupling

Hi @SarahAlidoost, @BSchilperoort, @yijianzeng In trying to adapt the code to calculate the groundwater recharge, I have uploaded a new branch and I have some questions

I have made two new functions (findPhreaticSurface and calculateGroundwaterRecharge). a- The findPhreaticSurface function finds the soil layer that includes the phreatic surface (saturated zone). The input to that function is a new variable called Sh. Sh is the soil matric potential (but from bottom to top, opposite of the variable h, following as what was done in STEMMUS_MODFLOW (see here)

b- The calculateGroundwaterRecharge function calculates the groundwater recharge. Here, I have some questions as follows:

1- I need to call the INBT variable, which is calculated already in the calculateBoundaryConditions function (line 34). In the main time, I had calculated the INBT again (lines 42-43), but I wonder if there is a way to call that INBT variable only without all inputs and outputs of the function?

2- A similar issue as in point 1, I need the QMT variable from the calculatesSoilWaterFluxes function. When I use the lines 65-66, I received an error that SAVE and hh variables are not definable. It is strange to me because in the function itself, SAVE and hh variables are not definable. How to call the QMT variable then?

3- Calculations of a new variable called RTOP are at lines 52-62, which is the sum of QL, and QV. Those variables are calculated at the energy/calcualteEnergyParameters function (lines 62-73). However, the calculations in the energy/calcualteEnergyParameters function are over the whole soil profile (see here, while here in the RTOP calculations, I need the variables at the saturated zone only, which is from the INBT to NL (so not from 1 to NL). What is the best way to do so? Do I need to repeat all the calculations in the calculateGroundwaterRecharge function again?

4- continue point 3, also something I noticed, QL is calculated for every soil layer (line 63), while Qv and Qa are only single values, how to sum them then?

5- My question in my last comment on the issue 211 is relevant. What does QMT mean?

6- I need to import new variables (SY and SS) from MODFLOW (similar to line 137 and lines 147-148 in STEMMUS_MODFLOW). Will we do this in BMI or can i make an input for this two variables from txt files?

7- I add the two functions in the soilmoisture folder. Is that the best location or maybe can I create a new folder called groundwater and maybe i can also move the readGroundwaterSettings to that new folder?

I then need to call those two functions in the main script. As I understand from Lianyu, those two functions need to run only with the MODFLOW timestep. However, as we agreed we will not include the MODFLOW timestep in the main code and we will do this in the BMI, then I am not sure where in the code I should call these two functions exactly. I followed what Lianyu did in STEMMUS_MODFLOW and based on what I understood so far as follows:

a- The initialization phase (lines 244-252 and lines 259-260). That was the same as lines 227-234 of STEMMUS_MODFLOW code. b- The update phase (lines 696-706). That was the same as lines 389-399 of STEMMUS_MODFLOW code.

I need confirmation if I insert the initialization and the update phases in the correct places in the code or not.
I also need a bit of explanation on how the update of the variables works (for example h, hh and Theta_L, Theta_LL)

Hi, I made a couple of changes and made a new branch (SSM_v.0.3.1), which are needed for the coupling as below:

A new folder called groundwater with 3 files The readGroundwaterSettings was moved from the io folder to the groundwater folder Two new files findPhreaticSurface and calculateGroundwaterRecharge are added

Changes in the readGroundwaterSettings:

A couple of variables will be called from MODFLOW (GroundwaterSettings.headBotmLayer, GroundwaterSettings.tempBotm, GroundwaterSettings.numAqL, GroundwaterSettings.aqLayers, GroundwaterSettings.SY, GroundwaterSettings.SS). The meaning of all of them are described in lines 11-30.
The GroundwaterSettings.indxBotmLayer will not be called from MODFLOW, but rather calculated (lines 80-99) inside the STEMMUS_SCOPE code (that was based on my better understanding of the Lianyu's STEMMUS_MODFLOW code).

The new function findPhreaticSurface: This function finds the soil layer where the groundwater recharge need to be exported and passed to MODFLOW

The new function calculateGroundwaterRecharge : This function calls the fluxes of interest to calculate the groundwater recharge. All details are described inside the function.

Modifications in other functions:

The main change after activating the groundwater coupling is that the soil profile is not fixed (default = 5 m) anymore and is changed per timestep based on the index of the bottom soil layer indxBotmLayer calculated in the readGroundwaterSettings (lines 80-99). The soil profile is discretized vertically in number of layers (start from index = 1 to index = NN, where NN is the number of soil nodes). Due to the coupling, the soil profile will not be from 1 to NN, but will be from indexBotmLayer to NN. That will affect several functions (almost all functions) in the soilmositure, energy and dryair folders.
An example is in lines 29-33 in the assembleCoefficientMatrices).
Also, when calling the RHS in any function, I change the RHS(1) to RHS(indexbotm). An example is in line 34 in the calcualteTimeDerivatives.
The soilmoisture/calculateBoundaryCondtions has some changes (here), and also the energy/calculateBoundaryCondtions has some changes (here)
An issue about the C6_a variable that was submitted here
I still have to confirm some issues with Yijian about the new changes, so that will be done soon.
When pushing up the code to the Github, the spaces in the lines are messed up. I think that will be cleaned after i use the miss_hit, so you may ignore these mistakes in the meantime
Two main Python scripts will be needed (one for interpolating the groundwater head received from MODFLOW, and one for aggregating the groundwater recharge to be passed to MODFLOW)

EcoExtreML / STEMMUS_SCOPE

Attempt v.0.3 to modify STEMMUS_SCOPE for MODFLOW coupling #226