A semi-analytical solution for the mean wind profile in the Atmospheric Boundary Layer : the convective case

A semi-analytical solution for the mean wind profile in the Atmospheric Boundary Layer: the convective case L. Buligon, G. A. Degrazia, O. C. Acevedo, C. R. P. Szinvelski, and A. G. O. Goulart Universidade Federal de Santa Maria, Departamento de Fı́sica,Santa Maria, RS, Brazil Universidade Federal do Pampa/UFSM, Centro de Ciências Exatas e Tecnológicas, Bagé, RS, Brazil Centro de Educação Superior do Alto Vale do Itajaı́, UDESC/CEAVI, Ibirama, Brazil Received: 26 June 2009 – Accepted: 25 August 2009 – Published: 23 September 2009 Correspondence to: G. A. Degrazia (degrazia@ccne.ufsm.br) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
The Navier-Stokes equations provide the framework for the interpretation of atmospheric boundary layer flows.However, their analytical solution requires approximations, which are, in many cases, idealized and distant from the physical reality.Particularly, the mean wind profile is a solution of the governing equations whose derivation can be applied to a wide variety of natural processes.
The classical Ekman expression for the mean wind profile is, probably, the most famous example of an analytical solution of the simplified Navier-Stokes equations (Sorbjan, 1989;Stull, 1988).Such solution needs, however, the strong, non-realistic assumption that the vertical eddy diffusivities (K) are constant with height.In fact, Grisogono (1995) argues that "...it is a complicated, nonlinear function of the flow structure and there is no explicit relation between the boundary-layer profiles and K".
The search for an analytical solution for the mean wind profiles in the atmospheric boundary layer under more realistic conditions has been a major focus of mathematical and physical research for a long time.For a given, imposed eddy Correspondence to: Gervásio A. Degrazia (degrazia@ccne.ufsm.br)diffusivity profile, solutions have been derived in a wide set of studies, reviewed by Zilitinkevich (1970), Monin and Yaglom (1971) and Grisogono (1995).
Other studies have expanded the solution to include baroclinic and advective conditions (Miles, 1994;Bannon and Salem, 1995).Berger and Grisogono (1998) extended the results obtained by Grisogono (1995) for the baroclinic case and a generic vertical eddy diffusivity profile.Tan (2001) proposed a semi-geostrophic Ekman layer solution for variable eddy diffusivities and baroclinicity.This model united the solutions presented by Wu and Blumen (1982) and Grisogono (1995), finding that the mean wind structure depends on the inertial acceleration, eddy diffusivity and baroclinic pressure gradient.The study concluded that anti-cyclonic wind shear accelerates the flow, while cyclonic shear has the opposite effect.Wilson and Flesch (2004) used a three-layer simplified model that provided a good comparison to observed wind profiles.
In the present study, using the Generalized Integral Transform Technique (GITT), we derive a semi-analytical solution of the Navier-Stokes equation to obtain the mean wind profile in the atmospheric boundary layer.Such technique is a hybrid numerical-analytical method applied to the treatment and solution of partial differential equations ( Özis ¸ik, 1993;Mikhailov and Özis ¸ik, 1984;Cotta, 1993).It provides a systematic, direct and efficient approximation to the solution of homogeneous and non-homogeneous, stationary and non-stationary, linear and non-linear boundary-value problems.The technique combines series expansion and an integration employing an inverse-transform pair.The PBL is discretized into N sub-intervals in such manner that inside each sub-region the eddy diffusivity is the average value (Moreira et al., 1999), that allows the use of realistic eddy diffusivity profiles, which depend on the physical characteristics of the energy-containing eddies.The nonlinear terms are written in terms of kinematical properties of the flow, such as divergence and vorticity, allowing the solutions to be interpreted in terms of large-scale synoptic conditions.The model results are compared to observed wind profiles obtained from the classical Wangara experiment (Clarke et al., 1971).

Basic equations
Considering that the turbulent momentum fluxes can be parameterized by a first order closure (K-theory), the mean horizontal wind spatial distribution is given by the Navier-Stokes equation in the following form: The equations above assume stationarity, no mean vertical motion, and that the molecular dissipation terms are neglectable.On the other hand, the flow is allowed to vary horizontally.
The geostrophic wind components in the baroclinic case are approximated by: where u g0 and v g0 are the surface geostrophic winds components and u T and v T are the thermal wind components (Sorbjan, 1989).The eddy diffusivities in each of the directions are represented by K x , K y and K z .
To realistically reproduce the wind profile, it is important to consider the vertical variation of the eddy diffusivities.As a consequence, in the present approach, the planetary boundary layer (PBL) is discretized into N sublayers (Vilhena and Barichello, 1991;Moreira et al., 1999;Degrazia et al., 2001).In each of the sublayers, the eddy diffusivities and horizontal wind components assume vertically averaged values.To overcome the difficulties that arise from the nonlinear character of equations (1a) and (1b), the advective terms are written in terms of the large-scale kinematical properties of the flow (Bluestein, 1992), namely, divergence and vorticity: Where the index n refers to the different sublayers considered.With the assumptions above, equations (1a) and (1b) can be written as: with z n ≤ z ≤ z n+1 , 0 < x < L x , 0 < y < L y and n = 1, 2, . . ., N .
Multiplying equation (4b) by i (i ∈ C), and adding term by term to equation (4a), yields where w n = u n + v n i, w gn = u gn + v gn i, z n ≤ z ≤ z n+1 , 0 < x < L x , 0 < y < L y and n = 1, 2, . . ., N .Equation ( 5) is a differential equation on the complex variable function w n = u n + v n i, whose solution provides an expression to the mean wind profile in terms of the flow divergence and vorticity.

Boundary and interface conditions
The horizontal wind speeds are assumed to be constant at the lower boundary at z = z 0 , and to be geostrophic at the upper boundary, the PBL top.Laterally, a horizontal domain is assumed, with dimensions L x × L y .At the lateral boundaries, the wind components are given by the imposed divergence and vorticity, so that: and At the interfaces between neighbor vertical layers, it is necessary to assume continuity of both the eddy fluxes and horizontal wind components.

Solution
Equation ( 5) can be solved using the Generalized Integral Transform Technique -GITT (Mikhailov and Özis ¸ik, 1984;Özis ¸ik, 1993;Cotta, 1993).In this method, the solution function is expanded in terms of the eigenfunctions corresponding to the auxiliary problem (Sturm-Liouville), associated with the original problem.The eigenfunction orthogonality condition is used to determine the expansion coefficients, hence originating the integral transform and its inverse.Applying the integral transform, the partial derivatives in relation to variables x and y are removed, reducing the problem to an ordinary second-order differential equation on variable z.Therefore, once the transformed problem is solved, the inverse formula is used to obtain a solution to the original problem.The truncation order is selected according to the desired precision.

The auxiliary problem
The auxiliary problem associated with w n is: with boundary conditions: ψ = 0 in y = 0, and y = L y . (10b) The solution to the problem above is (Mikhailov and Özis ¸ik, 1984;Özis ¸ik, 1993): ψ 2 (γ q , y) = sin γ q y for q = 1, 2, ... (11b) where β p = K zn K xn β p and γ q = K zn K yn γ q are the positive roots of equations sin β p L x = 0 and sin γ q L y = 0, respectively.Equations ( 11) and ( 12) are known, respectively, as the eigenfunctions and eigenvalues associated with the Sturm-Liouville problem.From the normalization integral: the corresponding norms can be obtained: 1

The transform problem
The integral transform is given by: (15)   with N (λ pq ) = N (β p ) N (γ q ).Applying the operator Lx 0 Ly 0 ψ (λ pq , x , y ) dy dx to equation (5) yields: with w n ≡ w n (x , y , z) and λ 2 pq = β 2 p + γ 2 q .Each term in equation ( 16) is solved by a different method.In the first integral, Lebniz rule is used.The second integral is solved by using Green s theorem, employing the eigenvalue problem (9), along with the boundary conditions (7).The third integral is directly substituted by definition (15).Finally, in the right-hand side, results (11) and ( 12) are used.This procedure yields in the following ordinary differential equation: where the constants are given in Appendix A.
The boundary and interface conditions associated with the transformed problem in z direction can be determined applying definition (15) to equations ( 6) and ( 8) given, respectively, by: w n = Gw 0 in z = z 0 and n = 1, (18a) w n = Gw gn in z = z i and n = N, (18b) with z = z n and n = 1, 2, ...(N − 1).

The inverse transform
From the GITT formalism, the solution to equation ( 5) is given by the expansion: The solution to equation ( 17) is given by: where A n , B n , r 1n , r 2n , w np ∈ C. The constants are given in Appendix A.
The boundary and interface conditions lead to the determination of constants A n and B n .Therefore, for each p and q, the resulting system is solved numerically.
Using the previous results in the equation ( 20), results where w n (λ pq , z) is given by the equation ( 21).Finally, we obtain the components of the average wind, u n and v n , from the fact that w n (x, y, z) = u n (x, y, z) + v n (x, y, z) i.In that case, where represents the real part of w n , and n represents the imaginary part of w n .

Eddy diffusivity
The eddy diffusivity vertical profiles employed in this study have been proposed by Degrazia et al. (2000).Such eddy diffusivities are based on Taylor s statistical diffusion theory, in which the shear buoyancy PBL spectra are modeled by means of a linear combination of the convective and mechanical forcings.Therefore, in the present case, such parameterization allowed reproducing the realistic case of a convective boundary layer where shear-generated turbulence occurs.The eddy diffusivities for such conditions are given by the following expression: where c i = α i α u (2πκ) −2/3 with α u = 0.5 ± 0.05 (Champagne et al., 1977;Sorbjan, 1989) and α i = 1, 4 3 , 4 3 for u,v and w components, respectively.w * = (u * ) 0   (Hφjstrup, 1982); L is the Monin-Obukov length and κ = 0.4 is the von Kármán constant; is the reduced frequency of the convective spectral peak, where (λ m ) i is the peak wavelength of the turbulent velocity spectra.According to Kaimal et al. (1976) and Degrazia and Anfonssi (1998) The reduced frequency of the neutral peak, (f m ) i , with

Analysis parameters
Equation ( 22) expresses, in terms of many independent parameters, the mean wind profiles.Among these parameters are: the size of the horizontal area defined by L x and L y , the thickness ∆z of the vertical sublayers over which the PBL was divided, the truncation order (p and q), and the values of large-scale divergence δ and vorticity ζ , which affect the lateral boundary conditions.
The size of the horizontal domain has an appreciable impact on the solution, but only for small areas (Figure 1).As L x and L y are successively increased from 1 km to 100 km, the solution becomes independent of the domain size for L x = L y ≥ 50 km .It means that the present solution is meaningful only over horizontal areas as large as 50 km.From this point on, the solutions shown were obtained with L x = L y = 50 km.
FIGURE 1 A similar convergence analysis was applied to the sublayer thickness ∆z and the truncation orders p and q, leading to the conclusion that ∆z = 5 m and p = q = 9 are values that warrant convergence of the mean wind profiles.The model results depend on the horizontal position within the domain, even when no large-scale divergence and vorticity are considered (Figure 2).The horizontal variation is larger close to the domain boundaries, so that there is a good portion of the domain, near its center, for which the wind profiles do not vary largely in the horizontal.The following analysis considers the vertical profiles at the domain center only.FIGURE 2

Comparison to observational data from the Wangara experiment
The Wangara experiment was conducted in Hay, New South Wales, Australia, from July to August, 1967 (Clarke et al., 1971).Wind profiles were obtained every hour up to a 2km height, using pilot balloons.A 16-m tower provided micrometeorological surface observations.In the present study, two convective days were chosen for comparison to the proposed model: days 33 and 40 (Table 1).
The wind components at the top of the domain are given by a thermal wind approximation (Equations 2), and both the surface geostrophic winds (u g0 and v g0 ) and the thermal wind magnitudes (u T and v T ) are given by Wangara observed values.
TABLE 1 The mean wind magnitudes simulated by the model are similar to the average magnitudes observed at Wangara (Figure 3).It is important to stress that such agreement concerns only the vertical overall average, but not the local maxima and minima observed at day 33, which characterize an unmixed wind profile.Indeed, such vertical variability is quite difficult to capture with a simplified model, as stated by Wyngaard (1988): " unfortunately, our knowledge of PBL physics does not yet allow us to calculate the wind profile from first principles . . .". Unmixed wind profiles, such as those observed at day 33, may be attributed to a number of reasons, such as local baroclinicity or vertical eddy diffusivity variability.Any of these reasons are, however, case-specific, and cannot be reproduced by a model where thermal wind is assumed to be constant.
The simulated values for different large-scale synoptic conditions (in terms of divergence and vorticity) cover a range of wind magnitudes, generally in agreement with the observations.Far from the surface, the condition without divergence and vorticity is the one which departs mostly from the measurements.This same condition, on the other hand, provides the best match to observations at the lowest levels (inlet).Different combinations of vorticity and divergence have been applied for the comparison.There is no clear distinction among most of them, as can be seen in Figure 3.However, an analysis based on statistical indices (Appendix B) reveals that, while the results are very similar for the wind magnitude (Table 2), the approximation for wind direction (Table 3) is improved when both the divergence and vorticity are positive.
Regarding the vertical profiles for day 33 (Figure 3), the analysis based on statistical indices shows that, when δ = ζ = 0 and δ = ζ = −f c , the model overestimates the mean observed wind magnitude (small negative values of F B). On the other hand, the statistical index F B shows that the horizontal wind direction is underestimated regardless of δ and ζ, meaning that the modeled winds are rotated counterclockwise with respect to the observations.The statistical index F S indicates that, except for the case δ = ζ = 0, the dispersion of the mean wind magnitude underestimated the experimental data.For the wind direction, this same index is negative in all cases, a consequence of the very small wind direction variability with height in the observed data, while the model results indicate a slight wind rotation with height.Other indices, such as N M SE, and F A2 are similar for all cases, and indicative of good agreement between model and observations.Finally, the correlation coefficient R was more variable, and therefore, serves as a measure of the best agreement in each case.

FIGURE 3 TABLE 2 TABLE 3
Thermal winds were observed only twice a day, at synoptic times, and those values were interpolated to 1500 LT .The large gradients near the top of the boundary layer arise from the assumed baroclinicity.For any case, the different modeled profiles agree to each other as a consequence of the top boundary conditions.They do not necessarily agree to the observed winds at the boundary layer top as a consequence of the interpolation used to calculate the thermal wind.This limitation has been noticed by Sorbjan (1989): "Finally, re-L.Buligon et al: Mean wind profile in the ABL sults of the Wangara experiment pointed out the difficulties and limitations of obtaining accurate measurements of thermal winds, vertical velocities, and representative spatially averaged fluxes".
Similar results were obtained for day 40 (Figure 4).Indeed, the statistical indices (Tables 4 and 5) indicate that the model reproduced the observations even better than for day 33.In this case, the condition without divergence and vorticity showed the largest departure from the observations for wind magnitude both at upper and lower levels.Again, the best representation offered by the model for both wind magnitude and direction occurred with positive values of divergence and vorticity.Furthermore, the solutions with negative vorticity provided the worst approximation for wind direction.A similar analysis of the statistical indices as that made for day 33 can be made for day 40 (Tables 4 and 5).

FIGURE 4 TABLE 4 TABLE 5
In both days 33 and 40, the best approximations to the observed profiles were obtained for positive vorticity and divergence, consistent with the occurrence of anti-cyclonic largescale flow at the period.Indeed, the synoptic surface pressure charts (Clarke et al., 1971) indicate the presence of a highpressure system at the region, for both days 33 and 40.Such consistency is further evidence that the model is able to reproduce the wind field realistically.

Conclusions
In the present study, a novel approach was used to obtain the average wind profile from the Navier-Stokes equations.The method is based on the Generalized Integral Transform Technique (GITT), applied to the convective boundary layer discretized in sublayers.Such discretization allows using eddy diffusivities that vary vertically.GITT is a procedure that combines series development and integral transforms, leading to a final solution (Equation 22) that contains the physical parameters determining the wind variability with height.Large-scale kinematical flow properties, such as divergence and vorticity are included in the solution, through the boundary conditions and the nonlinear advective terms of the original equations.
The results obtained are comparable to those found in the literature (Wilson and Flesch, 2004;Stull, 1988).Furthermore, the model provided a good comparison to observed data from Wangara experiment.The best approximations were obtained considering values of divergence and vorticity consistent with the synoptic charts from the experiment.
The main aim of this study is to establish an alternative method to determine the mean wind profiles.The method has been shown in detail, as well as its validation in comparison to observed data.From this point, it can be used for a further examination of a more generalized problem.As an example, in the approach taken here no temporal evolution is considered.The non-stationary problem can be solved using GITT along with Laplace transform applied to the time.Besides, the examples here were restrained to the convective case, but the development allows the use of the same approach for any stability condition.The use of appropriate eddy diffusivity profiles may lead to the determination of wind profiles under stable conditions as well.
velocity, in which α 1 = 1.7(Wyngaard et al., 1974); (u * ) 0 is the superficial friction velocity; z i is the convective PBL height; z is the height above the surface; ψ = b z i w 3 * is the adimensional dissipation rate functions, b = (0.75) the statistical indices used in this study are defined as:N M SE = (C o − C p ) o − C o C p − C p (σ o σ p )the analyzed amount and the subscript o and p refer to observed and predicted quantities, respectively, the over bar indicates an averaged value.The statistical index F B says if the predicted quantity underestimates or overestimates the average observed ones.The statistical index N M SE represents the quadratic error of the predicted quantities related to the observed ones.The statistical index F S indicates the as the model gets to simulate the dispersion of the observed data.The statistical index F A2 supply the fraction of the data (%) for the ones which 0, 5 ≤ C o C p ≤ 2. The best results are expected to have values near zero for the indices N M SE, F B and F S and near 1 in the indices R and F A2.

Fig. 1 .
Fig. 1.Simulated vertical profiles of a) wind magnitude and b) wind direction, for different domain sizes, as indicated in legend.Inlet shows the wind magnitude vertical profile for the lowest m in detail.

Fig. 2 .
Fig. 2. Simulated vertical profiles of a) wind magnitude and b) wind direction, for different positions within the domain, as indicated in legend.Inlet shows the wind magnitude vertical profile for the lowest 60 m in detail.