On Application of Artificial Neural Network Methods in Large-eddy Simulations with Unresolved Urban Surfaces

Micro-meteorological aspects of city comfort, land use management and air quality monitoring are rapidly growing areas of applications where environmental turbulence-resolving or large-eddy simulation (LES) models play the central role. Complex details of the urban surface morphology however remain unresolved or poorly resolved in the state-of-the-art LES due to severe limitations from computer facilities. The LES code LESNIC is applied in this study to simulate turbulent flow interaction with elements of the urban surface morphology. The study investigates a possibility to utilize a trained three layers’ artificial neural network (ANN) to parameterize the flow-to-surface interactions in the coarse LES where surface features are unresolved. It is concluded that the ANN can be a robust predictor for scalar concentrations and components of the surface stress tensor in the urban sub-layer with unresolved scalar sources and surface morphology. It has been noted however that dynamic training of the ANN may require more computational resources than adequate refinement of the LES resolution.


Introduction
Meteorological applications, e.g.city comfort, land use management and air quality monitoring, require development of turbulence-resolving large-eddy simulation (LES) models as well as approaches to utilize them in simulations of turbulent flows over poorly resolved complex surface features.Figure 1 gives an example of such an urban surface obtained from the MEGAPOLI database of the Paris surface morphology (Sievinen et al., 2009).Although the digital elevation model of Paris is sufficiently smooth and large-scale as it can be easily noticed in Figure 1, the building elevations on the top of it cannot be resolved with a moderately fine mesh in LES as the buildings have a lot of the elevation variability on scales 100 m and less.Elsewhere in the world, more complicated surface features may be required to incorporate in the model (e.g.Hj Kamaruzaman & Mohd Hasmadi, 2009).The under-resolved surface structures can be incorporated in the model boundary conditions in three ways: traditional parameterization methodologies (e.g.Britter and Hanna, 2003;Bou-Zeid et al., 2009); adjoint optimisation model (e.g.Losch and Wunsch, 2003); and artificial neural network technique (ANN).The latter approach is addressed in this study.We will investigate applicability and utility of the three layers' ANN in a priori numerical experiments.ANN is a statistical technique, which solves a set of parametric function equations, called static nonlinear maps, with empirically tuned coefficients or weights to relate the model variables to unresolved quantities of interest.The weights are obtained in so-called training process searching for the global minimum of a cost function in multi-dimensional parameter space.ANN were found useful in many closely related applications, e.g.turbulence control (Lee et al., 1997;Ayhan et al., 2004), turbulent flow simulations (Panigrahi, et al., 2003;Hocevar et al., 2004), air quality monitoring (Morabito and Versaci, 2002;Mammarella et al., 2006).In particular, Kukkonen et al. (2003) evaluated performance of five ANN methods, a linear regression statistical model and a hydrodynamic model using measurements of urban chemical concentrations at two stations in central Helsinki from 1996 to 1999.On average, the ANN methods demonstrated better agreement with data than other types of models.Mammarella et al. (2006) reported the use of the ANN in the fully operational air quality monitoring system "ATMOSFERA".
In the past, limited computational resources restricted the LES applications to the case of turbulent flows over homogeneous flat surfaces.More interesting applications like urban "hot-spots" of anthropogenic impact on the Earth's system were not fully addressed.In this study, the LES code LESNIC developed at the Nansen Environmental and Remote Sensing Center (Esau, 2004) will be utilized with the ANN application to elements of an urban surface.

Artificial neural network in modeling of surface stress
Limitations on computational resources and a lack of knowledge to parameterize sub-grid scale processes in coarsely resolved simulations have motivated a number of attempts to make shortcuts between practically important but unresolved quantities, like concentrations in the urban sub-layer, and more readily simulated ones, like the mean wind.The ANN application approach exploits persistent dynamical patterns on resolved scales created by the unresolved interactions between fixed surface features and the turbulent flow.Suppose a certain feature of boundary conditions causes a perturbation in modelled variables.For instance, a building, unresolved with the given model mesh, perturbs the flow.In this case, one can introduce a statistical correction included as an additional artificial force into the model.Since the exact structure of this force is unknown, it is approximated with a matrix of stochastic empirical coefficients and a set of linear equations that is in fact the ANN.The coefficients are chosen in such a way to minimize the difference between the model and the observations.Mathematically, the ANN can be considered as an optimal estimator of unresolved quantities using a finite time series of resolved-scale solutions (Moreau et al., 2006).In this study, a three-level ANN (one input layer of 16 neurons and two hidden layers of 4 an 1 neurons) with a back-propagation learning function and the tan-sigmoid transfer function is considered.The neuron distribution was selected optimal by trail-and-fail approach for the LES data.The ANN could be written as Here ) ( i x q is a predicted quantity, e.g. a chemical concentration or velocity stress, given at positions i is the input vector of predictors given at positions k x , which could be different from i x .Other symbols in Eq. ( 1) denote empirical weights j A , jk B , mi C and thresholds jm b and m a .The weights and thresholds are to be "trained", i.e. empirically fitted using one or another statistical approximation, on a sub-set of available data . Thus, the correspondence between i x and k x sets of positions will be implicitly built-in into the weights during this training process.This ANN feature will be used in its application to urban LES studies where the training process will implicitly account for the correspondence between the fixed positions of the unresolved surface features and the fixed positions of the mesh nodes in the model.In this particular application, the weights and thresholds are obtained by a conjugate-gradient method.We seek for minimization of the residual error on a sub-set of 5 min data sampled from the fine resolution LES run as the input vector ) ( k x X and the coarse resolution LES runs as the output vector . Finally, N M , are numbers of neurons on the hidden layers of the ANN and k N is the number of neurons linked to each input vector element.Actual calculations have been run in MATLAB using neural network toolbox (nnet) by Demuth and Beale.In this sense, the application of the ANN in this study is similar to the approach chosen by Manjunatha et al. (2010).The ANN application requires preconditioning through the model feedback, as large predictor vector dimension leads to a classical ill-conditioned problem of the linear algebra (Hsieh and Tang, 1998).
A correct way to introduce the unresolved boundary conditions into the LES model is through the sub-grid scale model.The sub-grid scale model describes the integral effect of unresolved flow features on the resolved scales.Sarghini et al. (2003) pioneered application of ANN to dynamic sub-grid scale modelling.LES equations are obtained through filtering of Navier-Stokes equations.The filtering result in appearance of a sub-grid stress term which must be modelled.Here the overbar with the superscript h denotes filtering operator at scale twice the grid scale while the grid filtering is done implicitly by the numerical scheme; is the mixing length scale; s C is the Smagorinsky constant; Δ is the LES mesh resolution.The strain rate tensor is ( ) where i ∇ denotes a gradient in i-th direction.As Eq. (2) reveals, s C is the parameter to be tuned to the specific unresolved features of the flow.Here we consider the following numerical experiments to test the ANN utility in the described problem of the flow around large unresolved obstacles below the first computational layer.
The LESNIC model run with developed atmospheric turbulent boundary layer with the geostrophic wind speed 4 m s -1 at latitude 45 degrees over homogeneous flat surface with the roughness length scale 0.01 m.The model domain was 1000 m by 750 m by 500 m in streamwise, spanwise and vertical directions with the mesh 64 by 48 by 64 grid nodes.The boundary layer evolved over one hour of model time.At the end of its evolution, the boundary layer occupied ½ of the domain depth.The last 3-dimensional state of the model was stored and later used as initial conditions to obtain additional 30 min runs both with homogeneous and heterogeneous surfaces.The simulations were run for truly neutral boundary layers with passive scalars without feedback effects on the flow.The lateral boundary conditions were periodic in all runs.At the surface, the turbulent flux of momentum was computed instantly and pointwise using the log-law and therefore did not account directly to the surface features.The ANN was utilized in the following way.The 3 layer back-propagation feeding ANN was trained on 5 min data from the 30 min additional LES runs described above.The input vector for the ANN was composed of components of the tensor in Eq. ( 3) and its module as

=
. (4) Since the tensor is symmetric, its off-diagonal components receive larger weights in the input vector.
Let us first consider the effect of the ANN in the Ekman layer simulations over flat, homogeneous surface.
Figure 2 compares histograms of s C computed with Eq. ( 2) only and with Eq. ( 2) following the ANN correction.Equation ( 2) represents a dynamic mixed model as described in Esau (2004).This model uses significant resolved scale information, obtained through an explicit discrete Gaussian filtering, to estimate the magnitude of the sub-filter stress.This so-called scale-similarity information becomes insufficient when a large fraction of the turbulent kinetic energy is unresolved.The mismatch between the required within the inertial interval of scales and the resolved spectral energy fluxes is compensated by lowering s C in the dynamic-mixed model in Eq. ( 2) (for details see Pope, 2004).However, too low s C results in numerical instability.Numerous ad hoc approaches have been proposed to address the low s C values.Desirable stabilizing modifications of the s C distribution introduced by the chosen ANN are clearly seen in Figure 2. The small values of s C , which destabilize the simulations, were almost eliminated and too large values ( s C > 0.4) occurred considerably less frequent in the flow during the analysed 30 min interval of time.Thus, we can conclude that sub-grid modelling with ANN correction may improve performance and stability of LES even in simulations over homogeneous surfaces.
The next step is to introduce a surface heterogeneity in form of rough spots that represent completely unresolved obstacles or just surface heterogeneity of any other kind.This experiment is LES H1 in Table 1.The surface properties remind unchanged during the simulations.Hence ANN can be trained to correct the simulations specific features of the perturbed flow.This training is achieved in the study trough the run with twice as fine resolution, which served as observations for the training process.The run LES-H1 had only surface roughness modified.The roughness was set to 5 m in a square of ¼ of the domain surface area; whereas the rest of the domain had roughness of 0.01 m. Figure 3 shows the mean values of s C in the coarse resolution run without ANN (Figure 3a) and with ANN (Figure 3b) correction.The added value of the ANN correction is clearly seen in localized modification of s C and therefore the surface sub-grid fluxes as related in Eq. (2). Figure 4 shows that ANN was able to improve the vertical profile of averaged s C as well.ANN has automatically introduced modifications, which were previously left for ad hoc fitting as described in Balaras et al. (1995).The observed changes in the mean Smagorinsky coefficient profiles as well as in their spatial distribution and clipping of the extreme values are welcome modifications made through ANN implementation to LES closure problem over roughness heterogeneous surface.
At the next step, a cubic obstacle was defined in the boundary conditions of the coarse resolution run (LES-H2).The obstacle was under-resolved in the sense that only the 1 st model level was modified to have immersed boundary conditions.The log-law was applied only at the upward facing facet of the obstacle, whereas the simple non-slip conditions were applied on the lateral facets.The obstacle in LES-H2 is of 12 m tall (2 computational levels) located in ¼ of the domain area.Thus, the LES-H1 run represented a resolved flow response on unresolved elements of urban morphology, and LES-H2 represents a corresponding response on semi-resolved elements.Figure 5 shows values of s C obtained as in Figure 3 but for LES-H2 experiment.An impenetrable obstacle indeed has strong impact on s C values and their spatial distribution.As in Glazunov (2009), s C is distributed asymmetrically around the obstacle with larger values in the front (upstream) facet and smaller values at the rear (downstream) facet.Remarkably, ANN correction was able to account for the localization of the smallest s C on the downstream sides of the obstacle and very close to the centre of the obstacle front facet.Those features were obtained in a fine resolution high quality LES by Glazunov (2009) but hardly reproduced by the dynamic closure in this LES.It is also worth noticing that the time needed to estimate s C with the ANN correction is about 2 orders of magnitude less than time used in the Glazunov's LES (personal communication) and an order of magnitude less than time used to compute it in the present high-resolution LES.

Artificial neural network in modeling of concentrations
One of the most important and interesting applications of ANN in conjunction with LES is the modelling of the scalar concentrations (pollution).The problem is to predict the concentration field ) , ( t x q r on the basis of the velocity field ) , ( t x u r r , which is known with a limited accuracy, and ) , ( , which is known only in a few positions k x .Using the control LES run, the application was computed for 1 min forecast of concentrations in the near field of a pollution source in the turbulent boundary layer.The forecast was computed 25 times over 25 min of model time.Initially, the ANN has been trained on 5 min data in order to produce ) , ( t x q r .To do the training, we use the same LES run assuming that the LES data supplant the real atmospheric observations.To avoid self-correlation, the data computed for 30 min before the forecast period have been used.In practice, the observations and the LES together should be used for the ANN training.The technology is known (e.g.Mamarella et al., 2006) but since it is not used here, we will not describe it.The input vector was composed of components of velocity and the historical value with t δ = 1 min.This choice is to very much degree arbitrary.
It could be motivated only by an expert opinion that it is difficult to get the observational data of higher frequency.Hence, the input vector was defined as for each mesh node k x .Figure 6 shows concentration field averaged over 30 min and within the layer 15 ± m above and below the injection source.The source was located at the height of 18 m above the surface so that the polluted wake meanders in three dimensions.The spots of high concentration are created by the wake meandering in the vertical direction.A nontrivial feature is a rarefied scalar concentration directly downstream from the source on the wake axes recognized also in Babatunde and Enger (2002).The ANN predicts robustly short-term fluctuations of the concentration but produce some noise in the absence of a real signal.The overall weighted correlation between the simulated concentration and the ANN averaged forecast was 0.92 in this experiment.The forecast interval of 1 min would probably seem too short for a practical use but the reader should keep in mind that we are analysing the ANN as an addition to the LES hydro-dynamical model.Thus, the ANN role is utilized only for predictions on unresolved intervals of time.The advection across larger distances and longer time scales are explicitly resolved in the model.Indeed over 1 min event with the wind speed of 4 m s -1 a polluted air volume would be transported over 240 m or more than 10 grid intervals in the present experiment.With typical wind of 10 m s -1 , the same volume will be transported over 600 m or 50 grid intervals.

Discussion: Compatibility of LES and ANN approaches
Artificial Neural Networks (ANN) were seen in several publications (e.g.Azamathulla et al., 2008) as a promising method to reduce dimensionality of models required to describe complex properties of the flow dynamics, dispersion in the flow and the flow interaction with the surface.These promises had not been met.One, probably insurmountable problem has been pointed out by Milano and Koumoutsakos (2002).To present it, these authors noticed analogy between the ANN and the Proper Orthogonal Decomposition (POD) techniques.
The linear ANN could be used for the POD analysis.Alternatively, the POD analysis could be considered as another training method to reduce dimensionality of the problem.In the case of POD, the training procedure is simply singular value decomposition retaining N first eigenvectors.As with POD, the dimension of the problem for any realistic flows is too large and untreatable.The flow domain is to be decomposed on small independent sections to circumvent the problem.If those sections are essentially interdependent, the decomposition becomes loose in the sense that essential part of the information will be lost.
In the ANN, the computational expenses are concentrated on calculation of the weight matrix W .Let us show that the POD and the linear ANN are equivalent.The POD could be written as where ) (t a is the eigenvalue vector of the size iables z y x n n n n var , Φ is the eigenvector matrix, and ) (t u is the transposed variables' vector.The ANN could be formulated as where 1 W and 2 W are the ANN weight matrices of the size iables z y x n n n n var to be found from a mini-max problem.Now obvious that while the POD problem can be hugely accelerated by choice of the orthogonal harmonic functions, the Fast Fourier Transformation (FFT) and the fast partial decomposition algorithms, the ANN problem should involve essentially slow iterative optimization problem solvers.The problem complexity increase with non-linear networks that can be expressed as where already 4 weight matrix are to be obtained in the training process.Milano and Koumoutsakos (2002) reported the largest trained network of just z y x n n n = 42 x 39 x 16 = 26 208 elements, which result in the size of the weight matrix 1 W of 6 710 272 elements.Large dimensionality of the optimization problem makes it not only computationally difficult but also ill-conditioned.Thus, any flow non-locality virtually prohibits the application of the ANN.Even neglecting non-locality, the problem remains untreatable.
There is however a serious deficiency in the Milano and Koumoutsakos analysis.They considered a global optimization problem in the sense that they reduce dimensionality within each of sub-domains solving the global mini-max problem.In the application approach proposed and tried in this study, the low-dimensional problem has been formulated for each mesh cell.In fact, the resolved-scale dimension of the problem here is remains O(10) or typically about O(100 m) in the length scale terms and therefore treatable.We simply use the scaling of the problem imposed by the specifics of the urban-type environment where the buildings are different but not particularly variable in size.Thus, the Milano and Koumoutsakos approach is insightful.But the ANN application for the considered urban problems is not obviously disadvantages.
Another approach to the ANN applications has been proposed by Jimenez et al. (1997).It refers to rigorous LES formalism, exemplified in e.g.Cook and Riley (1994) model.The formalism demands to express the sub-grid quantity ' q through resolved quantity q in a universal scale-similar manner.Cook and Riley exploited the fact that any quantity defined at [0;1], i.e. normalized q , and subjected to stochastic turbulent fluctuations u with the Gaussian distribution function should have the β-distribution where ) (a Γ is the gamma function of Euler With knowledge of the β-distribution, the true values of concentration and its ensemble-averaged value over a grid cell are related as In principle, Equations ( 8)-( 9) allow estimations of the concentration within the model grid cell with a minimal root mean square error.As Jimenez et al. (1997) argued, the envisaged model provides excellent results for a homogeneous flow even at high Re.The situation is however more complicated in heterogeneous urban sub-layer.Here, space-time filtering applied to compute 2 q σ and ultimately > < q is not valid any more due to presence of the vertical walls.The resolved variations of the concentration are far too small to represent scale similarity.

Conclusions
The micro-meteorological modelling tool, large-eddy simulation model, is well established for the research purposes.In order to utilize this tool in practical applications, unresolved surface and emission features with effect on the resolved turbulence have to be represented.The study addressed the problem investigating utilization of the artificial neural networks with the turbulence resolving large-eddy simulation (LES) model LESNIC.A priori analysis has been applied for idealized, relatively coarse resolution numerical experiments, namely, the experiments with turbulent neutrally stratified boundary layers over homogeneous surface, the surface with very large roughness patch and the surface with an obstacle.It is shown that the filtering property of the three layers' ANN is useful in improvement of numerical stability and physical consistency of the dynamic sub-grid closure in LES.The 3-layer tan-sigmoid ANN with back-propagation training over rather short interval of 5 min outperformed the usual engineering clipping approach to obtain the dynamic values of the Smagorinsky constant.Contrary to the scepticism induced by the earlier applications of the ANN to the sub-domain reduction in dimensionality, the ANN application for each of dynamical sub-regions (usually just of a few computational cells) in the LES domain was found promising and computationally efficient.The key to the problem is the fact that the small scale turbulent perturbations quickly decay.Thus the radius of locality of the problem is small and limits the size of the weight matrices in the ANN.The ANN is shown to be a robust predictor for the pollutant scalar concentration in the urban sub-layer with unresolved scalar sources.
Table 1.Domain-averaged correlation between actual, i.e. simulated in dynamic LES, sub-grid turbulent stress in the high resolution run and the stress in the control run with the ANN correction using the input vector X over post training period of 25 min.

Experiment
to filtering of too low s C values in ANN; fraction of s C < 0.15 is 48%.Heterogeneous surface roughness, LES-H1 0.595 Domain-averaged correlation improves as there is generally higher s C .Fraction of s C < 0.15 is 40%.Heterogeneous surface with an obstacle, LES-H2 0.581 Domain-averaged correlation improves as there is generally higher s C .Fraction of s C < 0.15 is 39%.

Figure 1 .
Figure 1.The surface elevation obtained for Paris area with the resolution of 10 m (1 m around Place Concordia).The gray scale indicates the elevations in meters.Courtesy MEGAPOLI team.

Figure 2 .Figure 4 .
Figure 2. Distribution of s C at the levels 1 through 6 in the LES-H1 without (squares) and with ANN (dots) correction.In ideal case, the distributions should be similar; the additional information resulting in the shown modification appears from the scale-similarity assumption, i.e. the averaged behavior of the fine-resolution run used for training remains valid for the coarse-resolution run.