Sea ice temperature advection for AWICM3

[JanStreffing]

@dsidoren pointed out that the way we have implemented sea ice temperature advection for the coupled AWICM3 setup does not look right to him. Rather than advection the heat energy we are advection the temperature directly. In case of convergence your temperatures are added up, in case of divergence they are reduced.

I would like to invite comments from @lzampier @qiangclimate

[lzampier]

I understand @dsidoren concerns. However, we do not advect the temperature, but temperature multiplied by the ice-area fraction (L189-193; ice_setup_step.F90 ). After advection, we convert the advected tracer back to temperature by dividing by sea ice concentration (L200-204). This is very similar to what happens for the sea-ice thickness, which is also not advected, to avoid the same issue. The quantity that you advect is the sea-ice volume per unit area (m_ice), which is thickness times ice-area fraction.

Sorry but today I do not have time for a more exhaustive answer... I will write a more detailed comment one of the next days. Or we can have a quick meeting to discuss this.

[qiangclimate]

Instead of using concentration, sea ice volume (mean thickness in the model) should be used. That is, advecting T*h.

[qiangclimate]

one can understand the issue in this way: assuming two grid cells with the same T, both with 100% sea ice and with the same ice thickness. If wind push ice from one cell to the other, the second cell will still have 100% sea ice. But the ice volume (mean thickness over the cell) is doubled in the second cell. If we advect T*h, its value will be doubled in the second cell, devide it with doubled h, T remains the same.

[qiangclimate]

following the talk today: if the temperature referred to on this page is "surface" (skin) temperature, I would suggest Not to advect it, but rather do vertical physics to get it. All other properties that are advected, including energey related, will be reflected in skin temperature when doing the ice thermodynamics after advection. (or, if we are dealing with a mean T in a layer with some thickness h, then it can be advected. But in this case, it is T*h that should be advected.)

[dsidoren]

one more example: skin temperature does not change due to ice divergence. if to advect it ... it will! if at all we shall advect only the total heat but 0 layer is assumed here and the advection is not necessary. (sorry I saw the discussion only now)

[dsidoren]

https://github.com/FESOM/fesom2/blob/b730dc8385132064017db0da1f1117a1fbbbac97/src/ice_thermo_cpl.F90#L471 shall not this be zicefl=con*(t-ctfreez)/zsniced ?

[JanStreffing]

https://github.com/FESOM/fesom2/blob/b730dc8385132064017db0da1f1117a1fbbbac97/src/ice_thermo_cpl.F90#L471

shall not this be zicefl=con*(t-ctfreez)/zsniced ?

If the line were wrong we would overestimate conductive heat flux an order of magnitude. Melting and freezing rates and the seasonal thickness cycle do not indicate this. Comparing to mo_surface_ice.f90 line 1160, it's implemented as in AWICM1&2.

I think the dependency on the temperature differential comes in here: https://github.com/FESOM/fesom2/blob/b730dc8385132064017db0da1f1117a1fbbbac97/src/ice_thermo_cpl.F90#L474-L476

The comment in line 471 is obviously not correct though. This is not actual conductive heat flux though sea ice. That would indeed require the formula you wrote. In mo_surface_ice.f90 this is done in line 1170 with the updated surface temperature.

[JanStreffing]

following the talk today: if the temperature referred to on this page is "surface" (skin) temperature, I would suggest Not to advect it, but rather do vertical physics to get it. All other properties that are advected, including energey related, will be reflected in skin temperature when doing the ice thermodynamics after advection. (or, if we are dealing with a mean T in a layer with some thickness h, then it can be advected. But in this case, it is T*h that should be advected.)

Hello Qiang, thank you for the feedback. If I were to implement this, this would entrail changing from a_ice to m_ice in https://github.com/FESOM/fesom2/blob/efae0221ecfc81650ec1e38322c089e685378ab2/src/ice_setup_step.F90#L187 https://github.com/FESOM/fesom2/blob/efae0221ecfc81650ec1e38322c089e685378ab2/src/ice_setup_step.F90#L198 and limiting m_ice to 0,05 meters max?

[qiangclimate]

hi Jan, if what you deal with here is skin temperature, it should not be advected. You get it just from ice thermodynamics locally (and do not advect it further!) in the case of 0 layer sea ice as in fesim. If sea ice has heat capacity as in ice pack, then energy or enthalpy will be advected. But skin temperature is still not to be advected in this case. Computing it from local thermodynamics will automatically contain the effect of advection of energy. I had a quick chat with Dima today and we share the same opinion.

[dsidoren]

https://github.com/FESOM/fesom2/blob/b730dc8385132064017db0da1f1117a1fbbbac97/src/ice_thermo_cpl.F90#L471

shall not this be zicefl=con*(t-ctfreez)/zsniced ?

If the line were wrong we would overestimate conductive heat flux an order of magnitude. Melting and freezing rates and the seasonal thickness cycle do not indicate this. Comparing to mo_surface_ice.f90 line 1160, it's implemented as in AWICM1&2.

I think the dependency on the temperature differential comes in here:

https://github.com/FESOM/fesom2/blob/b730dc8385132064017db0da1f1117a1fbbbac97/src/ice_thermo_cpl.F90#L474-L476

The comment in line 471 is obviously not correct though. This is not actual conductive heat flux though sea ice. That would indeed require the formula you wrote. In mo_surface_ice.f90 this is done in line 1170 with the updated surface temperature.

I believe there is everything correct in https://github.com/FESOM/fesom2/blob/b730dc8385132064017db0da1f1117a1fbbbac97/src/ice_thermo_cpl.F90#L437 although I cannot fully trace the approximations behind. I know you and Qiang can :).

we generally solve for skin T through the heat balance equation omitting dQ/dt and also apply several iterations to reach the convergence for skin T. At the same time we obtain all the fluxes we need. In this code we have the flux already but do not have neither atmospheric temperature nor specific humidity to complete the heat balance equation. A bit confusing.

[lzampier]

Here are my considerations about the advection issue, and sorry for the extreme delay in contributing to the discussion. To my understanding, the ECHAM6 implementation of the surface temperature calculation, which has been implemented in AWICM3, works as follows. The sea-ice temperature that we compute is representative of a thin, but finite, layer of ice and snow (if present) at the sea-ice/atmosphere interface. This thin superficial layer is the one responding most rapidly to changes in atmospheric temperatures, representing a compromise between a 0-layer and a more rigorous formulation of the sea-ice thermodynamics. In this respect, I agree with Dima: the choice is for sure suboptimal, but that is what we have, and we should make it work for the time being. Anyway, the thickness of the layer is defined as follows: 5cm of sea ice plus the snow thickness converted to snow-ice-equivalent. As the heat capacity of snow is much lower than that of sea-ice, the snow thickness component is small and furthermore homogeneous over multiple grid cells for typical climate simulations. The 5cm of sea ice are homogeneous by definition and, therefore, the thickness of this surface layer varies only slightly (depending on the snow cover) from one grid cell to another.

As this is not a simple skin temperature determined diagnostically at every time step, but a prognostic quantity representative of an energetic state, it needs advection. The correct quantity to advect would be the following: T_layer h_layer a_ice. Note that here I am using the real thickness of the layer and not a volume per unit area as in FESOM2. Hence, the multiplication by a_ice. As the variations of h_layer are almost negligible from one grid cell to another, as I argued before, I simplified the advection term to T_layer*a_ice, and that is currently implemented in the model. It would possible to be rigorous and advect also h_layer, it is not difficult but I don't expect this to have a major impact on the results.

Note that this advection is neglected altogether in ECHAM6 because the atmospheric model did not have any information about the ice velocity. Nevertheless, in that context, this represented a good approximation considering the low resolution of the model and the short timestep. With FESOM2, we need to be more careful as the model resolution can be substantially higher.

[qiangclimate]

Hi Lorenzo, to understand the error of the code, I repeat the example given above. Consider we have only two same grid cells with all variable the same. The same Tskin, a_ice, h_ice. And set a_ice=1 (100%). Now, blow ice from one grid cell to the other. What will happen in reality: sea ice ridges and h_ice is doubled; Tskin (almost) the same as before; a_ice is the same as before (100%). What the code as now will do: Tskin will be doubled! (a negative Tskin, if doubled, it is rendered to be too cold)

Independent on this issue, the way how to "thermodynamically" compute Tskin is discussed above as I can see. Doesnt matter how good or bad the way echam makes the approximation, it is just a way to approximate Tskin and has nothing to do with heat capacity advection. They did not choose to neglect the advection, but rather they cannot do advection, for the reason explained above.

[lzampier]

Yes, Qiang, I understand your point. Skin temperature, as you said, should not be advected in association with the thickness. I repeat however that we are not dealing with a skin temperature, nor with a temperature representative of the whole sea ice column. The temperature that we have in AWICM3 refers to a specific surface layer of ice of thickness h_layer = ~ 5cm, and needs to be advected accordingly.

I suggest we have a quick meeting to clarify this issue one of the next days, as it might be much quicker.

[qiangclimate]

on the contrary, T "should be" advected in association with the thickness, not "should not be". As T*h should be advected, while h is assumed to be fixed, then there is no sense to advect T. The Key point is that sea ice drift is NOT divergence free. If sea ice were non-compressible, as water, then it would be fine. Yes, we could have a quick chat on this :))

[dsidoren]

Here are my considerations about the advection issue, and sorry for the extreme delay in contributing to the discussion. To my understanding, the ECHAM6 implementation of the surface temperature calculation, which has been implemented in AWICM3, works as follows. The sea-ice temperature that we compute is representative of a thin, but finite, layer of ice and snow (if present) at the sea-ice/atmosphere interface. This thin superficial layer is the one responding most rapidly to changes in atmospheric temperatures, representing a compromise between a 0-layer and a more rigorous formulation of the sea-ice thermodynamics. In this respect, I agree with Dima: the choice is for sure suboptimal, but that is what we have, and we should make it work for the time being. Anyway, the thickness of the layer is defined as follows: 5cm of sea ice plus the snow thickness converted to snow-ice-equivalent. As the heat capacity of snow is much lower than that of sea-ice, the snow thickness component is small and furthermore homogeneous over multiple grid cells for typical climate simulations. The 5cm of sea ice are homogeneous by definition and, therefore, the thickness of this surface layer varies only slightly (depending on the snow cover) from one grid cell to another.

As this is not a simple skin temperature determined diagnostically at every time step, but a prognostic quantity representative of an energetic state, it needs advection. The correct quantity to advect would be the following: T_layer h_layer a_ice. Note that here I am using the real thickness of the layer and not a volume per unit area as in FESOM2. Hence, the multiplication by a_ice. As the variations of h_layer are almost negligible from one grid cell to another, as I argued before, I simplified the advection term to T_layer*a_ice, and that is currently implemented in the model. It would possible to be rigorous and advect also h_layer, it is not difficult but I don't expect this to have a major impact on the results.

Note that this advection is neglected altogether in ECHAM6 because the atmospheric model did not have any information about the ice velocity. Nevertheless, in that context, this represented a good approximation considering the low resolution of the model and the short timestep. With FESOM2, we need to be more careful as the model resolution can be substantially higher.

Hi Lorenzo, my concern is that if we choose to advect we must also advect this thin superficial layer (h_thin^{n}a_ice^{n}) together with (T^{n}h_thin^{n}*a_ice^{n}) and then extract the T^{n+1} otherwise it is purely diagnostic and the advection is not permitted.

[JanStreffing]

We had a short meeting on this issue. Protocol:

When we advect with convergence we divide by the sea ice concentration before the cuttoff. Therefore the overflown advected quantity is divided by the overflown concentration and we obtain the correct quantity in the target cell. The principle concern is resolved.

We further discussed about potaltial improvements to this very simple scheme.

1. We could iterate the prognostic equation
2. We could couple specific humidity as well as air temperature and use the more sophisticated therm_ice from standalone fesom2:

https://github.com/FESOM/fesom2/blob/b730dc8385132064017db0da1f1117a1fbbbac97/src/ice_thermo_oce.F90#L186-L213

The improvements may be made when ICEPACK coupled setups are created down the line.

I will leave this issue open a little while longer for final comments.

[JanStreffing]

Wolfgang Dorn reviewed part of our ice_thermo_cpl.F90 and commented:

As a sidenote: I find Lines 484-485: https://github.com/FESOM/fesom2/blob/5bfa6c9ba01d7947f288d4aec83ea2fb1a9b71b6/src/ice_thermo_cpl.F90#L479-L487 to be suspicious since we cannot heat or cool the entire snow layer in one timestep if the layer is thicker than a few centimeters. The term zcpdte can therefore be much too large, whereby the temperature changes can be too small. For a snow surface the heat capacity term should be smaller than for a pure ice surface. Here the opposite is the case.

I double checked with the AWI-CM1&2 routine and it seems to me like we adapted it 1-to-1 for AWI-CM3: https://gitlab.dkrz.de/modular_esm/echam6/-/blob/master/src/echam/mo_surface_ice.f90#L1162

That does not mean there is nothing wrong with it of course. I would like to hear your feedback.

[JanStreffing]

I plotted the local end of winter snow depths from /work/ollie/jstreffi/runtime/awicm-3.1/LED2/outdata/fesom

Is this enough in the SH? It's not clear if this is a cause or effect of having little sea ice in the SH.

[JanStreffing]

Make 0.05 m variable time step dependent.

[qiangclimate]

Even in Winton 3layer model, snow heat capacity is considered to be zero. So there is no reason to consider here. I would suggest the changes to the 3 lines as the following: zcpdte=hcapice/dt ! Energy required to change temperature of top ice "layer" [J/(sm²K)] ! zcprosn=rhowatcpsno/dt ! Specific Energy required to change temperature of 1m snow on ice [J/(sm³K)] ! zcpdte=zcpdt+zcprosnhsn ! Combined Energy required to change temperature of snow + 0.05m of upper ice

No idea whether this has any impact on the current issue. But doing this change is consistent with the overall model assumption.

[JanStreffing]

Hey Qiang, that's a good idea. It fits nicely with another test I'm running right now. Dmitry Sein suggested that the relevant ice thickness that should be considered (hardcoded to 0.05m at the moment) should depend on the timestep. Especially considering that FESOM2 with core2 mesh, as I used it for quick tuning runs, has much larger timesteps than ECHM6 T63 or T127 have.

So I'm running one test were I assume a much thicker layer of ice to be affected during one timestep. (0.25m). I will also run one test with 0.05m and without snow, then.

[qiangclimate]

forgot to mention, the following line (479) snicecond = con/consn*rhowat/rhosno ! equivalence fraction thickness of ice/snow should be snicecond = con/consn

[JanStreffing]

Hello @qiangclimate, I performed a few test runs.

• SWR4 is the default without any changes,
• spi1is a standalone FESOM2 run with the first ten years 2000-2009 and the last twenty years 1950-1969 forcing.
• ITH1 top layer of ice that is considered for heat capacity 0.05m -> 0.25m
• ITH2 no snow depth considered for heat capacity

As you can see the ITH1 test had a big impact on both hemispheres. Not exactly the way I wanted, but it's good to know for the future if we consider making this thickness timestep dependent.

The second test ITH2 however had virtually no impact. This surprised me. I would have thought that removing ~20cm water column mass equivalent of snow from the heat capacity of the SH ice surface layer might have an effect of similar magnitude as added 20 cm in ITH1 (albeit in the other direction). That's not the case. I have uploaded the changes I made for this test here: https://github.com/FESOM/fesom2/commit/b3ae997477ea1b058a380eef55196d090fd18e21. Pardon the misnomer of the branch.

[qiangclimate]

the results are not close to expected. Did tuning h0 help anyhow?

[JanStreffing]

Yes h0 tuning can get us quite far. With h0n=0.5 and h0s=1 we can reduce the ice volume bias by about 30-40% in both hemispheres. We are still looking for a root cause. Discussion is going in the internal issues. I've send you and invite.

https://github.com/AWI-CM3/project_management/issues/14 https://github.com/AWI-CM3/project_management/issues/18

[JanStreffing]

I will close this issue now. Should we revisit it later I will open a new one. Thank you for the code review and the fixes!