Predictive asymmetry analysis for Late PleistoceneEarth system dynamics.

(1)

time grid label for this notebook: NBR_postMPT_500

Predictive asymmetry analysis for Late Pleistocene Earth system dynamics.

Written by Maria Salem, in supplement to master thesis.

In this notebook we compute the normalized predictive asymmetry ( ) between empirical climate proxy records covering the last 480 kyrs leading up to the Holocene. The aim is to deduce causal coupling strength and directionality in the climate system between ice volume, northern hemisphere summer

insolation, atmospheric pCO2 concentration, and Southern Ocean Fe fertilization after the Mid-Pleistocene Transition.



Preamble: We import all necessary packages, wrangled time series and prededﬁned functions for analysis, by including notebook 3 MA_NB3_Toolbox.ipynb

(2)

ﬁle:///Users/maria/Jottacloud/Notebooks/Notebooks_PDF/NBR_postMPT_500.html 2/176

In [18]:

using NBInclude

@nbinclude "../NBRs_ePalus_ns20/NB3_Toolbox_ns20.ipynb"

(3)

The LR04 d18O record spans from -5319.0 to -1.0 kyrs BP.

The La2004 insolation time series (Ins) spans from -5000.0 to -0.0 k yrs BP.

Spratt Lisiecki GSL stack (SpraSL) spans from -797.0 to -1.0 kyrs B P.

The Elderfield GSL record (EldSL) spans from -1574.0 to -8.0 kyrs B P.

The Grant sea level record (GraSL) spans from -491.0 to -1.0 kyrs B P.

The Rohling sea level record (RohSL) spans from -5329.0 to -1.0 kyrs BP.

The Bereiter pCO2 record (BerCO2) spans from -803.0 to -2.0 kyrs BP.

The Chalk pCO2 record (ChaCO2) spans from -1240.0 to -1092.0 kyrs B P.

The Lambert dust record (IceDust) spans from -799.0 to -13.0 kyrs B P.

The Martinez-Garcia Fe flux record (MarFe) spans from -3999.0 to -2.

0 kyrs BP.

The higher resolution part of the MG record spans from -800.0 to -0.

5 kyrs BP.The record is now on the common time gridThe record is now on the common time gridResults are saved in the .jld2 file.

(4)

25.792495 seconds (133.00 M allocations: 9.590 GiB, 18.19% gc time)

(5)

The following records are used for analysis over this notebook's time interval (record labels in bold / code labels in italics):

Ins / La2004: Numerical solution for northern hemisphere insolation (top of atmosphere solar ﬂux mean for summer solstice at ), by Laskar et al. (2004).

LR04: global reference stack - a principal component analysis of 57 records, by Lisiecki &

Raymo (2005).

SpraSL / SL: Global sea level stack (principal component analysis) of 7 sea level records, whereof 4 spanning the last 800 kyrs, by Spratt & Lisiecki (2016).

EldSL / E: Temperature deconvolution of signal in marine sediment core from the Chatnam Rise (Paciﬁc Ocean), spanning the last 1.5 Myrs, by Elderﬁeld et al. (2012).

GraSL / G: High resolution record of the relative sea level at Bab el Mandab (Red Sea) spanning the last 500 kyrs, by Grant et al. (2014).

RohSL / R: Relative sea level at Straight of Gibraltar, spanning the last 5.3 Myrs (Rohling et al., 2014).

BerCO2 / B: Record of atmospheric pCO2 spanning the last 800 kyrs, from Epica Dome C ice core, East Antarctica, by Bereiter et al. (2015).

ChaCO2 / C: High resolution proxy record for pCO2, spanning a 150 kyr interval of the MPT, by Chalk et al. (2017).

IceDust / L: Record of Antarctic dust accumulation from Epica Dome C ice core, spanning the past 800 kyrs, by Lambert et al. (2008).

MarFe / MG: Marine sediment record of Fe accumulation in the Southern Ocean spanning the last 4 Ma, by Martinez-García et al. (2011).

𝑁 65

^∘

𝑂

𝛿

¹⁸

𝛿

¹⁸

𝑂

𝑂 𝛿

¹⁸

𝐵 𝛿

¹¹

Outline for this notebook:

1. Set the time series to a common time grid

Decide on a common time interval for analyses in this notebook.

Cut all records to the relevant time interval.

1. Run pairwise analyses for predictive asymmetry between the time series

Compute the normalized predictive asymmetry (function deﬁned in NB3) for each time series pair, and save the results in a .jld2 ﬁle.

Produce results plot of the normalized predictive asymmetry for each time series pair.

1. Produce overview plots of predictive asymmetry results

Produce the overview plot of the ensemble of predictive asymmetry results for the time period here studied. The ensemble plot will be use in the results chapter of main text.

1. Set the time series to a common time grid

For our method to work, it is important we have the data on the same time grid - covering the exact same time interval and with a regular time step. The time step was set by binsize of the grid in the

BinnedResampling in NB1, where we chose a regular grid with one value for every 1000 years. The common time interval, on the other hand, we set individually in each notebook.

(6)

- Decide on common time grid

Determine the common time interval for analyses in this notebook. We decide on a time interval for analysis according to where the records overlap. Additionally, we can delimit the time interval to a period of interest - in our case, we want to see if there are any changes in causal dynamics before, during and after the Mid-Pleistocene Transition.

In [3]:

# Determine the common time interval

#### We select a time interval according to the periods covered by the records

# the record with the "latest start" defines the start of the common grid gridstart = maximum([

tmin_LR04, tmin_La2004, tmin_SL, tmin_E,

tmin_G, # the Grant record is the shortest record and delimits the time interval analysed.

tmin_R, tmin_B,

# tmin_C, # the Chalk record is left out as it does not cover the post-M PT time interval

tmin_L, tmin_MG ])

# the record with the "earliest end" the end of the common grid gridend = minimum([

tmax_LR04, tmax_La2004, tmax_SL, tmax_E, tmax_G, tmax_R, tmax_B, # tmax_C, tmax_L, tmax_MG ])

print("Time interval for analyses in this notebook is from ", -gridstart, " ka B P to ", -gridend, " ka BP")

# Note: opposite to NB1, tmin and tmax here indicate the binmidpoints of the com mon grid.

Time interval for analyses in this notebook is from 491.0 ka BP to 1 3.0 ka BP

(7)

In [4]:

gridstart = maximum([tmin_B, tmin_G]) gridend = minimum([tmax_B, tmax_G])

Recall the time step for the common grid. This was given by the binsize in the grids for

BinnedResampling in NB1. We also recall time steps for the additional analyses of the high resolution records.

In [5]:

# Recall the binsize for the common grid

binsize_1 = 1 # 1 kyrs - is the default timestep on which all time seri es are binned (NB1)

# Additionally, we had prepared some records for higher resolution analyses binsize_hr125 = 0.125 # High resolution (hr) records are additionally binned o n a hr grid, 125 year timestep (Grant, Chalk, La2004)

binsize_hr500 = 0.5 # Martinez-García and La2004 are additionally binned on a g rid with a 500 year timestep

# (but the main timestep for all analyses will be 1000 years)

print("All records are on a regular grid with timestep of ", binsize_1, " kyr.

Additional high resolution analyses have timesteps of ", binsize_hr125 , " and "

,binsize_hr500, " kyrs.")

We can now deﬁne the time grid for analyses in this notebook. Objects associated with this grid are labeled with the grid suﬃx postMPT_500 .

Out[4]:

-2.0

All records are on a regular grid with timestep of 1 kyr.

Additional high resolution analyses have timesteps of 0.125 and 0.5 kyrs.

(8)

In [41]:

# the common grids for the time series are then defined by commongrid = gridstart : binsize_1 : gridend

Binmidpoints_commongrid =[commongrid[i] for i in 1:length(commongrid)]

# Additional high resolution grids

commongrid_hr125 = gridstart : binsize_hr125 : gridend

Binmidpoints_commongrid_hr125 = [commongrid_hr125[i] for i in 1:length(commongri d_hr125)]

commongrid_hr500 = gridstart : binsize_hr500 : gridend

Binmidpoints_commongrid_hr500 = [commongrid_hr500[i] for i in 1:length(commongri d_hr500)]

### But this is the main grid

print(gridstart : binsize_1 : gridend , " defines the main common grid for analy ses in this notebook.")

#= *Note*:

the common grid for the time series in this notebook is defined by tmin and tmax as the *bin midpoints*

(opposite to NB1, where tmin and tmax defined the *grid edges* for binned resamp ling.) =#

Binmidpoints_commongrid_hr500;

Binning of state space (ϵ) for transfer entropy estimation in the following analyses

The transfer entropy is a probability distribution, which we estimate by the visitation frequency test. Our estimation of the transfer entropy is thus sensitive to the binning resolution of state space (the amount of the bins we use to count visitations), which is given by . We choose according to the length of our time series, as given by the Palus horizon. This is all deﬁned inside our function_from_XY_to_normPA, see NB3, and NB2 for explanation for further information. For overview, let's check the used for the time series in this notebook.

𝜖 𝜖

𝜖

-491.0:1.0:-2.0 defines the main common grid for analyses in this no tebook.

(9)

In [42]:

# Palus horizon

N = length(Binmidpoints_commongrid) eD = 3

ϵ = Int(round( N^(1/(eD+1)) ))

print("binning resolution ϵ used to estimate transfer entropy is ", ϵ, " for the 1 kyr grid (time series length N = ", N,")."

)

# Palus horizon

N = length(Binmidpoints_commongrid_hr500) eD = 3

ϵ = Int(round( N^(1/(eD+1)) )) print("

", ϵ, " for the 500 yr grid (time series length N = ", N,")."

)

N = length(Binmidpoints_commongrid_hr125) eD = 3

ϵ = Int(round( N^(1/(eD+1)) )) print("

and ", ϵ, " for the 125 yr grid, (time series length N = ", N,")."

)

2. Cut the time series to the decided time interval.

We select the time series' common time array as following: all time values greater (younger) than, or equal to, the common grid's starting midpoint tmin, and all values smaller (older) than, or equal to, the common grid's ending midpoint tmax.

Insolation - La2004

La2004 is a numerical solution for insolation, with no associated uncertainty for the time interval here observed. Therefore no conﬁdence. interval on this record

Default version (La2004_insol_cut) is on a grid with 1 kyr timestep between each insolation value (as in the original dataset).

binning resolution ϵ used to estimate transfer entropy is 5 for the 1 kyr grid (time series length N = 490).

6 for the 500 yr grid (time series length N = 979).

and 8 for the 125 yr grid, (time series length N = 3913).

(10)

In [52]:

# cut out the relevant time interval (tmin:tmax) from the time array

La2004_t_cut = La2004_t_fullength[(La2004_t_fullength .>= gridstart) .& (La2004_

t_fullength .<= gridend)]

# cut out the relevant time interval from the insolation values array

La2004_insol_cut = La2004_insol65N_fullength[(La2004_t_fullength .>= gridstart) .& (La2004_t_fullength .<= gridend)]

# check that we have cut the correct time interval

print(La2004_t_cut[1] : La2004_t_cut[end] , " - the cut version of La2004 time s eries

", gridstart : gridend , " - is now the same as the common grid") # so all good

# define the time series plot plot_La2004 =

plot(La2004_t_cut, La2004_insol_cut, color = :orange,

label = "Ins",

xlabel = "Time [kyrs BP]", ylabel = string("Solar flux

", L"[W/m^{2}]"),

legend = :topleft, grid = false,

size = (800,200) );

High resolution (La2004_insol_cut_hr125)

Prepare a second version on the same grid as the Chalk pCO2 record / Grant GSL record, for high resolution analysis between the two (interpolated insolation values for every 125 years).

-491.0:1.0:-2.0 - the cut version of La2004 time series -491.0:1.0:-2.0 - is now the same as the common grid

(11)

In [7]:

# La2004 hr version with time step 125 years on the common grid

La2004_t_cut_hr125 = La2004_t_fullength_hr125[(La2004_t_fullength_hr125 .>= grid start) .& (La2004_t_fullength_hr125 .<= gridend)]

La2004_insol_cut_hr125 = La2004_insol65N_fullength_hr125[(La2004_t_fullength_hr1 25 .>= gridstart) .& (La2004_t_fullength_hr125 .<= gridend)]

# Check that the time series cut at the correct time interval (= tmin : tmax), a nd that it has the right timestep (= hr125)

timestep = La2004_t_cut_hr125[2] - La2004_t_cut_hr125[1]

record_timegrid = La2004_t_cut_hr125[1] : timestep : La2004_t_cut_hr125[end]

begin

if record_timegrid == commongrid_hr125

print("The record is now on the common time grid") else

print("Something's wrong") end

end

# define the time series plot plot_La2004_hr125 =

plot(La2004_t_cut_hr125, La2004_insol_cut_hr125, color = :orange,

label = "Ins_hr125",

xlabel = "Time [kyrs BP]",

ylabel = L"Solar \ flux \ [W/m^{2}]", legend = :topleft,

grid = false, size = (800,200) );

High resolution (La2004_insol_cut_hr500)

Prepare a third version on the same grid as the Martinez-García hr record, for high resolution analysis between the two (interpolated insolation values for every 500 years).

The record is now on the common time grid

(12)

In [8]:

# La2004 hr version with time step 125 years on the common grid

La2004_t_cut_hr500 = La2004_t_fullength_hr500[(La2004_t_fullength_hr500 .>= grid start) .& (La2004_t_fullength_hr500 .<= gridend)]

La2004_insol_cut_hr500 = La2004_insol65N_fullength_hr500[(La2004_t_fullength_hr5 00 .>= gridstart) .& (La2004_t_fullength_hr500 .<= gridend)]

# Check that the time series cut at the correct time interval (= tmin : tmax), a nd that it has the right timestep (= hr125)

timestep = La2004_t_cut_hr500[2] - La2004_t_cut_hr500[1]

record_time_grid = La2004_t_cut_hr500[1] : timestep : La2004_t_cut_hr500[end]

begin

if record_time_grid == commongrid_hr500

print("The record is now on the common time grid") else

end

# define the time series plot plot_La2004_hr500 =

plot(La2004_t_cut_hr500, La2004_insol_cut_hr500, color = :orange,

label = "Ins_hr500", xlabel = "Time [kyrs BP]",

ylabel = L"Solar \ flux \ [W/m^{2}]", legend = :topleft,

grid = false, size = (800,200) );

uivD cut troubleshooting

In [11]:

#= After several attempts, we see that the >= (greater than OR EQUAL TO) operato r does not work for uivDs (our time series format, which we want to cut).

We will therefore set the grid edges tmin and tmax a littlebit before and after the start and end of the common grid, so that we don't lose the datapoints on t he edges =#

tmin = gridstart - 0.001 tmax = gridend + 0.001

;

#= select the time series' common time array as following:

all values greater (younger) than ``tmin``, and smaller (older) than ``tmax`` =#

(13)

In [10]:

LR04 = LR04_binned_fullength_fullageunc

# Cut out the relevant time interval of the binned LR04 time series LR04_cut = UncertainIndexValueDataset(

LR04.indices[(LR04.indices .> tmin) .& (LR04.indices .< tmax)], # pick out a ll indices (ages) where ages are larger than tmin and smaller than tmax

LR04.values[(LR04.indices .> tmin) .& (LR04.indices .< tmax)] # pick out a ll values where the corresponding ages are larger than tmin and smaller than tma x

) Out[10]:

UncertainIndexValueDataset{UncertainIndexDataset,UncertainValueDatas et} containing 479 uncertain values coupled with 479 uncertain indic es

(14)

In [11]:

# check that the record was cut correctly and is now on the common time grid record_timestep = LR04_cut.indices[2].value - LR04_cut.indices[1].value

record_timegrid = LR04_cut.indices[1].value : record_timestep : LR04_cut.indices [end].value

if record_timegrid == Binmidpoints_commongrid print("The record is now on the common time grid") else

## or alternatively

binmidpoints_ts =[LR04_cut.indices[i].value for i in 1:length(LR04_cut.indices)]

if binmidpoints_ts == Binmidpoints_commongrid print("The record is now on the common time grid") else

#### Plot time series with the 95% confidence interval

# computing the median in each bin (0.5 quantile), and the confidence interval w e want to use (95%)

bin_median = quantile.(LR04_cut.values, 0.5)

bin_upperq = quantile.(LR04_cut.values, 0.975) .- bin_median bin_lowerq = bin_median .- quantile.(LR04_cut.values, 0.025)

;

plot_LR04 =

plot(binmidpoints_ts, bin_median, ribbon = (bin_lowerq, bin_upperq), fillalpha = 0.3,

color = :black, label = "LR04",

xlabel = "Time [kyrs BP]", ylabel = string(L"\delta{18}O","

",L"[\perthousand]"), grid = false, legend = :topleft, yflip = true

)

(15)

The record is now on the common time gridThe record is now on the co mmon time grid

Out[11]:

(16)

In [12]:

LR04 = LR04_binned_fullength_noageunc

# Cut out the relevant time interval of the binned LR04 time series LR04_cut_noageunc = UncertainIndexValueDataset(

LR04.indices[(LR04.indices .> tmin) .& (LR04.indices .< tmax)], # pick out a ll indices (ages) where ages are larger than tmin and smaller than tmax

LR04.values[(LR04.indices .> tmin) .& (LR04.indices .< tmax)]# pick out all values where the corresponding ages are larger than tmin and smaller than tmax )

# check that the record was cut correctly and is now on the common time grid binmidpoints_ts = [LR04_cut_noageunc.indices[i].value for i in 1:length(LR04_cut _noageunc.indices)]

begin

end

#### Plot time series with the 95% confidence interval

bin_median = quantile.(LR04_cut_noageunc.values, 0.5)

bin_upperq = quantile.(LR04_cut_noageunc.values, 0.975) .- bin_median bin_lowerq = bin_median .- quantile.(LR04_cut_noageunc.values, 0.025) plot_LR04_noageunc =

plot(

Binmidpoints_commongrid, bin_median, ribbon = (bin_lowerq, bin_upperq), fillalpha = 0.3, legend = :topleft, color = :black,

label = "LR04 - without

age model uncertainty",

ylabel = L"\delta{18}O \ [\perthousand]", grid = false,

yflip = true )

(17)

Sea level - SprattLisiecki GSL stack

The record is now on the common time grid Out[12]:

(18)

In [51]:

ts = SL_binned_fullength_noageunc # age uncertainty already included

# cut time series

ts_cut = UncertainIndexValueDataset(

ts.indices[(ts.indices .> tmin) .& (ts.indices .< tmax)], ts.values[(ts.indices .> tmin) .& (ts.indices .< tmax)])

# check that the record was cut correctly and is now on the common time grid binmidpoints_ts =[ts_cut.indices[i].value for i in 1:length(ts_cut.indices)]

begin

end

# save cut time series in an unambiguous name SL_cut = ts_cut

### Plot time series with the 95% confidence interval

bin_median = quantile.(ts_cut.values, 0.5)

bin_upper = quantile.(ts_cut.values, 0.975) .- bin_median bin_lower = bin_median .- quantile.(ts_cut.values, 0.025) plot_SL =

plot(binmidpoints_ts, bin_median, ribbon = (bin_lower, bin_upper), fillalpha = 0.3, legend = :topleft, color = :darkblue,

label = "SpraSL",

xlabel = "Time [kyrs BP]", ylabel = "GSL

[m]",

grid = false )

(19)

Sea level - Elderﬁeld record

(20)

In [14]:

ts = E_binned_fullength_ageunc

# cut time series

begin

end

# save cut time series in an unambiguous name E_cut = ts_cut

bin_upper = quantile.(ts_cut.values, 0.975) .- bin_median bin_lower = bin_median .- quantile.(ts_cut.values, 0.025)

;

binmidpoints_commongrid =[ts_cut.indices[i].value for i in 1:length(ts_cut.indic es)]

plot_E =

plot(binmidpoints_ts, bin_median, ribbon = (bin_lower, bin_upper),

fillalpha = 0.3, legend = :topleft, color = :royalblue, label = "EldSL",

xlabel = "Time [kyrs BP]", ylabel = "GSL

[m]",

grid = false )

(21)

Sea level - Grant record 1 kyr grid G_cut

(22)

In [19]:

ts = G_binned_fullength

# cut time series

begin

end

# save in an unambiguous name G_cut = ts_cut

;

plot_G =

fillalpha = 0.3, legend = :topleft, color = :cyan, label = "GraSL",

xlabel = "Time [kyrs BP]", ylabel = string(L"RSL_{BeM}", "

[m]"),

grid = false )

(23)

High resolution version with 500 yr timestep (G_cut_hr500), for hr analysis with Martinez-García The record is now on the common time grid

Out[19]:

(24)

In [30]:

ts = G_binned_fullength_hr500

# cut time series

begin

if binmidpoints_ts == Binmidpoints_commongrid_hr500 print("The record is now on the common time grid") else

end

# save in an unambiguous name G_cut_hr500 = ts_cut

;

plot_G_hr500 =

fillalpha = 0.3, legend = :topleft, color = :cyan, label = "GraSL_hr500",

ylabel = string(L"RSL_{BeM}", " [m]"), grid = false)

;

Mean resolution of Lambert record before interpolation was 155 years. We therefore consider it reasonable to run higher resolution analyses with this record, both with the Grant GSL record (time step 125 years), and with the Bereiter pCO2 and Martinez-García dust records (timestep 500 years)

High resolution version with 125 yr timestep (G_cut_hr125) The record is now on the common time grid

(25)

In [28]:

ts = G_binned_fullength_hr125

# cut time series

begin

end

# save in an unambiguous name G_cut_hr125 = ts_cut

;

plot_G_hr125 =

fillalpha = 0.3, legend = :topleft, color = :cyan, label = "GraSL_hr125",

ylabel = string(L"RSL_{BeM}", " [m]"), grid = false)

;

Sea level - Rohling record

(26)

In [18]:

ts = R_binned_full

# cut time series

begin

end

# save in an unambiguous name R_cut = ts_cut

;

plot_R =

fillalpha = 0.3, legend = :topleft, color = :skyblue, label = "RohSL",

xlabel = "Time [kyrs BP]", ylabel = string(L"RSL_{Gib}", "

[m]"),

grid = false )

(27)

pCO2 - Bereiter record 1 kyr grid (B_cut)

(28)

In [20]:

ts = B_binned_fullength_ageunc

# cut time series

begin

end

# save in an unambiguous name B_cut = ts_cut

# Plot time series with the 95% confidence interval

bin_upper = quantile.(ts_cut.values, 0.975) .- bin_median bin_lower = bin_median .- quantile.(ts_cut.values, 0.025) plot_B =

fillalpha = 0.3, legend = :topleft, color = :red, label = "BerCO2",

xlabel = "Time [kyrs BP]", ylabel = string(L"pCO_{2}", "

[ppm]"),

grid = false )

(29)

B_cut_hr500

(30)

In [21]:

ts = B_binned_fullength_ageunc_hr500

# cut time series

begin

end

# save in an unambiguous name B_cut_hr500 = ts_cut

bin_upper = quantile.(ts_cut.values, 0.975) .- bin_median bin_lower = bin_median .- quantile.(ts_cut.values, 0.025) plot_B_hr500 =

fillalpha = 0.3, legend = :topleft, color = :red, label = "BerCO2_hr500",

ylabel = string(L"pCO_{2}", " [ppm]"), grid = false

)

;

B_cut_hr125

(31)

In [27]:

ts = B_binned_fullength_ageunc_hr125

# cut time series

begin

end

# save in an unambiguous name B_cut_hr125 = ts_cut

bin_upper = quantile.(ts_cut.values, 0.975) .- bin_median bin_lower = bin_median .- quantile.(ts_cut.values, 0.025) plot_B_hr125 =

fillalpha = 0.3, legend = :topleft, color = :red, label = "BerCO2_hr125",

ylabel = string(L"pCO_{2}", " [ppm]"), grid = false

)

;

We also make a version without age uncertainty, to use for analysis with the Lambert dust record (both are from Epica Dome C ice core, and can thus be analysed by depth instead of age model)

B_cut_edc

(32)

In [22]:

# version without age uncertainty

ts = B_binned_fullength_noageuncEDC

# cut time series

begin

end

# Save with an unambiguous name B_cut_edc = ts_cut

bin_upper = quantile.(ts_cut.values, 0.975) .- bin_median bin_lower = bin_median .- quantile.(ts_cut.values, 0.025) plot_B_edc =

fillalpha = 0.3, legend = :topleft, color = :red, label = "BerCO2 -

without age model uncertainty",

xlabel = "Time [kyrs BP]", ylabel = string(L"pCO2","

[ppm]"),

grid = false )

(33)

B_cut_edc_hr500

(34)

In [23]:

# version without age uncertainty

ts = B_binned_fullength_noageuncEDC_hr500

# cut time series

begin

end

# Save with an unambiguous name B_cut_edc_hr500 = ts_cut

bin_upper = quantile.(ts_cut.values, 0.975) .- bin_median bin_lower = bin_median .- quantile.(ts_cut.values, 0.025) plot_B_edc_hr500 =

fillalpha = 0.3, legend = :topleft, color = :red,

label = "BerCO2_hr500 - without age model uncertainty", xlabel = "Time [kyrs BP]",

ylabel = string(L"pCO2", " [ppm]"), grid = false

)

;

Dust - Lambert record 1 kyr grid

(35)

In [24]:

ts = L_binned_full

# cut time series

begin

end

# Save the cut time series in an unambiguous name L_cut = ts_cut

;

plot_L =

fillalpha = 0.3, legend = :topleft, color = :lime, label = "IceDust",

xlabel = "Time [kyrs BP]", ylabel = "Dust conc.

[µg/kg]",

grid = false )

(36)

We also prepare a high resolution version, for analysis with the Grant record L_cut_hr125

(37)

In [25]:

ts = L_binned_full_hr125

# cut time series

begin

end

# Save the cut time series in an unambiguous name L_cut_hr125 = ts_cut

;

plot_L_hr125 =

fillalpha = 0.3, legend = :topleft, color = :lime, label = "IceDust_hr125",

xlabel = "Time [kyrs BP]", ylabel = "Dust conc. [µg/kg]", grid = false

)

;

L_cut_hr500

(38)

In [26]:

ts = L_binned_full_hr500

# cut time series

begin

end

# Save the cut time series in an unambiguous name L_cut_hr500 = ts_cut

;

plot_L_hr500 =

fillalpha = 0.3, legend = :topleft, color = :lime, label = "IceDust_hr500",

xlabel = "Time [kyrs BP]", ylabel = "Dust conc. [µg/kg]", grid = false

)

;

(39)

We also prepare a version without age uncertainty, for analysis with the Bereiter pCO2 record (Both records are ice cores from the same place, Epica Dome C. This allows us to analyse by depth instead of the

constructed age model, and thus avoiding the age uncertainty.

Help: some age uncertainty in the determination of lock-in time of gas in ice SHOULD be taken into account. However, I have not found how to do this (don't understand explanation in dataset (lidie?)

L_cut_edc

(40)

In [27]:

ts = L_binned_full_EDC

# cut time series

begin

end

# Save the cut time series in an unambiguous name L_cut_edc = ts_cut

bin_upperq = quantile.(ts_cut.values, 0.975) .- bin_median bin_lowerq = bin_median .- quantile.(ts_cut.values, 0.025)

;

plot_L_edc =

plot(binmidpoints_ts, bin_median, ribbon = (bin_lowerq, bin_upperq),

fillalpha = 0.3, legend = :topleft, color = :lime, label = "IceDust -

age uncertainty between gas/ice",

xlabel = "Time [kyrs BP]", ylabel = "Dust conc.

[µg/kg]",

grid = false )

(41)

L_cut_edc_hr500

(42)

In [28]:

ts = L_binned_full_hr500_edc

# cut time series

begin

end

# Save the cut time series in an unambiguous name L_cut_edc_hr500 = ts_cut

bin_upperq = quantile.(ts_cut.values, 0.975) .- bin_median bin_lowerq = bin_median .- quantile.(ts_cut.values, 0.025)

;

plot_L_edc_hr500 =

plot(binmidpoints_ts, bin_median, ribbon = (bin_lowerq, bin_upperq),

fillalpha = 0.3, legend = :topleft, color = :lime,

label = "IceDust_hr500 - age uncertainty between gas and ice", xlabel = "Time [kyrs BP]",

ylabel = "Dust conc. [µg/kg]", grid = false

);

dust - Martinez-García Fe MAR record

(43)

In [29]:

ts = MG_binned_fullength

# cut time series

begin

end

# save the cut time series uivD in an unambiguous name MG_cut = ts_cut

;

plot_MG =

fillalpha = 0.3, legend = :topleft, color = :green, label = "MarFe",

xlabel = "Time [kyrs BP]", ylabel = string("Fe MAR","

",L"[g/m^{2}/yr]"), grid = false )

(44)

We also prepare a high resolution version with timestep of 500 years, for analysis with La2004, Chalk and Grant hr records:

MG_cut_hr500

(45)

In [30]:

ts = MG_binned_postmpt_hr500

# cut time series

begin

end

# save the cut time series uivD in an unambiguous name MG_cut_hr500 = ts_cut

;

plot_MG_hr500 =

fillalpha = 0.3, legend = :topleft, color = :green, label = "MarFe_hr500",

ylabel = L"Fe \ MAR \ [g/m^{2}/year]", grid = false

)

;

Overview of time series

(46)

In [31]:

p1 = plot(plot_LR04, xlabel = "", xaxis = :off) p2 = plot(plot_SL, xlabel = "", xaxis = :off) p3 = plot(plot_E, xlabel = "", xaxis = :off) p4 = plot(plot_R, xlabel = "", xaxis = :off) p5 = plot(plot_G, xlabel = "", xaxis = :off) p6 = plot(plot_B, xlabel = "", xaxis = :off)

p6_edc = plot(plot_B_edc, xlabel = "", xaxis = :off) p7 = plot(plot_L, xlabel = "", xaxis = :off)

p7_edc = plot(plot_L_edc, xlabel = "", xaxis = :off) p8 = plot(plot_MG, xlabel = "", xaxis = :off)

p9 = plot(plot_La2004, xlabel = "", xaxis = :off)

p_xaxis = plot(xlabel = "Time [kyrs BP]", xaxis = :on, yaxis = :off, grid = :off )

A4 = (8.27*72, 11.69*72) # Dimensions of an A4 page (72 dots per inch) l = @layout [a;b;c;d;e;f;g;h;i;j;k;l{0.01h}]

po_alltimeseries = plot(

p1,p2,p3,p4,p5,p6,p6_edc,p7,p7_edc,p8,p9,p_xaxis,

#plot_LR04, plot_SL, plot_E, plot_G, plot_B, plot_B_edc, plot_L, plot_L_edc, plot_MG, plot_La2004,

layout = l, #grid(12,1), size = (500,1200), #size = A4,

xlims = (gridstart-50, gridend+250), xticks =(-500:100:0), #ylabel = "", ytickfont = false, yaxis = :off,

legend = :right, bg_legend = :white, ymirror = false

)

savefig("../../results_ePalus_ns20/timeseries/po_all_postMPT_500.pdf")

(47)

(48)

In [32]:

#plot(po_alltimeseries, size = A4, yaxis = :off, ytickfont = :off);

Compute the predictive asymmetry between the time series

See method explained in NB2 and function deﬁned in NB3

Compute the predictive asymmetry between the two time series

Normalize the results, and make a plot showing the 95% conﬁdence interval

(49)

Outline for analyses

1. - insolation LR04 - La2004 1. GSL - insolation

Spratt & Lisiecki - La2004 Elderﬁeld - La2004 Rohling - La2004 Grant - La2004

1. GSL / - pCO2 (post-MPT) Spratt & Lisiecki - Bereiter Elderﬁeld - Bereiter Rohling - Bereiter (..) Grant - Bereiter

1. GSL / - pCO2 (syn-MPT) Spratt & Lisiecki - Chalk Elderﬁeld - Chalk Rohling - Chalk Grant - Chalk 1. GSL / - dust

Spratt & Lisiecki - Lambert (..) Elderﬁeld - Lambert (..) Rohling - Lambert (...) Grant - Lambert (..)

Spratt & Lisiecki - Martinez-García Elderﬁeld - Martinez-García Rohling - Martinez-García Grant - Martinez-García 1. pCO2 - dust

Bereiter - Lambert (!)

Bereiter - Martinez-García (..) Chalk - Martinez-García 1. pCO2 - insolation

Bereiter - La2004 Chalk - La2004 1. dust - insolation

Lambert - La2004

Martinez-García - La2004

𝑂 𝛿

¹⁸

𝑂 𝛿

¹⁸

𝑂 𝛿

¹⁸

𝑂

𝛿

¹⁸

(50)

1. between and insolation

LR04 - La2004

 𝛿

¹⁸

𝑂

Let's start with LR04 and insolation - because the LR04 d18O stack is tuned to the insolation signal, and this

"bias" is carried on by tuning other gsl records to the LR04 stack, so it is interesting to see what to keep in mind from this.

(Both the Elderﬁeld and the Spratt-Lisiecki sea-level records are tuned to the LR04, whose age model is based on a lag of behind insolation. This is based on the understanding of the the Milankovitch cycles as a key driver behind the ice age cycles, but is however, a caveat for our method, which seeks to infer this causal connection, rather than a priori assume it. We will therefore compute the predictive asymmetry between the LR04 record and the La2004 northern hemisphere summer insolation record, to ... check for bias?

What would we expect? Drive insolation -> response . But we would expect this either way if it was caused by real signal or age model assumptions... Not sure what to make out as an argument here..)

𝑂 𝛿

¹⁸

𝑂 𝛿

¹⁸

𝑂

𝛿

¹⁸

(51)

In [33]:

# Plot overview of time series pair to analyse plot_overview_LR04_La2004 =

plot(layout = grid(2,1), size = (500, 400), plot_LR04,

plot_La2004 )

Out[33]:

(52)

Compute the predictive asymmetry between the time series pair. We compute the normalized predictive asymmetry using the function deﬁned in NB3.

In [ ]:

# Recall the time series on the common grid as X and Y X = LR04_cut

Y = La2004_insol_cut

# Compute the predictive asymmetry (function defined in NB3)

@time computePredictiveAsymmetries(X, Y, timestep = 1, ηmax = 20, ϵ = 5,

filepath = "../../results_ePalus_ns20/pa_jld2_files/postMPT_500/LR04_La2004n B.jld2")

(53)

In [34]:

# normalize the results for comparability

@load "../../results_ePalus_ns20/pa_jld2_files/postMPT_500/LR04_La2004nB.jld2"

normPA_XtoY = normalizePredictiveAsymmetry(rsTE_XtoY, rsPA_XtoY, ηmax = ηmax, f

= 1)

normPA_YtoX = normalizePredictiveAsymmetry(rsTE_YtoX, rsPA_YtoX, ηmax = ηmax, f

= 1);

#Calculate the quantiles for the confidence interval we want to plot (2σ, aka 9 5%).

### Calculate the quantiles for the confidence interval we want to plot (2σ, aka 95%)

# ... for PA from X to Y

normPA_XtoY_median = [quantile(normPA_XtoY[i,:], 0.5) for i in 1:ηmax]

# median

normPA_XtoY_upper = [quantile(normPA_XtoY[i,:], 0.975) for i in 1:ηmax] .- normP A_XtoY_median # upper quantile

normPA_XtoY_lower = normPA_XtoY_median .- [quantile(normPA_XtoY[i,:], 0.025) for i in 1:ηmax] # lower quantile

# ... for PA from Y to X

normPA_YtoX_median = [quantile(normPA_YtoX[i,:], 0.5) for i in 1:ηmax]

# median

normPA_YtoX_upper = [quantile(normPA_YtoX[i,:], 0.975) for i in 1:ηmax] .- normP A_YtoX_median # upper quantile

normPA_YtoX_lower = normPA_YtoX_median .- [quantile(normPA_YtoX[i,:], 0.025) for i in 1:ηmax] # lower quantile

;

#Define the results plotLR04_La2004.jld2

# define the predictive asymmetry plot plot_normPA_LR04_La2004 =

plot(#title = string("Normalized # define the predictive asymmetry plot", L"\mat hcal{A}", ") between ", L"\delta^{18}O", " (LR04) and northern hemisphere summer insolation (La2004)"),

xlims = (0, ηmax),

xticks = (0 : 5 : ηmax), ylims = (-4,2),

yticks = (-5:1:5), legend = :bottomleft,

xlabel = string(L"ηs"," [kyr]"), ylabel = L"\mathcal{A}",

size = (400,400), border = true, grid = true,

hline([1], line = (:dash, :black)), label = "" # significance treshold )

plot!(normPA_XtoY_median, # ...from X to Y

ribbon = (normPA_XtoY_lower, normPA_XtoY_upper), label = string("LR04", L"\rightarrow", "Ins"), fillalpha = 0.3, color = :black

)

plot!(normPA_YtoX_median, # ...from Y to X

ribbon = (normPA_YtoX_lower, normPA_YtoX_upper), label = string("Ins", L"\rightarrow", "LR04"),

(54)

)

# Plot overview of time series plot_overview_LR04_La2004 = plot(

layout= grid(2,1), plot_LR04,

plot_La2004)

# join plots of time series and pa results to a results subplot plot_results_LR04_La2004 =

plot(size = (500,250), layout = grid(1,2),

plot_overview_LR04_La2004, plot_normPA_LR04_La2004)

savefig("../../results_ePalus_ns20/pa_ResultPlots/postMPT_500/LR04_La2004nB.pdf"

)

Caption: Normalized predictive asymmetry between and northern hemisphere summer insolation.

Insolation (orange) is the La2004 numerical solution of insolation intensity on June at (Laskar et al., 2004). , shown in black, is the LR04 global stack by Lisiecki & Raymo (2005). Note that the age model of this record is built on the premiss that the \d18O signal follows northern hemisphere insolation intensity with an assumed lag of 300 years. This poses a caveat in the results of our analysis.

𝛿

¹⁸

𝑂

𝑏

31

^𝑠𝑡

65

^∘

𝑁 𝛿

¹⁸

𝑂

𝑏

2. between sea level and insolation 

SprattLisiecki - La2004

(55)

In [53]:

# First, recall the time series on the common grid as X and Y X = SL_cut

filepath = "../../results_ePalus_ns20/pa_jld2_files/postMPT_500/SL_La2004nB.

jld2")

Results are saved in the .jld2 file. 38.016789 seconds (215.90 M al locations: 15.629 GiB, 18.74% gc time)

(56)

In [54]:

# Load and normalize the predictive asymmetry results

# Load and normalize the random sequences results of predictive asymmetry

@load "../../results_ePalus_ns20/pa_jld2_files/postMPT_500/SL_La2004nB.jld2"

= 1)

= 1);

# calculate the quantiles for the 95% confidence interval # ... for PA from X to Y

# median

;

# defining the results plot plot_normPA_SL_La2004 = plot(#title =

xlims = (0, ηmax),

xticks = (0 : 5 : ηmax),

ylims = (-4,2), yticks = (-5:1:5), legend = :bottomleft,

xlabel = string(L"ηs"," [kyr]"), # prediction lags

ylabel = L"\mathcal{A}", # normalized predictive asymmetry...

plot!(normPA_XtoY_median, # ...from X to Y

ribbon = (normPA_XtoY_lower, normPA_XtoY_upper), label = string("SpraSL", L"\rightarrow", "Ins"), fillalpha = 0.3, color = :darkblue

)

plot!(normPA_YtoX_median, # ...from Y to X

ribbon = (normPA_YtoX_lower, normPA_YtoX_upper), label = string("Ins", L"\rightarrow", "SpraSL"), fillalpha = 0.3, color = :orange)

;

# join plots of time series and pa results to a results subplot plot_overview_SL_La2004 =

plot(layout = grid(2,1), size = (500,400), plot_SL,

plot_La2004

(57)

)

plot_results_SL_La2004 = plot(size = (500,250),

layout = grid(1,2), plot_overview_SL_La2004, plot_normPA_SL_La2004)

savefig("../../results_ePalus_ns20/pa_ResultPlots/postMPT_500/SL_La2004nB.pdf")

Elderﬁeld - La2004

In [ ]:

# Compute the predictive asymmetry between the time series

# First, recall the time series on the common grid as X and Y X = E_cut

@time computePredictiveAsymmetries(

X, Y, timestep = 1, ηmax = 20, ϵ = 5,

filepath = "../../results_ePalus_ns20/pa_jld2_files/postMPT_500/E_La2004nB.j ld2"

)

(58)

In [36]:

# Load and normalize the predictive asymmetry results

@load "../../results_ePalus_ns20/pa_jld2_files/postMPT_500/E_La2004nB.jld2"

= 1)

= 1);

# calculate the quantiles to plot the 95% confidence interval # ... for PA from X to Y

# median

;

# defining the results plot plot_normPA_E_La2004 = plot(xlims = (0, ηmax),

xticks = (0 : 5 : ηmax), ylims = (-4,3),

yticks = (-5:1:5), legend = :bottomleft,

xlabel = string(L"ηs"," [kyr]"), ylabel = L"\mathcal{A}",

plot!(normPA_XtoY_median,

ribbon = (normPA_XtoY_lower, normPA_XtoY_upper), label = string("EldSL", L"\rightarrow", "Ins"), fillalpha = 0.3, color = :royalblue

)

plot!(normPA_YtoX_median,

ribbon = (normPA_YtoX_lower, normPA_YtoX_upper), label = string("Ins", L"\rightarrow", "EldSL"), fillalpha = 0.3, color = :orange

)

;

# join plots of time series and pa results to a results subplot plot_overview_E_La2004 =

plot(layout = grid(2,1), size = (500,400), plot_E,

(59)

plot_La2004 )

plot_results_E_La2004 = plot(size = (500,250),

layout = grid(1,2), plot_overview_E_La2004, plot_normPA_E_La2004)

savefig("../../results_ePalus_ns20/pa_ResultPlots/postMPT_500/E_La2004nB.pdf")

Rohling - La2004

In [ ]:

# Recall the time series on the common grid as X and Y X = R_cut

filepath = "../../results_ePalus_ns20/pa_jld2_files/postMPT_500/R_La2004nB.j ld2")