Add CO2 flux decompisition to carbon submodule

[bradyrx]

See https://github.com/bradyrx/climpred/issues/5

[aaronspring]

I tried to do the decomposition. Its based on the HAMOCC output. Maybe we need to add a model flag to adapt it for different models.

I start with pCO2 decomp. of positive ENSO state: Freshwater contribution not implemented yet.
0. Load data ds_terms contains 'spco2','dissicos','talkos','tos','sos' from control run

ds_terms_anomaly = ds_terms-ds_terms.mean('time')

1. sensitivity:

def spco2_sensitivity(ds):
    DIC = ds['dissicos']
    ALK = ds['talkos']
    SALT = ds['sos']
    TEMP = ds['tos']
    pCO2 = ds['spco2']
    buffer_factor = dict()
    buffer_factor['ALK'] = -ALK**2 / ((2 * DIC - ALK) * (ALK - DIC))
    buffer_factor['DIC'] = (3*ALK*DIC - 2*DIC**2) / \
                           ((2 * DIC - ALK) * (ALK - DIC))
    # Compute sensitivities
    sensitivity = dict()
    sensitivity['tos'] = 0.0423
    sensitivity['sos'] = 1/SALT
    sensitivity['talkos'] = (1/ALK)*buffer_factor['ALK']
    sensitivity['dissicos'] = (1/DIC)*buffer_factor['DIC']
    sensitivity = xr.Dataset(sensitivity) * pCO2
    return sensitivity

sensitivity = spco2_sensitivity(ds_terms)

2. regression against enso(time)

enso = (ds_terms['tos']*nino34_mask).mean(['x', 'y']).squeeze()
def regression_against_index(ds, index, psig=None):
    terms = dict()
    for term in ds.data_vars:
        if term is 'spco2':
            continue
        print('Progress ...', term)
        # rewrote xr_linregress that you can input x,y
        reg = xr_linregress(index, ds[term], psig=psig)
        terms[term] = reg['slope']  # * sensitivity[term]
    terms = xr.Dataset(terms)
    return terms

terms = regression_against_index(ds_terms, enso)

3. Test whether decomposition seems useful: add up contributions

# shouldnt I also just take the sensitivities from ENSO-pos states?
(terms*sensitivity.mean('time')*1).to_array().sum('variable').plot(cmap='RdBu_r',
                                                                   vmin=-6, vmax=6, yincrease=False)

4. Compare with my ENSO composite.

def composite_analysis(field,
                       timeseries,
                       threshold=1,
                       plot=True,
                       ttest=True,
                       psig=0.05,
                       **plot_kwargs):
    index = standardize(timeseries)
    field = field - field.mean('time')
    composite = create_composites(field, index, threshold=threshold)
    if ttest:
        # test if pos different from none
        index = 'pos'
        m1 = composite.mean('time').sel(index=index)
        s1 = composite.std('time').sel(index=index)
        n1 = len(composite.groups[index])
        index = 'none'
        m2 = composite.mean('time').sel(index=index)
        s2 = composite.std('time').sel(index=index)
        n2 = len(composite.groups[index])
        # implemented ttest_ind_from_stats
        t, p = xr_ttest_ind_from_stats(m1, s1, n1, m2, s2, n2)
        comp_pos = composite.mean('time').sel(index='pos').where(p < psig)

        # test if neg different from none
        index = 'neg'
        m1 = composite.mean('time').sel(index=index)
        s1 = composite.std('time').sel(index=index)
        n1 = len(composite.groups[index])
        index = 'none'
        m2 = composite.mean('time').sel(index=index)
        s2 = composite.std('time').sel(index=index)
        n2 = len(composite.groups[index])

        t, p = xr_ttest_ind_from_stats(m1, s1, n1, m2, s2, n2)
        comp_neg = composite.mean('time').sel(index='neg').where(p < psig)

        composite = xr.concat([comp_pos, comp_neg], dim='index')
        composite['index'] = ['positive', 'negative']
    else:
        composite = composite.mean('time')
    if plot:
        composite.sel(index='pos').plot(**plot_kwargs)
        plt.show()
        composite.sel(index='neg').plot(**plot_kwargs)
    else:
        return composite

enso_comp = composite_analysis(
    pCO2, enso, plot=False, threshold=1, ttest=False).sel(index='pos')
enso_comp.plot(robust=True, vmin=-6, vmax=6, cmap='RdBu_r', yincrease=False)

5. Seems like a good fit. When 3 and 4 look similar, I am confident in the decomposition and the size of the terms.

(terms*sensitivity.mean('time')*1).to_array().plot(col='variable',
                                                   yincrease=False, robust=True, vmin=-6, vmax=6, cmap='RdBu_r')

As you can see in 5. I was a bit unsure about the scaling. (there 1).

Where I am really unsure how the whole method works for decomposing a timeseries and not just the ENSO-state. And then how to scale it.

[aaronspring]

Decomposing timeseries: (Somehow I though this will also involve xr_linregress but it doesnt)
0. Load data ds_terms contains 'spco2','dissicos','talkos','tos','sos' from control run

ds_terms_anomaly = ds_terms-ds_terms.mean('time')

1. scale anomalies with sensitivity

vmax = 40
norm = ds_terms_anomaly   #anomaly
#terms_in_pCO2_units = terms_spco2*sensitivity*norm
terms_in_pCO2_units = sensitivity*norm

2. Compare reconstructed pCO2 based on anomaly*sensitivity with the real pCO2 output.

x = terms_in_pCO2_units.to_array().sum('variable')
x.name = 'reconstr_spco2'
d = xr.merge([x, pCO2-pCO2.mean('time')])
d.isel(time=slice(0, 4)).to_array().plot(
    yincrease=False, robust=True, col='time', row='variable')

diff = (d['reconstr_spco2']-d['spco2'])
sum = (d['reconstr_spco2']+d['spco2'])
diff.isel(time=slice(0, 4)).plot(col='time', yincrease=False, robust=True)

Here there is large differences in the high-lats. Tropics seem fine.

terms_in_pCO2_units.isel(time=slice(0, 4)).to_array().plot(col='time', row='variable',
                                                           yincrease=False, robust=True, vmin=-vmax, vmax=vmax, cmap='RdBu_r')

[aaronspring]

And for CO2flux decomposition, there are less assumptions for the sensitivity derivation because we just take the partial derivatives.

So my questions are:

• Does the method work like this?
• What is your experience in analysing the mismatch between the reconstructed pCO2 and the real pCO2? Is there this kind of check?

I will add the composite method in that issue/PR also.