Alternative congelation formulation following Plante et al. 2024

columnphysics/icepack_parameters.F90

@@ -161,6 +161,11 @@ module icepack_parameters
                                   ! only for use with tr_aero or tr_zaero
          conserv_check = .false.  ! if true, do conservations checks and abort
 
+      character(len=char_len), public :: &
+         congel_freeze  = 'two-step'  ! congelation computation
+                                      ! 'two-step' = original formulation
+                                      ! 'one-step' = Plante et al, The Cryosphere, 2024
+
       character(len=char_len), public :: &
          tfrz_option  = 'mushy'   ! form of ocean freezing temperature
                                   ! 'minus1p8' = -1.8 C
@@ -465,7 +470,7 @@ subroutine icepack_init_parameters(   &
          atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, &
          atmiter_conv_in, calc_dragio_in, &
          tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, &
-         saltflux_option_in, &
+         saltflux_option_in, congel_freeze_in, &
          floeshape_in, wave_spec_in, wave_spec_type_in, nfreq_in, &
          dpscale_in, rfracmin_in, rfracmax_in, pndaspect_in, hs1_in, hp1_in, &
          bgc_flux_type_in, z_tracers_in, scale_bgc_in, solve_zbgc_in, &
@@ -575,15 +580,20 @@ subroutine icepack_init_parameters(   &
          phi_i_mushy_in           ! liquid fraction of congelation ice
 
       character(len=*), intent(in), optional :: &
-         tfrz_option_in              ! form of ocean freezing temperature
-                                     ! 'minus1p8' = -1.8 C
-                                     ! 'linear_salt' = -depressT * sss
-                                     ! 'mushy' conforms with ktherm=2
+         congel_freeze_in         ! congelation computation
+                                  ! 'two-step' = original formulation
+                                  ! 'one-step' = Plante et al, The Cryosphere, 2024
 
       character(len=*), intent(in), optional :: &
-         saltflux_option_in         ! Salt flux computation
-                                    ! 'constant' reference value of ice_ref_salinity
-                                    ! 'prognostic' prognostic salt flux
+         tfrz_option_in           ! form of ocean freezing temperature
+                                  ! 'minus1p8' = -1.8 C
+                                  ! 'linear_salt' = -depressT * sss
+                                  ! 'mushy' conforms with ktherm=2
+
+      character(len=*), intent(in), optional :: &
+         saltflux_option_in       ! Salt flux computation
+                                  ! 'constant' reference value of ice_ref_salinity
+                                  ! 'prognostic' prognostic salt flux
 
 !-----------------------------------------------------------------------
 ! Parameters for radiation
@@ -949,6 +959,7 @@ subroutine icepack_init_parameters(   &
       if (present(highfreq_in)          ) highfreq         = highfreq_in
       if (present(natmiter_in)          ) natmiter         = natmiter_in
       if (present(atmiter_conv_in)      ) atmiter_conv     = atmiter_conv_in
+      if (present(congel_freeze_in)     ) congel_freeze    = congel_freeze_in
       if (present(tfrz_option_in)       ) tfrz_option      = tfrz_option_in
       if (present(saltflux_option_in)   ) saltflux_option  = saltflux_option_in
       if (present(kitd_in)              ) kitd             = kitd_in
@@ -1192,7 +1203,7 @@ subroutine icepack_query_parameters(   &
          atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, &
          atmiter_conv_out, calc_dragio_out, &
          tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, &
-         saltflux_option_out, &
+         saltflux_option_out, congel_freeze_out, &
          floeshape_out, wave_spec_out, wave_spec_type_out, nfreq_out, &
          dpscale_out, rfracmin_out, rfracmax_out, pndaspect_out, hs1_out, hp1_out, &
          bgc_flux_type_out, z_tracers_out, scale_bgc_out, solve_zbgc_out, &
@@ -1311,16 +1322,21 @@ subroutine icepack_query_parameters(   &
          phi_i_mushy_out           ! liquid fraction of congelation ice
 
       character(len=*), intent(out), optional :: &
-         tfrz_option_out              ! form of ocean freezing temperature
-                                      ! 'minus1p8' = -1.8 C
-                                      ! 'constant' = Tocnfrz
-                                      ! 'linear_salt' = -depressT * sss
-                                      ! 'mushy' conforms with ktherm=2
+         congel_freeze_out         ! congelation computation
+                                   ! 'two-step' = original formulation
+                                   ! 'one-step' = Plante et al, The Cryosphere, 2024
+
+      character(len=*), intent(out), optional :: &
+         tfrz_option_out           ! form of ocean freezing temperature
+                                   ! 'minus1p8' = -1.8 C
+                                   ! 'constant' = Tocnfrz
+                                   ! 'linear_salt' = -depressT * sss
+                                   ! 'mushy' conforms with ktherm=2
 
       character(len=*), intent(out), optional :: &
-         saltflux_option_out         ! Salt flux computation
-                                     ! 'constant' reference value of ice_ref_salinity
-                                     ! 'prognostic' prognostic salt flux
+         saltflux_option_out       ! Salt flux computation
+                                   ! 'constant' reference value of ice_ref_salinity
+                                   ! 'prognostic' prognostic salt flux
 
 
 !-----------------------------------------------------------------------
@@ -1719,6 +1735,7 @@ subroutine icepack_query_parameters(   &
       if (present(highfreq_out)          ) highfreq_out     = highfreq
       if (present(natmiter_out)          ) natmiter_out     = natmiter
       if (present(atmiter_conv_out)      ) atmiter_conv_out = atmiter_conv
+      if (present(congel_freeze_out)     ) congel_freeze_out = congel_freeze
       if (present(tfrz_option_out)       ) tfrz_option_out  = tfrz_option
       if (present(saltflux_option_out)   ) saltflux_option_out = saltflux_option
       if (present(kitd_out)              ) kitd_out         = kitd
@@ -1926,6 +1943,7 @@ subroutine icepack_write_parameters(iounit)
         write(iounit,*) "  highfreq   = ", highfreq
         write(iounit,*) "  natmiter   = ", natmiter
         write(iounit,*) "  atmiter_conv = ", atmiter_conv
+        write(iounit,*) "  congel_freeze = ", trim(congel_freeze)
         write(iounit,*) "  tfrz_option= ", trim(tfrz_option)
         write(iounit,*) "  saltflux_option = ", trim(saltflux_option)
         write(iounit,*) "  kitd       = ", kitd

columnphysics/icepack_therm_vertical.F90

@@ -27,7 +27,7 @@ module icepack_therm_vertical
       use icepack_parameters, only: ustar_min, fbot_xfer_type, formdrag, calc_strair
       use icepack_parameters, only: rfracmin, rfracmax, dpscale, frzpnd, snwgrain, snwlvlfac
       use icepack_parameters, only: phi_i_mushy, floeshape, floediam, use_smliq_pnd, snwredist
-      use icepack_parameters, only: saltflux_option
+      use icepack_parameters, only: saltflux_option, congel_freeze
       use icepack_parameters, only: icepack_chkoptargflag
 
       use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso
@@ -46,6 +46,7 @@ module icepack_therm_vertical
       use icepack_mushy_physics, only: icepack_mushy_temperature_mush
       use icepack_mushy_physics, only: liquidus_temperature_mush
       use icepack_mushy_physics, only: icepack_enthalpy_mush, enthalpy_of_melting
+      use icepack_mushy_physics, only: enthalpy_mush_liquid_fraction, enthalpy_brine
 
       use icepack_aerosol, only: update_aerosol
       use icepack_isotope, only: update_isotope
@@ -85,7 +86,7 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
                                   fsnow,       fpond,     &
                                   fbot,        Tbot,      &
                                   Tsnice,      sss,       &
-                                  rsnw,                   &
+                                  sst,         rsnw,      &
                                   lhcoef,      shcoef,    &
                                   fswsfc,      fswint,    &
                                   Sswabs,      Iswabs,    &
@@ -164,6 +165,7 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
       real (kind=dbl_kind), intent(in) :: &
          fbot    , & ! ice-ocean heat flux at bottom surface (W/m^2)
          Tbot    , & ! ice bottom surface temperature (deg C)
+         sst     , & ! sea surface temperature (C)
          sss         ! ocean salinity
 
       ! coupler fluxes to atmosphere
@@ -403,6 +405,7 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
                              congel,      snoice,    &
                              mlt_onset,   frz_onset, &
                              zSin,        sss,       &
+                             sst,                    &
                              dsnow,       rsnw)
       if (icepack_warnings_aborted(subname)) return
 
@@ -1013,6 +1016,7 @@ subroutine thickness_changes (nilyr,     nslyr,    &
                                     congel,    snoice,   &
                                     mlt_onset, frz_onset,&
                                     zSin,      sss,      &
+                                    sst,                 &
                                     dsnow,     rsnw)
 
       integer (kind=int_kind), intent(in) :: &
@@ -1082,6 +1086,7 @@ subroutine thickness_changes (nilyr,     nslyr,    &
          zSin            ! ice layer salinity (ppt)
 
       real (kind=dbl_kind), intent(in) :: &
+         sst         , & ! sea surface temperature (C)
          sss             ! ocean salinity (PSU)
 
       ! local variables
@@ -1131,9 +1136,9 @@ subroutine thickness_changes (nilyr,     nslyr,    &
          qmlt            ! enthalpy of melted ice (J m-3) = zero in BL99 formulation
 
       real (kind=dbl_kind) :: &
-         qbotm       , &
-         qbotp       , &
-         qbot0       , &
+         qbotm       , & ! enthalpy of newly formed congelation ice
+         qbotp       , & ! enthalpy needed to grow new congelation ice (includes ocean enthalpy)
+         qbotw       , & ! enthalpy transferred to ocean during congelation freezing
          mass        , & ! total  mass from snow density tracers (kg/m^2)
          massi       , & ! ice mass change factor
          tmp1            ! temporary scalar
@@ -1254,15 +1259,22 @@ subroutine thickness_changes (nilyr,     nslyr,    &
 
       if (ktherm == 2) then
 
-         qbotm = icepack_enthalpy_mush(Tbot, sss)
-         qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy)
-         qbot0 = qbotm - qbotp
-
-         dhi = ebot_gro / qbotp     ! dhi > 0
-
+         if (congel_freeze == 'one-step') then
+            ! Plante et al., The Cryosphere, 18, 1685-1708, 2024
+            qbotm = enthalpy_mush_liquid_fraction(Tbot, phi_i_mushy)
+            qbotw = enthalpy_brine(sst)
+            qbotp = qbotm - qbotw
+            dhi = ebot_gro / qbotp     ! dhi > 0
+            hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss*phi_i_mushy
+         else ! two-step
+            qbotm = icepack_enthalpy_mush(Tbot, sss)
+            qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy)
+            qbotw = qbotm - qbotp
+            dhi = ebot_gro / qbotp     ! dhi > 0
+            hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss
+         endif
          hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm
-         hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss
-         emlt_ocn = emlt_ocn - qbot0 * dhi
+         emlt_ocn = emlt_ocn - qbotw * dhi
 
       else
 
@@ -2653,7 +2665,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
                                  fsnow=fsnow,         fpond=fpond,             &
                                  fbot=fbot,           Tbot=Tbot,               &
                                  Tsnice=Tsnice,       sss=sss,                 &
-                                 rsnw=rsnw,                                    &
+                                 sst=sst,             rsnw=rsnw,               &
                                  lhcoef=lhcoef,       shcoef=shcoef,           &
                                  fswsfc=fswsfcn  (n), fswint=fswintn      (n), &
                                  Sswabs=Sswabsn(:,n), Iswabs=Iswabsn    (:,n), &

configuration/driver/icedrv_init.F90

@@ -99,7 +99,7 @@ subroutine input_data
          natmiter, kitd, kcatbound
 
       character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
-         cpl_frazil, tfrz_option, saltflux_option, &
+         cpl_frazil, congel_freeze, tfrz_option, saltflux_option, &
          frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table
 
       logical (kind=log_kind) :: sw_redist, use_smliq_pnd, snwgrain, update_ocn_f
@@ -174,12 +174,12 @@ subroutine input_data
         update_ocn_f,    l_mpond_fresh,   ustar_min,       &
         fbot_xfer_type,  oceanmixed_ice,  emissivity,      &
         formdrag,        highfreq,        natmiter,        &
-        atmiter_conv,    calc_dragio,                      &
+        atmiter_conv,    calc_dragio,     congel_freeze,   &
         tfrz_option,     saltflux_option, ice_ref_salinity, &
         default_season,  wave_spec_type,  cpl_frazil,      &
         precip_units,    fyear_init,      ycycle,          &
         atm_data_type,   ocn_data_type,   bgc_data_type,   &
-        lateral_flux_type,                                &
+        lateral_flux_type,                                 &
         atm_data_file,   ocn_data_file,   bgc_data_file,   &
         ice_data_file,                                     &
         atm_data_format, ocn_data_format, bgc_data_format, &
@@ -225,6 +225,7 @@ subroutine input_data
            dSdt_slow_mode_out=dSdt_slow_mode, &
            phi_c_slow_mode_out=phi_c_slow_mode, Tliquidus_max_out=Tliquidus_max, &
            phi_i_mushy_out=phi_i_mushy, conserv_check_out=conserv_check, &
+           congel_freeze_out=congel_freeze, &
            tfrz_option_out=tfrz_option, saltflux_option_out=saltflux_option, &
            ice_ref_salinity_out=ice_ref_salinity, kalg_out=kalg, &
            fbot_xfer_type_out=fbot_xfer_type, puny_out=puny, &
@@ -791,6 +792,7 @@ subroutine input_data
          write(nu_diag,1005) ' hi_min                    = ', hi_min
          write(nu_diag,1030) ' fbot_xfer_type            = ', trim(fbot_xfer_type)
          write(nu_diag,1010) ' oceanmixed_ice            = ', oceanmixed_ice
+         write(nu_diag,1030) ' congel_freeze             = ', trim(congel_freeze)
          write(nu_diag,1030) ' tfrz_option               = ', trim(tfrz_option)
          write(nu_diag,*)    ' saltflux_option           = ', &
                                trim(saltflux_option)
@@ -979,6 +981,7 @@ subroutine input_data
            dSdt_slow_mode_in=dSdt_slow_mode, &
            phi_c_slow_mode_in=phi_c_slow_mode, Tliquidus_max_in=Tliquidus_max, &
            phi_i_mushy_in=phi_i_mushy, conserv_check_in=conserv_check, &
+           congel_freeze_in=congel_freeze, &
            tfrz_option_in=tfrz_option, saltflux_option_in=saltflux_option, &
            ice_ref_salinity_in=ice_ref_salinity, kalg_in=kalg, &
            fbot_xfer_type_in=fbot_xfer_type, &

configuration/scripts/icepack_in

@@ -113,6 +113,7 @@
     cpl_frazil      = 'fresh_ice_correction'
     update_ocn_f    = .false.
     l_mpond_fresh   = .false.
+    congel_freeze   = 'two-step'
     tfrz_option     = 'mushy'
     ice_ref_salinity = 4.0
     saltflux_option = 'constant'

configuration/scripts/machines/env.chicoma_intel

@@ -63,7 +63,7 @@ setenv ICE_MACHINE_WKDIR /lustre/scratch5/$user/ICEPACK_RUNS
 setenv ICE_MACHINE_INPUTDATA /usr/projects/climate/eclare/DATA/Consortium
 setenv ICE_MACHINE_BASELINE /lustre/scratch5/$user/ICEPACK_BASELINE
 setenv ICE_MACHINE_SUBMIT "sbatch "
-setenv ICE_MACHINE_ACCT t23_cice
+setenv ICE_MACHINE_ACCT t24_cice
 setenv ICE_MACHINE_QUEUE "debug"
 setenv ICE_MACHINE_TPNODE 128    # tasks per node
 setenv ICE_MACHINE_BLDTHRDS 12

configuration/scripts/options/set_nml.congel

@@ -0,0 +1 @@
+  congel_freeze      = 'one-step'

configuration/scripts/tests/base_suite.ts

@@ -23,6 +23,7 @@ smoke          col     1x1        debug,run1year,fluxopenw
 smoke          col     1x1        debug,run1year,dyn,fluxopenw
 smoke          col     1x1        debug,run1year,debug,snicartest
 smoke          col     1x1        debug,run1year,debug,snicar
+smoke          col     1x1        debug,run1year,debug,congel
 restart        col     1x1        debug
 restart        col     1x1        diag1
 restart        col     1x1        pondlvl
@@ -43,4 +44,5 @@ restart        col     1x1        modal
 restart        col     1x1        fluxopenw
 restart        col     1x1        snicartest
 restart        col     1x1        snicar
+restart        col     1x1        congel
 

doc/source/master_list.bib

@@ -32,7 +32,7 @@ @string{JGR
 @string{JGRO = {J. Geophys. Res. Oceans}}
 @string{JGRB = {J. Geophys. Res. Biogeo.}}
 @string{JGRA = {J. Geophys. Res. Atmos.}}
-@string{JCT = {J. Comput. Phys.}}
+@string{TC = {The Cryosphere}}
 @string{QJRMS = {Quart. J. Royal Met. Soc.}}
 @string{GRL = {Geophys. Res. Lett.}}
 @string{JAS = {J. Atmos. Sci.}}
@@ -244,7 +244,7 @@ @Article{Maykut95
 @Article{Murray96
   author = "R.J. Murray",
   title = "{Explicit generation of orthogonal grids for ocean models}",
-  journal = JCT,
+  journal = JCP,
   year = {1996},
   volume = {126},
   pages = {251-273},
@@ -805,13 +805,23 @@ @article{Roberts19
 @article{Dang19,
 author = {Dang, C. and Zender, C. S. and Flanner, M. G.},
 title = {Intercomparison and improvement of two-stream shortwave radiative transfer schemes in Earth system models for a unified treatment of cryospheric surfaces},
-journal = {The Cryosphere},
+journal = {TC},
 volume = 13,
 pages = {2325--2343},
 doi = {10.5194/tc-13-2325-2019},
 year = {2019}
 }
 
+@article{Plante24,
+author = {Plante, M. and J.-F. Lemieux and L. B. Tremblay and A. Tivy and J. Angnatok and F. Roy and G. Smith and F. Dupont and A. K. Turner},
+title = {Using Icepack to reproduce ice mass balance buoy observations in landfast ice: improvements from the mushy-layer thermodynamcis},
+journal = {TC},
+volume = 18,
+pages = {1685--1708},
+doi = {10.5194/tc-18-1685-2024},
+year = {20124}
+}
+
 %  **********************************************
 %  For new entries, see example entry in BIB_TEMPLATE.txt
 %  **********************************************

doc/source/science_guide/sg_thermo.rst

@@ -1737,9 +1737,21 @@ conductive heat flux at the bottom surface:
 
 If ice is melting at the bottom surface, :math:`q`
 in Equation :eq:`bottom-melting` is the enthalpy of the bottom ice layer. If
-ice is growing, :math:`q` is the enthalpy of new ice with temperature
-:math:`T_f` and salinity :math:`S_{max}` (``ktherm`` = 1) or ocean surface
-salinity (``ktherm`` = 2). This ice is added to the bottom layer.
+ice is growing, :math:`q` is the enthalpy of new ice added to the bottom layer with
+temperature :math:`T_f` and salinity :math:`S_{max}` (``ktherm`` = 1) for the Bitz99
+:cite:`Bitz99` formulation.  The original mushy thermodynamics
+(``ktherm`` = 2, ``congel_freeze`` = 'two-step') formed new ice in two steps, first
+moving the lower ice boundary into the ocean to form a mushy layer with an initial
+liquid fraction :math:`phi_{init}`, then freezing the new ice in the next step.  In
+this case, the freezing temperature is calculated using the ocean surface salinity
+:math:`SSS`.  A second mushy thermo method
+(``ktherm`` = 2, ``congel_freeze`` = 'one-step') freezes the ice immediately using
+the brine salinity :math:`phi_{init} * SSS` :cite:`Plante24`.
+In the two-step method, enthalpy not
+yet used to freeze the ice is returned to the ocean, causing frazil ice to form
+instead.  Together, the congelation and frazil ice add up to a similar total amount
+of new ice, but the differing freezing mechanism complicates comparisons with observational
+data and there is an unnecessary lag in ice-ocean coupling.
 
 In general, frazil ice formed in the ocean is added to the thinnest ice
 category. The new ice is grown in the open water area of the grid cell