[MITgcm-devel] Unrealistic low SST with SEAICE_GROWTH_LEGACY undef
Ian Fenty
ifenty at MIT.EDU
Fri Feb 3 15:20:12 EST 2012
Sea ice people,
After looking at Martin's output, I (tentatively) verified that actual
thermodynamical ice processes are working fine and that nothing
thermodynamical is driving T < T_freeze. However, I did find two
non-thermodynamical extraneous code blocks which can cause temps to fall
below T_f. Allow me to point them out.
Culprit 1) When ice is very thin, we melt whatever ice and snow volume
is present using ocean enthalpy, regardless of whether any energy is
actually available to melt:
429 C 1.25) treat the case of very thin ice:
438 if (HEFF(I,J,bi,bj).LE.siEps) then
439 tmpscal2=-HEFF(I,J,bi,bj)
440 tmpscal3=-HSNOW(I,J,bi,bj)
441 TICE(I,J,bi,bj)=celsius2K
442 endif
443 HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal2
444 d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
445 HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
446 d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
Culprit 2) Remove HEFF when it is too large by melting, regardless of
whether any energy is actually available to melt:
1446 #ifdef SEAICE_CAP_HEFF
1447 C This is not energy conserving, but at least it conserves
fresh water
1448 tmpscal0 = -MAX(HEFF(I,J,bi,bj)-MAX_HEFF,0. _d 0)
1449 d_HEFFbyNeg(I,J) = d_HEFFbyNeg(I,J) + tmpscal0
1450 HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal0
1451 #endif /* SEAICE_CAP_HEFF */
Both of these code blocks can cause heat loss from the surface grid cell
when T <= T_freeze.
With respect to the first culprit, I argued earlier that it wasn't
necessary but I failed to predict that it could actually lead to SST
pathologies. With respect to the second, I don't know 1) where it comes
from, 2) why excessive HEFF is a thermodynamic issue (i.e. why this cap
is in seaice_growth), or 3) why treating the problem of excessive HEFF
is done via "not energy conserving" code.
I also don't know whether these parts are responsible for the Martin's
or anyone else's problem, although SEAICE_CAP_HEFF is defined in Gael's
SEAICE_OPTIONS.h that he sent to the list.
Unless a cogent argument can be made for keeping the very thin ice
"pathology", I suggest that we remove it or at least make it an option
that must be explicitly specified. For both the thin ice removal and
the HEFF capping, I suggest that if either is done, we at least we do it
without removing imaginary ocean heat. If HEFF capping is related to
some issue of dynamics, then I suggest that it be done elsewhere
(perhaps a routine which prepares the ice field for the dynamics
solver). The more we dump code which isn't related to thermodynamic
growth and the expansion and contraction of ice cover into
seaice_growth, the harder and harder it becomes to read and debug.
-Ian
On 2/3/2012 7:31 AM, Martin Losch wrote:
> Hi Ian,
>
> I put 31 figures with the diagnosed fluxes (all in W/m^2, a modified version of your heat_flux.m that I used to make these plots) here:
> http://mitgcm.org/~mlosch/forIan.tar (with full range color scale)
> http://mitgcm.org/~mlosch/forIanFixed.tar (with fixed range color scale)
>
> I have a hard time interpreting the plots:
> As long as the ice cover is more or less closed, r_QbyATM_open is basically QSW_SURF, but once it is open, e.g starting around day 10, the SSTair-SSTocean difference causes a large heatflux into the ocean that should increase the (super-cool) SSTocean. Instead, in some places the surface temperature decreases even further, see grid point (61,148)
>
> d_HEFFbyOCNonICE is always negative, meaning a loss of heat from ocean to atmosphere (always melting) thus decreasing the ocean temperature. This term is only active when SSTocean>Tfreeze, as you explained earlier.
>
> Other terms do not seem to contribute significantly (especially not d_HEFF/SNWbyNEG).
>
> Do you have an explanation? Did I go wrong somewhere?
>
> Martin
>
> On Feb 3, 2012, at 2:01 AM, Ian Fenty wrote:
>
>> Martin,
>>
>> On Feb 2, 2012, at 5:02 AM, Martin Losch wrote:
>>
>>> Hi Ian,
>>>
>>> thanks for your explanation. So even without me reading it, writing it had some merit, eh? I, however, did read it, now I am curious where you think you found a bug.
>>>
>> Ultimately, I could find no bug. It was a wild goose chase that I could have avoided if I were able to maintain a complete working knowledge of the elaborate machinations and accounting required in the routine. That is to say, it's not that there isn't a bug, it's just that I haven't found it yet!
>>
>>> In the meantime, I can report, that there are
>>> 1. small negative contributions of pathological HEFF, order 1e-7 and smaller (corresponds to 0.025W/m^2, or 0.002degC/month), mostly near the ice edge, but also in individual points near the area of interest. I think that this is too small to explain the large negative temperatures.
>> I think the only way to nail this problem is to analyze each term contributing to QNET in seaice growth in the trouble area at the time when temperatures fall below T_freeze.
>>
>> Since QNET includes the contribution of shortwave radiation beneath the surface, we have to subtract out that bit: SWFRACB*(a_QSWbyATM_open + a_QSWbyATM_cover) to find the actual net heat flux operating on the surface cell. If that's negative, which if this routine is responsible for the problem it must be, then it is a simple matter of determining the responsible term and digging from there.
>>
>> I added the requisite diagnostics to seaice_growth and seaice_diagnostics_init and I will email it to you to separately along with the name of the three new diagnostic terms. Can't wait to find this bug.
>>
>> -Ian
>>
>>
>>
>>> Martin
>>>
>>> On Feb 1, 2012, at 9:32 PM, Ian Fenty wrote:
>>>
>>>> Martin et al,
>>>>> I get these very low SSTs (-20deg and less). This starts in the beginning of the "melting season" (month 6-7 of integration. I start the integration on Jan01 with no ice), when SIarea reduces from 1 to values<1 then recovers in the freezing season (months 9-10), In the next year the process starts earlier (ice cover reduces to values below<1).
>>>> This is a very important clue. The fact that your SSTs approach and then remain>= T_freeze throughout winter tells me two things. First, that the code is working as designed for flux divergence over open water and through ice cover since in winter air-sea heat fluxes divergence through both ice cover and leads are almost always positive and sea ice concentrations are almost never exactly 100%. Second, that the problem only presents in a narrow set of conditions related to cold atmospheres (even in July in the Arctic) and high QSW.
>>>>
>>>> Indeed, I have high suspicions for an issue related to QSW since it's only in June and July where the fluxes are great enough where temperatures could be driven so far down if somehow QSW fluxes were improperly handled
>>>>> Wenn I do not use MCPHEE_OCEAN_ICE_HEAT_FLUX, the surface temperatures are all near or above freezing (-1.9, actually this happens only with the non-legacy code. The legacy code has always temperature well below freezing, order -5deg in mycase). So my first idea is similar to Gael's: the ocean surface looses heat through openings (AREA<1), and then for some reason cannot use this negative heat to grow ice again.
>>>> I have to disagree. I don't think that heat loses through openings is the problem per se. Heat is lost through openings everywhere there is ice and your problem occurs only in a very narrow set of summertime conditions.
>>>>> You guesses:
>>>>> 1. pathological treatment:
>>>>> case 1: negative ice: if heff<0, this contributes to d_HEFFbyNEG>0 and d_HEFFbyNeg is substracted from QNET, making the heat loss smaller, not larger, so it tends to increase, not decrease theta. It also makes sense physically: removing ice (melting) extracts temperature (enthalpy), creating ice (freezing) should increase water temperature. Did I get this right?
>>>> You're right. Not enough coffee and I forgot the subtle arguments around negative ice. :)
>>>>
>>>> Negative ice --> positive d_HEFFbyNEG
>>>> HEFF = HEFF (negative) + d_HEFFbyNEG (positive) = 0
>>>> QNET = other terms - d_HEFFbyNEG
>>>> so QNET is decreased by d_HEFFbyNEG which is indeed is warming.
>>>>
>>>>> 2. Why should shortwave radiation change it's sign and nobody noticed? (needs to be checked)
>>>> It's not that QSW's sign changed and is now wrong everywhere, it's just that perhaps it is not being treated correctly where ice is present. I'm checking since I have a feeling a bug is hiding there.
>>>>
>>>>> 3. I think that advection is very unlikely the culprit for such large deviations (but I have not proven that, just my "feeling").
>>>> Oh good.
>>>>
>>>>> Now what's different between the MCPHEE and not MCPHEE (summarizing the code):
>>>> You've nailed the summary I'd say.
>>>>> If I understand your description, then you are saying that this mechanism (that a_QbyOCN contributes to freezing) is not physical (or not correct).
>>>> I'm not saying the mechanism is not correct exactly, I'm saying that unless we're in supercooling land, the temperature of an infinitesimal volume of liquid should approach the freezing point and then freeze when subjected to heat loss (excluding pressure and salt effects). While on the other hand, the bulk temperature of a much much larger liquid volume (like our model grid cells) can have ice and liquid together with a bulk liquid temperature above freezing. Removing enthalpy from that volume can generate ice without necessarily reducing the liquid bulk temperature. If we're careful, we will never need to "correct" the situation where we find T< T_f. It's just an alternative way of handling freezing that happens to work very well for the adjoint and is physical.
>>>>
>>>>> Physically, when water is cooled below freezing temperatures, it should freeze. Can you explain, where in the non-legacy code, this is happening? Or how/where does the non-legacy code make sure that the computed heat fluxes (by the seaice-package/seaice_growth) do not lead to temperatures below freezing?
>>>> I will explain it but you have to promise to read what I'm going to write!
>>>>
>>>> * Step 1 (Called Part 2 in the code):
>>>> First we calculate potential heff increments from turbulent air-sea fluxes across open water and the ice cover, shortwave fluxes across open water and ice cover, and turbulent fluxes between the ocean and any existing ice. These are respectively termed a_QbyATM_open, a_QbyATM_cover, a_QSWbyATM_open, a_QSWbyATM_cover, and a_QbyOCN. We then use these variables but with an "r_" instead of "a_" to mean "residual" for accounting purposes. Negative values mean potential melting.
>>>>
>>>> * Step 2 (Called Part 3 in the code):
>>>> Next we apply these potential heff increments to the variable HEFF, modifying them as necessary.
>>>>
>>>> ** sub step 1: delta HEFF from turbulent ocean-ice fluxes
>>>>
>>>> Consider the "normal" wintertime case where there is nonzero ice present, T>= T_freeze, but the ocean T is close to T_f such that ocean heat cannot melt all existing ice in one time step. Then the lines:
>>>> d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))
>>>> r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
>>>> HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)
>>>>
>>>> yields: d_HEFFbyOCNonICE = r_QbyOCN, and r_QbyOCN becomes 0, and HEFF is reduced leaving no more ocean heat available this time step to melt ice.
>>>>
>>>> ** sub step 2: delta HEFF from turbulent air-sea fluxes over ice cover.
>>>> consider the case of net divergence (r_QbyATM_cover> 0)
>>>>
>>>> tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
>>>> c Limit ice growth by potential melt by ocean
>>>> & AREApreTH(I,J) * r_QbyOCN(I,J))
>>>>
>>>> d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
>>>> r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
>>>> HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
>>>>
>>>> yelds: d_HEFFbyATMonOCN_cover = r_QbyATM_cover (since r_QbyOCN = 0), r_QbyATM_cover becomes 0, and HEFF is increased.
>>>>
>>>> ** sub step 3: delta HEFF from turbulent air-sea fluxes over open water,.
>>>> consider net divergence (r_QbyATM_open> 0) and no shortwave radiation.
>>>>
>>>> The argument is the same for sub step 2.
>>>> yelds: d_HEFFbyATMonOCN_open = r_QbyATM_open (since r_QbyOCN = 0), r_QbyATM_open becomes 0, HEFF is increased leaving no residual.
>>>>
>>>> * Step 3 (Called Part 7 in the code): determineoceanmodelforcing
>>>> Using these three terms described above, QNET is redefined as:
>>>>
>>>> QNET = r_QbyATM_cover(I,J)+r_QbyATM_open(I,J)- d_HEFFbyOCNonICE(I,J)
>>>>
>>>> which yields:
>>>> QNET = 0 + 0 - (a small, negative number) = a small positive number.
>>>>
>>>> Thus, only the heat which is taken from the ocean to melt ice reduces it's temperature. Heat taken out of the ocean from the atmosphere converts liquid water at the freezing point to ice leaving the model grid cell's bulk temperature unchanged but its bulk enthalpy reduced (if you include the sea ice in the grid cell). Since the amount of heat that is used to melt ice is a fraction of the total available heat, ocean temperatures should never fall below T_f.
>>>>
>>>> Anyhow, during the course of writing this email I think I found the bug. More later.
>>>> -Ian
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>> Martin
>>>>>
>>>>> On Feb 1, 2012, at 2:13 AM, Ian Fenty wrote:
>>>>>
>>>>>> Martin and whomever else is following,
>>>>>>
>>>>>> My formulation of seaice growth (i.e. seaice_growth with the flags I suggested) shouldn't generate SSTs below the freezing point from any non-pathological thermodynamic process. By design, all heat fluxes out of the ocean surface to the atmosphere are used to grow ice, not reduce seawater temperature.
>>>>>>
>>>>>> The normal way that seawater temperatures fall when ice is present with my flags occurs with the McPhee heat flux patermeterization [MHFP] in which a fraction (read: less than 100%) of the seawater enthalpy available to melt -- that is, (T_ocean - T_freezing)*cp*rho_ice -- is removed from the seawater to melt ice. When everything is working correctly, T_ocean quickly approaches, but never falls below, T_freezing when ice is present.
>>>>>>
>>>>>> Using the MHFP, the existence of SSTs below the freezing point (which again are never supposed to happen) are a warning sign that something else is going wrong that needs to be addressed, not an indication of a failure of seaice_growth.
>>>>>>
>>>>>> Therefore, I can make the following guesses about where your strange SSTs come from:
>>>>>> 1) treatment of pathological HEFF/AREA/SNOW values when we enter the subroutine
>>>>>> 2) error in shortwave radiation fields which somehow changes their sign
>>>>>> 3) errors from non-minima preserving theta advection
>>>>>>
>>>>>> Point 1: pathological treatment
>>>>>> ---------------------------------------------------
>>>>>> Problem: When seaice_growth begins, we handle some pathological cases in a way which can generate SSTs below T_freeze.
>>>>>>
>>>>>> case 1: negative heffs or negative hsnow
>>>>>> how we treat: extract seawater enthalpy to spontaneously create an equal and opposite quantity of ice such that the sum of old + new ice is zero regardless of whether there is any available enthalpy to do it.
>>>>>>
>>>>>> case 2: "very thin ice"
>>>>>> how we treat: extract seawater enthalpy to melt all of the "very thin ice" and also any snow that happens to be around, regardless of whether there is any available enthalpy to do it.
>>>>>>
>>>>>> Thus, if somehow negative or "very thin" ice keeps appearing in a grid cell, seawater temperatures can be driven far below T_freeze.
>>>>>>
>>>>>> Diagnosis: Make a diagnostic for d_HEFFbyNEG and d_HSNWbyNEG (they don't exist now) and see if they are nonzero near where you have the problem.
>>>>>>
>>>>>> Solution 1: fix sea ice advection so that negative values don't come out.
>>>>>> Solution 2: change pathological treatment such that you don't remove energy that isn't there!
>>>>>> Solution 3: stop caring about "very thin" ice (change siEps to zero) if that is the culprit
>>>>>>
>>>>>>
>>>>>> Point 2: QSW
>>>>>> ---------------------------------------------------
>>>>>> If somehow, a_QSWbyATM_cover or a_QSWbyATM_open are positive (i.e. indicate a shortwave radiative flux out of the ocean), heat will be removed from those grid cells "penetrated" by the erroneous QSW. This cooling will drive seawater temps far below the freezing point.
>>>>>>
>>>>>> Diagnosis: Since QSW should only ever be positive, it would be easy to diagnose this problem by checking a_QSWbyATM_cover and a_QSWbyATM_open and making sure they always have the right sign. Again, there don't seem to be built-in diagnostics for those terms so you'll have to add them.
>>>>>>
>>>>>> Solution: fix QSW sign.
>>>>>>
>>>>>> Point 3: Errors from non-minima preserving theta advection.
>>>>>> ---------------------------------------------------
>>>>>> I don't convert negative SSTs to ice so if the advection scheme doesn't preserve minima, SSTs below T_freeze could result. If there is a place where advection errors keep driving theta down, you could end up with places with very very negative SSTs.
>>>>>>
>>>>>> Diagnosis: Check theta pre and post advection to make sure T>= T_freeze going in and out.
>>>>>>
>>>>>> Solution: Put in special pathological case handling advection errors only.
>>>>>>
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> MITgcm-devel mailing list
>>>>>> MITgcm-devel at mitgcm.org
>>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>> _______________________________________________
>>>>> MITgcm-devel mailing list
>>>>> MITgcm-devel at mitgcm.org
>>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>>
>>>> _______________________________________________
>>>> MITgcm-devel mailing list
>>>> MITgcm-devel at mitgcm.org
>>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>>
>>> _______________________________________________
>>> MITgcm-devel mailing list
>>> MITgcm-devel at mitgcm.org
>>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>
>> _______________________________________________
>> MITgcm-devel mailing list
>> MITgcm-devel at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
>>
>>
>> _______________________________________________
>> MITgcm-devel mailing list
>> MITgcm-devel at mitgcm.org
>> http://mitgcm.org/mailman/listinfo/mitgcm-devel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mitgcm.org/pipermail/mitgcm-devel/attachments/20120203/d76d56f8/attachment-0001.htm>
More information about the MITgcm-devel
mailing list