Open SarahAlidoost opened 2 months ago
:red_circle: This issue should be fixed before releasing stemmus_scope, see details in https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/239
@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
?
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.
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.
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.
@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.
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.
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
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
@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?
I also noticed the implementation of
SoilVariable.TT
andTT
is not ideal inSTEMMUS_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.
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.
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
@SarahAlidoost @yijianzeng
After our meeting today, I test again by replacing the
tempBotm
in this line bytempBotm = SoilVariables.TT(indxBotm)
which is = 22.96 for the ZA-Kru and the problem is still there.
Hi Mostafa, please try the following:
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.
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
another test with changing the groundwater depth over time (~1m drop per 8 days)
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.
@MostafaGomaa93 here are the codes that we discussed at our meeting.
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 ofP_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
, anddryair
modules must be set manually to certain values (as we did already).3- (major)The soil temperature
SoilVariables.TT
should be set inupdate
loop here. But, due to current codes, the variableTT
should be set too.Okay, I add them here
4- (minor issue) RHS should be defined in module
+energy
asRHS = InitialValues.RHS;
inassembleCoefficientMatrices.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 thatSoilVariables.T
andT
are also related and should be set accordingly.I add the
SoilVariables.T
andT
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 useBoundaryCondition.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 asheadBotmLayer - topLevel + soilThick(indxBotmLayer_R)
. First of all the units don't match. Second, this equation might return zero.If the
BoundaryCondition.NBChB = 1
, then theRHS
is a head and has units of cm. If theBoundaryCondition.NBChB = 2
, then theRHS
is a flux and has units of cm/sec. In the case of groundwater coupling,BoundaryCondition.NBChB = 1
andRHS
is a head and has units of cm and theheadBotmLayer
,topLevel
,soilThick(indxBotmLayer_R)
are all in cm. and the value ofRHS(indexbotm)
is indeed close to zero4- (major)The condition in
energy/solveEnergyBalanceEquations.m
should be checked because it is usingForcingData.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 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
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 |
Another two examples at groundwater depth = 250
and (average condition) is not used
groundwater temperature = 23 |
groundwater temperature = 17 |
---|---|
@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.
@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"
@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.
@MostafaGomaa93 This is the output csv file for Fig 2 without setting groundwater_temperature:
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 Another physical condition that should be set correctly is
elseif groundWaterDepth > Tot_Depth
in calculateGroundWaterDepth. IfheadBotmLayer=1495.0
, thiselseif groundWaterDepth > Tot_Depth
will be True andindxBotmLayer
will be undefined in calculateIndexBottomLayer. This causes the model run to fail. There should be a condition to check whether the values ofheadBotmLayer
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'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.
Setting physical parameters: 1- (major)The coupling only sets soil temperature
SoilVariables.TT
. But, it seems thatSoilVariables.T
andT
are also related and should be set accordingly.I add the
SoilVariables.T
andT
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.
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 useBoundaryCondition.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.
4- (major)The condition in
energy/solveEnergyBalanceEquations.m
should be checked because it is usingForcingData.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.
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.
@MostafaGomaa93 Another physical condition that should be set correctly is
elseif groundWaterDepth > Tot_Depth
in calculateGroundWaterDepth. IfheadBotmLayer=1495.0
, thiselseif groundWaterDepth > Tot_Depth
will be True andindxBotmLayer
will be undefined in calculateIndexBottomLayer. This causes the model run to fail. There should be a condition to check whether the values ofheadBotmLayer
are in a valid range.Indeed, if
headBotmLayer
goes belowTot_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 theTot_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.'
@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):
_Originally posted by @BSchilperoort in https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/239#issuecomment-2242138630_