Surface soil thermal conductivity over snow-free areas (bug)
Summary
This issue was presented to the CLASSIC biweekly meeting agenda for May 8th, 2024 ⇾ #117 (closed)
See related talk: https://docs.google.com/presentation/d/1q9-5xaOgYiTXiobNIBiHBQYeAA44pbCmJRBcgpKBFkM/edit?usp=sharing
In summary, the soil surface thermal conductivity during periods of fractional snow cover (i.e., when snow depth is less than 10 cm for the current model version) was averaged between the first soil layer and the snowpack thermal conductivities, over all the subareas considered (i.e., with or without snow). This issue is due to an error of argument passed to the call of the soilHeatFluxPrep subroutine for the two non-snow-covered sub-areas (BARE GROUND and CANOPY OVER BARE GROUND) in energyBudgetDriver.f90
Details
The soilHeatFluxPrep.f90 subroutine includes a condition to compute the surface soil thermal conductivity over snow-covered and snow-free areas based on the snow depth:
if (ZSNOW(I) > 0.0) then
TCZERO = 1.0 / (0.5 / TCSNOW(I) + 0.5 / TCTOP(I,1))
else
TCZERO = TCTOP(I,1)
end if
However, when the snow depth is lower than 10 cm on average over the whole grid cell (i.e., for ZSNOW(I) < SNOLIM(I)
; see code below from radiationDriver.f90), the snow cover becomes fractional, and the snow depth is reset to 10 cm over the snow fraction area (SNOLIM is the limiting snow depth below which coverage is < 100%. It is defined as ZSNLROT = 0.10
in src/dynamicTiling.f90 and src/modelStateDrivers.f90.).
do I = IL1,IL2
if (SNO(I) > 0.0) then
ZSNOW(I) = SNO(I) / RHOSNO(I)
if (ZSNOW(I) >= (SNOLIM(I) - 0.00001)) then
FSNOW(I) = 1.0
else
FSNOW(I) = ZSNOW(I) / SNOLIM(I)
ZSNOW(I) = SNOLIM(I)
WSNOW(I) = WSNOW(I) / FSNOW(I)
end if
else
ZSNOW(I) = 0.0
FSNOW(I) = 0.0
end if
Therefore, when snow cover is fractional, ZSNOW(I) is always > 0 (as it is set for the fraction covered with snow), and the condition in soilHeatFluxPrep.f90 is not effective for snow-free areas. In turn, the soil surface conductivity is underestimated over snow-free areas (as snow thermal conductivity is usually lower than the soil), and the conductive flux too, which can lead to higher soil temperature at the beginning of the snow season when the atmosphere is cooler than the soil (as it won't cool "fast" enough) and the opposite in spring when the snow melt (i.e., the soil won't warm enough over the snow-free subareas).
This can be illustrated by printing TCZERO in the model:
as we can see, the soil surface thermal conductivity over the CANOPY OVER BARE GROUND (here as an example), i.e., snow-free subareas get a snow depth of 10 cm and therefore TCZERO has the same value as for the snow-covered subareas (as explained above).
Solution
An easy way to solve this issue is to pass ZERO
instead of ZSNOW
to both soilHeatFluxPrep calls over the snow-free subareas in energyBudgetDriver.f90 (as it is done in the energBalVegSolve and energBalNoVegSolve calls over the snow-free subareas. An alternative solution would be to use the flags ISNOW or IWATER, and modify the condition in soilHeatFluxPrep. The first solution will be proposed as a merge request to solve this issue.
This solution allow getting correct surface soil thermal conductivities over snow-free subareas:
Note that this issue only happens during fractional snow cover timings (but persisting soil temperature biases may occur all year round). Otherwise, when there is no snow, the snow depth is 0, and when it is > 10 cm, it covers the whole grid cell, which does not trigger this bug.
Implications
Here are 3 examples demonstrating the implication of this bug.
1D simulation at an Arctic site (Bylot Island)
This first example shows that depending on the years and how long the period of fractional snow cover lasts, the impact on soil temperatures can be negligible or important. On the second year (when the fractional snow lasts from September to October) the soil temperature decreases much faster with the bud corrected (orange line) leading to soil temperature differences of more than 5°C which attenuate during the snow season.
Spatial simulations over Canada
The figures below show similar effects but in a wider range. In November, we observed cooler temperatures for the simulation with the bug correction (Zero) compared to the control run, especially in the Arctic areas (where the snowpacks might still be thin and with large soil-atmosphere temperature gradients). These biases persist during winter and more or less cancel each other out in May (because the soil will warm faster), and warmer temperatures are simulated around the Rockies areas. The last figure show that those biases persist all year round beyond 3m soil depth.
More plots can be seen in the talk (link above in the summary).
Boreholes ⇾ permafrost
See the paper Melton et al. (2019) for more information. The bug correction induces slight degradation in the model skills.
For more information on those examples, please contact me (Mickaël Lalande), Libo Wang, and/or Joe Melton.
Discussion
The question remains if even over the snow-covered subareas, it is justified to use the weighted harmonic average between the first soil layer and the snowpack, as the snow is an explicit layer on top with its own snow thermal conductivity, therefore at the interface between the snow and the soil we should get:
for small dz. Further investigations will be carried out on that aspect to know if it is justified or not. Another question remains about the computation of the temperature at the bottom of the snowpack TSNBOT (https://gitlab.com/cccma/classic/-/blob/master/src/snowTempUpdate.f90?ref_type=heads#L139), but that will likely be covered in another issue if needed.
Miscellaneous
Model subareas (in each grid cell)
Illustration of the four subareas in CLASSIC. Each subarea evolves dynamically at each time step.
Snow cover fraction parameterization
Snow cover fraction parameterization in CLASSIC.
Example illustrating how the snow depth and liquid water content within the snowpack are treated when the average snow depth over the grid cell is below the SNOLIM threshold of 10 cm. In the above example, the 5 cm grid cell average snow depth (left) leads to a snow cover fraction of 0.5, and the snow depth is reset to 10 cm (right). The snow density stays unchanged, i.e., the snow water equivalent will be conserved (like the liquid water content in the snow). For example, if the snow density is 100 kg m-3, the snow water equivalent is 0. 05 × 100 = 5 kg m-2 (over the whole grid cell, i.e. fraction of 1; left), and it will be 10 kg m-2 of the 0.5 grid cell fraction (right).