Two questions about `atmo_boundary_layer`

[phil-blain]

Hi everyone, I'm working on upstreaming some changes we've made to subroutine atmo_boundary_layer in module icepack_atmo.

I have two questions about some part of the code that I'm not sure about:

1. In the iteration part of the subroutine, we compute the stability function in the stable case following Jordan et al 1999: https://github.com/CICE-Consortium/Icepack/blob/5cf223287d06167bb813d1fc2248258034511017/columnphysics/icepack_atmo.F90#L287-L294

However, at the end of the subroutine when we compute the diagnostic temperature Tref and humidity Qref, we instead use the formulation \psi_m = -5 \zeta (see for example Kauffman & Large p.40). See the first line of the following code fragment: https://github.com/CICE-Consortium/Icepack/blob/5cf223287d06167bb813d1fc2248258034511017/columnphysics/icepack_atmo.F90#L367-L374

Is this an implementation error ? I'm guessing it should not make a big impact but for consistency shouldn't we use the same formulation in both places ?

---

2. At the end of the subroutine, when we compute the diagnostic temperature Tref and humidity Qref, we compute the stability parameter zeta (hol in the code) at height zTrf by multiplying the existing hol value by zTrf / zlvl: https://github.com/CICE-Consortium/Icepack/blob/5cf223287d06167bb813d1fc2248258034511017/columnphysics/icepack_atmo.F90#L364 instead of reusing the same formula as in the iteration (with zTrf instead of zlvl): https://github.com/CICE-Consortium/Icepack/blob/5cf223287d06167bb813d1fc2248258034511017/columnphysics/icepack_atmo.F90#L278-L281

Both approaches do not give the same value because the former implicitly uses {u,q,t}star from the penultimate iteration instead of the last iteration (which would be the case if we would use the longer formula)...

Is the underlying assumption that if we are converged it should not matter? Because I've made a quick test and there can be some differences (in some cases it doubles the value of hol, others it's the same until the third decimal, sometimes it's at 1E-6, etc....) I'm not very familiar with the physics at play here so I don't know if in the end it would significantly change the model answers. I'm just sanity-checking what I'm understanding from my reading of the code...

@dabail10 @eclare108213 @proteanplanet (@JFLemieux73 )

[dabail10]

I guess I've never looked at this code very closely before. The while loop was added recently, but the default behaviour should still be the same where it always does 5 iterations. I'll have to do some digging. We have very similar code in the CESM coupler (shr_flux_mod.F90).

[eclare108213]

The original code came from CESM (a long time ago, before it was CESM), and I added the Jordan stability part later. I'll admit I've never paid much attention to the diagnostic fields Tref and Qref -- those were from CESM, for CESM, and I haven't gotten feedback on whether they're doing the right thing until now! They shouldn't matter for the solution because they're only diagnostic, right? The idea there was to move the output data to a standard reference level, maybe 2 m (I forget). There have been some modeling centers that decided to do more iterations than 5 in order to converge better, and Andrew added the wind dependence (along with the while loop) for RASM a few years ago. Some groups don't use this calculation at all (they input wind stress). It's good to look at it more closely - thanks.

[proteanplanet]

The reason this code is the way it is that the Community Atmosphere Model and other atmospheric models, too, require consistent state fields at the standard reference values for their own diagnostics and physics, and this must be calculated uniformly across all coupled components - the ocean, land and sea ice, and be consistent with the atmospheric model. Conversely, each of these components has their own physical requirements for passing the turbulent fluxes derived for their particular surfaces - the ocean with waves, land with vegetation and topography, and sea ice with surface heteogeneity, which are used by the atmosphere's physics package as flux fields. As a consequence, what looks like an inconsistency, is in fact physically correct according to the particular atmospheric model used - CAM in the case of CESM, and EAM in the case of E3SM, as two examples. Should the algorithm used be inconsistent for your particular application, then we could place flags around this code to signal preference of the scheme used, and you could add calculations consistent with the atmospheric coupling of your choosing.

[dabail10]

Thanks for pointing this out @proteanplanet. These are intended as diagnostics for the atmospheric component. They are not used elsewhere in the CICE code as far as I am aware. Here is the equivalent code in CESM shr_flux_mod.F90:

      !------------------------------------------------------------
      ! compute diagnositcs: 2m ref T & Q, 10m wind speed squared
      !------------------------------------------------------------
      hol = hol*ztref/zbot(n)
      xsq = max( 1.0_R8, sqrt(abs(1.0_R8-16.0_R8*hol)) )
      xqq = sqrt(xsq)
      psix2   = -5.0_R8*hol*stable + (1.0_R8-stable)*psixhu(xqq)
      fac     = (rh/loc_karman) * (alz + al2 - psixh + psix2 )
      tref(n) = thbot(n) - delt*fac
      tref(n) = tref(n) - 0.01_R8*ztref   ! pot temp to temp correction
      fac     = (re/loc_karman) * (alz + al2 - psixh + psix2 )
      qref(n) =  qbot(n) - delq*fac

[phil-blain]

Another thing I'm noticing just now is that the unit step function stable is not re-computed for the diagnostics, so it is still evaluated at the original (non-shifted) hol value...

[TillRasmussen]

I have a related question that is somewhat related. We have seen errors where icepack aborts due to vertical therm error (conservation of energy). This is with atmospheric forcing (ERA5). This is often triggered by unrealistic high latent heat (+1000 W/m^2) at low windspeeds (typically less than 1)and very small ice velocities. The dffference in humidity is around ~10^-3 but lhcoef becomes large. This is with formdrag and high frequent forcing turned on. One effect of adding the high frequent forcing is that the lower treshhold of the magnitude of wind (reduced by the ice speed) is set to 0.5 instead of 1. What is the reasoning behind this. Many of the high latent heat occurencies would be avoided by increasing this to one.

[mhrib]

Possible singularity / division by zero?? https://github.com/CICE-Consortium/Icepack/blob/5cf223287d06167bb813d1fc2248258034511017/columnphysics/icepack_atmo.F90#L298-L299

For psixh=alz+vonkar/rhn , the denominator will be zero and rh infinity. Thereby shcoef will get large(infinity) and so will fsensn: https://github.com/CICE-Consortium/Icepack/blob/5cf223287d06167bb813d1fc2248258034511017/columnphysics/icepack_therm_shared.F90#L150-L151 And similar for ren, re, lhcoef and flatn.

Thereby it seems, that it is possible to end up with unphysical huge sensible and latent heat fluxes for very special cases.

As I can see, this can happen if hol>0: Then psixh reduces to psimhs , as stable=1. for hol>=0

However, I do not have a (ready to use physical) solution to avoid it...

[proteanplanet]

@mhrib I'm struggling to understand why psixh=alz+vonkar/rhn would occur. Have you encountered this in your simulations?

[proteanplanet]

I have a related question that is somewhat related. We have seen errors where icepack aborts due to vertical therm error (conservation of energy). This is with atmospheric forcing (ERA5). This is often triggered by unrealistic high latent heat (+1000 W/m^2) at low windspeeds (typically less than 1)and very small ice velocities. The dffference in humidity is around ~10^-3 but lhcoef becomes large. This is with formdrag and high frequent forcing turned on. One effect of adding the high frequent forcing is that the lower treshhold of the magnitude of wind (reduced by the ice speed) is set to 0.5 instead of 1. What is the reasoning behind this. Many of the high latent heat occurencies would be avoided by increasing this to one.

@TillRasmussen High frequency ensures that stress can go to zero while the turbulent heat flux terms (sensible, latent) do not go to zero with no wind forcing. High frequency also ensures that the sea ice can complete inertial loops with the apparent wind speed, bearing in mind that during an inertial oscillation with a weak geostrophic flow, the ice will at times be traveling in exactly the opposite direction to the wind. This also depends on oceanic coupling, since the total ice-ocean mass transport is integral to a correct Ekman solution.

I would be testing this with and without form drag, as previous problems have been encountered with that parameterization.

[mhrib]

@proteanplanet No, we (Till and I) have not encountered this in our simulations. It is just considerations, but seems to be relevant to mention here.

We recently got unrealistic high fluxes of latent heat (flatn > +/- 1000 W/m2) once in a while using ERA5 forcing, and formdrag+highfreq turned on. We tried to debug writing out a lot of variables when it happens and found, that lhcoef takes very high values (order of +/- 10**7), and abs(re)>10. Most of the extreme values seems to be during low wind speed (minus ice velocity) as was also noted by @TillRasmussen above.

While debugging, we were looking into the equations and speculate, if the equations could lead to extreme high values in rare situations.

[proteanplanet]

@mhrib If time allows, we could discuss this at the end of the CICE call tomorrow. Previous application of form drag to RASM (fully coupled, regional model) proved difficult and was abandoned.

[TillRasmussen]

After the meeting yesterday I calculated some diagnostics of possible solutions for the variable re (line 298). x axis is ren which is the square root of the atmospheric drag coefficient (Cdn_atm). This is limited to .02 in the form drag calculations. hol is limited to [-10 10] (yaxis). All values where hol is negativ results in an instable atmsophere, which only happens when 2 meter atmospheric temperature is higher than the surface. Note that the color scheme indicates log(re). The white area is negative re, which result in latent heat flux into the ice, which seems unrealistic. Otherwise note how fast re increases in the unstable regime. In general we have issues when the drag coefficent is high and the winds are low (low ustar).

[eclare108213]

Returning to @phil-blain's questions about subroutine atmo_boundary_layer in module icepack_atmo.

1. In the iteration part of the subroutine, we compute the stability function in the stable case following Jordan et al 1999 However, at the end of the subroutine when we compute the diagnostic temperature Tref and humidity Qref, we instead use the formulation \psi_m = -5 \zeta (see for example Kauffman & Large p.40). Is this an implementation error ? I'm guessing it should not make a big impact but for consistency shouldn't we use the same formulation in both places ?

This sounds like an implementation error. If it only affects Tref and Qref then it doesn't matter for the (ice) solution, but it would better to be consistent. Changing it should be BFB in stand-alone configurations for prognostic variables.

1. At the end of the subroutine, when we compute the diagnostic temperature Tref and humidity Qref, we compute the stability parameter zeta (hol in the code) at height zTrf by multiplying the existing hol value by zTrf / zlvl instead of reusing the same formula as in the iteration (with zTrf instead of zlvl)

Both approaches do not give the same value because the former implicitly uses {u,q,t}star from the penultimate iteration instead of the last iteration (which would be the case if we would use the longer formula)...

Is the underlying assumption that if we are converged it should not matter? Because I've made a quick test and there can be some differences (in some cases it doubles the value of hol, others it's the same until the third decimal, sometimes it's at 1E-6, etc....) I'm not very familiar with the physics at play here so I don't know if in the end it would significantly change the model answers. I'm just sanity-checking what I'm understanding from my reading of the code...

Likewise, it would be better to be consistent. It's kind of you to offer the assumption of an underlying assumption -- in reality it's more likely that we (sea ice modelers) just haven't paid attention to Tref and Qref. I'd like @dabail10's input on this from the CESM point of view.

[eclare108213]

@TillRasmussen I'm not sure where to go next with your comments. What is your conclusion, or do you have a request for us to follow up?

[TillRasmussen]

I think that my conclusion for now is the same same @proteanplanet that the combination of formdrag and this calculation of latent and sensible heat provides unrealistic high fluxes in special cases, thus it should be used with cautious until it is improved or a potential bug is found. I will try to dig a bit further into this. We have seen this when running stand alone, coupled to nemo 4 and to HYCOM. All forced with ERA-5. A slight improvement was found when increasing the number of iterations to 50.

[dabail10]

Can we close this?

[proteanplanet]

Before we close this, I'll ask @erinethomas if she can post her findings from E3SM on stability, so it's stored away for future.

[erinethomas]

Here is a summary of the stability in the flux calculations within E3SM/MPAS-SI Column physics:

FluxConvergence_Summary_short.pptx

In short, E3SM flux calculations are stable in the sea ice column physics and converge very quickly (in fewer than 5 iterations).

[dabail10]

Any updates on this one? Can we close it?

[proteanplanet]

I think this can be closed.