GlobalMeanEBM

GlobalMeanEBM

Submodule of ConceptualClimateModels which provides processes useful in creating global mean energy balance models. It is inspired by Budyko-Sellers-Ghil type models (without the latitudinal dependence). The mean does not have to be global, it can be hemispheric, or for other smaller regions.

BasicRadiationBalance(; T, f, ASR, OLR, c_T = 5e8)

Create the equation

\[c_T \frac{dT}{dt} = ASR - OLR + f\]

representing the most basic radiative energy balance at the top of the atmosphere setting a global mean temperature, see e.g., any introductory article (North et al., 1981; Ghil, 1981). ASR is the absorbed solar radiation, which defaults to S*(1 - α) in the default processes. S is the received insolation, by default equal to solar_constant, but could e.g., be any astronomical forcing such as AstronomicalForcingDeSaedeleer. α is the albedo and f any additional forcing such as CO2Forcing. OLR defaults to A + B*T. c_T is the heat capacity of the system in J/K/m². However, for convenience, the parameter added to the final equation is τ_T which is the timescale in seconds, i.e., c_T/solar_constant.

ΔTLinearRelaxation(; ΔT, T, τ = 5e6, A = 36.53, B = 0.658)

Create the equation

\[\tau_{\Delta T}\frac{\Delta T}{dt} = \Delta T_{ref}(T) - \Delta T\]

which exponentially relaxes the equator-to-pole temperature difference ΔT to its reference value $\Delta T_{ref}(T) = A - B*(T - 275.15)$, i.e., it decreases linearly with global mean temperature $T$ (in Kelvin). The default values for A, B are obtained from Equation (2) of (Gaskell et al., 2022). We also fitted paleoclimate data of (Osman et al., 2021) and found very similar results, A = 35.8, B = -1.11 for north hemisphere and A = 27.4, b = -0.513 for south.

Here ΔT is defined as the temperature difference between average temperatures at (0, 30) and (60, 90) latitudes. The timescale is taken as 2 months, although if τ = 0 is given, the equation $\Delta T ~ \Delta T_{ref}(T)$ is created instead.

ΔTStommelModel(; ΔT=ΔT, ΔS=ΔS, η1 = 2, η2 = 1, η3 = 0.3)

Create the equations

\[\dot{\Delta T} = \eta_1 - \Delta T - |\Delta T - \Delta S| \Delta T \dot{\Delta S} = \eta_2 - \eta_3\Delta S - |\Delta T - \Delta S| \Delta S\]

which are the two equations of the Stommel box model for Atlantic thermohaline circulation (Stommel, 1961), here presented in nondimensionalized form (Lohmann et al., 2021), so that temperature and sality are normalized by their coefficients $a_T, a_S$ relating them to the density of water

\[\rho = \rho_0 [1 - a_T(T - T_0) + a_S(S-S_0)]\]

for some reference values.

LinearOLR(; OLR, T, A = -277.0, B = 1.8)

Create the equation OLR ~ A + B*T. This is a linearized outgoing longwave radiation (OLR), and is the same equation as (7) of (North et al., 1981): $OLR = A + BT$ with $T$ temperature in Kelvin and $A, B$ constants. However, default $A, B$ are fitted from current CERES all sky OLR and using ERA5 data for the 2-meter temperature. We assume T in Kelvin. This linear approximation is quite accurate for temporally averaged data $T \in (220, 280)$ however drops drastically in accuracy after that due to the nonlinear effects of clouds (as evident by observational data).

(Koll and Cronin, 2018) provide a "proof" of the linearity of the clear sky OLR due to spectral properties of water vapor.

We note a big difference between current CERES data and the values reported in (North et al., 1981): here A=214.67 (assuming $T$ in Celcius) and B=1.8 versus the values A=203.3 and B=2.09 in (North et al., 1981).

If instead of all sky, if we fit the clear sky CERES data, we get A = -326.0, B = 2.09. Interestingly, coefficient B here is the same as that reported by (North et al., 1981), but A=244.88 (assuming T in Celcius) is not.

LinearClearSkyOLR(; kw...)

Equivalent with LinearOLR(; A = A = -326.0, B = 2.09, kw...) and provided as a convenience for the clear sky fit to CERES data.

EmissivityStefanBoltzmanOLR(; ε, T)

Create the equation OLR ~ ε*σ*T^4 where σ is the Stefan Boltzmann constant and ε the effective emissivity, also known as the "grayness" of the system, or the deviation it has from being a perfect black body (Ghil, 1981). ε then needs to be parameterized itself to include greenhouse or other climate effects.

EmissivityFeedbackTanh(; ε, Τ, left = 0.5, right = 0.4, rate = 0.5, Tref = 288.0)

Create an equation that assigns emissivity ε to hyperbolic tangent of temperature T. This is an ad-hoc feedback that was used in (Bastiaansen et al., 2023), similar to EmissivitySellers1969. In essence this is a TanhProcess with the given keywords as parameters.

EmissivitySellers1969(; ε, T, m = 0.5)

Create the equation ε ~ 1 - m*tanh(19*T^6*1e-16), which was used originally in (Sellers, 1969) to represent the effect of "water vapor, carbon dioxide, dust and clouds on terrestrial radiation".

SoedergrenClearSkyEmissivity(; ε, T, CO2, RH = 0.8, H_H20 = 2.0)

Create Eq. 10 of (Södergren et al., 2018), which is the same as Eq. 21 of (Barker and Ross, 1999) for the effective emissivity of clear sky atmosphere:

\[\varepsilon = 1 - \exp(0.082 - (2.38*0.1*e_s*RH*H_{H2O} + 40.3*CO2*1e-6)^{0.294})\]

with $e_s$ the saturation_vapor_pressure. The equation assumes CO2 concentration is in ppm and vapor pressure in kPa hence the conversion factors 0.1 and 1e-6.

DirectAlbedoAddition(; α, α_bg = 0.1, other_albedo_variables = (α_ice, α_clouds))

Create the equation α ~ α_bg + other_albedo_variables..., meaning that planetary albedo α is a direct sum of all specified albedos.

CoAlbedoProduct(; α, albedo_variables = (α_ice, α_cloud))

Create the equation 1 - α ~ prod(a -> (1 - a), albedo_variables) meaning that the co-albedo is the product of the co-albedos of all albedo variables. This would be e.g., the planetary albedo if all components were uniform layers, while the bottom-most layer (surface) had perfect absorption and all other layers had 0 absorption and finite reflection.

SeparatedClearAllSkyAlbedo(; α, α_cloud, C, α_clr = 0.15)

Create the equation α ~ α_cloud*C + α_clr*(1 - C).

(Bender et al., 2017) argue that one can assume a separation between clear-sky and cloud albedo, so that α = α_cloud*C + α_clr*(1 - C) with C the cloud fraction and α_clr the clear sky albedo. They further cite (Cess, 1976) to facilitate the claim Additionally, Eq. (20) of (Barker and Ross, 1999) provides an identical expression.

In most cases you want to provide a variable with its own process for α_clr.

IceAlbedoFeedback(; T, α_ice,
    max = 0.45, min = 0.1, Tscale = 10, Tfreeze = 275.15, τ = 0
)

Create an equation that assigns ice albedo α_ice to a hyperbolic tangent of temperature T. This represents an approximately linear decrease with T, as ice melts over part of the earth, while it is constant for all T for which the earth would be either entirely ice covered (T < Tfreeze - scale) or ice free (T > Tfreeze).

In essence this is a TanhProcess with the given keywords as parameters with reference temperature Tref = Tfreeze - scale/2.

This albedo is the most common used large-scale feedback in energy balance models, e.g., (Ghil, 1981), although it is typically taken as a piece-wise linear function. There is little change with using a hyperbolic tangent instead, while the tanh offers a differentiable flow.

The timescale τ if not zero will make an ExpRelaxation process relaxing to the hyperbolic tangent.

saturation_vapor_pressure(T)

Given surface temperature (in K) return saturation pressure for water vapor (in kPa) using the Tetens-Murray formula from (Wikipedia contributors, 2023), which is A*exp(B*T/(C+T)). Different formula is used for when T is less than the freezing point.

AstronomicalForcingDeSaedeleer(; S = S, extensive = false)

Create the equation S ~ astronomical_forcing_desaedeleer(t, extensive) which is Eq. (1) of (de Saedeleer et al., 2013):

S = \sum_i s_i \sin(\omega_i t) + c_i \cos(\omega_i t)

where the values of $\omega_i, s_i, c_i$ come from (Berger, 1978) who performed a spectral expansion of the insolation. The validity range of this approximation is [-1, 0] Myr.

In the summation $i$ goes up to 35 if extensive, otherwise up to 8. The components are sorted according to magnitude of the spectral line, so the default version has only the 8 most important spectra lines.

Note that in contrast to Eq. (1) of (de Saedeleer et al., 2013) we do not normalize $f$ and its value is in W/m² (the mean value is still deducted). Additionally, the values of $\omega_i$ have been adjusted to expect time in units of seconds.

CO2Forcing(; f, CO2, CO2f = 3.7)

Create the equation $f ~ CO2f \log_2(CO2/400)$ which describes the forcing added to the TOA energy balance due to CO2 concentrations, assumming the OLR expression is calibrated for 0 added forcing at 400 ppm which is the default for OLR expressions provided by ConceptualClimateModels.jl.

The default value of $f$ comes from Eq. (3.2) of (Bastiaansen et al., 2023) which cites IPCC-5, while (Etminan et al., 2016) report practically the same value assuming a constant $f$ (note here the log is base 2). In reality $f$ depends on $CO2$ and other greenhouse gases concentrations due to spectral overlaps, see (Etminan et al., 2016) Sec. 4.

CloudAlbedoExponential(
    α_cloud, C, a = 2.499219232848238, b = 17.596369331717433
)

Create the equation α_cloud ~ sinh(a*C)/b relating cloud albedo to cloud fraction C. This equation is exponential and not linear, as in observations. (Engström et al., 2015) (and also (Bender et al., 2017)) discuss this exponential relation in detail, and provide as explanation that cloud effective albedo increases with latitude (due to solar zenith changes) while cloud fraction also increases with latitude.

Note that here however we modify the equation α_cloud ~ exp(a*C - b) of (Engström et al., 2015) to utilize the hyperbolic sine, so that α_cloud = 0 when C = 0 as is physically necessary. Then, a, b are extracted by fitting CERES data, using as α_cloud the energetically consistent cloud albedo as defined by (Datseris and Stevens, 2021), further yearly averaged and within latitudes (-60, 60) as in (Bender et al., 2017). This albedo can be directly added to the clear sky albedo to produce the planetary albedo.

CloudAlbedoLinear(; α_cloud, C, a = 0.2252861764703435)

Same as in CloudAlbedoExponential, but now the linear form α_cloud ~ a*C is returned, with a fitted from CERES data in the same way.

BudykoOLR(; OLR=OLR, T=T, C=C,
    BudykoOLR_A = -461.8068, BudykoOLR_B = 2.58978,
    BudykoOLR_Ac = -377.22741, BudykoOLR_Bc = 1.536171
)

Create the equation OLR ~ A + B*T - C*(Ac + Bc*T) for the dependence of OLR on both temperature and cloud fraction (in 0-1). This is the same as Eq. (1) of (Budyko, 1969). However, here T is expected in Kelvin, and the coefficients have been extracted by fitting into CERES data in the same way as in LinearOLR.