Question about nbedrock and excess-ice initialization in CTSM5.3.026

[dai lz]

Hi CTSM developers,

I have a question about how nbedrock is used when initializing excess ground ice in CTSM5.3.026.

In WaterStateType.F90, excess ice is initialized only when
j >= nexice .and. j < nbedrock
so the layer j = nbedrock is always excluded from excess-ice initialization.

However, my understanding is that nbedrock represents the last standard soil layer associated with the prescribed bedrock depth. For example, the code classifies
1:nbedrock
as LEVGRND_CLASS_STANDARD, while layers below nbedrock are classified as bedrock.

This leads to a few cases that I am having difficulty interpreting physically.

For the default configuration with 20 soil layers, the bedrock depth is set to the bottom interface of layer 20. Therefore, layer 20 is entirely above the bedrock interface and is still a standard soil layer. However, because the excess-ice condition is j < nbedrock, excess ice can only be initialized through layer 19, and layer 20 is excluded.

A similar issue occurs with a shallower prescribed bedrock depth. For example, if zbedrock = 3.2 m, this depth falls within layer 14 (approximately 2.72–3.26 m), so nbedrock = 14. Layer 14 is still classified as a standard layer, but excess ice is only allowed through layer 13. This means that the soil portion between about 2.72 m and 3.2 m cannot contain initialized excess ice.

An even clearer case is when the prescribed bedrock depth is exactly at the bottom interface of layer 14, approximately 3.26 m. In this case, layer 14 is entirely above bedrock, yet it is still excluded from excess-ice initialization because j < nbedrock

This also seems to introduce a discontinuity. For example, if the bedrock depth changes from exactly 3.26 m to slightly deeper than 3.26 m, nbedrock changes from 14 to 15, and layer 14 suddenly becomes eligible for excess ice. Thus, a very small change in bedrock depth can add or remove an entire soil layer of excess ice.

Could you please clarify the intended behavior here?

Specifically:

1. Is the exclusion of the nbedrock layer from excess-ice initialization intentional?
2. What is the intended physical interpretation of nbedrock in this context: the deepest soil layer, or the layer containing the soil-bedrock interface?
3. If the bedrock depth lies within the nbedrock layer, is the whole layer intentionally excluded to avoid assigning excess ice below the actual bedrock interface?
4. If the bedrock depth lies exactly at the bottom interface of the nbedrock layer, should that entire layer physically be allowed to contain excess ice?
5. Could the condition j < nbedrock be an off-by-one issue, or is this a deliberate discretization choice?

I would appreciate any clarification on the intended physical and numerical treatment of the deepest soil layer in the excess-ice initialization.

Thank you!

 

[mvdebolskiy]

Hi, when I was writing that initialization, since there is a lot of special treatment of the nbedrock layer in the soilwater hydrology and since actual vertical profile of excess ice is kind of arbitrary, I decided to exclude the last layer from getting any excess ice to not to get into any issues further down the line. Same way it is not added in the first layer and the excess_ice_coldstart_depth default value is 0.5.
Does this answer the question?

 

[dai lz]

Hi, when I was writing that initialization, since there is a lot of special treatment of the nbedrock layer in the soilwater hydrology and since actual vertical profile of excess ice is kind of arbitrary, I decided to exclude the last layer from getting any excess ice to not to get into any issues further down the line. Same way it is not added in the first layer and the excess_ice_coldstart_depth default value is 0.5.
Does this answer the question?

Yes, thank you! This answers my question and clarifies why the nbedrock layer was intentionally excluded.

I actually tested changing the condition from j &lt; nbedrock to j &lt;= nbedrock in my simulations, so that excess ice can also be initialized in the last soil layer. So far, the simulations have run without any apparent issues.

My bedrock depth is relatively shallow (around 3.2 m), so excluding the nbedrock layer can leave a relatively large portion of the soil profile without excess ice. This was the main reason I was interested in including that layer.

Given that the modified version appears to work normally, do you think there are any particular hydrological or subsidence-related issues that I should watch for when allowing excess ice in the nbedrock layer?

Thank you again for explaining the reasoning behind the original implementation.

 

[mvdebolskiy]

>My bedrock depth is relatively shallow (around 3.2 m), so excluding the nbedrock layer can leave a relatively large portion of the soil profile without excess ice. This was the main reason I was interested in including that layer.

Yes, your will get a bit more excess ice and subsidence (if the whole column thaws). I am not sure about the exact effects on soil hydrology but you should watch out for perched water and possibly nbedrock layer staying completely saturated After all excess ice melts compared to a simulation where there is no excess ice in the nbedrock layer. If you run modified and original up to the point where there is no more excess ice, can you post the timeseries of the difference in soil moisture, temperature and excess totexice_vol here?

 

[dai lz]

Yes, I will run both the original and modified cases and post the comparison here once I have the results. Thank you for the suggestion!