Values of soil temperature when groundwater coupling

[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:

Originally posted by @BSchilperoort in #239 (comment)

[SarahAlidoost]

🔴 This issue should be fixed before releasing stemmus_scope, see details in #239

[SarahAlidoost]

@BSchilperoort I found the reason for the artifact in your plot, the value of groundwater_head_bottom_layer is set to 1750, as model.set_value("groundwater_head_bottom_layer", np.array([2000-250.])) # 250 cm under ground surface, in your notebook see cell 5 in the notebook. if it is changed to default values in stemmus_scope i.e. 1950.0, the issue is gone

@MostafaGomaa93 is there any condition for the values of groundwater_head_bottom_layer?

[SarahAlidoost]

I also noticed the implementation of SoilVariable.TT and TT is not ideal in STEMMUS_SCOEP.m. However, this is a separate issue. It would be better if TT is replaced by SoilVariable.TT everywhere in the script, the same for these variables.

[BSchilperoort]

if it is changed to default values in stemmus_scope i.e. 1950.0, the issue is gone

The issue is still there, it occurs at the groundwater depth. See:

If the groundwater level was set to 4 m under the surface, the issue would occur there.

[SarahAlidoost]

These are the values of the test that creates the plots:

Location=ZA-Kru
StartTime=2001-01-01T00:00
EndTime=2001-01-01T10:30

model.set_value("groundwater_coupling_enabled", np.array([True]))
model.set_value("groundwater_elevation_top_aquifer", np.array([2000.]))
model.set_value("groundwater_head_bottom_layer", np.array([1950.]))  # 50 cm under ground surface
model.set_value("groundwater_temperature", np.array([23.])) # optional. 

[MostafaGomaa93]

@SarahAlidoost @yijianzeng

After our meeting today, I test again by replacing the tempBotm in this line by tempBotm = SoilVariables.TT(indxBotm) which is = 22.96 for the ZA-Kru and the problem is still there.

[yijianzeng]

          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:

This depends on the temperature value you set for groundwater temperature. See my other comments here.

[MostafaGomaa93]

          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:

This depends on the temperature value you set for groundwater temperature. See my other comments here.

Yes, I mean your suggestion by fixing the groundwater temperature = soil temperature wherever GW reaches a certain depth of the soil column. didn't solve the problem yet. So, it looks like there is a bug indeed

[MostafaGomaa93]

Please have a look at these two sheets at the column with first value = 230 (depth = 230), the depth above the groundwater table (groundwater depth = 250). In both cases, you will find the values at that column are lower than the two columns before and after that column (columns 220 and 245) which case these strange peaks in @BSchilperoort figures

case 1: setting bottom boundary = groundwater temperature = 23
Sim_Temp_gwtemp=23.csv

case 2: setting bottom temperature boundary = soil temperature, which is = 22.96
Sim_Temp_gwtemp=soiltemp.csv

[yijianzeng]

@MostafaGomaa93 is there any condition for the values of groundwater_head_bottom_layer?

From this plot, it seems to me at 50cm the groundwater temperature of 23 degree is used to replace soil temperature profile below 50cm?

[yijianzeng]

I also noticed the implementation of SoilVariable.TT and TT is not ideal in STEMMUS_SCOEP.m. However, this is a separate issue. It would be better if TT is replaced by SoilVariable.TT everywhere in the script, the same for these variables.

'TT' is the soil temperature updated after the current time step running. And i agree with you that this should be standardized.

[yijianzeng]

The issue is still there, it occurs at the groundwater depth. See:

I would say this is not the issue, but the artifact induced by manually setting the groundwater temperature.

[MostafaGomaa93]

I would say this is not the issue, but the artifact induced by manually setting the groundwater temperature.

The issue also occurs if we set groundwater temperature = soil temperature of bottom layer and the issue is there at any location of the groundwater table. see another example with a groundwater table is 3.5 m below the surface

[yijianzeng]

@SarahAlidoost @yijianzeng

After our meeting today, I test again by replacing the tempBotm in this line by tempBotm = SoilVariables.TT(indxBotm) which is = 22.96 for the ZA-Kru and the problem is still there.

Hi Mostafa, please try the following:

• first, check the soil temperature profile without coupling the groundwater flow.
• second, check the soil temperature profile after coupling the groundwater flow.
• third, try to see if you are sure that at the depth of the groundwater table, you did not change the soil temperature at that depth

[yijianzeng]

I would say this is not the issue, but the artifact induced by manually setting the groundwater temperature.

The issue also occurs if we set groundwater temperature = soil temperature of bottom layer and the issue is there at any location of the groundwater table. see another example with a groundwater table is 3.5 m below the surface

Another point i want to mention is that, in STEMMUS-SCOPE, the moisture and heat flow are coupled. So when coupled with groundwater flow, the moisture content of the bottom soil layer will increase, which will impact soil temperature simulation as well, see below:

Zhao, H., Zeng, Y., Lv, S., and Su, Z.: Analysis of soil hydraulic and thermal properties for land surface modeling over the Tibetan Plateau, Earth Syst. Sci. Data, 10, 1031–1061, https://doi.org/10.5194/essd-10-1031-2018, 2018.

[MostafaGomaa93]

I think the problem is relaxed with the time. Have a look at this figure of 8-days simulation
Location=ZA-Kru
StartTime=2001-01-01T00:00
EndTime=2001-01-08T00:00
tempBotm = SoilVariables.TT(indxBotm)

This solution is without the modification in this line (once agreed, we can delete it). Also, same results if we set GroundwaterSettings.tempBotm = 23

[MostafaGomaa93]

another test with changing the groundwater depth over time (~1m drop per 8 days)

[SarahAlidoost]

The current coupling implementation has several issues related to coding and setting physical parameters. Here is a summary:

Coding issues:
1-(major) introducing indxBotm to the modules +energy, +soilmoisture and +dryair is not needed at most parts and causes bugs for example in this line index of P_gg is still 2. Other examples of such bugs are in SV.T(2), C2(1, 2).

2- setting the values of RHS in modules +energy here, +soilmoisture here are not correct. Model should calculate those values.

3- (major)The soil temperature SoilVariables.TT should be set in update loop here. But, due to current codes, the variable TT should be set too.

4- (minor issue) RHS should be defined in module +energy as RHS = InitialValues.RHS; in assembleCoefficientMatrices.m

Setting physical parameters:

1- (major)The coupling only sets soil temperature SoilVariables.TT. But, it seems that SoilVariables.T and T are also related and should be set accordingly.

2- in coupling BoundaryCondition.NBChB = 1;, which is only valid for the sites with landcover 'Croplands'. But, for example ZA-Kru is not a 'Croplands'.

3- Values of RHS(1) in +soilmoisture is calculated as headBotmLayer - topLevel + soilThick(indxBotmLayer_R). First of all the units don't match. Second, this equation might return zero.

4- (major)The condition in energy/solveEnergyBalanceEquations.m should be checked because it is using ForcingData.Tmin. This might not be valid in the case of coupling.

[SarahAlidoost]

@MostafaGomaa93 here are the codes that we discussed at our meeting.

Codes for plotting the soil temperature values:

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import xarray as xr

# Read the CSV file
csv_path = "/path_to/output/ZA-Kru_2023-09-06-1228/Sim_Temp.csv"

# Read the CSV file
df = pd.read_csv(csv_path, header=None)

# Get the first row as depths
depths = df.iloc[0].values.astype(float) * -1

# Get the rest of the data as soil_temperature
soil_temperature = df.iloc[3:].values.astype(float)

# Create a time index
time = np.arange(1, len(df)-2)

# Create the xarray DataArray
da_t = xr.DataArray(
    data=soil_temperature,
    dims=("time", "depth"),
    coords={"time": time, "depth": depths},
)

da_t.isel(time=0).plot(y="depth")
da_t.isel(time=1).plot(y="depth")
da_t.isel(time=2).plot(y="depth")

Codes for setting values through BMI:

The notebook is available in pull request EcoExtreML/STEMMUS_SCOPE_Processing#101.

# as an example, every time step i, the value of groundwater_head_bottom_layer is set to [1950.-i*4])
model.set_value("groundwater_head_bottom_layer", np.array([1950.-i*4]))

Git commands for working with BMI implementation:

Based on the pystemmusscope documentation.

git clone https://github.com/EcoExtreML/STEMMUS_SCOPE_Processing.git
cd STEMMUS_SCOPE_Processing
git checkout bmi-groundwater-coupling  # switch to BMI branch
git checkout -b name_of_your_branch  # create your own branch
pip install -e .[dev]  # install pystemmusscope in development mode
pip install jupyterlab
jupyter lab  # start jupyter lab

[MostafaGomaa93]

Thanks a lot @SarahAlidoost for the comments. here are my replies

The current coupling implementation has several issues related to coding and setting physical parameters. Here is a summary:

Coding issues: 1-(major) introducing indxBotm to the modules +energy, +soilmoisture and +dryair is not needed at most parts and causes bugs for example in this line index of P_gg is still 2. Other examples of such bugs are in SV.T(2), C2(1, 2).

Yes, those 3 are my mistakes, I have corrected them in the new pull request. However, the changes doesn't change the soil temperature results because the three changes are in the condition when ModelSettings.Soilairefc = 1, and by default its = 0 (here) (the if condition related to the 3 mistakes are here, here and here).

2- setting the values of RHS in modules +energy here, +soilmoisture here are not correct. Model should calculate those values.

The model indeed should calculate the RHS for all layers. However, when we set the location of the groundwater table as the bottom boundary of the soil, then all layers below the groundwater table are saturated layers (full with water), so there is no soil anymore and STEMMUS is not adapted to calculate the variables below the groundwater table, and thats why we will couple STEMMUS_SCOPE to MODFLOW. So, the RHS of the saturated layers in the soilmoisture, energy, and dryair modules must be set manually to certain values (as we did already).

3- (major)The soil temperature SoilVariables.TT should be set in update loop here. But, due to current codes, the variable TT should be set too.

Okay, I add them here

4- (minor issue) RHS should be defined in module +energy as RHS = InitialValues.RHS; in assembleCoefficientMatrices.m

Okay, I have add this in this line. However, it doesn't change the soil temperature results

Setting physical parameters:

1- (major)The coupling only sets soil temperature SoilVariables.TT. But, it seems that SoilVariables.T and T are also related and should be set accordingly.

I add the SoilVariables.T and T here. Is that what we need only or we need further edits?

2- in coupling BoundaryCondition.NBChB = 1;, which is only valid for the sites with landcover 'Croplands'. But, for example ZA-Kru is not a 'Croplands'.

The BoundaryCondition.NBChB = 1 means that the bottom boundary is a head boundary and it is not related to the landcover type. In principle, if we know the head at the bottom of the soil (which is the case when we use MODFLOW), then we use BoundaryCondition.NBChB = 1. I think Yijian may help to answer why for the case of 'Croplands', BoundaryCondition.NBChB = 1

3- Values of RHS(1) in +soilmoisture is calculated as headBotmLayer - topLevel + soilThick(indxBotmLayer_R). First of all the units don't match. Second, this equation might return zero.

If the BoundaryCondition.NBChB = 1, then the RHS is a head and has units of cm. If the BoundaryCondition.NBChB = 2, then the RHS is a flux and has units of cm/sec. In the case of groundwater coupling, BoundaryCondition.NBChB = 1 and RHS is a head and has units of cm and the headBotmLayer, topLevel, soilThick(indxBotmLayer_R) are all in cm. and the value of RHS(indexbotm) is indeed close to zero

4- (major)The condition in energy/solveEnergyBalanceEquations.m should be checked because it is using ForcingData.Tmin. This might not be valid in the case of coupling.

I removed this condition while testing but found that no difference in the results. Maybe Yijian could clarify what is the purpose of this condition

[MostafaGomaa93]

@MostafaGomaa93 here are the codes that we discussed at our meeting.

Codes for plotting the soil temperature values:

Thanks so much for the codes. I have used the plotting code to observe the spikes in the figures below. However, the lines in my figures are sharp and not that smooth as in @BSchilperoort figures (don't know why)

In principle, I think as we increase the model run time, spikes will be minor and/or totally disappear and the problem is relaxed

[MostafaGomaa93]

For all the figures -> groundwater temperature = NaN, so groundwater temperature = soil temperature. Spikes for timestep = 1 or = 20 are marked with grey circle. Spikes for timestep = 60 or = 300 are marked with black circle

groundwater depth = 50 cm, (average condition) is not used	groundwater depth = 150 cm, average condition is not used
groundwater depth = 150 cm, average condition is used	groundwater depth = 250 cm, average condition is not used

[MostafaGomaa93]

Another two examples at groundwater depth = 250 and (average condition) is not used

groundwater temperature = 23	groundwater temperature = 17

[SarahAlidoost]

@MostafaGomaa93 here are the codes that we discussed at our meeting.
Codes for plotting the soil temperature values:

Thanks so much for the codes. I have used the plotting code to observe the spikes in the figures below. However, the lines in my figures are sharp and not that smooth as in @BSchilperoort figures (don't know why)

The code that I provided plots the entire depth profile (0, 5 m). In your plots here, the y axis is between (-20, 2.5 cm). Please make sure that you can create the exact same plots. So we are looking at the same things. You can also generate an exe file and run the BMI notebook to make sure that plots are the same.

[SarahAlidoost]

@MostafaGomaa93 the plots below are created using your changed code in #245, for two sites: ZA-Kru "Savannas" and DE-Geb "cropland". The plots show soil temperature and soil moisture in case no coupling and coupling with different temperatures and depths. The issues are indicated in each plot. it seems that results are better for DE-Geb "cropland" than ZA-Kru "Savannas". But the results are the same for DE-Geb with different temperatures 23 and 17, e.g. see Fig 3. It would be good to check the code for one more site with "cropland".

The issue of the first time step can be fixed by setting TT in the update section, see here.

As mentioned before, in stemmusscope (without coupling), BoundaryCondition.NBChB is set to 1 for sites with "cropland" whereas 3 for other sites. However, in groundwater coupling, BoundaryCondition.NBChB is set to 1 for all sites. @yijianzeng Is this a valid assumption? see definition here.

Fig 1: ZA-Kru "Savannas"

Fig 2: DE-Geb "cropland"

Fig 3: DE-Geb "cropland"

[SarahAlidoost]

@MostafaGomaa93 Another physical condition that should be set correctly is elseif groundWaterDepth > Tot_Depth in calculateGroundWaterDepth. If headBotmLayer=1495.0, this elseif groundWaterDepth > Tot_Depth will be True and indxBotmLayer will be undefined in calculateIndexBottomLayer. This causes the model run to fail. There should be a condition to check whether the values of headBotmLayer are in a valid range.

[SarahAlidoost]

@MostafaGomaa93 This is the output csv file for Fig 2 without setting groundwater_temperature:

• Location=DE-Geb,
• StartTime=2001-01-01T00:00,
• EndTime=2001-01-10T01:00,
• groundwater_elevation_top_aquifer=2000,
• groundwater_head_bottom_layer=1950,
• groundwater_temperature=Nan.

Sim_Temp.csv

[MostafaGomaa93]

The code that I provided plots the entire depth profile (0, 5 m). In your plots here, the y axis is between (-20, 2.5 cm). Please make sure that you can create the exact same plots. So we are looking at the same things. You can also generate an exe file and run the BMI notebook to make sure that plots are the same.

Thanks for the observation. I have updated my figures. I think the provided code may need this slight change
depths = df.iloc[0].index.astype(float) * -1 instead of depths = df.iloc[0].values.astype(float) * -1

After i have updated my figures, i think there are spikes in the first couple of time steps (time step = 1, and time step = 20, spikes are marked with grey circle). At longer timestep (timestep = 60, timestep = 300, spikes are marked with black circle), spikes are very minor or not there at all. So, the spikes problem relaxed with the time.

I also think that it is the same case with the figures that @SarahAlidoost posted

[MostafaGomaa93]

@MostafaGomaa93 Another physical condition that should be set correctly is elseif groundWaterDepth > Tot_Depth in calculateGroundWaterDepth. If headBotmLayer=1495.0, this elseif groundWaterDepth > Tot_Depth will be True and indxBotmLayer will be undefined in calculateIndexBottomLayer. This causes the model run to fail. There should be a condition to check whether the values of headBotmLayer are in a valid range.

Indeed, if headBotmLayer goes below Tot_Depth, then model will fail. The user should have a basic knowledge about the range of the groundwater depth of the site. If the groundwater depth is close to or larger than 5 m, then the Tot_Depth shoud be increased and be larger than 5 m. A solution is provided in this PR (hopefully we can check later after the soil temp issue)

[yijianzeng]

          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:

This depends on the temperature value you set for groundwater temperature. See my other comments here.

Yes, I mean your suggestion by fixing the groundwater temperature = soil temperature wherever GW reaches a certain depth of the soil column. didn't solve the problem yet. So, it looks like there is a bug indeed

Mostafa, as i mentioned the soil moisture content will impact soil temperature simulations. This is about soil physics.

[yijianzeng]

Setting physical parameters:
1- (major)The coupling only sets soil temperature SoilVariables.TT. But, it seems that SoilVariables.T and T are also related and should be set accordingly.

I add the SoilVariables.T and T here. Is that what we need only or we need further edits?

'SoilVariables.T' means the soil temperature calculated from the last time step. 'TT' means soil temperature calculated at this time step.

[yijianzeng]

2- in coupling BoundaryCondition.NBChB = 1;, which is only valid for the sites with landcover 'Croplands'. But, for example ZA-Kru is not a 'Croplands'.

The BoundaryCondition.NBChB = 1 means that the bottom boundary is a head boundary and it is not related to the landcover type. In principle, if we know the head at the bottom of the soil (which is the case when we use MODFLOW), then we use BoundaryCondition.NBChB = 1. I think Yijian may help to answer why for the case of 'Croplands', BoundaryCondition.NBChB = 1

I agree with Mostafa, that the bottom boundary condition is not related to landcover type.

[yijianzeng]

4- (major)The condition in energy/solveEnergyBalanceEquations.m should be checked because it is using ForcingData.Tmin. This might not be valid in the case of coupling.

I removed this condition while testing but found that no difference in the results. Maybe Yijian could clarify what is the purpose of this condition

I think Lianyu added this to allow the numerical model to converge better.

[yijianzeng]

As mentioned before, in stemmusscope (without coupling), BoundaryCondition.NBChB is set to 1 for sites with "cropland" whereas 3 for other sites. However, in groundwater coupling, BoundaryCondition.NBChB is set to 1 for all sites. @yijianzeng Is this a valid assumption? see definition here.

Most of time, we say the bottom boundary condition is free drainage. When it is coupled with groundwater model, we know the soil profile emerged in groundwater will be saturated, then we know the matric head will be '0'.

In short, yes, it is a valid assumption.

[yijianzeng]

@MostafaGomaa93 Another physical condition that should be set correctly is elseif groundWaterDepth > Tot_Depth in calculateGroundWaterDepth. If headBotmLayer=1495.0, this elseif groundWaterDepth > Tot_Depth will be True and indxBotmLayer will be undefined in calculateIndexBottomLayer. This causes the model run to fail. There should be a condition to check whether the values of headBotmLayer are in a valid range.

Indeed, if headBotmLayer goes below Tot_Depth, then model will fail. The user should have a basic knowledge about the range of the groundwater depth of the site. If the groundwater depth is close to or larger than 5 m, then the Tot_Depth shoud be increased and be larger than 5 m. A solution is provided in this PR (hopefully we can check later after the soil temp issue)

I would say, we shall provide a warning message saying something like 'Ground water table depth is larger than the total soil column thickness! Please enlarge the total depth of the soil column.'

[SarahAlidoost]

@BSchilperoort here are the results of the model for other sites, showing that there is an issue when the temperature is Nan and is set to soil temperature. This has been discussed with @MostafaGomaa93. The coupling results look okay for other sites except for some spikes in Fig 0-DE-Seh (croplands), no coupling which is not related to ground water coupling.

DE-Seh (croplands):

ES-LMa (Savannas):

AU-Stp (grassland):