Crash in extremely stable conditions - TKE production at the surface

[fjansson]

We have had a crash in the Dales surface scheme under stable conditions, over land. The symptoms are a sudden huge e12 - tendency at the surface, which causes the whole simulation to break and stop with a thermo crash a few steps later.
We run with isurf=4 - imposed surface fluxes. The crash occur in the evening when the heat flux is negative and the wind is very weak. The Obukhov length is then positive and small. A more stable surface layer -> smaller Obukhov length.

modsurface.f90, line 984:

phimzf = (1.+5.*zf(1)/obl(i,j))
...
dudz  (i,j) = ustar(i,j) * phimzf / (fkar*zf(1))*(upcu/horv)

if the Obukhov length is small here, dudz can be very large.
modsubgrid, line 437++: dudz is used to calculate the tke tendency at the surface, via the shear term sbshr. A large dudz results in the large tke tendency which leads to the crash. Note that this is in a situation with very weak wind, so there should be very little shear.

Another related issue is that ustar is forced to be >= 0.01. Small Obukhov length corresponds to a small ustar, so without the limitation a small ustar would to some extent compensate for the large phi. note however that ustar also is present in a denominator, so that tiny values could potentially give problems here.

dthldz(i,j) = - thlflux(i,j) / ustar(i,j) * phihzf / (fkar*zf(1))

Fix

phimzf is the phi function in Monin-Obukhov theory. It is used here to determine du/dz at the first layer at height zf(1). There is also a psi function, which determines the actual velocity as function of z. Psi is defined as an integral containing phi, and is used in determining the Obukhov length.

• a) The phi function used is not valid when z/obl is larger than ~1. See e.g. Optis 2015
• b) The phi and psi functions used in Dales are not consistent with each other - the psi function reported to have the correct z/obl -> infinity limit.

a fix for a) is to add a hard cap on both phi functions at z/obl >= 1. This is one of the simplest phi functions listed by Optis, and it solves the problem in our case. Proposed fix here: CloudResolvingClimateModeling@8469307

[Chiil]

Hi Frederik, I think that running the MO solver as it is in DALES with fixed fluxes has another complication, that might be worth to check if you doing SBLs with fixed fluxes. The bulk Richardson number that is used is locked at a certain temperature difference, whereas the solver should update this difference along the iterations. I have a fix for this, which I can show you in our meeting next month. I think the hard cap will keep the model alive, but you are entering the realm of stabilities where doing LES doesn’t make too much sense anymore. Cheers, Chiel

…

On 29 Jan 2018, at 16:53, Fredrik Jansson ***@***.***> wrote: We have had a crash in the Dales surface scheme under stable conditions, over land. The symptoms are a sudden huge e12 - tendency at the surface, which causes the whole simulation to break and stop with a thermo crash a few steps later. We run with isurf=4 - imposed surface fluxes. The crash occur in the evening when the heat flux is negative and the wind is very weak. The Obukhov length is then positive and small. A more stable surface layer -> smaller Obukhov length. modsurface.f90, line 984: phimzf = (1.+5.*zf(1)/obl(i,j)) ... dudz (i,j) = ustar(i,j) * phimzf / (fkar*zf(1))*(upcu/horv) if the Obukhov length is small here, dudz can be very large. modsubgrid, line 437++: dudz is used to calculate the tke tendency at the surface, via the shear term sbshr. A large dudz results in the large tke tendency which leads to the crash. Note that this is in a situation with very weak wind, so there should be very little shear. Another related issue is that ustar is forced to be >= 0.01. Small Obukhov length corresponds to a small ustar, so without the limitation a small ustar would to some extent compensate for the large phi. note however that ustar also is present in a denominator, so that tiny values could potentially give problems here. dthldz(i,j) = - thlflux(i,j) / ustar(i,j) * phihzf / (fkar*zf(1)) Fix phimzf is the phi function in Monin-Obukhov theory. It is used here to determine du/dz at the first layer at height zf(1). There is also a psi function, which determines the actual velocity as function of z. Psi is defined as an integral containing phi, and is used in determining the Obukhov length. a) The phi function used is not valid when z/obl is larger than ~1. See e.g. Optis 2015 <http://onlinelibrary.wiley.com/doi/10.1002/we.1883/full> b) The phi and psi functions used in Dales are not consistent with each other - the psi function reported to have the correct z/obl -> infinity limit. a fix for a) is to add a hard cap on both phi functions at z/obl >= 1. This is one of the simplest phi functions listed by Optis, and it solves the problem in our case. Proposed fix here: ***@***.*** <CloudResolvingClimateModeling@8469307> — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#30>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AAGu6bEQiZyYN4uY11783DsMCKHbUPomks5tPemQgaJpZM4Rw2Pv>.

[fjansson]

Hi Chiel,

Good point about the fluxes and Rib, I'm interested in the fix. I will not make it to the February meeting though - maybe you can discuss it with Pier, Stephan or Gijs.

Do you mean that LES is not very interesting in these situations, or that the results are not reliable? I get these problems when running Dales connected to a large-scale model, then Dales has to deal with whatever conditions the large-scale model gives it. Someone else, connecting Dales to another large-scale model, has also been suffering from crashes in stable situations, so this would be nice to get fixed.

Best,
Fredrik

[fjansson]

Fixed in version 4.2:
Commit 4098295 changes the phi function in the Monin-Obukhov similarity theory to behave better in stable conditions.
Commit db0c242 turns off the sgs_surface_fix by default. With the new phi definitions, the fix is no longer needed.