Gradient metrics

To isolate CO2 gradients driven by Southern Ocean fluxes, we examine CO2 anomalies relative to a local reference, using potential temperature (\(\theta\)) to delineate boundaries in the vertical. We define metrics quantifying the vertical and horizontal CO2 gradients

  • \(\Delta_{\theta}\ce{CO}_2\) is the difference between the median value of CO2 observed south of 45°S where \(\theta\) < 280 K and that in the mid- to upper-troposphere, \(\pu{295 K} < \theta < \pu{305 K}\).

  • \(\Delta_{y}\ce{CO}_2\) is the difference between CO2 averaged across stations in the core latitudes of summertime CO2-drawdown and that at the South Pole Observatory (SPO).

The atmospheric inversion models explicitly simulate CO2 tracers responsive only to ocean \((\ce{CO}_2^{ocn})\), land \((\ce{CO}_2^{lnd})\) and fossil fuel \((\ce{CO}_2^{ff})\) fluxes and subject to identical transport fields. We take advantage of these tracers to demonstrate that the gradient metrics (\(\Delta_{\theta}\ce{CO}_2\), \(\Delta_{y}\ce{CO}_2\)) are primarily responsive to local ocean influence.

Related analysis:

%load_ext autoreload
%autoreload 2
from itertools import product

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

import obs_aircraft
import emergent_constraint as ec
import datasets

import figure_panels
import util
clobber = False 
clobber_deep = False
constraint_type = 'ocean_constraint'
campaign_info = obs_aircraft.get_campaign_info(verbose=False)
model_input_lists = ec.get_model_tracer_lists(constraint_type)
model_input_lists
{'model_tracer_list': [('CT2017', 'CO2_OCN'),
  ('CT2019B', 'CO2_OCN'),
  ('CTE2018', 'CO2_OCN'),
  ('CTE2020', 'CO2_OCN'),
  ('MIROC', 'CO2_OCN'),
  ('CAMSv20r1', 'CO2_OCN'),
  ('s99oc_v2020', 'CO2_OCN'),
  ('s99oc_ADJocI40S_v2020', 'CO2_OCN'),
  ('s99oc_SOCCOM_v2020', 'CO2_OCN'),
  ('TM5-Flux-m0f', 'CO2_OCN'),
  ('TM5-Flux-mmf', 'CO2_OCN'),
  ('TM5-Flux-mrf', 'CO2_OCN'),
  ('TM5-Flux-mwf', 'CO2_OCN')],
 'model_tracer_ext_list': [('CT2017', 'CO2_LND+CO2_FFF'),
  ('CT2019B', 'CO2_LND+CO2_FFF'),
  ('CTE2018', 'CO2_LND+CO2_FFF'),
  ('CTE2020', 'CO2_LND+CO2_FFF'),
  ('CAMSv20r1', 'CO2_LND+CO2_FFF'),
  ('s99oc_v2020', 'CO2_LND+CO2_FFF'),
  ('s99oc_ADJocI40S_v2020', 'CO2_LND+CO2_FFF'),
  ('s99oc_SOCCOM_v2020', 'CO2_LND+CO2_FFF')],
 'model_list_sfco2_lnd': []}
air_parms = ec.get_parameters('default')
obj_srf = ec.whole_enchilada_srf(**model_input_lists)
obj_air = ec.whole_enchilada(clobber=clobber_deep, **model_input_lists)
ds_obs = datasets.aircraft_profiles_seasonal()
ds_obs
<xarray.Dataset>
Dimensions:              (season: 4, theta: 11)
Coordinates:
  * season               (season) <U3 'DJF' 'MAM' 'JJA' 'SON'
  * theta                (theta) float64 270.0 275.0 280.0 ... 310.0 315.0 320.0
Data variables: (12/135)
    ch4                  (season, theta) float64 -1.691 -1.112 ... -1.277 -8.285
    ch4_med              (season, theta) float64 -1.625 -0.805 ... 0.9865 -13.22
    ch4_med_std          (season, theta) float64 0.557 1.139 ... 11.36 19.76
    ch4_std              (season, theta) float64 -0.5766 -0.5704 ... 2.321 5.563
    ch4_std_std          (season, theta) float64 0.09758 0.2406 ... 2.37 4.173
    ch4mpanther          (season, theta) float64 nan 3.895 ... -1.017 -3.213
    ...                   ...
    sf6ucats_std_std     (season, theta) float64 nan 0.02498 ... 0.009239
    strat                (season, theta) float64 0.0 0.0 0.0 ... 0.2142 0.5099
    strat_med            (season, theta) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.5
    strat_med_std        (season, theta) float64 0.0 0.0 0.0 ... 0.3796 0.4841
    strat_std            (season, theta) float64 0.0 0.0 0.0 ... 0.1769 0.1614
    strat_std_std        (season, theta) float64 0.0 0.0 0.0 ... 0.2364 0.2356

Vertical gradients: ∆θCO2

ac = obj_air.get_ac(**air_parms, clobber=clobber)
ac.vg_ext
co2 gradient_mean gradient_std
model tracer campaign
TM5-Flux-mrf CO2_OCN HIPPO-1 -0.417654 -0.417654 NaN
HIPPO-2 -0.081392 -0.081392 NaN
HIPPO-3 -0.383548 -0.383548 NaN
HIPPO-5 -0.069285 -0.069285 NaN
ORCAS-J -0.560490 -0.560490 NaN
... ... ... ... ... ...
CAMSv20r1 CO2_OCN ORCAS-F -0.711700 -0.711700 NaN
ATom-1 -0.023400 -0.023400 NaN
ATom-2 -1.142900 -1.142900 NaN
ATom-3 -0.004600 -0.004600 NaN
ATom-4 -0.309100 -0.309100 NaN

450 rows × 3 columns

Horizontal gradients: ∆yCO2

data_a_clm = datasets.obs_surface_climatology('co2')
data_a_clm
<xarray.Dataset>
Dimensions:  (stncode: 8, month: 12)
Coordinates:
  * month    (month) int64 1 2 3 4 5 6 7 8 9 10 11 12
  * stncode  (stncode) object 'CRZ' 'MQA' 'DRP' 'PSA' 'SYO' 'CYA' 'MAA' 'HBA'
Data variables:
    lat      (stncode, month) float64 -46.43 -46.43 -46.43 ... -75.61 -75.61
    CO2      (stncode, month) float64 -0.4859 -0.3674 ... 0.009556 -0.1172
    CO2_std  (month) float64 0.1112 0.1557 0.09192 ... 0.0922 0.0766 0.09829
sc = obj_srf.get_sc(clobber=clobber)
sc.df_gradients_mon
season time_bounds gradient gradient_std
model tracer period month
obs CO2 2000-2004 1 DJF (2000, 2004) -0.262963 0.143441
2005-2009 1 DJF (2005, 2009) -0.191414 0.074473
2010-2014 1 DJF (2010, 2014) -0.297850 0.044947
2015-2019 1 DJF (2015, 2019) -0.346297 0.043477
1999-2020 1 DJF (1999, 2020) -0.254041 0.111184
... ... ... ... ... ... ... ...
CAMSv20r1 CO2_OCN 2010-2014 12 DJF (2010, 2014) -0.162536 0.029921
2015-2019 12 DJF (2015, 2019) -0.174630 0.027346
1999-2020 12 DJF (1999, 2020) -0.171655 0.055785
1999-2009 12 DJF (1999, 2009) -0.110739 0.055632
2009-2020 12 DJF (2009, 2020) -0.182731 0.048142

3864 rows × 4 columns

Seasonal aircraft profiles

def plot_seasonal_profiles(ds_obs, ax):
    theta = ds_obs.theta
    seasons = ds_obs.season.values
    field = 'co2'
    for n, season in enumerate(seasons):
        dsi = ds_obs.sel(season=season)
        x, xerr = dsi[f'{field}_med'], dsi[f'{field}_med_std']
        ax.errorbar(x, theta, 
                    xerr=xerr, 
                    marker='.', 
                    label=season,
                    lw=2, 
                    c=figure_panels.palette_colors[n],
                    markersize=10, 
                    zorder=100,
                   )
    ax.axvline(0., lw=0.5, c='k')    
    ax.legend(frameon=False)    
    ax.set_ylabel('$\\theta$ [K]')
    ax.set_xlabel('$\Delta$CO$_2$ [ppm]')
        
    xlm = ax.get_xlim()

fig, axs = util.canvas(1, 1, figsize=(4, 6))
plot_seasonal_profiles(ds_obs, axs[0, 0])

tbin = ac.theta_bins[0]
axs[0, 0].set_title(f'Aircraft obs: CO$_2$ minus ({tbin[0]:0.0f}-{tbin[1]:0.0f}K)');

util.savefig(f'seasonal-profiles-obs.pdf')
_images/gradients-main_13_0.png

Observed gradients

fig, axs = util.canvas(1, 1)
figure_panels.obs_theta_gradient(ac.flight_gradients, axs[0, 0], ac.theta_bins);

util.savefig(f'seasonal-DthetaCO2-obs.pdf')
_images/gradients-main_15_0.png

Simulated gradients

def vg_seasonal_cycle_model_flavors(ac, ax, co2_components=['CO2', 'CO2_OCN', 'CO2_LND', 'CO2_FFF',]):

    df = ac.vg_ext
    # 'CO2_LND+CO2_FFF']
   
    doy = []
    for c in campaign_info.keys():
        doy.append(util.day_of_year(np.array([campaign_info[c]['time']]))[0])

    data = {k: [] for k in co2_components}
    error = {k: [] for k in co2_components}
    for c in campaign_info.keys():    
        for tracer in co2_components:
            data[tracer].append(
                df.xs((c, tracer), level=('campaign', 'tracer')).gradient_mean.median()
            )
            error[tracer].append(
                df.xs((c, tracer), level=('campaign', 'tracer')).gradient_mean.std(ddof=1)
            )

    xhat = np.linspace(-30, 365+30, 365)

    for tracer in co2_components:
        x, y = util.antyear_daily(np.array(doy), np.array(data[tracer]))
        _, yerr = util.antyear_daily(np.array(doy), np.array(error[tracer]))
        ax.errorbar(x, y, yerr=yerr, 
                    linestyle='none', 
                    color=figure_panels.co2_colors[tracer], 
                    marker='.', 
                    markersize=10,
                    label=figure_panels.co2_names[tracer], 
                   )
        p, pcov = ec.curve_fit(ec.harmonic, xdata=x/365.25, ydata=y, sigma=yerr, absolute_sigma=True,)
        yhat = ec.harmonic(xhat/365.25, *p)
        ax.plot(xhat, yhat, color=figure_panels.co2_colors[tracer], lw=2)
        if tracer == 'CO2':
            yhat_out = yhat
            
    ax.axhline(0., lw=0.5, c='k')
    ax.legend(loc=(1.02, 0), ncol=1, fontsize=8);  
    ax.set_xlim((-10, 375))
    ax.set_xticks(figure_panels.bomday)
    ax.set_xticklabels([f'        {m}' for m in figure_panels.monlabs_ant]+['']);
    ax.legend(frameon=False)
    ax.set_ylabel('$\Delta_{ \\theta}$CO$_2$ [ppm]')
    return yhat_out

fig, axs = util.canvas(1, 1)
co2_air = vg_seasonal_cycle_model_flavors(ac, axs[0, 0]) 

theta_str = figure_panels.theta_bin_def(ac.theta_bins)
axs[0, 0].set_title(f'Models: {theta_str} CO$_2$ diff')
util.savefig(f'seasonal-DthetaCO2-model.pdf')
_images/gradients-main_17_0.png
fig, axs = util.canvas(1, 1)
co2_air = vg_seasonal_cycle_model_flavors(ac, axs[0, 0], co2_components=['CO2']) 

theta_str = figure_panels.theta_bin_def(ac.theta_bins)
axs[0, 0].set_title(f'Models: {theta_str} CO$_2$ diff')
util.savefig(f'seasonal-DthetaCO2-model-total.pdf')
_images/gradients-main_18_0.png
fig, axs = util.canvas(1, 1)
co2_air = vg_seasonal_cycle_model_flavors(ac, axs[0, 0], co2_components=['CO2', 'CO2_OCN']) 

theta_str = figure_panels.theta_bin_def(ac.theta_bins)
axs[0, 0].set_title(f'Models: {theta_str} CO$_2$ diff')
util.savefig(f'seasonal-DthetaCO2-model-total-ocn.pdf')
_images/gradients-main_19_0.png
fig, axs = util.canvas(1, 1)
co2_air = vg_seasonal_cycle_model_flavors(ac, axs[0, 0], co2_components=['CO2', 'CO2_LND', 'CO2_FFF']) 

theta_str = figure_panels.theta_bin_def(ac.theta_bins)
axs[0, 0].set_title(f'Models: {theta_str} CO$_2$ diff')
util.savefig(f'seasonal-DthetaCO2-model-total-lnd-fff.pdf')
_images/gradients-main_20_0.png

Observed gradients

fig, axs = util.canvas(1, 1)
figure_panels.obs_srf_seasonal(axs[0, 0], data_a_clm)
util.savefig(f'seasonal-DyCO2-obs.pdf')
_images/gradients-main_22_0.png
def hg_seasonal_cycle_model_flavors(sc, ax, period='1999-2020', 
                                    co2_components=['CO2', 'CO2_OCN', 'CO2_LND', 'CO2_FFF',]):

    df = sc.df_gradients_mon
    # 'CO2_LND+CO2_FFF']
   
    doy = util.doy_midmonth()

    data = {k: [] for k in co2_components}
    error = {k: [] for k in co2_components}
    for month in range(1, 13):
        for tracer in co2_components:
            gradient = df.xs((period, month, tracer), level=('period', 'month', 'tracer')).gradient
            data[tracer].append(gradient.median())
            error[tracer].append(gradient.std(ddof=1))
    
    xhat = np.linspace(-30, 365+30, 365)

    for tracer in co2_components:
        x, y = util.antyear_daily(np.array(doy), np.array(data[tracer]))
        _, yerr = util.antyear_daily(np.array(doy), np.array(error[tracer]))
        ax.errorbar(x, y, yerr=yerr, 
                    linestyle='none', 
                    color=figure_panels.co2_colors[tracer], 
                    marker='.', 
                    markersize=10,
                    label=figure_panels.co2_names[tracer], 
                   )
        p, pcov = ec.curve_fit(ec.harmonic, xdata=x/365.25, ydata=y, sigma=yerr, absolute_sigma=True,)
        yhat = ec.harmonic(xhat/365.25, *p)            
        ax.plot(xhat, yhat, color=figure_panels.co2_colors[tracer], lw=2)
        if tracer == 'CO2':
            yhat_out = yhat

    ax.axhline(0., lw=0.5, c='k')
    ax.legend(loc=(1.02, 0), ncol=1, fontsize=8);  
    ax.set_xlim((-10, 375))
    ax.set_xticks(figure_panels.bomday)
    ax.set_xticklabels([f'        {m}' for m in figure_panels.monlabs_ant]+['']);
    ax.legend(frameon=False, loc="lower left")
    ax.set_ylabel('$\Delta_{ y}$CO$_2$ [ppm]')
    return yhat_out 

fig, axs = util.canvas(1, 1)
co2_srf = hg_seasonal_cycle_model_flavors(sc, axs[0, 0]) 
axs[0, 0].set_title(f'Models: SO CO$_2$ minus SPO')
util.savefig(f'seasonal-DyCO2-model.pdf')
_images/gradients-main_23_0.png
fig, axs = util.canvas(1, 1)
co2_srf = hg_seasonal_cycle_model_flavors(sc, axs[0, 0], co2_components=['CO2']) 
axs[0, 0].set_title(f'Models: SO CO$_2$ minus SPO')
util.savefig(f'seasonal-DyCO2-model-total.pdf')
_images/gradients-main_24_0.png
fig, axs = util.canvas(1, 1)
co2_srf = hg_seasonal_cycle_model_flavors(sc, axs[0, 0], co2_components=['CO2', 'CO2_LND', 'CO2_FFF']) 
axs[0, 0].set_title(f'Models: SO CO$_2$ minus SPO')
util.savefig(f'seasonal-DyCO2-model-total-lnd-fff.pdf')
_images/gradients-main_25_0.png
fig, axs = util.canvas(1, 1)
co2_srf = hg_seasonal_cycle_model_flavors(sc, axs[0, 0], co2_components=['CO2', 'CO2_OCN']) 
axs[0, 0].set_title(f'Models: SO CO$_2$ minus SPO')
util.savefig(f'seasonal-DyCO2-model-total-ocn.pdf')
_images/gradients-main_26_0.png

How much bigger is the vertical gradient?

(co2_air.max() - co2_air.min())/(co2_srf.max() - co2_srf.min())
3.451083991601859

Fig 2: Seasonal evolution atmospheric CO2 gradients

fig = plt.figure(figsize=(14, 6)) #dpi=300)

# set up plot grid
gs = gridspec.GridSpec(
    nrows=2, ncols=2, 
    left=0.38, right=1,
    bottom=0., top=1,
    hspace=0.25, wspace=0.24,
)

nrow = 2
ncol = 2
axs = np.empty((nrow, ncol)).astype(object)
for i, j in product(range(nrow), range(ncol)):
    axs[i, j] = plt.subplot(gs[i, j])

ax_p = fig.add_axes([0.075, 0.08, 0.23, 0.82])
plot_seasonal_profiles(ds_obs, ax_p)
tbin = ac.theta_bins[0]
ax_p.set_title(f'Aircraft obs: CO$_2$ minus ({tbin[0]:0.0f}-{tbin[1]:0.0f}K)')
figure_panels.obs_theta_gradient(ac.flight_gradients, axs[0, 0], ac.theta_bins);
ylm = axs[0, 0].get_ylim()

theta_str = figure_panels.theta_bin_def(ac.theta_bins)
vg_seasonal_cycle_model_flavors(ac, axs[1, 0])
axs[1, 0].set_title(f'Models: {theta_str} CO$_2$ diff')

figure_panels.obs_srf_seasonal(axs[0, 1], data_a_clm)

hg_seasonal_cycle_model_flavors(sc, axs[1, 1])
axs[1, 1].set_title(f'Models: SO CO$_2$ minus SPO')

util.label_plots(fig, [ax_p]+[ax for ax in axs.ravel()], xoff=-0.02)

util.savefig(f'metrics-seasonal-cycle.pdf')
_images/gradients-main_30_0.png

Seasonal evolution of atmospheric CO2 over the Southern Ocean. (A) Vertical profiles of CO2 observations collected by aircraft south of 45°S, binned on 5 K potential temperature (θ) bins and averaged by season (whiskers show standard deviation). (B) The vertical gradient (\(\Delta_{\theta}\ce{CO}_2\)) in CO2 measured from aircraft south of 45°S. Small points show \(\Delta_{\theta}\ce{CO}_2\) for individual profiles; larger points show the median and standard deviation (whiskers) for each flight. The black line shows a two-harmonic fit to the flight-median points. (C) Monthly climatology (1999–2019) of the latitudinal gradient in CO2 measured by surface stations (Fig 1); the black line shows the station mean metric (\(\Delta_{y}\ce{CO}_2\)). Separate laboratory records at SYO and PSA have been averaged. The seasonal evolution of (D) \(\Delta_{\theta}\ce{CO}_2\) and (E) \(\Delta_{y}\ce{CO}_2\) simulated in a collection of atmospheric inverse models. The points show the median across the models, whiskers show the standard deviation; the colors correspond to the “total” CO2 (black), and CO2 tracers responsive to only ocean (blue), land (green), and fossil (red) surface fluxes. Note that the y-axis bounds differ by panel.