Merge branch 'thornhill-skeie'

CHANGELOG.rst

@@ -4,12 +4,23 @@ Changelog
 master
 ------
 
+v1.6.2
+------
+
+(`#99 <https://github.com/OMS-NetZero/FAIR/pull/99>`_) Add option to use ozone forcing dependent on N2O concentrations
+
+v1.6.1
+------
+
 (`#87 <https://github.com/OMS-NetZero/FAIR/pull/87>`_) Add Geoffroy temperature and prescribed forcing to inverse
 
 (`#84 <https://github.com/OMS-NetZero/FAIR/pull/84>`_) Fix concentration-driven runs for 45-species FaIR
 
 (`#83 <https://github.com/OMS-NetZero/FAIR/pull/83>`_) Support scmdata >= 0.7.1
 
+earlier
+-------
+
 (`#82 <https://github.com/OMS-NetZero/FAIR/pull/82>`_) Expand AR6 forcing diagnostics to get aerosol direct forcing by species
 
 (`#78 <https://github.com/OMS-NetZero/FAIR/pull/78>`_) Optimise ``fair.tools.scmdf.scmdf_to_emissions`` a bit

fair/forcing/ozone.py

@@ -0,0 +1,95 @@
+from __future__ import division
+
+import numpy as np
+from ..constants import cl_atoms, br_atoms, fracrel
+
+
+def thornhill_skeie(
+        emissions,
+        concentrations,
+        temperature=0,
+        feedback=-0.037,
+        beta=np.array([2.33379720e-04,  1.27179106e-03, -6.69347820e-05,
+                       1.14647701e-04,  5.14366051e-12,  3.78354423e-03]),
+        emissions_pi=np.zeros(40),
+        concentrations_pi=np.zeros(31),
+    ):
+    """Calculates total ozone forcing from precursor emissions and
+    concentrations based on AerChemMIP and CMIP6 Historical behaviour
+
+    Skeie et al. (2020)
+    Thornhill et al. (2021)
+
+    Unlike Stevenson CMIP5 no distinction is made for tropospheric and
+    stratospheric.
+
+    With this formulation, ozone forcing depends on concentrations of
+    CH4, N2O, ozone-depleting halogens, and emissions of CO, NVMOC and NOx,
+    but any combination of emissions and concentrations are allowed.
+
+    Inputs:
+        emissions: (nt x 40) numpy array in FaIR default units
+        concentrations: (nt x 31) numpy array of GHGs in FaIR default units
+        temperature: global mean surface temperature (for feedback)
+        feedback: temperature feedback on ozone forcing (W/m2/K) - set to zero
+            or False to turn off
+        beta: 6-element array of radiative efficiency coefficients in order of
+            CH4 concentrations,  W m-2 ppb-1 
+            N2O concentrations,  W m-2 ppb-1
+            ODS concentrations in EESC,  W m-2 ppt-1
+            CO emissions, W m-2 (Mt yr-1)-1
+            NMVOC emissions, W m-2 (Mt yr-1)-1
+            NOx emissions, W m-2 (MtN yr-1)-1
+
+        emissions_pi: pre-industrial/reference emissions
+        concentrations_pi: pre-industrial/reference concentrations
+
+    Outputs:
+        ozone ERF time series.
+    """
+
+    # we allow 2D output for quick calculation if feedbacks turned off
+    if emissions.ndim == 1:
+        nspec = len(emissions)
+        emissions = emissions.reshape((1, nspec))
+    if concentrations.ndim == 1:
+        nspec = len(concentrations)
+        concentrations = concentrations.reshape((1, nspec))
+
+    nt = emissions.shape[0]
+
+    # calculate EESC for halogens
+    cl = np.array(cl_atoms.aslist)
+    br = np.array(br_atoms.aslist)
+    fc = np.array(fracrel.aslist)
+
+    def eesc(c_ods, c_ods_pi):
+        return (
+            np.sum(cl * (c_ods-c_ods_pi) * fc/fc[0]) + 
+            45 * np.sum(br * (c_ods-c_ods_pi) * fc/fc[0])
+        ) * fc[0]
+
+
+    c_ch4, c_n2o = concentrations[:, [1, 2]].T
+#    delta_c_ods = eesc(concentrations[:,15:].T, concentrations_pi[None, 15:])
+    c_ods = concentrations[:,15:]
+    e_co, e_nmvoc, e_nox = emissions[:,[6, 7, 8]].T
+    c_ch4_pi, c_n2o_pi = concentrations_pi[[1, 2]]
+    c_ods_pi = concentrations_pi[15:]
+    e_co_pi, e_nmvoc_pi, e_nox_pi = emissions_pi[[6, 7, 8]]
+
+
+    forcing = np.zeros(nt)
+    if np.isscalar(temperature):
+        temperature = np.ones(nt) * temperature
+
+    for i in range(nt):
+        f_ch4   = beta[0] * (c_ch4[i] - c_ch4_pi)
+        f_n2o   = beta[1] * (c_n2o[i] - c_n2o_pi)
+        f_ods   = beta[2] * eesc(c_ods[i], c_ods_pi)
+        f_co    = beta[3] * (e_co[i] - e_co_pi)
+        f_nmvoc = beta[4] * (e_nmvoc[i] - e_nmvoc_pi)
+        f_nox   = beta[5] * (e_nox[i] - e_nox_pi)
+        forcing[i] = f_ch4 + f_n2o + f_ods + f_co + f_nmvoc + f_nox + feedback * temperature[i]
+
+    return forcing

fair/forcing/ozone_tr.py

@@ -87,7 +87,6 @@ def temperature_feedback(T, a=0.03189267, b=1.34966941, c=-0.03214807):
         F = F_CH4 + F_CO + F_NMVOC + F_NOx + temperature_feedback(T)
     else:
         F = F_CH4 + F_CO + F_NMVOC + F_NOx
-
     return F
 
 

fair/forward.py

@@ -8,8 +8,8 @@
 from .constants import molwt, lifetime, radeff
 from .constants.general import M_ATMOS, ppm_gtc
 from .defaults import carbon, thermal
-from .forcing import ozone_tr, ozone_st, h2o_st, contrails, aerosols, bc_snow,\
-                                         landuse
+from .forcing import (ozone, ozone_tr, ozone_st, h2o_st, contrails, aerosols,
+    bc_snow, landuse)
 from .gas_cycle.gir import calculate_alpha, step_concentration
 from .gas_cycle.fair1 import carbon_cycle
 from .forcing.ghg import co2_log, minor_gases
@@ -90,6 +90,7 @@ def fair_scm(
     ref_isSO2=True, # is Stevens SO2 emissions in units SO2 (T) or S (F)
     useMultigas=True,
     tropO3_forcing='stevenson',
+    ozone_feedback=-0.037, # W m-2 K-1, only for Thornhill-Skeie
     lifetimes=False,
     aerosol_forcing="aerocom+ghan",
     scaleAerosolAR5=True,
@@ -441,27 +442,38 @@ def fair_scm(
         # because SLCFs can still be given as emissions with GHGs as
         # concentrations
         if type(emissions) is not bool:
-            if tropO3_forcing[0].lower()=='s':
+            if tropO3_forcing[0].lower()=='s':  # stevenson
                 F[0,iF_tro3] = ozone_tr.stevenson(emissions[0,:], C[0,1],
                   T=np.sum(T_j[0,:]), 
                   feedback=useTropO3TFeedback,
                   fix_pre1850_RCP=fixPre1850RCP,
                   PI=pi_tro3)
-            elif tropO3_forcing[0].lower()=='c':
+            elif tropO3_forcing[0].lower()=='c':  # cmip6 stevenson
                 F[0,iF_tro3] = ozone_tr.cmip6_stevenson(emissions[0,:], C[0,1],
                   T=np.sum(T_j[0,:]),
                   feedback=useTropO3TFeedback,
                   PI=np.array([C_pi[1],E_pi[6],E_pi[7],E_pi[8]]),
                   beta=b_tro3)
-            elif tropO3_forcing[0].lower()=='r':
+            elif tropO3_forcing[0].lower()=='r':  # regression
                 F[0,iF_tro3] = ozone_tr.regress(emissions[0,:]-E_pi[:], beta=b_tro3)
+            elif tropO3_forcing[0].lower()=='t':  # thornhill-skeie
+                F[0,iF_tro3] = ozone.thornhill_skeie(
+                    emissions=emissions[0,:],
+                    concentrations=C[0,:],
+                    temperature=T[0],
+                    feedback=ozone_feedback,
+                    beta=b_tro3,
+                    emissions_pi=E_pi,
+                    concentrations_pi=C_pi,
+                )
             else:
                 F[0,iF_tro3] = F_tropO3[0]
         else:
             F[0,iF_tro3] = F_tropO3[0]
 
-        # Stratospheric ozone depends on concentrations of ODSs (index 15-30)
-        F[0,iF_sto3] = ozone_st.magicc(C[0,15:], C_pi[15:])
+        if tropO3_forcing[0].lower()!='t':
+            # Stratospheric ozone depends on concentrations of ODSs (index 15-30)
+            F[0,iF_sto3] = ozone_st.magicc(C[0,15:], C_pi[15:])
 
         # Stratospheric water vapour is a function of the methane ERF
         F[0,iF_ch4h] = h2o_st.linear(F[0,1], ratio=stwv_from_ch4)
@@ -711,9 +723,21 @@ def fair_scm(
                       beta=b_tro3)
                 elif tropO3_forcing[0].lower()=='r':
                     F[t,iF_tro3] = ozone_tr.regress(emissions[t,:]-E_pi, beta=b_tro3)
+                elif tropO3_forcing[0].lower()=='t':
+                    F[t,iF_tro3] = ozone.thornhill_skeie(
+                        emissions=emissions[t,:],
+                        concentrations=C[t,:],
+                        temperature=T[t-1],
+                        feedback=ozone_feedback,
+                        beta=b_tro3,
+                        emissions_pi=E_pi,
+                        concentrations_pi=C_pi,
+                    )
                 else:
                     F[t,iF_tro3] = F_tropO3[t]
-                F[t,iF_sto3] = ozone_st.magicc(C[t,15:], C_pi[15:])
+
+                if tropO3_forcing[0].lower()!='t':
+                    F[t,iF_sto3] = ozone_st.magicc(C[t,15:], C_pi[15:])
                 F[t,iF_ch4h] = h2o_st.linear(F[t,1], ratio=stwv_from_ch4)
 
                 # multiply by scale factors
@@ -824,11 +848,23 @@ def fair_scm(
                           beta=b_tro3)
                     elif tropO3_forcing[0].lower()=='r':
                         F[t,iF_tro3] = ozone_tr.regress(emissions[t,:]-E_pi, beta=b_tro3)
+                    elif tropO3_forcing[0].lower()=='t':
+                        F[t,iF_tro3] = ozone.thornhill_skeie(
+                            emissions=emissions[t,:],
+                            concentrations=C[t,:],
+                            temperature=T[t-1],
+                            feedback=ozone_feedback,
+                            beta=b_tro3,
+                            emissions_pi=E_pi,
+                            concentrations_pi=C_pi,
+                        )
                     else:
                         F[t,iF_tro3] = F_tropO3[t]
                 else:
                     F[t,iF_tro3] = F_tropO3[t]
-                F[t,iF_sto3] = ozone_st.magicc(C[t,15:], C_pi[15:])
+
+                if tropO3_forcing[0].lower()!='t':
+                    F[t,iF_sto3] = ozone_st.magicc(C[t,15:], C_pi[15:])
                 F[t,iF_ch4h] = h2o_st.linear(F[t,1], ratio=stwv_from_ch4)
 
                 # multiply by scale factors

tests/unit/unit_test.py

@@ -8,6 +8,7 @@
 from fair.defaults import carbon
 from fair.ancil import cmip5_annex2_forcing
 from fair.forcing.ozone_tr import regress
+from fair.forcing.ozone import thornhill_skeie
 from fair.inverse import inverse_fair_scm
 import numpy as np
 import os
@@ -523,3 +524,14 @@ def test_inverse_millar_fin():
         temperature_function = 'Millar',
         F_in = rcp85.Forcing.total
     )
+
+
+def test_thornhill_skeie():
+    forcing = thornhill_skeie(
+        emissions = rcp85.Emissions.emissions,
+        concentrations = rcp85.Concentrations.gases,
+        temperature = np.linspace(0,5,736),
+        feedback = -0.037,
+        emissions_pi = rcp85.Emissions.emissions[0,:],
+        concentrations_pi = rcp85.Concentrations.gases[0,:],
+    )