close all clear all %---------------------------------------------------- %This 'model' equation was composed by Lennart T. Bach for mesocosm export %datasets published in Global Biogeochemical Cycles. Please cite the Bach et %al., 2016 and Bach et al., 2019 studies when using this equation. Feel %free to modify and redistribute. %Formulation of a one dimensional export model using the new findings of %the dependence of sinking velocities and remineralization on plankton %community composition %basic equation reads as %POC(z)=POC(z-1) - POC(z-1)*C_model/SV_model*(z-(z-1)) %z>=1 %where z is a depth interval, C_model the remineralization equation, and SV_model the %sinking velocity equation. %---------------------------------------------------- %load data %temperature profile %The temperature profile was measured near Gran Canaria at the %European Station for Time-Series in the Ocean (ESTOC, 29.17 N; -15.50 W) %taken on 26 February 2014 during a cruise with R/V Poseidon (P465) %(°C). This profile was used in the Bach et al., 2019 study. %First column is depth, second column is temperature. depth_temp=csvread('depth_temp_high_res.csv'); %load temperature profile T=depth_temp(:,2); %(in °C) extract temperature vector from the file depth=depth_temp(:,1); %(in m) extract depth vector from the file %Sinking velocity and remineralization data from the Bach et al., 2019 %study. %First column is the day of experiment, second column is the mesocosm, %third column is C_remin, Fourth column is SV_insitu. (You will immediately %understand what these last two variables are when having a look in the %Bach et al., 2019 study. sink_resp=csvread('RLS_data_GC2_model.csv'); SV_insitu=sink_resp(:,4); %(in m/d) C_remin=sink_resp(:,3); %(in 1/d) %---------------------------------------------------- %calculation of C_model (i.e. depth dependent C_remin which declines relative %surface values do to decreasing temperature in the deeper water). %parameterization taken from Schmittner et al., (2008) in GBC %accounts for Q_10 factor %C_remin default value in Schmittner et al., (2008) is 0.048 d^-1 in the surface. %In this model we change this default value according to what has been measured %by Bach et al., (2019). (See C_remin in the sink_resp data file). b=1.038; %dimensionless scaling factor for this study. I tuned this from its %original value of 1.066 in Schmittner et al. down to 1.038 because more recent assessments have %shown that Martin's b is much lower (attenuation is lower) in the NE Atlantic %(Guidi et al., 2015, Henson et al., 2012). If I use the standard UVIC value of %1.066 from Schmittner et al., 2008 I get a too high attenuation %with the average C_remin and SV_insitu values measured in this study. %this loop profides the depth dependent remineralization in the model %(C_model) for all C_remin values. for i=1:1:length(C_remin); for j=1:1:length(T); C_model(i,j)=C_remin(i).*b.^T(j); %1/d (per day) end end C_model=C_model'; %need to transpose this matrix. %---------------------------------------------------- %calculation of SV_model. Note that SV_model depends on SV_insitu but %increases with depth as in the UVIC model. %SV_model = slope*depth+SV_insitu slope=0.04; %Uvic parameterization for the linear increase of sinking velocity % with depth as in the UVIC model (Schmittner et al., 2008). for i=1:1:length(SV_insitu); for j=1:1:length(depth); SV_model(i,j)=(slope.*depth(j)+SV_insitu(i)); %1/d (per day) end end SV_model=SV_model'; %need to transpose this matrix. %---------------------------------------------------- %calculation of attenuation profiles as in Bach et al., 2019. depth_initial = 15; % in meter. This is the depth where the the particle enters %the export pathway. In this case it is the bottom of the mesocosm in the %Bach et al., 2019 study where the mesocosms were 15 m deep. %Wd_control=(slope.*depth+64.37959792); %m/d %max=64.37959792 min=22.93814704 Average = 37.53427568 27.84224639 for i=1:1:length(C_remin); for z=find(depth==depth_initial):length(depth); %index with length of all the other vectors (e.g. T, depth...) if z==find(depth==depth_initial); %determining starting value for POC at depth of euphotic zone. In this case, the value is 1 meaning 100% POC(z,i)=1; elseif z>1; %below the initial depth level POC is allowed to decrease. POC(z,i)=POC(z-1,i)-(POC(z-1,i).*(C_model(z,i)./SV_model(z,i)).*(depth(z)-depth(z-1))); end end end %---------------------------------------------------- %---------------------------------------------------- %use this routine below if you want to test just one combination of C_remin %and SV_insitu. This is nice to play around with parameters. %JUST TAKE AWAY THE "%" BELOW %close all %clear all %load data %depth_temp=csvread('depth_temp_high_res.csv'); %load temperature profile %T=depth_temp(:,2); %(in °C) extract temperature vector from the file %depth=depth_temp(:,1); %(in m) extract depth vector from the file %depth_initial = 15; %(m) %slope=0.04; %SV_insitu=100; %(m/d) %C_remin=0.01; %(1/d) %b=1.038; %SV_model=(slope.*depth+SV_insitu); %m/d %C_model=C_remin.*b.^T; %1/d %for z=find(depth==depth_initial):length(depth); %index with length of all the other vectors (e.g. T, depth...) % if z==find(depth==depth_initial); %determining starting value for POC at depth of euphotic zone. In this case, the value is 1 meaning 100% % POC(z)=1; % elseif z>1; %below the initial depth level POC is allowed to decrease. % POC(z)=POC(z-1)-(POC(z-1).*(C_model(z)./SV_model(z)));%.*(depth(z)-depth(z-1) % % end %end %---------------------------------------------------- %---------------------------------------------------- %plotting POC(POC == 0) = NaN; % to exchange zeros with NaN in the output POC file when %the calcultion starts below the sea surface POC=POC.*100; %to give the relative change in percent plot(POC,depth,'LineStyle','-','Color','k','LineWidth',[1]) set(gca,'YDir','reverse') axis([0 100 0 1500]) xlabel('Remaining POC (%)') ylabel('Depth (m)')