Instability develops over long integrations

[p-hss]

Running SpeedyWeather v0.10.0 for more than a year leads to instabilities and a blow-up.

I ran Speedy with the following setup:

     spectral_grid = SpectralGrid(trunc=31,
                                   nlev=8,
                                   Grid=FullClenshawGrid)

     condensation = SpeedyWeather.ImplicitCondensation(spectral_grid)

     output = OutputWriter(spectral_grid,
                            PrimitiveWetModel,
                            #output_dt=Day(1),
                            output_dt=60*30,
                            path=path,
                            id=id)

     shortwave_radiation = NoShortwave() # enabling it also leads to instabilities

     model = PrimitiveWetModel(;spectral_grid,
                               large_scale_condensation=condensation,
                               shortwave_radiation=shortwave_radiation,
                               output=output)

     simulation = initialize!(model)

     run!(simulation,
          period=Day(num_days),
          output=true)

The spatial mean over time looks something like this before the blow-up:

And the vertically mean at the time step before the blow-up:

[milankl]

@p-hss Thanks for pointing this out!! (A bit embarrasing to be fair...) I've now moved all parameterizations from i to i-1 which should be more stable as also outlined in #445. Additionally I've put back in the older (harder) flux limiters, that prevent tendencies from being larger than some threshold. There's no good theory what these should be, it's just a dirty fix to prevent a blow up in some weather conditions. But because we have to guess/experiment these limits, it's hard to say whether they just prevent a blow up or also affect dynamics more seriously. I'd love to switch them off completely, but they're a preliminary solution ...

Could you repeat this analysis with the #mk/allprev branch? I.e. in the package manager do

] add https://github.com/SpeedyWeather/SpeedyWeather.jl#mk/allprev

[milankl]

Still blows up for me after ~1yr, I'm on it.

[milankl]

Of temperature and humidity you plotted the topmost (k=1) layer I assume?

[simone-silvestri]

I tried with better vertical numerics

spectral_grid = SpectralGrid(trunc=31,
                             nlev=8,
                             Grid=FullClenshawGrid)

condensation = SpeedyWeather.ImplicitCondensation(spectral_grid)

shortwave_radiation = NoShortwave() 

model = PrimitiveWetModel(;spectral_grid,
                          large_scale_condensation=condensation,
                          shortwave_radiation=shortwave_radiation,
                          vertical_advection=SpeedyWeather.WENOVerticalAdvection(spectral_grid))

simulation = initialize!(model, time=DateTime(1940,1,1))

run!(simulation,
     period = Day(1000),
     output=true)

and it seems to be stable on main (at least up to 1000 days). Without it, it crashes at 270 days

[milankl]

Yeah that doesn't surprise me, because on main the convection is still calculated at i (the current time step) leading to oscillations that WENO probably smoothes out. While I agree that WENO is obv the better scheme, I don't think the centred advection is the culprit here as there's a bunch of other stuff going on. The latest commit in #555 now takes out all flux limiters which shouldn't be necessary anyway (Fortran SPEEDY has them everywhere though...). Thanks for checking!!!

[simone-silvestri]

Makes sense. Indeed, limiters can often mask a more fundamental problem or bug (WENO is akin to a limiter).

[milankl]

There's now also another problem that the surface drag coefficient via the Bulk-Richardson number is in most areas 0, effectively disabling all surface fluxes leading to a cooling of the surface temperatures for example. E.g. -14C in the tropics here 😆

[p-hss]

@p-hss Thanks for pointing this out!! (A bit embarrasing to be fair...) I've now moved all parameterizations from i to i-1 which should be more stable as also outlined in #445. Additionally I've put back in the older (harder) flux limiters, that prevent tendencies from being larger than some threshold. There's no good theory what these should be, it's just a dirty fix to prevent a blow up in some weather conditions. But because we have to guess/experiment these limits, it's hard to say whether they just prevent a blow up or also affect dynamics more seriously. I'd love to switch them off completely, but they're a preliminary solution ...

Could you repeat this analysis with the #mk/allprev branch? I.e. in the package manager do

] add https://github.com/SpeedyWeather/SpeedyWeather.jl#mk/allprev

Thanks a lot for looking into this! I just finished a 10-year stable run with the new branch. :)

Mean precipitation looks maybe a bit worse than in an older version < 0.10.0:

[milankl]

Awesome and in the latest commit, I even took out all flux limiters. Amazing. It doesn't surprise me that the precipitation is still somewhat off given that after all it's highly dependend on the vertical stability which has to be off given we switched off shortwave radiation. Honestly, I'm actually quite surprised how comparably well the zonal averages still look!

[p-hss]

That's great. I will try another longer run. :)

[milankl]

I've just put the TransparentShortwave back in, which should make the tropics less stable (not numerical stability) in the vertical and trigger therefore more convection -- but this is still very crude as it doesn't go through a skin temperature but hits the lowermost layer air temperature directly via a flux.

[p-hss]

Okay, I will try it out!

[milankl]

Awesome thanks Philipp! When you check this could you also plot surface temperature at the last time step here? Cause at the moment there's no super strong control of sst / land surface temperatures on air temperature and given the bulk Richardson-based surface drag coefficient in stable conditions the surface fluxes are zero even. So without shortwave I found surface temperatures to easily drop to 10C in the tropics. Which is somewhat fine given that for 8 layers with k=8 it's air temperatures at ~930hPa technically but also unrealistically cold (maybe not surprising when you switch off the sun).

[p-hss]

Sure! I'm getting an ERROR: LoadError: type SurfaceHeatFlux has no field max_flux error on the current #mk/allprev branch. Is it the right branch?

[milankl]

Yes right branch and you should get that error (because I removed the flux limiter) but it shouldn't be triggered. Can you post the stacktrace?

@p-hss EDIT: found it, committed a bug fix. Should be resolved! (do julia> ] update SpeedyWeather if you have that branch installed)

[p-hss]

Perfect, thanks a lot! I ran the model, and it became unstable after around 1.5 years. The surface temperature on the last day looks like this:

[milankl]

Mhhh, surface temperature doesn't look to bad, it's clearly not running away to freezing as it did before. Could you throw in your previous analysis too?

[p-hss]

Sure, here is the hourly output as spatial means before the blow-up:

and the vertical mean at the hour before the blow-up:

It looks like the instability is originating at the bottom, maybe?

[milankl]

Okay there's some weird resonance with the daily cycle in the TransparentShortwave, which was never intended to be a radiation scheme but rather a quick hack to get some solar forcing in. Use either

• shortwave_radiation = NoShortwave() to disable the shortwave radiation and accept the lower temperatures at the surface (SST forcing only)
• planet = Earth(spectral_grid, daily_cycle=false) to remove the daily cycle in the solar radiation (daily averages instead)

[p-hss]

Okay thanks, I will try it out!

[milankl]

We can also remove the radiation scheme, i.e.

• shortwave_radiation = NoShortwave()
• longwave_radiation = NoLongwave()

but use Held-Suarez forcing instead

• temperature_relaxation = HeldSuarez(spectral_grid)

?

[p-hss]

Okay there's some weird resonance with the daily cycle in the TransparentShortwave, which was never intended to be a radiation scheme but rather a quick hack to get some solar forcing in. Use either

* `shortwave_radiation = NoShortwave()` to disable the shortwave radiation and accept the lower temperatures at the surface (SST forcing only)

* `planet = Earth(spectral_grid, daily_cycle=false)` to remove the daily cycle in the solar radiation (daily averages instead)

I ran some longer simulations, and disabling shortwave radiation made the simulation stable. :)

Removing the daily cycle leads to instabilities, as does using Held-Suarez forcing instead of the radiation schemes.

[milankl]

Awesome I'll take out the shortwave radiation and merge this. Meaning a better radiation should really be high up on my todo list! Thanks Phillip for checking this.

[p-hss]

Sounds good!