Lower boundary conditions (differences between CLM4.5 and CLM5)

[wvsi3w]

Hello,
Pardon my insufficient knowledge of the Fortran coding but I have to ask the following questions here as I only have this forum to check my inquiries with:

I have found the "Modified CLM4.5 files to change the bottom heat flux and add 24 bedrock layers of thickness 12.5 m at the bottom of the land column (to increase the soil depth from 42.1m to 342.1m)" = this link (Lower boundary conditions improvement in CLM4.5)

In CLM5, the file previously known as iniTimeConst.F90 in CLM4.5, which is used for modifying soil thermal properties, seems to have been renamed to SoilStateInitTimeConstMod.F90

I have located SoilStateInitTimeConstMod.F90 in the my_cesm_sandbox/components/clm/src/biogeophys directory and clm_varpar.F90 in the my_cesm_sandbox/components/clm/src/main directory of my ported CESM2.1.3 model. But, I can not find some of the required modifications explained in the above link, for instance:
1. iniTimeConst.F90 (Layers)
- Variable thick_equal changed form 0.2 to 12.5 (m)
- Original code for zsoi, zisoi and dzsoi changed to code in Appendix A

What is this thick_equal variable name in CLM5?
What is that Appendix A referring to?

Also, I noticed that the structure of those modified files in CLM4.5 (the link above) is almost different from the current equivalent files in CLM5. How should I modify them (in CLM5) without messing up the whole code?

I also need to do the rest of the modifications (regarding adding bottom heat flux) to CLM5, so your help would be tremendous.

P.S. I contacted the people who did this on CLM4.5 and unfortunately, they are unavailable or maybe unwilling to help.

 

[oleson]

Sorry, I'm not familiar with that work, so can't be of much help here. I will ask our permafrost expert here if they are aware of this work.
I assume Appendix A maybe refers to something in the author's paper? We don't have any Appendix A in our technical note.
thick_equal shows up in main/initVerticalMod.F90 in cesm2.1.3. It looks like if the soil layer structure is specified as 23SL_3.5m, then the soil layer structure "starts with standard exponential and then has several evenly spaced layers, then finishes off exponential. this allows the upper soil to behave as standard, but then continues with higher resolution to a deeper depth...". So thick_equal must refer to the thickness of the evenly spaced layers.

 

[oleson]

@slevis , perhaps you can provide some information on Monday.

 

[slevis]

At Dave L's recommendation, I will mention this capability as a possible path forward:
@wvs3iw consider using the soil_layerstruct_userdefined namelist option described in bld/namelist_files/namelist_definition_ctsm.xml. To use, enter a vector of soil thicknesses (dzsoi) in the namelist to define the soil structure and set soil_layerstruct_userdefined_nlevsoi accordingly. Please read the instructions carefully (in the same xml file) for details in the usage of these options.

 

[wvsi3w]

At Dave L's recommendation, I will mention this capability as a possible path forward:
@wvs3iw consider using the soil_layerstruct_userdefined namelist option described in bld/namelist_files/namelist_definition_ctsm.xml. To use, enter a vector of soil thicknesses (dzsoi) in the namelist to define the soil structure and set soil_layerstruct_userdefined_nlevsoi accordingly. Please read the instructions carefully (in the same xml file) for details in the usage of these options.

Thanks for the response,
I will read about that and try it and let you know.

The GOOD news is that I found another person who has worked on changing the lower boundary conditions of CTSM (not CESM) a few years ago, and in his work, I found the following:

He changed three main F90 scripts (clm_varpar, initVerticalMod, and namelist_definition_ctsm) which are attached here. I have ported and configured CESM and not CTSM on my system because I have to consider other components (like CAM) for my project, but for now, I have to use the land component. So, I checked the differences between the three F90 codes in the modified CTSM with what I have on CESM-CLM5 and the three "diff" files that are attached are the results.

I also checked his modifications on his version of CTSM with the newest CTSM available on github and except for some small differences most of them were similar.

The problem is the huge difference I found between these _DG modifications on CTSM and my CLM5. I know that I can do the same modifications he did on my version of CLM5 but I need to be sure that this is the right way to do it.

I mean, for those _DG modifications (_DG: Deep Ground) I should only look for the right place to define them in the three F90 codes that I have? and I don't need to be worried about other differences (not DG differences)?

Or, do you think it is better to port and configure the updated version of CTSM on my system and try those DG modifications with it?
Or do you think that the differences shown in those "diff" files are negligible for my purpose and I can continue these modifications with the CESM that I spend a lot porting and configuring?

Thanks for your time,

P.S. one last question: that namelist_definition file in CESM is named clm4_5 and I couldn't find any clm5 in that path. Is this ok?

 

[slevis]

I do not have a straightforward answer to your questions. You will have to make a decision and try it. If it doesn't go well, you may reconsider and try a different approach.

Regarding the namelist definition file, if the model has been running successfully, then it's correct.

 

[wvsi3w]

I do not have a straightforward answer to your questions. You will have to make a decision and try it. If it doesn't go well, you may reconsider and try a different approach.

Regarding the namelist definition file, if the model has been running successfully, then it's correct.

Dear @slevis thanks for the feedback.
I contacted the person who modified the code some years ago and he suggested to "definitely stick to the CESM version you have running on your system. Implementing the _DG code in your CESM version should be easy"

In his version, he managed to change the depth to 198m.
for instance, this is what he has for _DG in his clm-varpar file:

else if ( soil_layerstruct == '20SL_8.5m_DG' ) then
nlevsoi = 20
nlevgrnd = nlevsoi+15

And this is what he has in his initVerticalMod file:

else if ( soil_layerstruct == '20SL_8.5m_DG' ) then
do j = 1,4
dzsoi(j)= j*0.02_r8 ! linear increase in layer thickness of 2cm each layer
enddo
do j = 5,13
dzsoi(j)= dzsoi(4)+(j-4)*0.04_r8 ! linear increase in layer thickness of 4cm each layer
enddo
do j = 14,nlevsoi
dzsoi(j)= dzsoi(13)+(j-13)*0.10_r8 ! linear increase in layer thickness of 10cm each layer
enddo
do j = nlevsoi+1,nlevsoi+5 !bedrock layers
dzsoi(j)= dzsoi(nlevsoi)+(((j-nlevsoi)*25._r8)**1.5_r8)/100._r8 ! bedrock layers
enddo
do j = nlevsoi+6,nlevgrnd !bedrock layers
dzsoi(j)= dzsoi(nlevsoi)+(((5)*25._r8)**1.5_r8)/100._r8 ! bedrock layers
enddo

zisoi(0) = 0._r8
do j = 1,nlevgrnd
zisoi(j)= sum(dzsoi(1:j))
enddo

do j = 1, nlevgrnd
zsoi(j) = 0.5*(zisoi(j-1) + zisoi(j))
enddo

The bold lines are the differences between DG code and regular code.

For example here in the default soil_layerstruct == "20SL_8.5m" , we have:

else if ( soil_layerstruct == '20SL_8.5m' ) then
do j = 1,4
dzsoi(j)= j*0.02_r8 ! linear increase in layer thickness of 2cm each layer
enddo
do j = 5,13
dzsoi(j)= dzsoi(4)+(j-4)*0.04_r8 ! linear increase in layer thickness of 4cm each layer
enddo
do j = 14,nlevsoi
dzsoi(j)= dzsoi(13)+(j-13)*0.10_r8 ! linear increase in layer thickness of 10cm each layer
enddo
do j = nlevsoi+1,nlevgrnd !bedrock layers
dzsoi(j)= dzsoi(nlevsoi)+(((j-nlevsoi)*25._r8)**1.5_r8)/100._r8 ! bedrock layers
enddo

zisoi(0) = 0._r8
do j = 1,nlevgrnd
zisoi(j)= sum(dzsoi(1:j))
enddo

do j = 1, nlevgrnd
zsoi(j) = 0.5*(zisoi(j-1) + zisoi(j))
enddo

The question I have is how do you calculate the 20 soil layers. because when I do a simple calculation in excel, it is lower than that for the 20th layer, and if I continue the calculation for the next 5 layers (until nlevgrnd), I will see more than 8.5m of depth (15m). I attached the PNG of this calculation.

If I understand these simple calculations I would be sure of the198m modification too and I would probably be able to increase the depth to more than that, say 500m.

 

[wvsi3w]

Dear @slevis thanks for the feedback.
I contacted the person who modified the code some years ago and he suggested to "definitely stick to the CESM version you have running on your system. Implementing the _DG code in your CESM version should be easy"

In his version, he managed to change the depth to 198m.
for instance, this is what he has for _DG in his clm-varpar file:


And this is what he has in his initVerticalMod file:

The bold lines are the differences between DG code and regular code.

For example here in the default soil_layerstruct == "20SL_8.5m" , we have:



The question I have is how do you calculate the 20 soil layers. because when I do a simple calculation in excel, it is lower than that for the 20th layer, and if I continue the calculation for the next 5 layers (until nlevgrnd), I will see more than 8.5m of depth (15m). I attached the PNG of this calculation.

If I understand these simple calculations I would be sure of the198m modification too and I would probably be able to increase the depth to more than that, say 500m.

I checked the Fortran code I realized something (maybe I am wrong), it says "linear increase in layer thickness of Xcm each layer". This means that at the end we should sum up all the increments of the layers (as I depicted in the png attached here) to reach the actual depth.

# Python script to calculate the total depth based on layer thickness formulas

# Initialize variables
nlevsoi = 20  # Number of soil layers
dzsoi = [0.0] * 35  # Initialize a list for 35 layers (20 soil + 15 bedrock)

# Calculate thickness for layers 1 to 20
for j in range(1, 21):
    if j <= 4:
        dzsoi[j-1] = j * 0.02  # Thickness increases linearly by 2cm for each layer
    elif j <= 13:
        dzsoi[j-1] = 0.08 + (j - 4) * 0.04  # Thickness increases from the 4th layer by 4cm for each layer
    else:
        dzsoi[j-1] = 0.44 + (j - 13) * 0.10  # Thickness increases from the 13th layer by 10cm for each layer

# Calculate thickness for bedrock layers 21 to 35
for j in range(21, 36):
    if j <= 25:
        dzsoi[j-1] = dzsoi[19] + (((j - 20) * 25.0)**1.5) / 100.0  # Complex formula for layers 21 to 25
    else:
        dzsoi[j-1] = dzsoi[19] + ((5 * 25.0)**1.5) / 100.0  # Fixed thickness for layers 26 to 35

# Calculate total depth by summing up the thicknesses
total_depth = sum(dzsoi)

print(total_depth)  # Total depth in meters

 

[wvsi3w]

Hello dear @slevis and @oleson ,

I've extended the depth of the lower boundary in the CLM5 model from 50m to ~194m by modifying the layer structure based on a geometric progression as detailed in my previous question.



For layers 1st to 20th:
Δz_soil[j] = {
j × 0.02 if j ≤ 4
0.08 + (j - 4) × 0.04 if 5 ≤ j ≤ 13
0.44 + (j - 13) × 0.10 if j > 13
}

For layers 21st to 25th:
Δz_soil[j - 1] = Δz_soil[19] + ((j - 20) × 25.0^1.5) / 100.0

For layers 26th to 35th:
Δz_soil[j - 1] = Δz_soil[19] + ((5 × 25.0^1.5) / 100.0)


However, after running simulations, I'm observing that the active layer thickness (ALT) seems to reach up to the total depth of the model, which is unrealistic.
Default Simulation ALT Max value: 41.99 m
Modified Simulation ALT Max value: 193.15m

Could you provide insights into why the ALT might be erroneously reflecting the depth of the entire model? In extending the model depth, I have added 10 more bedrock layers following the layer thickness progression equations provided. I want to confirm if my approach is correct and aligns with the recommended practices for adjusting model depth in CLM5. Are there any additional considerations or typical mistakes to be aware of when implementing such modifications?

I didn't want to increase the ALT, I wanted to increase the bottom boundary condition of CLM5, I saw the results of ALT and plotted the difference between modified and default simulations and realized this significant increase in ALT for 40 years of test runs which was really questionable and odd. I think I did something wrong in the modification of depth. Please help.

 

[slevis]

It may help to look at the simulated soil temperatures and heat fluxes to determine whether it makes sense that the simulated ALT reaches the lower boundary of the model.

I do not remember, maybe you know: Does the model include a geothermal heat flux from the lower boundary? If not, this may be one explanation for why the model would behave one way and nature another way.

 

[wvsi3w]

Below is the table from two simulation I did (one with default depth (50m) and another with modified depth (194m)):

Variable​	Simulation​	Mean​	Median​	Standard Deviation​
ALT	Default	30.149​	41.998​	18.701​
	Modified	137.843​	193.153​	86.59​
CH4PROD	Default	9.69E-08​	0​	4.28E-07​
	Modified	9.69E-08​	0​	4.31E-07​
HEAT_CONTENT1	Default	-324,957,280​	683,516,544​	2,991,275,520​
	Modified	-2,076,964,608​	1,253,199,232​	9,961,728,000​
SOILC_CHANGE	Default	3.16E-06​	7.38E-11​	7.17E-06​
	Modified	3.16E-06​	7.37E-11​	7.17E-06​
STORVEGC	Default	200.292​	99.377​	249.385​
	Modified	200.282​	99.356​	249.384​
TSL	Default	283.189​	285.613​	17.481​
	Modified	283.185​	285.613​	17.485​
TSOI_10CM	Default	270.563​	274.332​	27.469​
	Modified	270.556​	274.332​	27.475​

As you can see the Heat content of the modified simulation (194m depth) has increased significantly and ALT also increased, but for 40 years of simulation with similar condition (except for the depth) this increase in ALT is really incorrect.

I also plotted the whole 40 years and by looking at the plots the Temperature of soil has decreased slightly. This part makes sense because by increasing the depth the propagated heat wont bounce back to the surface, BUT still I am confused by the ALT values.

I believe the way I increased the depth could be wrong and my question here is that do you know which part in the code I should change so that I only increase the bedrock and not ALT?

About the heat flux at the bottom, almost all the LSMs consider 0 heat flux from the bottom boundary and CLM is one of them. And as far as I read, the impact of heat flux at the bottom for these kinds of simulations is insignificant as the surface temperature is the main controller of surface layer energy.

 

[slevis]

I believe the way I increased the depth could be wrong and my question here is that do you know which part in the code I should change so that I only increase the bedrock and not ALT?

Correct me if I'm wrong, but I think ALT is a simulated quantity. If so, then you probably need to find the calculation and override it if you don't want ALT to change.

 

[wvsi3w]

Correct me if I'm wrong, but I think ALT is a simulated quantity. If so, then you probably need to find the calculation and override it if you don't want ALT to change.

Thank you @slevis for your response. The problem is we don't need ALT to increase as the bottom boundary does and ALT should be simulated normally because in our study the permafrost and ALT have a strong correlation and by hindering ALT from change we are sort of ruining the whole purpose of the study.

Where permafrost melts, ALT increases.

And by increasing the bottom boundary of CLM5 we aim to deepen the bedrock to a certain depth (say 500m) so that the signal (energy or T from the surface) propagates downward effectively (close to reality) without bouncing back from bottom. When the depth is around 50m this bounce happens and the simulation shows a lot of permafrost thaw but when the bottom boundary is more (in our case 200, 300 or 500m) this bounce of heat doesn't happen and the permafrost thaws closer to reality and ALT changes the same way.

The bottom heat flux also plays a role here but not as much as other stuff.

So I think the modification that I did to the CLM5 (increasing the depth to ~200m) is kind of wrong and I think there must be a line of code or some other modification to increase the bedrock layers and not the soil itself.

 

[slevis]

If you are concerned about an error in your modification, then you will need to troubleshoot by exploring the affected sections of code, possibly adding write statements to become familiar with the contents of different variables and so on.

 

[wvsi3w]

If you are concerned about an error in your modification, then you will need to troubleshoot by exploring the affected sections of code, possibly adding write statements to become familiar with the contents of different variables and so on.

Hello again @slevis

I checked the code and realized the soil layer structure that I modified is only increasing the nlevsoi variable (Maybe I am wrong and we need to turn on something related to bedrock in the simulation?????):

else if ( soil_layerstruct == '20SL_8.5m_DG' ) then
       do j = 1,4
          dzsoi(j)= j*0.02_r8          ! linear increase in layer thickness of 2cm each layer
       enddo
       do j = 5,13
          dzsoi(j)= dzsoi(4)+(j-4)*0.04_r8      ! linear increase in layer thickness of 4cm each layer
       enddo
       do j = 14,nlevsoi       
          dzsoi(j)= dzsoi(13)+(j-13)*0.10_r8     ! linear increase in layer thickness of 10cm each layer
       enddo
       do j = nlevsoi+1,nlevsoi+5 !bedrock layers
          dzsoi(j)= dzsoi(nlevsoi)+(((j-nlevsoi)*25._r8)**1.5_r8)/100._r8  ! bedrock layers
       enddo
       do j = nlevsoi+6,nlevgrnd !bedrock layers
          dzsoi(j)= dzsoi(nlevsoi)+(((5)*25._r8)**1.5_r8)/100._r8  ! bedrock layers
       enddo
      
       zisoi(0) = 0._r8
       do j = 1,nlevgrnd
          zisoi(j)= sum(dzsoi(1:j))
       enddo
      
       do j = 1, nlevgrnd
          zsoi(j) = 0.5*(zisoi(j-1) + zisoi(j))
       enddo


    end if

However, I noticed that there is another soil layer structure (23SL_3.5m) which is specifically for permafrost dynamic. I believe I should modify this one (?) but there is one question, in the name of soil layer it says 3.5m, which I guess is related to the total depth of this soil layer, right? But when I calculate the total depth it reaches to 116m !!!

do j = 1, toplev_equalspace
          zsoi(j) = scalez*(exp(0.5_r8*(j-0.5_r8))-1._r8)    !node depths
       enddo

       do j = toplev_equalspace+1,toplev_equalspace + nlev_equalspace
          zsoi(j) = zsoi(j-1) + thick_equal
       enddo

       do j = toplev_equalspace + nlev_equalspace +1, nlevgrnd
          zsoi(j) = scalez*(exp(0.5_r8*((j - nlev_equalspace)-0.5_r8))-1._r8) + nlev_equalspace * thick_equal
       enddo

       dzsoi(1) = 0.5_r8*(zsoi(1)+zsoi(2))             !thickness b/n two interfaces
       do j = 2,nlevgrnd-1
          dzsoi(j)= 0.5_r8*(zsoi(j+1)-zsoi(j-1))
       enddo
       dzsoi(nlevgrnd) = zsoi(nlevgrnd)-zsoi(nlevgrnd-1)

       zisoi(0) = 0._r8
       do j = 1, nlevgrnd-1
       zisoi(j) = 0.5_r8*(zsoi(j)+zsoi(j+1))         !interface depths
       enddo
       zisoi(nlevgrnd) = zsoi(nlevgrnd) + 0.5_r8*dzsoi(nlevgrnd)

and this is the python code for it:

import numpy as np

# Parameters
scalez = 1.0  # Scaling factor for exponential layers, adjust as needed
thick_equal = 0.1  # Thickness for each evenly spaced layer, adjust as needed
toplev_equalspace = 5  # Number of initial exponential layers, adjust as needed
nlev_equalspace = 15  # Number of evenly spaced layers, adjust as needed
nlevgrnd = toplev_equalspace + nlev_equalspace + toplev_equalspace  # Total number of layers

# Initialize arrays
zsoi = np.zeros(nlevgrnd)
dzsoi = np.zeros(nlevgrnd)
zisoi = np.zeros(nlevgrnd)

# Calculate zsoi for exponential and evenly spaced layers
for j in range(1, toplev_equalspace + 1):
    zsoi[j-1] = scalez * (np.exp(0.5 * (j - 0.5)) - 1)

for j in range(toplev_equalspace + 1, toplev_equalspace + nlev_equalspace + 1):
    zsoi[j-1] = zsoi[j-2] + thick_equal

for j in range(toplev_equalspace + nlev_equalspace + 1, nlevgrnd + 1):
    zsoi[j-1] = scalez * (np.exp(0.5 * ((j - nlev_equalspace) - 0.5)) - 1) + nlev_equalspace * thick_equal

# Calculate dzsoi for each layer
dzsoi[0] = 0.5 * (zsoi[0] + zsoi[1])
for j in range(2, nlevgrnd):
    dzsoi[j-1] = 0.5 * (zsoi[j] - zsoi[j-2])
dzsoi[nlevgrnd-1] = zsoi[nlevgrnd-1] - zsoi[nlevgrnd-2]

# Calculate zisoi for the interfaces
zisoi[0] = 0
for j in range(1, nlevgrnd):
    zisoi[j] = 0.5 * (zsoi[j-1] + zsoi[j])
zisoi[nlevgrnd-1] += 0.5 * dzsoi[nlevgrnd-1]

# Output (optional)
print("zsoi:", zsoi)
print("dzsoi:", dzsoi)
print("zisoi:", zisoi)

And when you run the python code this is the results:
zsoi values: [ 0.28 , 1.12 , 2.49 , 4.75 , 8.49 , 8.59 , 8.69 , 8.79 , 8.89 , 8.99 , 9.09 , 9.19 , 9.29 , 9.39 , 9.49 , 9.59 , 9.69 , 9.79 , 9.89 , 9.99 , 16.14 , 26.29 , 43.02 , 70.61 , 116.08 ]

dzsoi values: [ 0.70 , 1.10 , 1.82 , 3.00 , 1.92 , 0.10 , 0.10 , 0.10 , 0.10 , 0.10 , 0.10 , 0.10 , 0.10 , 0.10 , 0.10 , 0.10 , 0.10 , 0.10 , 0.10 , 3.13 , 8.15 , 13.44 , 22.16 , 36.53 , 45.48 ]

zisoi values: [ 0.00 , 0.70 , 1.80 , 3.62 , 6.62 , 8.54 , 8.64 , 8.74 , 8.84 , 8.94 , 9.04 , 9.14 , 9.24 , 9.34 , 9.44 , 9.54 , 9.64 , 9.74 , 9.84 , 9.94 , 13.07 , 21.22 , 34.66 , 56.81 , 116.08 ]

P.S. I have run the model with this soil layer structure and its ALT results are in reasonable range.

 

[Zh Chen]

Thanks for the response,
I will read about that and try it and let you know.

The GOOD news is that I found another person who has worked on changing the lower boundary conditions of CTSM (not CESM) a few years ago, and in his work, I found the following:

He changed three main F90 scripts (clm_varpar, initVerticalMod, and namelist_definition_ctsm) which are attached here. I have ported and configured CESM and not CTSM on my system because I have to consider other components (like CAM) for my project, but for now, I have to use the land component. So, I checked the differences between the three F90 codes in the modified CTSM with what I have on CESM-CLM5 and the three "diff" files that are attached are the results.

I also checked his modifications on his version of CTSM with the newest CTSM available on github and except for some small differences most of them were similar.

The problem is the huge difference I found between these _DG modifications on CTSM and my CLM5. I know that I can do the same modifications he did on my version of CLM5 but I need to be sure that this is the right way to do it.

I mean, for those _DG modifications (_DG: Deep Ground) I should only look for the right place to define them in the three F90 codes that I have? and I don't need to be worried about other differences (not DG differences)?

Or, do you think it is better to port and configure the updated version of CTSM on my system and try those DG modifications with it?
Or do you think that the differences shown in those "diff" files are negligible for my purpose and I can continue these modifications with the CESM that I spend a lot porting and configuring?

Thanks for your time,

P.S. one last question: that namelist_definition file in CESM is named clm4_5 and I couldn't find any clm5 in that path. Is this ok?

Hi, wvsi3w

I'm interested in the article you mentioned. If you don't mind, would you be willing to share it with me? I would greatly appreciate it. Thank you!

 

[wvsi3w]

Hi, wvsi3w

I'm interested in the article you mentioned. If you don't mind, would you be willing to share it with me? I would greatly appreciate it. Thank you!

Hello @Zh Chen

Unfortunately, this person did not publish the related article and I used my connections to go through his code modifications (contacted him personally). But if you are interested in similar papers this one (Lower boundary conditions in land surface models – effects on the permafrost and the carbon pools: a case study with CLM4.5) is much like what I have in mind but for previous version of the model.