Snow thermal conductivity changes: Sturm1997 and Jordan1991 (#101)

kvrigor · web-flow · commit 98054f5eb6df · 2026-03-25T15:13:59.000+01:00
This PR implements @AdrienDams snow thermal conductivity fixes in CLM5. > Damseaux, A., Matthes, H., Dutch, V. R., Wake, L., and Rutter, N.: **Impact of snow thermal conductivity schemes on pan-Arctic permafrost dynamics in the Community Land Model version 5.0**, The Cryosphere, 19, 1539–1558, https://doi.org/10.5194/tc-19-1539-2025, 2025. ## New namelist parameters in `lnd_in` These fixes allow changing the snow thermal conductivity per land unit type. New namelist switches must be added in the `lnd_in` with defaults set to [CLM6 defaults](https://github.com/ESCOMP/CTSM/blob/9f5c3602f9efea2ade07a8f5dcf9eb2b0425d540/bld/namelist_files/namelist_defaults_ctsm.xml#L576-L583). ```fortran &clm_inparm snow_thermal_cond_method = 'Sturm1997' snow_thermal_cond_glc_method = 'Sturm1997' snow_thermal_cond_lake_method = 'Jordan1991' ``` For reference the namelist parameters reproducing the **old defaults** ```fortran &clm_inparm snow_thermal_cond_method = 'Jordan1991' snow_thermal_cond_glc_method = 'Jordan1991' snow_thermal_cond_lake_method = 'Jordan1991' ``` ## Code changes Changes were based from these upstream PRs: - ESCOMP/CTSM#2148 - ESCOMP/CTSM#3072
diff --git a/src/clm5/biogeophys/LakeFluxesMod.F90 b/src/clm5/biogeophys/LakeFluxesMod.F90
@@ -8,6 +8,8 @@ module LakeFluxesMod
   ! !USES
   use shr_kind_mod         , only : r8 => shr_kind_r8
   use shr_log_mod          , only : errMsg => shr_log_errMsg
+  use abortutils           , only : endrun
+  use clm_varctl           , only : iulog
   use decompMod            , only : bounds_type
   use atm2lndType          , only : atm2lnd_type
   use EnergyFluxType       , only : energyflux_type
@@ -29,6 +31,9 @@ module LakeFluxesMod
   !
   ! !PUBLIC MEMBER FUNCTIONS:
   public :: LakeFluxes
+
+  character(len=*), parameter, private :: sourcefile = &
+  __FILE__
   !-----------------------------------------------------------------------
 
 contains
@@ -50,6 +55,7 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
     use clm_varcon          , only : hvap, hsub, hfus, cpair, cpliq, tkwat, tkice, tkair
     use clm_varcon          , only : sb, vkc, grav, denh2o, tfrz, spval, zsno
     use clm_varctl          , only : use_lch4
+    use clm_varctl          , only : snow_thermal_cond_lake_method
     use LakeCon             , only : betavis, z0frzlake, tdmax, emg_lake
     use LakeCon             , only : lake_use_old_fcrit_minz0
     use LakeCon             , only : minz0lake, cur0, cus, curm, fcrit
@@ -403,7 +409,20 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep,
             else
                !Need to calculate thermal conductivity of the top snow layer
                bw = (h2osoi_ice(c,jtop(c))+h2osoi_liq(c,jtop(c)))/dz(c,jtop(c))
-               tksur(c) = tkair + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkair)
+               select case (snow_thermal_cond_lake_method)
+               case ('Jordan1991')
+                  tksur(c) = tkair + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkair)
+               case ('Sturm1997')
+                  if (bw <= 156) then
+                     tksur(c) = 0.023 + 0.234*(bw/1000)
+                  else
+                     tksur(c) = 0.138 - 1.01*(bw/1000) +(3.233*((bw/1000)*(bw/1000)))
+                  end if
+               case default
+                  write(iulog,*) ' ERROR: unknown snow_thermal_cond_lake_method value: ', snow_thermal_cond_lake_method
+                  call endrun(msg=errMsg(sourcefile, __LINE__))
+               end select
+
                tsur(c) = t_soisno(c,jtop(c))
             end if
 
diff --git a/src/clm5/biogeophys/LakeTemperatureMod.F90 b/src/clm5/biogeophys/LakeTemperatureMod.F90
@@ -9,6 +9,7 @@ module LakeTemperatureMod
   !
   ! !USES
   use shr_kind_mod      , only : r8 => shr_kind_r8
+  use abortutils        , only : endrun
   use shr_log_mod       , only : errMsg => shr_log_errMsg
   use decompMod         , only : bounds_type
   use ch4Mod            , only : ch4_type
@@ -19,6 +20,7 @@ module LakeTemperatureMod
   use TemperatureType   , only : temperature_type
   use WaterfluxType     , only : waterflux_type
   use WaterstateType    , only : waterstate_type
+  use clm_varctl        , only : iulog
   use ColumnType        , only : col                
   use PatchType         , only : patch                
   !    
@@ -117,7 +119,7 @@ subroutine LakeTemperature(bounds, num_lakec, filter_lakec, num_lakep, filter_la
     use clm_time_manager   , only : get_step_size
     use clm_varcon         , only : hfus, cpliq, cpice, tkwat, tkice, denice
     use clm_varcon         , only : vkc, grav, denh2o, tfrz, cnfac
-    use clm_varctl         , only : iulog, use_lch4
+    use clm_varctl         , only : use_lch4
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds	  
@@ -1074,6 +1076,7 @@ subroutine SoilThermProp_Lake (bounds,  num_lakec, filter_lakec, tk, cv, tktopso
      use clm_varcon  , only : denh2o, denice, tfrz, tkwat, tkice, tkair
      use clm_varcon  , only : cpice,  cpliq, thk_bedrock, csol_bedrock
      use clm_varpar  , only : nlevsno, nlevsoi, nlevgrnd
+     use clm_varctl  , only : snow_thermal_cond_lake_method
      !
      ! !ARGUMENTS:
      type(bounds_type)      , intent(in)  :: bounds  
@@ -1170,7 +1173,20 @@ subroutine SoilThermProp_Lake (bounds,  num_lakec, filter_lakec, tk, cv, tktopso
              ! Only examine levels from snl(c)+1 -> 0 where snl(c) < 1
              if (snl(c)+1 < 1 .AND. (j >= snl(c)+1) .AND. (j <= 0)) then
                 bw = (h2osoi_ice(c,j)+h2osoi_liq(c,j))/dz(c,j)
-                thk(c,j) = tkair + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkair)
+                select case (snow_thermal_cond_lake_method)
+                case ('Jordan1991')
+                   thk(c,j) = tkair + (7.75e-5_r8 *bw + 1.105e-6_r8*bw*bw)*(tkice-tkair)
+                case ('Sturm1997')
+                   if (bw <= 156) then
+                      thk(c,j) = 0.023 + 0.234*(bw/1000)
+                   else
+                      thk(c,j) = 0.138 - 1.01*(bw/1000) +(3.233*((bw/1000)*(bw/1000)))
+                   end if
+                case default
+                   write(iulog,*) ' ERROR: unknown snow_thermal_cond_lake_method value: ', snow_thermal_cond_lake_method
+                   call endrun(msg=errMsg(sourcefile, __LINE__))
+                end select
+
              end if
 
           end do
diff --git a/src/clm5/biogeophys/SoilTemperatureMod.F90 b/src/clm5/biogeophys/SoilTemperatureMod.F90
@@ -47,7 +47,7 @@ module SoilTemperatureMod
   !     o The thermal conductivity of soil is computed from
   !       the algorithm of Johansen (as reported by Farouki 1981), and the
   !       conductivity of snow is from the formulation used in
-  !       SNTHERM (Jordan 1991).
+  !       Sturm (1997) or Jordan (1991) p. 18 depending on namelist option.
   !     o Boundary conditions:
   !       F = Rnet - Hg - LEg (top),  F= 0 (base of the soil column).
   !     o Soil / snow temperature is predicted from heat conduction
@@ -598,7 +598,7 @@ subroutine SoilThermProp (bounds,  num_nolakec, filter_nolakec, &
     use clm_varcon      , only : denh2o, denice, tfrz, tkwat, tkice, tkair, cpice,  cpliq, thk_bedrock, csol_bedrock
     use landunit_varcon , only : istice_mec, istwet
     use column_varcon   , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv
-    use clm_varctl      , only : iulog
+    use clm_varctl      , only : iulog, snow_thermal_cond_method, snow_thermal_cond_glc_method
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds                   
@@ -712,12 +712,50 @@ subroutine SoilThermProp (bounds,  num_nolakec, filter_nolakec, &
                endif
             endif
 
-            ! Thermal conductivity of snow, which from Jordan (1991) pp. 18
+            ! Thermal conductivity of snow
             ! Only examine levels from snl(c)+1 -> 0 where snl(c) < 1
             if (snl(c)+1 < 1 .AND. (j >= snl(c)+1) .AND. (j <= 0)) then  
                bw(c,j) = (h2osoi_ice(c,j)+h2osoi_liq(c,j))/(frac_sno(c)*dz(c,j))
-               thk(c,j) = tkair + (7.75e-5_r8 *bw(c,j) + 1.105e-6_r8*bw(c,j)*bw(c,j))*(tkice-tkair)
-            end if
+               l = col%landunit(c)
+
+               ! Select method over glacier land unit 
+               if (lun%itype(l) == istice_mec) then
+                  select case (snow_thermal_cond_glc_method)
+                  ! TODO, this code duplication isn't ideal and should likely be in it's own subroutine
+                  case('Jordan1991')
+                     thk(c,j) = tkair + (7.75e-5_r8 *bw(c,j) + 1.105e-6_r8*bw(c,j)*bw(c,j))*(tkice-tkair)
+                  case ('Sturm1997')
+                     if (bw(c,j) <= 156) then !LMW or 0.156 ?
+                        thk(c,j) = 0.023 + 0.234*(bw(c,j)/1000) !LMW - units changed by VRD
+                     else 
+                        thk(c,j) = 0.138 - 1.01*(bw(c,j)/1000) +(3.233*((bw(c,j)/1000)*(bw(c,j)/1000))) 
+                     end if
+                  case default
+                     write(iulog,*) ' ERROR: unknown snow_thermal_cond_glc_method value: ', snow_thermal_cond_glc_method
+                     call endrun(msg=errMsg(sourcefile, __LINE__))
+                  end select
+
+               else
+                  select case (snow_thermal_cond_method)
+                  case ('Jordan1991')
+                     thk(c,j) = tkair + (7.75e-5_r8 *bw(c,j) + 1.105e-6_r8*bw(c,j)*bw(c,j))*(tkice-tkair)
+                  case ('Sturm1997')
+                     ! Implemented by Vicky Dutch (VRD), Nick Rutter, and
+                     ! Leanne Wake (LMW)
+                     ! https://tc.copernicus.org/articles/16/4201/2022/
+                     ! See also https://tc.copernicus.org/articles/19/1539/2025/
+                     ! Code provided by Adrien Dams to Will Wieder
+                     if (bw(c,j) <= 156) then !LMW or 0.156 ?
+                        thk(c,j) = 0.023 + 0.234*(bw(c,j)/1000) !LMW - units changed by VRD
+                     else
+                        thk(c,j) = 0.138 - 1.01*(bw(c,j)/1000) +(3.233*((bw(c,j)/1000)*(bw(c,j)/1000)))
+                     end if
+                  case default
+                     write(iulog,*) ' ERROR: unknown snow_thermal_cond_method value: ', snow_thermal_cond_method
+                     call endrun(msg=errMsg(sourcefile, __LINE__))
+                  end select
+               end if ! close land unit if statement
+          end if
 
          end do
       end do
diff --git a/src/clm5/main/clm_varctl.F90 b/src/clm5/main/clm_varctl.F90
@@ -176,6 +176,11 @@ module clm_varctl
   ! use subgrid fluxes
   integer,  public :: subgridflag = 1                   
 
+  ! which snow thermal conductivity parameterization to use
+  character(len=25), public :: snow_thermal_cond_method
+  character(len=25), public :: snow_thermal_cond_glc_method
+  character(len=25), public :: snow_thermal_cond_lake_method
+
   ! true => write global average diagnostics to std out
   logical,  public :: wrtdia       = .false.            
 
diff --git a/src/clm5/main/controlMod.F90 b/src/clm5/main/controlMod.F90
@@ -204,6 +204,12 @@ subroutine control_init( )
          albice, soil_layerstruct, subgridflag, &
          irrigate, run_zero_weight_urban, all_active
 
+    ! Snow thermal conductivity
+    namelist /clm_inparm/  &
+         snow_thermal_cond_method, &
+         snow_thermal_cond_glc_method, &
+         snow_thermal_cond_lake_method
+
     ! vertical soil mixing variables
     namelist /clm_inparm/  &
          som_adv_flux, max_depth_cryoturb
@@ -718,6 +724,9 @@ subroutine control_spmd()
     call mpi_bcast (co2_ppmv, 1, MPI_REAL8,0, mpicom, ier)
     call mpi_bcast (albice, 2, MPI_REAL8,0, mpicom, ier)
     call mpi_bcast (soil_layerstruct,len(soil_layerstruct), MPI_CHARACTER, 0, mpicom, ier)
+    call mpi_bcast (snow_thermal_cond_method, len(snow_thermal_cond_method), MPI_CHARACTER, 0, mpicom, ier)
+    call mpi_bcast (snow_thermal_cond_glc_method, len(snow_thermal_cond_glc_method), MPI_CHARACTER, 0, mpicom, ier)
+    call mpi_bcast (snow_thermal_cond_lake_method, len(snow_thermal_cond_lake_method), MPI_CHARACTER, 0, mpicom, ier)
 
     ! snow pack variables
     call mpi_bcast (nlevsno, 1, MPI_INTEGER, 0, mpicom, ier)
