Synchronizing of fluxes especially when a point goes from ice free to ice covered.

[dabail10]

In ice_flux.F90 (and icedrv_flux.F90 in icepack) we set flwout to zero so that it can be accumulated in merge_fluxes. However, there is a timestepping bug here such that if you have a point that goes from ice free to ice covered in a single timestep, then flwout will not be computed in merge_fluxes (because aicen_init is all zero) and stays zero. However, the coupler fluxes are passed based on aice, so for a single timestep flwout is passed to the coupler with a zero value. This causes a crash in CAM and likely is a problem in other atmospheric models. I'm not sure why this was never a problem before.

The solution is to remove the line:

flwout = c0

in ice_flux.F90 and icedrv_flux.F90 and only do this just before the call to merge_fluxes based on aice_init. Working on it now.

[dabail10]

I have a solution. Let me know if you think this is ok. Basically, after merge_fluxes where lwout is computed, I make sure it never goes above the open ocean value (negative).

      ! Make sure flwout is not zero.

      TsfK = sst + Tffresh
      flwout = min(flwout,-stefan_boltzmann * TsfK**4)

[anton-seaice]

@kieranricardo - this might be relevant to you ?

[dabail10]

@peverwhee

[dabail10]

Actually I have this backwards. The reasonable range of temperatures in sea ice covered regions would be:

-80C to 0C or 193.15 to 273.15K.

So, this means a range of 79 W/m^2 to 316 W/m2. So, this should be:

flwout = min(flwout, -79)

since flwout is negative.

[dabail10]

Actually a better option is to add a check in scale_fluxes something like:

if (aice > puny .and. abs(flwout) < puny) then
   flwout = flwout_ocn
endif

This makes sure to catch the one special case of a point going from ice-free to ice-covered and the flwout is zero.

Dave

[dabail10]

I have a PR in for this now, but it occurs to me that there may be similar problems with fswabs, fsens, flat ... while a zero value is valid in these fields, we should technically be setting points that go from ice free to ice covered to have an open ocean value for these. For example, you can look in init_coupler_flux for some of these initial settings for different seasons.

[anton-seaice]

This doesn't solve you problem Dave, but I am slightly confused about the current behavior here.

In scale_fluxes, if there is some ice (aice>0) then the exported longwave is the radiation from the ice only and doesn't include the ocean. However if there is no ice, then the exported longwave is the radiation from the ocean.

These seems inconsistent. If the ice area is small, then the exported longwave is small (as it doesn't contain ocean longwave). If the ice area is zero, the longwave history output is large!

This is in scale_fluxes, so I think it only impacts the history output ? Why isn't

CICE/cicecore/cicedyn/general/ice_flux.F90

Lines 1263 to 1267 in 6aa677e

	fsens (i,j) = c0
	flat (i,j) = c0
	fswabs (i,j) = c0
	flwout (i,j) = -stefan_boltzmann *(Tf(i,j) + Tffresh)**4
	! to make upward longwave over ocean reasonable for history file

flwout just c0 ?

[dabail10]

This was super confusing! Took me a long time to figure this out. Let me try to explain.

1. In init_coupler_flux we have:

      flwout  (:,:,:) = -stefan_boltzmann*Tffresh**4
                        ! in case atm model diagnoses Tsfc from flwout

So, for all points this is set to the lwout equivalent to 0C (273.15K)

2. However in init_flux_atm:

flwout(:,:,:) = c0

This is everywhere so that flwout is computed as flwout = flwout + aicen(n)*flwoutn(n) for n=1,ncat in merge_fluxes. Since we are accumulating over subgridscale categories, this has to be initialized to zero.

3. The call to merge_fluxes is only done for aicen_init(n) > puny. So, flwout is only computed when there is ice at the beginning of the step.

4. However, it is possible that ice grows during the step and in scale_fluxes we do:

flwout = flwout / aice

where aice is the ice area at the end of the step. So, if flwout was not computed at the beginning of the step, then flwout is zero which is not physically possible.

Does this make sense?

[apcraig]

It seems like there might be two different issues.

1. What's on the history file. We can define this however we like. Ice free areas can have flwout whatever. This assumes the history file output is not being used to validate conservation. (see point 2).

2. What's being coupled. So, this seems important for conservation. Shouldn't the implementation be based entirely on conservation. The lw flux passed out at the ice surface should be consistent with the lw flux "used" in the ice model. If the ice fraction is zero, the lw flux value plays no role. But ultimately, the lw flux should be conserved. When ice goes from/to zero fraction, what lw flux is used in the ice model and are we conserving when we pass to the coupler? If lw flux is zero when aice > 0, and we fix that by setting lw flux to the ocean temperature lw flux, that seems wrong. If lw flux is non-zero when aice=0 (ice fraction -> 0), that seems wrong too. I wondered whether maybe the lags in the system were taking into account missing lw flux, but I don't think they can be. It could be that we have a conservation error that is reduced with this fix, but it still seems like we're not conserving. Is there an implementation that eliminates all conservation errors?

[dabail10]

This is not really a conservation issue. What the atmosphere is receiving is:

flwout_atm = aice * flwout_ice + (1-aice) * flwout_ocn

Because aice is non-zero, this means the atmosphere can receive a zero flwout from the ice meaning a temperature of zero degree Kelvin! This is completely unphysical. You might think that aice is very small as it is the growth over one timestep, but in CESM we are seeing the ice area go from zero to 0.8 in a single timestep over all of Hudson Bay!

[apcraig]

It is a conservation error as you describe it. If aice is 0.0001 when ice first grows, one could argue the error is small. If aice is 0.8, that's interesting and significant, but in either case, it seems like the system is not conserving. Changing the flw from zero to something associated with the ocean temperature may be reducing the error, but it's not conserving the heat. The flw*aice that the coupler receives should match exactly the long wave flux "used" in sea ice in order to conserve. Going to and from aice=0 may be a problem that we cannot overcome give the current implementation constraints, so minimizing the error may be the best option. But I'm pretty sure there is a conservation error.

I thought we were coupling the ice model each ice timestep to avoid these kinds of problems, in part by avoiding having to average (or similar) an ice fraction when some timesteps have zero ice fraction, but for other reasons related to ice state switching to/from zero ice fraction. I'm a little surprised this error has been in the implementation for so long. It must be quite small most of the time and especially when assessing conservation over longer periods.

[dabail10]

I agree this has been there a long time! They only recently added a check to the atmospheric radiation for this case. You do raise a good point in that the flwout should be computed based on the first temperature of the ice that forms in that timestep. Here we could use the Tsfc that is computed in the thermodynamics.

[NickSzapiro-NOAA]

When ice goes none->some or some->none, when is this formulated to occur in the timestep? Like, is it at the end of the timestep as a step function?

[dabail10]

When ice goes none->some or some->none, when is this formulated to occur in the timestep? Like, is it at the end of the timestep as a step function?

So, the idea is that if you have no ice, this doesn't matter for the sea ice model or the atmosphere model. The ocean upward longwave is computed separately. So, it is arbitrary what one sets the open ocean values to. This is for history purposes only. Then the thermodynamics (step_therm1) is called to compute the surface fluxes only when there is ice at the beginning of the timestep. These fluxes are computed over the ITD categories and then a weighted sum to get the grid cell value:

flwout = flwout + aicen_init(n)*flwoutn(n), n=1,ncat

where aicen_init is the category fraction at the beginning of the step. So, the fluxes (including flwout) must be zeroed out before we do this. However, the thermodynamics and then dynamics could create ice in the grid cell. The flux then is sent to the coupler / for history as:

flwout_cell = flwout / aice

where aice is the ice fraction at the end of the step. So, if aice_init is 0, but aice > 0, then you are never computing flwout and sending 0 to the coupler. This should actually be set to the open ocean value or computed using the updated surface temperature.

Now if we have the situation that aice_init > 0, but aice = 0, this is not a problem, because then the flwout is set to the open ocean value in scale_fluxes and is not used by the atmosphere model. So, there is a step function in the first case, but not really in the second except for history output.

Dave

[eclare108213]

I agree with @apcraig, conservation should guide what the correction should be. Moreover, the conservation check should not be for a particular point in time (beginning/end of time step), but rather across an entire coupling interval and including all components' LW fluxes. Is that what CESM is now doing, @dabail10?

The sea ice model was developed in isolation and so it conserves internally (or should), but there are time lags in the coupling that mean something is always going to be inconsistent in the coupled system, unless there's an iterated time barrier when they all resynchronize. E.g. we have what could be considered the first step of an iteration: we calculate albedo at both the beginning and the end of the ice model time step (also the coupling interval), first to have updated SW fluxes for the thermo calculation and then to ensure that the albedo sent to the coupler is consistent with the updated ice thickness distribution at the end of the time step.

Perhaps the LW value sent to the coupler could be an average between the beginning of the time step and the end, assuming that the ice area grows or shrinks uniformly through the time step and the surface temperature stays at the value computed in the thermo; then using 0.5*(aice - aice_init) as the scaling factor might be a better approximation. (This could be refined further by assuming that the temperature changes linearly and estimating LW throughout the timestep based on that temperature evolution.)

Regardless, I'd be interested in seeing a full conservation analysis to understand what's getting passed when in the coupled system(s), and where the inconsistencies are.

[dabail10]

There is already conservation across the CESM for lw fluxes. The lw that the atmosphere sees is:

flwup_atm = fice * flwup_ice + focn * flwup_ocn + flnd * flwup_lnd

This is the lower boundary condition for the longwave computation in the atmosphere. They do invert the LW up to check the surface temperature. This was where the issue was discovered. The temperature was far too low because the flwup_ice term was zero.

Let's think about the coupling fluxes then.

Step 1 (t = 0)

fice = 0., focn = 1.

flwup_atm = flwup_ocn

Step 2 (t = 30min)

fice = 0.6, focn = 0.4

flwup_atm = 0.6*flwup_ice + 0.4 * flwup_ocn

In this case currently, flwup_ice = 0 which is unphysical. We have a couple of choices.

• flwup_ice = flwup_ocn and then flwup_atm is exactly the same as it was at step one.

• flwup_ice = flwup_ice*

where flwup_ice* could be based on the freezing point, the current sst (same as flwup_ocn), or an ice surface temperature of 0C (273.15K), or some average of these.

Dave

[dabail10]

I guess, I am not thinking of this as a conservation issue. One expects d(flwout) / dt to change from step to step. The atmosphere will use whatever it is given. I guess the "smoothest" transition is to go from SST to Tsfc on the ice?

[eclare108213]

I do think there's a conservation issue. Consider
flwup_atm = 0.6*flwup_ice + 0.4 * flwup_ocn
I believe this is calculated as
flwup_atm = sum(aice(n)) * (sum(aice_init(n)*flwup_ice(n)) / sum(aice(n)) + (1-sum(aice(n))) * flwup_ocn
for n=1,ncat (please correct me if I'm wrong!).
The first sum is what's multiplied by the coupler. The second sum is what we compute in the thermo and then (consistently) aggregate across the initial categories. The third sum is our scaling so that when the coupler multiplies by the first sum, the flux is what we calculated over the ice. But then what the atmosphere is really seeing is
flwup_atm = sum(aice_init(n)*flwup_ice(n)) + (1-sum(aice(n))) * flwup_ocn.
Since the area fractions are inconsistent, there's a conservation problem. I think this should be
flwup_atm = sum(aice_init(n)*flwup_ice(n)) + (1-sum(aice_init(n))) * flwup_ocn.

I reached this conclusion through a derivation from a completely different direction that considered three area fractions, the ice area that doesn't change during the time step (if any), the ocean area that doesn't change (if any), and a transition area where either ice becomes open water or open water becomes ice. What I came up with is that the coupler should be multiplying the component fluxes by the initial surface area fractions (i.e. at the beginning of the time step) rather than the final fractions, for both ice and ocean. I can type up the derivation if you want to see it.

An alternative would be to calculate
flwup_atm = sum(aice(n)) * (sum(aice(n)*flwup_ice(n)) / sum(aice(n)) + (1-sum(aice(n))) * flwup_ocn(n)
(or simplify it by canceling aice) and maybe that's what the coupler / atmosphere are expecting. This would be consistent for ice and ocean but still wouldn't fix the problem in this issue, when aice_init is zero and aice is not. Using the initial area fractions in the coupler is the right way to handle it, I think.

Let me know if I'm doing something wrong here...

[dabail10]

Well, I should really have @duvivier and @marikaholland chime in here. The reason we divide by aice (and not aice_init) is so that the computation of:

flwup_atm = aice * lwup_ice + (1-aice)*lwup_ocn

where lwup_ice = sum(aicen_init(n)*flwupn_ice(n)) / aice.

where aice is cancelled out in the first term as you pointed out. The ice fraction aice is the only fraction the atmosphere sees and is the only thing at a given step that can be used in the coupler. The flwup_ice is computed at the beginning of the timestep in the sea ice before the ice state changes. It is the surface energy budget that drives the thermodynamic changes to aice. We always have this timestep offset with the fluxes and the ice state at the end of the step. This is why the some of the fluxes used to get passed in the middle of the timestep before the ice state changed. However, all the coupling is done at the end of the timestep now.

[NickSzapiro-NOAA]

If the longwave the atmosphere receives needs to be lost by the ice (at the beginning of each ice timestep), it seems like aice_init is the relevant coupling ice_fraction (for longwave)

Doesn't aice*fill_value send energy to atmosphere without any energy balance from ice or ocean?

[eclare108213]

I see a mistake in my derivation, so let me think about this some more.
Regardless, this is still true in the coupled system at the moment, since the ice model divides the aggregated ice flux by aice to scale it up to the grid cell and the coupler then multiplies by aice:
flwup_atm = sum(aice_init(n)*flwup_ice(n)) + (1-sum(aice(n)))*flwup_ocn.
Do you agree that's wrong? Or at least problematic?

[apcraig]

From my perspective, conservation means that the energy needs to be exactly consistent across components. So that means Joules (or Watts) needs to be conserved. But we are coupling fluxes, so there is an area involved. There are an infinite number of ways to conserve by varying fraction and flux at any given coupling period. But there is only one value of Joules (or Watts) that will conserve. If the sea ice model uses or releases some energy and the ice fraction is zero, then there is no way to conserve from what I can tell. If you are computing a flux and you know the energy used, there is only one fraction that will conserve. If you know the ice fraction and the energy, then there is only one value of flux that will conserve. This isn't a case of "guessing" a flux associated with a fraction.

Having said all that, if the lw flux up from the ice is simply associated with the surface temperature and the value of that flux is just an ice diagnostic, then maybe we can set it to whatever we want understanding 0 degK is a particularly bad choice. But it's also a little hard for me to believe that the "right" surface temperature is arbitrary. I don't think the problem is how the coupler merges the ice and ocean fluxes, the problem is what do we set ice fraction and ice flux to in these edge cases.

Finally, I think it's worth taking a slightly bigger step out considering all fluxes. I've always worried a little about whether the time varying ice fraction is really being taken into account properly with respect to conservation, especially on these edge cases where the ice was or becomes zero. In that case, there is no area to apply the fluxes, so it doesn't matter what the flux value is. But if there is some energy that needs to be conserved across the system, how is that handled?

[eclare108213]

Setting Tsfc to SST for computing lwup seems like a reasonable stop-gap for the particular case prompting this issue, but it doesn't solve the problem for other fluxes. I think what's needed is to recompute all of the thermo fluxes at the end of the time step, to match the final ITD. That would take care of the aice_init=0 issue and the difference in aice time levels used for the atmosphere fluxes at the end of the time step. There might still be weird inconsistencies regarding time level through the time step, though, and figuring out what the total/net flux to the coupler should be might not be trivial. It might just move the time lag error from the coupler into the ice model...

I do agree with Tony re the need for a deeper review of the coupling in light of conservation requirements.
Do any of the other centers use a different process order during the time step?

An interim step could be to just recompute the fluxes only for the new ice area immediately after the call to add_new_ice, and include it in merge_fluxes. That would help with the aice_init=0 issue for frazil growth but not for ice transported into previously empty cells, and it would be less accurate than a full thermo recompute because the fluxes over the existing thickness categories wouldn't have been updated using their new thicknesses. Even so, it might help with some of our other ice edge issues. I haven't really thought through this, so feel free to shoot it down.

[dabail10]

All good thoughts. I do think it is impossible to get rid of the time lags here. I know GFDL does something where it splits the "fast" and "slow" thermodynamic processes. The fast ones are computed in the atmospheric boundary layer and the slow ones are computed in the oceanic boundary layer. Similarly for the dynamics. Hadley Centre does a similar thing with the surface fluxes. Maybe @proteanplanet has some thoughts here as well.

Think of this as a timestep with a forward Euler discretization:

Flwout_ice(t+1) = Flwout_ice(t) - d(stefan_boltzmann*Tsfc**4)/dt

where

d/dt of the Tsfc term is:

4*stefan_boltzmann*(Tsfc(t+1)**3 - Tsfc(t)**3) * (Tsfc(t+1)-Tsfc(t)) / dt

So, this is a case where Tsfc(t) is the open ocean SST and Tsfc(t+1) has not been computed yet as ice formed at the end of the step. The thermodynamics (step_therm1) is only called when aicen_init > puny. So, I think we have to assume that the surface temperature of this newly formed ice is the same as the open ocean. Hence, Flwout_ice(t+1) = Flwout_ice(t). This is all in the sea ice model. The atmosphere uses whatever it is given and passes back a new incoming longwave for the next step.

[proteanplanet]

Flagging @erinethomas who is leading our collab with Lili Manzo at UCI to couple emissivity similar to shortwave albedo.

[dabail10]

If the longwave the atmosphere receives needs to be lost by the ice (at the beginning of each ice timestep), it seems like aice_init is the relevant coupling ice_fraction (for longwave)

Doesn't aice*fill_value send energy to atmosphere without any energy balance from ice or ocean?

I'm not sure what you are asking here?

[NickSzapiro-NOAA]

Sorry. So I agree that the atmosphere can use/conserve whatever flux boundary condition the ice provides.

But doesn't this (aice x flwup_ice*) term never enter the "top surface flux balance"? Then mass+energy conservation is broken in the coupled system as it's broken within the ice model

[dabail10]

I think we can agree to disagree on this. Regardless, this is a bug to have zero longwave out being passed to the coupler. I think the most consistent thing is to set the temperature (and hence longwave) to open ocean values in the polar regions. The other fluxes such as sensible and latent heat are a larger discussion.

[eclare108213]

But doesn't this (aice x flwup_ice*) term never enter the "top surface flux balance"?

@NickSzapiro-NOAA, it is included (the sigma *T^4 term in the figure above). I think the confusing part is that the surface temperature / vertical thermo calculation is strictly at a point, and that point is over ice in the grid cell. Using the ice thickness distribution, we do the thermo calculation over each of the thickness categories, as if they cover the whole grid cell, then aggregate the results into one set of fluxes from the ice, using the area fractions of each thickness category to weight them. The ocean calculates what its fluxes would be over open water, and then the coupler merges the ice and ocean parts using the ice and ocean fractions. Is that what you are asking about?

There are actually 2 related questions now, first (Dave's) what to do when the initial ice fraction is zero but ice appears during the time step, since the ice model's flux calculation is based on the initial ice fraction, and second (mine) what area fractions the coupler is using to merge the ice and ocean fluxes. We can stick to Dave's question in this issue, but I think the second question needs to be looked into also.

[dabail10]

I have added a fix for the flwout = 0 problem. A more important issue is the timestepping of the fluxes versus the ice state.

1. In init_flux_atm, flwout, fswabs, fsens, flat, etc are are initialized to zero.

• this is because they are accumulated in merge_fluxes to get the "gridcell category weighted fluxes".

2. When there is ice a the beginning of a step, then flwoutn, fsensn, flatn, are computed in step_therm1. If not, then these are all zero.

3. At the end of the step, the radiation code is called and the albedos and fswabsn are updated based on the ice state at the end of the step.

4. So, this means when a point goes from ice free to ice covered, the flwout, fsens, and flat are not up to date with the latest ice state.

5. This is all complicated by the fact (in CESM) that the albedos and fluxes are sent half hourly to the coupler and yet the atmospheric radiation is computed hourly. That means the sea ice receives the same incoming radiation for two time steps. This might be impossible to synchronize.

[duvivier]

@dabail10 Sorry it took me a while to comment (blame AGU).

We have talked about this and my understanding is that this inconsistency comes because of the use of aice_init in the initial fluxes and then aice_final in the grid cell averaging. Basically, those two values should be consistent. I think we all agree with that. But as @apcraig and @eclare108213 mention, the conservation issue is that then CICE passes flwout=0, so no energy loss at all from the ice (because there is no ice). But of course over the course of the timestep when aice_init goes to aice_final there actually is energy lost as soon as the ice forms.

Honestly, I odn't know what the best solution is, but it is causing issues with CESM at least for the coupling. My inclination is to do something like @eclare108213 suggested:
flwup_atm = sum(aice_init(n)*flwup_ice(n)) + (1-sum(aice(n)))*flwup_ocn
Where the flux passed uses a consistent aice value as that used to compute the flux. From what @dabail10 says, I believe this would apply to longwave up, sensible, and latent heat fluxes. Maybe shortwave out (I wasn't clear if that was included in @dabail10 message, above).

@dabail10, I don't see a PR with your suggested change for the fix. Can you point me to the update?

[dabail10]

The LW fix was brought in with PR #999. Thanks for your thoughts @duvivier. I have changed this to be a larger issue now. See my comment above for the latest.

[duvivier]

Ah, I missed it because I was looking for open PRs. Thanks!