Skip to content

Example Eddy Covariance Flux Calculation

This Matlab® 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 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)

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.

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:

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

Boudreau, B.P. and Jorgensen, B.B. eds., 2001. The benthic boundary layer: Transport processes and biogeochemistry. Oxford University Press.

Aubinet, M., Vesala, T. and Papale, D. eds., 2012. Eddy covariance: a practical guide to measurement and data analysis. Springer Science & Business Media.

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.

Created by: Matthew Long 5 June 2013 Modified by: Brian Collister 7 August 2017

Contents

The Eddy Covariance Method

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:

Define File Structures

First, define the input/output files to be used for each calculation. 

infile = 'outfile1.dat'; %input filename with instrument counts corrected to standard units
outfile2 = 'eddyoutfile2.dat'; %output filename for fluxes (mmol m-2 h-1)and burst means
outfile3 = 'eddyoutfile3.dat'; %output filename with cumulative fluxes (mmol m-2 H-1)
outfile4 = 'outfile2O2.dat'; %O2 mean for block averaging

Set Instrument Parameters and Constants

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.

mheight = 0.30; % measuring height (m)
hz = 8; % frequency of measured data (Hz)
windowsize = 2401; % averaging window for running mean
B = 4; % number of bursts per hour (e.g. 15 minute bursts = 4)

% input timestamps for beginning and end of each burst
timestamps = [

13.759	13.9999
14.009	14.2499
14.259	14.4999
14.509	14.7499
14.759	14.9999
15.009	15.158
15.259	15.4999
15.509	15.7499
15.759	15.9999
16.009	16.2499
16.259	16.4999
16.509	16.7499
16.759	16.9999
17.009	17.2499
17.259	17.4999
17.509	17.7499

];
n = length(timestamps);

Import Data File

Now, import the calibrated instrument data file (defined by the variable "infile")

data = importdata(infile,' ',1); % load the data from the user defined input file

% data columns
time = data.data(:,1); % instrument time (hr)
vx = data.data(:,2); vy = data.data(:,3); vz = data.data(:,4); % velocity measurements in the x,y,z direction (ms^-1)
o2 = data.data(:,5);  % oxygen (uM)
pres = data.data(:,6); % CTD pressure measurements (dbar)
SNR = data.data(:,7); % Signal:Noise
correlation = data.data(:,8);
clear data;

Preallocate Flux Variables

% Timestamp Indices
j1 = zeros(n,1); j2 = zeros(n,1);


% Measured Parameter Variables
vxmean = zeros(n,1); vymean = zeros(n,1); vzmean = zeros(n,1);  % average velocity in the (x,y,z) directions
o2mean = zeros(n,1); % average oxygen concentration (uM)
presmean = zeros(n,1); % average CTD pressure (dbar)
SNRmean = zeros(n,1); correlationmean = zeros(n,1); % waveheight (m)
flux1o2 = zeros(n,1);  % O2 flux using method 1
flux2o2 = zeros(n,1);  % O2 flux using method 2
flux3o2 = zeros(n,1);  % O2 flux using method 3
timemean = zeros(n,1); % average time ????? why?

% Find the index of burst begin and end (using timestamps)
for i = 1:n;
    j1(i,1) = find(time > timestamps(i,1),1);
    j2(i,1) = find(time < timestamps(i,2),1,'last');
end


% Calculate means for each burst
for i = 1:n;
    timemean(i) = mean(time(j1(i):j2(i)));
    vxmean(i) = mean(vx(j1(i):j2(i)));
    vymean(i) = mean(vy(j1(i):j2(i)));
    vzmean(i) = mean(vz(j1(i):j2(i)));
    o2mean(i) = mean(o2(j1(i):j2(i)));
    presmean(i) = mean(pres(j1(i):j2(i)));
    SNRmean(i) = mean(SNR(j1(i):j2(i)));
    correlationmean(i) = mean(correlation(j1(i):j2(i)));
end



% Calculate the mean velocity (vmean) for each burst
vmean = (vxmean.^2+vymean.^2+vzmean.^2).^(1/2); % mean velocity (ms^-1)

Calculate Flux by Method #1:

%calculate fluxes by subtracting the mean for each time burst from the measurements along that time burst

vxprime1 = zeros(length(time),1);
vyprime1 = zeros(length(time),1);
vzprime1 = zeros(length(time),1);
o2prime1 = zeros(length(time),1);


for i = 1:n;
	vxprime1(j1(i):j2(i)) = vx(j1(i):j2(i))-vxmean(i);
	vyprime1(j1(i):j2(i)) = vy(j1(i):j2(i))-vymean(i);
	vzprime1(j1(i):j2(i)) = vz(j1(i):j2(i))-vzmean(i);
	o2prime1(j1(i):j2(i)) = o2(j1(i):j2(i))-o2mean(i);
	flux1primeo2 = vzprime1.*o2prime1;
	flux1o2(i) = mean(flux1primeo2(j1(i):j2(i)))*60*60;
	end

Calculate Flux by Method #2:

%calculate deviations as the residuals of the linear regression of x vs. time for each time burst

vxprime2 = zeros(length(time),1);
vyprime2 = zeros(length(time),1);
vzprime2 = zeros(length(time),1);
o2prime2 = zeros(length(time),1);
o2fit = zeros(n,2);
vxfit = zeros(n,2); vyfit = zeros(n,2); vzfit = zeros(n,2);


for i = 1:n;
	o2fit(i,:) = polyfit(time(j1(i):j2(i)),o2(j1(i):j2(i)),1);
	vxfit(i,:) = polyfit(time(j1(i):j2(i)),vx(j1(i):j2(i)),1);
	vyfit(i,:) = polyfit(time(j1(i):j2(i)),vy(j1(i):j2(i)),1);
	vzfit(i,:) = polyfit(time(j1(i):j2(i)),vz(j1(i):j2(i)),1);
	vxprime2(j1(i):j2(i)) = vx(j1(i):j2(i))-(vxfit(i,1)*time(j1(i):j2(i))+vxfit(i,2));
	vyprime2(j1(i):j2(i)) = vy(j1(i):j2(i))-(vyfit(i,1)*time(j1(i):j2(i))+vyfit(i,2));
	vzprime2(j1(i):j2(i)) = vz(j1(i):j2(i))-(vzfit(i,1)*time(j1(i):j2(i))+vzfit(i,2));
	o2prime2(j1(i):j2(i)) = o2(j1(i):j2(i))-(o2fit(i,1).*time(j1(i):j2(i))+o2fit(i,2));
	flux2primeo2 = vzprime2.*o2prime2;
	flux2o2(i) = mean(flux2primeo2(j1(i):j2(i)))*60*60;
	flux2primevx(j1(i):j2(i)) = vzprime2(j1(i):j2(i)).*vxprime2(j1(i):j2(i));
	flux2primevy(j1(i):j2(i)) = vzprime2(j1(i):j2(i)).*vyprime2(j1(i):j2(i));

end

Calculate flux by Method 3:

calculate xprime by subtracting the data smoothed by running average of size ("windowsize') from the raw data

vxprime3 = zeros(length(time),1);
vyprime3 = zeros(length(time),1);
vzprime3 = zeros(length(time),1);
o2prime3 = zeros(length(time),1);


for i = 1:n;
	vxprime3(j1(i):j2(i)) = vx(j1(i):j2(i))-smooth(vx(j1(i):j2(i)),windowsize);
	vyprime3(j1(i):j2(i)) = vy(j1(i):j2(i))-smooth(vy(j1(i):j2(i)),windowsize);
	vzprime3(j1(i):j2(i)) = vz(j1(i):j2(i))-smooth(vz(j1(i):j2(i)),windowsize);
	o2prime3(j1(i):j2(i)) = o2(j1(i):j2(i))-smooth(o2(j1(i):j2(i)),windowsize);
	flux3primeo2 = vzprime3.*o2prime3;
	flux3o2(i) = mean(flux3primeo2(j1(i):j2(i)))*60*60;
	flux3primevx(j1(i):j2(i)) = vzprime3(j1(i):j2(i)).*vxprime3(j1(i):j2(i));
	flux3primevy(j1(i):j2(i)) = vzprime3(j1(i):j2(i)).*vyprime3(j1(i):j2(i));

end

Calculate Cumulative Fluxes

cum1o2 = zeros(length(time)+n,1);
cum2o2 = zeros(length(time)+n,1);
cum3o2 = zeros(length(time)+n,1);
time1 = zeros(length(time)+n,1);

k = 1;
for i = 1:n;
    for j = j1(i):j2(i);
        if j == j1(i);
            cum1o2(k) = flux1primeo2(j);
            cum2o2(k) = flux2primeo2(j);
            cum3o2(k) = flux3primeo2(j);
            time1(k) = time(j);
            k = k+1;
        else
        time1(k) = time(j);
        sum1o2 = cum1o2(k-1)+flux1primeo2(j);
        cum1o2(k) = sum1o2;
        sum2o2 = cum2o2(k-1)+flux2primeo2(j);
        cum2o2(k) = sum2o2;
        sum3o2 = cum3o2(k-1)+flux3primeo2(j);
        cum3o2(k) = sum3o2;
        k = k+1;
        end
        if j == j2(i);
            cum1o2(k) = NaN;
            cum2o2(k) = NaN;
            cum3o2(k) = NaN;
            time1(k) = NaN;
            k = k+1;
        end
    end
end

Write Flux Calculations to Files

%Convert uMol cm L-1 s-1 at 16Hz to mmol m-2
cum1o2b = cum1o2.*0.01./hz.*B;
cum2o2b = cum2o2.*0.01./hz.*B;
cum3o2b = cum3o2.*0.01./hz.*B;

%write data to files
y1 = [timemean, vmean, o2mean, presmean, flux1o2, flux2o2, flux3o2, SNRmean, correlationmean, vxmean, vymean, vzmean];
y1 = y1';

fid1 = fopen(outfile2, 'w');


fprintf(fid1, 'timemean vmean o2mean presmean flux1o2 flux2o2 flux3o2 SNRmean correlationmean vxmean vymean vzmean \r\n');
fprintf(fid1,'%f %f %f %f %f %f %f %f %f %f %f %f \r\n', y1);
fclose(fid1);

Write Cumulative Flux Data to Files

y2 = [time1, cum1o2b, cum2o2b, cum3o2b];
k = find(time1,1,'last');
y2 = y2(1:k,:);
y2 = y2';

fid2 = fopen(outfile3, 'w');

fprintf(fid2, 'time1 cum1o2b cum2o2b  cum3o2b \r\n');
fprintf(fid2,'%f %f %f %f\r\n', y2);
fclose(fid2);

Published with MATLAB® R2016b