Project Topics
Get Complete Project Material File(s) Now! »
Methods for the numerical resolution of the advection equation
In its general form, the advection equation does not have a straightforward analytical solution. Even though analytical solutions may be obtained for some particular cases and under specific restrictions (e.g. [Perez Guerrero et al., 2009], [Zoppou and Knight, 1997], [ECMWF, 2002]), these methods are not of general and practical use. As such, there is a strong need for numerical solutions of discretized forms of the advection equation, either in the advective (Eq. 1.1) or flux form (Eq. 1.2).
Precursors of modern calculus, Gottfried von Leibniz and Isaac Newton, were the first to address the fact that some differential equations lack analytical solutions. Consequentially, they were also the first to introduce different approaches to deal with this problem. While Leibniz resorted to geometry to find new forms of analytical solutions by introducing novel trigonometric functions, such as the inverse tangent (Leibniz [1684, 1693]), these newly introduced functions only suited some specific forms of differential equations. Newton then discussed approximate solutions to some of these equations by series expansion in his treatise on differential calculus [Newton, 1736]. These ideas, along with with the work of Brook Taylor—the Taylor series—set the theoretical bases for the numerical solution of differential equations, which then awaited the era of computers. Focusing on equations whose solution depends on space and time, such as the advection equation, we can be more precise regarding the discretization process.
Usually—if we omit spectral methods based on Fourier (or spherical harmonic) expansions and the domain of mesh-free methods [Fasshauer, 2007]—obtaining approximate numerical solutions for a partial differential equation (PDE) such as the advection equation (Eq. 1.1 or 1.2) requires the domain of definition to be divided into small sub-domains. Accordingly, we suppose that the studied variable (in our case, the mixing ratio ) can be correctly approximated in each of the sub-domains by a low-order multivariate polynomial expansion (e.g. constant, linear, quadratic, etc.). Each of these elementary spatial subdomains is called a cell, and together they form a discretization of the domain called a mesh, or grid. Numerical methods can be designed to build approximations of the exact solution in each cell in such a way that the merging of these finitesimal solutions converges towards the exact solution of the PDE as the grid becomes finer. Such numerical methods are called numerical schemes.
Having discussed the discretization in terms of space, we now turn to time as the advection equation also incor-porates the time dimension. In such equations, both time and space are discretized, with one to three dimensions of space discretized following a mesh as discussed above, and one dimension of time discretized by a partition into finite intervals. The latter intervals often have a (conditionally) constant length called a time step.
Numerical strategies for solving the advection equation therefore involve two essential steps which we examine below: designing a mesh that fits the geometry of the domain, and designing the numerical schemes which permit to obtain approximate solutions of the equation for each mesh cell and each time step. In the present section, we first examine different existing strategies to build meshes covering an entire sphere, and then we move on to some of the main numerical schemes designed to obtain numerical solutions of the advection equation on these meshes.
Cartesian grid formulations on the sphere
Whether Rene´ Descartes was the first to propose the classical Cartesian coordinate system in his book La Geom´etrie´ [Descartes, 1637], or simply developed its use in localization therein, it came to carry his name. From this coordinate system derives a type of mesh called the Cartesian mesh or—more prevalently—the Cartesian grid. These grids are made up of the set of squares or rectangles obtained by splitting each of the two coordinate axes of the plane into intervals and—for each axis—constructing the lines parallel to the other axis, passing through interval boundaries. The same method is applied for the meshing of n-dimensional spaces as well.
Even though this type of mesh is not directly applicable to a sphere—which does not admit a Cartesian coordinate system, as its points cannot be described as a pair of coordinates along two orthogonal axes—meshes with a structure comparable to the Cartesian meshes of the plane can be built on the sphere. They are traditionally called Cartesian, or rectilinear, grids and can cover the entire sphere or a part of it. Rectilinear grids can be obtained readily from any map projection of the sphere: if we choose a map projection associating a point X on the sphere to a point (x1; x2) of the plane5, then any Cartesian mesh of the plane can be mapped onto a corresponding rectilinear mesh covering the sphere (or a part of it). The properties of such rectilinear meshes depend on the properties of the map projection.
Variable-resolution Voronoi-mesh generation
Similarly to the quasi-uniform Voronoi meshes, a variable resolution Voronoi mesh is built upon a Delaunay triangulation. However, to account for the variable resolution, the starting point cannot, in this case, be a regular icosahedron. Instead, a set of points is randomly distributed on the surface of a unit sphere, with a higher probability in regions where higher resolution is sought. A Delaunay triangulation is formed using these points, and then the Voronoi dual is taken in the same way as in the previous case. Also, the target resolution is taken into account during the Lloyd iteration process [Du et al., 1999]. It is important to note that, with variable resolution, the Voronoi mesh is no longer hexagonal-pentagonal: at least some cells have seven edges or more.
Edge values of the mixing ratio
The mixing ratio at the midpoints of the cell edges is denoted by e. This edge value is needed in the numerical schemes detailed in Section 2.3. Its evaluation depends on the chosen numerical scheme, starting from the value at the respective cell center. This step of the numerical formulation of schemes is referred to as reconstruction [Dubey et al., 2015]. To transport the tracer field from one cell to one of its neighbouring cells, we reconstruct the tracer mixing ratio at the edge common to the two cells.
Discretized velocity field
The velocity field, u, on the other hand, is a vector field, and a different approach is used to discretize it. To each edge, of index e, we attach a real number, ue, which we define as the velocity vector—at the edge midpoint—projected onto the normal unit vector of the edge.
Discretized gradient
This quantity is useful for the second order numerical schemes requiring a discrete gradient (Section 2.3.3). This is the discrete counterpart of the continuous gradient of the continuous mixing ratio field. In our study, we go by the work of Dubey et al. [2015], who demonstrate that a first order approximation of the gradient is sufficient to maintain second order scheme accuracy. Following their first order approximation, a first-order estimate of the gradient, rv is obtained at a each dual cell (with index v) as the dual-cell-averaged gradient: Z 1 rv = r dA (2.10).
where Av is the area of the dual cell v, and Vv is the dual cell itself (i.e. a triangular subset of the sphere). The vertices of these dual cells correspond to the centers of the the primal cells, where the mixing ratio can be estimated with second-order accuracy as discussed above. Using the Green-Gauss theorem, the above area integral is converted into a line integral: ZVv r dA = Z nd (2.11). where is the perimeter of the vth cell, and n is the normal relative to the infinitesimal edge, d , directed outwards.
Discretizing the above expression using second-order-accurate estimates of the mixing ratio at primal-cell centers, we obtain the discrete gradient of the continuous scalar field at the dual cell v: 1 X (2.12).
where nv;e is the vector normal to the dual edge e, and mv;e = 1 is the factor ensuring the vector is directed out of dual cell v (analogous to the nk;e of primal cells). de is the length of this dual edge, v;e is the half sum of estimated at the extremities of the dual edge (which are cell centers) and v is the equal-weight average of at dual-cell vertices, which are also cell centers. We apply the Green-Gauss formula (2.11) to v, which has the same gradient as . This ensures that rv is zero whenever the mixing ratio is uniform. In a second step, rk is computed for each cell k by performing a dual-cell-area-weighted average of rv for all dual cells v intersecting the primal cell k.
One-dimensional transport schemes
We consider finite-volume schemes of first-, second- and third-order accuracy. The first-order scheme is the original Godunov scheme [Godunov, 1959], which the Voronoi upwind scheme (Section 2.3.3) generalizes to Voronoi meshes. The second-order scheme is the Van Leer method [Van Leer, 1977]. The SLFV scheme can be regarded as an extension of this scheme to Voronoi meshes. More specifically, the Van Leer method relies on a linear reconstruction in each cell and a semi-Lagrangian approach to time integration. As with SLFV, monotonicity is enforced by a slope-limiting technique.
The Piece-wise Parabolic Method (PPM) scheme [Colella and Woodward, 1984], [Carpenter et al., 1990] is also used in our study. The PPM scheme implements a semi-Lagrangian design similar to Van Leer [1977] but with a better accuracy due to a parabolic (instead of linear) reconstruction of in each cell. While formally third-order away from extrema in the advected field, the slope-limiting procedure presented in Colella and Woodward [1984] generates order-1 error terms in the immediate vicinity of extrema, preventing the scheme from formally reaching third-order accuracy [Colella and Sekora, 2008]. The PPM scheme—as described by Colella and Woodward [1984]—is therefore order-2, but generally more accurate than Van Leer [1977] (see, e.g., Mailler et al. [2021]). As discussed above, solvers for the two- and three-dimensional advection equations are obtained by combining these schemes with time splitting with the above-described strategies of Lie (order-1) and Strang (order-2).
Double cosine-bell distribution
The cosine-bell is a smooth distribution with a sinusoidal variation between a background value and a maximum value (Eq. 3.1). In its double form, two cosine-bell-shaped hills, are allocated to two distinct parts of the domain.
In 2D, cosine bells actually have a circular shape, with a maximum concentration at a central point, decreasing in a smooth sinusoidal manner to reach a minimum value at the edge, matching the background value. Following Lauritzen and Thuburn [2012], we use cosine bells characterized by a common radius R = 1=2 (in radians) and by the coordinates ( n; n) of their centers (n = 1; 2 for two cosine bells): (;)= 8 b + c h rn(R ; if rn( ; ) < R with n = 1; 2: (3.1) < ; ) b; otherwise: :
where b and c are constants, rn( ; ) is the great-circle distance from the center of cosine bell n, found as: rn( ; ) = arccos [sin n sin + cos n cos cos (n)] (3.2) and h(r) is the shape function : 1 h(r) = [1 + cos ( r)] : (3.3) b is the background mixing ratio, and b + c is the maximum mixing ratio. c; b 2 [0; 1] are chosen so that b + c = 1, and the mixing ratio, , remains in the interval [0; 1]. When computing entropy (Section 3.4.4), a background value of 0 is necessary to eliminate any contribution to entropy coming from regions of the full hexagonal-pentagonal mesh beyond the limits of the Cartesian domain. Values b = 0 and c = 1 are used in this case. This pair of values is also used to investigate the positivity of the schemes, which is a prerequisite for studying the preservation of non-linear tracer-relationships, using, again, the same values. Otherwise, we use b = 0:1 and c = 0:9.
Single cosine-bell distribution
A simplified form of the above distribution, it is represented by only one of the two cosine bells, the location of which can be varied. Such a distribution is particularly useful for tracking a maximum mixing-ratio throughout space. With only one peak, the location of the maximum is unique. Similar to other smooth distributions, the single cosine-bell is useful when investigating a property such as convergence, where sharp discontinuities would prevent a useful comparison of the various schemes(such as the slotted-cylinder distribution used by Nair and Lauritzen [2010]). In this case, we use the parameters set as b = 0:1 and c = 0:9.
Tilted solid-body-type periodic oscillation
Such a wind field is simple in its similarity to the ordinary solid-body field. However, it endows some challenges to the schemes through an inclination, while still avoiding strong deformations that increase numerical diffusion. A sinusoidal time-dependence guarantees the return of the scalar field to its initial position after a full period (Figure 3.3). These properties make this wind field a reasonable candidate for studying convergence, which requires a comparison between the numerical and exact solutions, the latter identical to the initial field distribution.
At the same time, by introducing a tilt of the wind field relative to the zonal direction, we reduce the advantage endowed to the Cartesian mesh by a solid-body wind-field aligned with the Cartesian grid. The tilt and wind speed are set in a way that ensures that no tracer mass is transported beyond the boundaries of the Cartesian domain.
To define this field, we begin with the classical solid-body rotation around the axis passing through the two poles. A tilt is made to the axis of rotation, such that it is skewed away from the pole axis by an angle . If we define the sphere of radius a as the surface x2 + y2 + z2 = a2 in a three-dimensional Cartesian system (x; y; z), this is represented by a skew from the z-axis, within the z and x plane, with no change along y. The tilted axis has the ( ; ) coordinates (0; 2 ), where is the tilt angle in radians. In terms of 3D (x; y; z) coordinates, the tilted axis can be expressed as follows: which falls back to the pole axis ez when the tilt is zero. At any point (a cos cos ; a cos sin ; a sin ) on the surface of the sphere, the wind-field stream-function is: ( ; ; t) = u0a cos 2T er e (3.11) t with: er e = a (cos cos sin + sin cos ).
Table of contents :
1 Introduction
1.1 Context
1.1.1 Sharp plumes
1.1.2 The advection equation
1.2 Methods for the numerical resolution of the advection equation
1.2.1 Cartesian grid formulations on the sphere
1.2.2 Triangular and polygonal meshes
1.2.3 Numerical schemes
1.2.4 Applications
1.3 Research questions
1.4 General outline
2 Discretization
2.1 Mesh generation and descriptions
2.1.1 Quasi-uniform Voronoi-mesh generation
2.1.2 Variable-resolution Voronoi-mesh generation
2.1.3 Voronoi and Delaunay mesh elements
2.2 Discrete representation of fields
2.2.1 Discrete degrees of freedom
2.2.2 Edge values of the mixing ratio
2.2.3 Discretized velocity field
2.2.4 Discretized gradient
2.3 Numerical solution of the advection equation
2.3.1 Time discretization
2.3.2 Method of lines: time-integration schemes
2.3.3 Method of lines: spatial schemes
2.3.4 Coupled time and space scheme
2.4 Slope limiting
2.5 Cartesian meshes
2.5.1 Cartesian-mesh generation and elements
2.5.2 Operator splitting
2.5.3 One-dimensional transport schemes
3 Two-Dimensional Test Cases
3.1 Tracer-concentration distributions
3.1.1 Uniform distribution
3.1.2 Double cosine-bell distribution
3.1.3 Single cosine-bell distribution
3.2 Wind fields
3.2.1 Solid-body rotation
3.2.2 Tilted solid-body-type periodic oscillation
3.2.3 Dual vortex with solid-body rotation (NL2010)
3.3 Limiting to the Cartesian band
3.4 Metrics
3.4.1 Stability
3.4.2 Monotonicity
3.4.3 Convergence
3.4.4 Numerical diffusion
4 Findings on the Numerical Mesh-Scheme Pairs in 2D
4.1 Stability
4.2 Monotonicity
4.3 Numerical diffusion
4.3.1 Entropy
4.3.2 Preservation of non-linear relationships between tracers
4.4 Convergence
4.4.1 Zonal solid-body-rotation convergence
4.4.2 Tilted periodic solid-body convergence
4.4.3 NL2010 low-shear convergence
4.4.4 NL2010 intermediate-shear convergence
4.4.5 NL2010 default-shear convergence
4.5 Partial conclusions
5 Puyehue-Cord´on Caulle 3D Simulation Results
5.1 Puyehue-Cord´on Caulle eruption
5.2 3D simulations
5.3 Trajectory of the plume
5.3.1 SO2 concentration at 330 hPa
5.3.2 Column-integrated SO2 concentration
5.3.3 Plume trajectory
5.4 Summary
6 Conclusion
6.1 Main results
6.2 Perspectives for chemistry-transport modeling
7 Sommaire en franc¸ais
7.1 R´ esultats principaux
7.2 Perspectives pour la mod´ elisation chimie-transport
