Towards a new multiscale air quality transport model using the fully unstructured anisotropic adaptive mesh technology of Fluidity (version 4.1.9)

An integrated method of advanced anisotropic hr-adaptive mesh and discretization numerical techniques has been, for first time, applied to modelling of multiscale advection–diffusion problems, which is based on a discontinuous Galerkin/control volume discretization on unstructured meshes. Over existing air quality models typically based on static-structured grids using a locally nesting technique, the advantage of the anisotropic hr-adaptive model has the ability to adapt the mesh according to the evolving pollutant distribution and flow features. That is, the mesh resolution can be adjusted dynamically to simulate the pollutant transport process accurately and effectively. To illustrate the capability of the anisotropic adaptive unstructured mesh model, three benchmark numerical experiments have been set up for two-dimensional (2-D) advection phenomena. Comparisons have been made between the results obtained using uniform resolution meshes and anisotropic adaptive resolution meshes. Performance achieved in 3-D simulation of power plant plumes indicates that this new adaptive multiscale model has the potential to provide accurate air quality modelling solutions effectively.


Introduction
It is well-known that the interaction of multiscale physical processes in atmospheric phenomena poses a formidable challenge for numerical modelling (Kühnlein, 2011).Largescale processes can trigger small-scale features that again have an important influence/feed back to the large scale (Behrens, 2007).For example, the processes of tropical cyclone involve a range over a continuous spectrum of scales from the large-scale flow environment ∼ O(10 6 -10 7 ) m, tropical cyclone itself ∼ O(10 5 -10 6 ) m, embedded eye wall and rainbands ∼ O(10 3 -10 4 ) m, as well as down to microscales of the boundary layer turbulence ∼ O(10-10 2 ) m (Kühnlein, 2011).For air pollution, the dynamic and chemical processes also involve a wide range of scales.The initial transformation of emissions from urban and industrial centres or dispersion of plumes from large power plant stacks occur on relatively small scales, but would be engaged to much larger scales after long range transport.It is a gargantuan computational challenge to modelling large regions with uniform resolution at the finest relevant scale.Therefore, mesh adaptation may be the only effective way to encompass different scales (e.g. local, urban, regional, global) Published by Copernicus Publications on behalf of the European Geosciences Union.

J. Zheng et al.:
A new multiscale air pollutant transport model using anisotropic adaptive mesh methods in a unified modelling system to better capture the interactions among the processes relevant at each scale (Garcia-Menendez and Odman, 2011;Kühnlein, 2011;Weller et al., 2010;Nikiforakis, 2009).
So far, the accurate numerical modelling of advection (or transport) remains a central problem for many applications such as air pollution, atmospheric chemistry, meteorology and other physical sciences.There have been many studies on the numerical advection schemes (e.g.PPM, Bott and Walcek etc.) that have been used in many air quality models (e.g.CMAQ, CMAx, NAQPMS etc.) (Colella and Woodward, 1984;Bott, 1989;Walcek and Aleksic, 1998).These advection algorithms were implemented based on a fixed uniform mesh system.The successive global refinement can be used to capture the details of small-scale flow features, but is prohibitively expensive and not feasible for practical applications.Alternatively, the nesting technique, placing finer meshes within coarser meshes, is often used for achieving local higher resolution in many air quality models (Garcia-Menendez and Odman, 2011;Frohn et al., 2002;Wang, 2001).In static mesh nesting, the solutions obtained from the global coarse mesh model provide the boundary conditions for the nested mesh regional model; in turn, the solutions in the global model are updated with the high resolution solutions.However this may lead to spurious oscillations at the interface between the coarse mesh and nested fine mesh, especially when concentration gradients are large cross the interface.Although the numeral techniques such as blending, nudging and selective damping approaches can be used to remove these oscillations, the small-scale features on the fine meshes may be damped (Garcia-Menendez and Odman, 2011;Zhang et al., 1986;Debreu and Blayo, 2008;Alapaty et al., 1998).Moreover, due to highly unsteady atmospheric flows, it is almost impossible to construct a static optimal nested mesh suitable for an accuracy simulation over a long time period.The use of dynamically adaptive mesh techniques can therefore be considered in such as way that the mesh resolution can be adjusted locally in response to the evolution of the flow and passive tracer (Piggott et al., 2009;Behrens, 2007).
In contrast to locally nested mesh techniques, adaptive mesh techniques can not only resolve multiscale processes in a consistent way, but also enable to follow and capture the features of flows as time evolves.Dynamic mesh adaptation can be achieved, either by re-locating mesh nodes or by locally increasing (and decreasing) the number of nodes in time and space.The former, known as mesh movement (i.e.r-adaptive technique), can be used to improve the accuracy of solutions by optimally re-locating mesh nodes to resolve the small-scale features of interest (Garcia-Menendez and Odman, 2011;Srivastava et al., 2000;Lagzi et al., 2009;Kühnlein et al., 2012;Nikiforakis, 2009).However, the accuracy of solutions using r-adaptivity is restricted a priori for achieving an optimal dynamic mesh (where the total number of nodes is fixed).The latter, known as mesh en-richment (i.e.h-adaptive technique), can guarantee a minimum solution accuracy level by providing sufficient resolution where and when it is needed (Baker et al., 2013;Constantinescu et al., 2008;Piggott et al., 2005).Various hadaptive techniques based on structured meshes as well as the r-adaptive techniques on unstructured/structured meshes have been explored in atmospheric modelling; furthermore some of these techniques have been applied to air quality models (Garcia-Menendez and Odman, 2011).Recently, significant research efforts have been focused on application of this new adaptive mesh techniques in ocean modelling (Pain et al., 2005;Piggott et al., 2009Piggott et al., , 2008a, b), b).
This article applies a new anisotropic hr-adaptive mesh technique into air quality transport (advection) modelling.This adaptive unstructured mesh technique provides the dynamic spatial and temporal resolution to capture moving features, e.g.moving fronts or power plant plumes.Using the hr-adaptive technique, existing elements can be split (hadaptive) or element vertices can be moved (r-adaptive), to periodically modify the mesh geometry.Hence, the purpose of this article is to demonstrate, through example problems, the capability of anisotropic mesh adaptivity for modelling of multiscale transport phenomena.
The remaining structure of this article is as follows: Sect. 2 describes numerical methods, including discontinuous Galerkin (DG) and control volume (CV) methods based on unstructured meshes.Section 3 covers the topics of mesh adaptivity, error measures and interpolation.Section 4 introduces three-dimensional (3-D) unstructured anisotropic adaptive mesh model (Fluidity).Section 5 discusses its performance in three benchmark advection problems and a model problem for dispersion of power plant plumes.Conclusions are drawn in Sect.6.

Numerical methods for transport equation
As a model problem, we consider the generic transport equation for a scalar quantity, c, is given in conservative form by where u = (u, v, w) T is the velocity vector, κ is the diffusivity (tensor) and s represents any source or reaction terms.If κ = 0 and s = 0, Eq. (1) reduces to the advection equation

Spatial discretization
Integrating Eq. ( 2) by part over the computational domain , its weak form can be written (3)

Discontinuous Galerkin discretization
As a locally conservative, stable and high-order accurate method, the discontinuous Galerkin methods can easily construct discontinuous approximations on unstructured meshes to capture highly complex solutions and are well-suited for hr-adaptivity and parallelization (Cockburn et al., 2000;Cockburn and Shu, 2001;Flaherty et al., 2002;Hesthaven and Warburton, 2007).Moreover DG methods, as a generalization of finite volume methods, can directly make numerical fluxes and slope limiters available in the finite element framework (Burbeau et al., 2001;Hoteit et al., 2004;Krivodonova, 2007;Krivodonova et al., 2004).Integrating Eq. ( 2) over a single element and summing over all elements, we obtain where the hatted term represents fluxes across the element facets.If κ = 0 and s = 0, Eq. (4) becomes a pure advection equation Due to the discontinuous nature of fields, there is no unique value for the flux term; however the requirement that c is a conserved quantity, does demand that adjacent elements make a consistent choice for the flux between them.In this work, two advective flux schemes, the upwind and local Lax-Friedrichs flux methods, are used to represent n • u c for DG methods (AMCG, 2014).In n • u c, the advecting velocity u can be calculated by either averaging it on each side of the face or applying a Galerkin projection to project the velocity onto a continuous basis.
In the upwind flux formulation, the value of c at each quadrature point on the face is taken to be the upwind value; that is, if fluid is into/out of the element then it is the value on the exterior/interior side of the face.Integrating the advection term by parts twice, then Eq. ( 5) becomes (AMCG, 2014): where û represents the flux velocity and a weakly imposed boundary condition c = c b is applied on the inflow part of boundaries; c ext and c int are the values on the exterior and interior side of the face respectively.
In local Lax-Friedrichs flux formulation, the tracer advection is given by where for each facet s ⊂ ∂e: Here, "sup" is the abbreviation of supremum.
To ensure nonlinear stability and effectively suppress spurious oscillations, the slope limiting techniques are used here (Kuzmin, 2010;Cockburn and Shu, 2001;Luo et al., 2007).

Control volume discretization
The control volume discretization uses a dual mesh constructed around the nodes of the parent finite element mesh.Once the dual control volume mesh has been defined, it is possible to discretize the transport Eq. (1) using piecewise constant shape functions within each volume, v. Integrating Eq. ( 1) by parts within a volume, v and summing over all volumes, we obtain For the flux term n • u c, the velocity is well-defined since the control volume facets are in the centre of the elements of the parent mesh where it is continuous.The face value of c k is computed at each quadrature point of the facet k using the finite element interpolation approach, i.e. interpolating it using the finite element basis functions on the parent mesh.Usually the first-order quadrature is performed on the control volume facets; however if higher-order control volume facet quadrature is selected then k refers to each quadrature point on the facet.To avoid spurious oscillations, the CV-TVD (control volume -total variation diminishing) limiter is used to make the solutions total variation diminishing (Sweby, 1984;AMCG, 2014).
For diffusion term n • κ∇c, ∇c is treated with controlvolume element-based gradients, equal-order Bassi-Rebay and staggered mesh Bassi-Rebay discretization (for details, see AMCG, 2014).the maximum mesh size is set to be 0.2).

Time discretization
The semi-discrete matrix form of Eq. ( 3) can be written as in which the vector c = (c 1 , . .., c N ) T contains the solution of variable c at nodes (N is the number of nodes), M is the mass matrix, A(u) is the advection operator, K is the diffusion operator, and r is the right-hand side vector containing boundary, source and absorption terms, where for continuous Galerkin discretization: The time derivative term at time level n + 1 is treated using the θ method to yield where θ ∈ [0, 1] and the terms c n+θ are given by Equation ( 12) can be rearranged for unknown vector c n+1 : Equation ( 14) can be solved in two stages: where r in Eq. ( 14) is split into Dirichlet r D and Neumann boundary components r N , and a source component r s .
For discontinuous Galerkin discretization, the explicit Euler scheme (θ = 0) is used in Eq. ( 15).An advection subcycling method based upon a CFL criterion or a fixed number of subcycles is adopted in modelling advection flows; that is, the time step t is split to N sub-time-step t sub = t N to satisfy the specified Courant number: To guarantee a bounded solution, the slope limiter is applied to c new after each sub-time-step.Note that the matrix M − ( t/N )A (u n ) is constant within one time step.Therefore the process of solving Eq. ( 17) only involves the matrixvector multiplication, thus reducing a large amount of the  CPU time required for assembling the matrices, especially when unstructured meshes are used.For control volume discretization, an explicit scheme is simple but strictly limited by the CFL number, which can be restrictive on adaptive meshes as the minimum mesh size can be very small.Here, we adopt a new time stepping θ scheme based on traditional Crank-Nicolson scheme (θ aim = 1/2) because of its robustness, unconditional stability and secondorder accurate in time (Pavlidis et al., 2015;Versteeg and Malalasekera, 2007;Donea and Huerta, 2003).For the given time step, the value of θ min can be estimated at each CV face based on the satisfaction of a total variation diminish-ing (TVD) criterion.Therefore, for each control volume v, we can choose θ v ∈ [θ min , 1] to be as close to θ aim as possible; that is, θ v = max{θ min , θ aim }.In this way, it can eliminate the local time step restriction for physically realistic and bounded solution although it may be at a cost of losing some local accuracy (for details, see Pavlidis et al., 2015).

Mesh adaptivity
The optimization-based adaptivity technique, developed by the Applied Modelling and Computation Group (AMCG) at Imperial College London (AMCG, 2014), is introduced The key objective of using adaptive mesh methods is to reduce the overall computational cost in achieving an error goal; thus ensuring that fine resolution is used only when and where it is needed (Fang et al., 2010;Pain et al., 2001).A error metric tensor to guide an adaptive meshing algorithm can be defined (Pain et al., 2001): where is the required level of errors and γ a scalar constant (here, γ = 1) and H is the Hessian matrix of variable fields (here, the tracer concentration c): The absolute value of the symmetric Hessian matrix is defined as (Pain et al., 2001) where the matrices V and contain the eigenvectors e i and eigenvalues λ i of the Hessian matrix H, respectively.The required edge length in the direction e i to achieve the required level of errors can be obtained (Piggott et al., 2009): The rotation matrix V in combination with can be used to adapt the original element to an anisotropic element required  for the given level of errors.To bound the aspect ratio of elements in physical space, the eigenvalues of the metric can be modified (Pain et al., 2001): where where a is the a given aspect ratio of elements, and h min and h max are the minimum and maximum sizes of elements, respectively.
To represent small-scale dynamics, a relative error metric formulation is suggested: where c is the field under consideration, is now a relative tolerance and min is the minimum tolerance used to ensure that the denominator never becomes zero.
To guide refinement/coarsening of the mesh, the maximum and minimum mesh sizes are set to allow one to  impose different limits in different directions (for details, see AMCG, 2014).Assuming that these directions are aligned with the coordinate axes allows one to define diagonal tensors.The maximum and minimum number of nodes are also set for mesh adaptivity.This is effected by computing the expected number of nodes from the given metric.If the expected number of nodes is greater than the maximum number of nodes, the metric resolution is homogeneously decreased so that the expected number of nodes is the maximum number of nodes.
Another key issue of mesh adaptivity is to interpolate any necessary data from the previous mesh to the adapted one.The consistent interpolation is often adopted in mesh adaptivity.However, the consistent interpolation can introduce a suboptimal interpolation error, unsuitability for discontinuous fields, and lack of conservation.An alternative conservative interpolation approach, the Galerkin projection is proposed for discontinuous fields.A supermeshing algorithm (Farrell et al., 2009) is used for implementation of the Galerkin projection.are available for mesh-to-mesh interpolations between adaptations.
Fluidity is parallelized using MPI and is capable of scaling to many thousands of processors.It has a userfriendly GUI and a python interface that can be used to calculate diagnostic fields, set prescribed fields or set user-defined boundary conditions (for details see https://www.imperial.ac.uk/engineering/departments/ earth-science/research/research-groups/amcg/).It is noted that version 4.1.9 of Fluidity is not necessarily required, but older versions might work as well.

Numerical examples
To illustrate the efficiency and accuracy of anisotropic adaptive schemes, four benchmark problems have been adopted,   which are representative and challenging enough to predict how the new adaptive multiscale model would behave in future real-life applications (LeVeque, 1996;Kuzmin, 2009;Staniforth et al., 1987;Walcek and Aleksic, 1998;Bott, 1989Bott, , 1993Bott, , 2010)).
In the following comparative study, we consider FEM_Fix and FEM_Adapt schemes (FEM represents CV or DG) based on the control volume and discontinuous Galerkin discretization.The CV_Fix_L and DG_Fix_L schemes use fixed uniform triangular meshes, while the CV_Adapt_L and DG_Adapt_L schemes use adaptive meshes (where L represents the different mesh schemes, as shown in Table 1).For CV discretization, a finite element interpolation is used at the control volume faces with a CV-TVD limiter to bound the solution.The time discretization used here is the new time stepping θ scheme based on the Crank-Nicolson scheme.For DG discretization, the upwind flux is chosen in combination with vertex-based slope limiter.The slope limiter used with the discontinuous Galerkin formulation only guarantee a bounded solution in conjunction with an explicit advection scheme.Therefore, advection subcycling based upon a CFL criterion is necessary for DG discretization (AMCG, 2014).Equation ( 14) is solved by the generalized minimum residual method (Saad, 1993).The successive over-relaxation precondition is invoked to speed up convergence at large time steps.
To assess the difference between the analytical solution c and its numerical approximation c h , we introduce the error norms: The order of accuracy in modelling is used to assess the numerical convergence rate: where h is the mesh size.All computations were performed on a workstation using the Gfortran Compiler for Linux.The simulation workstation has 8 processors and a 4 GB random access memory (RAM).The processor used in workstation is Intel(R) Core(TM) i7-2600 CPU @ 3.40 GHz.A single processor with frequency of 3.40 GHz was used since the test cases were simulated in serial.

Case one: solid body revolution
A standard test problem applied to the advection Eq. (2) in 2-D is solid body revolution (LeVeque, 1996;Kuzmin, 2009).The incompressible velocity field in the domain = which corresponds to a anticlockwise rotation around the centre (0.5, 0.5) of .Following LeVeque (1996), we consider a slotted cylinder, a sharp cone and a smooth hump as the initial solid bodies defined within the circle centred at each reference point (x 0 , y 0 ): where r 0 = 0.15.After each full revolution (t = 2π k), the exact solution is returned to the initial distribution as depicted in Fig. 1.For the slotted cylinder, the reference point is (x 0 , y 0 ) = (0.5, 0.75) and The cone is centred at (x 0 , y 0 ) = (0.5, 0.25) and its geometry is given by The peak of the smooth hump is located at (x 0 , y 0 ) = (0.25, 0.5) and the shape function is In the rest of the domain , the solution of Eq. ( 2) is initialized by zero.The challenge of this numerical test case is to preserve the shape of the rotating bodies as time evolves.
The mesh size used for the FEM_Fix_L schemes and the FEM_Adapt_L schemes are listed in Table 1.The time step is set to t = 0.01 for all different mesh schemes.For DG discretization, the explicit advection subcycling scheme with a tight CFL criterion (here 0.1) is used to make sure that the simulation is converging as the mesh is refined.For CV discretization, although the time stepping θ scheme based on the Crank-Nicolson scheme can maintain high accuracy, the subcycling number is set to be {2, 4} for h = {1/400, 1/800} respectively such that the sub-time-step is small enough to guarantee convergence and higher accuracy.
Figure 2 shows the errors of results at t = 2π (one full revolution) and the CPU time required.It can be seen that compared with the CV method, the DG method is more accurate but requires more computer memory and CPU time.For the CV method, the accuracy of results using the adaptive mesh scheme is very close to that using the fixed mesh (global mesh refinement) scheme, while the CPU time required by the adaptive mesh scheme is reduced by 75 %.For the DG method, to achieve a given level of accuracy of results, for example, E 1 = 0.0055 and E 2 = 0.035, by using adaptive meshes, the CPU time can be reduced by 45 % of that required using fixed meshes.With increasing mesh resolution, the CPU time for the adaptive schemes increases at a much slower rate than those for the fixed (global mesh refinement) approach (see Fig. 3 and Table 1).Compared with that in the fixed mesh (global mesh refinement) schemes, the problem size is reduced by 68-97.7 % using the adaptive mesh schemes.Hence, the use of adaptive meshes provides an efficient approach to lower the storage requirement,  thus leading to the reduction of the overall computing time, while remaining the accuracy of numerical results.To estimate the rate of convergence, the order of accuracy is calculated in Eq. ( 26) (here, h = 1/200).The order of accuracy is {0.83, 0.54, 0.95, 0.72} for {CV_Fix, CV_Adapt, DG_Fix, DG_Adapt} schemes, respectively.It is argued that non-smooth profiles in the complex problems presented here lead to a low order of accuracy, that is, a low convergence rate.If we only consider the hump-smooth profile as the initial data, the order of accuracy can increase to be {1.98,1.52, 1.54, 1.13}.Figure 4 shows the numerical results at t = 2π (after one full revolution) using the adaptive and fixed mesh schemes.For comparison purpose, the FEM_Fix_1 and FEM_Adapt_4 schemes are chosen since the number of nodes in these two mesh schemes are almost the same, where N = 10 201 for FEM_Fix_1 scheme while N ≈ 11 500 for FEM_Adapt_4 scheme.The solutions of CV_Fix_1 and DG_Fix_1 are computed on a structured uniform mesh of triangular elements with mesh size h = 1/100 and t = 0.01.It can be seen that there is severe erosion of the slotted cylinder when the fixed mesh scheme is adopted.The adaptive mesh scheme provides an improvement in accuracy of results.It is shown that with use of adaptive meshes (especially DG Adapt_Adapt_4), the initial shape of bodies is well preserved.

Case two: swirling flow
The capability of the adaptive mesh model has been further demonstrated in modelling swirling flow phenomena.The set-up of the simulation in this case is similar with case one; however the velocity field is provided by the formula (LeVeque, 1996;Kuzmin, 2009): where g(t) = cos(π t/T ) on the time interval 0 ≤ t ≤ T (here T = 1).
The initial mass distribution will be deformed by the timedependent velocity field, which gradually slows down to zero and reverses its direction at t = T /2.Thus, the initial profile will be reproduced at the final time t = T as depicted in Fig. 1.Due to the flow here being time variable, the time step is set small enough to be t = 0.00125 for all different mesh schemes.A comparison of results using fixed and adaptive meshes is illustrated in Figs.5-8.Again, it can be observed that by using the adaptive mesh scheme in the model, both the CPU time and number of nodes required are significantly reduced for a given level of accuracy of results (see Fig. 6).To improve the stability of solutions when the mesh resolution is increased, the explicit advection subcycling based upon a CFL criterion is used for DG discretization, and the Crank-Nicolson scheme is used for CV discretization.In this case, the order of accuracy is {0.81, 0.68, 0.92, 0.79} for {CV_Fix, CV_Adapt, DG_Fix, DG_Adapt} schemes respectively.Again the convergence rate is low due to non-smooth profiles in solutions.
The numerical solutions in Figs.7 and 8 (at time levels t = T /2 and T ) were computed by different fixed and adaptive mesh schemes.Again adaptive mesh modelling is able to present better deformation of shapes at t = T /2 (Fig. 7) and preserve the initial shape after one full revolution (t = T ) much better than fixed mesh modelling Fig. 8.
Figure 9 displays the change of adaptive meshes as time evolves.It is observed the dynamic mesh adaptation algorithm is capable of following the evolution details of transient flows.As the simulation progresses, the mesh has to be adapted not only to the current solution profile but also to its expected shape in the future.It can be seen that the mesh is adapted to capture the details of local flows, i.e. increasing the resolution around the shape's boundary with anisotropic elements and then capturing the shape of deformed bodies.

Case three: swirling deformation
A comparison of the anisotropic adaptive mesh schemes with the Walcek (or Bott) scheme (Walcek and Aleksic, 1998;Bott, 2010) adopted by many air quality models has been undertaken in this section.The case used here was described in Staniforth et al. (1987).In this case, we only focus on the subdomain [0.24, 0.76]×[0.12,0.88] of .A cone is initially centred at (x 0 , y 0 ) = (0.5, 0.5) with a negative (−0.2) background as shown in Fig. 10 and its geometry is given by The velocity field defined by the following stream function (Staniforth et al., 1987): with u(x, y) = (−ψ y , ψ x ) = ( Ak sin(kx) sin(ky), Ak cos(kx) cos(ky)), where A = 0.08, k = 4π .Staniforth et al. (1987) defined two flow regimes (short time periods and long time periods) that have different evaluation criteria for the numerical advection schemes.Here, we focus on the evaluation of the first regime (short time periods) so that the numerical solutions should be compared with the analytical solutions in a qualitative manner.Figures 11 and 12 show the comparison of three different schemes' results with the analytical solution at time t = 3T /20 and t = T /5, where T = 2.6376.The solutions of the Walcek scheme were computed on a structured uniform mesh with h = 1/200 and t = 0.003297.For FEM_Adapt_128 schemes (see Table 1), they were computed on dynamic adaptive mesh with constant t = 0.006594.The minimum mesh size is 7.8125 × 10 −6 while the maximum mesh size is 0.2.
It can be observed the initial c field is split into two rotations within the areas of the two central vertices as time evolves.Since the spatial gradient of solutions increases as time evolves especially at the boundaries of the central vortices, high resolution of meshes around the boundaries is needed to present the sharp shape accurately.Due to lack of high resolution of meshes, the solutions using the Walcek scheme fail to represent the analytical one and maintain the shape distribution.For the Walcek scheme, at time t = 3T /20 (see Fig. 11), the gradients of numerical distribution begin to disappear at the upper and middle boundaries of the central vortices and nearly completely disappear at time t = T /5 (Fig. 12).By adapting the mesh in time and space, the mesh resolution increases around the boundaries of each vertices, thus improving the accuracy of results.There is close agreement between the adaptive mesh modelling results and the analytical ones although the gradients for CV_Adapt_128 scheme are not as strong as the exact solution.
The sequence of triangulations presented in Fig. 13 demonstrates that the dynamic mesh adaptation algorithm succeeds in locally refining the mesh in the vicinity of steep fronts; therefore reducing the amount of numerical diffusion and follow steep fronts as time evolves.To further reduce the number of elements, the anisotropic adaptive algorithm has been used for all the above adaptive mesh scheme, allowing the mesh to be adapted along different directions.As shown in Fig. 14, which depicts a close up view of locally adapted mesh, the adapted mesh size across the boundaries is small enough to capture the sharp fronts while large along the boundaries since the c field does not change much.Therefore, the mesh sizing desired in anisotropic adaptive algorithm is not only a function of space but also a function of direction.At a given point, the desired mesh sizing differs in different directions.
Figure 15 shows the number of nodes required for CV_Adapt_128 scheme is less than the node number (15 808) for fixed Walcek scheme during most of the simulation period.However, as local mesh resolution increases with time, the max CFL number of CV_Adapt_128 scheme exceeds 10 and even reach 80.To keep the stabilization of solutions, the time stepping θ scheme is used to eliminate the time-step restrictions and maintain high accuracy as far as possible, where θ v (1/2 ≤ θ v ≤ 1) is chosen at 0.5 for most of the elements while being big enough (close to 1) for a small fraction of individual elements with a large CFL number (see Fig. 16).In this way, the use of a large time step is acceptable when applying adaptive mesh techniques into comprehensive air quality models, which can make the computation much more efficient.As shown in Figs.11-12, in combination with the time stepping scheme, the adaptive mesh CV modelling solutions can stay stable and accurate without reducing the time step size, even if the max CFL number of CV_Adapt_128 exceeds 80.All of these can further illustrate the efficiency and the potential of dynamic mesh adaptation for future real applications in air quality model.

Case four: power plant plumes
In this case, the anisotropic adaptive mesh model is applied to an advection-diffusion problem (Eq.2): atmospheric dispersion of emissions from power plants.This is a first step towards applying the adaptive mesh model to realistic cases.The SO 2 emission of power plants was obtained from the Regional Emission inventory in ASia (REAS 2.1) data developed by National Institute of Environmental Sciences of Japan.As shown in Fig. 17a, the simulated domain covers the whole Shanxi-Hebei-Shandong-Henan region of China with 1090 km × 1060 km, and there are about 100 power plants in this area.The meteorological fields are provided by the mesoscale meteorological model WRF (v3.5) with a hor-izontal resolution of 5 km × 5 km and 20 vertical sigma layers.The simulation started at 00:00 UTC on 10 January 2013 and ran through to the 15 January 2013.In this case, the CV method is used for simulation of power plant plumes.
We started with a numerical investigation of a simplified 2-D test.The mixing layer height is 600 m and the turbulent horizontal diffusivity is 100 m 2 s −1 .The horizontal wind fields are obtained by averaging the lowest five layers of WRF's meteorological fields and stored at hourly intervals during 5-day period.For fixed mesh schemes, three mesh resolution levels in horizontal are used: 10 km×10 km (level 1), 5 km×5 km (level 2) and 2.5 km×2.5 km (level 3).For coarse meshes (level 1), there are 110 × 107 nodes and 23 108 elements.The total number of fixed elements increases by a factor of 4 when doubling the horizontal mesh resolution.For adapt mesh schemes, the minimum (maximum) mesh size is set to be 2 km (30 km), and the maximum number of nodes is set to be 12 000, which is the same as that of the fixed coarse mesh scheme (level 1).To represent the emission sources accurately, the fixed mesh with a high resolution of 2 km is used around the power plant points within a radius of 6 km (see Fig. 17b).
Figure 18 shows SO 2 concentrations at 21:00 UTC on 12 January after spin-up of simulations.An artificial dilution effect can be seen when coarse meshes are used in modelling.This can be improved by increasing the mesh resolution or applying an adaptive mesh scheme.The results using adaptive meshes are in agreement with those using fixed meshes with a high mesh resolution of 2.5 km while the number of nodes decreases by a factor of 16 with use of adaptive meshes.The evolution of adaptive meshes displayed in Fig. 19 illustrates that the adaptive algorithm is able to capture not only the detailed small-scale plume structures near the point sources, but also the regional high concentrations at large downwind distances.
To further demonstrate the adaptive mesh model's ability in 3-D modelling, we extended the above 2-D case to 3-D dispersion of plumes.According to the terrain data of the modelling domain, the initial 2-D adaptive mesh (see Fig. 17b) can be extruded to create a layered 3-D mesh from the top 20 km (above sea level) to the terrain surface, with 11 terrain-following layers.There are seven layers within the lowest 1 km above the terrain surface (see Fig. 20a).The power plant emissions were injected into the third layer about 200 m above the surface.Similarly, the 3-D velocity fields produced by WRF were interpolated from the fixed mesh in WRF onto the adaptive mesh.The vertical eddy diffusivity is parameterized based on a scheme by Byun and Dennis (1995).Figure 20 shows the evolution of 3-D SO 2 concentration visualization, which includes surface concentrations and the corresponding adaptive mesh, as well as the 3-D pollutant plumes defined as a constant concentration surface for concentrations greater than 100 µg m −3 .It can be seen that full 3-D mesh adaptivity has been used to improve the ability of the model to capture the details of flow dynamics and follow the evolution of power plant plumes.

Conclusions
In this paper, a new anisotropic adaptive mesh technique has been introduced and applied to modelling of multiscale transport phenomena, which is a central component in air quality modelling systems.The first two benchmark test cases using the fixed mesh and adapted mesh schemes have been set up to illustrate the efficiency and accuracy of anisotropic adaptive mesh technique, which is an important means to improve the competitiveness of unstructured mesh air quality models.The third case presents the irreplaceable advantage of this new adaptive mesh method to reveal detailed small-scale plume structure (large gradients) that cannot be resolved with static grids, using comparable computational resources.Dispersion of power plant plumes, as a real model problem, has been simulated in the last case to illustrate that the adaptive algorithm is able to capture the detailed small-scale plume structures near each point source as well as the regional high concentrations at large downwind distances.
It is demonstrated that the dynamic anisotropic adaptive mesh technique can be used to automatically adapt the mesh resolution to follow the evolving pollutant and transient flow features in time and space, thus reducing the CPU time and memory requirement significantly.In combination with the time stepping θ scheme based on the Crank-Nicolson method, the adaptive mesh air pollution model is able to maintain the stability and accuracy of results without reducing the time step size when the minimum mesh size is getting smaller.This is of great significance for the future applications in multiscale modelling.
The third test case serves as a proof-of-concept to further illustrate the capability of anisotropic mesh adaptivity techniques.In this case, the swirling deformation flow exhibits very high aspect ratios (1000, for example), which means that the pollutant distribution can possess very strong anisotropies as time evolves.Hence, the anisotropic mesh adaptation provides a very useful and effective way to simulate and represent this special atmospheric phenomena.
In summary, the results obtained in this work show the capability and potential of adaptive mesh methods to simulate multiscale air pollutant transport problems (spanning a range of scales) with higher numerical accuracy.The mesh adaptation can be used to improve the mesh resolution when and where it is needed without performing successive global refinement, which is prohibitively expensive, and therefore, not feasible for realistic applications.Future work will consider chemical reactions to further demonstrate the capability of dynamic adaptive mesh techniques.

Figure 2 .
Figure2.Case one -solid body revolution: the errors in the c field solutions and the CPU time (as a function of the mesh size h) required for one revolution (t = 2π), where h is the mesh size for FEM_Fix_L schemes while the minimum mesh size for FEM_Adapt_L schemes.

Figure 3 .
Figure 3. Case one -solid body revolution: the evolution of number of nodes for (a) CV_Adapt, (b) DG_Adapt.

Fig. 4 .Figure 4 .
Fig. 4. Case one -solid body revolution: the results from the fixed and adaptive mesh schemes using almost the same nodes number N , where N = 10201 for FEM Fix 1 scheme while N ≈ 11500 for FEM Adapt 4 scheme, at t = 2π.

Figure 5 .
Figure5.Case two -swirling flow: the errors in the c field solutions and the CPU time (as a function of the mesh size h) required for one revolution, where h is the mesh size for FEM_Fix_L schemes while the minimum mesh size for FEM_Adapt_L schemes, using the same t = 0.00125.

Figure 6 .
Figure 6.Case two -swirling flow: the evolution of number of nodes for (a) CV_Adapt, (b) DG_Adapt.

Fig. 7 .Figure 7 .
Fig. 7. Case two -swirling flow: the results from the fixed and adaptive mesh schemes using almost the same nodes number N , where N = 10201 for FEM Fix 1 scheme while N ≈ 12000 for FEM Adapt Figure 7. Case two -swirling flow: the results from the fixed and adaptive mesh schemes using almost the same node number N, where N = 10 201 for FEM_Fix_1 scheme while N ≈ 12 000 for FEM_Adapt_4 scheme, at t = T /2(= 0.5).

4
Introduction of a 3-D unstructured anisotropic adaptive mesh model (Fluidity, version 4.1.9)Thenew multiscale air quality transport model has been developed with a 3-D unstructured and adaptive mesh model (Fluidity; developed by the AMCG at Imperial College London).Fluidity, an open-source LGPL model, numerically solves the 2-D/3-D Navier-Stokes equation (being nonhydrostatic, to model dense water formation and flows over steep topography) and field equations with a range of control volume and finite element discretization methods.It includes a number of novel, advanced methods based upon adapting and moving anisotropic unstructured meshes, advanced finite element and control volume discretization, and a range of numerical stabilization and large-eddy simulation (LES) turbulence models.Among existing unstructured mesh models, Fluidity is the only model that can simultaneously resolve both small-and large-scale fluid flows while smoothly varying resolution and conforming to complex topography.The model employs 3-D anisotropic mesh adaptivity to resolve and reveal fine-scale features as they develop while reducing resolution elsewhere.A number of interpolation methods (e.g.non-conservative pointwise and conservative methods)

Fig. 8 .Figure 8 .
Fig. 8. Case two -swirling flow: the results from the fixed and adaptive mesh schemes using almost the same nodes number N , where N = 10201 for FEM Fix 1 scheme while N ≈ 12000 for FEM Adapt Figure 8. Case two -swirling flow: the results from the fixed and adaptive mesh schemes using almost the same node number N, where N = 10 201 for FEM_Fix_1 scheme while N ≈ 12 000 for FEM_Adapt_4 scheme, at t = T (= 1).

Fig. 9 .Figure 9 .
Fig. 9. Case two -swirling flow: the evolution of the adaptive mesh colored with tracer value c, where DG Adapt Figure 9. Case two -swirling flow: the evolution of the adaptive mesh coloured with tracer value c, where DG_Adapt_4 scheme is used.

Figure 10 .
Figure 10.Case three -swirling deformation: initial distribution and velocity field.

Figure 11 .
Figure 11.Case three -swirling deformation: comparison of the analytical solution with the results from different schemes using almost the same number of nodes N ≈ 15 000, at t = 3T /20, where T = 2.6376.

Figure 12 .
Figure 12.Case three -swirling deformation: comparison of the analytical solution with the results from different schemes using almost the same number of nodes N ≈ 15 000, t = T /5, where T = 2.6376.

Figure 13 .
Figure 13.Case three -swirling deformation: the evolution of the adaptive mesh coloured with tracer value c, where DG_Adapt_128 scheme is used.

Figure 15 .
Figure 15.Case three -swirling deformation: the evolution of (a) number of nodes, (b) max local and (c) integral of CFL number for CV_Adapt_128 schemes.

Figure 19 .
Figure 19.Case four -power plant plumes: the evolution of the adaptive mesh coloured with SO 2 concentrations (µg m −3 ), using the CV_Adapt scheme.

Figure 20 .
Figure 20.Case four -power plant plumes: the evolution of 3-D plumes visualization, surface SO 2 concentrations (µg m −3 ) and the corresponding adaptive mesh.