This Matlab\u00ae code represents an example code for processing and analyzing oxygen flux data from an Eddy Covariance Oxygen System. I recommend that everyone develop their own code because:<\/p>\n
1.) Every covariance system and application are different<\/p>\n
2.) Corrections like rotations, time-lag corrections, storage, coordinate transformations, and additional scalars are not included here as these are again system specific but also complicate the underlying principles (I am happy to share my most recent version of the code that includes these – mlong@whoi.edu<\/a>)<\/p>\n 3.) The small underwater eddy covariance community is constantly improving their capabilities and accuracies – being able to adapt, modify, and test your code is imperative.<\/p>\n For detailed descriptions of the eddy covariance method and its application to measurements of fluxes at the sediment\/water interface, please refer to the following peer-review articles:<\/p>\n Baldocchi, D. D. 2003. Assessing the eddy covariance technique for evaluating carbon dioxide exchange rates of ecosystems: Past, present and future. Glob. Chang. Biol. 9: 479-492. doi:10.1046\/j.1365-2486.2003.00629.x<\/p>\n Boudreau, B.P. and Jorgensen, B.B. eds., 2001.\u00a0The benthic boundary layer: Transport processes and biogeochemistry<\/i>. Oxford University Press.<\/p>\n Aubinet, M., Vesala, T. and Papale, D. eds., 2012.\u00a0Eddy covariance: a practical guide to measurement and data analysis<\/i>. Springer Science & Business Media.<\/p>\n Long, Matthew H., et al. “Oxygen metabolism and pH in coastal ecosystems: Eddy Covariance Hydrogen ion and Oxygen Exchange System (ECHOES).” Limnology and Oceanography: Methods 13.8 (2015): 438-450.<\/p>\n Created by: Matthew Long 5 June 2013 Modified by: Brian Collister 7 August 2017<\/p>\n The basis for the eddy covarianve technique is that turbulent mixing, caused by the interaction of current flow with the interaction of current flow with the benthic surface, is the dominant vertical transport process in the ocean. Therefore, vertical fluxes accross the sediment-water interface can be derived from high-resolution measurements of the covariation between vertical velocity and a solute concentration. The average solute flux (C) is calculated as:<\/p>\n\n First, define the input\/output files to be used for each calculation.\u00a0<\/p>\ninfile = ‘outfile1.dat’; %input filename with instrument counts corrected to standard units\noutfile2 = ‘eddyoutfile2.dat’; %output filename for fluxes (mmol m-2 h-1)and burst means\noutfile3 = ‘eddyoutfile3.dat’; %output filename with cumulative fluxes (mmol m-2 H-1)\noutfile4 = ‘outfile2O2.dat’; %O2 mean for block averaging\n\n The user must define the instrument parameters necessary to correcting for instrument geometry and sampling scheme. In this case, the EC system used a flow through pH sensor, so knowledge of the flow time through the sensor was necessary.<\/p>\nmheight = 0.30; % measuring height (m)\nhz = 8; % frequency of measured data (Hz)\nwindowsize = 2401; % averaging window for running mean\nB = 4; % number of bursts per hour (e.g. 15 minute bursts = 4)\n% input timestamps for beginning and end of each burst\ntimestamps = [\n13.759\t13.9999\n14.009\t14.2499\n14.259\t14.4999\n14.509\t14.7499\n14.759\t14.9999\n15.009\t15.158\n15.259\t15.4999\n15.509\t15.7499\n15.759\t15.9999\n16.009\t16.2499\n16.259\t16.4999\n16.509\t16.7499\n16.759\t16.9999\n17.009\t17.2499\n17.259\t17.4999\n17.509\t17.7499\n];\nn = length(timestamps);\n\n Now, import the calibrated instrument data file (defined by the variable “infile”)<\/p>\ndata = importdata(infile,’ ‘,1); % load the data from the user defined input file\n% data columns\ntime = data.data(:,1); % instrument time (hr)\nvx = data.data(:,2); vy = data.data(:,3); vz = data.data(:,4); % velocity measurements in the x,y,z direction (ms^-1)\no2 = data.data(:,5); % oxygen (uM)\npres = data.data(:,6); % CTD pressure measurements (dbar)\nSNR = data.data(:,7); % Signal:Noise\ncorrelation = data.data(:,8);\nclear data; \n calculate xprime by subtracting the data smoothed by running average of size (“windowsize’) from the raw data<\/p>\nvxprime3 = zeros(length(time),1);\nvyprime3 = zeros(length(time),1);\nvzprime3 = zeros(length(time),1);\no2prime3 = zeros(length(time),1);\nfor i = 1:n;\n\tvxprime3(j1(i):j2(i)) = vx(j1(i):j2(i))-smooth(vx(j1(i):j2(i)),windowsize);\n\tvyprime3(j1(i):j2(i)) = vy(j1(i):j2(i))-smooth(vy(j1(i):j2(i)),windowsize);\n\tvzprime3(j1(i):j2(i)) = vz(j1(i):j2(i))-smooth(vz(j1(i):j2(i)),windowsize);\n\to2prime3(j1(i):j2(i)) = o2(j1(i):j2(i))-smooth(o2(j1(i):j2(i)),windowsize);\n\tflux3primeo2 = vzprime3.*o2prime3;\n\tflux3o2(i) = mean(flux3primeo2(j1(i):j2(i)))*60*60;\n\tflux3primevx(j1(i):j2(i)) = vzprime3(j1(i):j2(i)).*vxprime3(j1(i):j2(i));\n\tflux3primevy(j1(i):j2(i)) = vzprime3(j1(i):j2(i)).*vyprime3(j1(i):j2(i));\nend\n\n Example Eddy Covariance Flux Calculation This Matlab\u00ae code represents an example code for processing and analyzing oxygen flux data from an Eddy Covariance Oxygen System. I recommend that everyone develop their own code because: 1.) Every covariance system and application are different 2.) Corrections like rotations, time-lag corrections, storage, coordinate transformations, and additional scalars are…<\/p>\n","protected":false},"author":84,"featured_media":0,"parent":202,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":[],"_links":{"self":[{"href":"https:\/\/www2.whoi.edu\/staff\/mlong\/wp-json\/wp\/v2\/pages\/302"}],"collection":[{"href":"https:\/\/www2.whoi.edu\/staff\/mlong\/wp-json\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/www2.whoi.edu\/staff\/mlong\/wp-json\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/www2.whoi.edu\/staff\/mlong\/wp-json\/wp\/v2\/users\/84"}],"replies":[{"embeddable":true,"href":"https:\/\/www2.whoi.edu\/staff\/mlong\/wp-json\/wp\/v2\/comments?post=302"}],"version-history":[{"count":3,"href":"https:\/\/www2.whoi.edu\/staff\/mlong\/wp-json\/wp\/v2\/pages\/302\/revisions"}],"predecessor-version":[{"id":579,"href":"https:\/\/www2.whoi.edu\/staff\/mlong\/wp-json\/wp\/v2\/pages\/302\/revisions\/579"}],"up":[{"embeddable":true,"href":"https:\/\/www2.whoi.edu\/staff\/mlong\/wp-json\/wp\/v2\/pages\/202"}],"wp:attachment":[{"href":"https:\/\/www2.whoi.edu\/staff\/mlong\/wp-json\/wp\/v2\/media?parent=302"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}Contents<\/h2>\n
\n
The Eddy Covariance Method<\/h2>\n
Define File Structures<\/h2>\n
Set Instrument Parameters and Constants<\/h2>\n
Import Data File<\/h2>\n
Preallocate Flux Variables<\/h2>\n% Timestamp Indices\nj1 = zeros(n,1); j2 = zeros(n,1);\n% Measured Parameter Variables\nvxmean = zeros(n,1); vymean = zeros(n,1); vzmean = zeros(n,1); % average velocity in the (x,y,z) directions\no2mean = zeros(n,1); % average oxygen concentration (uM)\npresmean = zeros(n,1); % average CTD pressure (dbar)\nSNRmean = zeros(n,1); correlationmean = zeros(n,1); % waveheight (m)\nflux1o2 = zeros(n,1); % O2 flux using method 1\nflux2o2 = zeros(n,1); % O2 flux using method 2\nflux3o2 = zeros(n,1); % O2 flux using method 3\ntimemean = zeros(n,1); % average time ????? why?\n% Find the index of burst begin and end (using timestamps)\nfor i = 1:n;\n j1(i,1) = find(time > timestamps(i,1),1);\n j2(i,1) = find(time < timestamps(i,2),1,’last’);\nend\n% Calculate means for each burst\nfor i = 1:n;\n timemean(i) = mean(time(j1(i):j2(i)));\n vxmean(i) = mean(vx(j1(i):j2(i)));\n vymean(i) = mean(vy(j1(i):j2(i)));\n vzmean(i) = mean(vz(j1(i):j2(i)));\n o2mean(i) = mean(o2(j1(i):j2(i)));\n presmean(i) = mean(pres(j1(i):j2(i)));\n SNRmean(i) = mean(SNR(j1(i):j2(i)));\n correlationmean(i) = mean(correlation(j1(i):j2(i)));\nend\n% Calculate the mean velocity (vmean) for each burst\nvmean = (vxmean.^2+vymean.^2+vzmean.^2).^(1\/2); % mean velocity (ms^-1)\n\n
Calculate Flux by Method #1:<\/h2>\n%calculate fluxes by subtracting the mean for each time burst from the measurements along that time burst\nvxprime1 = zeros(length(time),1);\nvyprime1 = zeros(length(time),1);\nvzprime1 = zeros(length(time),1);\no2prime1 = zeros(length(time),1);\nfor i = 1:n;\n\tvxprime1(j1(i):j2(i)) = vx(j1(i):j2(i))-vxmean(i);\n\tvyprime1(j1(i):j2(i)) = vy(j1(i):j2(i))-vymean(i);\n\tvzprime1(j1(i):j2(i)) = vz(j1(i):j2(i))-vzmean(i);\n\to2prime1(j1(i):j2(i)) = o2(j1(i):j2(i))-o2mean(i);\n\tflux1primeo2 = vzprime1.*o2prime1;\n\tflux1o2(i) = mean(flux1primeo2(j1(i):j2(i)))*60*60;\n\tend\n\n
Calculate Flux by Method #2:<\/h2>\n%calculate deviations as the residuals of the linear regression of x vs. time for each time burst\nvxprime2 = zeros(length(time),1);\nvyprime2 = zeros(length(time),1);\nvzprime2 = zeros(length(time),1);\no2prime2 = zeros(length(time),1);\no2fit = zeros(n,2);\nvxfit = zeros(n,2); vyfit = zeros(n,2); vzfit = zeros(n,2);\nfor i = 1:n;\n\to2fit(i,:) = polyfit(time(j1(i):j2(i)),o2(j1(i):j2(i)),1);\n\tvxfit(i,:) = polyfit(time(j1(i):j2(i)),vx(j1(i):j2(i)),1);\n\tvyfit(i,:) = polyfit(time(j1(i):j2(i)),vy(j1(i):j2(i)),1);\n\tvzfit(i,:) = polyfit(time(j1(i):j2(i)),vz(j1(i):j2(i)),1);\n\tvxprime2(j1(i):j2(i)) = vx(j1(i):j2(i))-(vxfit(i,1)*time(j1(i):j2(i))+vxfit(i,2));\n\tvyprime2(j1(i):j2(i)) = vy(j1(i):j2(i))-(vyfit(i,1)*time(j1(i):j2(i))+vyfit(i,2));\n\tvzprime2(j1(i):j2(i)) = vz(j1(i):j2(i))-(vzfit(i,1)*time(j1(i):j2(i))+vzfit(i,2));\n\to2prime2(j1(i):j2(i)) = o2(j1(i):j2(i))-(o2fit(i,1).*time(j1(i):j2(i))+o2fit(i,2));\n\tflux2primeo2 = vzprime2.*o2prime2;\n\tflux2o2(i) = mean(flux2primeo2(j1(i):j2(i)))*60*60;\n\tflux2primevx(j1(i):j2(i)) = vzprime2(j1(i):j2(i)).*vxprime2(j1(i):j2(i));\n\tflux2primevy(j1(i):j2(i)) = vzprime2(j1(i):j2(i)).*vyprime2(j1(i):j2(i));\nend\n\n
Calculate flux by Method 3:<\/h2>\n
Calculate Cumulative Fluxes<\/h2>\ncum1o2 = zeros(length(time)+n,1);\ncum2o2 = zeros(length(time)+n,1);\ncum3o2 = zeros(length(time)+n,1);\ntime1 = zeros(length(time)+n,1);\nk = 1;\nfor i = 1:n;\n for j = j1(i):j2(i);\n if j == j1(i);\n cum1o2(k) = flux1primeo2(j);\n cum2o2(k) = flux2primeo2(j);\n cum3o2(k) = flux3primeo2(j);\n time1(k) = time(j);\n k = k+1;\n else\n time1(k) = time(j);\n sum1o2 = cum1o2(k-1)+flux1primeo2(j);\n cum1o2(k) = sum1o2;\n sum2o2 = cum2o2(k-1)+flux2primeo2(j);\n cum2o2(k) = sum2o2;\n sum3o2 = cum3o2(k-1)+flux3primeo2(j);\n cum3o2(k) = sum3o2;\n k = k+1;\n end\n if j == j2(i);\n cum1o2(k) = NaN;\n cum2o2(k) = NaN;\n cum3o2(k) = NaN;\n time1(k) = NaN;\n k = k+1;\n end\n end\nend\n\n
Write Flux Calculations to Files<\/h2>\n%Convert uMol cm L-1 s-1 at 16Hz to mmol m-2\ncum1o2b = cum1o2.*0.01.\/hz.*B;\ncum2o2b = cum2o2.*0.01.\/hz.*B;\ncum3o2b = cum3o2.*0.01.\/hz.*B;\n%write data to files\ny1 = [timemean, vmean, o2mean, presmean, flux1o2, flux2o2, flux3o2, SNRmean, correlationmean, vxmean, vymean, vzmean];\ny1 = y1′;\nfid1 = fopen(outfile2, ‘w’);\nfprintf(fid1, ‘timemean vmean o2mean presmean flux1o2 flux2o2 flux3o2 SNRmean correlationmean vxmean vymean vzmean rn’);\nfprintf(fid1,’%f %f %f %f %f %f %f %f %f %f %f %f rn’, y1);\nfclose(fid1);\n\n
Write Cumulative Flux Data to Files<\/h2>\ny2 = [time1, cum1o2b, cum2o2b, cum3o2b];\nk = find(time1,1,’last’);\ny2 = y2(1:k,:);\ny2 = y2′;\nfid2 = fopen(outfile3, ‘w’);\nfprintf(fid2, ‘time1 cum1o2b cum2o2b cum3o2b rn’);\nfprintf(fid2,’%f %f %f %frn’, y2);\nfclose(fid2);\n\n\n\n","protected":false},"excerpt":{"rendered":"