% Adapt Johnson 1D dissolved Fe model Excel to MATLAB, % based on Johnson et al.,1999 Marine Chemistry % % Adapted by Mak Saito, November 26, 2005 clear % Set Metal Constants and Boundary conditions soly=0.6 ; %nM Ligand concentration FeToC=5; %umol/mol Metal to carbon Ratio in sinking material KmO2=0; % Half saturation constant for O2 dependent scavenging k=0.005; %/yr Scavenging rate constant Fei= 0.0500; % initial metal condition for profile BoundaryMetalConc=0.6; % Boundary condition % % This version examines variation in the Flux100 variable and is modified % to cobalt concentrations % Adapt April, 25, 2009 - add in MnO2 scavenging model to this model % Need to parse the k prime by depth steps and assign them to this model, % will be step like though % April 30, 2009, O2 scavenging added, but only a % May 1, 2009 - add OH and revise constants. % Adapt September 27, 2009 % Renamed Hybrid1.m - make universal for Fe, Mn, and Co. % Change O2 sensitivity to Michealis Menten formulation instead of third % order kinetic parameter % Input Oxygen datasets for three stations % Process O2 data S5 KM0405; %load stn5O2UP.mat; load stn502.dat; load P16O2.dat; % from 18.3998 -154.4737, depth, O2 stn502=flipud(stn502); data=ones(2,2291)'; data2=ones(2,80)'; for i= 1:2291; if stn502(i,1) <= 4001; data(i,1)=stn502(i,1); data(i,2)=stn502(i,2); else data(i)= NaN; end end for j=1:80; for i= 1:2291; data2(j,1)=j*50; if data(i,1) <=(j)*50; data2(j,2)=data(i,2); end end end S5O2=data2(:,2); data=ones(2,2291)'; data2=ones(2,80)'; for i= 1:2291; if stn502(i,1) <= 4001; data(i,1)=stn502(i,1); data(i,2)=stn502(i,3); else data(i)= NaN; end end for j=1:80; for i= 1:2291; data2(j,1)=j*50; if data(i,1) <=(j)*50; data2(j,2)=data(i,2); end end end S5Temp=data2(:,2); % this could be cleaned up combining with stn5O2 data process in three column matrix % Process WOCE P16 O2 data (Near Station 9 KM0405) data=ones(2,36)'; data2=ones(2,80)'; for i= 1:36; if P16O2(i,1) <= 4001; data(i,1)=P16O2(i,1); data(i,2)=P16O2(i,2); else data(i)= NaN; end end for j=1:80; for i= 1:36; data2(j,1)=j*50; if data(i,1) <=(j)*50; data2(j,2)=data(i,2); end end end P16O2=data2(:,2); % Process O2 data from Vertex T-5 % Martin 1989 Vertex Station T-5 39.60N 140.77W % Depth Oxygen Mn Co Fe % Note The data from Martin and the data in Johnson's figure are quite % different. The upper water column has about 195 O2 not 247. % O2 and pH data are from Martin 1985 even in Johnson's paper, whereas % metal data are from Martin 1989 T5=[100 247 0.43 13.4 0.06 150 217 0.3 32 nan 200 211 0.19 38.1 0.2 250 204 0.23 36 nan 300 181 0.21 36 0.27 400 143 0.2 32.6 0.4 500 100 0.36 36.8 0.46 600 65 0.47 37.4 0.68 690 40 0.54 49.4 nan 780 24 0.53 32.5 0.72 900 18 0.44 29.8 0.71 1000 16 0.47 31.2 0.76 1240 17 0.38 23.8 0.66 1500 28 0.26 22.4 0.66 1750 50 nan nan nan 2000 75 nan nan nan]; % add 1750m and 2000m estimated values from Johnson 1996 figure, change 100 to 101 data=ones(2,16)'; data2=ones(2,80)'; for i= 1:16; if T5(i,1) <= 4001; data(i,1)=T5(i,1); data(i,2)=T5(i,2); else data(i)= NaN; end end for j=1:80; for i= 1:16; data2(j,1)=j*50; if data(i,1) <=(j)*50; data2(j,2)=data(i,2); end end end data2(1,2)=data2(2,2); % copy 100m O2 to 50m O2, check T5O2=data2(:,2); % Cobalt model - % Figure 1 Cobalt scenarios % Metal constants now at beginning of script % Set non-metal Constants B=0.858; %/m Depth exponent in carbon flux eqn Kz=3153.6; %m2/yr Vertical eddy diffusion coefficient Fluxcon=223.0782825; %Constant part of carbon flux eqn FluxconV=FeToC*B/(100^-B); %Allows This constant to peak dependent of Fe quota value delt=0.2; %yr Time step delz=50; %m Depth step reset=1; %if = 0, reinitialize all iron depths to 0.05nM iterations = 5000 ; % Set Non-metal Boundary conditions MaxD = 4000; % NumDepSteps=(MaxD)/delz; % Note the model cannot start above 100m or it will be unstable % Initialize arrays z=((1:NumDepSteps).*delz)'; Fet=(ones(1,NumDepSteps).*Fei)'; Diff=(zeros(1,NumDepSteps))'; Prod=(zeros(1,NumDepSteps))'; Scav=(zeros(1,NumDepSteps))'; Fetplus1=(ones(1,NumDepSteps).*Fei)'; Fetplus1(80)=BoundaryMetalConc; % Run Model and Plot figure(1) clf for a=1:4; for f=1:10 Flux100=f; %mol/m2/yr Export flux of carbon at 100 m for r=1:5000 % iterations Fet=Fetplus1; for i=2:(NumDepSteps-1) Diff(i)=Kz*(Fet(i+1)-2*Fet(i)+Fet(i-1))/delz^2; % note there is a if statement not included here yet Prod(i)=FluxconV*Flux100*z(i)^(-(1+B)); if Fet(i)soly; if a==1 Scav(i)=k*(soly-Fet(i))*(P16O2(i)/(KmO2+P16O2(i))); % new form 9.27.09 monod O2 form end if a==2 Scav(i)=k*(soly-Fet(i))*(S5O2(i)/(KmO2+S5O2(i))); % new form 9.27.09 monod O2 form end if a==3 Scav(i)=k*(soly-Fet(i))*(T5O2(i)/(KmO2+T5O2(i))); % new form 9.27.09 monod O2 form end if a==4 % k=0.052; %/yr Scavenging rate constant % Scav(i)=k*(soly-Fet(i)); %no O2 effect, original Johnson form end end Fetplus1(i)=Fet(i)+delt*(Diff(i)+Prod(i)+Scav(i)); end end if a==1 subplot(3,2,1); plot(Fet,-z); title('Stn 9 Metal'); xlabel('nM');ylabel('Depth (m)'); hold on subplot(3,2,2); plot(Scav,-z); title('Stn 9 Scavenging term') ; xlabel('');ylabel('Depth (m)');hold on end if a==2 subplot(3,2,3); plot(Fet,-z); title('Stn 5 Metal'); xlabel('nM');ylabel('Depth (m)'); hold on subplot(3,2,4); plot(Scav,-z); title('Stn 5 Scavenging term') ; xlabel('');ylabel('Depth (m)');hold on end if a==3 subplot(3,2,5); plot(Fet,-z); title('Stn T5 Metal'); xlabel('nM');ylabel('Depth (m)'); hold on subplot(3,2,6); plot(Scav,-z); title('Stn T5 Scavenging term') ; xlabel('');ylabel('Depth (m)');hold on end end end figure(2) subplot(3,1,1); plot(P16O2, -z); title('Stn 9 O2'); ylabel('Depth (m)'); xlabel('uM'); axis([ 0 300 -4000 0]); subplot(3,1,2); plot(S5O2, -z); title('Stn 5 O2'); ylabel('Depth (m)'); xlabel('uM');axis([ 0 300 -4000 0]); subplot(3,1,3); plot(T5O2, -z); title('Vertex T5 O2'); ylabel('Depth (m)'); xlabel('uM');axis([ 0 300 -4000 0]);