Mass–Consistent Wind Field Models: Numerical Techniques by L2–Projection Methods

For several meteorological problems and a large number of applications, the knowledge of the 3–D wind field over a region is required. Examples include prediction of the transport, diffusion and dispersion of air pollutants in the atmosphere (Finardi et al., 2010; Sherman, 1978), realization of wind maps for the design of different urban and general projects (Castino et al., 2003), and the effect of wind on structures and fire spreading (Potter & Butler, 2009), among others. Moreover, meteorological wind fields are also required inputs for air quality models. In practice, usually limited horizontal wind fieldmeasurements are available, and therefore the calculation of the vertical motion must be predicted or calculated. Several methods and strategies, with various levels of complexity, have been proposed to address this problem. They can be included into two general model types: prognostic models and diagnostic models. Prognostic models are complex time–dependent hydrodynamic models governing air flow, including thermal effects, density variation and turbulent interaction. While these models are “realistic”, they are expensive to operate, need extensive computer facilities, and require specialized training for their operation. On the other hand, diagnostic wind models do not require the integration of the non–linear hydrodynamic equations. Instead, available interpolated data is used to generate wind fields, which satisfy some physical or dynamical constraints. For instance, to assure mass conservation, a simplified steady–state version of the continuity equation is imposed, and the resulting model is then called a mass–consistent model. A review of these models is available in Ratto et al. (1994) and Ratto (1996). We focus in a variational mass–consistent model which is based in the original formulation by Sasaki (Sasaki, 1958). This approach has been used for a variety of meteorological problems (Castino et al., 2003; Pennel, 1983; Sherman, 1978; Wang et al., 2005). Mass–consistent models are attractive because of their simplicity, and because they are easy and economical to operate. In some applications, these models outperform the more sophisticated and expensive dynamical models (Ratto et al., 1994). However, mass–consistent models have some disadvantages, because they are based on incomplete or idealized models and have difficulty 2


Introduction
For several meteorological problems and a large number of applications, the knowledge of the 3-D wind field over a region is required. Examples include prediction of the transport, diffusion and dispersion of air pollutants in the atmosphere (Finardi et al., 2010;Sherman, 1978), realization of wind maps for the design of different urban and general projects (Castino et al., 2003), and the effect of wind on structures and fire spreading (Potter & Butler, 2009), among others. Moreover, meteorological wind fields are also required inputs for air quality models. In practice, usually limited horizontal wind field measurements are available, and therefore the calculation of the vertical motion must be predicted or calculated. Several methods and strategies, with various levels of complexity, have been proposed to address this problem. They can be included into two general model types: prognostic models and diagnostic models. Prognostic models are complex time-dependent hydrodynamic models governing air flow, including thermal effects, density variation and turbulent interaction. While these models are "realistic", they are expensive to operate, need extensive computer facilities, and require specialized training for their operation. On the other hand, diagnostic wind models do not require the integration of the non-linear hydrodynamic equations. Instead, available interpolated data is used to generate wind fields, which satisfy some physical or dynamical constraints. For instance, to assure mass conservation, a simplified steady-state version of the continuity equation is imposed, and the resulting model is then called a mass-consistent model. A review of these models is available in Ratto et al. (1994) and Ratto (1996). We focus in a variational mass-consistent model which is based in the original formulation by Sasaki (Sasaki, 1958). This approach has been used for a variety of meteorological problems (Castino et al., 2003;Pennel, 1983;Sherman, 1978;Wang et al., 2005). Mass-consistent models are attractive because of their simplicity, and because they are easy and economical to operate. In some applications, these models outperform the more sophisticated and expensive dynamical models (Ratto et al., 1994).
However, mass-consistent models have some disadvantages, because they are based on incomplete or idealized models and have difficulty representing flows accurately in data-sparse regions as mountains or oceans. Despite these limitations, mass-consistent models are a valuable tool for air quality applications and consequently several developments have taken place over last decades (Ferragut et al., 2010;Ratto, 1996;Ratto et al., 1994;Ross et al., 1988;Wang et al., 2005). Most of the results presented in this chapter has been published in the last few years (Flores et al., 2010;Núñez et al., 2007;2006), but we also include some additional ideas and recent results. The variational method proposed by Sasaki uses the continuity equation ∇·u = 0, where u is the wind velocity vector field on a given domain Ω. The method is based on the minimization of the functional L defined by where u I is an initial observed wind field, λ is a Lagrange multiplier and S is a diagonal matrix with weighting parameters α i > 0, i = 1, 2, 3, called Gaussian precision moduli, related to the scales of the respective components of the velocity field. The vertical component of the initial wind field is taken as zero because meteorological stations usually do not measure this component. The Euler-Lagrange equations of (1) are: Usually u is obtained from (2), after λ is computed. Since ∇·u = 0, then from (2) we obtain the elliptic equation −∇ · S −1 ∇λ = ∇·u I ,f r o m w h i c h λ is obtained. To complement (close) this equation, two types of boundary conditions are commonly used: homogeneous Dirichlet boundary conditions, λ = 0, for open or "flow through" boundaries (like truncated boundaries), and Neumann boundary conditions, ∂λ/∂n = 0,forclosedor"noflowthrough" boundaries (like the surface terrain or topography). Many authors have been used and recommend these boundary conditions (Kitada et al., 1986;1983;Ratto et al., 1994;Sherman, 1978). However they are physically and mathematically inconsistent as we will show in this work. Even though, there have been several sophisticated developments in the numerical simulations of this model as, for instance, the application of multigrid methods (Wang et al., 2005), and the application of genetic algorithms to estimate parameters (Montero et al., 2005), it seems that the analysis of boundary conditions has not attracted the attention of the community in meteorology.
In this work we study how boundary conditions affect solutions of the elliptic equation for λ. We show that the application of incorrect boundary conditions may degrade the solutions several orders of magnitude, and we propose some strategies to overcome this problem. In particular, we introduce a new approach based on the saddle-point formulation of the constrained least squares formulation of the problem, which allows the introduction of successful techniques from computational fluid dynamics. This new approach does not require boundary conditions for the multiplier. It produces much better results, and it also helps us to establish more consistent boundary conditions on truncated nonphysical boundaries. We also explore other boundary conditions for the multiplier better suited for artificial truncated boundaries. Furthermore, we present some preliminary numerical results using a meshfree method based on a radial basis function collocation method.

Mathematical formulation of the problem
Let Ω be an open, simply connected and bounded region in R d (d = 2o r3 ) with Lipschitz boundary ∂Ω = Γ N ∪ Γ D ,w h e r eΓ N is the part of the boundary associated to the surface terrain (topography), and Γ D is the rest of the boundary (artificial vertical boundaries and top boundary), as shown in Figure 1. Given an initial vector field u I in Ω (which can be obtained by interpolating atmospheric data, or by other means), our goal is to generate a solenoidal field u -called adjusted field-as close to u I as possible in a sense that will be clarified below, such that u · n = 0onΓ N . We define the following vector function spaces: L 2 (Ω)=L 2 (Ω) d and H(div; Ω)= { v ∈ L 2 (Ω) : ∇·v ∈ L 2 (Ω) }. Then, the adjusted wind field u must belong to the normed closed space with the norm · S,Ω associated to the inner product u, v S = Ω Su · v dx,w h e r ev · w = ∑ d 1 v i w i is the usual scalar product in R d . We can now formulate the problem as a least squares projection problem. For this purpose, we define a convex quadratic functional J : V → R as Then, our problem can be stated as follows: Due to the properties of this functional, u ∈ V is a minimizer of J if and only if it is a stationary point of J: Will-be-set-by-IN-TECH

Derivation of the elliptic problem
The first approach is based on a Helmholtz-type decomposition of the Hilbert vector space L 2 (Ω), and it reduces to the traditional approach used by meteorologists. Proposition 1 The orthogonal complement in L 2 (Ω) of the closed subspace V is An argument very similar to that given by Girault and Raviart (Girault & Raviart, 1986), shows that this decomposition is valid (details are given in (Núñez et al., 2007)). Therefore we get from (6) that S u − u I = ∇λ,withλ in With the above properties, we obtain a saddle-point problem for u and λ (left), as well as the correspondent elliptic problem for λ (right): Su −∇λ = Su I ,a n d∇·u = 0i n Ω, −∇ · S −1 ∇λ = ∇·u I in Ω,( 8 ) To obtain the elliptic problem, we eliminate u from the saddle-point problem using that u belongs to V.O n c eλ is calculated from (8)-(10), the adjusted field is recovered from (2). Equation (8) has traditionally been used by meteorologists. However, this equation is generally introduced from a discussion in which it is not clear how to establish the proper boundary conditions for λ. The crucial argument in our study is the decomposition of L 2 (Ω) in orthogonal subspaces V and V ⊥ , from which the boundary conditions for λ arises in a natural way, from the mathematical point of view. We would like to mention that the boundary condition (10) has already been used in recent research (Ferragut et al., 2010).

Finite element solution of the elliptic problem
The variational formulation of the elliptic problem (8)- (10) is Here, we consider the two-dimensional case. Let T h be a finite element triangulation of Ω ⊂ R 2 (Ciarlet, 2002), where h is taken as the space discretization step. Let's denote by P 1 the space of polynomials of degree less or equal than 1. Then, L 2 (Ω) and H 1 D (Ω) are approximated by the finite dimensional spaces respectively. Thus, the finite element algorithm is: where u I h ∈ L h is the interpolant of the given initial velocity field u I . We obtain λ h after solving the resulting system of linear equations, and the numerical approximation u h of u is computed by the weak version of (2) as follows: From now on, we identify the algorithm (14)-(15) as the E1-algorithm. Example 1. We consider the solenoidal vector field u( Assuming that we have u I (x, z)=(x,0 ) as an initial horizontal wind field, we want to apply the E1-algorithm to see how much we can recover of the vertical component of u. For this numerical calculation, Ω is divided into a 80 × 80 triangular mesh, and we choose the following values for the Gaussian Precision moduli: α 1 = 1andα 3 = 0.001. Figure 2 shows the exact field in red and the computed adjusted field in blue. Both fields agree fairly well almost everywhere, except on the vertical artificial boundaries x = 1andx = 2.  Fig. 2. Exact field u =(x, −z) in red, adjusted field obtained by the E1-algorithm in blue.
The relative error and the mean divergence of the computed solution are defined as e r = ||u − u h || 2 ||u|| 2 ,a n dmdiv = mean respectively. The point-wise divergence is computed in a weak sense, as follows where φ i is the piece-wise linear base function associated to vertex node x i .F o rt h ep r e s e n t example, we obtain e r = 1.9 × 10 −2 and mdiv = 4.1 × 10 −2 . The values for the Gaussian precision moduli were chosen based on numerical performance. Table 1 shows the behavior of e r and mdiv for different values of α 3 when α 1 is kept constant and equal to one. Clearly the best results were obtained with α 3 = 0.001. We will explain this behavior later on. But, for the moment, we want to emphasize that this algorithm produces satisfactory results almost everywhere, except on the boundary Γ D , where homogeneous Dirchlet boundary conditions were imposed.
On the other hand, one of its major drawbacks is that inconsistent or incorrect boundary conditions, on truncated artificial boundaries, degrade the accuracy of the solution. In the rest of the chapter, we introduce some alternatives to overcome these problems.

A saddle-point formulation and the conjugate gradient algorithm 4.1 Derivation of the formulation
The second approach to solve the problem (5), or equivalently problem (6), is based on the usual methodology to solve constrained optimization problems. That is, we introduce the space of vector functions together with the Lagrangian L defined on V N × L 2 (Ω) as A stationary point (u, λ) of L solves the following saddle-point problem where λ need not satisfy boundary conditions.Thesolutionu is the minimizer of J, and now it is obtained from the enlarged space V N where free divergence is not required. Instead, the condition ∇·u = 0 is relaxed by the introduction of the Lagrange multiplier λ so that u must satisfy the weaker condition (20). To solve (19)- (20) we introduce a method which has shown to be very effective for solving Stokes problems in computational fluid dynamics (Glowinski, 2003). The idea is as follows: assuming that (u, λ) is a solution of the problem (19)-(20), the vector field u is decomposed as u = u I + u λ ,w h e r eu I is the given initial vector field, and Furthermore, u λ must satisfy (20) which has the following equivalent strong version A key point is that problem (21)-(22) can be formulated as a functional equation. For this we introduce the linear operator A from L 2 (Ω) into L 2 (Ω) defined by where With this definition, it is clear, from (21)-(22), that the multiplier λ satisfies the functional equation

Conjugate gradient algorithm
Operator A is selfadjoint, and strongly elliptic, since from (23) and (24) we have Therefore, the following iterative conjugate gradient algorithm may be used to solve the infinite dimensional problem (25): 3. If g k , g k ≤ε g 0 , g 0 , take λ = λ k+1 and u = u k+1 and stop. Otherwise, compute

29
Mass-Consistent Wind Field Models: Numerical Techniques by L2-Projection Methods www.intechopen.com Above, ·, · indicates the usual scalar product in L 2 (Ω). Observe that the adjusted field u is also computed as an intermediate step in the algorithm. In this algorithm, no boundary conditions are imposed on λ, contrary to what it was done in the first approach. This fact has a very important effect in the numerical calculation.

A mixed finite element method
To approximate the functions in V N and L 2 (Ω), considered in the previous algorithm, we use the Bercovier-Pironneau finite element approximation (Bercovier & Pironneau, 1979). Functions in L 2 (Ω) are approximated by continuous piecewise linear polynomials over a triangulation T h of Ω, while the elements in V N are also approximated by linear polynomials but now over a twice finer triangulation T h/2 of Ω. The fine triangulation T h/2 is obtained from a regular subdivision of each triangle T ∈T h . Then, the functional spaces V N and L 2 (Ω) will be approximated, respectively, by the finite dimensional spaces We apply this mixed method, particularly to solve the integral equations in steps 1 and 2, as well as for the calculation of the weak divergence to obtain g 0 in step 1 and w k in step 2. Those calculations require this mixed method, or any other stable finite element pair, to avoid instabilities in the numerical solution. Actually, the main cost of this algorithm is the solution at each iteration of the integral equation to get u k and the calculation of w k . However, if the trapezoidal rule is applied to approximate the left hand side of the integral equations, we obtain a system of algebraic equations with diagonal matrix, and the cost to solve them is just a vector multiplication. We call this new algorithm the CG-algorithm. Example 2. We consider again the initial horizontal field u I =( x,0),a si nE x a m p l e1t o test the performance of the CG-algorithm. In order to compare the numerical results with those obtained with the E1-algorithm,w ec h o s eh = 1/40 and h/2 = 1/80 in this case.
To stop the iterations we choose ε = 10 −8 at step 3. Figure 3 shows the exact and the adjusted wind fields. The agreement is excellent this time, even at the vertical boundaries x = 1a n dx = 2. The relative error and the average divergence are e r = 5.9 × 10 −4 and mdiv = −5.3 × 10 −12 , respectively. Note that we got a significant improvement: nearly two orders of magnitude better on the relative error, and about ten orders of magnitude better on the average divergence. The improvement of the relative error is mainly due to the reduction of the error on truncated boundaries, while the enhancing of average divergence is mainly due to the iterative method, because it stops when it reaches the tolerance (i.e. when the norm of the divergence is small enough).
To test further the CG-algorithm we consider two, more "realistic", additional examples. The first one includes a domain with a topography of a cosine-shape, and the second one includes a domain with a real topography. In both cases, the "exact" wind field was obtained with a Stokes solver using the methodology described in (Glowinski, 2003). The initial wind field u I was obtained dropping the vertical component of the "exact" one in both cases. Then, the vector wind field is recovered using the same discretization parameters as in example 2.  Fig. 3. Exact field u =(x, −z) in red, adjusted field obtained by the CG-algorithm in blue.
The "exact" wind field satisfies ∇·u = 1.2 × 10 −16 . Figure 4 shows the adjusted and "exact" wind fields. Example 4. Terrain elevation from real data. In this case, the domain is defined as where h(x) is a function constructed via cubic splines, which interpolate discrete data over 10 Km of real topography of a certain region in Mexico, contained in a database (GTOPO, 1997). The "exact" wind field satisfies ∇·u = 6.1 × 10 −16 . Figure 5 shows the adjusted and "exact" wind fields. We have an excellent agreement in all cases, even on truncated artificial boundaries. The relative error and the computed mean divergence are about the same order as in example 2.

Preconditioned conjugate gradient method
The CPU time to solve the problem with the CG-algorithm, at the level of accuracy shown in Table 2, is about twice the CPU time needed to solve the problem with the E1-algorithm.
In order to make this algorithm more reliable we need to speed up the iterative algorithm to get at least a comparable computational efficiency. Fortunately, we have found a good preconditioner for the iterative algorithm. This preconditioner is an optimal one, and we are presently working in its computer implementation, so we only describe here the main ideas without presenting numerical results yet. Let B : L 2 (Ω) → L 2 (Ω) be an operator defined by Bq = φ q ,w h e r eφ q solves : Operator B is self-adjoint and elliptic, and satisfies A (Bq)=q,insideΩ, for every q ∈ L 2 (Ω). An easy way to see these properties is considering the differential form of operators A and B: Bq= φ q = −[∇·(S −1 ∇)] −1 q,s i n c e −∇·(S −1 ∇φ q )=q in Ω.
Then, from (27)-(28) we obtain This shows that B can be used as an optimal preconditioner. Therefore, the additional cost of the preconditioned conjugate gradient algorithm is the solution of an elliptic problem at each iteration. However, this additional cost is offsetted by two nice properties: a) the preconditioning must reduce drastically the number of iterations (from about 1000 to less than 20, based on previous experience in CFD); b) there is a significant reduction of degrees of freedom in the discrete version of the elliptic problem associated to operator B. This elliptic problem is solved in a coarse mesh, and it is four times smaller than the elliptic problem for the multiplier λ in the 2-D case, and about eight times smaller in 3-D problems.

Some extensions and future research
In this section, we present some additional alternatives to look at the problem. We first consider a different set of boundary conditions on vertical truncated boundaries for the multiplier λ, and we show that it produces better results than the traditional approach. We also show that if we introduce ghost nodes on the truncated artificial boundaries, we get even a better improvement. Finally, we introduce radial basis functions to solve the elliptic problems for the multiplier, and show that this is a promising alternative for 3-D wind fields.

Alternative boundary conditions for the elliptic problem
From equations (19)- (20), we obtain The boundary integral in (30) vanishes in two cases, namely: when v · n = 0orwhenλ = 0 on Γ Γ N . The first case is not possible since it holds only on Γ N , and the second case is not a good choice on vertical boundaries as we have already seen in Section 3. However, there is a possibility: decompose Γ D as the union of the vertical boundaries, Γ V ,andthetopboundary , Γ T .N o w ,a tΓ T we still impose λ = 0, and on Γ V we impose the new boundary condition u · n = u I · n. This new boundary condition is reasonable, since we assume that u I is the horizontal part of u. Therefore, with this choice, we obtain the saddle-point problem (left) and its corresponding elliptic problem (right): Su −∇λ = Su I ,a n d∇·u = 0i n Ω, u · n = 0o n Γ N . −S −1 ∇λ · n = u I · n on Γ N . 12

Will-be-set-by-IN-TECH
The finite element algorithm for the elliptic problem is: where H h is defined as in (13), but with q = 0o nΓ T instead of Γ D . Equation (36) differs from equation (14) only by the boundary integral on Γ V . We call (36), together with (15), the E2-algorithm. Example 5. Let us consider one more time the problem introduced in example 1 and in example 3, with the same discretization parameter, h = 1/80. The recovered wind field for both cases is better than the one obtained with the E1-algorithm, since the vertical component is recovered fairly well, not only in the interior of the domain but also at the vertical boundaries. Figure 6 shows the exact and recovered wind fields. Table 3 shows a summary of the results obtained in examples 1, 2, 3 and 5. For the case with exact wind field u =(x, −z) the immediate effect of this improvement is the reduction of the relative error by two orders of magnitude. However, we do not obtain a comparable reduction of the mean divergence. For the problems with cosine topography occurs the opposite. Actually, the numerical results show that the most effective algorithm to reduce both, the relative error and the average divergence is the CG-algorithm. 3 cosine topography CG-algorithm 3.6 × 10 −6 9.8 × 10 −9 955 3.30 5 cosine topography E2-algorithm 9.2 × 10 −2 3.6 × 10 −4 -2 . 0 8 Table 3. Summary of the numerical results obtained in examples 1, 2, 3 and 5.

Ghost nodes
Given that u belongs to V and satisfies (2), then λ satisfies the equations −∇ · (S −1 ∇λ)=∇·u I ,i nΩ, −(S −1 ∇λ) · n = u I · n,o nΓ N . Instead of looking for boundary conditions on Γ D , we may enforce mass conservation by asking any solution of (37)- (38) to satisfy Equations (37)-(39) imply the identity Ω ∇·u I dx = Γ u I · n dΓ. Actually, this is the compatibility condition associated to the the above Poisson-Neumann-like problem. Therefore, this problem has a unique solution λ ∈ H 1 (Ω)/R, and its variational formulation is: Given u I ∈ L 2 (Ω),findλ ∈ H 1 (Ω)/R such that Observe that when q = 1inH 1 (Ω)/R, we recover (39). However, the computational solution of this problem is not trivial, since the matrix associated to the discrete version is semidefinite.
On the other hand, the symmetry of the matrix is lost because the boundary integral in the right-hand side has the unknown λ. A way to overcome this computational problem is to introduce "ghost nodes", around and beyond of the nonphysical truncated boundary Γ D . Then, we may impose λ h = 0a n d / o r∂λ h /∂n = 0 on the outer layer of those ghost nodes. At the end, we discard the solution on the ghost nodes, and we only keep the solution values on the actual nodes. Actually, this is a well-known way to deal with differential equations in domains with truncated boundaries. Example 6. We consider one more time the problem from example 1 with the same discretization parameters. We incorporate two layers of ghost nodes and impose λ = 0o n the outer layer. The recovered wind field obtained is such that the relative error and average weak divergence are e r = 2.1 × 10 −5 and mdiv = 1.6 × 10 −6 , respectively. The figure with the comparison of the adjusted wind field and the exact wind field is not shown, because it is very similar to Figure 5. Instead, we summarize in Table 4 the results for this example with the different algorithms.
Case E1-algortihm E2-algortithm Ghost-Nodes CG-algorithm e r 1.9 × 10 −2 4.0 × 10 −4 2.1 × 10 −5 5.4 × 10 −4 mdiv 4.1 × 10 −2 1.8 × 10 −2 1.6 × 10 −6 −5.2 × 10 −12 Table 4. Comparison of numerical solutions obtained with different algorithms. Table 4 shows how boundary conditions degrade numerical calculations. It is observed that the solution improves each time the Dirichlet boundary condition λ = 0 is applied to a smaller section of the non-physical boundary. This is not surprising, since this boundary condition introduces a large artificial gradient, mainly on vertical truncated boundaries, when calculating the term ∇λ in order to get u = u I + S −1 ∇λ at the corresponding boundary nodes. At this time, and taking in account the performance of every algorithm, we may recommend to use either the classical approach with ghost nodes or the saddle point problem approach with the conjugate gradient algorithm, specially if we do not have enough information at truncated boundaries.
14 Will-be-set-by-IN-TECH

Approximation with radial basis functions
Af u n c t i o nΦ : R d → R is called a radial basis function (RBF) if Φ(x)=φ( x ) where the kernel φ is a scalar function φ : R + → R,andd the spatial dimension. Usually x is denoted by r, and typical functions used in applications are, among others: 1. Multiquadrics, φ(r)= √ c 2 + r 2 .
The constant value c es called the shape parameter. The radial basis function method was first introduced in the 1970s for multivariate scattered data approximation, (Hardy, 1971). This interpolation problem is defined as follows: Given a set of points {x j } n j=1 ⊂ Ω ⊂ R d , approximate the function f (x) from the set of values where p is a polynomial which depends on the specific RBF. Then the interpolation condition gives an algebraic system of equations for λ = {λ i } k j=1 . However, the corresponding matrix could be ill conditioned and, in some cases, even rank deficient, and special techniques are needed, like preconditioning and least squares (Buhmann, 2003), (Wendland, 2005). In the last two decades, the main focus of the applications seems to have slowly shifted from scattered data approximation to the numerical solution of PDE. Radial basis function collocation methods for solving PDE are truly meshfree algorithms, in the sense that collocation points can be chosen freely and no connectivity between the points is needed or used (Kansa, 1990), (Narcowich & Ward, 1994). The main attraction of RBF collocation method to solve PDE is that it can be extended directly to solve 3-D problems. Moreover, due to the absence of a grid, these techniques are better suited than classical methods to cope with problems having complex boundaries. So, we think that the RBF collocation method is a good choice to study our problem, and we want to explore this alternative. Let us denote by L the linear elliptic differential operator for the multiplier λ,a n dB the boundary operator. Suppose that we want to solve the problem We consider a set of collocation points {x i } n i=1 ,w i t hn i points in the interior and n b points on the boundary, so that n = n i + n f . We look for an approximate solution λ h (x)= Therefore the recovered wind field u h is given by Example 8. As a last example we include a 3-D numerical calculation using radial basis functions. The exact syntectic wind field for this example is u =( x, y, −2z). We dropped the vertical component so that u I =( x, y,0). A multiquadric kernel with c = 11.33 was used, and we chose α 1 = α 2 = α 3 = 1. The collocation points were obtained by a 5 × 5 × 5 regular subdivision of Ω =( 1, 2) × (0, 1) × (0, 1). Figure 7 shows the collocation points and the comparison between the exact and recovered wind field. The agreement is excellent, we obtained a relative error e r = 1.98 × 10 −4 and mean divergence mdiv = −5.59 × 10 −6 .

Conclusions
We studied the problem of generating an adjusted wind field from horizontal wind data by different numerical techniques. We have shown that boundary conditions can significantly affect numerical solutions depending of how we treat artificial truncated boundaries. The usual methodology (E1-algorithm) does not produce satisfactory results close to the vertical boundaries due to the high gradients introduced by the term S −1 ∇λ in (2), when homogeneous Dirichlet boundary conditions are imposed there. The formulation of the problem as a saddle-point one, with a functional equation that has a self-adjoint and strongly elliptic operator, allows the use of an iterative conjugate gradient algorithm (the CG-algorithm). This new methodology, in the context of mass consistent models, produces much better results, it does not involve the solution of differential equations, and it does not require boundary conditions for the multiplier. An optimal preconditioner for the conjugate gradient algorithm was introduced, and we hope, based on previous experience with the solution of the Stokes problem, to reduce the number of iterations of nearly a thousand to afewtens.
In an attempt to improve the numerical results obtained with the traditional approach, we introduced new boundary conditions for the multiplier on vertical boundaries. These boundary conditions are demonstrated to be physically and mathematically consistent. The numerical reconstruction of the wind field was improved by two orders of magnitude, but a comparable reduction on the weak divergence is not always obtained. However, the introduction of "ghost nodes" produces more satisfactory results, reducing both the relative error and the mean weak divergence by two or more orders of magnitude. On the other hand, we have just started to explore meshfree methods. In particular, radial basis collocation methods seem to be a very simple reliably alternative to the reconstruction of three-dimensional wind fields, according to the preliminary numerical results shown in this work.
The application of the different alternatives and methodologies presented here to the more realistic three-dimensional case is a continuation of the present work. Another interesting issue is the potential extension and application of these methodologies to other experimental fields, such as fluid dynamics and computer vision. In particular, the reconstruction of solenoidal velocity fields from experimental data, obtained through the particle image velocimetry technique, is an important issue, (Adrian, 2005). Its relation with computer vision is established by optical flow estimation, (Ruhnau & Schnorr, 2007). The content of this book covers several up-to-date topics in fluid dynamics, computational modeling and its applications, and it is intended to serve as a general reference for scientists, engineers, and graduate students. The book is comprised of 30 chapters divided into 5 parts, which include: winds, building and risk prevention; multiphase flow, structures and gases; heat transfer, combustion and energy; medical and biomechanical applications; and other important themes. This book also provides a comprehensive overview of computational fluid dynamics and applications, without excluding experimental and theoretical aspects.

How to reference
In order to correctly reference this scholarly work, feel free to copy and paste the following: wind-field-models-numerical-techniques-by-l2-projection-methods