Ice deposition (the inverse of sublimation) in coupled simulations

Versions v4.2.x
1

Hello,

We have NEMO4.2 coupled to our in-house atmosphere through our in-house coupler. (CanAM, CanCPL)

We are having some issues getting a realistic sea ice cover in coupled simulations, (our sea ice volume explodes 20 fold in 100 years) while forced simulations are giving a fairly good sea ice cover.

So we are investigating the fields going through our coupler and compare them to ones from the forced simulation (ERA5) through the NEMO bulk scheme “sbcblk”.

I noticed that the evap_ice field from our coupler is positive and negative while in the forced simulation, it is only positive. I then noticed the sentence in the SI3 documentation saying that

Sublimation converts the evaporation mass flux evap_ice into a snow loss
assuming constant density. The energy cost of sublimation is included in the latent
heat flux and needs not being accounted for. Evaporation removes snow mass, and
is handled layer by layer. Deposition of snow if any should be routed through the
solid precipitation mass flux. (p54, sect 6.4.3)

So I assume that we would move all deposition (negative evap_ice) to the sprecip while keeping the sublimation (positive evap_ice) in `evap_ice in our coupler

My issue here is that sprecip is also for fluxes going into the ocean, so I am not sure how to proceed here because it is realistic for snowfall but not for deposition, we should not deposit ice onto the ocean. We use the conservative setting for sn_rcv_emp in namsbc_cpl in namelist_cfg.

Also some comments in the code would suggest that deposition could be handled by SI3 through evap_ice? Has anyone experimented with this?

Thanks in advance!

2

Reporting here if someone else is interested:

After a careful reading of the code, I think that SI3 can handle deposition.

If there is some, a corresponding amount of snow is added to the (0) snow layer where new snow is also added/removed from snowfall/sublimation. This new snow is added at the same temperature as the current snow, or 0degC if no snow. If snow is left in this (0) layer at the end of sbc timestep (after surface melt), the snow is moved to the “normal” snow layers ( layer (1) and above (if more than one layer is used) ).

I did not see any other mass or energy issue that would be created by having evap_ice < 0. The associated energy and mass flux should be taken into account by hfx_sub and wfx_snw_sub.

I will test to see if what is given by our coupler corresponds to what is in the sublimation diagnostics.

Maybe I missed or forgot something, if anyone has any idea, please let me know =)

I annotated the code for my analysis, details below:

Summary

My comments after !==>DR:

Set in sbccpl from coupler values

            zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) !==>DR: <0 if deposition
            zevap_ice_total(:,:) = zevap_ice(:,:,1) !DR: Same
            DO jl = 2, jpl
               zevap_ice(:,:,jl) = zevap_ice(:,:,1) !DR just the same thing again
            END DO
[...]
      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
         ! For conservative case zemp_ice has not been defined yet. Do it now.
         zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:) !==>DR: incuded in `emp_ice` but `emp_ice` does not do much in SI3, it seems to be used in CICE
      ENDIF

      ! --- evaporation over ocean (used later for qemp) --- !
      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:) !==>DR: evap_tot = evap_oce * (1-a) + evap_ice *a, so if evap_ice <0 then evap_oce > evap_tot 

      evap_ice(:,:,:) = zevap_ice(:,:,:) !==>DR: temp. var. to global var.

[...]

      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea' , zevap_ice_total(:,:) * picefr(:,:) * tmask(:,:,1)     )  ! Sublimation over sea-ice (cell average) !==>DR: Also deposition
      IF( iom_use('evap_ao_cea') )   CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
         &                                                         - zevap_ice_total(:,:) * picefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average) 
      IF (        iom_use('hflx_evap_cea') )    &                                                    ! heat flux from evap (cell average)
         &   CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:) )  &
         &                                  * zcptn(:,:) * tmask(:,:,1) )
      IF(         iom_use('hflx_subl_cea') )    &                                                    ! heat flux from sublimation
         &   CALL iom_put('hflx_subl_cea' ,   SUM( qevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) * tmask(:,:,1) )

in ice_thd

         CALL tab_2d_1d( npti, nptidx(1:npti), evap_ice_1d   (1:npti), evap_ice(:,:,kl)     ) !==>DR: From 2D to 1D

in ice_thd_dh

      ! Snow sublimation
      !-----------------
      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
      !    comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
      zdeltah   (1:npti) = 0._wp ! total snow thickness that sublimates, < 0 !==>DR: > 0 if deposition
      zevap_rema(1:npti) = 0._wp
      DO ji = 1, npti ! horizontal grid in 1D
         zdeltah   (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rDt_ice, - h_s_1d(ji) )   ! amount of snw that sublimates, ==>DR: evap_ice_1d<0 so zdeltah>0
         zevap_rema(ji) = evap_ice_1d(ji) * rDt_ice + zdeltah(ji) * rhos               ! remaining evap in kg.m-2 (used for ice sublimation later on) ==>DR: because no limit up, zevap_rema = 0
      END DO
      
      DO jk = 0, nlay_s ! snow layers just one but 0 for the current snowfall 1 for the layer
         DO ji = 1, npti ! horizontal grid in 1D
            zdum = MAX( -zh_s(ji,jk), zdeltah(ji) ) ! snow layer thickness that sublimates, < 0 ==> DR: zh_s>0 so zh_s(ji,jk) <0 so zh_s(ji,jk) < zdeltah(ji) always so zdum = zdeltah(ji) > 0
            !
            hfx_sub_1d    (ji) = hfx_sub_1d    (ji) + ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice  ! Heat flux of snw that sublimates [W.m-2], < 0 ==>DR: >0 for depostion, 
            wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos        * zdum * a_i_1d(ji) * r1_Dt_ice  ! Mass flux by sublimation ==>DR: <0 for depostion negative substraction means addition
            
            ! update thickness
            h_s_1d(ji)    = MAX( 0._wp , h_s_1d(ji)    + zdum ) ! ==>DR: Adding total thickness zdum >0
            zh_s  (ji,jk) = MAX( 0._wp , zh_s  (ji,jk) + zdum ) ! ==>DR: Adding layer thickness zdum >0
            !!$ IF( zh_s(ji,jk) == 0._wp )   ze_s(ji,jk) = 0._wp
            
            ! update sublimation left for the next layer
            zdeltah(ji) = MIN( zdeltah(ji) - zdum, 0._wp ) !==>DR: must be zero for the deposition at the first layer, as  zdum = zdeltah(ji)
         END DO
      END DO
      
      [...]

      !                       ! ============ !
      !                       !     Ice      !
      !                       ! ============ !        
     
      [...]

      DO jk = 1, nlay_i ! ice layers
         DO ji = 1, npti ! horizontal grid in 1D
            [...]
            ! Ice sublimation
            ! ---------------
            zdum               = MAX( - zh_i(ji,jk) , - zevap_rema(ji) * r1_rhoi ) ! ==>DR: is zero for deposition, what comes after is all zero
            !
            hfx_sub_1d(ji)     = hfx_sub_1d(ji)     + e_i_1d(ji,jk) * zdum              * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0
            wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi          * zdum              * a_i_1d(ji) * r1_Dt_ice ! Mass flux > 0
            sfx_sub_1d(ji)     = sfx_sub_1d(ji)     - rhoi          * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0
            !                                                                                                      clem: flux is sent to the ocean for simplicity
            !                                                                                                            but salt should remain in the ice except
            !                                                                                                            if all ice is melted. => must be corrected
            ! update remaining mass flux and thickness
            zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi
            zh_i(ji,jk)    = MAX( 0._wp, zh_i(ji,jk) + zdum )
            h_i_1d(ji)     = MAX( 0._wp, h_i_1d(ji)  + zdum )
            dh_i_sub(ji)   = dh_i_sub(ji) + zdum
            
            ! update heat content (J.m-2) and layer thickness
            eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk)
            h_i_old (ji,jk) = h_i_old (ji,jk) + zdum

         END DO
      END DO
3

Hi Damien,

sorry for this late answer.

We did not pay too much attention to deposition actually. For forced simulations, we ensure that the latent heat qla_ice is always positive, hence no deposition is possible. For coupled simulation, there is no such a thing and you can have deposition indeed if evap_ice < 0. In that case, the mass flux seems to be taken into account ok. But then I notice a problem with the heat flux. They are two cases. 1) if there is snow precipitation, then ice “enthalpy” ze_s(0) takes the value of the precip which can be considered as ok. 2) if no snow fall, then ice “enthalpy” is 0 which is not ok. It should be equal to rhos*(Lfus-cpi*Ts). As it is written, it means that the temperature of this deposition is equal to Lfus/cpi, hence something like 150°C.

Solution:
One has to set a value for ze_s(0) when deposition AND no precip. It can just be rhos*Lfus if you consider that snow deposit is at 0°C. Or it can be ze_s(1). I am not sure about the best pick.

Thanks for pointing that out. I will add this item on my to-do list. Would you mind keep me in touch if you adopt one solution or the other and which one seems the most convincing?

Clem

4

Hi Clem,

Thanks a lot for your answer,

I see! Effectively, I did not look at the deposition of energy in the layer. Actually @cpelletier told me something similar privately and I did not have the time to go back to it.

Ok I will try to implement a fix, and of course I will keep you updated of the results.

Damien