Semi-Lagrangian Subgrid Reconstruction for Advection-Dominant Multiscale Problems

We introduce a new framework of numerical multiscale methods for advection-dominated problems motivated by climate sciences. Current numerical multiscale methods (MsFEM) work well on stationary elliptic problems but have difficulties when the model involves dominant lower order terms. Our idea to overcome the assocociated difficulties is a semi-Lagrangian based reconstruction of subgrid variablity into a multiscale basis by solving many local inverse problems. Globally the method looks like a Eulerian method with multiscale stabilized basis. We show example runs in one and two dimensions and a comparison to standard methods to support our ideas and discuss possible extensions to other types of Galerkin methods, higher dimensions and nonlinear problems.


Introduction.
1.1. Motivation and Overview. Simulating complex physical processes at macroscopic coarse scales poses many problems to engineers and scientists. Such simulations strive to reflect the effective behavior of observables involved at large scales even if the processes are partly driven by highly heterogenous micro scale behavior. On the one hand hand resolving the microscopic processes would be the safest choice but such a strategy is prohibitive since it would be computationally expensive. On the other hand microscopic processes significantly influence the macroscopic behavior and can not be neglected.
Incorporating micro scale effects into macro simulations in a mathematically consistent way is a challenging task. There exist many scenarios in different disciplines of science that are faced with such challenges. In fully coupled paleo climate simulations, i.e., climate simulations over more than hundred thousand years a typical grid cell has edge lengths around 200 kilometres and more. Consequently, subgrid processes such as heterogeneously distributed and moving ice shields are not or just insufficiently resolved [32]. These subgrid processes are usually taken care of by so-called parameterizations. One can imagine this as small micro scale simulations that are then coupled to the prognostic variables such as wind speed, temperature and pressure on the coarse grid (scale of the dynamical core). This coupling from fine to coarse scales is being referred to as upscaling and is unfortunately often done in rather heuristic ways. This leads to wrong macroscopic quantities like wrong pressures and eventually even to wrong wind directions and more undesired effects such as phase errors.
The increasing complexity of earth system models (ESMs), in general, demands for mathematically consistent upscaling of parametrized processes that occur or are modelled on very different scales relative to the dynamical core. A first example was already mentioned: In global climate simulations a sea ice model computes the average ice cover for each coarse cell. This information then enters heat fluxes between ocean and atmosphere since sea ice forms an interface between them and is hence modeled as an averaged diffusive process. It is known from homogenization theory that simple averaging of heterogenous diffusive quantities does not reflect the correct coarse scale diffusion [8].
Another, quite contrary, approach to perfom multiscale simulations is adaptive mesh refinement (AMR), see [7]. The idea of AMR is to initially run a coarse scale simulation and then assess the quality of the solution locally by means of mathematically or physically based a posteriori estimators. The coarse mesh is then refined in regions where the quality of the solution is not sufficient or coarsened where the solution is smooth. This approach is different from the above upscaling since micro scales are locally resolved -it is therefore sometimes referred to as downscaling. Such an approach is attractive and is being used in practice but it is not feasible in situations in which heterogeneities in the solution are expected to occur globally.
One early attempt to use multiscale methods are Eulerian Lagrangian localized adjoint methods (ELLAM), which constitute a space-time finite element framework, see [9,20], and form the basis of their multiscale version MsELLAM [36,10]. Such methods rely on operator splitting and basis functions are required to satisfy an adjoint equation. For a review see [34].
There exist many other multiscale methods. Homogenization is an originally analytical tool to find effective models of otherwise heterogenous models in the limit of a large scale separation [5,8,27]. Unfortunately it is often too difficult to find such an effective equation but there exist numerical techniques that aim at designing numerical algorithms to effectively capture the behavior of the solution to the unknown homogenized problem.
The heterogenous multiscale method (HMM) was pioneered by E and Enquist [12,37] and refers to a rich population of variants [3,18,17]. This method emphasizes important principles in the design of multiscale methods such as the choice of macroand micro-models and their communication. For reviews and further references see [2,11] and [1] for a discussion of the HMM with special focus on PDEs.
Variational multiscale methods (VMM) have been developed in the 1990s by Hughes and collaborators, see [24,25]. The spirit of the method lies in a decomposition of the solution space and in the design of variational forms that reflect the relevant scale interactions between the spaces. Many versions of it exist and we point the reader to recent reviews [4,33] and the vast literature therein.
The present work is inspired by multiscale finite element methods (MsFEM) which can be seen as a special variant of the VMM. The idea of this method is to introduce subgrid variations into basis functions and can be dated back to works by Babuška, Caloz and Osborn [6] and shares ideas with the partition of unity method [31]. The MsFEM in its current form was introduced in [23,22,14]. The essential idea of the method is to capture the local asymptotic structure of the solution through adding problem dependent bubble correctors to a standard basis and use these as ansatz and trial functions. Many variations of this method exist and refer the reader to [13,16] for a review.
Many of the afore mentioned methods have the advantage that they work well for elliptic or parabolic problems and that they are accessible to an analysis. The difficulty in many applications on the other hand is their advection or reaction dominated character, i.e., the dynamics is often driven by low order terms. This poses major difficulties to numerical multiscale methods. Multiscale finite element methods naively applied will not converge to any reasonable solution since basis functions will exhibit artificial boundary layers that are not present in the actual physical flow. Ideas to tackle this problem are based on combining transient multiscale methods with Lagrangian frameworks [35] or with stabilization methods for stationary problems, see [28] for an overview. A HMM based idea for incompressible turbulent flows can be found in [29]. For a method based on the VMM see [30].

Contribution.
Our main contribution is a framework of numerical methods for advection-dominated flows which by reconstructing subgrid variations on local basis functions aims at reflecting the local asymptotic structure of solutions correctly. The idea combines multiscale Galerkin methods with semi-Lagrangian methods by locally solving an inverse problem for the basis representation of solutions that is adapted to the actual flow scenario. We demonstrate the idea on one and twodimensional advection-diffusion equations with heterogenous background velocities and diffusivities in both non-conservative and conservative form.
Applying standard MsFEMs directly leads to failure since boundary conditions (BCs) for the modified basis functions can not be prescribed arbitrarily. Suitable BCs on the subgrid scale are an essential ingredient for many multiscale methods. A wrong placement of Dirichlet BCs, for example, leads to boundary layers in basis functions that are not there in the real large scale flow, i.e., the way information propagates in advection-dominant flow needs to be respected and not artificially blocked. Other choices of BCs usually complicate the enforcement of conformity or obfuscate additional assumptions on the problem structure (such as local periodicity).
Conformal MsFEM techniques for advection-dominated tracer transport were already explored in our previous work in one spatial dimension on a transient advectiondiffusion equation [35]. The finding is that one has to follow a Lagrangian point of view on coarse scales such that flow is "invisible". On fine scales one can then simplify Lagrangian transforms in order to make advective effects locally milder without going to a fully Lagrangian setting. This amounts to prescribing Dirichlet BCs on coarse flow characteristics and has the effect that basis functions do not develop spurious boundary layers not present in the actual flow. While this work gave some useful insights it is unfortunately not feasible for practical applications since it suffers from several weaknesses. First, it is not directly generalizable to higher dimensions. Secondly, it needs assumptions on the background velocity that are not necessarily fulfilled in practical applications to ensure that coarse scale characteristics do not merge.
In order to circumvent these problems we suggest here a new idea based on a semi-Lagrangian framework that locally in time constructs a multiscale basis on a fixed Eulerian grid. The construction is done in a semi-Lagrangian fashion on the subgrid scale whereas the macroscopic scale is conveniently treated fully Eulerian. This is in complete contrast to our previous work but still respects that information in advection-dominant flows is "mostly" propagated along flow characteristics.
The construction of the basis in each cell at time t n+1 is done by tracing each Eulerian cell back to time t n . This yields a distorted cell. A basis on this distorted cell is then reconstructed by solving an inverse problem for the representation of the global solution at the previous time step. Then the local representation of the solution on this cell is propagated forward in time to get a basis on the Eulerian grid at t n+1 instead of propagating the solution itself.
All reconstructions of the basis in coarse cells are independent from each other and can be performed in parallel. The global time step with the modified basis is not critical since the coarse problem is small and since we can use algebraic constructs to make the assembly of the global problem efficient. Note that although we formulate the algorithm globally in an implicit form it can well be formulated explicitly which will allow for computational efficiency for the global step.
Our new approach has several advantages. First, we can effectively incorporate subgrid behavior of the solution into the multiscale basis via the solution of local inverse problems. Secondly, the method can handle advection-dominated flows in a parallel fashion and, furthermore, the idea works in any dimension. Numerical tests show that it is accurate in both L 2 and H 1 since it represents subgrid variability correctly. The method can handle problems that involve an additional reaction term. This consequently includes conservation problems. Furthermore, the idea is generic and can be used for vector valued problems and problems that involve subgrid information coming from actual data.
2. The Semi-Lagrangian Basis reconstruction in One and Two Spatial Dimensions. In this section we will outline our ideas on an advection-diffusion equation (ADE) with periodic boundaries as a model problem. We will outline all ideas for didactical purposed first in d = 1 and in d = 2 dimensions on where c δ (x, t) is the background velocity, A ε (x, t) is a matrix-valued diffusivity, and f and u 0 are some smooth external forcing and initial condition. We will use bold letters for the vectors and tensors independent of the dimension.
The indices δ > 0 and ε > 0 indicate that both quantities may have large variations on small scales that are not resolved on coarse scales H > 0 of our multiscale method. We will also work locally on a scale h H that can resolve the variations in the coefficients. Furthermore, we assume that c δ A ε (see remark below) and that c δ is well-behaved, for example C 1− in space and continuous in time. This assumption is often satisfied in practice for example in climate simulations where c δ is given at the nodes of a coarse grid. Depending on the application variations in c δ may be resolved on scale H or not. The diffusivity tensor is assumed to be positive definite (uniformly in ε and point-wise in x) with A ε ∈ L ∞ t (L ∞ x ) d×d , and is often derived from parametrized processes such as a varying sea ice distribution, convection, topographical features, or land use patterns. Note, that (2.1) does not conserve the tracer u in contrast to equation (2.2) provided f = 0.
Remark 2.1. Advection-dominance of a flow (what we sloppily expressed by c δ A ε ) is usually expressed by a dimensionless number -the Péclet number Pe which is essentially the ratio between advective and diffusive time scales. There exist several versions of this number [26]. Since for large variations of the coefficients on the subgrid scale advection-dominance is a very local property we need be more precise with what me mean by that. Here we assume that Pe is high on average, i.e., Pe = c δ L 2 L Aε L 2 where L is a characteristic length. We take L to be the length of the computational domain.
Although the general idea of reconstructing subgrid variability in basis functions is not substantially different for (2.1) and (2.2) there are some differences in the propagation of local boundary information described below that strongly influence the accuracy. We will explicitly describe these differences.
We start with outlining our method in one dimension since the idea is very simple and avoids complications that arise in higher dimensions. The idea is to represent non-resolved fine scale variations of the solution locally on a set of non-polynomial basis functions in each coarse cell. That is we fix two (Eulerian) meshes: A coarse mesh T H of width H and on each cell K ∈ T H of the coarse mesh we have a fine mesh T K h of width h H (we assume that h does not depend on K). On the coarse mesh we have a multiscale basis ϕ H,ms where N H is the number of nodes of T H . This basis depends on space and time and will be constructed so that we will obtain a spatially H 1 -conformal (multiscale) finite element space. Note that this is a suitable space for problems (2.1) and (2.2) since they have [15]. The initial condition u 0 can therefore be assumed to be in L 2 (T d ) but in later experiments we will choose it to be smooth. The fine mesh on each cell K ∈ T H is used to represent the basis locally, see Figure 1.
Standard MsFEM methods used for porous media flows are designed in such a way that the coarse scale basis solves the PDE model locally with the same boundary conditions as standard FEM basis functions. Note that the global coarse scale MsFEM solution is a linear combination of modified local basis functions that do resolve the fine scale structure induced by the problem's heterogenities and therefore do represent the asymptotics correctly. This works for stationary elliptic problems and for parabolic problems as long as there is no advective term involved or even dominant. The reason is that advective terms as they appear in our model problems prevent a basis constructed by a standard MsFEM technique to represent the correct asymptotics. This is since flow of information is artificially blocked at coarse cell boundaries. Due to the incorrect boundary conditions for the multiscale basis artificial steep boundary layers form in basis functions that do not exist in the actual global flow. For transient problems another difficulty ist that the local asymptotics around a point depends on the entire domain of dependence of this point and hence subscale information represented by a base function must contain memory of the entire history of the local domain of that base function.
A first attempt to bypass these difficulties was to pose boundary conditions for the basis on suitable space-time curves, across which naturally no information is propagated on the coarse scales. Such curves are Lagrangian paths. Therefore, a Lagrangian framework on the coarse scale and an "almost Lagranian" setting on fine scales was proposed by us earlier [35]. Unfortunately, this method is not feasible since it can not be generalized to higher dimensions and needs restrictive assumptions on the flow field.
Nonetheless, the results of [35] suggest that a coarse numerical splitting of the domain should correspond to a resonable physical splitting of the problem. Instead of a fully Lagrangian method on coarse scales semi-Lagrangian techniques to build the basis can circumvent the difficulties of Lagrangian techniques. But these are only local in time and therefore they do not take into account the entire domain of dependence of a point. We show how to deal with this in the following. First, we start with the global problem.
2.1. The global time step in 1D/2D. Suppose we know a set of multiscale basis functions in a conformal finite element setting, i.e., we approximate the global solution at each time step in a spatially coarse subspace V H (t) ⊂ H 1 (T d ) in which the solution u is sought (almost everywhere). We denote this finite-dimensional subspace as First, we expand the solution u H (x, t) in terms of the basis at time t ∈ [0, T ], i.e., Then we test with the modified basis and integrate by parts. Therefore, the spatially discrete version of both problem (2.1) and (2.2) becomes the ODE for (2.1) and The mass matrix is given by f H (t) contains forcing and boundary conditions and the initial condition u H,0 is the projection of u 0 ∈ L 2 (T d ) onto V H (0). Note that (2.5) contains a derivative of the mass matrix: known can be reconstructed unknown This is necessary since the basis functions depend on time and since we discretized in space first. The reader will notice that using Rothe's method of lines, i.e., discretizing in time first, it is not clear what basis function to use for testing and this a priori leads to a different linear system. For the time discretization we simply use the implicit Euler method. The discrete ODE then reads Other time discretization schemes, in particular, explicit schemes are of course possible. For didactic reasons we choose to present the algorithm in an implicit version. We will come back to that later. The next step is to show how to construct the multiscale basis. Convention. In both 1D and 2D we use bold letters like x, A to express potentially vector valued quantities although they would be scalar in one dimension in all formulas and figures. We do so since we would like to avoid confusing the reader with different notations that are merely due to different dimensionality although the basic ideas are the same. Difficulties that arise in more than one dimension will be pointed out explicitly. For the sake of consistency we also use ∇ instead of ∂ x in one dimension. Quantities marked with a tilde likex signalize (semi-)Lagrangian quantities.

The Reconstruction Mesh in
1D/2D. Our idea combines the advantage of both semi-Lagrangian and multiscale methods to account for dominant advection. The reconstruction method is based on the simple observation that local information of the entire domain of dependence is still contained in the global solution at the previous timestep. This can be used to construct an Eulerian multiscale basis: we trace back an Eulerian cell K ∈ T H at time t n+1 on which the solution and the basis are unknown to the previous time step t n . This gives a distorted cellK over which the solution u n is known but not the multiscale basisφ i , i = 1, 2.
In order to find the points where transported information originates we trace back all nodes in T K H from time t n+1 to t n . For this one simply needs to solve an ODE with the time-reversed velocity field that reads for each x l and then takex l =x l (−t n ), see Figures 3 and 2 for an illustration. This procedure is standard in semi-Lagrangian schemes and can be parallelized.

Basis Reconstruction in 1D.
After tracing back each point x l of K ∈ T H to its originx l in a distorted coarse cellK ∈ T H we need to reconstruct a local representation of the (known) solution u n onK: wherex j andx j+1 are the boundary points ofK. In one dimension one can of course choose a representation using the standard basis of hat functions but this would not incorporate subgrid information at step t n at all. Even in the forward advection step (explained below) that is being carried out in the next step subgrid information on the basis will not contain the correct subgrid variability because the information contained in the basis does not take into account the entire domain of dependence of K. We choose to bypass this problem by solving an inverse problem for the basis to adjusts the representation. The idea is to fit a linear combination of the basis locally such that u n is optimally represented, i.e., we solve The functions R i : C 0 (K) → R are regularizers weighted by positive numbers α i ∈ R. A simple regularizer that we found useful in one spatial dimension is a penalization of the deviation of the modified basis function from the standard linear basis function with the same boundary values in the quadratic mean, i.e., we use whereφK ,i denotes the t-th standard (linear) basis onK. This amounts to a system of linear equations that can be computed explicitly. In a spatially discrete version this system will be small and cheap to solve. A suitable choice of a regularizer depends on the problem at hand. For problems in two dimensions we choose a different one (see below). Figure 4 illustrates the effect of a local reconstruction of a basis compared to a representation with a standard basis.

Figure 5
The basis reconstructed according to (2.13) at time t n is propagated forward to time t n+1 according to (2.15) or (2.16).
reconstructed basis reconstructed basis after propagation 2.4. Basis Propagation in 1D. After having reconstructed a suitable basis on each coarse cellK we have an H 1conformal basis. This basis, however is a basis at time step t n and does not live on the coarse Eulerian grid T H that we initially fixed. The step to take now is to construct a basis at t n+1 on T H . This is done similarly to [35], i.e., we evolve the basis according to the model at hand with a vanishing external forcing. Note, however, that we compute the basis at t n+1 along Lagrangian trajectories starting from t n , i.e., we need to transform the original model. Equation (2.1) becomes (2.16) Note that these evolution equations are solved onK, i.e., on the element K ∈ T H traced back in time. Advection is "invisible" in these coordinates. The end state Trace back each node in K one time step from t n+1 to t n according to (2.11) 10 Reconstruct the optimal basis representation of u n (x)|K according to (2.13) 11 Propagate the optimal basis forward onto K according to (2.15)  Postprocess the solution 17 Return 18 end ϕ K,i (x, t n+1 ) onK can then be transformed onto the Eulerian element K ∈ T H to obtain the desired basis function ϕ K,i (x, t n+1 ) ∼ ϕ n+1 K,i (x) at the next time step. Corresponding basis functions in neighboring cells can then be glued together to obtain a modified global basis ϕ H,ms i , i = 1, . . . , N H . This way we get a basis of a subspace of H 1 that is neither a partition of unity nor is it necessarily positive. Nonetheless, it is adjusted to the problem and the data at hand. The propagation step is illustrated in Figure 5.
Using our method we reconstruct and advect the representation of the global solution first and then the solution itself using the modified representation. The global step is completely Eulerian while the local reconstruction step is semi-Lagrangian in contrast to [35] where the global step is Lagrangian and and the local step is "almost"-Lagrangian. The essential steps described above are summerized in Algorithm 2.1. Note that the steps to reconstruct the multiscale basis are embarrassingly parallel and consist of small local problems.

Basis Reconstruction in 2D.
The basis reconstruction step in two dimensions is slightly different. We separate the description of the procedure from the reconstruction in one dimension to make the reader aware of difficulties that arise when increasing the space dimension. Keeping the difficulties in mind one can then easily generalize the procedure even to 3D. We will comment on that again later. It is worth mentioning that these difficulties do not need to arise in non-conformal settings but since we agreed on a conformal method we will be consistent here.
As the reader may have guessed to construct a H 1 -conformal basis we need to ensure that the reconstructed global basis is continuous across coarse cell boundaries. This can be achieved by first looping over all coarse edges of the traced back mesh and reconstructing the solution at the previous time step t n with a basis representation on each edge, i.e., we solve first an inverse problem of the kind (2.13). (2.17) Note that the regularizer (2.14) needs to be replaced since the edge can be curved (it is unclear whatφ 0 K,i in this case is). We use a harmonic prior for the reconstructed basis with a low weight α i in (2.13). The operator ∆ g(Γ) denotes the Laplace-Beltrami operator induced by the standard Laplace operator with the trace topology of the respective (traced back) edgeΓ, i.e., g(Γ) is the metric tensor. We pass on providing details here. The edge reconstructed basis then serves as boundary value for the cell basis reconstruction. The optimization problem to solve on the traced back cellK then reads x j is corner point of (traced back) cellK ϕK ,j |Γ l = ϕΓ ,j l .
(2.19) However, the essential task is to ensure conformity of the global basis by first reconstructing representations on all edges and then inside the cellsK. In three dimensions it would be necessary to first reconstruct edge then face and only then the interior representations. This sounds potentially expensive but it is embarrassingly parallel since all reconstructions are independent.
2.6. Basis Propagation in 2D. Again, for didactical reasons we make the reader aware of differences to the propagation of the basis in contrast to one spatial dimension. As in Section 2.5 we need to preserve conformity of the global basis. This can easily be done by fixing the boundary values in the propagation step of the basis functions similar to the propagation in 1D, see Section 2.4.
However, this strategy can result in numerical errors from two sources. Recall that the goal of our multiscale method is to represent the local asymptotic structure of the solution correctly by putting subgrid variations into the basis. This can be interpreted as adding a non-polynomial corrector function in each coarse cell to a standard FEM basis which changes the boundary conditions. These boundary conditions are themselves subject to evolution and keeping them fixed in time will not reflect any system intrinsic dynamics.
The first source of a numerical error is a resonance error. This error usually becomes dominant if the scale of oscillations in the diffusion term is close to the scales resolved by the coarse grid. This is well-documented in the literature on multiscale FEMs for stationary elliptic problems that are usually not dominated by lower order terms. Note, that in practical problems it is often unclear if a scale separation exists or it can not be quantified which makes it difficult to identify a resonance regime [19,23].
Another source for not capturing the asymptotic structure with fixed boundary values even if there was a scale separation between coarse mesh and diffusion oscillations (which can not be guaranteed) appears in the conservative case (2.1). This is due to the reactive term (2.20) ∇ ·c δ ϕ K,i that appears also in (2.16). This term compensates divergence, i.e., it triggers a reaction on the boundary wherever the divergence of the background velocity does not vanish. Essentially that means that this term is very local but can be quite strong, i.e., it is of the order of δ −1 if c is of order 1 but oscillates on a scale of order δ. This term must be taken into account when evolving a reconstructed basis. Fortunately, both errors can be dealt with by employing two strategies. The first one is oversampling, i.e., the local reconstruction and basis propagation on a coarse cell K is performed on a slightly larger domain than K os ⊃ K. After that the basis is being truncated to the relevant domain K. Using this strategy is generally possible but results in a non-conformal technique and hence we will not use it here [19]. Instead we will solve a reduced evolution problem on the edges that predicts the boundary values and then evolve the reconstructed interior basis for which we then know the boundary values in time. The edge evolution strategy is in general only necessary in the case of a missing scale separation for the diffusion as explained above or in the conservative case. Note that both errors do not appear in one dimension. We summarize this discussion below.
Strategy without Edge Evolution.. In this case we keep the boundary values of the basis that we reconstructed at the previous time step fixed. Similarly to (2.15) this amounts to solving the evolution problem for the i-th basis function.
Strategy with Edge Evolution.. This method requires a prediction of the boundary values. The prediction is done by solving a reduced evolution problem Reconstruct the optimal basis representation of u 0 (x)| K according to (2.13) 24 end /* Now the time stepping starts. We compute the solution at t n+1 . */ 25 for n = 0 to n ≤ N steps do 26 for K ∈ T H do in parallel 27 Trace back each node in K one time step from t n+1 to t n according to (2.11) 28 Reconstruct the boundary condition of optimal basis representation of u n (x)|K according to (2.17), see Section 2.5 29 Reconstruct the optimal basis representation of u n (x)|K according to (2.19), see Section 2.5 30 if Problem is not in conservation form then 31 Propagate the optimal basis forward onto K according to (2.21) andÃ red ε only contain the tangential parts of their counterpartsc δ andÃ ε . Once this evolution problem is solved for each edgeΓ the boundary values for the full evolution problem onK are known and one can solve d dt to get the i-th basis.
In both cases note again that the evolution equations (2.21) and (2.23) are solved onK. The end state ϕ K,i (x, t n+1 ) onK can again be mapped onto the Eulerian element K ∈ T H to obtain the desired basis function ϕ K,i (x, t n+1 ) ∼ ϕ n+1 K,i (x) at the next time step.
The described modifications for 2D problems are summarized in Algorithm 2.2.
3. Numerical Examples in 1D. We will show several 1D examples in a nonconservative and conservative setting according to (2.1) and (2.2), respectively. For all 1D tests we use a Gaussian with variance σ = 0.1 centered in the middle of the domain T 1 , i.e., µ = 0.5. The end time is set to T = 1 with a time step δt = 1/300. We show our semi-Lagrangian multiscale reconstruction method (SLMsR) with a coarse resolution H = 10 −1 in comparison to a standard FEM with the same resolution and high order quadrature. As a reference method we choose a high-resolution standard FEM with h ref = 10 −3 .
For the multiscale method we choose a fine mesh T K h with h = 10 −2 in each coarse cell K ∈ T H . Test 1.. We start with four tests showing the capability of the SLMsR to capture subgrid variations correctly. Note that the coarse standard FEM has as many cells the the SLMsR has coarse cells and that it does not capture subgrid variations in the following tests. The resolution for the reference solution resolves all subgrid variations but the reader should keep in mind that practical applications do not allow the application of high-resolution methods. The coefficients    The first test demonstrates how the SLMsR behaves when all coefficients of the equation are resolved by the SLMsR, i.e., Hh ≤ ε, δ but not by the low resolution FEM that has the same resolution as the SLMsR's coarse resolution, i.e., H ε, δ. For the coefficients we chose A ε (x, t) = 10 −3 + 9 · 10 −4 cos(174πx) . Error plots of the results are shown in the first row of Figure 8. The plots indeed show that the SLMsR captures the reference solution much more accurately than the standard FEM in both L 2 and H 1 even for small H. Only when decreasing H to a resolution that resolves all coefficients the FEM is able to represent the solution correctly and starts converging. The reader should keep in mind that in real world applications the standard FEM will be too expensive to compute while the SLMsR allows for massive parallelization. This is the regime that is of practical interest. Increasing the fine resolution in each coarse cell on the other hand does not improve  A ε (x, t) = 10 −3 + 9 · 10 −4 cos(16πx) .

(3.7)
and refined with the sequence H = 1 16 , 1 32 , 1 64 , 1 128 while we we fixed the number of fine cells in each coarse cell to n f = 32. Note that we start with a coarse resolution regime that almost resolves the data. The error plots for the SLMsR and the FEM indeed indicate that as H decreases the SLMsR does not dramatically increase its accuracy while the FEM with linear basis functions increases its accuracy according to the standard estimates. This can be observed in the second row of Figure 8. The relative errors at t = 1/2 and t = 1 for both test problems (3.6) and (3.7) are summarized in Tables 2 and 3 in the appendix together with the estimated order of convergence (EOC) in space. We recall that the EOC can be computed as follows: Let e H : R n + → R + be an (error) function that maps a high dimensional vector to a positive real for each H > 0. Then the EOC w.r.t. H can then be computed as Similar results in terms of accuracy can be obtained in case of a conservative equation. We pass on showing this here. Test 3.. This test shows an example where both diffusion and background velocity are generated randomly. We intend to show an example of the SLMsR behaves when data is involved that does not exhibit a clear scale separation. For this we initially generate fixed mesh based functions with random nodal coefficients. In each mesh cell the functions are interpolated linearly. Note that this is not to simulate a sampled stochastic process. We simply intend not to create any scale or symmetry bias when constructing coefficient functions. The results look appealing and show a clear advantage of the SLMsR, see Figure 9.

Numerical Examples in 2D.
This section is to experimentally demonstrate that the SLMsR can handle conservative and non-conservative advection-diffusion equations with dominant advection term in higher dimensions. All tests are being carried out on the torus T 2 (periodic unit square) in the time interval t ∈ [0, 1] and use the same initial condition. As initial value we chose a normalized super-position of two non-isotropic Gaussians All tests of the SLMsR are are perfomed on a coarse unstructured uniform triangular Delaunay mesh with n c = 62 coarse cells, i.e., for our triangulation H ∼ 0.3 (mean diameter of circumcircle of a cell). We compare the SLMsR to a standard low resolution FEM with the same resolution and to a standard high resolution FEM with approximately n f = 63K cells. To get a fine mesh on each coarse cell of the SLMsR we created a triangulation such that the sum of all fine cells over all coarse cells is approximately n f to get a fair comparison of the SLMsR to the low resolution standard FEM with respect to the reference solution that resolves all coefficients involved. Test 1.. Here we test our multiscale reconstruction method with a solenoidal field c δ described by the stream function (4.3) ψ(x, t) = sin(2π(x 1 − t)) sin(2πx 2 ) so that (a) (b) (c) Figure 11 Background velocity for Test 1 and Test 3. Four vortices moving through the domain from left to right and come back to their starting points at T = 1.
This background velocity desribes four vortices moving in time through the (periodic) domain from left to right and get back to their starting point at T = 1. Note that this velocity field involves scales that are resolved by the coarse mesh as well as scales that are not resolved, see Figure 11. Note that since ∇ · c δ = 0 equation (2.1) and (2.2) are (analytically) identical and hence we only solve (2.1).
The diffusion tensor is chosen to be (4.5) Note that in this case advection dominance is a local property and Péclet numbers are ranging from Pe = 0 to Pe ∼ 6 · 10 6 . Snapshots of the solutions are shown in Figure 10. It can be observed that the low resolution FEM does not capture the effective solution well since it diffuses too strongly while the SLMsR reasonably caputeres the effective behavior of the solution and even the fine scale structure. This can be seen also in the error plots, Figure 16.

Figure 12
Background velocities for Test 2a and Test 2b given by (4.6) and (4.7), respectively. Both fields are divergent.
Test 2a.. For this test we use two divergent background velocities and hence we split it into two parts to distiguish between the diffusive transport problem (2.1) and the conservation law (2.2). We start with showing how the SLMsR behaves on equation (2.1). The velocity field is a modification of (4.4) and is given by .
while the diffusion tensor is also given by (4.5). Here the standard solution does not converge to any reasonable approximation of the effective solution and no valuable understanding of the dynamics can be drawn from it. The SLMsR shows a suprisingly good approximation of the reference solution, see Figure 13. Test 2b.. Here we solve equation (2.2), i.e., a conservation problem. The divergent velocity field for this test describes regions of fast and slow flow with two separatrices across which there is no flow moving in time from left to right once through the domain during the time interval. The field is given by (4.7) c δ (x, t) = 2π 5 sin(2π(x 1 − t)) 2 cos(π(x 2 − 0.5)) 2 2 sin(2π(x 1 − t)) cos(2π(x 1 − t)) cos(π(x 2 − 0.5)) 2 .   Both velocity fields for Test 2a and Test 2b are shown in Figure 12. The diffusion tensor is again given by (4.5). The difference to the previous test is that we have quite a large additional reaction term to handle. This will be taken care of in the basis representation as described in Sections 2.5 and 2.6 as well as in the gloabl weak form, see equation (2.7). Note that the difficulty we focus on here is not the actual conservation (although this is an important issue) but the fact that we do have another lower order (reaction) term in the equation with a multiscale character. Fortunately this term is very local and is handled nicely by the SLMsR since it scales with its magnitude. The results are shown in Figure 14 where we also show how the quality of the SLMsR solution is negatively affected when the reconstructed basis boundary conditions are not propagated. As the reader can see the SLMsR without edge evolution shows a slight grid imprinting of the coarse mesh. This effect can be more pronounced depending on the strength of the reactive term since it is due to keeping the boundary conditions fixed in time in the local basis evolution step. This in turn increases the error in the evolution of the asymptotic structure near coarse grid boundaries due to the additional reaction term.
We pass on showing a comparison of all snapshots to the low resolution FEM here (which does not behave well) for the sake of brevity. Error plots including the low resolution FEM are again summarized in Figure 16.
Test 3.. This test has a similar intention as Test 3 in Section 3. We intend to show a reasonable behavior of the SLMsR if the coefficients do not show a clear scale separation on equation (2.1). For this we again used (4.4) as background velocity and Figure 17 Magnitude of randomly generated (but fixed) diffusion tensor for Test 3. generated a diagonal diffusion tensor that is constant in each cell of a mesh generated to represent solely the diffusion. The cell based constants are random and are fixed at the beginning of the simultion. The diffusion tensor was then scaled to contain values in the range of 10 −5 to 10 −1 , see Figure 17. Snapshots of the solution can be found in Figure 15 and error plots in Figure 16.
5. Summary and Discussion. In this work we introduced a multiscale method for advection-diffusion equations that are dominated by the advective term. Such methods are of importance, for example, in reservoir modeling and tracer transport in earth system models. The main obstacles in these applications are, first, the advection-dominance and, secondly, the multiscale character of the background velocity and the diffusion tensor. The latter makes it impossible to simulate with standard methods due to computational constraints while simulating using standard methods with lower resolution that does not resolve variations in the coefficients leads to wrong effective solutions.
Our idea to cope with these difficulties is inspired by ideas for semi-Lagrangian methods, ideas based on "convected fluid microstucture" as described in [21], inverse problems and multiscale finite elements [13]. At each time step we reconstruct fine scale information from the solution at the previous time step. This fine scale information enters the local representation of the solution in each coarse cell, i.e., it is added as a corrector to the local basis such that the basis representation is optimal in some sense. Due to advection-dominance the optimal basis needs to be reconstructed in its departure area, i.e., we locally trace back information like in a semi-Lagrangian method for the reconstruction. The reconstruction is performed by solving an inverse problem with a suitable regularizer. It constructs a basis that does not constitute a partition of unity (PoU) and is tailored for the actual problem at hand. The idea of adding prior knowledge about the solution to a local representation in PoU methods, however, is similar, see for example [31,17]. After reconstructing the basis at the previous time step it is evolved with suitable boundary conditions to the time level, where it is needed; i.e., we evolve the local representation of the solution rather that the solution itself. Note that the global framework of the SLMsR is completely Eulerian while only the local reconstruction step in each coarse cell is semi-Lagrangian.
The SLMsR also shares shallow similarities with stabilized methods such as streamline upwind Petrov-Galerkin (SUPG) methods in which "advected" basis functions are used as correctors for test functions to add a small diffusion in the direction of the flow. This is however not the same since our basis also corrects for the fine scale asymptotic structure induced by diffusion and other low order terms.
Numerical experiments that we carried out in one and two dimensions for conservative and non-conservative advection diffusion equations show a clear advantage in terms of accuracy for the SLMsR at low (coarse) resolutions. We formulated the method globally as a conformal FEM to show how the idea works in a quite restrictive setting. Other frameworks are possible. A non-conformal FEM, for example, would put less restrictions on the basis reconstruction at the price of having to compute a boundary integral in the global assembly for each cell.
One of the main features of the SLMsR is its scalability: Although it sounds expensive to trace back each coarse cell, then solve an inverse problem and then solve a PDE at each time step (the so-called offline phase) we would like to point out that these local problems are independent and usually small and therefore the offline phase is embaressingly parallel, although we did not take advantage of that in our implementation yet. The global time step (online phase) also consists of a small problem and matrix assembly procedures can be made very efficient by using algebraic expressions, see [23,13].
We would like to further emphasize the flexibility of the SLMsR. Here we presented an implicit version but explicit time stepping is possible. Also, the method can be transferred to higher dimensions as well as it can be extended to deal with advectiondiffusion-reaction problems. Furthermore, the use of inverse problems in the local steps to adjust the basis makes it generally possible to incorporate knowledge coming from measurement data. For this a thorough understanding of the data is necessary (as for any other assimilation method). We would like to explore that opportunity in the future.