An Iterative Method for Calculation of Wind Profiles at the Mesoscale and Microscale

This paper presents the variational diagnostic model and iterative procedure, which enables the wind field in subdomains to be adjusted. Diagnostic models are not time dependent. Consideration of more complex features of the thermodynamic structure requires models with high resolution, which require large calculation times. The model presented applies the variational approach and enables topographical complexity of the terrain to be considered. The problem of adjusting the wind field is solved in two steps. The first step adjusts the initial wind field by means of experimental measurements or a prognosis in the larger domain, which includes smaller domains. Then the results obtained are used as the initial wind field when the grid refinement in the smaller domain is performed. This allows more precise mapping of the terrain and its architecture. Nevertheless the algorithm proposed ensures a considerable reduction in calculation time. This approach also allows us to eliminate the problem of the lack of initial data when the number of meteorological stations in the smaller domain is insufficient. The algorithm is described and validated, and numerical simulations for pollutant dispersion for a chosen town are described, followed by discussion of the iterative procedure.

numerical efficiency and accuracy, even when limited data are available, are widely used in dispersion models and research on wind energy.
An essential feature of a kinematic-diagnostic model is the lack of time-dependent equations as opposed to dynamic-prognostic models. The latter use a continuity equation and the Navier-Stokes momentum equation but also state and thermodynamic equations. Sometimes, surface warming, radiative effects, condensation, and evaporation are also considered. These models assess the temporal change of atmospheric variables and they can be used for prediction purposes. Most dynamic models require data referring to the thermodynamic structure of the atmosphere, which means that more data and calculation time is necessary than for kinematic models.
Thus, kinematic models are suitable for defining the wind profile when data necessary for dynamic models are not available. Also, when computational resources are limited, kinematic models are used. Much research (Ratto et al. 1994;Stage et al. 1999;Scire et al. 2000;Burlando et al. 2007a, b;Singh et al. 2008;Lundquist and Katopodes 2013;Gopalaswami et al. 2015) stress the fact that diagnostic models take into account measurements and terrain topography, and they should be equipped with algorithms for parametrization of the flow field. Diagnostic models are not used for predictive purposes. However, a comparison of the use of diagnostic and predictive models (Cox et al. 2003) indicates that both models allow for similar compatibility with the results of the experiment in which concentrations of pollutants were determined. Robe and Scire (1998) concluded that the kinematic model California Meteorological (CALMET) gave results compatible with those from the dynamic model MM5 (Stage et al. 1999;Lundquist and Katopodes 2013). Gopalaswami et al. (2015) note that there are models combining both kinematic and dynamic approaches (Robe and Scire 1998;Barna et al. 2000;Chandrasekar et al. 2003;Ghannam and El-Fadel 2013), and following such an approach they propose combining the MM5 model (dynamic model) with CALMET (prognostic model).
The dynamic models are used in data-assimilation approaches presented by, among others, Evensen (2009), Gao et al. (2013), Lin and Wang (2013) as well as Chandramouli et al. (2020). The data-assimilation approach, which is applied in relatively large areas for weather forecasting, combines calculations from computational fluid dynamics (CFD) with experimental fluid dynamics. The calculation process from CFD is controlled by the results of measurements. Two approaches can be encountered: segmental data assimilation with the most popular ensemble Kalman filter method (Evensen 2009;Gao et al. 2013;Lin and Wang 2013;Hamrud et al. 2015) or variational data assimilation (VDA) (Lorenc 2003;Chandramouli et al. 2020) with 3DVar and 4DVar formulations.
The diagnostic model and the procedure presented here are similar to the diagnostic version of the WindNinja model (Forthofer et al. 2014) used, inter alia, in Sanjuan et al. (2014Sanjuan et al. ( , 2015 and Wagenbrenner et al. (2016).
A similar procedure can also be found in the Quick Urban and Industrial Complex dispersion model (QUIC), specifically in the QUIC module for wind-field determination (Singh et al. 2008;Girard et al. 2018). This model allows the wind field to be estimated in urbanized areas similarly to that presented here. Because the diagnostic models approximate measurement data, they are critically dependent on the quality and density of these data.
The procedure for calculating the wind field in diagnostic models is based on adjustment of the wind field resulting from approximation by means of some corrections so that the final wind field fulfils the criterion of mass consistency. Homicz (2002) distinguishes four methods of solving the problem, each using one of the following algorithms: direct-differencing, pointiterative, hybrid, and variational calculus.
The formulation of the variational-diagnostic model used here is presented by Montero et al. (1998Montero et al. ( , 2005, Montero and Sanín (2001), Sanín and Montero (2007) and Brzozowska (2015a, b). In this formulation, the problem of minimizing the deviation of the wind field and satisfying the continuity equation for initial velocities (obtained by approximation of measurements and parametrization) involves solving the differential equation with partial derivatives. The main difference between the model presented and the WindNinja and QUIC models lies in different approaches to the solution of the differential equations. We use the equidistant finite-difference method, while the finite-element method is applied in WindNinja. Our approach enables us to use a kinematic model, avoiding the situation where the data from meteorological stations, or more precisely, the density of the stations, are too small to be fully used.
Only a short description of the method is given here together with a modification, which enables us, using the finite-difference method, to adjust the wind field into subdomains, by means of a variational method, to meteorological data with consideration of the topographical complexity of the terrain analyzed. In the first step of the algorithm, the wind field, adjusted with experimental measurements to the large domain 1 including domain 2 , is calculated. The results obtained are treated as the boundary conditions for subdomain 2 . Then, the problem of adjusting the wind field is solved again, but for the new conditions of the initial wind field in subdomain 2 . The algorithm enables us to perform grid refinement in domain 2 in comparison with the grid in 1 , and thus more precise mapping of the terrain and its architecture. The task of the adjustment of the wind field in subdomain 2 can be solved several times until the results of calculations become stable.
The main elements of the variational algorithm and description of conformal transformations are presented in Sect. 2. Section 3 presents the solution of an elliptical equation by means of the finite-difference method. The iterative procedure and an application of the algorithms developed for the solution of the task in domain 1 and subdomain 2 are described subsequently in Sects. 4 and 5. Section 6 discusses validation of the algorithm proposed and the results of simulations for a selected town. Conclusions regarding the iterative procedure of wind field calculations are presented in Sect. 7.

Variational Method, Conformal Variables
The first step is concerned with the definition of the initial wind field. The procedure used in this step is only outlined in broad terms. The details can be found in Brzozowska (2015a, b). This field, which can be obtained either by measurements, parametrization, and approximation (Girard et al. 2018;Bauer 2019;Barbano et al. 2020) or by earlier calculations, is denoted as u 0 (x 1 , x 2 , x 3 ). Usually it does not ensure mass consistency in the domain analyzed. The next step of the adjustment procedure consists in finding such a wind field u(x 1 , x 2 , x 3 ), which being close to field u 0 , fulfils the continuity equations in the following form in domain ( Fig. 1 as well as boundary conditions where n is a vector normal to the terrain surface. Using Lagrange multipliers, the problem of adjusting the wind field becomes the following minimization problem where α (1) = α (2) , α (3) are Gaussian precision moduli, φ x (1) , x (2) , x (3) is the Lagrange multiplier, and K is a set of permissible functions (Montero et al. 2005).
In order to solve the task, the following elliptic equation is obtained (Brzozowska 2015a where T (l) = 1/ 2α (l) 2 are components of the diagonal transmission tensor with the assumption that if α (1) = α (2) then T (1) = T (2) = T (1,2) Boundary conditions for (2) are in the form of Dirichlet and Neumann conditions. Three-dimensional domain , for which the solution is sought, can be reduced by means of conformal coordinates matched to the unit cube , in which the terrain surface is a horizontal plane. Such an approach simplifies the grid development process but at the same time complicates the discrete representation of Eq. 2. The procedure is presented in Fig. 2 (Brzozowska 2013;Brzozowska 2015a, b).
The following transformation of the spatial coordinates is assumed: where m are maximal spatial dimensions of the domain (length, width), x m is the maximum height of the domain (the height of the mixing layer), and s x (1) , x (2) is the

. Boundary conditions can be written in the form
Differential Eq. 4 with boundary conditions [(5a)-(5d)] can be solved using spatial discretization. To this end, the finite-element method, finite-volume method (Montero et al. 1998;Cascón and Ferragut 2007), or the finite-difference method can be used. Here, the latter method is used for the discretization (Brzozowska et al. 2009;Brzozowska 2015a, b).

Discretization of Domain Ä
According to the idea of the finite-difference method, we can assume that Eq. 4 is fulfilled if it is fulfilled at internal points of the domain , which are ξ It is also assumed that discretization points can be spread unevenly (Fig. 3).
Partial derivatives occurring in (2) can be approximated using the differences of the second order. As a result we obtain a system of m = n (1) n (2) n (3) linear algebraic equations with constant coefficients, which can be written in the matrix form as where is the vector of values of function inside the domain , φ is the vector of values of function φ at the boundaries of domain . Matrix A is a sparse matrix. Values of φ from 15 nodes of the discretization grid occur in (6).
All the equations obtained from boundary conditions using finite differences can be written as Finally, the set of equations discretized by finite differences takes the form Equation 7 forms a set of n = n (1) n (2) n (3) + 2 algebraic linear equations.

Iterative Procedure
Domain 1 with subdomain 2 ⊂ 1 analyzed is presented in Fig. 4. Let us consider a general problem of calculating the wind field in 2 when values of wind velocity are measured in a number of meteorological stations located in domain 1 . Montero et al. (1998), Brzozowska (2015b), and Brzozowska et al. (2009) present a procedure that enables the approximate values of initial wind field u 0 from (2), (3) to be calculated at nodes of the grid. The values depend on, among others, the distance of the grid nodes from measurements stations, terrain topography, the stability class of the atmosphere, and roughness of the ground. In order to capture the flow features, the parametrization is necessary (Girard et al. 2018;Oettl 2019a, b;Barbano et al. 2020). Having calculated the wind field u 0 , Eq. 4 with boundary conditions, (5a-5d) can be solved. Then the solution of algebraic linear Eq. 7 gives us the adjusted wind field u fulfilling the continuity equation at the grid nodes. The calculation time depends on the dimension of the problem. According to (7), the number of unknowns in the system is n = n (1) n (2) n (3) + 2 and the calculation time for the solution of such a system is proportional to n α where α > 1. Then it may be necessary to use a dense mesh in domain 2 ⊂ 1 , but the use of such a mesh in the whole domain 1 is not numerically effective. On the other hand, limitation of calculations to domain 2 with the use of measurements from meteorological stations situated in 2 or even in domain 1 can lead to considerable errors in calculations of the initial wind field u 0 . Thus we propose the following procedure: Stage Z(1) Calculation of wind field u 0 in domain 1 using the initial wind field obtained by approximating data from measurement stations and parametrization. The domain is discretized into n (l) + 1 intervals along axes x (l) l = 1, 2, 3. The resulting wind field is u 1 . Stage Z(2) Calculation of the wind field in domain 2, assuming that the initial wind field is obtained by interpolating wind field u 1 obtained as a solution of task Z(1). The wind field obtained in domain 2 is denoted as u 2 . Domain 2 is discretized into n (l) 2 + 1 intervals along axes x ( j) j = 1, 2 and n , which means that the number and the length of intervals along vertical axis x (3) is in task Z(2) the same as in task Z(1). Stage Z(i) i = 3, 4, . . . Calculation of wind field in domain 2 , assuming that the initial wind field in task Z(i) is the solution of task Z(i-1) and n The iterative process described in task Z(i) is carried out until the wind field in domain 2 stabilizes or other conditions for the accuracy of results are met.

Formulation of the Problem in Domain
The solution of the task formulated above in domain 1 requires an assumption of boundary conditions for ξ 2 ), the influence of the boundary conditions on the results can be considerable. Also the procedure of calculating the initial wind field u 0 with a small number of meteorological stations and a large grid can be the source of errors in the calculated wind field. Thus, let us assume that subdivision 2 ⊂ 1 is placed at a certain distance from borders ( 1 ) a of domain 1 (Fig. 5). The solution of the task in domain 1 leads us to the determination of values 1 + 1, j = 0, . . . , n 1 + 1. With these values the task of the adjustment of the wind field in domain 2 (defined above as Z(2)), can be formulated as follows. The initial conditions for the iteration process in domain 2 are assumed as where u 1 * is a vector of velocities and the function φ at all points of 2 (together with boundaries) is obtained either directly from u 1 , having solved the task in domain 1 , or from the interpolation formulae of the second order in the case when the grid in 2 is more dense than in 1 . Domain 2 is divided into n 2 +1 subintervals along directions x (1) , x (2) and it is assumed that n 1 . This means that the division in the direction of the z axis in domain 2 is the same as in domain 1 . When the grid in 2 is more dense than the one used in 1 , the respective values of the wind field u 0 can be calculated from the following interpolation formulae k+1 are nodes of the grid in domain 1 , It is important to note that the problem of adjustment of the wind field in domain 2 assumes that: (1) Discretization in the vertical direction is the same in both domains 1 and 2 . Nodes The task for 2 can be solved several times using the same grid as before and assuming The iteration process can be carried out until the wind field converges. In simulations presented in the next section, it is assumed that the wind field in 2 converges when , and ε is a relative error.
In the sub-areas, the mass continuity equation is also fulfilled. Results of numerical simulations using the iterative process described are presented and discussed in the next sections. They indicate efficiency and numerical effectiveness of the algorithm proposed as well as a considerable reduction of calculation time.

Results of Numerical Simulations
In order to validate the model of the wind field, the usual metrics are used (COST Action 710 1998; Britter and Schatzmann 2007): − the normalized-mean-square error − the proportion of calculated values within a factor-of-two of the measured values where w is the value calculated, w e is the value measured, n is the number of measurement data; − the hit rate q (VDI 2005) where D w is the normalized allowed range (assumed 0.25 for the analysis carried out), W w is the threshold value accounting for relative uncertainty of the measurements assumed as 8 × 10 -3 for the main component of the velocity vector. The lowest acceptable value of the hit rate q for prognostic models is 0.66 (VDI 2005).

Validation of the Model
Validation of the model is performed for a domain in which there are obstacles disturbing the airflow. Results of numerical simulations obtained using the authors' own model are compared with the laboratory-scale experimental measurements. The experiment was carried out in a wind tunnel in Hamburg. Geometrical conditions of the Mock Urban Setting Test (MUST) experiment (Harms et al., 2011;Bezpalcová et al. 2006;Bezpalcová 2007;Leitl et al. 2007) have been reflected on a 1:75 scale. Although the experiment was conducted in a wind tunnel, all data provided hereafter relate to the MUST test site, so the length units are not divided by 75. The authors are aware that the mere mapping of real conditions in the tunnel may cause some errors. These issues are considered in the research to which we refer. One hundred and twenty solid blocks imitating real containers (119 cuboid containers with size 12.2 × 2.42 × 2.54 m and one cuboid container of 6.09 × 2.43 × 3.51 m in the middle) were placed in the domain as shown in Fig. 6. The wind field was calculated and compared with the experimental measurements for a flow direction of 45°with respect to the position of the containers. During the experiment the flow velocity was registered for two components of the air velocity during three measurement series. The total number of measurement points was 1721 (series 1-498 points, series 2-847 points, series 3-376 points). The number of measuring profiles (planes) and their location on the x (3) axis is as follows: In order to reflect the airflow in the domain (Fig. 9), at the border of the domain (for x (1) = −140 m, x (2) = 0 and x (1) = 135 m, x (2) = 0) the input wind field was calculated for the reference speed u re f = 5.5 m s −1 (registered at x (3) re f = 7.29 m) according to Macdonald (2000) for The calculations were carried out on a discretized grid with 15 layers in the vertical direction based on trigonometric distribution (Chebyshev nodes). Such a distribution of nodes provides better convergence of the iterative algorithm for the wind field than the parabolic distribution or equal distance method (Brzozowska 2015a). Table 1 presents the results of the quantitative analysis of the main component of the flow velocity vector for each series for different density of the calculation grid in directions x (1) and x (2) . These are the results of stage Z(1) described in Sect. 4.
Analysis of the metrics presented in Table 1 shows that better results are obtained for measurement series I and II. The change of density of the grid in directions x (1) and x (2) below 1 m between the neighbouring nodes does not result in improvement of the quality of the results (reduction of error values). Domain 1 is defined by ( Below, the results of calculations of the wind field in three smaller subdomains 2 , separated from domain 1 in such a way that all points of the same measurement series belong to one subdomain, are presented. The results of calculations for domain 1 with a discretized grid of resolution 5 m are assumed as input data for subdomain 2 . The results of calculations are again compared with the results of the experiment, in this case focusing only on subdomain 2 . The results presented in the Table 2 refer to the first iteration. Even in this case, a considerable improvement in the accuracy can be seen.
The errors obtained are similar to those obtained when calculations are carried out for the whole domain. It is important to note that the application of the iterative procedure improved the metrics for series 3. The reason can be seen in the fact that subdomain 2 lies far from the edges < −150; 150 > and the initial solution for 1 with the 5-m grid gives a better approximation of the initial wind field in 2 than the initial field adopted for the domain 1 with the 1-m grid according to the parametrization.
A comparison of the calculation times for the analyzed grids for domain 1 and subdomains 2 is presented in Table 3.
Having analyzed data in Table 3, it can be seen that the calculation time necessary for the subdomains is considerably smaller. When the wind field is determined using the division of the domain into subdomains, the calculation time in subdomain 2 is several times smaller than the calculation time in domain 1 , which demonstrates the efficiency of the algorithm proposed. Table 4 shows the results for the iterative process of wind field stabilization. The wind field calculated at a given calculation step is used as the initial data for the solution of the task in the next step. The discretization grid is the same at all steps in the iterative procedure.
The results presented indicate that the conformity metrics stabilize after one or two iterations. Subsequent iterations do not significantly improve the metrics.

Simulations for a Built-Up Area in a Medium-Size Town
In this section, the algorithms presented are used in numerical simulations concerning a real area. Further analysis is concerned with part of a medium-size town (170,000 inhabitants, Bielsko-Biała, Poland) lying near mountains, which means the area has a diversified terrain. The area simulated is shown in Fig. 10.
The area presented in Fig. 10a forms domain 1 , and Fig. 10b shows subdomains 2 modelled in which P1, P2, and P3 denote three cases analyzed. It is assumed that the flow is from the north-east. Table 1 Metrics   Grid resolution in directions A logarithmic wind profile is assumed at nine points of the area marked in Fig. 10b as S 1 , . . . , S 9 (components in directions x (1) , x (2) of the wind speed at h = 10 m are u x (1) = 4 m s −1 , u x (2) = 3 m s −1 ). The coordinates of the points are as follows: Data for the elevation and hydrodynamic roughness of the terrain are obtained from GIS (geographical information system) (1:10,000 scale topographical map). The parametrization is adopted from Montero et al. (1998Montero et al. ( , 2005, Macdonald (2000), and Montero and Sanín (2001). Calculations in area 1 are carried out with the grid resolution x (1) = x (2) = 10 m (notation 10 1 ) or the grid resolution x (1) = x (2) = 50 m (notation 50 1 ). In direction x (3) , Chebyshev nodes are used with greater density towards the ground, assuming that n (3) + 1 = 16. The grid with x (1) = x (2) = 10 m is used for all subdomains in 10 2 . Sub-areas P1, P2, and P3 are denoted as 10 2,1 , 10 2,2 , and 10 2,3 respectively. Table 5 presents the results of simulations.
The results for sub-areas P j ( j = 1, 2, 3) in the second iteration (assuming input data from 50 1 ) are denoted as 10 2, j .
where u i, p , u i−1, p are absolute velocities of the wind in the i-th and i − 1 iteration and N is the total number of points (nodes of the grid in the area analyzed).
The admissible values of the metrics are presented by Chang and Hanna (2004) and the value of the average percentage error is calculated as p are velocities calculated in area (10)  1 and sub-area (10) 2, j respectively, N (2) is the number of points (nodes of the grid in area (10) 2, j ). If we assume that the input wind field u 0 2, j analyzed are presented in Fig. 11. Analysis of the results indicates that the first iteration step with the grid of 50 m in the horizontal plane already provides correct results (with respect to 10 m mesh in the whole area).
The calculation time of wind fields in each sub-area is considerably shorter. The average percentage error for wind fields calculated in area 1 and sub-areas 2,1 , 2,2 , 2,3 with the same mesh does not exceed 1%.

Conclusions
The method of subdomains enables us to considerably shorten the calculation time. Validation of the model presented in Sect. 6.1 concerns a standard experiment widely used for the validation of wind field models for built-up areas (with complicated topography). The analysis indicates that the results already converge after the first iteration of the algorithm. The method applied in areas with complex topography makes it possible to obtain more accurate results with a considerably reduced time of calculations. It can be successfully used for typical meteorological conditions in a given area/town (using the database of the wind field) when the wind field and pollutant concentration have to be defined in the microscale. The calculations presented in Sect. 6.2 used a mesh not smaller than 10 m, which is the typical grid size considered in a mesoscale. However, the model, as shown in Brzozowska et al. (2007) and Brzozowska (2016), also enables modelling of the wind field and dispersion of pollutants in the microscale with good accuracy.
In general, it should be noted that dynamic, prognostic models of the CFD class controlled by measurements have the greatest accuracy-as in the data-assimilation approach. However, they require detailed data and high computing power for performing calculations. Diagnostic, stationary models that fulfil the continuity equation, sometimes equivalent to the momentum equation, are less accurate, as in the NinjaWind or QUIC-URB packages.
An interesting comparison of diagnostic and CFD class models (Reynolds-averaged Navier−Stokes and large-eddy simulation) is presented by Hayati et al. (2019). It is shown that diagnostic models (as the model presented in this paper) are definitely cheaper in terms of costs and calculation times. The use of parallel computations can significantly improve computational efficiency. Here, parallel computations were used by means of multithreading. As demonstrated by Bozorgmehr et al. (2021), an additional improvement in the efficiency of computations can be obtained by applying graphic processing units.
The model presented can be used successfully when the prognostic model is not applicable (too much input data is required), when there are no meteorological stations in the analyzed area, or the grid of available data points from meteorological stations is so sparse that it does not cover the area analyzed. The use of the method proposed allows for the initial calculations of the air velocity field and then the calculation of the wind field in any sub-area, at the same time refining the discretization grid. If, additionally, a database of wind fields is created for typical meteorological conditions in a given region, the primary calculations can be avoided each time, and only performed for selected sub-areas. This considerably reduces the time needed to perform simulations, for example, in the case of sudden releases of hazardous substances or episodes related to low emission.