WW3 coupling CICE

[anton-seaice]

Hi @ezhilsabareesh8

To turn on wave interactions with the sea-ice, it looks like we need:

• set 'wav_coupling_to_cice' in nuopc.runconfig, which will add Si_thick and Si_floediam to the advertised fields sent to the mediator from CICE
• turn on floe size distribution in ice_in. Set tr_fsd to true, and then set nfsd and wave_spec_type (see Icepack Docs). These last two fields we will need to experiment with a bit / talk to the end users. It also looks like if we use initial ice conditions (like we do currently) then the initial floe-size distrubution would be based on measurements from the Arctic, so we might want to experiment with setting ice_ic=default to see if the floe-sizes stabilise quicker.

That's would start us (I think) to get the sea-ice to reduce wave penetration.

I am not sure about getting the wave field to impact sea-ice breakup / floesize. That might be this 'Sw_elevation_spectrum' field, which would looks to be exported from WW.

I suggest we start by copying the CICE config from the 'MOM6-CICE6' repo, and then modify from there? (i.e. this commit)

EDIT: Also might need to set nfreq to the number of frequencies in ocean surface wave spectral forcing

[anton-seaice]

@ofa001 - You might have some thoughts?

[ofa001]

Hi @anton-seaice @ezhilsabareesh8 , thanks for the quick summary of the basic steps to start the coupling of the WW3 and CICE6 codes. Noah Day has been using 12 fsd categories, he did experiment with more, Lettie Roach alos used 12 so that would be a good number to choose and is the default, so they should be in the code or I can get him to send the boundaries values to you.

I think you are right the Sw_elevation spectum is the field exported, it should have a full directional spectrum, I will check and get back to you. We have been using a simpler approach in the wave data forced set up.

[ezhilsabareesh8]

Thanks @anton-seaice and @ofa001 . The wave damping and scattering by sea ice is not included in the current MOM6-CICE6-WW3 configuration, which can be enabled by ICx and ISx switches in WW3. There are five different
versions of dissipation processes activated with the switches IC1, IC2, IC3, IC4 and IC5 that can be combined with two different versions of scattering effects IS1 and IS2. For these models, ice-related inputs like ice thickness, density, effective shear modulus and viscosity are required. The spatial and temporal variability of ice-related inputs are only available for ice-concentration and ice-thickness which can be given as a input from the model and it is the most commonly used ice input parameters in dissipation models. The other parameters can be under option of ww3_shel.nml.

The Ice floe mean diameter is specifically required for IS2, a floe-size dependent scattering and dissipation model. I am checking whether this parameter can be provided to WW3 from the ice model through a coupler, or if only homogeneous input options are available.

[ofa001]

Hi @anton-seaice In relation to the earlier comment on the flds_wave option on the wave coupling is enabled in the CICE portion of the nuopc cap the SW_elevation_specrtum is the incoming wave field. Whilst Si_thick and SI_floedaim are the two outgoing ones in the ice_import_export fields that are adveriesed to the mediator. But as it says in the comments.
"the following are advertised but might not be connected if they are not advertised in the
! in the cmeps esmFldsExchange_xxx_mod.F90 that is model specific
! from wave"
All the fields might not be available for coupling.
This may be the case for @ezhilsabareesh8 next issue when looking at the different scattering problems, some of them need ice thickness to be passed from the sea ice model only one needs thickness and floe diameter. The ice model can be run without the FSD switched on and it sets a constant floe diameter of 300m. You don't get meaningful FSD distributions without waves.

[anton-seaice]

For these models, ice-related inputs like ice thickness, density, effective shear modulus and viscosity are required.

I believe ice density is fixed in CICE and we use the default of 917kg/m3. Ill need to do some more reading on shear modulus and viscosity.

The Ice floe mean diameter is specifically required for IS2, a floe-size dependent scattering and dissipation model. I am checking whether this parameter can be provided to WW3 from the ice model through a coupler, or if only homogeneous input options are available.

Isn't this the Si_floediam field ? It looks to me like it is in both the CICE and WW caps. As Siobhan says, we want to turn on floesize distributions if we have wave input to CICE, even if the wave model isn't using it.

[dougiesquire]

Whilst Si_thick and SI_floedaim are the two outgoing ones in the ice_import_export fields that are adveriesed to the mediator. But as it says in the comments.
"the following are advertised but might not be connected if they are not advertised in the
! in the cmeps esmFldsExchange_xxx_mod.F90 that is model specific
! from wave"
All the fields might not be available for coupling.

Si_thick and Si_floediam are handled in esmFldsExchange_cesm_mod.F90 (which we are using currently for ACCESS-OM3) so long as wav_coupling_to_cice=.true.

[ezhilsabareesh8]

I believe ice density is fixed in CICE and we use the default of 917kg/m3.

Thanks Anton, the ice density is set to the same value in WW3 regtests.

Si_thick and Si_floediam are handled in esmFldsExchange_cesm_mod.F90 (which we are using currently for ACCESS-OM3) so long as wav_coupling_to_cice=.true.

It is the same in WW3 wav_import_export.F90 as well.

[ofa001]

Thanks @dougiesquire esmFieldsExchange_xxx_mod.F90 is not one a look at as but I will check now I guess its on step closer to the mediator and I have been only checking the CICE code, I was just pointing out that we can run with variable ice thickness and a constant floe diameter if no FSD is switched on, but if we want to start engaging the wave coupling we need to switch all on. @ezhilsabareesh8 also has options in the wave code that only couples ice thickness not floe diameter but using both are newer schemes.

[ezhilsabareesh8]

@ezhilsabareesh8 also has options in the wave code that only couples ice thickness not floe diameter but using both are newer schemes.

Thanks @ofa001 I think floe diameter can also be passed to WW3 by enabling wav_coupling_to_cice=.true.. In wav_import_export.F90, Si_thick is passed to ICEP1 which is ice-thickness (refer here) and Si_floediam is passed to ICEP5 which is floe-diameter (refer here) in WW3.

[anton-seaice]

For reference, these are the switches in WW:

Selection of term for damping by sea ice:
ic0 No damping by sea ice.
ic1 Simple formulation.
ic2 Liu et al. formulation.
ic3 Wang and Shen formulation.
ic4 Frequency-dependent damping by sea ice.

Selection of term for scattering by sea ice:
is0 No scattering by sea ice.
is1 Diffusive scattering by sea ice (simple).
is2 Floe-size dependent scattering and dissipation.

They are described here in Section 2.4: Source terms for wave-ice interactions

https://raw.githubusercontent.com/wiki/NOAA-EMC/WW3/files/manual.pdf

The more I think about it the more I think its likely users will want to be able to configure these switches themselves.

[anton-seaice]

To summarise:

These are the cice config changes discussed in this thread:
fe3569c, this runs but I haven't looked at the output or the WW3 config changes needed.

A couple of notes for my reference later:

From Noah:

I agree that setting the wave_spec_type to random is best, the other settings are meant for testing.

The other thing wave_spec_type changes is in the generation of the sea surface height (SSH) field (icepack_wavefracspec.F90 L494). So, as you said Anton, you should get bit-for-bit with constant, but I’m not sure how ‘realistic’ that is.

and (re: nfsd)

The problem with increasing the maximum floe size is that if the wave spectra are too small the new floe size parameterisation (Shen et al. Annals, 2001) is turned off and CICE will assign the new ice to the largest floe size category (skewing the FSD even further). Hopefully, with a coupled wave model, some wave energy enters coastal polynyas, limiting the likelihood of these very large and thin-sheet ice.

From Ezhil

Regarding the wave frequency in WW3, the following parametrisation is used

&SPECTRUM_NML
SPECTRUM%XFR = 1.1. ! frequency increment
SPECTRUM%FREQ1 = 0.04118 ! first frequency (Hz)
SPECTRUM%NK = 25 ! number of frequencies (wavenumbers)
SPECTRUM%NTH = 24 ! number of direction bins
/

which is the same as the hardcoded following wave frequency in icepack_wavefracspec.F90 in CICE:

! hardwired for wave coupling with NIWA version of Wavewatch
! From Wavewatch, f(n+1) = C*f(n) where C is a constant set by the user
! These freq are for C = 1.1
wavefreq = (/0.04118, 0.045298, 0.0498278, 0.05481058, 0.06029164, &
0.06632081, 0.07295289, 0.08024818, 0.08827299, 0.09710029, &
0.10681032, 0.11749136, 0.1292405, 0.14216454, 0.15638101, &
0.17201911, 0.18922101, 0.20814312, 0.22895744, 0.25185317, &
0.27703848, 0.30474234, 0.33521661, 0.36873826, 0.40561208/)

If we decide we want to change these values, we should make them a configuration option and only have to set them in one place!

Similarly nfreq is 25 by default in CICE, which matches our WW3 config, but also needs configuring in two places to match.

[ezhilsabareesh8]

Exciting News: By activating the floe size distribution calculation in CICE through the configuration of tr_fsd = .true. and specifying nfsd = 12, the Sw_elevation_spectrum is exchanged from WW3 to CICE. Additionally, the reciprocal transfer of floe size diameter and ice thickness from CICE to WW3 is also operational and the WW3 output shows the variable floe size distribution. The following observations are made by reviewing the CICE codebase and the output of the simulations:

1. WW3 passes the 1-D power spectral density of surface elevation, E(f) (m^2/s), as Sw_elevation_spectrum to CICE. No phase information is passed along with the spectral density.

2. In CICE, the wave_spec_type option is limited to only two valid choices: none or random while passing spectrum from WW3 to CICE. All other specified options are not utilized in the FSD calculation.

3. When wave_spec_type = none is selected in CICE, the wave elevation spectrum from WW3 is used directly.
   A constant phase of phi = 2*pi*0.5 is added for all spatial locations. refer here

4. When wave_spec_type = random is chosen in CICE, the wave elevation spectrum from WW3 is used and a random phase of phi = 2*pi*rand_array is added for all spatial locations as defined here. The fracture code is then run to convergence by generating multiple realizations of sea surface height and adding resulting fractures to a histogram until successive histograms converge within a small error tolerance.

5. Workflow of Sw_elevation_spectrum Passing:

• Sw_elevation_spectrum from WW3 is assigned to the wave_spectrum variable in the ice_import_export.f90 file of CICE. Refer here
• The wave_spectrum variable is passed through the ice_arrays_column module by other subroutines in CICE. Refer here
• The step_dyn_wave subroutine is called during the advancement of the CICE model, which further calls icepack_step_wave_fracture to evolve the floe size distribution. Refer here
• wave_frac is called to calculate fracture histogram during which constant or random phase is being added.Refer here
• Within the wave_frac subroutine, the only utilized wave_spec_type option is random, triggering the convergence-seeking fracture code. Refer here

In the current workflow, only the random option for wave_spec_type is actively utilized. Routines associated with other wave_spec_type options, such as constant, are not invoked when passing the spectrum from WW3 to CICE.

The following are the outputs that I obtained while adding constant phase, i.e., wave_spec_type=none and random phase wave_spec_type=random

Experiment: 1 deg MOM6-CICE6-WW3 config, with IC1 dissipation and IS2 scattering model in WW3.

For wave_spec_type = none (constant phase)

For wave_spec_type = random (random phase)

The use of a random phase appears to introduce some unrealistic outcomes, notably evident in depicting complete ice coverage in the Sea of Okhotsk near Japan and in the Antarctic region. This might be attributed to exceptionally small wave spectra. Consequently, the activation of the new floe size parameterization (as proposed by Shen et al. in Annals, 2001) is disabled. In this scenario, CICE resorts to assigning the new ice to the largest floe size category, thereby exacerbating the skewness of the floe size distribution (FSD) as pointed out by Noah Day.

@ofa001, @aekiss and @NoahDay any thoughts on adding the random phase?

[anton-seaice]

Well done Ezhil!

I think we would have to look at sea ice concentration (aice) to know if that is complete ice coverage - or if the ice coverage is appropriate and only the FSD is incorrect? i.e. the Antarctic sea ice looks too low in the top left plots for wave_spec_type=constant

I will investigate the output afsd, it looks like it is weighted by the sea ice concentration (aice) (link), in CICE. I don't know if the mediator or WW3 corrects for that.

[NoahDay]

That's looking like great progress @ezhilsabareesh8!

I would be surprised if you see a large difference is sea ice concentration from the setting of wave_spec_type, since the new floe size parameterisation only determines which floe size category the new ice is allocated to (i.e., not the amount of new ice formed).

In my standalone CICE6 runs I have seen large floes occur in the open ocean < 1% concentrations, but not to this extent. It looks like you've plotted fsdrad from CICE which is independent of ice concentration (see here). If this is the case, I'm not sure why the floe diameter increases from 800m for wave_spec_type = none to 1500m for wave_spec_type = random. With the default bins and nfsd=12, the maximum floe size radius (fsdrad) should be ~850m.

@anton-seaice the output of afsd is weighted by the ice concentration in each thickness category, aicen (see here).

[ofa001]

Well done for getting a run underway @ezhilsabareesh8 I guess its for January which makes the result for the random setting look a little strange particularly in Antarctica and in the peripheral seas in the Arctic. I am assuming you only running for a few rime steps with 5 day in the heading. So it should be worth looking as @NoahDay says at why its jumped up the size of floe size maximum has jumped, after 1 day interval, or even checking after every coupling time step to check its not related to anything from the coupler side, as that's what new to us here.

[ezhilsabareesh8]

Thanks @anton-seaice, @NoahDay and @ofa001 .

It looks like you've plotted fsdrad from CICE which is independent of ice concentration

I plotted floe diameter (ICEF) from WW3 output. In the WW3 output, the floe diameter (Si_floediam) is derived from CICE6 (ice_import_export.F90 source). This computation in CICE involves the ice concentration, as outlined here. Subsequently, it gets updated within the IS2 scattering model of WW3.

With the default bins and nfsd=12, the maximum floe size radius (fsdrad) should be ~850m.

I anticipate that the maximum floe diameter is updated by IS2 model in WW3.

I guess its for January which makes the result for the random setting look a little strange particularly in Antarctica and in the peripheral seas in the Arctic. I am assuming you only running for a few rime steps with 5 day in the heading.

@ofa001 The plots correspond to the month of May. The coupled model was run for ten months, starting from January 1999 to October 1999.

[anton-seaice]

There are several namelist history parameters I missed in the first pass of looking at this:

https://github.com/CICE-Consortium/CICE/blob/7a4b95e6deec0ec72c1da35a23ae1eb3ffe3d077/configuration/scripts/ice_in#L710

There two for average radius and perimeter

         f_fsdrad      = 
         f_fsdperim    = 

are probably worth turning on

There are also these 3, which at some point we would want to figure out why these exist (and are different to the fields exported in ice_import_export.F90 and fsdrad/aice/hice)

            f_aice_ww  = 'x'
            f_diam_ww  = 'x'
            f_hice_ww  = 'x'

[anton-seaice]

For these models, ice-related inputs like ice thickness, density, effective shear modulus and viscosity are required.

@ofa001 Do you know what the default values of Ice elasticity (Young's Modulus) and Viscosity are? Do they vary with time?

It looks like with ktherm=1 (Bitz99) thermodynamics, CICE doesn't do anything with viscosity?

[anton-seaice]

@ezhilsabareesh8

This is what I was talking about with the three definitions for floe size diam/radius in cice.

In ice_import_export.F90:

 do j = jlo,jhi
      do i = ilo,ihi
       ...
            workx = c0
            worky = c0
            do n = 1, ncat
               do k = 1, nfsd
                  workx = workx + floe_rad_c(k) * aicen_init(i,j,n,iblk) * trcrn(i,j,nt_fsd+k-1,n,iblk)
                  worky = worky + aicen_init(i,j,n,iblk) * trcrn(i,j,nt_fsd+k-1,n,iblk)
               end do
            end do
            if (worky > c0) workx = c2*workx / worky
            floediam(i,j,iblk) = MAX(c2*floe_rad_c(1),workx)

In ice_history_fsd.F90:

  if (f_fsdrad(1:1) /= 'x') then
     do j = 1, ny_block
     do i = 1, nx_block
        worka(i,j) = c0
        if (aice(i,j,iblk) > puny) then
         do k = 1, nfsd_hist
            do n = 1, ncat_hist
              worka(i,j) = worka(i,j) &
                           + (trcrn(i,j,nt_fsd+k-1,n,iblk) * floe_rad_c(k) &
                           * aicen(i,j,n,iblk)/aice(i,j,iblk))
             end do
          end do
        endif
     end do
     end do
     call accum_hist_field(n_fsdrad, iblk, worka(:,:), a2D)

and:

if (f_diam_ww (1:1) /= 'x') then
     do j = 1, ny_block
     do i = 1, nx_block
        worka(i,j) = c0
        workb      = c0
        do n = 1, ncat_hist
        do k = 1, nfsd_hist
           workc = aicen(i,j,n,iblk)*trcrn(i,j,nt_fsd+k-1,n,iblk) &
                 / (c4*floeshape*floe_rad_c(k)**2)
           ! number-mean radius
           worka(i,j) = worka(i,j) + workc * floe_rad_c(k)
           ! normalization factor
           workb      = workb      + workc
        end do
        end do
        ! number-mean diameter, following Eqn. 5 in HT2017
        if (workb > c0) worka(i,j) = c2*worka(i,j) / workb
        worka(i,j) = MAX(c2*floe_rad_c(1),worka(i,j)) ! >= smallest resolved diameter
     end do
     end do
     call accum_hist_field(n_diam_ww, iblk, worka(:,:), a2D)
  endif

[aekiss]

Also see @dabail10's slides and presentation on the fields they couple between WW3 and CICE6.

[ezhilsabareesh8]

This issue has not been updated on GitHub for a while, so I am consolidating some of our email communications here for documentation purposes.

The initial choice of wave attenuation model is IC3, where ice is treated as a viscous elastic layer, combined with the IS2 floe-size dependent scattering model. This setup shows a variable floe size distribution in the WW3 output. Notably, switching to the random scheme of wave_spec_type in CICE while using IC3 and IS2 resulted in improvements in the Arctic region, as shown in the plots. This scheme prevents the infilling of the default high floe diameter at low ice concentrations near the ice edge, allowing for greater wave penetration further north. However, a significant issue with this combination is the wave-fracture algorithm in CICE6, which relies on wave energy. This algorithm faces convergence issues and does not function as intended, resulting in the warning: warning: step_wavefracture struggling to converge.

IC4M2 Wave Attenuation Model

Subsequently, the empirical wave attenuation model IC4M2 with the coefficient values from Meylan et al. 2014 was tested while maintaining the IS2 scattering model. In this setup, waves with small amounts of energy are being cut off, and the wave attenuation is too strong. Additionally, thin ice is being categorised as having the largest floe size, unlike in the IC3 model. This issue arises from the implementation of the Shen "pancake" theory in CICE6, where new ice consolidates into a continuous sheet if the wave energy in the cell is below a certain threshold. The attached output demonstrates that significant wave heights are being cut off within the ice cover.

The spectral components plot indicates that higher frequency waves dissipate first, likely due to wave scattering. This is evident along a transect where the significant wave height decreases sharply from 10^−2 to 10^−15 around -59° S. In comparison, the expected behavior is a continuous decrease in significant wave height as waves encounter ice, which is not observed in the current simulation.

Spectral components plot:

Test experiments with IC4M2 model alone, without scattering:

To address these issues, test experiments using the IC4M2 wave attenuation model alone were conducted, similar to the wave attenuation and scattering model choice of E3SM. However, this did not yield significant improvements. Transect analysis revealed a rapid decrease in wave energy, represented by normalized significant wave height, as soon as waves encounter sea ice.

Thanks @NoahDay and Luke, for your valuable inputs, analysis, and particularly for the plots of transects and spectral components. Thanks @aekiss, @ofa001 and @anton-seaice for the valuable insights.

[ezhilsabareesh8]

Recent developments discussed after the CICE consortium workshop:

1. Most groups (CESM, UFS, E3SM) are working towards the Meylan 2021 method, referred to as IC4M8 (but different from IC4M8 in the main WW3 repository). Some of the significant works in improving the wave attenuation scheme in WW3 are Cooper et al. 2022 and Meylan et al. 2021

2. Numerics issue highlighted by Cecilia (CC):

   • An issue for it is available here: Improve numerical solution for waves in sea ice NOAA-EMC/WW3#738
   • The fix for the numerics is in Dave’s branch and also Meylan et al. (2021) formulation (i.e., IC4M8). This will be integrated into the main WW3 fork over the next few months. This fix applies to all sea ice attenuation schemes.

3. Anton noted that the switch choice only impacts sea-ice attenuation. He also recommended to validate the open ocean first, and confirm that the model choices for the other “Source terms” are reasonable.

Response:

I did a initial wave buoy data validation for the open ocean, and it showed a deviation from the measurements for higher wave heights. However, this comparison is not sufficient as we also need to validate with satellite data. Validating with wave buoy data is a complex process that requires additional processing, as highlighted by Stefan.

Below shows the wave buoy and model comparison for 46059 (LLNR 382) - WEST CALIFORNIA, Year 1999.

Moreover, the choice of other source terms was consulted with Stefan, and a reasonable parameterization has been made. Currently, we are using ST6 with the recommended parameters by Stefan. Additionally, other source terms like the Unresolved Obstacles Source Term (S_uo) need to be considered. This term accounts for unresolved bathymetric and coastal features, such as cliffs, shoals, and small islands, which are major sources of local error in spectral wave models. Their dissipating effects can accumulate over long distances, and neglecting them can compromise the simulation skill over large portions of the domain.

Next steps:

• Set up the code for IC4M8 so that it can be tested and fix numerics issue.
• Validate the open ocean with satellite data.
• Consider including other source terms such as unresolved obstacles source term S_uo.
• Currently, we are using the first order propagation scheme since only this scheme is suitable for tripolar grids. We may either have to switch to an unstructured grid or use the rotated lat-lon grid to move the north pole to Greenland as done in E3SM to use higher order propogation schemes. Anything to add here @aekiss?

[dabail10]

I will also add here, that there is a heat, freshwater, salt conservation issue in CICE with the current FSD scheme. I am working with Lettie Roach on this. This leads to conservation check failures in add_new_ice. (If conserv_check = .true.)

[aekiss]

Related: https://forum.access-hive.org.au/t/access-om3-evaluation/1462/13