In situ determination of porewater gases by ... - Wiley Online Library

Download PDF
6 downloads 16092 Views 545KB Size Report
Comment
set(0,'defaultAxesFontWeight', 'bold')%change default figure font weight ... height]; % as per: http://www.aslo.org/lomethods/instructions/manuscript.html#figures.
Contents Primary Script ................................................................................................................................................................................................................. 2 Flux_Fig.m .................................................................................................................................................................................................................... 20 GasSol.m ...................................................................................................................................................................................................................... 26 InterpTthenD.m ........................................................................................................................................................................................................... 33 PlotContour.m.............................................................................................................................................................................................................. 36

Primary Script % % % % %

POREWATERMIMSDATAWORKUP.m Ryan Bell © SRI International, 2011 Currently at Vancouver Island University permanent email: [email protected] Tested in R2010b

% % % %

This script works up MIMS data collected at depth intervals in sediment porewater. It takes in MIMS data, CTD data and various other meteological data and produces publication quality figures for L&O:Methods

%**** “RawData.xlsx” must be put in the same folder as this script to run****

%%%%% Dependencies:%%%%%% % % Function: Flux_Fig (Included in this Package) % Subtoolbox: \shared\optimlib\lsqcurvefit.m (Mathworks) % Subfunction: tight_subplot.m (Mathworks File Exchange) % Function: GasSol.m (Included in this Package) % Subfunction: CO2SYS.m (http://cdiac.ornl.gov/oceans/co2rprt.html) % Function: InterpDthenT.m (Included in this Package) % Function: PlotContour.m (Included in this Package) % Function: CO2SYS.m (http://cdiac.ornl.gov/oceans/co2rprt.html) % Function: herrorbar.m (Mathworks File Exchange) % Public Toolbox: CSIRO Seawater Tools (http://www.cmar.csiro.au/datacentre/ext_docs/seawater.htm) % Mathworks Toolbox: \curvefit\curvefit\smooth.m (Mathworks) % Mathworks Toolbox: \stats\stats\norminv.m (Mathworks) % Mathworks Toolbox: \stats\stats\regress.m (Mathworks) % Mathworks Toolbox: \stats\stats\zscore.m (Mathworks) % DISCLAIMER: % This software is provided "as is" without warranty of any kind. %% Start with couple global parameters fprintf('Start!...\n'); clear; %Find the folder this .m file is in. Assumes raw data is also in this folder folder = fileparts(mfilename('fullpath')); XLFilename = 'RawData.xlsx';%The file name of the raw data printfigs = 'yes'; %'yes' will cause the figures to be saved into 'folder'

%% Ensure a couple defaults are set up set(0,'defaultAxesFontName', 'times')%change default figure fonts set(0,'defaultTextFontName', 'times')%change default figure fonts set(0,'defaultAxesFontWeight', 'bold')%change default figure font weight set(0,'defaultTextFontWeight', 'normal')%change default figure font weight warning('off', 'MATLAB:interp1:NaNinY');

%the error would otherwise pop up all the time

%% Collect a sample time series for demo purposes MIMSInFile = fullfile(folder,XLFilename); [Data,Headers] = xlsread(MIMSInFile,'ExampleTimeseries'); ExMSTime ExMass32 ExMass34 ExMass40 ExMass44

= = = = =

Data(:,strcmp('Timestamp(days-local)',Headers)); Data(:,strcmp('Mass( 32 )(amps)',Headers)); Data(:,strcmp('Mass( 34 )(amps)',Headers)); Data(:,strcmp('Mass( 40 )(amps)',Headers)); Data(:,strcmp('Mass( 44 )(amps)',Headers));

%% Demostate 32 vs 34 to ID any sulfide presence figure(101);clf; set(gcf,'Name','Mass 32 vs Mass34'); fontsize =12; height = 4; width = 7.25; set(gcf,'units','inches',... 'position',[.05 1 width height],... 'Color','w',... 'InvertHardCopy', 'off',... 'PaperUnits', 'inches',... 'PaperSize', [width height],... 'PaperPositionMode', 'manual',... 'PaperPosition', [0 0 width height]); BG32 = mean(ExMass32(1:100)); BG34 = mean(ExMass34(1:100)); [bhat] = regress((ExMass32-BG32)*1000000000, (ExMass34-BG34)*1000000000); plot(ExMSTime,bhat*(ExMass34-BG34)*1000000000,'b','linewidth',1.5);hold on plot(ExMSTime,(ExMass32-BG32)*1000000000,'k','linewidth',1.5); datetick('x',16); set(gca,'fontsize',fontsize,'linewidth',1.5); set (gca, 'Box', 'off', 'YAxisLocation', 'left'); set(gca, 'tickdir','out'); xlabel('Time (EDT) 08/29/2008','fontsize',fontsize,'fontweight','bold'); ylabel('Ion Current (nA)'); legend('Mass( 34 )','Mass( 32 )','location','NorthWest');

%% BBL Aug 08 Example Data ExampleFig = figure(102);clf; set(gcf,'Name','Aug 08 Example Data'); fontsize =10; height = 2; width = 3.5; %2 cols: 7.25 inches 1 col: 3.5inches FigDim =[.05 1 width height]; % as per: http://www.aslo.org/lomethods/instructions/manuscript.html#figures set(gcf,'units','inches','position',FigDim,'Color','w',... 'InvertHardCopy', 'off'); set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [width height]); set(gcf, 'PaperPositionMode', 'manual'); set(gcf, 'PaperPosition', [0 0 width height]); Range = 4600:5600; plot(ExMSTime,ExMass40*1000000000,'k','linewidth',1.5);hold on plot(ExMSTime,ExMass44*1000000000,'b','linewidth',1.5); plot(ExMSTime,ExMass32*1000000000,'m','linewidth',1.5); datetick('x',16); set(gca,'fontsize',fontsize,'fontweight','bold','linewidth',1.5); set (gca, 'Box', 'off', 'YAxisLocation', 'left'); set(gca, 'tickdir','out'); xlabel('Time (EDT) 08/29/2008','fontsize',fontsize,'fontweight','bold'); ylabel('Ion Current (nA)'); % ind = [0 1 3 4 5 6]; val = [20 154 470 610 742 870]; for n=1:6 text(ExMSTime(val(n)), ExMass32(val(n))*1000000000+1.9, '\downarrow');%,'fontweight','bold') end shg %res: >350dpi min for b/w; 1200pdi for line art; >350 for greyscale if strcmp(printfigs,'yes') print (ExampleFig,fullfile(folder,'Figure3 Exampleseries.tiff'),... '-cmyk','-dtiff','-r400'); end %% Get the raw MIMS Data for individual measuremnts from probe fprintf('Collecting MS Data...\n'); MIMSInFile = fullfile(folder,XLFilename); [Data,Headers] = xlsread(MIMSInFile,'MSIntervalData'); MSTime = MSBin =

Data(:,strcmp('Timestamp(days-local)',Headers)); Data(:,strcmp('Bin',Headers));

BG = MSInitial2Probe = Mass14 = Mass15 = Mass17 = Mass28 = Mass32 = Mass34 = Mass40 = Mass44 =

Data(:,strcmp('Background(amps)',Headers)); Data(:,strcmp('ProbePosition (cm)',Headers)); Data(:,strcmp('Mass( 14 )(amps)',Headers)); Data(:,strcmp('Mass( 15 )(amps)',Headers)); Data(:,strcmp('Mass( 17 )(amps)',Headers)); Data(:,strcmp('Mass( 28 )(amps)',Headers)); Data(:,strcmp('Mass( 32 )(amps)',Headers)); Data(:,strcmp('Mass( 34 )(amps)',Headers)); Data(:,strcmp('Mass( 40 )(amps)',Headers)); Data(:,strcmp('Mass( 44 )(amps)',Headers));

%% Background subract Everything %this is NOT handled in the calibration further down Mass14=Mass14-BG; Mass15=Mass15-BG; Mass17=Mass17-BG; Mass28=Mass28-BG; Mass32=Mass32-BG; Mass34=Mass34-BG; Mass40=Mass40-BG; Mass44=Mass44-BG; %% Collect the ctd data fprintf('Collecting CTD Data...\n'); InFile = fullfile(folder,XLFilename); [Data,Headers] = xlsread(InFile,'CTDData'); CTDTime = CTDT = CTDCond = CTDPress =

Data(:,strcmp('Timestamp(days-local)',Headers)); Data(:,strcmp('Temp (°C)',Headers)); Data(:,strcmp('Cond(mS/cm)',Headers)); Data(:,strcmp('Pressure (dB)',Headers));

%% Derive CTD related parameters %Calcualte Gas solubilities CTDDepth = sw_dpth(CTDPress, 31.375); %http://www.cmar.csiro.au/datacentre/ext_docs/seawater.htm CTDS =real(sw_salt(CTDCond/sw_c3515,CTDT,CTDPress)); %http://www.cmar.csiro.au/datacentre/ext_docs/seawater.htm [CTDCH4,~,~,~,~] = GasSol(CTDS,CTDT,0.000001745*ones(size(CTDT)),0); [~,CTDN2 ,~,~,~] = GasSol(CTDS,CTDT,0.781*ones(size(CTDT)),0); [~,~,CTDO2 ,~,~] = GasSol(CTDS,CTDT,0.209*ones(size(CTDT)),0); [~,~,~,CTDAr ,~] = GasSol(CTDS,CTDT,0.00934*ones(size(CTDT)),0); [~,~,~,~,CTDCO2] = GasSol(CTDS,CTDT,0.000385*ones(size(CTDT)),0); %interpolate the important CTD data to the same timestamps as the MIMS Data

AmbD = interp1(CTDTime,CTDDepth,MSTime); AmbT = interp1(CTDTime,CTDT,MSTime); AmbS = interp1(CTDTime,CTDS,MSTime); AmbP = interp1(CTDTime,CTDPress,MSTime); AmbAr = interp1(CTDTime,CTDAr,MSTime); %These data are not used, but are an interesting comparison to MISM Data AmbCH4Sat = interp1(CTDTime,CTDCH4,MSTime); AmbN2Sat = interp1(CTDTime,CTDN2,MSTime); AmbO2Sat = interp1(CTDTime,CTDO2,MSTime); AmbCO2Sat = interp1(CTDTime,CTDCO2,MSTime); %Define a couple other necessary parameters AmbMIMST = 45*ones(size(AmbT)); %membrane temperature was 45C AmbALK = 2389*ones(size(AmbT)); %Alkalintiy Measured from SSW collected ar R2 on this day

%% Obtain the Calibration coefficients InFile = fullfile(folder,XLFilename); [Data,~,Raw] = xlsread(InFile,'MIMSCalData'); Headers = Raw(1,2:end);% remove first row to properly match with Data b14 b15 b17 b28 b32 b40 b44

= = = = = = =

Data(:,strcmp('Nitrogen(14)',Headers)); Data(:,strcmp('Methane',Headers)); Data(:,strcmp('Water',Headers)); Data(:,strcmp('Nitrogen(28)',Headers)); Data(:,strcmp('Oxygen',Headers)); Data(:,strcmp('Argon',Headers)); Data(:,strcmp('Carbon Dioxide',Headers));

%% B0 Baseline Correction _ m/z 17 fprintf('Performing Water correction on baseline...\n'); %create acorrected baseline for each data point H2OCorr = b17(1)./Mass17; %predicted water over measured water b14_0 b15_0 b28_0 b32_0 b40_0 b44_0

= = = = = =

b14(1)./H2OCorr; b15(1);%./H2OCorr; b28(1)./H2OCorr; b32(1)./(H2OCorr); b40(1);%./H2OCorr; b44(1)./H2OCorr;

%m/z 15 not affected by water %m/z 40 not affected by water

%% B1 Sensitivity Correction_ m/z 40 fprintf('Argon Correcting Data...\n');

%obtain Argon Data from MS that is comparable to Ambient saturated Ar %use this to create a correction factor, then correct the other gases %Calculate the Argon conceatrions measured by MIMS ArConc = (Mass40-b40_0)./b40(2); %obtain the argon measurements of topmost determinations (above sediment) AmbMSAr = ArConc(MSInitial2Probe==0); %obtain time of measurements AmbMSTime = MSTime(MSInitial2Probe==0); %obtain calculated saturation Argon concentrations from CTD AmbArSol = AmbAr(MSInitial2Probe==0); %calculate Argon correction coefficitn: predicted Argon over measured Argon ArCorr = AmbArSol./AmbMSAr; %interpolate data to match time stamps of all measuremns, including the %'below sediment' measuremtns ArCorr = interp1(AmbMSTime,ArCorr,MSTime,'cubic','extrap'); % plot er up for QA/QC figure(3);clf; plot(MSTime,ArCorr); Title('Argon Correction Factor'); shg; %create a variable for corrected sensitivity for each data point b14_1 = b14(2)./ArCorr; b15_1 = b15(2)./ArCorr; b28_1 = b28(2)./ArCorr; b32_1 = b32(2)./ArCorr; %b34_1 = b34(2)./ArCorr; b40_1 = b40(2)./ArCorr; b44_1 = b44(2)./ArCorr; %% B2 correction

(interference factors)

%B2 does not need a correction, it's intrinsic to the ion source %create a variable each data point for consistencies sake b14_2 = b14(3)*ones(size(b14_1)); b15_2 = b15(3)*ones(size(b14_1)); b28_2 = b28(3)*ones(size(b14_1)); b32_2 = b32(3)*ones(size(b14_1)); b40_2 = b40(3)*ones(size(b14_1)); b44_2 = b44(3)*ones(size(b14_1)); %% Calibrate the data

N2Conc1Corr = (((Mass14-b14_2.*Mass15))-b14_0)./(b14_1); %include interfence factor from maas15 CH4ConcCorr = (((Mass15-b15_2.*Mass14))-b15_0)./(b15_1); %include interfence factor from maas14 N2Conc2Corr = (((Mass28-b28_2.*Mass44))-b28_0)./b28_1; %include interfence factor from maas44 O2ConcCorr = ((Mass32)-b32_0)./b32_1; % by definition, topmost Ar measurement will be equal to the calcuated Ar saturation conc, ArConcCorr = ((Mass40)-b40_0)./b40_1; CO2ConcCorr = ((Mass44)-b44_0)./b44_1; N2overAr1 = N2Conc1Corr./ArConc; N2overAr2 = N2Conc2Corr./ArConc; %% Correction CO2 for from in MIMS T to situ T %!!Assumes T and S are the same in the porewater as in the ambient water % potentially minor differences will have a very small effect (negible) fprintf('Correction CO2 for Temp shift...\n'); %example calculation at correct conditions to obatin Henry's Coefficient [Data,~,~]=CO2SYS(380,AmbALK,5,1,AmbS,AmbMIMST,AmbMIMST,AmbP,AmbP,0,0,1,4,1); % Henry's Coeffient K0 = Data(:,5)./Data(:,8); %calcualte K0 = fug(uatm)/co2*(umol/kg) %fugacity of measured CO2 measured at the MIMS interface fCO2 = K0.*CO2ConcCorr; %calculate fCO2 %use fugacity measurement to calculate [CO2] in situ [Data,~,~]=CO2SYS(fCO2,AmbALK,5,1,AmbS,AmbMIMST,AmbT,AmbP,AmbP,0,0,1,4,1); %Calculate Co2 in MIMS then in situ CO2starin2 = Data(:,8) ; %mims (same as above) CO2starout2 = Data(:,23); %in situ %redefine [Co2] for the in situ conditions CO2ConcCorr = CO2starout2; %% Remove Negative Concentrations N2Conc1Corr(N2Conc1Corr(n)Boundary(n) %take all previous data, first point of this profile, %boudnary of this profile, then rest of this profile DepthPro{m}= [Boundary(n); Depth(n:n+6)]; MSTimePro{m} = [interp1(Depth(n-1:n+1), MSTime(n-1:n+1), Boundary(n)); MSTime(n:n+6)]; MSBinPro{m} = [MSBin(n) ; MSBin(n:n+6)]; N2Conc1CorrPro{m} = [N2Conc1Corr(n-1);N2Conc1Corr(n:n+6)]; CH4ConcCorrPro{m} = [CH4ConcCorr(n-1);CH4ConcCorr(n:n+6)]; N2Conc2CorrPro{m} = [N2Conc2Corr(n-1);N2Conc2Corr(n:n+6)]; O2ConcCorrPro{m} = [O2ConcCorr(n-1);O2ConcCorr(n:n+6)]; ArConcCorrPro{m} = [ArConcCorr(n-1);ArConcCorr(n:n+6)]; CO2ConcCorrPro{m} = [CO2ConcCorr(n-1);CO2ConcCorr(n:n+6)];

else %take all previous %point of profile, %of the profile DepthPro{m}= MSTimePro{m} = MSBinPro{m} = N2Conc1CorrPro{m}= CH4ConcCorrPro{m}= N2Conc2CorrPro{m}= O2ConcCorrPro{m} = ArConcCorrPro{m} = CO2ConcCorrPro{m}=

data, first point of profile, second boundary point of profile, then the rest [Boundary(n); Depth(n+1:n+6)]; [interp1(Depth(n:n+2), MSTime(n:n+2), Boundary(n)); MSTime(n+1:n+6)]; [MSBin(n) ; MSBin(n+1:n+6)]; [N2Conc1Corr(n); N2Conc1Corr(n+1:n+6)]; [CH4ConcCorr(n);CH4ConcCorr(n+1:n+6)]; [N2Conc2Corr(n); N2Conc2Corr(n+1:n+6)]; [O2ConcCorr(n); O2ConcCorr(n+1:n+6)]; [ArConcCorr(n); ArConcCorr(n+1:n+6)]; [CO2ConcCorr(n); CO2ConcCorr(n+1:n+6)];

end end end %% Check out some fits fprintf('Fit the profiles to a model...\n'); ConcData.Time=MSTimePro; ConcData.Depth=DepthPro; ConcData.CH4=CH4ConcCorrPro; ConcData.N2=N2Conc1CorrPro; ConcData.O2=O2ConcCorrPro; ConcData.CO2=CO2ConcCorrPro; %may need to commented off for licencing issues fontsize=9; height = 2; width = 3.5; %2 cols: 7.25 inches 1 col: 3.5inches %optimization toolbox required: [DepResFig FitFig CH4Beta,O2Beta]=Flux_Fig(ConcData,fontsize,height,width); FileName='Modelled Profile'; %res: >350dpi min for b/w; 1200pdi for line art; >350 for greyscale if strcmp(printfigs,'yes') print (FitFig,'-dtiff','-r400','-cmyk',... fullfile(folder,'Figure7 PooledData.tiff')); print (DepResFig,'-dtiff','-r400','-cmyk',... fullfile(folder,'Figure8 Production.tiff')); end %% A plot indicating Statistcal relevance of the methane maximum %Is methane data >4cm statically higher than below, what about maxicline?

n=9; m =numel(DepthInt)/n; depthArr = reshape(DepthInt,n,m); DepthHat=mean(depthArr,2); CH4ConcCorrIntArr = reshape(CH4ConcCorrInt,n,numel(CH4ConcCorrInt)/n); CH4Hat = mean(CH4ConcCorrIntArr,2); CH4std = std(CH4ConcCorrIntArr,0,2); confidence = .95; p = max(norminv([(1-confidence)/2, 1-((1-confidence)/2)])); %95 condifence

figure(99);clf; set(gcf,'name','CH4 statistical analysis') title('Methane Profile with Errorbars') %From the file exchange: http://www.mathworks.com/matlabcentral/fileexchange/3963-herrorbar herrorbar(CH4Hat,DepthHat,(p*CH4std)./(sqrt(m)),'k'); hold on set(gca,'ydir','reverse') ylabel('Depth (cm)');xlabel('[CH_4]'); plot(CH4ConcCorrInt,DepthInt,'.','color',[.7 .7 .7]) %% Interpolate the 1D data timeseries to 2D for contour plot %cubic interpolation used, linear give super boring difference data %(ie first derivative of a linear line is constant) tsteps=size(MSTimeInt,1)*5; dsteps=2*180; %methods: 'nearest','linear','spline','pchip','cubic','v5cubic' DepthType = 'cubic'; TimeType = 'linear'; fprintf('Manually Interpolating the N2-14 Data...\n'); [N2_1 , ~, ~] = InterpDthenT(MSTimeInt, DepthInt, N2Conc1CorrInt,... MSBinInt, tsteps, dsteps, DepthType, TimeType); fprintf('Manually Interpolating the CH4 Data...\n'); [CH4, ~, ~] = InterpDthenT(MSTimeInt, DepthInt, CH4ConcCorrInt, ... MSBinInt, tsteps, dsteps, DepthType, TimeType); fprintf('Manually Interpolating the N2-28 Data...\n'); [N2_2, ~, ~] = InterpDthenT(MSTimeInt, DepthInt, N2Conc2CorrInt, ... MSBinInt, tsteps, dsteps, DepthType, TimeType); fprintf('Manually Interpolating the O2 Data...\n'); [O2, ~, ~] = InterpDthenT(MSTimeInt, DepthInt, O2ConcCorrInt, ... MSBinInt, tsteps, dsteps, DepthType, TimeType); fprintf('Manually Interpolating the Ar Data...\n'); [Ar, ~, ~] = InterpDthenT(MSTimeInt, DepthInt, ArConcCorrInt, ... MSBinInt, tsteps, dsteps, DepthType, TimeType);

fprintf('Manually Interpolating the CO2 Data...\n'); [CO2 T D] = InterpDthenT(MSTimeInt, DepthInt, CO2ConcCorrInt, ... MSBinInt, tsteps, dsteps, DepthType, TimeType); %% Find Chemoclines fprintf('Find Chemoclines...\n'); %location of mean O2 concentration for n = 1:size(O2,2) meanO2 = (max(max(O2))+min(min(O2)))/2; [g i]=min(abs(O2(:,n)-meanO2)); OxyMean(n)=interp1(O2(i-1:i+1,n),D(i-1:i+1,n),meanO2); end %oxygen chemocline O2Diff = diff(O2,1,1); Ddiff = D(1:end-1,:) + diff(D,1,1); %centered location of each differential [~,j] = max(abs(O2Diff),[],1); O2Chemocline = Ddiff(j); %CO2 Chemocline CO2Diff = diff(CO2,1,1); Ddiff = D(1:end-1,:) + diff(D,1,1); %centered location of each differential [~,j] = max(abs(CO2Diff),[],1); CO2Chemocline = Ddiff(j); %CO2Maxocline [~,i] = max(CO2,[],1); CO2Maxocline = D(i); % figure(10);clf;hold on set(gcf,'position',[20 20 600 400],'name','Chemocline demo') plot(T(1,:),O2Chemocline,'.') %plot(T(1,:),smooth(O2Chemocline,20),'.k') plot(T(1,:),CO2Chemocline,'.g') plot(T(1,:),CO2Maxocline,'.r') plot(T(1,:),OxyMean,'.m') legend('O2Chemocline','CO2Chemocline','CO2Maxocline','OxyMean') set(gca,'ydir','reverse') datetick('x',2) shg fprintf(1, 'The mean O2 Chemocline is at %g cm\n', mean(O2Chemocline)); fprintf(1, 'The mean CO2 Chemocline is at %g cm\n', mean(CO2Chemocline));

fprintf(1, 'The mean difference between the two is: %g cm\n', mean(O2Chemocline-CO2Chemocline)); fprintf(1, 'The mean Depth of the CO2 maximum is %g cm\n', mean(CO2Maxocline)); %%

Some time series plots for chemocline comparisons

fprintf('Plot Timeseries model...\n'); Chemofig=figure(6);clf set(gcf,'Name','BBL Aug 08 ChemoCline'); height = 2.5; width = 6; %2 cols: 7.25 inches 1 col: 3.5inches FigDim =[.5 1 width height]; % as per: http://www.aslo.org/lomethods/instructions/manuscript.html#figures set(gcf,'units','inches','position',FigDim,'Color','w', 'InvertHardCopy', 'off'); set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [width height]); set(gcf, 'PaperPositionMode', 'manual'); set(gcf, 'PaperPosition', [0 0 width height]); FontWeight = 'bold'; FontName = 'times'; FontMag = 10; %%%%%%%%%%%%%%%%% Set up the x axes %%%%%%%%%%%%%%%% %first the times timeint = 1/3; %time intervals for x labels in days timemax=timeint*(floor((max(T(1,:)))/timeint)); timemin=timeint*(ceil((min(T(1,:)))/timeint)); xlim=[(min(T(1,:))) (max(T(1,:)))]; xtick = timemin:timeint:timemax; xticklabel = datestr(xtick, 15); %%%%%%%%%Set up the data%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Prpeare some data for multivariate analysis %use a linearly space time range, over a limited range normtime = linspace(xlim(1), xlim(2), 1000); %Assignt the variables to be analyzed, and interpolate them on normtime interface =interp1(FritTime,Initial2Sand,normtime,'linear','extrap'); chemo = interp1(T(1,:),O2Chemocline,normtime,'linear','extrap'); %chemo = interp1(T(1,:),OxyMean,normtime,'linear','extrap'); depth= interp1(CTDTime(1:10:length(CTDTime)),CTDDepth(1:10:length(CTDTime)),normtime,'linear','extrap'); Par3 = interp1(ParTime,Par,normtime,'linear','extrap'); %interpolate over NANs interface = interp1(normtime(~isnan(interface)),interface(~isnan(interface)),normtime,'linear','extrap'); nans

%remove

chemo = interp1(normtime(~isnan(chemo)),chemo(~isnan(chemo)),normtime,'linear','extrap')'; %remove nans depth = interp1(normtime(~isnan(depth)),depth(~isnan(depth)),normtime,'linear','extrap'); %remove nans Par3 = interp1(normtime(~isnan(Par3)),Par3(~isnan(Par3)),normtime,'linear','extrap'); %remove nans %Standardize interface = zscore(interface); chemo = zscore(chemo); depth = zscore(depth); Par3 = zscore(Par3); %Detrend %uncomment this section to observe multiregression for removeing %long-term trends (ie constant plunging of the probe) % interface = detrend(interface); % chemo = detrend(chemo); % depth = detrend(depth); % Par3 = detrend(Par3); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% make the first plot%%%%%%%%%%%%%%%%%%%% color1='k'; ah1=axes; set(ah1,'fontsize',FontMag,'fontweight',FontWeight,'fontname',FontName); p1=plot(ah1,normtime,chemo,'color',color1,'linewidth',3);hold on; ax1pos = get(ah1,'position'); set(ah1,'position',[ax1pos(1)+.15 ax1pos(2) ax1pos(3)*.8 ax1pos(4)]); set (ah1, 'Box', 'off', 'Color', 'none', 'YAxisLocation', 'right','linewidth',1.5,'fontweight',FontWeight,'fontname',FontName); set(ah1, 'tickdir','out','ycolor',color1,'YMinorTick','on');%,'ydir','reverse'); l5 = ylabel('O_2 Chemocline (cm)','color',color1,'fontsize',FontMag,'fontweight',FontWeight,'fontname',FontName); set(ah1,'xlim',xlim,'xtick',xtick,'xticklabel',xticklabel); xlabel('Time (EDT)'); % and now the dates timeint=1; %Time interval of label - 1 equals once a day timemax=timeint*(floor((max(T(1,:)))/timeint)); %max label timemin=timeint*(ceil((min(T(1,:)))/timeint)); %min label xtickb = timemin:timeint:timemax ; %create the ticks xticklabelb = datestr(xtickb, 2) ; %create the labels ylimits = get(ah1,'ylim'); ydateposition = ylimits(1)- .152*abs(ylimits(1)-ylimits(2)); for n=1:size(xticklabelb,1) text(xtickb(n),ydateposition ,datestr(xtickb(n), 2),'units','data',... 'fontsize',FontMag*1.1,'fontweight',FontWeight,'HorizontalAlignment','center','fontname',FontName) end %%%%%%%%%%%%% a tricky plot

%%%%%%%%%%%%%%%%%%%%%%

ax1pos = get(ah1,'position'); ah3=axes('position',ax1pos); set(ah3,'fontsize',FontMag,'fontweight',FontWeight,'fontname',FontName); p3=plot(ah3,normtime,Par3,'r','linewidth',1.5);hold on set(ah3,'xlim',xlim,'visible','off','YLimMode','manual'); %set (ah3, 'Box', 'off', 'Color', 'none', 'XTick', [], 'YAxisLocation', 'left', 'tickdir','out','linewidth',1.5,'fontweight',FontWeight,'fontname',FontName); ylim3=get(ah3,'ylim'); ah3b=axes('position',[ax1pos(1)-.1 ax1pos(2) ax1pos(3)*.8 ax1pos(4)]); set (ah3b, 'Box', 'off', 'XTick', [], 'YAxisLocation', 'left', 'tickdir','out','linewidth',1.5,'fontweight',FontWeight,'fontname',FontName); set(ah3b,'xlim',xlim,'Color', 'none'); set(ah3b,'ycolor','r','ylim',ylim3,'YMinorTick','on','xcolor','w'); set(ah3b,'Layer','bottom') ylabel('APAR (\muE/m^2s)','color','r','fontsize',FontMag,'fontweight',FontWeight,'fontname',FontName); %%%%%%%%%%%%% another tricky plot %%%%%%%%%%%%%%%%%%%%%% color4 = [0 .8 0]; ax1pos = get(ah1,'position'); ah4=axes('position',ax1pos); set(ah4,'fontsize',FontMag,'fontweight',FontWeight,'fontname',FontName); p4=plot(ah4,normtime,interface,'color',color4,'linewidth',1.5);hold on set(ah4,'xlim',xlim,'visible','off','YLimMode','manual','ydir','reverse'); %set (ah4, 'Box', 'off', 'Color', 'none', 'XTick', [], 'YAxisLocation', 'left', 'tickdir','out','linewidth',1.5,'fontweight',FontWeight,'fontname',FontName); ylim4=get(ah4,'ylim'); ah4b=axes('position',[ax1pos(1)-.2 ax1pos(2) .01 ax1pos(4)]); set (ah4b, 'Box', 'off', 'XTick', [], 'YAxisLocation', 'left', 'tickdir','out','linewidth',1.5,'fontweight',FontWeight,'fontname',FontName); set(ah4b,'xlim',xlim,'Color', 'none'); set(ah4b,'ycolor',color4,'ylim',ylim4,'YMinorTick','on','xcolor','w','ydir','reverse'); ylabel('Sediment Interface (cm)','color',color4,'fontsize',FontMag,'fontweight',FontWeight,'fontname',FontName); %%%%%%%%%%%%%now the next plot %%%%%%%%%%%%%%%%%%%%%% color2='b'; ax1Pos = get(ah1,'position'); ah2=axes('position',ax1pos); set(ah2,'fontsize',FontMag,'fontweight',FontWeight,'fontname',FontName); p2=plot(ah2,normtime,depth,'color',color2,'linewidth',1.5);hold on set(ah2,'xlim',xlim); set (ah2, 'Box', 'off', 'Color', 'none', 'XTick', [], 'YAxisLocation', 'left', 'tickdir','out','linewidth',1.5,'fontweight',FontWeight,'fontname',FontName); set(ah2,'ycolor',color2,'YMinorTick','on'); ylabel('Depth (m)','color',color2,'fontsize',FontMag,'fontweight',FontWeight,'fontname',FontName);

%res: >350dpi min for b/w; 1200pdi for line art; >350 for greyscale if strcmp(printfigs,'yes') print (Chemofig,fullfile(folder,'Figure6 Multivariate.tiff'),... '-cmyk','-dtiff','-r400'); end %%

multivariate Analysis

%Assign predictors Y=-chemo; X=[ones(size(Par3)); interface; Par3; depth]'; %perform multiregression [B,BINT,R,RINT,stats] = regress(Y,X); %announce results fprintf(1, 'Interface Beta is %g cm\n', B(2)); fprintf(1, 'Par Beta is %g cm\n', B(3)); fprintf(1, 'Depth(tide) Beta is %g cm\n', B(4)); %% Set up the figure Plot Contour ContourFig=figure(5);clf; set(gcf,'Name','BBL Aug 08 Contour Plots'); height = 6;%7.5; width = 5;%6; %2 cols: 7.25 inches 1 col: 3.5inches % as per: http://www.aslo.org/lomethods/instructions/manuscript.html#figures FigDim =[.05 1 width height]; set(gcf,'units','inches','position',FigDim,'Color','w', 'InvertHardCopy', 'off'); set(gcf, 'PaperUnits', 'inches'); set(gcf, 'PaperSize', [width height]); set(gcf, 'PaperPositionMode', 'manual'); set(gcf, 'PaperPosition', [0 0 width height]);

C.N2_1=N2_1; C.CH4=CH4; C.N2_2=N2_2; C.O2=O2; C.Ar=Ar; C.CO2=CO2; PlotContour(C,T,D,MSTime,Depth,MSTime,Boundary,ContourFig); %res: >350dpi min for b/w; 1200pdi for line art; >350 for greyscale if strcmp(printfigs,'yes')

print(ContourFig,fullfile(folder,'Figure5 Contours.tiff'),... '-cmyk','-dtiff','-r400'); end %% Done fprintf('Done!\n'); beep

Flux_Fig.m function [DepResFig FitFig CH4Beta O2Beta]=Flux_Fig(ConcData,fontsize,height,width) % Flux_Fig performs production analysis and plotting on oxygen and methane % depth profile data % Ryan Bell © SRI International, 2011 % Currently at Vancouver Island University % permanent email: [email protected] % Works in R2010b % DISCLAIMER: % This software is provided "as is" without warranty of any kind. %========================================================================= %% %% setup the figure for saving etc, FitFig = figure(7);clf; set(FitFig,'Name','BBL Aug 08 Fit Methane and O2'); % height = 3; % width = 6; % %FigDim =[16 -1.2 width height]; % locate it on second screen on right FigDim =[1 1 width height]; % locate it on first screen set(FitFig,'units','inches','position',FigDim,'Color','w', 'InvertHardCopy', 'off'); set(FitFig, 'PaperUnits', 'inches'); set(FitFig, 'PaperSize', [width height]); set(FitFig, 'PaperPositionMode', 'manual'); set(FitFig, 'PaperPosition', [0 0 width height]); %%

Sort out Data, and create some globally used variables

D = vertcat(ConcData.Depth{:}); methane = vertcat(ConcData.CH4{:}); oxygen = vertcat(ConcData.O2{:}); j = [D, oxygen, methane]; k=sortrows(j,1); D=-k(:,1); % depth vector positive up oxygen = k(:,2); methane = k(:,3);

linewidth=2; DPred=0:-.1:-15; %depth data onto which to predict concentration data phi = .37; % Rao thetasq = 3; % Bourdreau KmolO2= 9.9E-6; %cm^2/sec -- NOTE UNITS KmolO2 = KmolO2 * 60*60*24; %cm^2/day alpha0 = 2; % units are day-1 HalvingDepth = 2; % Jahnke best fit had a "halving depth" of 2 cm alphaZ = alpha0*2.^(DPred/HalvingDepth); %% Setup subplots %from file exchange: %http://www.mathworks.com/matlabcentral/fileexchange/27991-tight-subplot ha = tight_subplot(1,2,[.00 0.07],[.23 .03],[.15 .03]); %file exchange ax1=ha(1);ax2=ha(2); %% Oxygen %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ax1=axes('position',[.13 .22 .33 .7]); LB=[0 0 2]; %set lower bound limits for the fit to avoid some BS due to data sparsity UB =[inf inf inf]; %set upper bound limits for the fit to avoid some BS due to data sparsity beta0=[0.004 0.001 3]; % initial guess beta1 = lsqcurvefit(@ProfileDataModelO2,beta0,D,oxygen,LB,UB); O2Beta = beta1; %store for output %Calculate preduicted profiles OxyHat = ProfileDataModelO2(beta1,DPred); OxyHat1st = ProfileDataModelO2_1stDer(beta1,DPred); OxyHat2nd = ProfileDataModelO2_2ndDer(beta1,DPred); axes(ax1); plot(ax1,oxygen,D,'+k','markersize',2,'color',[.6 .6 .6]);hold on plot(ax1,OxyHat,DPred,'k','linewidth', linewidth);hold on %plot(OxyHat1st,DPred,'r','linewidth', linewidth) %plot(OxyHat2nd,DPred,'g','linewidth', linewidth) %Calculate resulting flux KO2 = (phi/thetasq)*KmolO2*OxyHat2nd; % O2secder units = (umol/L)/cm^2 %overall units = (umol/L)/day XO2Z = (OxyHat(1) - OxyHat); %units = (umol/L) nonlocalO2 = XO2Z.*alphaZ ; %alphaZ units = 1/day; overall (umol/L)/day NLO2 = -(KO2+nonlocalO2); %(umol/L)/day %Claculate net O2 consumption O2cons = NLO2*phi*0.001; %mmol/(m^2 * 0.1 cm)/day [~,i]=min(O2cons); [~,j]=max(OxyHat1st);

fprintf(1, 'The max O2 consumption depth is %g cm\n', DPred(i)) fprintf(1, 'The average O2 chemocline depth is %g cm\n', DPred(j)) netO2cons = sum(O2cons); %mmol O2/m^2/day (over the top 15 cm)) fprintf(1, 'The total O2 production

is %g mmol m-2 d-1 (over the top 15 cm)\n', netO2cons)

ylabel(ax1,'Depth (cm)','fontname','times','fontweight','bold','fontsize',fontsize); xlabel(ax1,'O_{2} \mumol/kg','fontname','times','fontweight','bold','fontsize',fontsize); set(ax1,'fontname','times','fontweight','bold','fontsize',fontsize); %% Methane %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ax2=axes('position',[.63 .22 .33 .7]); axes(ax2) %obtain fit coeffiecents for average of all profiles beta0=[0.0051 0.0014 0.0017 0.27 0.048]; % good intial guesses for coeffieicits beta1 = lsqcurvefit(@ProfileDataModel,beta0,D,methane); %obtain best fit CH4Beta=beta1; %report the best fit (and store in new variable)

%obtain predicted data from best fit MethHat = ProfileDataModel(beta1,DPred); %obtain predicted data MethHat1st=ProfileDataModel_1stDer(beta1,DPred);%obtain predicted differential data MethHat2nd=ProfileDataModel_2ndDer(beta1,DPred);%obtain predicted differential data KCH4 = (phi/thetasq)*KmolO2*MethHat2nd; % O2secder units = (umol/L)/cm^2 %overall units = (umol/L)/day XCH4Z = (MethHat(1) - MethHat); %units = (umol/L) nonlocalCH4 = XCH4Z.*alphaZ; %alphaZ units = 1/day; overall (umol/L)/day NLCH4 = -(KCH4+nonlocalCH4); %(umol/L)/day %Claculate net O2 consumption CH4cons = NLCH4*phi*1000; %nmol/(m^2 * 0.1 cm)/day [~,i]=max(CH4cons); [~,j]=min(MethHat1st); fprintf(1, 'The max CH4 consumption depth is %g cm\n', DPred(i)) fprintf(1, 'The average CH4 chemocline depth is %g cm\n', DPred(j)) netCH4cons = sum(CH4cons); %nmol CH4/m^2/day (over the top 15 cm) fprintf(1, 'The total CH4 production is %g nmol m-2d-1 (over the top 15 cm)\n', netCH4cons)

%plot the data plot(ax2, methane,D,'+','color',[.6 .6 .6],'markersize',2);hold on plot(ax2, MethHat,DPred,'k','linewidth', linewidth)

%ylabel(ax2, 'Depth (cm)','fontname','times','fontweight','bold','fontsize',fontsize); xlabel(ax2, 'CH_{4} \mumol/kg','fontname','times','fontweight','bold','fontsize',fontsize); set(ax2,'fontname','times','fontweight','bold','fontsize',fontsize); set(ax2,'yticklabel',[]) shg %% Another Figure for the Depth resolved prod/cons % setup the figure for saving etc, DepResFig = figure(8);clf; set(DepResFig,'Name','BBL Aug 08 Depth Resolved Methane and O2'); % % %

% %

height = 3; width = 6; %FigDim =[16 -1.2 width height]; % locate it on second screen on right FigDim =[1 1 width height]; % locate it on first screen set(DepResFig,'units','inches','position',FigDim,'Color','w', 'InvertHardCopy', 'off'); set(DepResFig, 'PaperUnits', 'inches'); set(DepResFig, 'PaperSize', [width height]); set(DepResFig, 'PaperPositionMode', 'manual'); set(DepResFig, 'PaperPosition', [0 0 width height]); ha = tight_subplot(1,2,[.00 0.07],[.23 .03],[.15 .03]); for ii = 1:6; axes(ha(ii)); plot(randn(10,ii)); end set(ha(1:4),'XTickLabel',''); set(ha,'YTickLabel','') %ax1=axes('position',[.15 .2 .3 .7]); ax1=ha(1);ax2=ha(2); plot(ax1,O2cons,DPred,'k','linewidth', linewidth) ylabel(ax1,'Depth (cm)','fontname','times','fontweight','bold','fontsize',fontsize); xlabel(ax1,'O_{2} mmolm^-^2d^-^1','fontname','times','fontweight','bold','fontsize',fontsize); set(ax1,'fontname','times','fontweight','bold','fontsize',fontsize);

% ax2=axes('position',[.65 .2 .3 .7]); plot(ax2,CH4cons,DPred,'k','linewidth', linewidth) %ylabel('Depth (cm)','fontname','times','fontweight','bold','fontsize',fontsize); xlabel(ax2,'CH_{4} nmolm^-^2d^-^1','fontname','times','fontweight','bold','fontsize',fontsize); set(ax2,'fontname','times','fontweight','bold','fontsize',fontsize); set(ax2,'yticklabel',[]) end

function yhat=ProfileDataModel(beta, x) a=beta(1); b=beta(2); c=beta(3);

d=beta(4); e=beta(5); yhat=(a+b*x + c*x.^2)./(1+d*x +e*x.^2); end function yhat=ProfileDataModel_1stDer(beta, x) % %diff of: yhat=(a+b*x + c*x.^2)./(1+d*x +e*x.^2); % to obtain this type: % syms x a b c d e; % result = diff((a+b*x+c*x^2)/(1+d*x+e*x^2),x) a=beta(1); b=beta(2); c=beta(3); d=beta(4); e=beta(5); yhat=(b + 2*c*x)./(e*x.^2 + d*x + 1) - ((d + 2*e*x).*(c*x.^2 + b*x + a))./(e*x.^2 + d*x + 1).^2; end function yhat=ProfileDataModel_2ndDer(beta, x) %diff of: yhat=(b + 2*c*x)./(e*x.^2 + d*x + 1) - ((d + 2*e*x).*(c*x.^2 + b*x + a))./(e*x.^2 + d*x + 1).^2; %to obtain this type: %syms x a b c d e; %result = diff(yhat=(b + 2*c*x)./(e*x.^2 + d*x + 1) - ... %((d + 2*e*x).*(c*x.^2 + b*x + a))./(e*x.^2 + d*x + 1).^2,x) a=beta(1); b=beta(2); c=beta(3); d=beta(4); e=beta(5); yhat = ((2*b*e^2-2*c*d*e).*x.^3+(6*a*e^2-6*c*e)*x.^2+(6*a*d-6*b)*e.*x-2*a*e+2*a*d^2-2*b*d+2*c)... ./(e^3*x.^6+3*d*e^2*x.^5+(3*e^2+3*d^2*e).*x.^4+(6*d*e+d^3)*x.^3+(3*e+3*d^2).*x.^2+3*d*x+1); end function yhat=ProfileDataModelO2(beta, x) a=beta(1); b=beta(2); c=beta(3); yhat=1./(a + b.*((-x).^c)); end function yhat=ProfileDataModelO2_1stDer(beta, x) % %differential of ProfileDataModel2 % to obtain this type: % syms x a b c; % result = diff(1/(a + b*(x^c))) a=beta(1);

b=beta(2); c=beta(3); yhat=(b.*c.*(-x).^(c - 1))./(a + b.*(-x).^c).^2; end function yhat=ProfileDataModelO2_2ndDer(beta, x)% 2nd differential of ProfileDataModel2 % to obtain this type: % syms x a b c; % result = diff(-(b*c*x^(c - 1))/(a + b*x^c)^2) a=beta(1); b=beta(2); c=beta(3); %yhat=(2.*b.^2.*c.^2.*x.^(2.*c - 2))./(a + b.*x.^c).^3 - ... %(b.*c.*x.^(c - 2).*(c - 1))./(a + b.*x.^c).^2; yhat = (2*b.^2.*c.^2.*(-x).^(2.*c - 2))./(a + b.*(-x).^c).^3 - ... (b.*c.*(-x).^(c - 2).*(c - 1))./(a + b.*(-x).^c).^2; end

GasSol.m function [ConcCH4,ConcN2,ConcO2,ConcAr, ConcCO2] = GasSol(S,TC,pGdry,units) % GasSol Solubility of a gas in sea water %========================================================================= % Recent Additions:2011,Ryan Bell SRI International, [email protected] % % Insipired/modeled from: % Arsol Version 1.X 7/13/2004 % Author: Roberta C. Hamme (Scripps Inst of Oceanography) % % USAGE: [ConcCH4,ConcN2,ConcO2,ConcAr, ConcCO2] = GasSol(S,T,pGdry,Gas,units) % % DESCRIPTION: % Solubility (saturation) of a Gas (Ch4,N2,O2,Ar and CO2) in sea water % % % INPUT: (if S and T are not singular they must have same dimensions) % S = salinity [PSS] % T = temperature [degree C] % pGdry = dry mole of a gas fraction (0.00934 Ar in dry atm air ) % units = 0,1,2,3 for Molinity,Molarity,Molality,Molar Fraction % - implemented but disengaged for the moment due to periodic errors % - wikipedia's conentration article has a good discussion % on these differences % OUTPUT: % [ConcCH4,ConcN2,ConcO2,ConcAr, ConcCO2] = solubility of % CH4,N2,O2,Ar,CO2 in umole/kg % % DISCLAIMER: % This software is provided "as is" without warranty of any kind. %========================================================================= %---------------------% CHECK INPUT ARGUMENTS %---------------------if nargin ~=4 error('GasSol.m: Must pass 4 parameters') else

% Check S,T,pGdry dimensions and verify consistent [ms,ns] = size(S); [mt,nt] = size(TC); [mp,np] = size(pGdry); [mu,nu] = size(units); % Check that T,S&P have the same shape or are singular if ((ms~=mt) || (ns~=nt) || np~=ns || mp~=ms) error('GasSol: S, T, pGdry and units must have same dimensions') end if (mu~=1) || (nu~=1) error('GasSol: units must be singular (0,1,2,or3)') end end if units~=0 error('GasSol: units outside of molinity not yet debugged') end %-----% BEGIN %-----%% Temperature conversions % convert T to scaled temperature temp_S = log((298.15 - TC)./(273.15 + TC)); % convert T to Kelvin T = 273.15 + TC; %% For Methane % REFERENCES: % Wiesenburg, D. A. 1979 % "Equilibrium Solubilties of Methane, Carbon Monoxide, and Hydrogen in Water and Seawater." % Journal of Chemical and Engineering Data, 24(4), 1979,pg 356-360. % constants for Eqn (7) of Wiesenburg 1979 a1 = -417.5053; a2 = 599.8626; a3 = 380.3636; a4 = -62.0764; %This term corrects for the use of dry air b1 = -0.064236; b2 = 0.034980; b3 = -0.0052732; % Eqn (7) of Wiesenburg 1979 lnC = log(pGdry) + a1 + a2.*(100./T) + a3.*log(T./100) + a4.*(T./100) + ... S.*( b1 + b2.*(T./100) + b3.*((T./100).^2) );

c = exp(lnC); %reported in nmol/kg by Wiesenburg ConcCH4=c/1000; %Converted to umol/kg by Ryan %% Do the Calcualtion for Nitrogen % REFERENCE: % Roberta Hamme and Steve Emerson, 2004. % "The solubility of neon, nitrogen and argon in distilled water and seawater." % Deep-Sea Research I, in press. % constants from Table 4 of Hamme and Emerson 2004 A0_n2 = 6.42931; A1_n2 = 2.92704; A2_n2 = 4.32531; A3_n2 = 4.69149; B0_n2 = -7.44129e-3; B1_n2 = -8.02566e-3; B2_n2 = -1.46775e-2; % Eqn (1) of Hamme and Emerson 2004 conc_N2_air = exp(A0_n2 + A1_n2*temp_S + A2_n2*temp_S.^2 + A3_n2*temp_S.^3 +... S.*(B0_n2 + B1_n2*temp_S + B2_n2*temp_S.^2)); %umole/kg of gas from a moist atmosphere ConcN2 = conc_N2_air.*pGdry.*(1-vpress(S,TC))./((0.78084.*(1-vpress(S,TC)))); %% Do the Calcualtion for Oxygen % REFERENCE: % Hernan E. Garcia and Louis I. Gordon, 1992. % "Oxygen solubility in seawater: Better fitting equations" % Limnology and Oceanography, 37, pp. 1307-1312. % constants from Table 1 of Garcia & Gordon for the fit to Benson and Krause (1984) A0_o2 = 5.80871; A1_o2 = 3.20291; A2_o2 = 4.17887; A3_o2 = 5.10006; A4_o2 = -9.86643e-2; A5_o2 = 3.80369; B0_o2 = -7.01577e-3; B1_o2 = -7.70028e-3; B2_o2 = -1.13864e-2; B3_o2 = -9.51519e-3; C0_o2 = -2.75915e-7; % Corrected Eqn (8) of Garcia and Gordon 1992 conc_O2_air = exp(A0_o2 + A1_o2*temp_S + A2_o2*temp_S.^2 + A3_o2*temp_S.^3 +...

A4_o2*temp_S.^4 + A5_o2*temp_S.^5 +... S.*(B0_o2 + B1_o2*temp_S + B2_o2*temp_S.^2 + B3_o2*temp_S.^3) + C0_o2*S.^2); ConcO2 = conc_O2_air.*pGdry.*(1-vpress(S,TC))./(0.20946.*(1-vpress(S,TC)));

%% Do the Calcualtion for Argon % REFERENCE: % Roberta Hamme and Steve Emerson, 2004. % "The solubility of neon, nitrogen and argon in distilled water and seawater." % Deep-Sea Research I, in press. % constants from Table 4 of Hamme and Emerson 2004 A0_ar = 2.79150; A1_ar = 3.17609; A2_ar = 4.13116; A3_ar = 4.90379; B0_ar = -6.96233e-3; B1_ar = -7.66670e-3; B2_ar = -1.16888e-2; % Eqn (1) of Hamme and Emerson 2004 conc_Ar_air = exp(A0_ar + A1_ar*temp_S + A2_ar*temp_S.^2 + A3_ar*temp_S.^3+... S.*(B0_ar + B1_ar*temp_S + B2_ar*temp_S.^2)); ConcAr = conc_Ar_air.*pGdry.*(1-vpress(S,TC))./(0.00934.*(1-vpress(S,TC))); %% Do the Calcualtion for Carbon Dioxide %From Weiss % a1 = -60.2409; % a2 = 93.4517; % a3 = 23.3585; % b1 = 0.023517; % b2 = -0.023656; % b3 = 0.0047036; % % % Eqn (12) of Weiss 1974 % lnK = a1 + a2.*(100./T) + a3.*log(T./100) + ... % S.*( b1 + b2.*(T./100) + b3.*((T./100).^2) ); % % K = exp(lnK); % %convert dry pCO2 to wet pCO2 % pCO2wet = pGdry.*(1-vpress(S,TC)); % %convert to fugacify % fCO2=pCO2tofCO2(pCO2wet,TC); % %fCO2=pCO2;

% ConcCO2 = K.*fCO2.*1000000; %

%in umol/kg-soln

%from http://cdiac.ornl.gov/oceans/co2rprt.html pCO2wet = pGdry.*(1-vpress(S,TC)); pCO2wet=pCO2wet*1000000; %to uatm [Data,~,~]=CO2SYS(pCO2wet,2000,4,1,S,TC,TC,0,0,0,0,1,4,1); ConcCO2 = Data(:,8) ;

%Calculate the [CO2]

%% Some code I'm working on, it works, but has an error the pops up periodically % % switch units % case 0 % ConcAr = Molinity; % case 1 % [Molarity Molality MolarFraction] = sw_ConcSolver(S,TC,Molinity); % ConcAr = Molarity; % case 2 % [Molarity Molality MolarFraction] = sw_ConcSolver(S,TC,Molinity); % ConcAr = Molality; % case 3 % [Molarity Molality MolarFraction] = sw_ConcSolver(S,TC,Molinity); % ConcAr = MolarFraction; % end end %% the vpress function from Hamme function [vapor_press] = vpress(S,T) % vpress Vapor pressure of sea water %========================================================================= % vpress Version 1.0 8/30/2004 % Author: Roberta C. Hamme (Scripps Inst of Oceanography) % % USAGE: vapor_press = vpress(S,T) % % DESCRIPTION: % Vapor press of sea water % % INPUT: (if S and T are not singular they must have same dimensions) % S = salinity [PSS] % T = temperature [degree C] % % OUTPUT:

% vapor_press = vapor pressure of seawater [atm] % % AUTHOR: Roberta Hamme ([email protected]) % % REFERENCE: % Vapor pressure of pure water: D. Ambrose and I.J. Lawrenson % "The vapour pressure of water" % Journal of Chemical Thermodynamics, v.4, p. 755-671. % Correction for seawater: Frank J. Millero and Wing H. Leung % "The thermodynamics of seawater at one atmosphere" % American Journal of Science, v. 276, p. 1035-1077.

%---------------------% Check input parameters %---------------------if nargin ~=2 error('vpress.m: Must pass 2 parameters') end %if % CHECK S,T dimensions and verify consistent [ms,ns] = size(S); [mt,nt] = size(T); % Check that T&S have the same shape or are singular if ((ms~=mt) | (ns~=nt)) & (ms+ns>2) & (mt+nt>2) error('vpress: S & T must have same dimensions or be singular') end %if %-----% BEGIN %-----%Calculate temperature in Kelvin and modified temperature for Chebyshev polynomial temp_K = T+273.15; temp_mod = (2*temp_K-(648+273))/(648-273); %Calculate value of Chebyshev polynomial Chebyshev = (2794.0144/2)+(1430.6181*(temp_mod))+... (-18.2465*(2*temp_mod.^2-1))+... (7.6875*(4*temp_mod.^3-3*temp_mod))+... (-0.0328*(8*temp_mod.^4-8*temp_mod.^2+1))+... (0.2728*(16*temp_mod.^5-20*temp_mod.^3+5*temp_mod))+... (0.1371*(32*temp_mod.^6-48*temp_mod.^4+18*temp_mod.^2-1))+... (0.0629*(64*temp_mod.^7-112*temp_mod.^5+56*temp_mod.^3-7*temp_mod))+... (0.0261*(128*temp_mod.^8-256*temp_mod.^6+160*temp_mod.^4-32*temp_mod.^2+1))+...

(0.02*(256*temp_mod.^9-576*temp_mod.^7+432*temp_mod.^5-120*temp_mod.^3+9*temp_mod))+... (0.0117*(512*temp_mod.^10-1280*temp_mod.^8+1120*temp_mod.^6-400*temp_mod.^4+50*temp_mod.^2-1))+... (0.0067*(1024*temp_mod.^11-2816*temp_mod.^9+2816*temp_mod.^7-1232*temp_mod.^5+220*temp_mod.^3-11*temp_mod)); %Vapor pressure of pure water in kiloPascals and mm of Hg vapor_0sal_kPa = 10.^(Chebyshev./temp_K); vapor_0sal_mmHg = vapor_0sal_kPa*1000*0.00750062; %Correct vapor pressure for salinity vapor_press =(vapor_0sal_mmHg+(S.*(-0.0023311+(-0.00014799*T)+(-0.00000752*T.^2)+... (-0.000000055185*T.^3)))+S.^1.5.*(-0.00001132+(-0.0000087086*T)+... (0.00000074936*T.^2)+(-0.000000026327*T.^3)))*0.001315789; end

InterpTthenD.m

function [C T D] = InterpDthenT(Time, Depth, Data, Bin, tsteps, dsteps, method1, method2) % InterpDthenT interpolates data for the contour plots % % % %

Ryan Bell © SRI International, 2011 Currently at Vancouver Island University permanent email: [email protected] Works in R2010b

% USAGE: [C T D] = InterpDthenT(Time, Depth, Data, Bin, tsteps, dsteps, method1, method2) % % DESCRIPTION: % Data in each Bin is interpolated from it associated Depth to new % depths (n= dsteps) across the same range using method1. Data is then % interpolated acrosss Time to new time (where n=tsteps) using method2. % Error checking has not been implented % % INPUT: % Time - scalar array of times associtaed with Data % Depth - scalar array of depths associtaed with Data % Data - scalar array of concentrations % Bin - scalar array of integers indicating the bin(depth profile) % associated concentrations belong to. % tsteps - the number of interpolated time intervals % dsteps - the number of interpolated depth intervals % method1 - the interpolation method used for the interpolation along depth % method2 - the interpolation method used for the interpolation along time % OUTPUT: % [C,T,D] = 2D matrices of Concentration, time and depth % % DISCLAIMER: % This software is provided "as is" without warranty of any kind. %=========================================================================

%% Interpolate across all points of equal bin di=linspace(min(Depth), max(Depth), dsteps);

%Sort on Bin (ie all the points of equal bin (ie 1 probe into sediment) AllData = horzcat(Bin,Depth,Time,Data); AllData =sortrows(AllData,1); bj = AllData(:,1); dj = AllData(:,2); %depth sorted depth array tj = AllData(:,3); %depth sorted time array cj = AllData(:,4); %depth sorted conc array n=1; m=1; %m will track the colomns j=1; %j will track start position, for i=1:size(bj,1) %i will track end position of each bin if i==size(bj,1) %if we're at the end Cj(:,m) = interp1(dj(j:i),cj(j:i),di,method1,NaN); warning off last Tj(:,m) = (mean(tj(j:i))).*ones(size(di)); Dj(:,m) = di; m=m+1;%count the columns elseif bj(i)~=bj(i+1) %if the next bin value is different, interpoate values for all values from i up to j Cj(:,m) = interp1(dj(j:i),cj(j:i),di,method1,NaN); %interpolate accross all the data at this depth warning off last Tj(:,m) = (mean(tj(j:i))).*ones(size(di)); Dj(:,m) = di; m=m+1;%count the columns j=i+1; else n=n+1;

%count the rows

end end

%% Interpolate across all points of equal depth

ti=linspace(min(Time), max(Time), tsteps); Ck=ones(size(Cj,1),numel(ti)); %preallocate, this makes a big differnce in theis case Dk=ones(size(Cj,1),numel(ti)); %preallocate Tk=ones(size(Cj,1),numel(ti)); %preallocate for i=1:size(Cj,1)

if sum(isfinite(Cj(i,:)))>1 Ck(i,:)= interp1(Tj(i,:),Cj(i,:),ti,method2,'extrap'); else %other methods get tripped up on sparse data Ck(i,:)= interp1(Tj(i,:),Cj(i,:),ti,'linear','extrap'); end warning off last Dk(i,:) = Dj(i,1)*ones(size(ti)); Tk(i,:) = ti; end warning on last C=Ck; T=Tk; D=Dk;

% % % % %

figure(11);clf set(gcf,'name','test contour','position',[20 20 400 200]); [p1 h1]=contourf(gca,T',D',C'); set(gca,'ydir','reverse') shg

PlotContour.m function PlotContour(Conc,T,D,t1,d1,t2,d2,PlotFig) % PlotContour is a custom plotting function for Depth-Time contours % Ryan Bell © SRI International, 2011 % Currently at Vancouver Island University % permanent email: [email protected] % Works in R2010b % % DISCLAIMER: % This software is provided "as is" without warranty of any kind. %========================================================================= figure(PlotFig) %% D=-D; d1=-d1; d2=-d2; Dmax = ceil(max(max(D))); Dmin = floor(min(min(D))); FontMag=10; Fontweigh ='bold'; fontname = 'Times New Roman'; xtag =0.01; ytag = .82; colortag='k'; ytickpos = 3*floor(Dmin/3):3:3*ceil(Dmax/3); ylim = [Dmin Dmax]; Tags = {'a) O_{2}','b) CO_{2}','c) CH_{4}',' d) N_{2}/Ar'}; %% Define x axis and Generate X axis labels for hours on contour %[yr mon day hr mn sec]=datevec(T); timeint = 1/4; %time intervals for x labels in days timemax=timeint*(floor(max(max(T))/timeint)); timemin=timeint*(ceil(min(min(T))/timeint)); xlim=[min(min(T)) max(max(T))]; xtick = timemin:timeint:timemax; xticklabel = datestr(xtick, 15);

%% Lets start from the bottom, C=Conc.N2_2./Conc.Ar;

- plot n2/ar

res =256; %define the color resolution colormap(jet(res)); Cmax =max(max(C)); %I'd like to put some code here, but its a little complex with the colorbar Cmin =min(min(C)); colormap_index = fix(((C-Cmin)/(Cmax-Cmin))*(res-1))+1; %this is the y data rounded to value capable of the resolution yplotted = ((colormap_index-1)*(Cmax-Cmin))/(res-1) + Cmin;

%% Set up axes xoff = 0.1; offset = .08; Hcorr = 1.15; yoff = offset+(Hcorr*0)/4; w = .86; h = 0.2; ah1=axes('position',[xoff yoff w h]); %integers between 1 and res represent color, from the colormap [~, h1]=contourf(gca,T',D',colormap_index',res); set(h1,'linestyle','none'); set(gca,'ydir','normal','linewidth',1.5,'ylim',ylim,'ytick',ytickpos); set(gca,'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname,'xlim',xlim,'xtick',xtick); set (gca, 'Box', 'off', 'YAxisLocation', 'left','xticklabel', xticklabel, 'tickdir','out'); xlabel('Time (EDT)'); l1=ylabel('Depth (cm)'); %Title(CTitle); hold on plot(t1,d1,'+k','linewidth',1.5,'markersize',2); plot(t2,d2,'-k','linewidth',1.0); %plot(CarboTime, Carbocline,'-k','linewidth',1.5); %text(.015,.9,'e) CO_2 \mumole/kg','color','w','units','normalized','fontsize',FontMag,'fontweight',Fontweigh); text(xtag,ytag,Tags{4},'color',colortag,'units','normalized','fontsize',... FontMag,'fontweight',Fontweigh,'fontname',fontname); %% Put in some Dates too timeint=1; timemax=timeint*(floor(max(max(T))/timeint)); timemin=timeint*(ceil(min(min(T))/timeint));

%Time interval of label - 1 equals once a day %max label %min label

xtickb = timemin:timeint:timemax ; xticklabelb = datestr(xtickb, 2) ;

%create the tick x location %create the labels

for n=1:size(xticklabelb,1) text(xtickb(n),Dmax+5 ,datestr(xtickb(n), 2),'units','data','fontsize',... FontMag,'fontweight',Fontweigh,'HorizontalAlignment','center','fontname',fontname) end %%

Plot Color bar

cb=colorbar; %title(cb,'CO_2 \mumole/kg','fontweight',Fontweigh,'fontname',fontname,'fontsize',FontMag); set(cb,'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname,'linewidth',1.5); %How many markers on the colorbar? ticknum = 5; tickspot = linspace(1, res-1, ticknum); %Space out the markers linearly between 1 and res-1 %relabel them to match the actual concentrations ticklabels =linspace(min(min(yplotted)), max(max(yplotted)), ticknum);

%handle significatn figures if max(ticklabels)>1000 ticklabels2 = (round(ticklabels*.1))/.1; elseif max(ticklabels)>=100 ticklabels2 = (round(ticklabels*1))/1; elseif max(ticklabels)>=10 ticklabels2 = (round(ticklabels*10))/10; elseif max(ticklabels)>=1 ticklabels2 = (round(ticklabels*100))/100; else ticklabels2 = (round(ticklabels*1000))/1000; end set(cb,'YTick',tickspot); set(cb,'YTickLabel',ticklabels2);

%% %create a box,(it had to be removed to handle left/right tick issues ah0 = axes ('Position', get(ah1,'position')); set (ah0, 'Box', 'on', 'Color', 'none', 'XTick', [], 'YTick', [],'linewidth',1.5);

% % % % % % % % % % % % % % % % % % % % % % % %

% % % % % % % % % % % % %%%%%%%%%%%%%%%%Now on to Metrhane%%%%%%%%%%%%%%%%%%%%%% % % % % % % % % % % % % %

C=Conc.CH4; res =256; %define the color resolution colormap(jet(res)); Cmax =max(max(C)); Cmin =min(min(C)); %I'd like to put some code here, but its a little complex with the colorbar colormap_index = fix(((C-Cmin)/(Cmax-Cmin))*(res-1))+1; %this is the y data rounded to value capable of the resolution yplotted = ((colormap_index-1)*(Cmax-Cmin))/(res-1) + Cmin;

%% Set up axes yoff = offset+(Hcorr*1)/5; ah1=axes('position',[xoff yoff w h]); [~, h1]=contourf(gca,T',D',colormap_index',res); %integers between 1 and res represent color, from the colormap set(h1,'linestyle','none'); set(gca,'ydir','normal','linewidth',1.5,'ylim',ylim,'ytick',ytickpos); set(gca,'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname) set(gca,'xlim',xlim,'xtick',xtick,'xticklabel', []); set (gca, 'Box', 'off', 'YAxisLocation', 'left', 'tickdir','out'); %xlabel('Time (EDT)'); l2=ylabel('Depth (cm)'); %Title(CTitle); hold on plot(t1,d1,'+k','linewidth',1.5,'markersize',2); plot(t2,d2,'-k','linewidth',1.0); %text(.015,.9,'d) Oxygen \mumole/kg','color',colortag,'units','units','normalized',... 'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname); text(xtag,ytag,Tags{3},'color',colortag,'units','normalized',... 'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname); %%

Plot Color bar

cb=colorbar; %title(cb,'N_2 \mumole/kg','fontweight',Fontweigh,'fontname',fontname,'fontsize',FontMag); set(cb,'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname,'linewidth',1.5); %How many markers on the colorbar? ticknum = 5;

tickspot = linspace(1, res-1, ticknum);%Space out the markers linearly between 1 and res-1 %relabel them to match the actual concentrations ticklabels =linspace(min(min(yplotted)), max(max(yplotted)), ticknum); if max(ticklabels)>1000 ticklabels2 = (round(ticklabels*.1))/.1; elseif max(ticklabels)>=100 ticklabels2 = (round(ticklabels*1))/1; elseif max(ticklabels)>=10 ticklabels2 = (round(ticklabels*10))/10; elseif max(ticklabels)>=1 ticklabels2 = (round(ticklabels*100))/100; else ticklabels2 = (round(ticklabels*1000))/1000; end set(cb,'YTick',tickspot); set(cb,'YTickLabel',ticklabels2);

%% %create a box,(it had to be removed to handle left/right tick issues ah0 = axes ('Position', get(ah1,'position')); set (ah0, 'Box', 'on', 'Color', 'none', 'XTick', [], 'YTick', [],'linewidth',1.5);

% % % % % % % % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %%%%%%%%%%%%%%%%Now on to CO2%%%%%%%%%%%%%%%%%%%%%% % %

C=Conc.CO2; res =256; %define the color resolution colormap(jet(res)); Cmax =max(max(C)); Cmin =min(min(C)); colormap_index = fix(((C-Cmin)/(Cmax-Cmin))*(res-1))+1; %this is the y data rounded to value capable of the resolution yplotted = ((colormap_index-1)*(Cmax-Cmin))/(res-1) + Cmin;

%% Set up axes yoff = offset+(Hcorr*2)/5; ah1=axes('position',[xoff yoff w h]); %integers between 1 and res represent color, from the colormap [~, h1]=contourf(gca,T',D',colormap_index',res); set(h1,'linestyle','none'); set(gca,'ydir','normal','linewidth',1.5,'ylim',ylim,'ytick',ytickpos); set(gca,'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname); set(gca,'xlim',xlim,'xtick',xtick,'xticklabel', []); set (gca, 'Box', 'off', 'YAxisLocation', 'left', 'tickdir','out'); %xlabel('Time (EDT)'); l3 = ylabel('Depth (cm)'); %Title(CTitle); hold on plot(t1,d1,'+k','linewidth',1.5,'markersize',2); plot(t2,d2,'-k','linewidth',1.0); %text(.015,.9,'c) O_2 \mumole/kg','color',colortag,'units','units','normalized',... 'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname); text(xtag,ytag,Tags{2},'color','w','units','normalized','fontsize',FontMag,... 'fontweight',Fontweigh,'fontname',fontname); %%

Plot Color bar

cb=colorbar; %title(cb,'N2 \mumole/kg','fontweight',Fontweigh,'fontname',fontname,'fontsize',FontMag); set(cb,'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname,'linewidth',1.5); ticknum = 5; %How many markers on the colorbar? tickspot = linspace(1, res-1, ticknum); %Space out the markers linearly between 1 and res-1 %relabel them to match the actual concentrations ticklabels =linspace(min(min(yplotted)), max(max(yplotted)), ticknum); if max(ticklabels)>1000 ticklabels2 = (round(ticklabels*.1))/.1; elseif max(ticklabels)>=100 ticklabels2 = (round(ticklabels*1))/1; elseif max(ticklabels)>=10 ticklabels2 = (round(ticklabels*10))/10; elseif max(ticklabels)>=1 ticklabels2 = (round(ticklabels*100))/100; else ticklabels2 = (round(ticklabels*1000))/1000; end set(cb,'YTick',tickspot);

set(cb,'YTickLabel',ticklabels2);

%% %create a box,(it had to be removed to handle left/right tick issues ah0 = axes ('Position', get(ah1,'position')); set (ah0, 'Box', 'on', 'Color', 'none', 'XTick', [], 'YTick', [],'linewidth',1.5);

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

% % % %

Conc.O2(end,end)=0;

% % % %

% % % %%%%%%%%%%%%%%%%Now on to Oxygen %%%%%%%%%%%%%%%%%%%%%% % % %force one of the point to zero for scaling purposese(cheat a little)

C=Conc.O2; res =256; %define the color resolution colormap(jet(res)); Cmax =max(max(C)); Cmin = min(min(C)); %I'd like to put some code here, but its a little complex with the colorbar colormap_index = fix(((C-Cmin)/(Cmax-Cmin))*(res-1))+1; %this is the y data rounded to value capable of the resolution yplotted = ((colormap_index-1)*(Cmax-Cmin))/(res-1) + Cmin;

%% Set up axes yoff = offset+(Hcorr*3)/5; ah1=axes('position',[xoff yoff w h]); %integers between 1 and res represent color, from the colormap [~, h1]=contourf(gca,T',D',colormap_index',res); set(h1,'linestyle','none'); set(gca,'ydir','normal','linewidth',1.5,'ylim',ylim,'ytick',ytickpos); set(gca,'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname) set(gca,'xlim',xlim,'xtick',xtick,'xticklabel', []); set (gca, 'Box', 'off', 'YAxisLocation', 'left', 'tickdir','out'); %xlabel('Time (EDT)'); l4 = ylabel('Depth (cm)'); %Title(CTitle); hold on plot(t1,d1,'+k','linewidth',1.5,'markersize',2);

plot(t2,d2,'-k','linewidth',1.0); %text(.015,.9,'b) CH_4 \mumole/kg','color',colortag,'units','units','normalized',... 'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname); text(xtag,ytag,Tags{1},'color','w','units','normalized','fontsize',FontMag,... 'fontweight',Fontweigh,'fontname',fontname); %%

Plot Color bar

cb=colorbar; %title(cb,'CH_4 \mumole/kg','fontweight',Fontweigh,'fontname',fontname,'fontsize',FontMag); set(cb,'fontsize',FontMag,'fontweight',Fontweigh,'fontname',fontname,'linewidth',1.5); ticknum = 5; %How many markers on the colorbar? tickspot = linspace(1, res-1, ticknum); %Space out the markers linearly between 1 and res-1 %relabel them to match the actual concentrations ticklabels =linspace(min(min(yplotted)), max(max(yplotted)), ticknum); %ticklabels =linspace(0, max(max(yplotted)), ticknum); %relabel them to match the actual concentrations

if max(ticklabels)>1000 ticklabels2 = (round(ticklabels*.1))/.1; elseif max(ticklabels)>=100 ticklabels2 = (round(ticklabels*1))/1; elseif max(ticklabels)>=10 ticklabels2 = (round(ticklabels*10))/10; elseif max(ticklabels)>=1 ticklabels2 = (round(ticklabels*100))/100; else ticklabels2 = (round(ticklabels*1000))/1000; end set(cb,'YTick',tickspot); set(cb,'YTickLabel',ticklabels2);

%% %create a box,(it had to be removed to handle left/right tick issues ah0 = axes ('Position', get(ah1,'position')); set (ah0, 'Box', 'on', 'Color', 'none', 'XTick', [], 'YTick', [],'linewidth',1.5); % % % % % % % % % % %%%%%%%%%%%%%And last and defintley least:%%%%%%%%%%%%%%55 %% set label positions p1 = get(l1,'position');p2 = get(l2,'position');p3 = get(l3,'position'); p4 = get(l4,'position');%p5 = get(l5,'position'); ylabelypos = min([p1(1) p2(1) p3(1) p4(1)]);

set(l1,'position',[ylabelypos p1(2) p1(3)]); set(l2,'position',[ylabelypos p2(2) p2(3)]); set(l3,'position',[ylabelypos p3(2) p3(3)]); set(l4,'position',[ylabelypos p4(2) p4(3)]); %set(l5,'position',[ylabelypos p5(2) p5(3)]); %% Put on icons % [X2,map2] = imread('sri_logo_1_blu_000.gif'); % ah4=axes('position',[.91 0.005 .08 .08]); % axes(ah4); % subimage(X2,map2); %need to use subimage to avoid colormap issues % %%%%imshow('sri_logo_1_blu.jpg') % set(ah4,'xtick',[],'ytick',[]);