MATLAB Examples

Unstable and neutral stratifications

The present example illustrates the calculation of the Obukhov length using the eddy-covariance method. The input is made of time series of the along-wind and vertical wind velocity components, as well as the potential temperature in three different positions on a "simulated" met mast.

A comparison between the target and simulated spectra is provided, where the target spectra is taken from Højstrup [1], which is used to generate the wind spectra in the surface layer for unstable and neutral stratifications. No model for the stable stratification is provided, and this case is therefore not investigated in the following.

[1] Højstrup, J. (1981). A simple model for the adjustment of velocity spectra in unstable conditions downstream of an abrupt change in roughness and heat flux. Boundary-Layer Meteorology, 21(3), 341-356.

Contents

Notations

  1. theta: potential temperature (Kelvin)
  2. u: along-wind wind component (m/s)
  3. v: across-wind wind component (m/s)
  4. w: vertical wind component (m/s)
  5. u_star: friction velocity (m/s)
  6. T_star: scaling temperature (K)
  7. phiM: similarity function for the momentum flux (wind shear) - no dimensions
  8. zi: Lowest inversion height (m)
  9. t: time vector (s)

Assumptions

  1. The humidity flux is not modelled here, so that the humidity is assumed negligeable. In other words, virtual potential temperature = potential temperature.
  2. The turbulent flow is a stationary, homogeneous ergodic random process.
  3. The Monin-Obukhov Similarity Theory (MOST) is assumed to be verified.

Definitions

Obukhov length [2] :

$L = -\frac{\overline{\theta_v} u^3_{*_0}}{g \kappa (\overline{w^{\prime}\theta_v^{\prime}})_0}$

Non-dimensional stability parameter (zL):

$\zeta = \frac{z}{L}$

Friction velocity (u_star):

$u_{*} = \left( \overline{u^{\prime}w^{\prime}}^{2}+\overline{u^{\prime}v^{\prime}}^{2}\right)^{1/4}$

Scaling temperature (T_star):

$T_{*} = -\frac{\overline{w^{\prime}\theta_{v}^{\prime}}}{u_*}$

Non-dimensional velocity profile:

$\phi_m = \frac{\kappa z}{u_*}\frac{\partial\overline{u}}{\partial z}.$

Empirical formulation of $\phi_m$ [3,4]:

$\mathrm{For} \hspace{1em} -2 \leq \zeta \leq 0$:

$\phi_m = \left(1+15.2 |\zeta|\right)^{-1/4}$

$\mathrm{For} \hspace{1em} 0 \leq \zeta \leq 1$:

$\phi_m =  1 + 4.8\zeta$

Ref:

[2] Monin, A. S., & Obukhov, A. M. F. (1954). Basic laws of turbulent mixing in the surface layer of the atmosphere. Contrib. Geophys. Inst. Acad. Sci. USSR, 151(163), e187.

[3] Högström, U. L. F. (1988). Non-dimensional wind and temperature profiles in the atmospheric surface layer: A re-evaluation. Boundary-Layer Meteorology, 42(1), 55-78.

[4] Dyer, A. J. (1974). A review of flux-profile relationships. Boundary-Layer Meteorology, 7(3), 363-372.

Constants

clearvars;close all;clc;
kappa = 0.4; % von karman constant
COLOR = {'k','r','b'};

Load simulated data

Each sample contains 6 variables: t,u,w,nodes,theta and zi. The
variable "nodes" is a structure containing information about the
simulated mast. The across-wind component is neglected
data = dir('Sample_L_*');

clear z meanU L wT *_star U
v = []; % The across-wind component is neglected
for ii=1:numel(data),
    load(data(ii).name);
    z(ii,:) = nodes.Z;
    U(ii,:) = nanmean(u);
    Theta(ii,:) = nanmean(theta);
    fprintf(['Summary of data set ',num2str(ii),': \n']);
    [L(ii,:),wT(ii,:),T_star(ii,:),u_star(ii,:)] = getStability(u,w,theta,nodes.Z,'v',v); % the 'v' parameter is optional

end

z = mean(z);% z is the same for every samples
zL = bsxfun(@times,z,1./L);
% According to the MOST, wT and u_star are constant in the sruface layer:
u_star = nanmean(u_star,2);
T_star = nanmean(T_star,2);
wT = nanmean(wT,2);
Summary of data set 1: 

Measured T* = -0.016 K 
Measured u* = 0.62 m/s 
Measured heat flux wT = 0.01 
Measured L = -1663 m 

Summary of data set 2: 

Measured T* = -1.7 K 
Measured u* = 0.63 m/s 
Measured heat flux wT = 1.1 
Measured L = -16 m 

Summary of data set 3: 

Measured T* = -0.78 K 
Measured u* = 0.63 m/s 
Measured heat flux wT = 0.49 
Measured L = -36 m 

Summary of data set 4: 

Measured T* = -3.1 K 
Measured u* = 0.6 m/s 
Measured heat flux wT = 1.8 
Measured L = -8 m 

Overview of the "nodes" on the simulated met-mast

disp(nodes)

figure
hold on, box on;
plot([0,0],[0,15],'color',[0.5,0.5,0.5],'linewidth',2);
for ii=1:numel(z),
    h(ii) = plot(nodes.Y(ii),nodes.Z(ii),'ko','markerfacecolor',COLOR{ii});
end
legend(h,nodes.name{:})
hold on
ylim([0,18])
ylabel('$z$ (m)','interpreter','latex')
xlabel('$y$ (m)','interpreter','latex')
set(gcf,'color','w')
title('Location of the 3 nodes simulating the anemometers in N1, N2 and N3')
annotation('textarrow',[0.3,0.51],[0.55,0.55],'String','met-mast ','interpreter','latex')
       U: [5.9317 6.8365 7.1847]
       Y: [0 0 0]
       Z: [3 9 15]
    name: {'N1'  'N2'  'N3'}

Comparison with similarity functions for momentum

The reference similarity functions are taken from Högström [3], based on the
works of Dyer et al. [4]
figure
plotSimilarityFun(z,zL,u_star,T_star,U,Theta);

Overview of the target spectra for the along-wind velocity component

Each panel corresponds to a different stability, characeterized by a given value of the Obukhov length L.

Normalization = 1; % We want the normalized spectra (otherwise Normalization =  0)
f = logspace(-3,0.5,100);
Su = zeros(numel(data),numel(z),numel(f));
Sw = zeros(numel(data),numel(z),numel(f));


for ii=1:numel(data),
    for jj=1:numel(z),
        fr = f*z(jj)/U(ii,jj);
        fi = f*zi/U(ii,jj); % zi is the inversion height
        Su(ii,jj,:) = longitudinalSpectrum_Hojstrup(u_star(ii),L(ii),f,fr,fi,zi,Normalization);
        Sw(ii,jj,:) = verticalSpectrum_Hojstrup(u_star(ii),z(jj),L(ii),f,fr,Normalization);
    end
end

clf;close all
figure('units','pixels','outerposition',  [560 10 560 1000])
for ii=1:numel(data),
    subplot(4,1,ii)
    for jj=1:numel(z),
        hold on;box on;
        plot(f*z(jj)/U(ii,jj),squeeze(Su(ii,jj,:)),COLOR{jj});
        xlabel('$nz/\overline{u}$','interpreter','latex')
        ylabel('$n S_u / u^2_*$','interpreter','latex')
        grid on
        ylim([0.1,5]);xlim([1e-4,2]);
        title(['L = ',num2str(round(L(ii))),' m'],'FontSize',8)
    end
    set(gca,'Xscale','log','Yscale','log');
    if ii==4,
        legend('z = 40 m','z = 60 m','z = 80 m','location','south')
    end
end

set(gcf,'color','w')

Overview of the target spectra for the vertical wind velocity component

COLOR = {'k','r','b'};
clf;close all
figure('units','pixels','outerposition',  [560 10 560 1000])
for ii=1:numel(data),
    subplot(4,1,ii)
    for jj=1:numel(z),
        hold on;box on;
        plot(f*z(jj)/U(ii,jj),squeeze(Sw(ii,jj,:)),COLOR{jj});
        xlabel('$nz/\overline{u}$','interpreter','latex')
        ylabel('$nS_w/u^2_*$','interpreter','latex')
        grid on
        ylim([0.0005,2]);xlim([1e-4,2])
        title(['L = ',num2str(round(L(ii))),' m'],'FontSize',8)
    end
    set(gca,'Xscale','log','Yscale','log');
    if ii==4,
        legend('z = 40 m','z = 60 m','z = 80 m','location','south')
    end
end

set(gcf,'color','w')

Computed vs target spectra

Comparison of the target and computed spectra, for each stability class.

clf;close all;
for ii=1:numel(data),

    load(data(ii).name,'u','w','t');
    fs = 1./median(diff(t));
    figure
    plotVelocitySpectra(u,w,fs,squeeze(Su(ii,1,:))',squeeze(Sw(ii,1,:))',f,u_star(ii))
    title(['L = ',num2str(round(L(ii))),' m'],'FontSize',8)
end