No clouds are forming in the bomex validation case

[mmr0]

We are attempting to use BOMEX observations/LES as a test case for Breeze (latest in this branch). The experimental design is from Siebesma et al. (2003). Briefly, these are 6-hour runs of marine shallow cumulus convection, with a well-mixed 500m-thick "dry" convective layer topped by a 1000m thick conditionally-unstable cloud layer. These clouds do not form in our Breeze run, and I'm not sure why.

Model setup

The model is forced with geostrophic winds, prescribed heat and vapor fluxes at the air-sea interface, and prescribed friction velocity. Other forcings include pressure gradient forcing for geostrophic balance, prescribed radiative cooling and drying profiles, as well as a subsidence forcing that introduces tendencies in u, v, pot. temp., and humidity.

Initial conditions:

The subsidence profile, and "radiative cooling" & "advective drying" forcings:

The tendency terms due to subsidence (as calculated at the first timestep):

These all look reasonable. The v tendency is not shown as it is zero at the first timestep since v=0. The subsidence leads to warming/drying as expected. I haven't shown the geostrophic forcing here, but since the wind is steady in time (outside of the actively-mixing lower layer) I'm confident that this is working as expected.

Some actual output

This is plane- and time-averaged output (time-averaging over hours 3-6) of: potential temperature; specific humidity (vapor); velocity; specific humidity (liquid).

Compared to the Siebesma et al paper:

• Our mixed layer (ML) is deeper
• We have no liquid water (i.e. no clouds)
• Our ML is too moist/cloud layer is too dry
• Velocity profiles are somewhat different, likely due to deeper ML/less mixing in cloud layer

Here is a movie:

bomex_xz.mp4

[Edit 21/07] As a quick check that the heat/vapor fluxes were doing the right thing, I removed the geostrophic forcing & velocity IC to run convection-only cases. Some relevant profiles here:

Why no clouds? Some ideas....

1. The problem setup isn't quite right

I have growing confidence in the problem setup based on the profiles above + the tendency terms which appear to be doing what they're supposed to do, according to descriptions in Siebesma et al 2003.

2. The thermodynamics isn't quite right

We know we can get liquid water (e.g. the free convection example), so that's something. I feel very ill-equipped to check the code, but am happy to try, or to run any tests we can think of!

3. The LES approach isn't appropriate

We have tried both implicit LES (using a WENO advection scheme) and an Anistotropic Minimum Dissipation (AMD) closure. AMD resulted in lots more mixing above the ML but did not help with clouds. It could be that neither of these approaches is appropriate for this coarse boundary layer setup.
[Edit 18/7] - Found a nice reference for implicit LES of clouds which indicates WENO for momentum/scalar advection is the best choice. I note that there is an eddy diffusivity model applied to the surface layer though which we could implement easily https://pmc.ncbi.nlm.nih.gov/articles/PMC5586241/

Some things to try/check

• Compare our turbulent flux profiles to Siebesma et al.
• Also compare them to the flux BC's we impose at the lower boundary
• Sensitivity to resolution. In theory, if the closure (or lack thereof) is having an unintended effect we should at least see some change to the soln as we change resolution. I have done this previously, but could repeat more carefully now that we are confident in the problem setup.
• ???

[glwagner]

@trontrytel

[mssingh]

Hi Maddie,

As you say, you seem to get strong dry convection and turbulence in the mixed layer, but no condensation, and therefore no ability to penetrate the conditionally unstable layer above.

Can't say exactly what is going on here, but here are a couple of questions/things to try.

What is the saturation specific humidity formulation in Breeze? Can you calculate the relative humidity? If clouds formed, one would expect you would get RH=1 at the top of the mixed layer in regions of upward motion. How close is the model to this?

I agree about plotting the flux profiles. You should be imposing a latent heat flux of something like 150 W/m^2 - where is this water going?

If everything looks right and you just don't seem to be getting enough moisture to saturate, perhaps reduce or remove the drying. This should moisten the domain pretty quickly and you should get some sort of clouds - at least to check that clouds do form in this kind of setup.

[mmr0]

Thanks for the suggestions Marty - will look into them. I also found the configuration & output from DALES for this case, which might be a useful point of reference: https://github.com/dalesteam/dales/tree/main/cases/bomex

[mssingh]

15 years ago! I am very pleased someone decided to upload this in 2010.

[mmr0]

I know I did laugh when I looked at the dates

[trontrytel]

Nice to meet you all! Apologies I can't help in the next 7 days - I'll be traveling for a conference. But I would be very happy to help work on this when I'm back. (And then add all our microphysics schemes as options).

I did plot RH at some point. It looked like the plumes were getting really close to saturation at a correct cloud base height, and then getting "smeared" by the flow instead of producing clouds. But at that point I was running the simulation more as a 2D case for a faster turn-around time. And was wondering if me taking shortcuts in the setup was causing the behavior.

Btw, the LES case is based on measurements done in 1969. BOMEX is the T-REX of cloud LES cases.

[mmr0]

You too @trontrytel! That would be amazing, thank you :) I will update this issue if I learn anything new. I wonder if that fact the the mixed layer is too deep & coincides with the shear in the velocity profile @ 700m is causing the "smearing" that you describe...

[glwagner]

Ok, here are a few preliminary results. Adding RH to this case, I find that it reaches 90-92% when I use the resolution reported in Siebesma et al 2003, plus WENO(order=9) and no explicit subgrid closure. Interestingly, I did a quick test with a reduced resolution of Nx = Ny = 32 (rather than Nx = Ny = 64) and found somewhat higher values of RH, up to 94% (after 1.5 hours).

Now, reading Siebesma et al 2003, I find this table which details the numerical methods + closures used in various models:

I'm noting that the codes with an "all or nothing" (AN) assumption for subgrid saturation (eg, computing saturation solely in terms of the resolved humidity/temperature with no specific subgrid saturation model) mostly use centered second-order advection with SmagorinskyLilly closure. So I added support for SmagorinskyLilly to MoistAirBuoyancy in #22. This change produced slightly larger RH up to 96%, but no saturation within 40 minutes of simulation (we can return to this with longer simulations if we like). Note, SmagorinskyLilly runs significantly slower (right now) because have not implemented a way to store the computation of temperature / saturation adjustment result in the Boussinesq code (we do support this in AtmosphereModel, which we aren't testing yet). As a result we have to recompute temperature to get N^2 for the Smag Lilly scheme.

Finally, I noticed that we were not prescribing initial perturbations as prescribed:

I fixed this, but it seems to only have a minor effect.

[glwagner]

Here is something more interesting perhaps. Siebesma et al state that the sea surface potential temperature is 299.1 K, which should correspond to sea surface temperature 300.4 K:

I don't understand this. What reference pressure are they using then? I am honestly a little confused by this, because then they have a non-unity Exner function at the sea surface? Our thermodynamics is right now formulated so that the pressure profile is given by:

Breeze.jl/src/Thermodynamics/reference_states.jl

Line 55 in 931dedd

return p₀ * (1 - g * z / (cᵖᵈ * θᵣ))^(Rᵈ / cᵖᵈ)

where p₀ is the "base pressure" (sometimes this is written p00) and is equivalent to the pressure at z=0. So we simply cannot describe a non-unity Exner function at z=0 with this formulation. Not sure it matters but it means that our temperature and saturation specific humidity are not 300.4 K and 0.02245 at the surface:

julia> T[1, 1, 1]
298.7584708527437

julia> qᵛ★[1, 1, 1]
0.02013170843706923

[mmr0]

Not sure it matters but it means that our temperature and saturation specific humidity are not 300.4 K and 0.02245 at the surface

I don't think it does - my read on that section of the paper is that that's the SST to prescribe if you're calculating the fluxes rather than prescribing them as in this experiment.

Nothing to add to your technical question though. Thanks for looking into this!

[trontrytel]

The potential temperature by definition uses p=1000 hPa as reference. So unless our surface pressure is also 1000 hPa, we should see non-unity Exner function at the surface

See for example AMS Glossary: https://glossary.ametsoc.org/wiki/Potential_temperature

[glwagner]

The potential temperature by definition uses p=1000 hPa as reference. So unless our surface pressure is also 1000 hPa, we should see non-unity Exner function at the surface

See for example AMS Glossary: https://glossary.ametsoc.org/wiki/Potential_temperature

Ok, we need to add an extra parameter for the Exner value.

[glwagner]

This may have been the bug we were looking for: 133f008