MATLAB Examples

# INFAUNAL: INdividual Foraminiferal Approach UNcertainty AnaLysis

## Contents

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: 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 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.