Attempt v.0.1 to modify STEMMUS_SCOPE for MODLFOW coupling

[MostafaGomaa93]

The Lianyu's STEMMUS_MODFLOW code includes 4-5 modules that need to be included in STEMMUS_SCOPE for the coupling with MODFLOW. I am trying to include those modules one by one. The first one is called "TIME_INTERPOLATION.m". I have created a new branch called "SSM_v.0.1". The changes I did so far in this branch is the following:

1- Add the TIME_INTERPOLATION.m 2- Changes in the "STEMMUS_SCOPE.m" file to read the variables got from the new module "TIME_INTERPOLATION.m" as follows: 2.1- One constant at define_constant.m -> line 44 2.2- Additional model settings at getModelSettings -> lines [82-93] 2.3- additional settings at lines [49-61] 2.4- update variables at lines 343-350. The variables here should be received from MODFLOW. In the meantime, I have just assigned them as arbitrary values.

I would need, please a- @Lianyu-Yu and @yijianzeng: to check if the "TIME_INTERPOLATION.m" module is doing its job with the arbitrarily assigned values so far or we can change the values to test. I am still not sure about the variables meaning, so can't judge alone so far. In the meantime, no changes are expected in the results because the outputs from the function are not passed anywhere yet, but will do in the next step b- Then, @SarahAlidoost and @BSchilperoort to check the variable names and the variable locations

Please correct me if I miss anything or if anything is not clear

[MostafaGomaa93]

I copy paste the meaning of the variables I got from Lianyu from issue 212, so to help to track if those variables are already in the code and if their names are changed after refactorating

TEND -> end of the current time step; TIME -> start of the current time step; IP0STm -> indicate the number of STEMMUS columns for MODFLOW; BOTm -> elevation of the bottom layer of MODFLOW; HPILR1N -> water heads at the end of the current time step; HPILR0 -> water heads at the start of the current time step; DELT -> time interval of the current time step; XElemnt -> layer thickness of STEMMUS soil layers; NN -> number of nodes (number of soil layers in STEMMUS + 1); Q3DF -> indicator for the quasi-3d simulation =1, yes; =0, no; ADAPTF -> indicator for the adaptive lower boundary setting, =1, the moving lower boundary scheme is used; =0, the fixed lower boundary scheme is used.

My expectations are: 1- TEND, TIME, DELT, IP0STm, BOTm., XElemnt, NN are there in the code, but could be with different names 2- HPILR1N, HPILR0 are new variables

[MostafaGomaa93]

Update Testing the code with the new changes gave an error in the calculations (the arbitrary values here didn't work and have to be changed). I would appreciate the help of @Lianyu-Yu for this

[MostafaGomaa93]

Update 2 @Lianyu-Yu sent me by email the reason why the code was not working and I fixed the error. Now I have updated the description of the changes (my first comment above). @SarahAlidoost, could you please check the changes and let me know if I need to make a pull request? or I continue with a new issue and new branch since no changes to the outputs so far because the new variables are not used in the code yet

[SarahAlidoost]

@SarahAlidoost, could you please check the changes and let me know if I need to make a pull request? or I continue with a new issue and new branch since no changes to the outputs so far because the new variables are not used in the code yet

It seems that these changes be applied to the version of STEMMUS_SCOPE model that lives in the main branch. Please first update your local main branch by doing git pull, then make a new branch, then implement those changes and then submit a pull request. Here some comments/suggestions:

1- Add the TIME_INTERPOLATION.m

The script should be refactored to improve readability and code quality. For example, documentation and perhaps a reference should be added, proper meaningful variable names should be chosen, ...

2- Changes in the "STEMMUS_SCOPE.m" file to read the variables got from the new module "TIME_INTERPOLATION.m" as follows: 2.1- Some constants at lines 64-86

Are these really constants or settings. Constant should be added to define_constants.m whereas settings to getModelSettings.m

2.2- update variables at lines 353-360. The variables here should be received from MODFLOW. In the meantime, I have just assigned them as arbitrary values. Some variables like (TEND, TIME, NN, XElement) should be already there in the STEMMUS_SCOPE code so far, but their names may be changed after refactoring. So, this has to be checked also

adding variables should be fine as STMMUS_SCOPE does not use those variables. However, they should be written to some files. Also TEND, TIME, NN are already there except XElement. Please add inline comments to explain why you make those changes.

When you submit a pull request, look at the checklist in the pull request. Every item should be checked accordingly. Hope it helps.

[MostafaGomaa93]

@Lianyu-Yu and @yijianzeng Could you please check the table below and correct and/or fill in the description of the variables

Variable	Variable Description	Location
IPOSTm	number of STEMMUS columns for MODFLOW	STEMMUS_SCOPE.m & getModelSettings.m
BTOm	elevation of the bottom layer of MODFLOW	STEMMUS_SCOPE.m & getModelSettings.m
HPILR1N	water heads (m) at the end of the current time step, should be recieved from MODFLOW	STEMMUS_SCOPE.m & TIME_INTERPOLATION.m
HPILR0	water heads (m) at the start of the current time step, should be recieved from MODFLOW	STEMMUS_SCOPE.m & TIME_INTERPOLATION.m
Q3DF	indicator for quasi-3d simulation, 1 means yes; 0 means no	STEMMUS_SCOPE.m & TIME_INTERPOLATION.m
ADAPTF	indicator for adaptive lower boundary setting, 1 means moving lower boundary; 0 means fixed lower boundary	STEMMUS_SCOPE.m & TIME_INTERPOLATION.m
HPUNIT	= 100, unit conversion from m to cm (MODFLOW values in m)	STEMMUS_SCOPE.m & getModelSettings.m & define_constants.m
XElemnt	layer thickness of STEMMUS soil layers	STEMMUS_SCOPE.m
NN	total number of soil layers in STEMMUS	STEMMUS_SCOPE.m & getModelSettings.m
NPILR
ITERQ3D
IGRID
RELAXF
hBOT	bottom water heads	STEMMUS_SCOPE.m
IBOTM	index of bottom soil layer	STEMMUS_SCOPE.m

[MostafaGomaa93]

@Lianyu-Yu Could you also explain the following: 1- what is the purpose of the following lines? why there are two XElenmt and XElemnt1?

XElemnt=zeros(NN,1); XElemnt1=zeros(NN,1); XElemnt1(1)=0; TDeltZ=flip(DeltZ); for ML=2:NL XElemnt1(ML)=XElemnt1(ML-1)+TDeltZ(ML-1); end XElemnt1(NN)=Tot_Depth; TOPELEV=[3840.000 3717.500 3645.161 3586.464 3560.000 3611.383 3520.358 3483.765 3454.116 3427.815 3397.764 3414.816 3397.583 3429.178 3385.020 3393.255 3810.000 3743.777 3821.000 3772.387 3665.000 3649.507 3567.249 3545.872 3502.934 3474.846 3566.573 3599.505 3559.742 3511.335 3459.323 3395.369 3394.994 3386.328 3861.316 3775.572 3662.498 3496.611 3431.786 3461.962 3409.241 3399.618 3395.802 3388.953 ].*HPUNIT; % elevation at the top surface, added by LY XElemnt=XElemnt1;%TOPELEV(IP0STm)

2- In the TIME_INTERPOLATION.m -> the function is as : function [HPILLAR,IBOT]=TIME_INTERPOLATION(T1,T0,IP,BOT,HPILLAR1N,HPILLAR0,T,DELT,X,NUMNP,Q3DF,ADAPTF)

later when you use the function in the Mainloop11.m, you did the following: [HPILLAR,IBOT]=TIME_INTERPOLATION(TEND1,TTIMEx86400,IP0STm,BOTm.xHPUNIT,HPILR1NxHPUNIT,HPILR0xHPUNIT,TIME,TTENDx86400,XElemnt,NN,Q3DF,ADAPTF)

There are a lot of time variables (T1, T0, T, TEND, TEND1, TTIME, TIME). what are the differences between them?

3- In the TIME_INTERPOLATION.m, I noticed that T1 is not used in any of the calculations. Maybe we can remove it?

[BSchilperoort]

Hi Mostafa,

Are all of these changes essential for the coupling? In Lianyu's case he had to determine how to couple the models together, as well as exchange the information. When you do the coupling through BMI, a lot of this code is not required in Matlab anymore, as you'd code it in Python (where you have the BMI of both models.).

The water head inside STEMMUS_SCOPE has to be connected to the head of MODFLOW somehow (for a two way coupling). So perhaps it would be more efficient to focus on that. I.e. answer:

• How is the water head modeled and updated in STEMMUS_SCOPE
  • in which variable(s) is the state stored?
  • can we easily change this value?
• How is the water head modeled and updated in MODFLOW
  • can we also easily change this value
• As the time steps are different, how to interpolate between the two.

[MostafaGomaa93]

Hi Bart

Are all of these changes essential for the coupling? In Lianyu's case he had to determine how to couple the models together, as well as exchange the information. When you do the coupling through BMI, a lot of this code is not required in Matlab anymore, as you'd code it in Python (where you have the BMI of both models.).

I am currently investigating and doing my best to differentiate between what has to be inside the code and what can be from outside using the BMI.

The water head inside STEMMUS_SCOPE has to be connected to the head of MODFLOW somehow (for a two way coupling). So perhaps it would be more efficient to focus on that. I.e. answer:

Yes, that's indeed what I am starting to add now to STEMMUS_SCOPE code. There is nowhere in the code so far to read the MODFLOW information (water level and at which layer it exists) and act upon it

• How is the water head modeled and updated in STEMMUS_SCOPE

Honestly, I don't know so far.

• in which variable(s) is the state stored?

I am not 100% sure yet but I think it is the "hh" variable

• can we easily change this value?

it is calculated through the code and not an input (as I understand)

• How is the water head modeled and updated in MODFLOW
  • can we also easily change this value

We will not change the water head in MODFLOW. The water head is updated in MODFLOW through the Recharge. The Recharge is the variable that MODFLOW needs from STEMMUS_SCOPE. Recharge is the summation of qh, qa, and qv variables. Those variables are already in the code. So, I think getting the recharge from STEMMUS_SCOPE will be through the BMI and not that many (or almost not at all) modifications in the code are needed.

• As the time steps are different, how to interpolate between the two.

It depends on the study case. Most likely, yes, because STEMMUS_SCOPE is half-hourly and MODFLOW is usually daily.

[MostafaGomaa93]

Hi @BSchilperoort, @SarahAlidoost , @yijianzeng

I had a meeting with Lianyu and got some useful feedback regarding my questions. There are 3 points I want to discuss

1. Spatial averaging: In the STEMMUS_MODFLOW of Lianyu: a number of MODFLOW cells are connected to one soil column. i do think its better we set a 1 to 1 connection (every MODFLOW cell is connected to 1 STEMMUS soil column) inside the STEMMUS_SCOPE code. Then, the spatial averaging can be done better in the BMI.

2. Temporal linear interpolation: Because MODFLOW is usually daily time step and STEMMUS is half hourly, the STEMMUS_MODFLOW code of Lianyu is doing linear interpolation in the TIME_INTERPOLATION.m function that i am trying to add now in the STEMMUS_SCOPE. However, i again think it's better this interpolation is done in the BMI and not in the code.

3- The data inputs for MODFLOW (e.g. number of MODFLOW layers, top and bottom elevations of each layer) are assigned inside the STEMMUS_MODFLOW code of Liaynu. This data will change per study area, so i think, in STEMMUS_SCOPE, its better we have a new txt file(s) to import the MODFLOW data inputs.

Based on your replies, i may need to do slight changes in the SSM_v.0.1 branch, and the pull request

[yijianzeng]

1. Spatial averaging: In the STEMMUS_MODFLOW of Lianyu: a number of MODFLOW cells are connected to one soil column. i do think its better we set a 1 to 1 connection (every MODFLOW cell is connected to 1 STEMMUS soil column) inside the STEMMUS_SCOPE code. Then, the spatial averaging can be done better in the BMI.

Hi Mostafa, first of all, i would suggest not to set a one-to-one connection, main reasons are following:

• STEMMUS-SCOPE's running can take a long time for a catchment, which will become a major bottleneck;
• lumping a number of MODFLOW cells based on the zonation concept is reasonable, considering groundwater table distribution will not differ much from one cell to another cell, in the same zone (this is also the criteria to do zonation)
• lumping a number of MODFLOW cells and connect then to one soil column can be also done in the BMI

[yijianzeng]

2. Temporal linear interpolation: Because MODFLOW is usually daily time step and STEMMUS is half hourly, the STEMMUS_MODFLOW code of Lianyu is doing linear interpolation in the TIME_INTERPOLATION.m function that i am trying to add now in the STEMMUS_SCOPE. However, i again think it's better this interpolation is done in the BMI and not in the code.

Time Interpolation should be certainly done in BMI, and this is the strength of BMI and also discussed with Bart last time.

[yijianzeng]

3- The data inputs for MODFLOW (e.g. number of MODFLOW layers, top and bottom elevations of each layer) are assigned inside the STEMMUS_MODFLOW code of Liaynu. This data will change per study area, so i think, in STEMMUS_SCOPE, its better we have a new txt file(s) to import the MODFLOW data inputs.

Yes, indeed.

[MostafaGomaa93]

1. Spatial averaging: In the STEMMUS_MODFLOW of Lianyu: a number of MODFLOW cells are connected to one soil column. i do think its better we set a 1 to 1 connection (every MODFLOW cell is connected to 1 STEMMUS soil column) inside the STEMMUS_SCOPE code. Then, the spatial averaging can be done better in the BMI.

Hi Mostafa, first of all, i would suggest not to set a one-to-one connection, main reasons are following:

• STEMMUS-SCOPE's running can take a long time for a catchment, which will become a major bottleneck;
• lumping a number of MODFLOW cells based on the zonation concept is reasonable, considering groundwater table distribution will not differ much from one cell to another cell, in the same zone (this is also the criteria to do zonation)
• lumping a number of MODFLOW cells and connect then to one soil column can be also done in the BMI

Hi Yijian, yeah I am aware of the computational running time problem. My point is that we set a 1-1 connection inside the code and we do the spatial averaging outside in the BMI. It will be more flexible to change the number of MODFLOW cells that are connected to one STEMMUS soil column. For ex, Lianyu set 44 MODFLOW cells that is connected to one soil column. I think this 44 can be changed based on the application conditions (e.g. size of the study area, size of pixel, and so on), so I think it's better to make these settings in the BMI and not in the code.

[yijianzeng]

Hi Yijian, yeah I am aware of the computational running time problem. My point is that we set a 1-1 connection inside the code and we do the spatial averaging outside in the BMI. It will be more flexible to change the number of MODFLOW cells that are connected to one STEMMUS soil column. For ex, Lianyu set 44 MODFLOW cells that is connected to one soil column. I think this 44 can be changed based on the application conditions (e.g. size of the study area, size of pixel, and so on), so I think it's better to make these settings in the BMI and not in the code.

Hi Mostafa, i am not sure if we are talking about the same thing here now. Particularly, when you mentioned "44 MODFLOW cells that is connected to one soil column", this is not the case. '44 MODFLOW' zones mean 44 soil profiles to be simulated by STEMMUS for the Maqu catchment. And, for each 'MODFLOW' zone, there are columns and rows of cells, the GW table and elevation of these cells are then averaged into one value to be communicated further with ONE soil profile to be simulated by STEMMUS.

[MostafaGomaa93]

Hi Yijian, yeah I am aware of the computational running time problem. My point is that we set a 1-1 connection inside the code and we do the spatial averaging outside in the BMI. It will be more flexible to change the number of MODFLOW cells that are connected to one STEMMUS soil column. For ex, Lianyu set 44 MODFLOW cells that is connected to one soil column. I think this 44 can be changed based on the application conditions (e.g. size of the study area, size of pixel, and so on), so I think it's better to make these settings in the BMI and not in the code.

Hi Mostafa, i am not sure if we are talking about the same thing here now. Particularly, when you mentioned "44 MODFLOW cells that is connected to one soil column", this is not the case. '44 MODFLOW' zones mean 44 soil profiles to be simulated by STEMMUS for the Maqu catchment. And, for each 'MODFLOW' zone, there are columns and rows of cells, the GW table and elevation of these cells are then averaged into one value to be communicated further with ONE soil profile to be simulated by STEMMUS.

Okay, thanks for the clarification. So again this averaging can be different from one to another study case (based on size of the study area, size of model cells, and so on). So, averaging them in the BMI will be more flexible in my opinion. Also if the study area is so small, we may not need any averaging

[yijianzeng]

Okay, thanks for the clarification. So again this averaging can be different from one to another study case (based on size of the study area, size of model cells, and so on). So, averaging them in the BMI will be more flexible in my opinion. Also the study area is so small, we may not need any averaging

Yes, this is correct. We should do it in BMI. The size of study area and how the zonation is done is another topic though,

[SarahAlidoost]

@yijianzeng The discussion and action points: As seen in the schematic diagram of models coupling here, there are two exchanges of information:

• From modeflow to stemmus: water table depth for all stemmus soil profiles: this can be done via set_value bmi function
• From stemmus to modflow: Recharge fluxes. this can be done via get_value bmi function

However, the variable water table depth is not in the stemmus_scope code. We need to figure out which variables are affected by this variable and then expose them to the bmi interface. Another issue is the difference between the temporal resolution of the two models (modflow resolution is daily and stemmus_scope is half hourly). The interpolation can be implemented in a python script and then called where bmi interface is used to interact with models. Therefore, if variables are needed for interpolation and are in stemmus_scope (even if they are static), they should be exposed to the bmi interface, too.

Here is the list of action points:

• Defining the variable names to be exposed like this table
• Developing a python script for temporal interpolation and other conversions
• Creating a simple workflow of coupling e.g. in a notebook
• Running the model coupling for a bigger scope spatially and temporally

@MostafaGomaa93 @BSchilperoort please feel free to edit this comment.

[MostafaGomaa93]

Hi @SarahAlidoost, @BSchilperoort and @yijianzeng

By Following Lianyu's STEMMUS_MODFLOW code, and trying to know which variables will be affected by the water heads (received from MODFLOW), Please have a look below:

1- The received MODFLOW water heads will be used to calculate new variables "hBOT" and "IBOTM" (TIME_INTERPOLATION.m) and here). We agreed that these calculations will be shifted into a Python script (outside the code).

2- These new variables ("hBOT" and "IBOTM") are then used to define the boundary conditions. previously was called "h_BC.m", and now refactored to "calculateBoundaryConditions.m". Those calculations are in lines (22-34) of Linayu's code.

3- The comparison of the boundary conditions calculations is as follow (can also be found here)	code	Lianyu's STEMMUS_MODFLOW	STEMMUS_SCOPE
Equivalent lines	lines 9-19	lines 18-27 with refactorated names
non-equivalent lines	lines 22-34	doesn't exist

4- For the equivalent lines between the two codes, MODFLOW variables are not used (default). For the non-equivalent lines, MODFLOW variables are used. The affected variables are "RHS", "C4", "C4_a". Honestly, so far, I don't know what does "RHS", "C4", and "C4_a" mean. But what I understand is that those variables need to be exposed to the BMI, right?

5- Then those "RHS", "C4", "C4_a" shoud be added to the table, right?

A further question, what do @BSchilperoort and @SarahAlidoost propose about these changes in the calculations of RHS and C4? Do we need to add those calculations in the code itself (in the calculateBoundaryConditions.m)? or we can do it also in a separate Python script?

[MostafaGomaa93]

I would like also to add that the same variables "hBOT" and "IBOTM" are used in another file "Diff_Moisture_Heat.m" of Lianyu's code, and seems like it also do changes for the "RHS" and "C4" variables (lines 232-242). However, this "Diff_Moisture_Heat.m" file is not in the STEMMUS_SCOPE code

[SarahAlidoost]

Hi @SarahAlidoost, @BSchilperoort and @yijianzeng

By Following Lianyu's STEMMUS_MODFLOW code, and trying to know which variables will be affected by the water heads (received from MODFLOW), Please have a look below:

1- The received MODFLOW water heads will be used to calculate new variables "hBOT" and "IBOTM" (TIME_INTERPOLATION.m) and here). We agreed that these calculations will be shifted into a Python script (outside the code).

yes, all the inputs of the function TIME_INTERPOLATION.m as you implemented here are either model settings or modflow variables except variable TIME which is stemmus_scope variable and it gets updated in the main loop. So it should be exposed to BMI. As we discussed this function TIME_INTERPOLATION.m can be implemented outside of the model repository.

2- These new variables ("hBOT" and "IBOTM") are then used to define the boundary conditions. previously was called "h_BC.m", and now refactored to "calculateBoundaryConditions.m". Those calculations are in lines (22-34) of Linayu's code.

Looking at the script, IBOTM is not used anywhere. So, let's ignore it. The two variables hBOT and IBOT are used in h_BC.m.

3- The comparison of the boundary conditions calculations is as follow (can also be found here)

code Lianyu's STEMMUS_MODFLOW STEMMUS_SCOPE Equivalent lines lines 9-19 lines 18-27 with refactorated names non-equivalent lines lines 22-34 doesn't exist

Looking at the implementation h_BC.m, I would say these calculations do not exist in the calculateBoundaryConditions.m. note if ~Q3DF and index of RHS.

4- For the equivalent lines between the two codes, MODFLOW variables are not used (default). For the non-equivalent lines, MODFLOW variables are used. The affected variables are "RHS", "C4", "C4_a". Honestly, so far, I don't know what does "RHS", "C4", and "C4_a" mean. But what I understand is that those variables need to be exposed to the BMI, right?

The meaning of RHS is defined here. Variables "RHS", "C4", and "C4_a" are used in many places in the code. I would say exposing hBOT and IBOT are more straightforward as they are only used in this script.

5- Then those "RHS", "C4", "C4_a" shoud be added to the table, right?

no, it is better to add hBOT and IBOT. see my previous comment.

A further question, what do @BSchilperoort and @SarahAlidoost propose about these changes in the calculations of RHS and C4? Do we need to add those calculations in the code itself (in the calculateBoundaryConditions.m)? or we can do it also in a separate Python script?

The changes are about NBChB and calculating RHS, C4, ... Now, we have two different calculations:

• Case A: when there is modflow- consider these lines in h_BC,
• Case B: when there is no modflow- consider these lines, in calculateBoundaryConditions.

I would say keep the implementation in MATLAB. There are two ways to implement this:

1- replace case B completely with case A: then variables Q3DF, hBOT and IBOT should be defined perhaps by some default values when there is no modflow. 2- or we keep both cases with/without modflow. @yijianzeng any suggestions on which way to choose?

[SarahAlidoost]

I would like also to add that the same variables "hBOT" and "IBOTM" are used in another file "Diff_Moisture_Heat.m" of Lianyu's code, and seems like it also do changes for the "RHS" and "C4" variables (lines 232-242). However, this "Diff_Moisture_Heat.m" file is not in the STEMMUS_SCOPE code

The Diff_Moisture_Heat.m script is only used if CPLD==0 see here. CPLD is defined here that does not exist in STEMMUS_SCOPE. @yijianzeng can you please help to address this issue.

Meantime, the calculations in Diff_Moisture_Heat.m script are very similar to calculateHydraulicConductivity.

[MostafaGomaa93]

The changes are about NBChB and calculating RHS, C4, ... Now, we have two different calculations:

• Case A: when there is modflow- consider these lines in h_BC,
• Case B: when there is no modflow- consider these lines, in calculateBoundaryConditions.

I would say keep the implementation in MATLAB. There are two ways to implement this:

1- replace case B completely with case A: then variables Q3DF, hBOT and IBOT should be defined perhaps by some default values when there is no modflow.

@SarahAlidoost I am working now on updating the calculateBoundaryConditions.m

Two questions

1. For the Q3DF: a- I would like to change it from "Q3DF" to "MODFLOW" (more meaningful) b- Can we move this new variable "MODFLOW" to the config_file? with Default value = 0 -> means no MODFLOW coupling and "MODFLOW" = 1, for MODFLOW coupling
2. Can we make a new function (called for e.g. loadMODFLOWData) with a new input txt file? In that txt file, the user will define all related MODFLOW settings (e.g. a number of MODFLOW layers, top and bottom of each layer, etc)?

[SarahAlidoost]

The changes are about NBChB and calculating RHS, C4, ... Now, we have two different calculations:

• Case A: when there is modflow- consider these lines in h_BC,
• Case B: when there is no modflow- consider these lines, in calculateBoundaryConditions.

I would say keep the implementation in MATLAB. There are two ways to implement this: 1- replace case B completely with case A: then variables Q3DF, hBOT and IBOT should be defined perhaps by some default values when there is no modflow.

You are going to replace case B completely with case A, right?

@SarahAlidoost I am working now on updating the calculateBoundaryConditions.m

Two questions

1. For the Q3DF: a- I would like to change it from "Q3DF" to "MODFLOW" (more meaningful)

Variable names should be meaningful but also self-explanatory as much as possible for example ModflowCoupling. Considering the implementing case B, it is even possible to ignore the if~Q3DF statement and define default values for hBOT and IBOT in case of no modflow coupling. Later, these two variables can be set by BMI outside the model code for coupling purposes. Therefore, the question is: what are the values of hBOT and IBOT in case of no modflow coupling?

b- Can we move this new variable "MODFLOW" to the config_file? with Default value = 0 -> means no MODFLOW coupling and "MODFLOW" = 1, for MODFLOW coupling

no need for that. If needed, I suggest adding for example, a variable ModflowConfigPath to the config_file where ModflowConfigPath is a path to modflow config file, see below.

2. Can we make a new function (called for e.g. loadMODFLOWData) with a new input txt file? In that txt file, the user will define all related MODFLOW settings (e.g. a number of MODFLOW layers, top and bottom of each layer, etc)?

It depends what the function does: loading data or reading settings/configurations? Following your explanations, I think you want a function only for all related modflow settings. Therefore, you can define a function e.g. readModflowConfigs.m where it reads the text file Modflow_config.txt and sets all required variables. Here is an example function that reads a text file read_config.m and sets some variables. However, the question is where modflow settings (e.g. a number of MODFLOW layers, top and bottom of each layer, etc) are used in stemmus_scope?

[MostafaGomaa93]

You are going to replace case B completely with case A, right?

Yes

Variable names should be meaningful but also self-explanatory as much as possible for example ModflowCoupling.

Okay, thanks for the suggestion

Considering the implementing case B, it is even possible to ignore the if~Q3DF statement and define default values for hBOT and IBOT in case of no modflow coupling. Later, these two variables can be set by BMI outside the model code for coupling purposes. Therefore, the question is: what are the values of hBOT and IBOT in case of no modflow coupling?

If no MODFLOW coupling, hBOT and IBOT won't be used. so, whatever values assigned, they won't be used

b- Can we move this new variable "MODFLOW" to the config_file? with Default value = 0 -> means no MODFLOW coupling and "MODFLOW" = 1, for MODFLOW coupling

no need for that. If needed, I suggest adding for example, a variable ModflowConfigPath to the config_file where ModflowConfigPath is a path to modflow config file, see below.

My idea is that from the beginning, the user is aware if he needs MODFLOW coupling or not. So, the user won't bother later by searching inside the code to activate MODFLOW coupling, and simply activate MODFLOW coupling from the config file.

2. Can we make a new function (called for e.g. loadMODFLOWData) with a new input txt file? In that txt file, the user will define all related MODFLOW settings (e.g. a number of MODFLOW layers, top and bottom of each layer, etc)?

It depends what the function does: loading data or reading settings/configurations? Following your explanations, I think you want a function only for all related modflow settings. Therefore, you can define a function e.g. readModflowConfigs.m where it reads the text file Modflow_config.txt and sets all required variables. Here is an example function that reads a text file read_config.m and sets some variables. However, the question is where modflow settings (e.g. a number of MODFLOW layers, top and bottom of each layer, etc) are used in stemmus_scope?

Yes, indeed I need a function to define all MODFLOW related settings. For example, this "BOTm", which is MODFLOW data input, refers to bottom level of all layers, and will be used in the calculations of the "RHS". if no MODFLOW coupling, this new function won't be used at all. So, I will prepare a template for this function which works only if ModflowCoupling = True

[yijianzeng]

Hi @SarahAlidoost, @BSchilperoort and @yijianzeng By Following Lianyu's STEMMUS_MODFLOW code, and trying to know which variables will be affected by the water heads (received from MODFLOW), Please have a look below: 1- The received MODFLOW water heads will be used to calculate new variables "hBOT" and "IBOTM" (TIME_INTERPOLATION.m) and here). We agreed that these calculations will be shifted into a Python script (outside the code).

yes, all the inputs of the function TIME_INTERPOLATION.m as you implemented here are either model settings or modflow variables except variable TIME which is stemmus_scope variable and it gets updated in the main loop. So it should be exposed to BMI. As we discussed this function TIME_INTERPOLATION.m can be implemented outside of the model repository.

Hi Sarah, Mostafa, to me, 'TIME_INTERPOLATION.m' indeed can be implemented outside of the model repository. On top of what has been discussed, i think TIME_INTERPOLATION.m can be coded directly in BMI, i think it is in that sense, Sarah mentioned that "Developing a python script for temporal interpolation and other conversions"

[yijianzeng]

3- The comparison of the boundary conditions calculations is as follow (can also be found here) code Lianyu's STEMMUS_MODFLOW STEMMUS_SCOPE Equivalent lines lines 9-19 lines 18-27 with refactorated names non-equivalent lines lines 22-34 doesn't exist

Looking at the implementation h_BC.m, I would say these calculations do not exist in the calculateBoundaryConditions.m. note if ~Q3DF and index of RHS.

What is mentioned here by Mostafa is to incorporate 'IBOT' & 'hBOT' into the calculateBoundaryConditions.m.

[yijianzeng]

The changes are about NBChB and calculating RHS, C4, ... Now, we have two different calculations:

• Case A: when there is modflow- consider these lines in h_BC,
• Case B: when there is no modflow- consider these lines, in calculateBoundaryConditions.

I would say keep the implementation in MATLAB. There are two ways to implement this:

1- replace case B completely with case A: then variables Q3DF, hBOT and IBOT should be defined perhaps by some default values when there is no modflow. 2- or we keep both cases with/without modflow. @yijianzeng any suggestions on which way to choose?

Great thanks, Sarah. I fully agree with you observations! And i would suggest to keep both cases, e.g., we can switch between with/without modflow.

[yijianzeng]

I would like also to add that the same variables "hBOT" and "IBOTM" are used in another file "Diff_Moisture_Heat.m" of Lianyu's code, and seems like it also do changes for the "RHS" and "C4" variables (lines 232-242). However, this "Diff_Moisture_Heat.m" file is not in the STEMMUS_SCOPE code

The Diff_Moisture_Heat.m script is only used if CPLD==0 see here. CPLD is defined here that does not exist in STEMMUS_SCOPE. @yijianzeng can you please help to address this issue.

Meantime, the calculations in Diff_Moisture_Heat.m script are very similar to calculateHydraulicConductivity.

The 'Diff_Moisture_Heat.m' is a simplified 'STEMMUS' for numerically solving 'Richards' equation only. That is why 'CPLD==0' (meaning no coupling between soil water and heat transport). We need to confirm with Lianyu, if the 'Diff_Moisture_Heat.m' is used in STEMMUS-MODFLOW or not.

@Lianyu-Yu please respond here.

[yijianzeng]

1. For the Q3DF: a- I would like to change it from "Q3DF" to "MODFLOW" (more meaningful)

Hi Mostafa, 'Q3DF' could be referring to 'Quasi-3D Flow'. @Lianyu-Yu to confirm.

[yijianzeng]

@SarahAlidoost I am working now on updating the calculateBoundaryConditions.m Two questions

1. For the Q3DF: a- I would like to change it from "Q3DF" to "MODFLOW" (more meaningful)

Variable names should be meaningful but also self-explanatory as much as possible for example ModflowCoupling. Considering the implementing case B, it is even possible to ignore the if~Q3DF statement and define default values for hBOT and IBOT in case of no modflow coupling. Later, these two variables can be set by BMI outside the model code for coupling purposes. Therefore, the question is: what are the values of hBOT and IBOT in case of no modflow coupling?

Hi Sarah, if no coupling, 'hBOT' and 'IBOT' are not needed, and can be set to zero or null.

And i agree 'ModflowCoupling' is more meaningful.

[yijianzeng]

However, the question is where modflow settings (e.g. a number of MODFLOW layers, top and bottom of each layer, etc) are used in stemmus_scope?

The modflow settings are not used in STEMMUS-SCOPE.

[MostafaGomaa93]

@yijianzeng and @Lianyu-Yu

I have uploaded the example (mh_9_1_1) of Linayu STEMMUS_MODFLOW in a new branch because I have some related questions. The answers will affect the proposed changes in STEMMUS_SCOPE code (SSM_v.0.2)

1. In the h_BC.m (line 19) -> what does it mean by INBT ?
2. In the example “mh_9_1_1”, IPOSTm is not assigned I can't find that in the MainLoop00). Based on a previous discussion, issue 212, IPOSTm means -> number of STEMMUS columns for MODFLOW. Still not clear to me what does it mean by IPOSTm and what is the value of IPOSTm?
3. What does the BOTM mean? In the model example “mh_9_1_1”, there is a conflict about the value of BOTM as follows: o In the main code MainLoop00 (line 141) -> BOTM is retrieved from the external “Bot_nodes.txt” o In the h_BC.m (line 20) -> BOTM = XElemnt(NN)

[SarahAlidoost]

@MostafaGomaa93 I have some questions about TIME_INTERPOLATION.m. When you have time, please have a look, thanks:

1. in this line gwTable = botmLayerLevel(1, nSoilColumns) - headBotmLayer, what are botmLayerLevel and nSoilColumns?
2. in this line indexBotmLayer = indexBotmLayer + 3, what is the reason for adding 3?
3. is there any reference in literature for this code block?