MATLAB Examples

INFAUNAL: INdividual Foraminiferal Approach UNcertainty AnaLysis

Contents

1. About

INFAUNAL: INdividual Foraminiferal Approach UNcertainty AnaLysis

This code constructs synthetic time series with varying ENSO and seasonal cycle amplitudes for a point location in the tropical Pacific Ocean and produces an IFA probability contour plot.

Written by Kaustubh Thirumalai (kau@ig.utexas.edu) of the University of Texas Institute for Geophysics in August 2013. Latest Update: 08/28/2014

Citation: Thirumalai, K., J. W. Partin, C. S. Jackson, T. M. Quinn (2013) Statistical constraints on El Niño Southern Oscillation reconstructions using individual foraminifera: a sensitivity analysis, Paleoceanography, 28, 401-412; doi: 10.1002/palo.20037

2. Input

Input parameters for uncertainty model:

  • num - Number of foraminifera
  • tsl - Number of years represented by sediment sample i.e. length of synthetic time series
  • anerr - Analytical Uncertainty
  • run - Number of Monte Carlo runs for picking loop (Start low [100] and end high [1000-5000])
  • mcs - No. of iterations for overall exercise (mcs = [1-5] - 5 produces a smooth plot; start at 1 and slowly increase; 5 takes time!)

Location of core:

  • lat - Latitude for SST gridpoint (Format: XY.5)
  • latS - Latitude for SSS gridpoint (Format: XY)
  • lon - Longitude for gridpoint (Format: XY.5 - same for SST & SSS)

Domain: [30°N-30°S, 120.5°E-70.5°W; Convention: -ve for S,W])

clear;
clc;
tic
%---- Input Parameters ----

num = 70;
tsl = 100;
anerr = 0.1;
run = 500;
mcs = 1;

%---- Location of Core ----

% Domain: [30°N-30°S, 120.5°E-70.5°W; Convention: -ve for S,W])

lat = -1.5;                     % Data point for VM21-30 (EEP)
latS = -1;
lon = -89.5;

3. Instrumental Data Extraction

Extract SST from HadISST1 (Rayner et al. [2003]) and SSS from the LEGOS dataset (Delcroix et al. [2011]).

TempSal is a workspace file containing the above datasets and it is availble for download here: https://db.tt/Zx4gvxVE

Alternatively, you can load any (equal-length) timeseries for temperature in 'T' and salintiy in 'S' (ex. thermocline) and bypass (or comment out) this block and jump directly to the next.

load('TempSal.mat');            % Load any T/S data through a workspace here

%---- HadISST lat/lon ----

[latco,lonco] = hadco(lat,lon);

%---- Time constraints ----

lyear = 2009;                   % Last year of time series
years = 39;                     % Length of time series: 1970(Jan)-2008(D)
syear = lyear - years;          % Starting year of time series
len = years*12;                 % No. of months
las = (lyear-1870)*12;           % Last month
fir = las - len + 1;            % First month
tmslice = tsl*12;               % Time slice in months

%---- Declaration ----

T = zeros(len,1);               % SST data
S = zeros(len,1);               % SSS data
Y = zeros(len,1);               % 'Y' contains all the months

ys = syear;
c = 0;
for i=1:len                     % Converting the unit of time to years
    if (c<=11)
        Y(i) = ys + (c/12);
        c = c+1;
    else
        c = 0;
        ys = ys+1;
        Y(i) = ys + (c/12);
        c = c+1;
    end
end

%---- Extracting temperature from site ----

for i=fir:las
    T(i-(las-len)) = sstd(lonco,latco,i);
end

%---- GSSS Lat/Lon ----

[latco,lonco] = gssco(latS,lon);

%---- Time constraints ----

las = gssyear(lyear);
fir = gssyear(syear)+1;

%---- Extracting salinity from site ----

k = 1;
for i=fir:las
    S(k) = sssd(lonco,latco,i);
    k = k + 1;
end

4. Foraminiferal Forward Model

Transfer SST and SSS to Pseudo-$\delta^{18}O$ space:

  • Foraminiferal isotopic equilibrium based on low-light equation from Bemis et al. [1998]
  • VSMOW to VPDB-CO2 conversion from Hut [1987]
  • SSS: $\delta^{18}O$ relationship for tropical Pacific is from Legrande & Schmidt [2006]

Change equations according to depth/location (ex. Benway and Mix, 2004)

d18Osw = zeros(len,1);
d18O = zeros(len,1);

for i=1:len
    d18Osw(i) = 0.27*S(i) - 8.88;
    d18O(i) = (d18Osw(i)-0.27) + ((16.5-T(i))/4.8);
end

5. Climatology

Extract seasonal cycle, La Niña and El Niño climatology (beginning with the month of April) from instrumental SST, SSS and $\delta^{18}O$ data.

If your SST and SSS time series is longer than 1970-present, you have to manually add in El Niño years and La Niña years for the record. Alternatively, you can code a function that defines ENSO from the Niño3.4 box and compute particular El Niño and La Niña years.

%---- ENSO Events in Record ----

elyear = [1972;1982;1986;1987;1991;1994;1997;2002;2004;];
layear = [1970;1971;1973;1974;1975;1984;1988;1998;1999;2000;2007;];

elen = (elyear-syear)*12 + 4;           % + 4 => Begin w April
llen = (layear-syear)*12 + 4;

elin = (elyear - syear) + 1;            % Year for injecting ENSO event
lain = (layear - syear) + 1;

%---- Declaration ----

mel = zeros(12,length(elyear));
melT = zeros(12,length(elyear));
melS = zeros(12,length(elyear));

mla = zeros(12,length(layear));
mlaT = zeros(12,length(layear));
mlaS = zeros(12,length(layear));

Cel = zeros(12,1);
Celsd = zeros(12,1);
CelT = zeros(12,1);
CelS = zeros(12,1);

Cla = zeros(12,1);
Clasd = zeros(12,1);
ClaT = zeros(12,1);
ClaS = zeros(12,1);

CT = zeros(12,1);
CS = zeros(12,1);
CC = zeros(12,1);

%---- Yearly Data Extraction ----

%-- El Niño Events --

for i = 1:length(elen)
    mel(:,i) = d18O(elen(i):(elen(i)+11));
    melT(:,i) = T(elen(i):(elen(i)+11));
    melS(:,i) = S(elen(i):(elen(i)+11));
end
elt = reshape(mel,length(elyear)*12,1);
eltT = reshape(melT,length(elyear)*12,1);
eltS = reshape(melS,length(elyear)*12,1);

%-- La Niña Events --

for i = 1:length(llen)
    mla(:,i) = d18O(llen(i):(llen(i)+11));
    mlaT(:,i) = T(llen(i):(llen(i)+11));
    mlaS(:,i) = S(llen(i):(llen(i)+11));
end
lats = reshape(mla,length(layear)*12,1);
latsT = reshape(mlaT,length(layear)*12,1);
latsS = reshape(mlaS,length(layear)*12,1);

%-- Mean Climatologies --

for month=1:12
    Cel(month) = mean(elt(month:12:end));
    Cla(month) = mean(lats(month:12:end));

    Celsd(month) = std(elt(month:12:end));
    Clasd(month) = std(lats(month:12:end));

    CelT(month) = mean(eltT(month:12:end));
    ClaT(month) = mean(latsT(month:12:end));

    CelS(month) = mean(eltS(month:12:end));
    ClaS(month) = mean(latsS(month:12:end));
end

%---- Climatology ----

l = 1;
for month=4:15              % 4-15 => April to March
    CC(l) = mean(d18O(month:12:end));
    CT(l) = mean(T(month:12:end));
    CS(l) = mean(S(month:12:end));
    l = l + 1;
end

%---- ENSO Parameters ----

del = Cel - CC;
dla = Cla - CC;

%---- ENSO Relationship Check ----

if (mean(Cel) < mean(Cla))  % eg. in SouthWestPacific, El Nino => Cold, Dry
    wp = 0;
else
    wp = 1;
end

6. Synthetic Time Series: Modeled Modern Variability

Constructs synthetic time series mimicking 'normal' modern variability. Red and white noise is added and altered according to a while loop check with corresponding tolerances.

[rratai,rfw,rpwf,rwp,rrp] = powrat(d18O);  % Function that extracts spectral parameters (including noise)

%---- Annual:Interannual Power Tolerance ----

aitol = 30;             % Percent tolerance of annual:interannual
mdif = 2*aitol;         % Initial condition
mag = 1;                % Magnitude of ENSO event to ENSO event variability

%---- Noise Model Parameters ----

mnt = 0.75;             % Multiplication component
rnt = 0.25;             % Initial condition for red noise parameter
wnt = 0.25;             % Initial condition for white noise parameter

%---- Noise Model Tolerance ----

wth = 30;               % White noise threshold
rth = 30;               % Red noise threshold
wnp = 2*wth;
rnp = 2*rth;

%---- Modeled Time Series ----

X = 0;
SCheck = zeros(6,5);

basin1 = zeros(12,years);
for i = 1:years
    basin1(:,i) = (CC(:));
end

while ((mdif > aitol) || (wnp > wth) || (rnp > rth))

    X = X + 1;
    conel1 = basin1;

    %---- Injecting ENSO Events ----

    for i = 1:length(elin)
        if (wp == 0)
            conel1(:,elin(i)) = Cel - rand*Celsd*mag;
        else
            conel1(:,elin(i)) = Cel + rand*Celsd*mag;
        end
    end
    for i = 1:length(lain)
        if (wp == 0)
            conel1(:,lain(i)) = Cla + rand*Clasd*mag;
        else
            conel1(:,lain(i)) = Cla - rand*Clasd*mag;
        end
    end

    el18O1 = reshape(conel1,len,1);

    %---- Mimicking Reality ----

    fel = el18O1;

    %-- Noise Addition --

    for i=2:len
        fel(i) = mnt*fel(i) + rnt*fel(i-1) + wnt*(rand);
    end

    %---- Spectral Parameters ----

    [modratai,mfw,mpwf,mwp,mrp] = powrat(fel);
    mdif = (abs(modratai-rratai)/rratai)*100;

    %---- Noise Parameters ----

    rnp = abs((mrp-rrp)/rrp*100);
    wnp = abs((mwp-rwp)/rwp*100);

    %---- Storage ----

    SCheck(X,1) = mdif;
    SCheck(X,2) = wnt;
    SCheck(X,3) = wnp;
    SCheck(X,4) = rnt;
    SCheck(X,5) = rnp;

    %---- Tweaking Noise ----

    if (rnp > rth && ((mrp-rrp)/rrp*100) < -50)     % Red noise check
        rnt = rnt + 0.05;
    elseif (rnp > rth && ((mrp-rrp)/rrp*100) > 50)
        rnt = rnt - 0.05;
    end

    if (wnp > wth && ((mwp-rwp)/rwp*100) < -50)     % White noise check
        wnt = wnt + 0.05;
    elseif (wnp > wth && ((mwp-rwp)/rwp*100) > 50)
        wnt = wnt - 0.05;
    end
end

%---- Making full-length time series ----

mfull = zeros(tmslice,1);

if (tmslice ~= len)
    for i = 1:tmslice
        mfull(i) = el18O1(mod(i,len)+1);
    end
else
    mfull = el18O1;
end

%---- Noise Model Tolerance ----

wnt = SCheck(X,2);
rnt = SCheck(X,4);

for i=2:tmslice     % 'mfull' is modeled modern variability time series
    mfull(i) = mnt*mfull(i) + rnt*mfull(i-1) + wnt*(rand);
end

7. Synthetic Time Series: Modeled Altered (Past) Variability

Construction of synthetic time series with altered ENSO/Seasonal Cycle amplitude: This section also compares picked statistics of altered vs. modeled modern variability and provides a probability value depending on the significance threshold defined for each set of populations.

probR = zeros(20,20,mcs);           % Probability for Range
probS = zeros(20,20,mcs);           % Probability for Std. Dev.

%---- MonteMonte Carlo ----

for mc = 1:mcs
    r = 1;                          % c/w => ENSO; r/p => SeasonalCycle;

    for p=-90:10:100                % SC = kA; k = [.1,2]
        c = 1;

        for w=-90:10:100            % E = kA; k = [.1,2]

            %---- Seasonal Cycle Change ----

            sc = (CC - mean(CC))*(1+(p/100));   % Change in seasonal cycle
            nCC = mean(CC) + sc;                % New (altered) seasonal cycle

            %---- ENSO Amplitude Change ----

            amp = (w/100)+1;            % Change in ENSO amplitude
            nCel = nCC + amp.*del;      % New (altered) El Niño
            nCla = nCC + amp.*dla;      % New (altered) La Niña

            %---- Altered Sin Wave Generation ----

            basin2 = zeros(12,years);
            for i = 1:years
                basin2(:,i) = (nCC(:));
            end
            conel2 = basin2;

            %---- Injecting ENSO Events ----

            for i = 1:length(elin)
                if (wp == 0)
                    conel2(:,elin(i)) = nCel - rand*Celsd*mag;
                else
                    conel2(:,elin(i)) = nCel + rand*Celsd*mag;
                end
            end
            for i = 1:length(lain)
                if (wp == 0)
                    conel2(:,lain(i)) = nCla + rand*Clasd*mag;
                else
                    conel2(:,lain(i)) = nCla - rand*Clasd*mag;
                end
            end

            %---- Full-length Altered Time Series ----

            el18O2 = reshape(conel2,len,1);
            afull = zeros(tmslice,1);

            if (tmslice ~= len)
                for i = 1:tmslice
                    afull(i) = el18O2(mod(i,len)+1);
                end
            else
                afull = el18O2;
            end

            for i=2:tmslice
                afull(i) = mnt*afull(i) + rnt*afull(i-1) + wnt*(rand);
            end

            %---- Monte Carlo Simulations ----

            ts1v = zeros(3,run);
            ts2v = zeros(3,run);

            for m = 1:run                   % Monte Carlo Loop
                MC1 = zeros(num,1);
                MC2 = zeros(num,1);
                for j=1:num                 % Picking Loop
                    k1 = ceil((tmslice)*rand);
                    MC1(j) = mfull(k1) + anerr*randn;
                    k2 = ceil((tmslice)*rand);
                    MC2(j) = afull(k2) + anerr*randn;
                end
                ts1v(1,m) = mean(MC1);
                ts1v(2,m) = std(MC1);
                ts1v(3,m) = max(MC1) - min(MC1);
                ts2v(1,m) = mean(MC2);
                ts2v(2,m) = std(MC2);
                ts2v(3,m) = max(MC2) - min(MC2);
            end

            %---- True SD & Range ----

            o_r1 = Range(mfull);
            o_s1 = std(mfull);

            o_r2 = Range(afull);
            o_s2 = std(afull);

            %---- Extract Monte Carlo SD & Range ----

            sm1 = 2*std(ts1v(1,:));
            ss1 = 2*std(ts1v(2,:));
            sr1 = 2*std(ts1v(3,:));

            sm2 = 2*std(ts2v(1,:));
            ss2 = 2*std(ts2v(2,:));
            sr2 = 2*std(ts2v(3,:));

            ms1 = mean(ts1v(2,:));
            ms2 = mean(ts2v(2,:));
            mr1 = mean(ts1v(3,:));
            mr2 = mean(ts2v(3,:));

            %---- Significance Test ----

            R = 0;
            SD = 0;
            sr3 = sqrt(sr1^2 + sr2^2);  % Significance threshold (Range)
            ss3 = sqrt(ss1^2 + ss2^2);  % Significance threshold (SD)
            for xx=1:run
                for xy=1:run
                    if (abs(ts1v(3,xx)-ts2v(3,xy)) > sr3)
                        R = R + 1;
                    end
                    if (abs(ts1v(2,xx)-ts2v(2,xy)) > ss3)
                        SD = SD + 1;
                    end
                end
                probR(r,c,mc) = R/(run^2)*100;
                probS(r,c,mc) = SD/(run^2)*100;
            end
            c = c +1;
        end
        r = r + 1;
    end
end

8. Output

Probability contour plot for concurrent ENSO vs. Seasonal Cycle amplitude changes.

%---- Averaging MonteMonte Carlo Output ----

pR = zeros(20,20);
pS = zeros(20,20);

for i=1:20
    for j = 1:20
    pR(i,j) = mean(probR(i,j,:));
    pS(i,j) = mean(probS(i,j,:));
    end
end

%---- SigmaPlot Software Output ----

sigmaR = pR';        % Range
sigmaS = pS';        % Std. Dev

pp = (-90:10:100)';
ww = (-90:10:100)';

figure(1);
title('Probability Contour Plot');
clf;

r = subplot(1,2,1);
[cs,h]=contourf(pp,ww,pR,0:10:100,'edgecolor',[0.8 0.8 0.8]);  %[1 0.6 0]
clabel(cs,h,'FontName','Helvetica Neue','fontsize',10,'color',[1 1 1]);
hold on;
h1 = plot(0,0,'*y');
xlim([-90 100]);
ylim([-90 100]);
caxis([0 100]);
colormap(gray);

xlabel('Change in ENSO Amplitude (%)','FontName','Helvetica Neue','fontsize',14);
ylabel('Change in Seasonal Cycle (%)','FontName','Helvetica Neue','fontsize',14);
title(r,'P(Detection)_{Range}','FontName','Helvetica Neue','fontsize',16,'Fontweight','bold');

sd = subplot(1,2,2);
[cs,h]=contourf(pp,ww,pS,0:10:100,'edgecolor',[0.8 0.8 0.8]);
clabel(cs,h,'FontName','Helvetica Neue','fontsize',10,'color',[1 1 1]);
hold on;
plot(0,0,'*y');
xlim([-90 100]);
ylim([-90 100]);
caxis([0 100]);
colormap(gray);

xlabel('Change in ENSO Amplitude (%)','FontName','Helvetica Neue','fontsize',14);
ylabel('Change in Seasonal Cycle (%)','FontName','Helvetica Neue','fontsize',14);
title(sd,'P(Detection)_{Standard Deviation}','FontName','Helvetica Neue','fontsize',16,'Fontweight','bold');
toc

% Yellow star at (0,0) in the contour plot is modeled modern variability.
Elapsed time is 721.720150 seconds.