Examples 3D

The following showcases several examples of SpeedyWeather.jl simulating the Primitive equations with and without humidity and with and without physical parameterizations.

See also Examples 2D for examples with the Barotropic vorticity equation and the shallow water model.

Jablonowski-Williamson baroclinic wave

using SpeedyWeather
spectral_grid = SpectralGrid(trunc=31, nlayers=8, Grid=FullGaussianGrid, dealiasing=3)

orography = ZonalRidge(spectral_grid)
initial_conditions = InitialConditions(
    vordiv = ZonalWind(),
    temp = JablonowskiTemperature(),
    pres = ZeroInitially())

model = PrimitiveDryModel(spectral_grid; orography, initial_conditions, physics=false)
simulation = initialize!(model)
run!(simulation, period=Day(9))

The Jablonowski-Williamson baroclinic wave test case^[JW06] using the Primitive equation model particularly the dry model, as we switch off all physics with physics=false. We want to use 8 vertical levels, and a lower resolution of T31 on a full Gaussian grid. The Jablonowski-Williamson initial conditions are ZonalWind for vorticity and divergence (curl and divergence of $u, v$), JablonowskiTemperature for temperature and ZeroInitially for pressure. The orography is just a ZonalRidge. There is no forcing and the initial conditions are baroclinically unstable which kicks off a wave propagating eastward. This wave becomes obvious when visualised with

using CairoMakie

vor = simulation.diagnostic_variables.grid.vor_grid[:, end]
heatmap(vor, title="Surface relative vorticity")

Held-Suarez forcing

using SpeedyWeather
spectral_grid = SpectralGrid(trunc=31, nlayers=8)

# construct model with only Held-Suarez forcing, no other physics
model = PrimitiveDryModel(
    spectral_grid,

    # Held-Suarez forcing and drag
    temperature_relaxation = HeldSuarez(spectral_grid),
    boundary_layer_drag = LinearDrag(spectral_grid),

    # switch off other physics
    convection = NoConvection(),
    shortwave_radiation = NoShortwave(),
    longwave_radiation = NoLongwave(),
    vertical_diffusion = NoVerticalDiffusion(),

    # switch off surface fluxes (makes ocean/land/land-sea mask redundant)
    surface_wind = NoSurfaceWind(),
    surface_heat_flux = NoSurfaceHeatFlux(),

    # use Earth's orography
    orography = EarthOrography(spectral_grid)
)

simulation = initialize!(model)
run!(simulation, period=Day(20))

The code above defines the Held-Suarez forcing ^[HS94] in terms of temperature relaxation and a linear drag term that is applied near the planetary boundary but switches off all other physics in the primitive equation model without humidity. Switching off the surface wind would also automatically turn off the surface evaporation (not relevant in the primitive dry model) and sensible heat flux as that one is proportional to the surface wind (which is zero with NoSurfaceWind). But to also avoid the calculation being run at all we use NoSurfaceHeatFlux() for the model constructor. Many of the NoSomething model components do not require the spectral grid to be passed on, but as a convention we allow every model component to have it for construction even if not required.

Visualising surface temperature with

using CairoMakie

temp = simulation.diagnostic_variables.grid.temp_grid[:, end]
heatmap(temp, title="Surface temperature [K]", colormap=:thermal)

Aquaplanet

using SpeedyWeather

# components
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
ocean = AquaPlanet(spectral_grid, temp_equator=302, temp_poles=273)
land_sea_mask = AquaPlanetMask(spectral_grid)
orography = NoOrography(spectral_grid)

# create model, initialize, run
model = PrimitiveWetModel(spectral_grid; ocean, land_sea_mask, orography)
simulation = initialize!(model)
run!(simulation, period=Day(20))

Here we have defined an aquaplanet simulation by

• creating an ocean::AquaPlanet. This will use constant sea surface temperatures that only vary with latitude.
• creating a land_sea_mask::AquaPlanetMask this will use a land-sea mask with false=ocean everywhere.
• creating an orography::NoOrography which will have no orography and zero surface geopotential.

All passed on to the model constructor for a PrimitiveWetModel, we have now a model with humidity and physics parameterization as they are defined by default (typing model will give you an overview of its components). We could have change the model.land and model.vegetation components too, but given the land-sea masks masks those contributions to the surface fluxes anyway, this is not necessary. Note that neither sea surface temperature, land-sea mask or orography have to agree. It is possible to have an ocean on top of a mountain. For an ocean grid-cell that is (partially) masked by the land-sea mask, its value will be (fractionally) ignored in the calculation of surface fluxes (potentially leading to a zero flux depending on land surface temperatures).

Now with the following we visualize the surface humidity after the 50 days of simulation. We use 50 days as without mountains it takes longer for the initial conditions to become unstable. The surface humidity shows small-scale patches in the tropics, which is a result of the convection scheme, causing updrafts and downdrafts in both humidity and temperature.

using CairoMakie

humid = simulation.diagnostic_variables.grid.humid_grid[:, end]
heatmap(humid, title="Surface specific humidity [kg/kg]", colormap=:oslo)

Aquaplanet without (deep) convection

Now we want to compare the previous simulation to a simulation without deep convection, called DryBettsMiller, because it is the Betts-Miller convection but with humidity set to zero in which case the convection is always non-precipitating shallow (because the missing latent heat release from condensation makes it shallower) convection. In fact, this convection is the default when using the PrimitiveDryModel. Instead of redefining the Aquaplanet setup again, we simply reuse these components spectral_grid, ocean, land_sea_mask and orography (because spectral_grid hasn't changed this is possible).

# Execute the code from Aquaplanet above first!
convection = DryBettsMiller(spectral_grid, time_scale=Hour(4))

# reuse other model components from before
model = PrimitiveWetModel(spectral_grid; ocean, land_sea_mask, orography, convection)

simulation = initialize!(model)
run!(simulation, period=Day(20))

humid = simulation.diagnostic_variables.grid.humid_grid[:, end]
heatmap(humid, title="No deep convection: Surface specific humidity [kg/kg]", colormap=:oslo)

But we also want to compare this to a setup where convection is completely disabled, i.e. convection = NoConvection() (many of the No model components don't require the spectral_grid to be passed on, but some do!)

# Execute the code from Aquaplanet above first!
convection = NoConvection(spectral_grid)

# reuse other model components from before
model = PrimitiveWetModel(spectral_grid; ocean, land_sea_mask, orography, convection)

simulation = initialize!(model)
run!(simulation, period=Day(20))

humid = simulation.diagnostic_variables.grid.humid_grid[:, end]
heatmap(humid, title="No convection: Surface specific humidity [kg/kg]", colormap=:oslo)

And the comparison looks like

Large-scale vs convective precipitation

using SpeedyWeather

# components
spectral_grid = SpectralGrid(trunc=31, nlayers=8)
large_scale_condensation = ImplicitCondensation(spectral_grid)
convection = SimplifiedBettsMiller(spectral_grid)

# create model, initialize, run
model = PrimitiveWetModel(spectral_grid; large_scale_condensation, convection)
simulation = initialize!(model)
run!(simulation, period=Day(10))

We run the default PrimitiveWetModel with ImplicitCondensation as large-scale condensation (see Implicit large-scale condensation) and the SimplifiedBettsMiller for convection (see Simplified Betts-Miller). These schemes have some additional parameters, we leave them as default for now, but you could do ImplicitCondensation(spectral_grid, relative_humidity_threshold = 0.8) to let it rain at 80% instead of 100% relative humidity. We now want to analyse the precipitation that comes from these parameterizations

using CairoMakie

(; precip_large_scale, precip_convection) = simulation.diagnostic_variables.physics
m2mm = 1000     # convert from [m] to [mm]
heatmap(m2mm*precip_large_scale, title="Large-scale precipiation [mm]: Accumulated over 10 days", colormap=:dense)

Precipitation (both large-scale and convective) are written into the simulation.diagnostic_variables.physics which, however, accumulate all precipitation during simulation. In the NetCDF output, precipitation rate (in mm/hr) is calculated from accumulated precipitation as a post-processing step. More interactively, you can also reset these accumulators and integrate for another 6 hours to get the precipitation only in that period.

# reset accumulators and simulate 6 hours
simulation.diagnostic_variables.physics.precip_large_scale .= 0
simulation.diagnostic_variables.physics.precip_convection .= 0
run!(simulation, period=Hour(6))

# visualise, precip_* arrays are flat copies, no need to read them out again!
m2mm_hr = (1000*Hour(1)/Hour(6))    # convert from [m] to [mm/hr]
heatmap(m2mm_hr*precip_large_scale, title="Large-scale precipiation [mm/hr]", colormap=:dense)
heatmap(m2mm_hr*precip_convection, title="Convective precipiation [mm/hr]", colormap=:dense)

As the precipitation fields are accumulated meters over the integration period we divide by 6 hours to get a precipitation rate $[m/s]$ but then multiply with 1 hour and 1000 to get the typical precipitation unit of $[mm/hr]$.

References