Rm unused variables and equations

[SarahAlidoost]

Description

relates #240

🔴 exe should be updated by merging the main into this branch.

Checklist

• Add a reference to related issues.
• @mentions of the person or team responsible for reviewing proposed changes.
• This pull request has a descriptive title.
• Code is written according to the guidelines.
• The checks by MISS_HIT style checker and linter, below the pull request, are successful (green).
• Documentation is available.
• Add changes to the changelog file under section Unreleased.
• Model runs successfully.
• Ask a meinatainer to re-generate exe file if matlab codes are changed. About how to create an exe file, see exe readme.

[SarahAlidoost]

@MostafaGomaa93 in this PR, unused variables and equations are commented out, only botmSoilLevel is removed. see Files changed.

@BSchilperoort the exe file is re-genertaed. please let me know if more fixes are needed.

[BSchilperoort]

I found it difficult to understand the current code.

Why not reduce the initialization to:

function GroundwaterSettings = readGroundwaterSettings()
    % Activate/deactivate Groundwater coupling
    GroundwaterSettings.GroundwaterCoupling = 0; % (value = 0 -> deactivate coupling, or = 1 -> activate coupling); default = 0, update value to = 1 -> through BMI

    % Initialize the variables (head, temperature, air pressure) at the bottom boundary (start of saturated zone)
    GroundwaterSettings.tempBotm = 0; % groundwater temperature (C), received from MODFLOW through BMI

    % Call MODFLOW layers information (number of aquifer layers and their elevations, etc)
    GroundwaterSettings.gw_Dep = 1.0; % elevation of the top surface aquifer layer
    %     to avoid model crashing, assign minimum value of 1 cm
end

Move the soil thickness calculation to a more appropriate place. It does nothing with the groundwater inputs:

    % Calculate soil layers thickness (cumulative layers thickness; e.g. 1, 2, 3, 5, 10, ......., 480, total soil depth)
    soilThick = zeros(NN, 1); % cumulative soil layers thickness
    soilThick(1) = 0;
    DeltZ = ModelSettings.DeltZ;
    DeltZ_R = ModelSettings.DeltZ_R;

    % Load model settings
    ModelSettings = io.getModelSettings();
    NN = ModelSettings.NN; % Number of nodes
    NL = ModelSettings.NL; % Number of layers

    for i = 2:NL
        soilThick(i) = soilThick(i - 1) + DeltZ_R(i - 1);
    end
    soilThick(NN) = ModelSettings.Tot_Depth; % total soil depth

---

Also, the indxBotmLayer thing is a bit problematic. It is calculated in readGroundwaterSettings.m based on a groundwater depth value, but this is only run once at initialization, even though it's supposed to skip saturated layers:

        % index of bottom layer after neglecting saturated layers (from bottom to top)
        indxBotm = GroundwaterSettings.indxBotmLayer;

So, shouldn't the following code be run every time step...?

    for i = 1:NL
        midThick = (soilThick(i) + soilThick(i + 1)) / 2;
        if gw_Dep >= soilThick(i) && gw_Dep < soilThick(i + 1)
            if gw_Dep < midThick
                indxBotmLayer_R = i;
            elseif gw_Dep >= midThick
                indxBotmLayer_R = i + 1;
            end
            break
        elseif gw_Dep >= soilThick(i + 1)
            continue
        end
    end

---

If these things are fixed, then the BMI implementation will be simple: only three variables (coupling toggle, temperature of the groundwater, groundwater head at top aquifer) have to be implemented.

[SarahAlidoost]

@MostafaGomaa93 can you check the calculation of indxBotmLayer and if it should be calculated every time step?

[MostafaGomaa93]

I found it difficult to understand the current code.

Why not reduce the initialization to

Move the soil thickness calculation to a more appropriate place. It does nothing with the groundwater inputs:

Hi @BSchilperoort , its fine for me to split this function into two as you suggested if that will increase the readability of the code

[MostafaGomaa93]

So, shouldn't the following code be run every time step...?

    for i = 1:NL
        midThick = (soilThick(i) + soilThick(i + 1)) / 2;
        if gw_Dep >= soilThick(i) && gw_Dep < soilThick(i + 1)
            if gw_Dep < midThick
                indxBotmLayer_R = i;
            elseif gw_Dep >= midThick
                indxBotmLayer_R = i + 1;
            end
            break
        elseif gw_Dep >= soilThick(i + 1)
            continue
        end
    end

Yes, gw_Dep has to be calculated everytime step based on the updated GroundwaterSettings.headBotmLayer

[MostafaGomaa93]

@MostafaGomaa93 can you check the calculation of indxBotmLayer and if it should be calculated every time step?

Yes, indxBotmLayer has to be updated everytime step of MODFLOW

[SarahAlidoost]

@MostafaGomaa93 can you please review the changes in this pull request? Feel free to suggest better names for functions/variables, if needed. and also double-check the logic in this code block. Thanks.

[MostafaGomaa93]

I did several tests to track where the problem was. First of all, this strange temperature behaviour happens in the soil layer above the groundwater table and only happens when we set groundwater temperature > temperature of the layer above the groundwater table. So, when I was testing with groundwater temperature = 17 C, there was no problem. However, when @BSchilperoort changed the value to 23 C, we observed this strange behaviour. The same happens when the groundwater temperature = NaN because in this case, the groundwater temperature is replaced by the bottom temperature boundary which = 22.96 C (same case as when groundwater temperature = 23 C).

So, in the code, this strange behaviour is related to the RHS(indxBotm + 1) which is calculated and updated based on different variables in several lines (here, here and here).

Honestly, I tried several things but doesn't solve the problem. Also, the solution suggested by @BSchilperoort doesn't solve it. I came up with an end-up solution to solve the problem (suggested here, not the ideal one so far). A better solution will need a lot of tests and take more time.

From the physics point of view, an important question now is, Does it happen in reality that groundwater temperature > soil temperature? I was thinking that this is not the case and that groundwater is always cooler than the soil. But I need to discuss this further with @yijianzeng. I also sent to one of the MODFLOW development team members about this and will update you.

This is now the temperature figure with groundwater temperature = 23 C

[SarahAlidoost]

@MostafaGomaa93 can you please address the comments above? Then I can create the exe file and test whether the model runs with/without coupling.

[MostafaGomaa93]

@MostafaGomaa93 can you please address the comments above? Then I can create the exe file and test whether the model runs with/without coupling.

Done

[BSchilperoort]

I'm afraid this did not fully fix the issue... While there is no negative spike anymore, a growing positive spike does form over time:

[SarahAlidoost]

I'm afraid this did not fully fix the issue... While there is no negative spike anymore, a growing positive spike does form over time:

@BSchilperoort and @MostafaGomaa93 It seems that fixing this issue of temperature needs more work and tests and this pull request has expanded beyond just addressing unused variables. Since the BMI functionality is now working, I prefer to merge this pull request and the BMI one in Pystemmusscope. The other issue should be resolved in a separate pull request.

[MostafaGomaa93]

Since the BMI functionality is now working, I prefer to merge this pull request and the BMI one in Pystemmusscope. The other issue should be resolved in a separate pull request.

I agree with @SarahAlidoost

[BSchilperoort]

I prefer to merge this pull request and the BMI one in Pystemmusscope. The other issue should be resolved in a separate pull request.

We can merge this PR, however, I don't think it's a good idea to make a release yet until the issue has been fixed.

We cannot merge a working PR on PyStemmusScope without there being a new STEMMUS_SCOPE container image (which ideally requires a new release).

[yijianzeng]

From the physics point of view, an important question now is, Does it happen in reality that groundwater temperature > soil temperature? I was thinking that this is not the case and that groundwater is always cooler than the soil. But I need to discuss this further with @yijianzeng

Dear all, thanks a lot for your great efforts into this. Indeed, physically speaking for the soil bottom when it is emerged in groundwater, the soil temperature should be equal to the groundwater temperature. As such, if you have an initial soil temperature like 17-19 degree or so, the groundwater temperature should be set equal to that value.

To consider the heat flow of groundwater is indeed another topic to tackle with, since this is linked to MODFLOW's heat flow version, which is not yet released yet as i understood.

For now, I would suggest fixing the groundwater temperature = soil temperature wherever GW reaches a certain depth of the soil column.