function doall=box_geochem % A 3 box THC model with geochemistry. Written as a practice problem % to study basic geochemistry. See documantation in enclosed latex % notes, file paleo_co2_tutorial.tex. % Eli Tziperman clear all; close all % Solution of the full model: global Tplot1 Splot1 uplot1 % profile on % Plotted circulation corresponds to a positive THC: % % -------------------------------------- % | | % | atmosphere | % | | % -------------------------------------- % -------------------------------------- % | | | % | 1 ======> 2 | % | ^ | | | % | | |------------|-----------| % |-----|------ | | % | | v | % | | % | | % | | % | 3 | % | | % -------------------------------------- % S. Pole Equator % % set model parameters and initial conditions: [L,H,W,V,Area,Sv,year,u0,alpha,beta,gamma_T,gamma_S ... ,k_h,k_v_highlat,k_v_lowlat,PV ... ,FW,Tstar,Sstar,delta_t,N_max,ubar,T,S,ALK,TCO2,PO4,Hp ... ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4 ... ,solve_physics,solve_geochemistry] ... =set_parameters; % run the model to a steady state: % run the box model: [var_T,var_S,var_THC,Tbar,Sbar,ubar... ,ALKbar,PO4bar,TCO2bar,CO2_atm_bar ... ,Tplot1,Splot1,uplot1]= ... box_model(T,S,ALK,TCO2,PO4,Hp,u0,alpha,beta ... ,k_h,k_v_highlat,k_v_lowlat,PV ... ,gamma_T,gamma_S,year,FW,delta_t,N_max,Tstar,Sstar ... ,V,H,L,Area,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4 ... ,solve_physics,solve_geochemistry); %print -depsc figure1.eps disp(' Steady state solution is:') Tbar Sbar ubar ALKbar PO4bar TCO2bar CO2_atm_bar % profile report %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [var_T,var_S,var_THC,Tbar,Sbar,ubar... ,ALKbar,PO4bar,TCO2bar,CO2_atm_bar ... ,Tplot,Splot,uplot] ... =box_model (T,S,ALK,TCO2,PO4,Hp,u0,alpha,beta ... ,k_h,k_v_highlat,k_v_lowlat,PV ... ,gamma_T,gamma_S,year,FW,delta_t,N_max,Tstar,Sstar ... ,V,H,L,Area,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4 ... ,solve_physics,solve_geochemistry) % Initialize: ALKnew=ALK; ALKnow=ALK; ALKold=ALK; PO4new=PO4; PO4now=PO4; PO4old=PO4; Snew=S; Snow=S; Sold=S; TCO2new=TCO2; TCO2now=TCO2; TCO2old=TCO2; Tnew=T; Tnow=T; Told=T; % Initialize arrays containing plot data: Tplot =zeros(3,N_max); Splot =zeros(3,N_max); ALKplot =zeros(3,N_max); PO4plot =zeros(3,N_max); TCO2plot=zeros(3,N_max); uplot =zeros(1,N_max); % some silly consmetics stuff for indicating progress of run: status_bar='[ ]'; i_status_bar=1; fprintf(1,'%s',['running the box model to steady state: ',status_bar]) % main time loop: for n=1:N_max % some silly consmetics stuff for indicating progress of run: N_print=floor(N_max/20); if n == floor(n/N_print)*N_print i_status_bar=i_status_bar+1; for ii=1:i_status_bar status_bar(ii+1)='.'; end % Delete old status bar and print updated one: fprintf(1,'\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%s',status_bar) end if n == N_max disp(' ') end % Calculate densities: % rho=linear_density(alpha,beta,Tnow,Snow); pressure=1.0; for i=1:3 rho(i)=density(Tnow(i),Snow(i),pressure)/1000.0-1.0; end % Calculate transport: u=transport(u0,rho); if (solve_geochemistry) % Solve for carbonate system: Hp_init_guess=Hp; [CO2_g,H2CO3,OHm,Hp,pH,HCO3m,CO3m2,B_OH_3,B_OH_4m,B_T]= ... solve_carbonate_system(ALKnow,TCO2now,Tnow,Snow,Hp_init_guess); if n == N_max % print details of carbonate system on final time step: disp('The solution of carbonate system on final time step is') ALKnow,TCO2now,Tnow,Snow,Hp_init_guess [CO2_g,H2CO3,OHm,Hp,pH,HCO3m,CO3m2,B_OH_3,B_OH_4m,B_T]= ... solve_carbonate_system(ALKnow,TCO2now,Tnow,Snow,Hp_init_guess) disp('End of carbonate system solution.') disp('export production at end of run:') EP1=export_production(PO4now(1),hPO4,light_factor(1),rPO4(1)); EP2=export_production(PO4now(2),hPO4,light_factor(2),rPO4(2)); % our units are PO4/cm^2/sec. Need to convert to gr C/year and to m^2: export_production_gr_C_per_m2=... [EP1*(122+Rain(Tnow(1))),EP2*(122+Rain(Tnow(2)))] ... *(12.011*10^-6)*(365*86400)*100 global_export_production_Gt_C_per_yr=... [EP1*(122+Rain(Tnow(1)))*Area(1),EP2*(122+Rain(Tnow(2)))*Area(2)] ... *(12.011*10^-6)*(365*86400)*1.e-15 end % solve for atmospheric co2 concentration: CO2_atm=(Area(1)*CO2_g(1)+Area(2)*CO2_g(2))/(Area(1)+Area(2)); end % Leap frog time stepping: if (solve_physics) Tnew=Told+2*delta_t*time_rate_of_change_T ... (n,Tnow,u,Tstar,gamma_T,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat); Snew=Sold+2*delta_t*time_rate_of_change_S ... (n,Snow,u,Sstar,gamma_S,FW,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat); end if (solve_geochemistry) TCO2new=TCO2old+2*delta_t*time_rate_of_change_TCO2 ... (TCO2now,PO4now,Tnow,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat,PV,CO2_atm,CO2_g ... ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4); ALKnew=ALKold+2*delta_t*time_rate_of_change_ALK ... (ALKnow,PO4now,Tnow,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat ... ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4); PO4new=PO4old+2*delta_t*time_rate_of_change_PO4 ... (PO4now,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat ... ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4); end % Robert filter: robert_filter_coeff=0.25; if (solve_physics) Tnow = Tnow + robert_filter_coeff * ( Tnew - 2*Tnow + Told ); Snow = Snow + robert_filter_coeff * ( Snew - 2*Snow + Sold ); end if (solve_geochemistry) ALKnow = ALKnow + robert_filter_coeff * ( ALKnew - 2*ALKnow + ALKold ); PO4now = PO4now + robert_filter_coeff * ( PO4new - 2*PO4now + PO4old ); TCO2now = TCO2now + robert_filter_coeff * ( TCO2new - 2*TCO2now + TCO2old ); end if (solve_physics) Told=Tnow; Tnow=Tnew; Sold=Snow; Snow=Snew; end if (solve_geochemistry) ALKold=ALKnow; ALKnow=ALKnew; PO4old=PO4now; PO4now=PO4new; TCO2old=TCO2now; TCO2now=TCO2new; end % plotting data: Tplot(:,n)=Tnow'; Splot(:,n)=Snow'; ALKplot(:,n)=ALKnow'; PO4plot(:,n)=PO4now'; TCO2plot(:,n)=TCO2now'; uplot(n)=u/1.e12; end time=(1:n)*delta_t/year; figure axes('position', [0.1,0.7,0.8,0.2]); % [left, bottom, width, height]) title(' results of full box model') hold on plot(time,Tplot(1,:),time,Tplot(2,:),time,Tplot(3,:)) hold on legend('T_1','T_2','T_3') ylabel('T'); axes('position', [0.1,0.4,0.8,0.2]); % [left, bottom, width, height]) plot(time,Splot(1,:),time,Splot(2,:),time,Splot(3,:)) hold on legend('S_1','S_2','S_3') ylabel('S'); axes('position', [0.1,0.1,0.8,0.2]); % [left, bottom, width, height]) plot(time,uplot(:)) hold on ylabel('THC'); figure axes('position', [0.1,0.7,0.8,0.2]); % [left, bottom, width, height]) title(' results of full box model') hold on plot(time,ALKplot(1,:),time,ALKplot(2,:),time,ALKplot(3,:)) hold on legend('ALK_1','ALK_2','ALK_3') ylabel('ALK'); axes('position', [0.1,0.4,0.8,0.2]); % [left, bottom, width, height]) plot(time,TCO2plot(1,:),time,TCO2plot(2,:),time,TCO2plot(3,:)) hold on legend('TCO2_1','TCO2_2','TCO2_3') ylabel('TCO2'); axes('position', [0.1,0.1,0.8,0.2]); % [left, bottom, width, height]) plot(time,PO4plot(1,:),time,PO4plot(2,:),time,PO4plot(3,:)) legend('PO4_1','PO4_2','PO4_3') hold on ylabel('PO4'); Tbar=Tnow; Sbar=Snow; ALKbar=ALKnow; PO4bar=PO4now; TCO2bar=TCO2now; ubar=u; CO2_atm_bar=CO2_atm; var_THC=var(uplot(round(n/2):n)); var_T=var(Tplot(:,round(n/2):n)'); var_S=var(Splot(:,round(n/2):n)'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function u=transport(u0,rho) u=-20.0e12; %dbgz: need to make this a function of meridional temp gradient... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function EP=export_production(PO4,hPO4,light_factor,rPO4) % In units of micro mol PO4 per cm^2 per sec: EP=rPO4*light_factor*PO4^2/(hPO4+PO4); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function Rain_ratio=Rain(Temperature) Rain_ratio=61.0/(1+exp(0.1*(10-Temperature))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dT_dt=time_rate_of_change_T ... (n,T,u,Tstar,gamma_T,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat) % rate of change of temperature: % beware of Matlab mess: make sure each of the elements of the vector % dX_dt are sorrounded by (). Otherwise continuation character at end % of line makes matlab think that each line is a different vector % element: if (u>=0) dT_dt=[(gamma_T*(Tstar(1)-T(1)) + u*(T(3)-T(1))/V(1) ... +k_v_highlat*(T(3)-T(1))/H(1)+k_h*(T(2)-T(1))/L(1) ) ... ,(gamma_T*(Tstar(2)-T(2)) + u*(T(1)-T(2))/V(2) ... +k_v_lowlat*(T(3)-T(2))/H(2)+k_h*(T(1)-T(2))/L(2) ) ... ,(u*(T(2)-T(3))/V(3) ... -(k_v_lowlat*(T(3)-T(2))+k_v_highlat*(T(3)-T(1)))/H(3) ) ... ]; else dT_dt=[(gamma_T*(Tstar(1)-T(1)) - u*(T(2)-T(1))/V(1) ... +k_v_highlat*(T(3)-T(1))/H(1)+k_h*(T(2)-T(1))/L(1) ) ... ,(gamma_T*(Tstar(2)-T(2)) - u*(T(3)-T(2))/V(2) ... +k_v_lowlat*(T(3)-T(2))/H(2)+k_h*(T(1)-T(2))/L(2) ) ... ,(-u*(T(1)-T(3))/V(3) ... -(k_v_lowlat*(T(3)-T(2))+k_v_highlat*(T(3)-T(1)))/H(3) ) ... ]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dS_dt=time_rate_of_change_S ... (n,S,u,Sstar,gamma_S,FW,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat) % rate of change of salinity: % beware of Matlab mess: make sure each of the elements of the vector % dX_dt are sorrounded by (). Otherwise continuation character at end % of line makes matlab think that each line is a different vector % element: if (u>=0) dS_dt=[ ((FW(1) ... +gamma_S*(Sstar(1)-S(1))*V(1) +u*(S(3)-S(1)) )/V(1) ... +k_v_highlat*(S(3)-S(1))/H(1)+k_h*(S(2)-S(1))/L(1) ) ... ,((FW(2) ... +gamma_S*(Sstar(2)-S(2))*V(2) +u*(S(1)-S(2)) )/V(2) ... +k_v_lowlat*(S(3)-S(2))/H(2)+k_h*(S(1)-S(2))/L(2) ) ... ,(u*(S(2)-S(3))/V(3) ... -(k_v_lowlat*(S(3)-S(2))+k_v_highlat*(S(3)-S(1)))/H(3) ) ... ]; else dS_dt=[ ((FW(1)... +gamma_S*(Sstar(1)-S(1))*V(1) -u*(S(2)-S(1)))/V(1) ... +k_v_highlat*(S(3)-S(1))/H(1)+k_h*(S(2)-S(1))/L(1) ) ... ,((FW(2)... +gamma_S*(Sstar(2)-S(2))*V(2) -u*(S(3)-S(2)))/V(2) ... +k_v_lowlat*(S(3)-S(2))/H(2)+k_h*(S(1)-S(2))/L(2) ) ... ,(-u*(S(1)-S(3))/V(3) ... -(k_v_lowlat*(S(3)-S(2))+k_v_highlat*(S(3)-S(1)))/H(3) ) ... ]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dALK_dt=time_rate_of_change_ALK ... (ALK,PO4,T,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat... ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4) % rate of change of alkalinity: % need to take care of units here, because biochemistry is in units of % micro mol per liter, and export production in micro mol/ (cm^2 sec). cm3_per_liter=1000.0; S_ALK_surface_high=-Area(1)*(2*Rain(T(1))-1.0/Redfield_PN) ... *cm3_per_liter*export_production(PO4(1),hPO4,light_factor(1),rPO4(1)); S_ALK_surface_low=-Area(2)*(2*Rain(T(2))-1.0/Redfield_PN) ... *cm3_per_liter*export_production(PO4(2),hPO4,light_factor(2),rPO4(2)); S_ALK_deep=-(S_ALK_surface_high+S_ALK_surface_low); % beware of Matlab mess: make sure each of the elements of the vector % dX_dt are sorrounded by (). Otherwise continuation character at end % of line makes matlab think that each line is a different vector % element: if (u>=0) dALK_dt=[ ((S_ALK_surface_high ... +u*(ALK(3)-ALK(1)) )/V(1) ... +k_v_highlat*(ALK(3)-ALK(1))/H(1)+k_h*(ALK(2)-ALK(1))/L(1) ) ... ,((S_ALK_surface_low ... +u*(ALK(1)-ALK(2)) )/V(2) ... +k_v_lowlat*(ALK(3)-ALK(2))/H(2)+k_h*(ALK(1)-ALK(2))/L(2) ) ... ,((S_ALK_deep ... +u*(ALK(2)-ALK(3)))/V(3) ... -(k_v_lowlat*(ALK(3)-ALK(2))+k_v_highlat*(ALK(3)-ALK(1)))/H(3) ) ... ]; else dALK_dt=[ ((S_ALK_surface_high ... - u*(ALK(2)-ALK(1)))/V(1) ... +k_v_highlat*(ALK(3)-ALK(1))/H(1)+k_h*(ALK(2)-ALK(1))/L(1) ) ... ,((S_ALK_surface_low ... - u*(ALK(3)-ALK(2)))/V(2) ... +k_v_lowlat*(ALK(3)-ALK(2))/H(2)+k_h*(ALK(1)-ALK(2))/L(2) ) ... ,((S_ALK_deep ... -u*(ALK(1)-ALK(3)))/V(3) ... -(k_v_lowlat*(ALK(3)-ALK(2))+k_v_highlat*(ALK(3)-ALK(1)))/H(3) ) ... ]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dTCO2_dt=time_rate_of_change_TCO2 ... (TCO2,PO4,T,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat,PV,CO2_atm,CO2_g ... ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4) % rate of change of total co2: % need to take care of units here, because biochemistry is in units of % micro mol per liter, and export production in micro mol/ (cm^2 sec). cm3_per_liter=1000.0; S_TCO2_surface_high=-Area(1)*(1.0/Redfield_PC+Rain(T(1)))... *cm3_per_liter*export_production(PO4(1),hPO4,light_factor(1),rPO4(1)); S_TCO2_surface_low=-Area(2)*(1.0/Redfield_PC+Rain(T(2)))... *cm3_per_liter*export_production(PO4(2),hPO4,light_factor(2),rPO4(2)); S_TCO2_deep=-(S_TCO2_surface_high+S_TCO2_surface_low); % beware of Matlab mess: make sure each of the elements of the vector % dX_dt are sorrounded by (). Otherwise continuation character at end % of line makes matlab think that each line is a different vector % element: if (u>=0) dTCO2_dt=[ ((S_TCO2_surface_high + Area(1)*PV*(CO2_atm-CO2_g(1)) ... +u*(TCO2(3)-TCO2(1)) )/V(1)... +k_v_highlat*(TCO2(3)-TCO2(1))/H(1)+k_h*(TCO2(2)-TCO2(1))/L(1) ) ... ,((S_TCO2_surface_low + Area(2)*PV*(CO2_atm-CO2_g(2)) ... +u*(TCO2(1)-TCO2(2)) )/V(2)... +k_v_lowlat*(TCO2(3)-TCO2(2))/H(2)+k_h*(TCO2(1)-TCO2(2))/L(2) ) ... ,((S_TCO2_deep ... +u*(TCO2(2)-TCO2(3)))/V(3) ... -(k_v_lowlat*(TCO2(3)-TCO2(2))+k_v_highlat*(TCO2(3)-TCO2(1)))/H(3) ) ... ]; else dTCO2_dt=[ ((S_TCO2_surface_high + Area(1)*PV*(CO2_atm-CO2_g(1)) ... -u*(TCO2(2)-TCO2(1)))/V(1) ... +k_v_highlat*(TCO2(3)-TCO2(1))/H(1)+k_h*(TCO2(2)-TCO2(1))/L(1) ) ... ,((S_TCO2_surface_low + Area(2)*PV*(CO2_atm-CO2_g(2)) ... -u*(TCO2(3)-TCO2(2)))/V(2) ... +k_v_lowlat*(TCO2(3)-TCO2(2))/H(2)+k_h*(TCO2(1)-TCO2(2))/L(2) ) ... ,((S_TCO2_deep ... -u*(TCO2(1)-TCO2(3)))/V(3) ... -(k_v_lowlat*(TCO2(3)-TCO2(2))+k_v_highlat*(TCO2(3)-TCO2(1)))/H(3) ) ... ]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function dPO4_dt=time_rate_of_change_PO4 ... (PO4,u,V,Area,H,L,k_h,k_v_highlat,k_v_lowlat ... ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4) % rate of change of nutrient: % need to take care of units here, because biochemistry is in units of % micro mol per liter, and export production in micro mol/ (cm^2 sec). cm3_per_liter=1000.0; S_PO4_surface_high=-Area(1)*cm3_per_liter*export_production(PO4(1),hPO4,light_factor(1),rPO4(1)); S_PO4_surface_low=-Area(2)*cm3_per_liter*export_production(PO4(2),hPO4,light_factor(2),rPO4(2)); S_PO4_deep=-(S_PO4_surface_high+S_PO4_surface_low); % beware of Matlab mess: make sure each of the elements of the vector % dX_dt are sorrounded by (). Otherwise continuation character at end % of line makes matlab think that each line is a different vector % element: if (u>=0) dPO4_dt=[ ((S_PO4_surface_high ... +u*(PO4(3)-PO4(1)) )/V(1)... +k_v_highlat*(PO4(3)-PO4(1))/H(1)+k_h*(PO4(2)-PO4(1))/L(1) ) ... ,((S_PO4_surface_low ... +u*(PO4(1)-PO4(2)) )/V(2)... +k_v_lowlat*(PO4(3)-PO4(2))/H(2)+k_h*(PO4(1)-PO4(2))/L(2) ) ... ,((S_PO4_deep ... +u*(PO4(2)-PO4(3)))/V(3) ... -(k_v_lowlat*(PO4(3)-PO4(2))+k_v_highlat*(PO4(3)-PO4(1)))/H(3) ) ... ]; else dPO4_dt=[ ((S_PO4_surface_high ... -u*(PO4(2)-PO4(1)))/V(1) ... +k_v_highlat*(PO4(3)-PO4(1))/H(1)+k_h*(PO4(2)-PO4(1))/L(1) ) ... ,((S_PO4_surface_low ... -u*(PO4(3)-PO4(2)))/V(2) ... +k_v_lowlat*(PO4(3)-PO4(2))/H(2)+k_h*(PO4(1)-PO4(2))/L(2) ) ... ,((S_PO4_deep ... -u*(PO4(1)-PO4(3)))/V(3) ... -(k_v_lowlat*(PO4(3)-PO4(2))+k_v_highlat*(PO4(3)-PO4(1)))/H(3) ) ... ]; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [CO2_g,H2CO3,OHm,Hp,pH,HCO3m,CO3m2,B_OH_3,B_OH_4m,B_T]= ... solve_carbonate_system(ALK,TCO2,T,S,Hp_init_guess) % Inputs: % ------- % alkalinity, total co_2, temperature, salinity, pressure. % Outputs: % -------- % CO2_g = co_2(g) % H2CO3 = H_2co_3^* = H_2co_3+co_2(aq) % OHm = OH^- % Hp = H^+ % pH = -log10(Hp) % HCO3m = HCO_3^- % CO3m2 = CO_3^= % B_OH_3 = B(OH)_3 % B_OH_4m = B(OH)_4^- % solve for all three boxes: for i=1:3 % calculate dissolution coefficients to save operations during % the iterations below: K0(i)=K_0(T(i),S(i)); K1(i)=K_1(T(i),S(i)); KB(i)=K_B(T(i),S(i)); K2(i)=K_2(T(i),S(i)); Kw(i)=K_w(T(i),S(i)); % calculate total Borate from salinity; multiply by 1.e6 to get % it in micromole: rho_0 = 1.024; B_T(i) = 1.e6*1.212*10^-5*S(i)*rho_0; % initialize a=1/[H^+] and ALK_C from previous time step: a=1/Hp_init_guess(i); max_carbonate_iterates=20; min_carbonate_iterates=5; iterate=1; epsilon_a=1.0; while iterate 1.0e-6 ) iterate=iterate+1; %% Calculate carbonate alkalinity: ALK_C=ALK(i)-( (a*KB(i)*B_T(i))/(1+a*KB(i)) + a*Kw(i) - 1/a ); %% gamma: gamma=ALK_C/TCO2(i); %% solve for a=1/[H^+]: a_new=(-(1-gamma)*K1(i)... +sqrt((1-gamma)^2*K1(i)^2 ... + 4*(2-gamma)*K1(i)*K2(i)*gamma)... )/(2*(2-gamma)*K1(i)*K2(i)); %% evaluate relative increment to a in this iteration: epsilon_a=abs(a_new-a)/a; a=a_new; end if iterate == max_carbonate_iterates fprintf(1,['reached max carbonate system iterates,' ... ' t=%g, iterate=%d, Alk=%8.2f, DIC=%8.2f\n'],t,iterate,ALK,TCO2); end %% Now solve for all other carbonate system variables: Hp(i) = 1/a; % add 6 to pH because units of Hp are micro mole/liter: pH(i) = -log10(Hp(i))+6; OHm(i) = Kw(i)*a; H2CO3(i) = (1/(a*K1(i)+a^2*K1(i)*K2(i)+1))*TCO2(i); HCO3m(i) = ((a*K1(i))/(a*K1(i)+a^2*K1(i)*K2(i)+1))*TCO2(i); CO3m2(i) = ((a^2*K1(i)*K2(i))... /(a*K1(i)+a^2*K1(i)*K2(i)+1))*TCO2(i); CO2_g(i) = H2CO3(i)/K0(i); B_OH_4m(i) = ((a*KB(i))/(1+a*KB(i)))*B_T(i); B_OH_3(i) = B_OH_4m(i)/(a*KB(i)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [L,H,W,V,Area,Sv,year,u0,alpha,beta,gamma_T,gamma_S ... ,k_h,k_v_highlat,k_v_lowlat,PV ... ,FW,Tstar,Sstar,delta_t,N_max,ubar,T,S,ALK,TCO2,PO4,Hp ... ,Redfield_PC,Redfield_PN,hPO4,light_factor,rPO4 ... ,solve_physics,solve_geochemistry] ... =set_parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% cgs!! %% %% ===== %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% solve_physics=0; solve_geochemistry=1; % box geometry: % north-south size: L(1)=5000.0e5; L(2)=1000.0e5; L(3)=L(1); % depth: H(1)=1.0e5; H(2)=H(1); H(3)=4.0e5; if (H(1) ~= H(2)) disp(['need to check that advection/ diffusion equation are' ... ' consistent with the case of h1!=h2!!']) quit end % width: W=4000.0e5; % volumes: V(1)=L(1)*H(1)*W; V(2)=L(2)*H(2)*W; V(3)=L(3)*H(3)*W; % areas: Area(1)=L(1)*W; Area(2)=L(2)*W; % parameters for linearized model matrix: Sv=1.e12; year=3600.0*24.0*365.0; % parameter for the relation between u (THC) and the density % gradients: u0=16*Sv/0.0004; % 16*Sv/0.0005; % parameters for linearized equation of state: alpha=1668.0e-7; beta=0.781e-3; % set SST restoring time in days per 50 meters upper layer % thickness: days_per_50m_mixed_layer=50.0; gamma_T=1.0/(300.0*86400.0); gamma_S=0.0/((days_per_50m_mixed_layer*86400.0)*(H(1)/50.e2)); FW(1)=-1.1*(100.0/year)*35.0*Area(1); FW(2)=-FW(1); % the temperature to which the model surface temperature is restored: Tstar=[0.0,22.0]; % the salinity to which the model surface salinity is restored (only % if gamma_S is not zero): Sstar=[34.0,36.0]; % parameters for numerically solving the steady state problem: % time step: delta_t=5*86400.0; % integration time: integration_time=200.0*year; % number of time steps: N_max= floor((integration_time/delta_t)/100)*100; % Initial conditions: ubar= 20*Sv; T = [1.0,20.0 ,2.0]; S = [34.0,35.5,34.5]; PO4 = [2.09,2.09,2.09]; ALK = [2371,2371,2371]; TCO2 = [2258,2258,2258]; Hp=[10^-8.2,10^-8.2,10^-8.2]; % some biochemistry parameters: Redfield_PC=1.0/122; Redfield_PN=1.0/16; % paraemters for export production parameterization: hPO4=1.0; light_factor=[350.0,200.0]; rPO4=1.0e-10*[0.25,250.0]; % piston velocity for air-sea co2 exchanges PV=1.0*0.59e-2 ; % \approx 5 m/day % mixing coefficients: k_h = 1.0e7/(0.5*(L(1)+L(2))); k_v_highlat= 200.0e0/(0.5*(H(1)+H(3))); k_v_lowlat = 0.0e0/(0.5*(H(1)+H(3))); % Test the density calculation: if 0 pressure=1.0; for i=1:3 rho_linear(i)=linear_density(alpha,beta,T(i),S(i)); rho(i)=density(T(i),S(i),pressure)/1000.0-1; end disp(' test T,S,density:') T,S,rho disp(' test T,S,linear density:') T,S,rho_linear end % Test the geochemical coefficient calculation and the solution of the % carbonate system: if 0 disp('Testing geochemical coefficient calculation:') T1=18,S1=35 Kw=K_w(T1,S1),k0=K_0(T1,S1),k1=K_1(T1,S1),k2=K_2(T1,S1),kb=K_B(T1,S1) disp('Testing solution of carbonate system:') Hp_init_guess=Hp; T,S [CO2_g,H2CO3,OHm,Hp,pH,HCO3m,CO3m2,B_OH_3,B_OH_4m,B_T]= ... solve_carbonate_system(ALK,TCO2,T,S,Hp_init_guess) disp('Testing export production calculation (printed in units of gr C/yr/m^2):') EP1=export_production(PO4(1),hPO4,light_factor(1),rPO4(1)); EP2=export_production(PO4(2),hPO4,light_factor(2),rPO4(2)); % our units are PO4/cm^2/sec. Need to convert to gr C/year and to m^2: export_production_gr_C_per_m2=... [EP1*(122+Rain(Tnow(1))),EP2*(122+Rain(Tnow(2)))] ... *(12.011*10^-6)*(365*86400)*100 export_production_Gt_C_per_yr=... [EP1*(122+Rain(Tnow(1)))*Area(1),EP2*(122+Rain(Tnow(2)))*Area(2)] ... *(12.011*10^-6)*(365*86400)*1.e-15 disp('Testing rain ratio paramterization:') T,rain_T1=Rain(T(1)),rain_T2=Rain(T(2)),rain_T3=Rain(T(3)) end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function rho=linear_density(alpha,beta,T,S) rho=-alpha*T+beta*S; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function rsw=density(t,s,pp) %------------------------------------------------------------------- % PH_DENSITY $Revision: 1.0 $ $Date: 1998/12/07 $ % Copyright (C) GEOMAR, Roger Luff 1998. % % USAGE: rsw = ph_density(t,s,pp) % % DESCRIPTION: % density calculates the density of seawater using the unesco % (1983) formulae, which are valid for a temperature range % of 0-40 deg c and 0.5-43 ppt % % INPUT: % t : temperature of seawater at the specific location [deg C] % s : salinity of seawater at the specific location [PSU] % pp : pressure at the specific location [atm, (with 1 atm % at sea surface)] % % OUTPUT: % rsw : density of a "seawater" [kg/m^3] % % ORIGINAL AUTHOR: % Bernard Boudreau % Dept. Oceanography % Dalhousie University % Halifax, NS B3H 4J1,Canada % bernie@boudreau.ocean.dal.ca % MATLAB VERSION % Roger Luff 98-12-07 % GEOMAR Research Centre for Marine Geosciences % Wischhofstr. 1-3, D-24148 Kiel, Germany % visumod@technologist.com % % DISCLAIMER: % This software is provided "as it is" without warranty of any kind. % % REFERENCES: % Boudreau B. P. (1996) "A method-of-lines code for carbon and % nutrient diagenesis in aquatic sediments". Computers & Geosciences % 22(5), 479-496. % Luff, Haeckel and Wallmann (submitted) "Robust and fast FORTRAN % and MATLAB libraries to calculate pH distributions in % marine systems", Computers & Geosciences. % Fofonoff P. and Millard R. C. (1983) "Algorithms for computation % of fundamental properties of seawater." UNESCO Technical Paper % in Marine Science 44, 53. %------------------------------------------------------------------- %------------------------------------------------------------------- % CHECK INPUT ARGUMENTS %------------------------------------------------------------------- if nargin ~=3 error('ph_density.m: Must pass 3 parameters') end %------------------------------------------------------------------- % convert pressure to bars of sea pressure only %------------------------------------------------------------------- p = (pp-1.0)*1.013; %------------------------------------------------------------------- % density of pure reference water %------------------------------------------------------------------- a0 = 999.842594; a1 = 6.793952d-02; a2 = - 9.095290d-03; a3 = 1.001685d-04; a4 = - 1.120083d-06; a5 = 6.536332d-09; rw = ((((a5*t + a4)*t + a3)*t + a2)*t + a1)*t + a0; %------------------------------------------------------------------- % density at t, s and 1 atm %------------------------------------------------------------------- b0 = 8.24493d-01; b1 = - 4.0899d-03; b2 = 7.6438d-05; b3 = - 8.2467d-07; b4 = 5.3875d-09; c0 = - 5.72466d-03; c1 = 1.0227d-04; c2 = - 1.6546d-06; d0 = 4.8314d-04; rts = rw + ((((b4*t + b3)*t + b2)*t +b1)*t + b0)*s + ((c2*t + c1)*t + c0)*sqrt(s*s*s) + d0*s*s; if p==0.0 rsw = rts; else %------------------------------------------------------------------- % density at t, s and p %------------------------------------------------------------------- e0 = 19652.21; e1 = 148.4206; e2 = - 2.327105; e3 = 1.360477d-02; e4 = - 5.155288d-05; kw = (((e4*t + e3)*t + e2)*t + e1)*t + e0; h0 = 3.239908; h1 = 1.43713d-03; h2 = 1.16092d-04; h3 = - 5.77905d-07; aw = ((h3*t + h2)*t + h1)*t + h0; k0 = 8.50935d-05; k1 = - 6.12293d-06; k2 = 5.2728d-08; bw = (k2*t + k1)*t + k0; m0 = - 9.9348d-07; m1 = 2.0816d-08; m2 = 9.1697d-010; b = bw + ((m2*t + m1)*t + m0)*s; i0 = 2.2838d-03; i1 = - 1.0981d-05; i2 = - 1.6078d-06; j0 = 1.91075d-04; a = aw + ((i2*t + i1)*t + i0)*s + j0*sqrt(s*s*s); f0 = 54.6746; f1 = - 0.603459; f2 = 1.09987d-02; f3 = - 6.1670d-05; g0 = 7.944d-02; g1 = 1.6483d-02; g2 = - 5.3009d-04; kst0 = kw + (((f3*t + f2)*t + f1)*t + f0)*s + ((g2*t + g1)*t + g0)*sqrt(s*s*s); kstp = kst0 + a*p + b*p*p; rsw = rts/(1.0 - p/kstp); end; % These functions calculate the dissociation coefficients needed for % calculation of pCO2, and the solubility of oxygen. function output=K_0(T,S) % solubility of CO2 in seawater % [CO2*] % K_0=-------- % [CO2g] % from Weiss (1974) % Eli: looks like units are mole CO2/(kg Sea Water*atm) % which is the same as the units used here: micro mole CO2/(kg SW*ppm) T=T+273.15; output=exp( ... (9345.17/T) -60.2409+23.3585*log(T/100.) ... +S*( 0.023517-0.00023656*T + ... 0.00000047036*T*T )); function output=K_1(T,S) % [H+][HCO3-] % K_1=----------- % [CO2*] % from Millero (1987) T=T+273.15; output=10.^(-((6320.18/T)-126.3405+19.568*log(T)+ ... (19.894-(840.39/T)-3.0189*log(T))* ... sqrt(S)+0.00678*S) ); % convert units to micromole: output=output*1.e6; function output=K_2(T,S) % [H+][CO3--] % and K_2=----------- from Millero (1987) % [HCO3-] T=T+273.15; output=10.^(-((5143.69/T)-90.1833+14.613*log(T)+ ... (17.176-(690.59/T)-2.6719*log(T))* ... sqrt(S)+0.02169*S) ); % convert units to micromole: output=output*1.e6; function output=K_B(T,S) % for Boric acid, [H+][B(OH4)-] % K_B=-------------- is from Johansson & Wedborg (1982) % [B(OH)3 T=T+273.15; output=10.^(-((1030.5/T)+5.5076-0.015469*S+ ... 1.5339d-4*S*S ) ); % convert units to micromole: output=output*1.e6; function output=K_w(T,S) % for water. K_w=[h+]*[OH-], from Dickson & Riley (1979) % pKw=3441.0*(1./T)+2.241-0.09415*sqrt(S) T=T+273.15; output=10.^(-((3441.0/T)+2.241-0.09415*sqrt(S) )); % convert units to micromole: output=output*1.e12;