Surface fluxes

The surfaces fluxes in SpeedyWeather represent the exchange of momentum, heat, and humidity/moisture between ocean and land as surface into or out of the lowermost atmospheric layer (see Vertical coordinates for how vertical layers are indexed). Surface fluxes of momentum represent a drag that the boundary layer wind experiences due to friction over more or less rough ground on land or over sea. Surface fluxes of heat represent a sensible heat flux from a warmer or colder ocean or land into or out of the surface layer of the atmosphere. Surface fluxes of humidity represent evaporation of sea water or condensation of humidity onto the ocean surface. Similarly, humidity fluxes over land with a given soil moisture and vegetation's evapotranspiration.

Surface flux implementations

Currently implemented surface fluxes of momentum are

using SpeedyWeather
subtypes(SpeedyWeather.AbstractSurfaceMomentumFlux)

1-element Vector{Any}:
 SurfaceMomentumFlux

Interdependence of surface flux computations

model.boundary_layer.surface_condition extrapolates atmospheric variables (wind, temperature, humidity) to the surface. model.boundary_layer computes centrally a drag coefficient for surface momentum, heat and humidity fluxes. Setting model.boundary_layer.surface_condition = nothing or model.boundary_layer = nothing will therefore disable all other surface fluxes unless they have use_boundary_layer_drag = false which lets them use independently their own drag coefficient.

with more explanation below. The surface heat fluxes currently implemented are

subtypes(SpeedyWeather.AbstractSurfaceHeatFlux)

5-element Vector{Any}:
 PrescribedLandHeatFlux
 PrescribedOceanHeatFlux
 SurfaceHeatFlux
 SurfaceLandHeatFlux
 SurfaceOceanHeatFlux

and the surface humidity fluxes (this does not include Convection or Large-scale condensation those are precipitation) implemented are

subtypes(SpeedyWeather.AbstractSurfaceHumidityFlux)

5-element Vector{Any}:
 PrescribedLandHumidityFlux
 PrescribedOceanHumidityFlux
 SurfaceHumidityFlux
 SurfaceLandHumidityFlux
 SurfaceOceanHumidityFlux

The calculation of thermodynamic quantities at the surface (air density, temperature, humidity) and the computation of drag coefficients (which is used by default for the surface fluxes above) are handled through the model.boundary_layer where currently implemented are

subtypes(SpeedyWeather.AbstractBoundaryLayer)

1-element Vector{Any}:
 BoundaryLayer

See Parameterizations for the general parameterization interface these implement.

Boundary layer

The model.boundary_layer is a composite of four components that run one after another for every grid point, each writing variables the next one (and the surface fluxes) read

BoundaryLayer
├── surface_condition    # 1. surface air temperature, density, wind speed
│   ├── wind_slowdown    #    ratio of near-surface to lowest-layer wind
│   └── wind_gust        #    sub-grid scale gusts added to the surface wind speed
├── neutral_wind_speed   # 2. stability-corrected (neutral) surface wind speed
├── surface_roughness    # 3. roughness lengths for momentum, heat, humidity over land and ocean
│   ├── land_surface_roughness
│   └── ocean_surface_roughness
│       ├── charnock
│       ├── sea_ice_momentum_roughness
│       └── max_ocean_momentum_roughness
└── drag                 # 4. drag coefficients for momentum, heat, humidity over land and ocean

Circularity of wind slowdown and ocean roughness

The LearnedWindSlowdown (default) reduces the lowermost-layer wind towards the surface depending on the momentum roughness length : rougher surfaces slow the wind more. Over the ocean, however, the roughness is itself a function of the wind: the LearnedOceanRoughness (default) computes from the neutral wind speed, which is derived from the surface wind. So the surface wind depends on which depends on the surface wind (surface_condition → neutral_wind_speed → surface_roughness → surface_condition). Over land there is no such circularity, as the land roughness depends on vegetation, snow and soil only for now.

There are several ways to resolve this circularity, each with a downside:

• Lag: use the roughness from the previous time step. This introduces a hidden state that is not part of the prognostic variables (it is not, e.g., written to restart files), needs an ad-hoc roughness on the first time step, and can oscillate in time at low wind speeds, where the learned ocean roughness drops by many orders of magnitude.

• Iterate: repeat surface condition → neutral wind speed → surface roughness within a time step until converged. Stateless and consistent, but this multiplies the cost of these parameterizations and converges slowly at the same low-wind points.

• Decouple: use fixed roughness lengths in the wind slowdown only. Stateless and cheap, but the wind slowdown does not see the actual roughness.

The default is to lag. We set LearnedWindSlowdown with use_computed_roughness = true which uses the roughness from the previous timestep, and we fall back to the fixed roughness length on the first time step.

To decouple, one sets use_computed_roughness = false, and the LearnedWindSlowdown uses fixed roughness lengths over land and ocean (land_roughness = 0.1 m, ocean_roughness = 1e-4 m, weighted by the land fraction) so that the surface condition depends only on the current state. The computed roughness is still used for the drag coefficients.

The iterative solution is currently not implemented.

The subcomponents of SurfaceCondition are

subtypes(SpeedyWeather.AbstractWindSlowdown)

2-element Vector{Any}:
 ConstantWindSlowdown
 LearnedWindSlowdown

subtypes(SpeedyWeather.AbstractWindGust)

2-element Vector{Any}:
 BeljaarsWindGust
 ConstantWindGust

the neutral wind speed

subtypes(SpeedyWeather.AbstractNeutralWindSpeed)

1-element Vector{Any}:
 NeutralWindSpeed

which can also be set to nothing (BoundaryLayer(spectral_grid, neutral_wind_speed = nothing)). Parameterizations that require a neutral wind speed (e.g. LearnedOceanRoughness) then use the (non-neutral) surface wind speed instead, i.e. without stability correction.

the surface roughness, with its land and ocean subcomponents (used by LearnedSurfaceRoughness)

subtypes(SpeedyWeather.AbstractSurfaceRoughness)

4-element Vector{Any}:
 ConstantSurfaceRoughness
 LearnedSurfaceRoughness
 SpeedyWeather.AbstractLandRoughness
 SpeedyWeather.AbstractOceanRoughness

subtypes(SpeedyWeather.AbstractLandRoughness)

1-element Vector{Any}:
 LearnedLandRoughness

subtypes(SpeedyWeather.AbstractOceanRoughness)

1-element Vector{Any}:
 LearnedOceanRoughness

and the functors used within the learned ocean roughness

subtypes(SpeedyWeather.AbstractRoughnessFunctor)

3-element Vector{Any}:
 Charnock
 MaxOceanMomentumRoughness
 SeaIceMomentumRoughness

and finally the drag coefficients

subtypes(SpeedyWeather.AbstractBoundaryLayerDrag)

2-element Vector{Any}:
 BulkRichardsonDrag
 ConstantDrag

All of these take the spectral_grid as first argument and their parameters (or subcomponents) as keyword arguments, so that they can be composed together and passed on to the model constructor, e.g.

spectral_grid = SpectralGrid(truncation = 32, nlayers = 8)

# constant wind slowdown and gusts instead of the defaults
surface_condition = SurfaceCondition(
    spectral_grid,
    wind_slowdown = ConstantWindSlowdown(spectral_grid, wind_slowdown = 0.9),
    wind_gust = ConstantWindGust(spectral_grid, wind_gust = 1),
)

# learned roughness but with a different Charnock parameterization
ocean_surface_roughness = LearnedOceanRoughness(
    spectral_grid,
    charnock = Charnock(spectral_grid, charnock_min = 1.0e-4),
)
surface_roughness = LearnedSurfaceRoughness(spectral_grid; ocean_surface_roughness)

# or constant roughness lengths instead
# surface_roughness = ConstantSurfaceRoughness(spectral_grid, momentum_roughness_ocean = 2e-4)

boundary_layer = BoundaryLayer(
    spectral_grid;
    surface_condition,
    surface_roughness,
    drag = BulkRichardsonDrag(spectral_grid, critical_Richardson = 5),
)

model = PrimitiveWetModel(spectral_grid; boundary_layer)
model.boundary_layer

BoundaryLayer{SurfaceCondition{ConstantWind...} <: SpeedyWeather.AbstractBoundaryLayer
├ surface_condition::SurfaceCondition{ConstantWindSlowdown{Float32}, Constant...
├ neutral_wind_speed::NeutralWindSpeed{Float32}
├ surface_roughness::LearnedSurfaceRoughness{LearnedLandRoughness{Float32}, L...
└ drag::BulkRichardsonDrag{Float32}

Fluxes to tendencies

In SpeedyWeather.jl, parameterizations can be defined either in terms of tendencies for a given layer or as fluxes between two layers including the surface flux and a top-of-the-atmosphere flux. The upward flux out of layer   into layer (vertical indexing increases downwards) is (  for half, as the flux sits on the cell face, the half-layer in between and  ) and similarly is the downward flux. For clarity we may define fluxes as either upward or downward depending on the process although an upward flux can always be regarded as a negative downward flux. The absorbed flux in layer is

A quick overview of the units

Quantity	Variable	Unit	Flux unit
Momentum	Velocity
Heat	Temperature
Moisture	Specific humidity

The time-stepping in SpeedyWeather.jl (see Time integration) eventually requires tendencies which are calculated from the absorbed fluxes of momentum or and moisture as

with gravity and layer-thickness (see Sigma coordinates) so that the right-hand side divides the absorbed flux by the mass of layer (per unit area). Tendencies for equivalently with their respective absorbed fluxes.

The temperature tendency is calculated from the absorbed heat flux as

with heat capacity of air at constant pressure with a typical value of . Because we define the heat flux as having units of the conversion includes the division by the heat capacity to convert to a rate of temperature change.

Time integration of surface fluxes

Surface fluxes as a relaxation

The bulk surface fluxes of momentum, heat and humidity below all relax a variable of the lowermost layer ( or ) towards a surface value ( for momentum, the skin or sea surface temperature for heat, the saturation humidity at the surface for humidity). Each of them can be written as

with the conductance    in . is a mass exchange rate: the mass of air per square metre and second that effectively exchanges its momentum, heat or moisture with the surface. contains the surface air density , the (dimensionless) drag coefficient and the surface wind speed , it is additionally scaled by for momentum, and divided by the sea-ice or snow insulation factors for heat and humidity. All fluxes are positive upward.

Flux			units
Momentum	[]
Heat	[]	or SST	, see below	/ insulation
Humidity	[]	or		/ insulation

with the saturation humidity and the soil moisture availability over land. Note that only scales the surface value , not , so it does not enter .

A flux through the bottom of layer with pressure thickness changes the content of that layer. The layer has a mass of per square metre, so its momentum or moisture content is per square metre and therefore

The relaxation time scale is the time it takes to exchange the mass of the lowermost layer with the surface at a rate , which is the same for momentum, heat and humidity if they use the same drag coefficient. is short for a thin lowermost layer (small ), strong winds or high drag. With FriersonSigmaCoordinates for example, the lowermost layer is about 3 times thinner than with equi-distant levels, which can reduce below the time step.

Heat fluxes and the heat capacity

The sensible heat flux is defined in as    , so that   in , not . However, the heat content of the lowermost layer is per square metre, so the temperature tendency is

and cancels: it appears both in the flux and in the heat content of the layer. The relaxation time scale   is therefore the same as above and does not depend on . In the code the heat flux is accordingly calculated without , i.e. as     in  , for the time integration and only multiplied by when stored in .

Exact solution and flux damping

Holding and fixed over a time step ( for leapfrog, as the tendency is applied from   to  ), the relaxation equation has the exact solution of an exponential decay towards the surface value

How does this compare to what the model actually does? The time stepping in SpeedyWeather is explicit: all tendencies, including those from the surface fluxes, are calculated from the state at the beginning of the time step and then used to step forward. For the surface flux alone this is a forward Euler step

with   the value at the beginning of the time step and the numerical approximation of the exact  . Using     this can be written as a change in the deviation from the surface value

So within one time step the explicit step multiplies the deviation   by  , whereas the exact solution multiplies it by . The two agree for small , as     , but they differ more and more for larger : for   the factor   is negative, so overshoots past (e.g. the wind reverses direction or the air becomes colder than the surface below it), and for   we have   , so the deviation grows with every time step: the explicit step is unstable even though the true solution simply decays.

Can we get the exact solution but keep the explicit time stepping? Compare the change in over one time step from both

With   both changes are proportional to the same deviation  , they only differ by the factor

which only depends on . So if we replace the flux by the explicit step yields the exact solution,   , without changing the time stepping itself:

There is also a physical interpretation of : as relaxes towards during the time step, the flux decays exponentially too,  . Its average over the time step is

So is the time-mean flux over the time step, whereas the explicit step uses the flux at the beginning of the time step, which is too large when the flux decays quickly. is therefore a damping factor   , close to when the time step resolves the relaxation ( ) and small when it does not ( ). This approach of integrating the linear part of an equation exactly is known as exponential time differencing. It has two practical advantages here: is calculated pointwise from , so no iteration or system of equations has to be solved, and the damped flux is a flux itself, which can be passed on to the ocean and land models (see below).

For small we have     , and   for large , meaning is relaxed fully to within one time step.

However, evaluating requires an exponential for every grid point, every surface flux and every time step. This is relatively expensive: at T32L8 using ExponentialFluxes for all surface fluxes increases the time per time step by about 5%, compared to about 1% for ImplicitFluxes. We therefore look for cheaper approximations of that keep its essential properties:   for small (the flux is unchanged if the time step resolves the relaxation) and   for large ( relaxes to but does not overshoot it). The other time integration schemes correspond to such approximations of the exact :

• Explicit,  , only approximates the exact solution to first order in .   replaces , which overshoots the surface value for   and is unstable (  ) for  .

• Implicit (backward Euler) evaluates the flux with at the end of the time step,      , so that   . This is unconditionally stable and monotone, and has the correct behaviour for large , but is only first-order accurate,   , so it overdamps for small .

• Rational,     , matches the expansion of the exact to second order,   , and has the same behaviour for large , at the cost of the implicit scheme (one division) instead of an exp.   becomes negative for   with a minimum of about near  , meaning that overshoots the surface value by at most 10% of the initial deviation, then tends to zero for large .

Scheme			Single step	Cost
ExplicitFluxes()			overshoots for  , unstable for	none
ImplicitFluxes()			monotone, stable, first order	1 division
ExponentialFluxes()			monotone, stable, exact for fixed	exp, ~5% of a time step at T32L8
RationalFluxes() (default)			stable, ≤10% overshoot for	1 division

Because is applied to the flux (not only to the tendency) before it is stored, the fluxes passed on to the ocean and land models and to the output are consistent with those the atmosphere receives, conserving momentum, energy and water between the components. For the land and ocean fluxes, is weighted by the land or ocean fraction respectively, as only that fraction of the flux enters the lowermost layer of the atmosphere. The drag coefficient , the wind speed and the surface value are held fixed over the time step (linearisation), which makes and constant within a time step. Prescribed fluxes (e.g. PrescribedOceanHeatFlux) are applied as given.

Choosing the time integration

The scheme is chosen per flux, for example

spectral_grid = SpectralGrid(truncation = 31, nlayers = 8)
surface_momentum_flux = SurfaceMomentumFlux(spectral_grid, time_integration = ExponentialFluxes())
surface_heat_flux = SurfaceHeatFlux(
    spectral_grid,
    ocean = SurfaceOceanHeatFlux(spectral_grid, time_integration = ImplicitFluxes()),
    land = SurfaceLandHeatFlux(spectral_grid, time_integration = ExplicitFluxes()),
)
model = PrimitiveWetModel(spectral_grid; surface_momentum_flux, surface_heat_flux)

and similarly for SurfaceOceanHumidityFlux and SurfaceLandHumidityFlux.

Bulk Richardson-based drag coefficient

All surface fluxes depend on a dimensionless drag coefficient which we calculate as a function of the bulk Richardson number (see also Vertical diffusion for a similar Richardson-number-based diffusion coefficient) following with some simplification as outlined below. We use different drag coefficients for momentum, heat and moisture fluxes. The bulk Richardson number at the lowermost model layer   of height is

with   the Geopotential at  ,    the virtual dry static energy, the Virtual temperature and the surface wind speed including gusts (see Surface momentum fluxes).

The virtual dry static energy at the surface   is

computed separately over land and ocean. Over ocean, is the sea surface temperature and  . Over land, is the top soil layer temperature and  , with the soil moisture availability as in the land humidity flux. In the dry model  .

Then the drag coefficient follows as

with the neutral drag coefficient and the minimum drag defined below, and roughness lengths computed in the surface roughness scheme boundary_layer.surface_roughness, per surface and per flux type. There is a maximum drag for negative bulk Richardson numbers but the drag drops to for bulk Richardson numbers being larger than a critical   with a smooth transition in between. The height of the -th model layer is   with the geopotential

which depends on the temperature of that layer. This model layer height above the surface is recalculated on every time step for every grid cell as it depends on temperature. The maximum drag coefficient then follows as

with   the von Kármán constant, the momentum roughness length and the roughness length of momentum, heat or humidity, so the three coefficients share the stability dependence but differ in their roughness lengths. This simplifies the drag coefficient calculation to

with   the clamped which is at least and at most , and   (drag_min) a minimum drag that avoids vanishing surface fluxes in strongly stable conditions.

Surface momentum fluxes

The surface momentum flux is calculated from the surface wind velocities

meaning it is scaled down by   (Fortran SPEEDY default; ()) from the lowermost layer wind velocities . A wind speed scale accounting for gustiness with   (Fortran SPEEDY default; ()) is then defined as

such that for low wind speeds the fluxes are somewhat higher because of unresolved winds on smaller time and length scales. The momentum flux is then

with   the surface air density calculated from surface pressure and lowermost layer temperature . Better would be to extrapolate to a surface air temperature assuming adiabatic descent but that is currently not implemented. How this flux is integrated in time is described in Time integration of surface fluxes.

Surface heat fluxes

The surface heat flux is proportional to the difference of the surface air temperature and the land or ocean skin temperature . We currently just approximate as the lowermost layer temperature although an adiabatic descent from pressure to surface pressure would be more accurate. We also use land and sea surface temperature to approximate although future improvements should account for faster radiative effects on compared to sea and land surface temperatures determined by a higher heat capacity of the relevant land surface layer or the mixed layer in the ocean. We then compute

in , positive up. How this flux is integrated in time is described in Time integration of surface fluxes.

Surface humidity fluxes

The surface humidity flux, i.e. evaporation/condensation of soil moisture over land and evaporation of sea water over the ocean or condensation of humidity on the ocean surface is proportional to the difference of the surface specific humidity and the saturation specific humidity given and surface pressure . This assumes that a very thin layer of air just above the ocean is saturated but over land this assumption is less well justified as it should be a function of the soil moisture and how much of that is available to evaporate given vegetation, hence we include a soil moisture availability. We again make the simplification that  , i.e. specific humidity of the surface is the same as in the lowermost atmospheric layer above.

The surface humidity flux in , positive up, is then

with the saturation specific humidity calculated from the skin temperature and surface pressure . The availability of soil water over land is represented by (over the ocean  )

following the Fortran SPEEDY documentation (), which follows . The variables (or spatially prescribed arrays) are water content in the top soil layer and the root layer below (see Land soil moisture for how these are modelled) using the vegetation fraction    composed of a (dimensionless) high and low vegetation cover per grid cell . The constants are depth of top soil layer  , depth of root layer   (defaults subject to change, check model.land.geometry), soil wetness at field capacity (volume fraction)  , and soil wetness at wilting point (volume fraction)  . How this flux is integrated in time is described in Time integration of surface fluxes.

Land-sea mask

SpeedyWeather uses a fractional land-sea mask, i.e. for every grid-point

• 1 indicates land

• 0 indicates ocean

• a value in between indicates a grid-cell partially covered by ocean and land

The land-sea mask determines solely how to weight the surface fluxes coming from land or from the ocean. For the sensible heat fluxes this uses land and sea surface temperatures and weights the respective fluxes proportional to the fractional mask. Similar for humidity fluxes. You can therefore define an ocean on top of a mountain, or vice versa. Let be the fluxes coming from land and sea, respectively. Then the land-sea mask   weights them as follows for the total flux

Note that this weighting is purely multiplicative: it does not protect against NaN propagating into should or be NaN. Land and sea surface temperatures (and other quantities entering the flux calculation, such as soil moisture) must therefore be finite everywhere, even where the land-sea mask is 0 or 1, i.e. a NaN land or sea surface temperature is not supported anymore. Use a fallback temperature for grid cells without "real" land or ocean instead (e.g. SlabOcean's land_temperature, or SeasonalLandTemperature's ocean_temperature).

Setting the land-sea mask to ocean therefore will disable any fluxes that may come from land, and vice versa. However, with an ocean-everywhere land-sea mask you must still also define sea surface temperatures everywhere, otherwise the fluxes in those regions will be incorrect rather than simply zero.

For more details see The land-sea mask implementation section.