Preamble

This notebook deals with the comparison between inventories and fluxes of nitrate on one hand, and different aspects of primary production on the other, such as net primary production, new production, and export fluxes.

All figures exported from this notebook are prefixes by FIGURE_PP_.

In [ ]:
%load_ext autoreload
%autoreload 2
%run imports.py

dims = dict(
    FN=hv.Dimension('FN', label='Nitrate flux', unit='mmol N m⁻² d⁻¹', range=(.005, 10)),
    pp_no3_eq=hv.Dimension('pp_no3_eq', label='NPP', unit='mmol N m⁻² d⁻¹'),
    flux=hv.Dimension('flux', label='N flux', unit='mmol N m⁻² d⁻¹', range=(.01, 10)),
)

Net primary production

First, we visualize annual net primary production inferred using ocean colour [@arrigo2015continued] as an average summer value during the last two decades. This map can then be compared with observed upward fluxes of nitrate during winter. It will not come as a surprise that NPP several times larger than the upward nitrate fluxes as it includes a considerable fraction regenerated production.

FIGURE_PP_map_NPP-FN

Annual average maps

Data were acquired as binary files through personal contact with G. van Dijken and K. Arrigo and converted to the netcdf format.

In [ ]:
if os.path.exists('../data/primprod/'):
    ds = xr.open_mfdataset(
        '../data/primprod/*summer_05-09_mean.nc',
        concat_dim='date', combine='nested').compute()
    ds = ds.mean(dim='date')
    df = ds.to_dataframe()
    df['pp_no3_eq'] = df.pp / 12 * 16 / 106
    # df.pp /= 1e3

    def logmean(x):
        return np.nanmedian(np.log10(x))

    options = (
        opts.HexTiles(aggregator=logmean, projection=ccrs.NorthPolarStereo(), gridsize=30,
                      tools=['hover'], colorbar_opts=dict(major_label_text_font_size='12pt'),
                      colorbar=True, width=570, height=510, cmap=default_cmap,
                      hooks=[logcolor_ticks([.4, 1., 3, 10, 30, 60])]),
        opts.Feature(scale='50m'),
    )

    pp = gv.HexTiles(df, ['lon', 'lat'], 'pp_no3_eq')
    l = pp * land
    l = l.redim(**dims).opts(*options)
    hv.save(l, '../nb_fig/FIGURE_PP_map', fmt='html')
else:
    print('Data not available')

Pan-Arctic patterns of annual net primary production show a very productive Atlantic sector, somewhat lower values in the Beaufort Sea, and a very oligotrophic central basin, similar to winter fluxes of nitrate. The main difference difference is that net primary production is at least an order of magnitude larger than any of the fluxes, so it seems hard to get more out of this than qualitative agreement at best.

The reason is of course that "net primary production" includes regenerated production, meaning the same nitrogen is just recycled over and over again, and can mostly do without input of new nitrogen.

In [ ]:
display(HTML('../nb_fig/FIGURE_PP_map.html'))

Pan-Arctic comparison of winter nitrate fluxes with annual new and export production

First, we need to compile pan-Arctic estimates of annual new

New production

Sakshaug's 2004 article in "The Organic Carbon Cycle in the Arctic Ocean" by Stein & MacDonald (eds.) has a good compilation. We use the following values:

In [ ]:
gC_per_yr_to_mmolN_per_day = 1. /12 /106*16 * 1e3/365

data = {
    # Area: new production estimate in [gC / yr]
    'Barents Sea': [8., 100],
    'Chukchi Sea': [5, 160],
    'Beaufort Sea': [7, 17],
    'Baffin Bay': [25, 50],
    'Central basin': [.5]
}

df = pd.DataFrame(zip(*[(k, v) for k, vs in data.items() for v in vs]), index=['Area', 'newp']).transpose()
df.newp = gC_per_yr_to_mmolN_per_day * df.newp.astype(float)

df = df.groupby('Area').mean().reset_index()

df_newp = df.copy()

Export flux of particulate organic carbon

The basis of our brief discussion of POC export fluxes are the values compiled by I. Wiedmann's PhD thesis (2015, Fig. 2), supplemented with a number of additional references to increase data coverage in selected regions.

Ingrid Wiedmann, Fig. 2 PhD thesis

These values are all valid for a depth of 200 m.

Because the raw data to produce the graph was not available at the time of writing, we did some screen grabbing using a custom python script. This script is saved as grab_data.py, alongside the original and the reproduced figure.

In [ ]:
# extracted raw coordinates relative to the png
xy = pd.read_csv('../data/wiedmann2015_export_fluxes/coords.csv', index_col=[0])[['x', 'y']]

# x0, ...: figure coordinates
# x0d, ...: data coordinates. y: logscale!
x0, x1, y0, y1 = xy.loc[0, 'x'], xy.loc[1, 'x'], xy.loc[0, 'y'], xy.loc[2, 'y']
x0d, x1d, y0d, y1d = 1., 12., 0., 3.

def convertx(c):
    return round((c-x0)/(x1-x0) * (x1d-x0d) + x0d)

def converty(c):
    return (c-y0)/(y1-y0) * (y1d-y0d) + y0d

df = dd.read_csv('../data/wiedmann2015_export_fluxes/wiedmann2015_fig2/*.csv').compute()

df['month'] = df.x.apply(convertx)
df['exp'] = 10**df.y.apply(converty)

df['exp_n_eq'] = df.exp /12 /106 *16
df = df.drop(columns=['x', 'y'])

# Sorting data

df['Season'] = 'Winter'
df.loc[df.month.between(4, 9), 'Season'] = 'Summer'

df['Area'] = ''

df.loc[df.ref.str.contains('barents'), 'Area'] = 'Barents Sea'
df.loc[df.ref.str.contains('nsvalbard'), 'Area'] = 'Barents Sea'
df.loc[df.ref.str.contains('beaufort-amundsen'), 'Area'] = 'Amundsen Gulf'
df.loc[df.ref.str.contains('baffin|nwp'), 'Area'] = 'Baffin Bay'
df.loc[df.ref.str.contains('kara'), 'Area'] = 'Kara'
df.loc[df.ref.str.contains('nwp'), 'Area'] = 'Northwater'
df.loc[df.ref.str.contains('greenland'), 'Area'] = 'Greenland Sea'
In [ ]:
l = hv.Dataset(df, ['month', 'exp', 'ref']).to(hv.Points, groupby='ref').overlay().opts(
    opts.Points(legend_position='right', width=700, height=500, size=10, logy=True, padding=.1, jitter=.1)
)
hv.save(l.opts(toolbar=None), '../data/wiedmann2015_export_fluxes/Wiedmann-Figure2-recreated.png')

Some additional export estimates

These all pertain to the Canadian Basin.

In [ ]:
#Honjo et al. 2010
flux = [0.05, 0.65] # gC m-2 yr-1

#Cai et al. 2010
flux.append(0.9) # gC m-2 yr-1

# convert to mmol N m-2 d-1
flux =  np.array(flux) / 12/106*16 *1e3/365

df = df.append(pd.DataFrame(dict(exp_n_eq=flux, Season='Summer', Area='Canadian Basin')), sort=False)

df_exp = df.copy()

TBD

miquel2015downward:

0.2–2.5 mg N m−2 d−1

In [ ]:
print(0.2/14)
print('--')
print(2.5/14)
In [ ]:
500/12/106*16

Plot

In [ ]:
options = opts.BoxWhisker(
    xrotation=45, xlabel='', logy=True, padding=.1, toolbar=None, yticks=[0.1, 1, 5, 10],
    height=500, width=500,
)

exp = hv.BoxWhisker(df_exp, ['Season', 'Area'], hv.Dimension('exp_n_eq', label='C export, Redfield N-equiv.', unit='mg N m-2 d-1'))
exp.opts(options)

Nitrate fluxes

In [ ]:
df_fn = pd.read_csv('../data/fn-compilation.csv')

fn = hv.BoxWhisker(df_fn, ['Season', 'Area'], hv.Dimension('FN', label='FN', unit='mg N m-2 d-1'))
fn.opts(options)

Merge compilations of export production, new production, and nitrate fluxes

In [ ]:
df_fn = df_fn.replace({'Perennial':'Winter', 'Barents Sea, AABC': 'Barents Sea',
                       'N Svalbard/Fram Strait': 'Barents Sea', 'Canada Basin': 'Central basin', 
                       'Makarov Basin': 'Central basin', 'Amundsen Gulf': 'Beaufort Sea'})

df_fn = df_fn.groupby(['Season', 'Area']).mean()['FN'].loc['Winter'].reset_index()

df_exp = df_exp.replace({'Amundsen Gulf': 'Beaufort Sea', 'Canadian Basin': 'Central basin'})

df_exp = df_exp.groupby(['Season', 'Area']).mean()['exp_n_eq'].loc['Summer'].reset_index()

df = df_fn.merge(df_exp, how='outer').merge(df_newp, how='outer').groupby('Area').mean()

df = df.stack().rename_axis(['Area', 'kind']).rename('flux').to_frame().reset_index()

df = df.replace({'exp_n_eq': 'Vertical export', 'FN': 'NO₃ flux', 'newp': 'New production'})

df.head()
# Enter uncertainties manually from each paper # Bourgault et al, 2011 df_exp.loc[('Amundsen Gulf', 'NO3 flux'), 'flux_std'] = 0.04 # unpublished df_exp.loc[('Baffin Bay', 'NO3 flux'), 'flux_std'] = 0.5 # Randelhoff et al., 2015 df_exp.loc[('Barents Sea', 'NO3 flux'), 'flux_std'] = 0.5
In [ ]:
hv.renderer('bokeh').theme = theme

opts_bk = [opts.Bars(
    color=hv.Cycle(['#AA5F77', '#2DA3A2', '#D6AE4A']),
    backend='bokeh',
    xrotation=45, xlabel='', logy=True, yticks=[.01, 0.1, 1, 10],
    width=800, height=400, 
    fontsize={v: 14 for v in ['ylabel', 'xlabel', 'yticks', 'xticks']},
    line_width=3, yformatter='%g', show_legend=True,
    tools=['hover'],
    toolbar=None,
)]
opts_mpl = translate_options(opts_bk, bokeh2mpl)

l = hv.Bars(df, ['Area', 'kind'], 'flux').opts(*opts_bk, *opts_mpl).redim(**dims)
l = l[['Baffin Bay', 'Barents Sea', 'Beaufort Sea', 'Central basin']]
fname = '../nb_fig/FIGURE_PP_FN_NEWP_EXP_regional'
hv.save(l, fname, fmt='png')
hv.save(l, fname, fmt='html')
save_bokeh_svg(l, fname+'.svg')
l

Single case studies

First, manually enter the data from the studies.

In [ ]:
## Randelhoff et al 2016 JGR

df = pd.DataFrame(dict(
    FN=[1.2, 0.6, 0.3, 1.1, 0.3, 0.1], 
    newP=[2.6, 3.1, 8.4, 0.015, 0.018, 0.048],
    #station=['P1', 'P3', 'P4', 'P5', 'P6', 'P7'], 
    season=['spring', 'spring', 'spring', 'summer', 'summer', 'summer']
))

df = df.set_index('season').stack().rename('flux')

df = df.rename_axis(['season', 'kind']).reset_index()

df = df.replace(dict(spring='Spring', summer='Summer', newP='New prod.', FN='NO3 flux')).assign(ref='aRandelhoff et al. (2016)')

## Nishino et al 2018

df = df.append(
    pd.DataFrame({'New prod.': 0.24, 'NO3 flux': .19,}, index=[0]).stack()
    .rename('flux').rename_axis(('','kind')).reset_index(1)
    .assign(season='Summer', ref='bNishino et al. (2018)'),
    sort=False, ignore_index=True
)

# Also append our Young Sound data

def YS_newP():
    """
    Load Young Sound new production estimates
    """
    df = pd.read_csv('../data/youngsound2015/YS2015_N15_inc.csv')
    df = df.rename(columns={
        'newProd_uM_N_d-1': 'newP', 'regProd_uM_N_d-1': 'regP', 
        'depth_m': 'depth', 'NOx_init_uM': 'NOx',
    })
    df['fratio'] = df.newP / (df.newP + df.regP)
    
    # It is difficult to estimate integrated new production 
    # from measurements at only two depths. Here we give two very rough estimates:
    # One is conservative:
    # by assuming the 5 m value holds in the [0-5 m] interval and the 20 m value holds over [5, 20].
    newPint_lo = (
        df.set_index(['depth', 'station']).newP.to_xarray()
        * xr.DataArray(data=[5, 15], coords=dict(depth=[5, 20]), dims=['depth'])
    ).sum(dim='depth').rename('newPint').to_dataframe().assign(estimate='lo')
    # the other one is presumably much larger than the real value: We assume the 5-m value holds down to 20m, 
    # and the 20-m value holds down to 40 m
    newPint_hi = (
        df.set_index(['depth', 'station']).newP.to_xarray()
        * xr.DataArray(data=[20, 20], coords=dict(depth=[5, 20]), dims=['depth'])
    ).sum(dim='depth').rename('newPint').to_dataframe().assign(estimate='hi')
    return pd.concat([newPint_hi, newPint_lo])

def YS_FN():
    """
    Load Young Sound nitrate fluxes
    """
    df = xr.load_dataset('../data/youngsound2015/YS2015_FN.nc').to_dataframe()
    on_sill = df.station.isin(['YS3.02', 'GH01', 'Tyro07']) | df.station.str.contains('SILL')
    df.loc[on_sill, 'station_type'] = 'sills'
    df.loc[~on_sill, 'station_type'] = 'interior'
    return df.assign(kind='NO₃ flux, '+df.station_type)[['FN', 'kind']]

df = df.append(
    pd.concat([
        YS_newP().newPint.rename('flux').to_frame()
        .assign(kind='New prod.'), 
        YS_FN().rename(columns={'FN': 'flux'}),
    ], sort=False).assign(ref='cYoung Sound', season='Summer'),
    sort=False
)
In [ ]:
options = (
    opts.BoxWhisker(logy=True, xlabel='', xrotation=45, yformatter='%g',
                    frame_width=600, frame_height=200, box_fill_color='grey',
                    padding=.1, tools=['hover']),
)

l = hv.BoxWhisker(
    df.replace('NO3 flux', 'NO₃ flux'), 
    [hv.Dimension('ref', value_format=lambda s: s[1:]), 'season', 'kind'], 
    hv.Dimension('flux', label='N flux', range=(.01, 10), unit=u'mmol N m\u207B\u00B2 d\u207B\u00B9')
)
l = l.opts(*options)
panelA = l.clone().opts(title='A')

Seasonal cycle of DON and PON (Paulsen et al., 2019)

First, enter the values given in their Table 1.

In [ ]:
df = pd.DataFrame(dict(
    Month=[1, 3, 5, 8, 11],
    DON=   [84., 51, 35, 109, 81],
    DONerr=[9, 32, 22, 36, 10],
    PON=   [4., 8, 46, 30, 3],
    PONerr=[1, 5, 41, 17, 1]
))

Then define an NdOverlay of Spread elements:

In [ ]:
griddict = {pool: hv.Spread(df, 'Month', [pool, pool+'err'])
            for pool in ['PON', 'DON']}

options = [
    opts.Curve(padding=.1, line_width=3, width=500, tools=['hover'], align='end',
               xticks=[(1, 'Jan'), (3, 'Mar'), (5, 'May'), (8, 'Aug'), (11, 'Nov')]),
    opts.NdOverlay(legend_position='top_left'),
]

err = hv.HoloMap(griddict)
conc = err.map(hv.Curve, hv.Spread)

l = err*conc
l = l.overlay()
l.opts(*options)
for c, g in zip(['k', 'g'], ['PON', 'DON']):
    l[g].opts(opts.Curve(color=c), opts.Spread(color=c))

l = l.redim(DON=hv.Dimension('DON', label='DON, PON', unit='µM'))
panelB = l.clone().opts(title='B')
In [ ]:
l = (panelA + panelB).cols(1)
In [ ]:
fname = '../nb_fig/FIGURE_PP_incubations_vs_flux'
hv.save(l, fname, fmt='html')
hv.output(l)
l = l.opts(toolbar=None)
hv.save(l, fname, fmt='png')
save_bokeh_svg_multipanel(l, fname+'.svg', 'v', align='start')