Cyril Voyant
email: cyrilvoyant@gmail.com

Fabrice Motte

Alexis Fouilloy

Gilles Notton

Christophe Paoli

Marie-Laure Nivet

Kalman filtering and classical time series tools for global radiation prediction

Keywords: Time Series forecasting, model free forecasting, Artificial Neural Networks, ARMA, Kalman filter

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

For nowcasting and short term forecasting of salar irradiation, the usual technics are based on machine learning predictions such as Artificial Neural Network (ANN) [START_REF] Voyant | Multi-horizon solar radiation forecasting for Mediterranean locations using time series models[END_REF], Support Vector Machines (SVM) [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF], AutoRegressive-Moving-Average (ARMA) models [START_REF] Hamilton | Time series analysis[END_REF], etc. A significant inconvenience of these methods is related to the large historic data set required during the training phase of the predictors; thus, in this work, we propose a simple methodology able to predict a global radiation time series without the need of historical data, making the method easily applicable for poor instrumented areas. We suggest to call these intuitive methods in the following "training-less" methods. The accuracy of these methods will be compared against other classical prediction methods, taking into account the time horizon of the prediction.

Data

The solar data used in the models are global horizontal irradiations (GHI) measured at the meteorological station of Ajaccio (Corsica Island, France, 41°55 N, 8°44 E, 4m asl) with pyranometers (CM 11 Kipp & Zonen). Data are related to an hourly basis from 1998 to 2009 (11 complete years).

Models

Data driven methods

-Linear model: AutoRegressive process (AR). In an AR model, the future value of a variable namely 𝐺𝐻𝐼 ̂(𝑡 + ℎ) is assumed to be a linear combination of several past observations (𝑡 -𝑖). In this study, the complexity of the model depends on the autoregressive order p which is optimized using the auto-mutual information factor.

-Non-linear model: neural network model (ANN). It is the predominant method in the domain of time series forecasting. Indeed, the availability of meteorological historical data databases and the fact that ANN are data driven approaches capable of performing a non-linear mapping between sets of input and output variables make this modelling tool very attractive. In the present study, the ANN model has been computed with the Matlab © software and its Neural Network toolbox. The Levenberg-Marquardt learning algorithm with a max fail parameter before stopping training equal to 3 was used to estimate the ANN model's parameters. The max fail parameter corresponds to a regularization tool limiting the learning steps after a characteristic number of predictions failures. Consequently it is a means to control the model complexity [START_REF] Paoli | Forecasting of preprocessed daily solar radiation time series using neural networks[END_REF].

3.2.

Training-less methods -The persistence. This naïve method is the most cost-effective forecasting model, it provides a reference against which more sophisticated models can be compared. Using the naïve approach produced forecasts are equal to the last observed value. It simply states that future GHI values will be equal to observed GHI at time t (i.e. the atmospheric conditions and solar irradiation remain unchanged between current time t and future time t+h).

-Scaled persistence. Using the clear sky model (simplified solis model [START_REF] Mueller | Rethinking satellite-based solar irradiance modelling: The SOLIS clear-sky module[END_REF]), this training-free prediction model can be defined by:

𝐺𝐻𝐼 ̂𝑆𝑃 (𝑡 + 1) = 𝐺𝐻𝐼(𝑡). 𝐺𝐻𝐼 𝑐𝑙𝑠𝑘 (𝑡+1) 𝐺𝐻𝐼 𝑐𝑙𝑠𝑘 (𝑡) (1)
-Kalman filter. It is a recursive estimator. This means that only the estimated state from the previous time step and the current measurement are needed to compute the estimate for the current state. The Kalman filter can be written as a single equation, however it is most often constructed with two distinct phases: prediction and update. The prediction phase uses the state estimated from the previous timestep (t-1) to produce an estimation of the state at the current timestep (t). This predicted state estimated is also known as the a priori state estimated because, although it is an estimation of the state at the current timestep, it does not include observation information from the current timestep. In the update phase, the current a priori prediction is combined with current observation information to refine the state estimated. This improved estimation is termed the a posteriori state estimated. With this methodology, the forecasting algorithm becomes [START_REF] Grewal | Kalman filtering[END_REF]:

𝐺𝐻𝐼(𝑡 + 1) = 𝐴(𝑡). 𝐺𝐻𝐼(𝑡) + 𝜔(𝑡) (2)
With 𝜔 a multivariate normal distribution with covariance Q (=ℵ(0,Q)) and 𝐴(𝑡) = 𝐺𝐻𝐼 𝑐𝑙𝑠𝑘 (𝑡+1)

𝐺𝐻𝐼 𝑐𝑙𝑠𝑘 (𝑡) . At time t, an observation (or measurement) z(t) of the true state GHI(t) is made according to:

𝑧(𝑡) = 𝐻(𝑡). 𝐺𝐻𝐼(𝑡) + 𝑣(𝑡) (3)
Where v is the observation noise which is assumed to be zero mean Gaussian white noise with covariance R (=ℵ (0,R)) [START_REF] Julier | New extension of the Kalman filter to nonlinear systems[END_REF]. The initial state, and the noise vectors at each step are all assumed to be mutually independent. In our problem, the prediction defining the state vector is defined by:

𝐺𝐻𝐼 ̂(𝑡|𝑡 -1) = 𝐺𝐻𝐼 ̂(𝑡 -1|𝑡 -1). 𝐴(𝑡 -1) (4)
𝐺𝐻𝐼 ̂(𝑡|𝑡) is the predicted value of GHI given the measured value of GHI at time (t). It is in fact, an a posteriori state estimated at time t given observations up to and including time t. Then the a posteriori error covariance matrix P (a measure of the estimated accuracy of the state estimated) is calculated [START_REF] Shumway | An approach to time series smoothing and forecasting using the EM algorithm[END_REF].

𝑃(𝑡|𝑡 -1) = 𝐴(𝑡 -1). 𝑃(𝑡 -1|𝑡 -1). 𝐴(𝑡 -1) 𝑇 + 𝑄

From this last equation, we define the filter gain K which is then computed:

𝐾(𝑡) = 𝑃(𝑡|𝑡 -1). 𝐻(𝑡). (𝐻(𝑡). 𝑃(𝑡|𝑡 -1) + 𝑅) -1 (6)

A correction factor is then introduced and defined by [START_REF] Harvey | Forecasting, structural time series models and the Kalman filter[END_REF]:

𝐺𝐻𝐼 ̂(𝑡|𝑡) = 𝐺𝐻𝐼 ̂(𝑡|𝑡 -1) + 𝐾(𝑡). (𝑧(𝑡) -𝐻(𝑡). 𝐺𝐻𝐼 ̂(𝑡|𝑡))

𝑃(𝑡|𝑡) = 𝑃(𝑡|𝑡 -1) -𝐾(𝑡). 𝐻(𝑡). 𝑃(𝑡|𝑡 -1) (8)

The prediction for the horizon 1 becomes:

𝐺𝐻𝐼 ̂𝑘𝑎𝑙𝑚𝑎𝑛 (𝑡 + 1|𝑡) = 𝐴(𝑡). 𝐺𝐻𝐼 ̂(𝑡|𝑡) =

𝐺𝐻𝐼 𝑐𝑙𝑠𝑘 (𝑡+1) 𝐺𝐻𝐼 𝑐𝑙𝑠𝑘 (𝑡) . 𝐺𝐻𝐼 ̂(𝑡|𝑡) [START_REF] Harvey | Forecasting, structural time series models and the Kalman filter[END_REF] Note that the approach is easily generalizable for other horizons h (𝐺𝐻𝐼 ̂𝑘𝑎𝑙𝑚𝑎𝑛 (𝑡 + ℎ|𝑡)).

Results

This section details the main results we obtained on the previous presented models for time granularity 1 hour and for 10 horizons (t+1 to t+10). In the We can see that the best predictors (in the Annual case) are AR and NN but these are also the more elaborated methods. It is essentially during the second and last quarters that errors are very large, during these seasons it is ANN and AR which give the best results. Because the prediction of the global forecasting is not only interesting 1 hour in advance, we will propose to study the impact of the considered horizon on the predictors. Figure 1 shows the nRMSE related to the horizon. Previously (1 hour horizon case) we've shown that trainless model (scaled persistence and Kalman filter) were similar to NN and AR models. When horizon increases, this effect becomes false, from the 2 hours horizon, the two last one are very better than the first ones. Moreover the NWP is better than NN and AR from the 6 hours horizon and better than the order_1 and Kalman filter from the 3 hours horizon.

Conclusion

In the scope of satisfying electrical operator needs, several kinds of short term prediction methodologies: a naïve model, a linear model, non-linear models and models without training phase have been described. The difference of performance being very small, it is very difficult to compare all the presented predictors. All of them seem interesting to predict the global solar radiation depending on use situation: timestep, horizon, size of the training set, etc. Applying the parsimony concept, a model using the recursive Kalman method of scaled persistence is used performing the prediction from an ad-hoc filtering. The Kalman filter is a recursive estimator, this means that only the estimated state from the previous time step and the current measurement are needed to compute the estimate for the current state. In contrast to batch estimation techniques, no history of observations and/or estimates is required. Kalman filter is applied here for short-term forecasting for both one hour and one minute data sets timestep. The presented approach presents interesting results as it allows to improve quasisystematically the order_1 prediction and Kalman model performances are, for one hour's horizons, competitive with much more complicated models such as ANN which require both consistent historical data sets (at least 200 days) and computational resources (time consuming, Matlab toolboxes, etc.). This method is easily applicable for poor instrumented or isolated sites but it will be tested on several distinct geographical spots.

Figure 1 .

 1 Figure 1. nRMSE evolution in term of time horizons

Table 1 .

 1 Table 1 are shown results of prediction for all the predictors described previously.

	Annual	Quarter 1	Quarter 2	Quarter 3	Quater 4

Results (nRMSE) in the hourly case of the 8 predictors among the 3 types of forecasters (models, trainless models and recursive trainless model), in bold the best results for each column.


Richard E Ewing
email: richard-ewing@tamu.edu

Hong Wang

A summary of numerical methods for time-dependent advection-dominated partial di erential equations

Keywords: Advection-di usion equations, Characteristic methods, Eulerian methods, Numerical simulations

We give a brief summary of numerical methods for time-dependent advection-dominated partial di erential equations (PDEs), including ÿrst-order hyperbolic PDEs and nonstationary advection-di usion PDEs. Mathematical models arising in porous medium uid ow are presented to motivate these equations. It is understood that these PDEs also arise in many other important ÿelds and that the numerical methods reviewed apply to general advection-dominated PDEs. We conduct a brief historical review of classical numerical methods, and a survey of the recent developments on the Eulerian and characteristic methods for time-dependent advection-dominated PDEs. The survey is not comprehensive due to the limitation of its length, and a large portion of the paper covers characteristic or Eulerian-Lagrangian methods.

Mathematical models

We present mathematical models arising in subsurface porous medium uid ow (e.g. subsurface contaminant transport, reservoir simulation) to motivate time-dependent advection-dominated PDEs. These types of PDEs also arise in many other important ÿelds, such as the mathematical modeling of aerodynamics, uid dynamics (e.g. Euler equations, Navier-Stokes equations) [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF][START_REF] Smoller | Shock Waves and Reaction-Di usion Equations, 2nd Edition[END_REF], meteorology [START_REF] Salmon | Lectures on Geophysical Fluid Dynamics[END_REF], and semiconductor devices [START_REF] Markowich | Semiconductor Equations[END_REF]. 1

Miscible ows

A mathematical model used for describing fully saturated uid ow processes through porous media is derived by using the mass balance equation for the uid mixture [START_REF] Bear | Hydraulics of Groundwater[END_REF][START_REF] Ewing | The Mathematics of Reservoir Simulation[END_REF]

@ @t ( ) -• K (p -g) = q; x ∈ ; t ∈ [0; T ]: (1.1)
Here is the physical domain, u, p, and are the Darcy velocity, the pressure, and the mass density of the uid, K (x) is the absolute permeability of the medium, is the dynamic viscosity of the uid, g is the acceleration vector due to gravity, and q represents the source and sink terms, which is often modeled via point or line sources and sinks.

The transport of a speciÿc component in the uid mixture is governed by the mass conservation for the component and is expressed as

@( c) @t + • (uc) -• (D(u)c) = cq; x ∈ ; t ∈ [0; T ]: (1.2)
Here c, a fraction between 0 and 1, represents the concentration of the component, is the porosity of the medium, c(x; t) is either the speciÿed concentrations of the injected uids at sources or the resident concentrations at sinks, and D(u) is the di usion-dispersion tensor.

Multiphase ows

When either air or a nonaqueous-phase liquid (NAPL) contaminant is present in groundwater transport processes, this phase is immiscible with the water phase and the two phases ow simultaneously in the ow process. Likewise, in the immiscible displacement in petroleum production, the oil phase and the water phase are immiscible. In both cases, there is no mass transfer between the two phases and so the following equations hold for each phase [START_REF] Bear | Hydraulics of Groundwater[END_REF][START_REF] Chavent | A new formulation of diphasic incompressible ows in porous media[END_REF][START_REF] Chavent | Mathematical Models and Finite Elements for Reservoir Simulation[END_REF][START_REF] Ewing | The Mathematics of Reservoir Simulation[END_REF]: @ @t ( j S j ) -• j K k rj j p j = j q j ; x ∈ ; t ∈ [0; T ]:

Here S j , u j , j , p j , k rj , j , and q j are the saturation, velocity, density, pressure, relative permeability, viscosity, and source and sink terms for the phase j. The indices j=n and w stand for the nonwetting and wetting phases, respectively. The saturations S n and S w satisfy the relation S n + S w = 1.

Eqs. (1.3) may be rearranged in a form that resembles Eqs. (1.1) and (1.2) by letting S n = 1 -S w . The pressure between the two phases is described by the capillary pressure p c (S w ) = p n -p w . The global pressure p and total velocity u of a two-phase ow model is given by the following equations [START_REF] Chen | Comparison of various formulations of three-phase ow in porous media[END_REF]:

S n C n
Dp Dt -• (K p) = q(x; S w ; p); x ∈ ; t ∈ [0; T ]; (1.4) where (D=Dt) = (@=@t) + (u n =S n ) • and p = 1 2 (p n + p w ) + 1 2 S Sc (( nw )= )(dp c =d ) d with p c (S c ) = 0. The total mobility = n + w , the phase mobility j = k rj = j , and the compressibility C j = (1= j )(d j =dp j ) are functions of time and space.

The governing equation for the wetting phase now has a form @S w @t + • (f(S w )u -D(S w )S w ) = q w ; x ∈ ; t ∈ [0; T ]; (1.5) where c = nw . The capillary di usion D(S w ) = -K n f w (dp c =dS w ) and the fractional ow functions f j = j = .

In practice, the di usion term in Eq. (1.2) or (1.5) is often a small phenomenon relative to advection. Hence, these equations are time-dependent advection-di usion partial di erential equations (PDEs) in terms of the concentration c or the saturation S. In particular, Eq. (1.5) has an S-shaped nonlinear ux function f and a degenerate capillary di usion term [START_REF] Chavent | Mathematical Models and Finite Elements for Reservoir Simulation[END_REF][START_REF] Ewing | The Mathematics of Reservoir Simulation[END_REF]. Sometimes the di usion phenomenon is so small that its e ect is neglected. In this case, Eq. (1.2) or (1.5) is reduced to a ÿrst-order hyperbolic PDE. Finally, initial and boundary conditions also need to be speciÿed to close the system (1.1) -(1.2) or (1.4) -(1.5).

Conventional ÿnite di erence and ÿnite element methods

We carry out a brief historical review of classical numerical methods in this section and a survey of the recent developments on the Eulerian and characteristic methods in the next section primarily for time-dependent advection-dominated PDEs, including ÿrst-order hyperbolic PDEs and nonstationary advection-di usion PDEs. Because of the extensive research carried out in these areas, it is impossible to describe adequately all these developments in the space available. Hence, this review is not comprehensive in that we try to describe and review only some representatives of the huge amount of works in the literature. Notice that since relatively more references and survey papers can be found on the Eulerian methods for unsteady state advection-dominated PDEs, we intend to use a relatively large portion to cover characteristic or Eulerian-Lagrangian methods for advection-dominated PDEs. Finally, we refer interested readers to the works of Morton [START_REF] Morton | Numerical Solution of Convection-Di usion Problems[END_REF] and Roos et al. [START_REF] Roos | Numerical Methods for Singular Perturbed Di erential Equations[END_REF] for detailed descriptions on the recent developments for the numerical methods for stationary advection-di usion PDEs.

It is well known that advection-dominated PDEs present serious numerical di culties due to the moving steep fronts present in the solutions of advection-di usion transport PDEs or shock discontinuities in the solutions of pure advection PDEs or advection-di usion PDEs with degenerate diffusion. Additional di culties include the strong couplings and nonlinearities of advection-dominated PDE systems, the e ect of the singularities at point=line sources and sinks, the strong heterogeneity of the coe cients, anisotropic di usion-dispersion in tensor form, and the enormous sizes of ÿeld-scale applications.

Finite di erence methods (FDMs)

Due to their simplicity, FDMs were ÿrst used in solving advection-dominated PDEs. For convenience, of presentation, we consider the one-dimensional constant-coe cient analogue of Eq. (1.2)

@c @t + V @c @x -D @ 2 c @x 2 = 0; x ∈ (a; b); t ∈ [0; T ] (2.1)
and assume a uniform spatial and temporal partition x i = a +i x for i =0; 1; : : : ; I with x =(b -a)=I and t m = m t for m = 0; 1; : : : ; M with t = T=M . We deÿne the Courant number Cr = V t= x and the Peclet number Pe = V x=D. It is known that the solution to the space-centered explicit scheme

c m+1 i -c m i t + V c m i+1 -c m i-1 2 x -D c m i+1 -2c m i + c m i-1 ( x) 2 = 0 (2.2)
does not oscillate only when the Peclet number Pe62 and the CFL condition (Cr61) is satisÿed [START_REF] Courant | Uber die partiellen di erenzen-gleichungen der mathematisches physik[END_REF][START_REF] Price | Applications of oscillation matrices to di usion-convection equations[END_REF]. For Pe ¿ 2, damped oscillations occur with nonreal eigenvalues [START_REF] Finlayson | Numerical Methods for Problems with Moving Fronts[END_REF][START_REF] Price | Applications of oscillation matrices to di usion-convection equations[END_REF]. Furthermore, for the linear hyperbolic PDE @c @t

+ V @c @x = 0; x ∈ (a; b); t ∈ [0; T ]; (2.3)
which can be viewed as a limiting case of D → 0 in Eq. (2.1), the corresponding scheme to scheme (2.2)

c m+1 i -c m i t + V c m i+1 -c m i-1 2 x = 0 (2.4)
is unconditionally unstable [START_REF] Finlayson | Numerical Methods for Problems with Moving Fronts[END_REF][START_REF] Strikwerda | Finite Di erence Schemes and Partial Di erential Equations[END_REF].

The upwind FDM (UFDM) uses a one-sided ÿnite di erence in the upstream direction to approximate the advection term in the transport PDE (2.1) and can be expressed as follows (assuming V ¿0):

c m+1 i -c m i t + V c m i -c m i-1 2 x -D c m i+1 -2c m i + c m i-1 ( x) 2 = 0: (2.5)
The Lax-Friedrichs scheme

c m+1 i -(c m i+1 + c m i-1 )=2 t + V c m i+1 -c m i-1 2 x -D c m i+1 -2c m i + c m i-1 ( x) 2 = 0 (2.6)
is obtained by replacing c m i in the ÿrst term in Eq. (2.2) by its mean value (c m i+1 + c m i-1 )=2.

Remark 1. Schemes (2.5) and (2.6) eliminate the nonphysical oscillations present in Scheme (2.2), and generate stable solutions even for very complicated multiphase and multicomponent ows. It can be shown that the UFDM scheme is actually a second-order approximation to Eq. (2.1) with a modiÿed di usion D(1 + (Pe=2)(1 -Cr)), while the Lax-Friedrichs scheme is a second-order approximation to Eq. ( 2.1) with an extra numerical di usion (( x) [START_REF] Ewing | The Mathematics of Reservoir Simulation[END_REF][START_REF] Hedstrom | Models of di erence schemes for u t + ux = 0 by partial di erential equations[END_REF][START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. Hence, these methods introduce excessive numerical di usion and the numerical solutions are dependent upon grid orientation. Detailed description on the theory and the use of modiÿed equations can be found in [START_REF] Hedstrom | Models of di erence schemes for u t + ux = 0 by partial di erential equations[END_REF][START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF][START_REF] Warming | The modiÿed equation approach to the stability and accuracy analysis of ÿnite di erence methods[END_REF].

2 =2 t)(1-Cr 2 )
The Lax-Wendro scheme is based on the Taylor series expansion and Eq. (2.3)

c(x; t m+1 ) = c(x; t m ) + t @c(x; t m ) @t + ( t) 2 2 @ 2 c(x; t m ) @t 2 + O(( t) 3 ) = c(x; t m ) -V t @c(x; t m ) @x + (V t) 2 2 @ 2 c(x; t m ) @x 2 + O(( t) 3 ): (2.7)
Dropping the O(( t) 3 ) term in Eq. (2.7) and using centered di erences to approximate the spatial derivatives yields the Lax-Wendro scheme

c m+1 i = c m i - Cr 2 (c m i+1 -c m i-1 ) + Cr 2 2 (c m i+1 -2c m i + c m i-1 ); (2.8)
which is a second-order scheme. The Beam-Warming scheme is a one-sided version of the Lax-Wendro scheme. It uses secondorder accurate one-sided di erences to approximate the spatial derivatives in Eq. (2.7)

c m+1 i = c m i - Cr 2 (3c m i -4c m i-1 + c m i-2 ) + Cr 2 2 (c m i -2c m i-1 + c m i-2 ): (2.9)
Remark 2. The Lax-Wendro scheme and the Beam-Warming scheme give third-order approximations to the modiÿed advection-dispersion equation

@c @t + V @c @x -ÿ @ 3 c @x 3 = 0; x ∈ (a; b); t ∈ [0; T ] with ÿ = (V ( x) 2 =6)(Cr 2 -1) for (2.8
) and (V ( x) 2 =6)(2 -3Cr + Cr 2 ) for (2.9). The theory of dispersive waves and its utility in the study of numerical methods are covered in [START_REF] Trefethen | Group velocity in ÿnite di erence schemes[END_REF][START_REF] Whitham | Linear and Nonlinear Waves[END_REF], which show that the Lax-Wendro scheme tends to develop oscillations behind shock fronts while the Beam-Warming scheme tends to develop oscillations in front of shock fronts.

Remark 3. Solving Eq. (2.3) yields c(x; t m+1 ) = c(x -V t; t m ). When the CFL condition is satisÿed, the UPFD or the Lax-Friedrichs scheme can be viewed as an linear interpolation of c(x -V t; t m ) by the nodal values c(x i-1 ; t m ) and c(x i ; t m ), or c(x i-1 ; t m ) and c(x i+1 ; t m ), respectively. This explains why these schemes are free of oscillations and introduce smearing from another point of view [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. Second, because UFDM takes advantage of upstream information, it is slightly more accurate than the Lax-Friedrichs scheme. On the other hand, the latter is symmetric and can be easily implemented, which is an important feature for nonlinear hyperbolic conservation laws. In contrast, the Lax-Wendro scheme (2.8) or Beam-Warming scheme (2.9) can be viewed as a quadratic interpolation of c(x -V t; t m ) by the nodal values c(x i-1 ; t m ), c(x i ; t m ), and c(x i+1 ; t m ), or c(x i ; t m ), c(x i+1 ; t m ), and c(x i+2 ; t m ). This is why they introduce oscillations across shock discontinuities.

The leap-frog scheme for Eq. (2.3) is obtained by replacing the forward di erence in time in (2.4) by a centered di erence

c m+1 i -c m-1 i t + V c m i+1 -c m i-1 2 x = 0: (2.10)
Scheme (2.10) has an improved truncation error of O(( x) 2 + ( t) 2 ), but it is a multi-level scheme. This leads to increased computational storage, a particular disadvantage for large multi-dimensional nonlinear systems. These methods can be extended to solve nonlinear hyperbolic conservation laws

@c @t + @f(c) @x = 0; x ∈ (a; b); t ∈ [0; T ] (2.11)
and their viscous analogue [START_REF] Finlayson | Numerical Methods for Problems with Moving Fronts[END_REF][START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. For example, a large class of upwind schemes have been developed, based on the Godunov scheme [START_REF] Godunov | A di erence scheme for numerical computation of discontinuous solutions of uid dynamics[END_REF]; they have often been presented in terms of Riemann solvers. The Lax-Friedrichs scheme is a basis for the development of nonoscillatory central schemes (see e.g. [START_REF] Nessyahu | Non-oscillatory central di erencing for hyperbolic conservation laws[END_REF]).

Galerkin and Petrov-Galerkin ÿnite element methods (FEMs)

Many FEM schemes have been developed in parallel. For instance, the Galerkin and Petrov-Galerkin FEMs that are analogues to Scheme (2.2) and the UFDM (2.5) for Eq. (2.1) can be uniformly written as follows:

b a c(x; t m+1 )w i (x) dx - b a c(x; t m )w i (x) dx + t b a D @c(x; t m+1 ) @x @w i (x) @x d x + b a V @c(x; t m+1 ) @x w i (x) dx = -(1 -) t b a D @c(x; t m ) @x @w i (x) @x d x + b a V @c(x; t m ) @x w i (x) dx : (2.12)
Here c(x; t m+1 ) is a piecewise-linear trial function. In the linear Galerkin FEM, the test functions w i (x) are standard hat functions centered at the node x i and correspond to the space-centered scheme (2.2) (see e.g. [START_REF] Ewing | The Mathematics of Reservoir Simulation[END_REF][START_REF] Finlayson | Numerical Methods for Problems with Moving Fronts[END_REF]).

In the quadratic Petrov-Galerkin FEM (QPG), the test functions are constructed by adding an asymmetric perturbation to the original piecewise-linear hat functions [START_REF] Barrett | Approximate symmetrization and Petrov-Galerkin methods for di usion-convection problems[END_REF][START_REF] Celia | A new numerical approach for the advective-di usive transport equation[END_REF][START_REF] Christie | Finite element methods for second order di erential equations with signiÿcant ÿrst derivatives[END_REF]

w i (x) =              x -x i-1 x + (x -x i-1 )(x i -x) x 2 ; x ∈ [x i-1 ; x i ]; x i+1 -x x - (x -x i )(x i+1 -x) x 2 ; x ∈ [x i ; x i+1 ]; 0 otherwise:
With a choice of =3, the QPG reproduces the UFDM. With an optimal choice of =3[coth(Pe=2)-2=Pe], the QPG is reduced to the optimal FDM of Allen and Southwell [START_REF] Allen | Relaxation methods applied to determining the motion in two dimensions of a uid past a ÿxed cylinder[END_REF]. For a stationary analogue of Eq. (2.1), the QPG method yields solutions that coincide with the exact solution at the nodal points, and minimizes the errors in approximating spatial derivatives [START_REF] Barrett | Approximate symmetrization and Petrov-Galerkin methods for di usion-convection problems[END_REF][START_REF] Celia | A new numerical approach for the advective-di usive transport equation[END_REF]. However, the QPG is susceptible to strong time truncation errors that introduce numerical di usion and the restrictions on the size of the Courant number, and hence tends to be ine ective for transient advection-dominated PDEs.

In the cubic Petrov-Galerkin FEM (CPG), the test functions are deÿned as the original piecewiselinear hat functions with a symmetric cubic perturbation added to each nonzero piece [START_REF] Bouloutas | An improved cubic Petrov-Galerkin method for simulation of transient advectiondi usion processes in rectangularly decomposable domains[END_REF][START_REF] Westerink | Consider higher degree Petrov-Galerkin methods for the solution of the transient convection-di usion equation[END_REF]

w i (x) =              x -x i-1 x + (x -x i-1 )(x i -x)(x i-1 + x i -2x) x 3 ; x ∈ [x i-1 ; x i ]; x i+1 -x x - (x -x i )(x i+1 -x)(x i + x i+1 -2x) x 3 ; x ∈ [x i ; x i+1 ]; 0 otherwise:
Here = 5Cr 2 . The CPG intends to use nonzero spatial error to cancel the temporal error to improve the overall accuracy. In these treatments the e ects on mass balance come from spatial dependence of test functions in the ÿrst terms on both the sides of Eq. (2.12). Detailed descriptions of the FDMs and FEMs that have been used in the petroleum industry can be found in [START_REF] Ewing | The Mathematics of Reservoir Simulation[END_REF][START_REF] Russell | Finite element and ÿnite di erence methods for continuous ows in porous media[END_REF].

Corresponding to the Lax-Wendro scheme (2.8) and the leap-frog scheme (2.10) is the Taylor-Galerkin scheme

c m+1 i+1 -c m-1 i+1 6 + 2(c m+1 i -c m-1 i ) 3 + c m+1 i-1 -c m-1 i-1 6 = - Cr 2 (c m i+1 -c m i-1 ) + Cr 2 2 (c m i+1 -2c m i + c m i-1 );
and the leap-frog Galerkin scheme

c m+1 i+1 -c m-1 i+1 6 t + 2(c m+1 i -c m-1 i ) 3 t + c m+1 i-1 -c m-1 i-1 6 t + V c m i+1 -c m i-1 2 x = 0:
In addition, a wide variety of other methods can be devised for advection-dominated transport PDEs by using di erent FDM and FEM approximations, or Taylor expansions. Many large-scale simulators use fully implicit discretization so that large time steps can be allowed. However, in implicit methods, the temporal error and the spatial error add together. Hence, increasing the size of time steps can signiÿcantly reduce the accuracy of the solutions [START_REF] Ewing | The Mathematics of Reservoir Simulation[END_REF]. This is also observed computationally [START_REF] Wang | A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations[END_REF].

In contrast, in explicit schemes the temporal error and the spatial error cancel each other. Hence, reducing the time step size further with ÿxed spatial step size will actually reduce the accuracy of the numerical solutions. The sizes of spatial grids and temporal steps have to be reduced simultaneously to improve the accuracy of the solutions, leading to signiÿcantly increased overall computational and storage cost [START_REF] Wang | A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations[END_REF].

Recent developments for advection-di usion PDEs

Recent developments in e ectively solving advection-di usion PDEs have generally been along one of two approaches: Eulerian or characteristic Lagrangian methods. Eulerian methods use the standard temporal discretization, while the main distinguishing feature of characteristic methods is the use of characteristics to carry out the discretization in time.

Eulerian methods for advection-di usion PDEs

Many methods directly apply to a nonconservative analogue of Eq. (1.2)

@c @t + u • c -• (Dc) = cq; x ∈ ; t ∈ [0; T ]: (3.1)
3.1.1. The streamline di usion ÿnite element method (SDFEM)

The SDFEM directly applies to Eq. (3.1). It is based on the framework of space-time FEMs on the space-time strip × [t m ; t m+1 ], and uses continuous and piecewise polynomial trial and test functions in space as standard FEM but a discontinuous Galerkin approximation in time at time level t m and t m+1 such that

t m+1 t m @c @t + u • c -• (Dc) w + @w @t + u • w dx dt + t m+1 t m w • (Dc) dx dt + c(x; t m + )w(x; t m + ) dx = t m+1 t m cq w + @w @t + u • w dx dt + c(x; t m -)w(x; t m + ) dx:
Here w(x; t m + )=lim t→t m ; t¿t m w(x; t) and w(x; t m -)=lim t→t m ; t¡t m w(x; t). At the initial time step, c(x; t 0 -)= c 0 (x) is the prescribed initial condition. The second term on the left-hand side is carried out elementwise, since it is not well deÿned for continuous and piecewise polynomials. The parameter , which determines the amount of numerical di usion introduced, is typically chosen to be of order O( ( x) 2 + ( t) 2 ).

The SDFEM was ÿrst proposed by Hughes and Brooks [START_REF] Hughes | A multidimensional upwind scheme with no crosswind di usion[END_REF]. Since then, various SDFEM schemes have been developed and studied extensively by Brooks and Hughes and Hughes [START_REF] Brooks | Streamline upwind Petrov-Galerkin formulations for convection dominated ows with particular emphasis on the incompressible Navier-Stokes equations[END_REF][START_REF] Hughes | Multiscale phenomena: Green functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles, and the origins of stabilized methods[END_REF] and Johnson et al. [START_REF] Hansbo | A velocity-pressure streamline di usion ÿnite element method for the incompressible Navier-Stokes equations[END_REF][START_REF] Johnson | An analysis of discontinuous Galerkin methods for a scalar hyperbolic equation[END_REF]. The SDFEM adds a numerical di usion only in the direction of streamlines to suppress the oscillation and does not introduce any crosswind di usion. However, the undetermined parameter in the SDFEM scheme needs to be chosen very carefully in order to obtain accurate numerical results. An optimal choice of the parameter is heavily problem-dependent. We refer readers to the work of Shih and Elman on the study of the choice in the SDFEM formulation and the related numerical experiments [START_REF] Shih | Modiÿed streamline di usion schemes for convection-di usion problems[END_REF].

While the SDFEM can capture a jump discontinuity of the exact solution in a thin region, the numerical solution may develop over-and under-shoots about the exact solution within this layer. A modiÿed SDFEM with improved shock-capturing properties was proposed [START_REF] Hughes | A new ÿnite element formulation for computational uid dynamics: III, The general streamline operator for multidimensional advective-di usive systems[END_REF][START_REF] Johnson | On the convergence of shock-capturing streamline di usion ÿnite element methods for hyperbolic conservation laws[END_REF], which consists of adding a "shock-capturing" term to the di usion by introducing a "crosswind" control that is close to the steep fronts or "shocks". This modiÿed SDFEM performs much better in terms of catching the steep fronts or the jump discontinuities of the exact solutions. However, the modiÿed SDFEM is a nonlinear scheme and involves another undetermined parameter.

Total variation diminishing (TVD) methods

Notice that when oscillations arise, the numerical solutions will have larger total variation. TVD methods are designed to yield well-resolved, nonoscillatory shock discontinuities by enforcing that the numerical schemes generate solutions with nonincreasing total variations. One approach is to take a high-order method and add an additional numerical di usion term to it. Since this numerical di usion is needed only near discontinuities, one wants it to vanish su ciently quickly so that the order of accuracy of the method on smooth regions of the solutions is retained. Hence, the numerical di usion should depend on the behavior of the solutions, being larger near shock regions than in smooth regions. This leads to a nonlinear method even for the linear advection equation (2.3). The idea of adding a variable amount of numerical di usion dates back to some of the earliest work on the numerical solution of uid dynamics [START_REF] Courant | On the solution of nonlinear hyperbolic di erential equations by ÿnite di erences[END_REF][START_REF] Lax | Systems of conservation laws[END_REF][START_REF] Von Neumann | A method for the numerical calculation of hydrodynamic shocks[END_REF]. The di culty with this approach is that it is hard to determine an appropriate amount of numerical di usion that introduces just enough dissipation without causing excessive smearing.

For this reason, the high-resolution methods developed more recently are based on fairly di erent approaches, including ux-and slope-limiter approaches that impose the nonoscillatory requirement more directly. In the ux-limiter approach, one ÿrst chooses a high-order numerical ux F H (c; i) = F H (c i-lH ; c i-l+1 ; : : : ; c i+rH ) that generates accurate approximations in smooth regions and a low-order numerical ux F L (c; i) = F L (c i-lL ; c i-l+1 ; : : : ; c i+rL ) that yields nonoscillatory solutions near shock discontinuities. One then combines F H and F L into a single numerical ux F, e.g. in the form of

F(c; i) = F L (c; i) + (c; i)(F H (c; i) -F L (c; i)); (3.2)
such that F reduces to F H in smooth regions and to F L in shock regions. Here (c; i), the ux limiter, should be near one in smooth regions and close to zero near shock discontinuities.

The ux-corrected transport (FCT) method of Boris and Book can be viewed as one of the earliest ux limiter methods [START_REF] Boris | Flux-corrected transport I, SHASTA, a uid transport algorithm that works[END_REF][START_REF] Boris | Flux-corrected transport, III, Minimal-error FCT algorithms[END_REF][START_REF] Zalesak | Fully multidimensional ux-corrected transport algorithms for uids[END_REF]. In the FCT method, an anti-di usive term (i.e., the correction term in (3.2)) is added to reduce the excessive numerical di usion introduced by the lower-order ux F L as much as possible without increasing the total variation of the solution.

Sweby studied a family of ux-limiter methods in [START_REF] Sweby | High resolution schemes using ux limiters for hyperbolic conservation laws[END_REF]. By choosing

F L (c m ; i) = Vc m i and F H (c m ; i) = Vc m i + 1 2 V (1 -Cr)(c m i+1 -c m i )
; to be the ÿrst-order upwind ux in (2.5) and the Lax-Wendro ux in (2.8) and using (3.2), a family of ux-limiter methods can be deÿned

c m+1 i = c m i - t x [F(c m ; i) -F(c m ; i -1)] (3.3)
with the ux F(c m ; i) being given by

F(c m ; i) = Vc m i + (c m ; i) 2 V (1 -Cr)(c m i+1 -c m i ):
One way to measure the smoothness of the solution is to look at the ratio of consecutive gradients and to deÿne the ux limiter accordingly

(c m ; i) = ( i ) with  i = c m i -c m i-1 c m i+1 -c m i :
Sweby obtained algebraic conditions on the limiter functions that guarantee second-order accuracy and the TVD property of the derived methods [START_REF] Sweby | High resolution schemes using ux limiters for hyperbolic conservation laws[END_REF]. Harten proved a su cient condition on that can be used to impose constraints on [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF][START_REF] Harten | On ÿnite-di erence approximations and entropy conditions for shocks[END_REF]. Among the di erent choices of limiters are the "superbee" limiter of Roe [START_REF] Roe | Some contributions to the modeling of discontinuous ows[END_REF] (Â) = max{0; min{1; 2Â}; min{Â; 2}} and a smoother limiter by van Leer [START_REF] Van Leer | Towards the ultimate conservative di erence scheme II, Monotonicity and conservation combined in a second order scheme[END_REF] (Â) = |Â| + Â 1 + |Â| :

The extension of ux limiter methods to nonlinear conservation laws and numerical comparisons can be found in [START_REF] Colella | The piecewise-parabolic method (PPM) for gasdynamical simulations[END_REF][START_REF] Sweby | High resolution schemes using ux limiters for hyperbolic conservation laws[END_REF][START_REF] Zalesak | A preliminary comparison of modern shock-capturing schemes: linear advection equation[END_REF].

Another approach is to use slope limiters. These intend to replace the piecewise-constant representation of the solutions in Godunov's method by more accurate representations, and can be expressed in the following steps for Eq. (2.11):

(i) Given the piecewise-constant cell-average representation { c m i } i=+∞ i=-∞ of the solution at time level t m , deÿne a (e.g., piecewise-linear) reconstruction at time t m by ĉ(x;

t m ) = c m i + m i (x -x i ) (3.4) for x on the cell [x i-1=2 ; x i+1=2 ].
Here m i is a slope on the ith cell that is based on the data { c m i }. (ii) Solve Eq. (2.11) with the data ĉ(x; t m ) at time t m to obtain the solution c(x; t m+1 ) at time t m+1 . (iii) Compute the cell average { c m+1 i } i=+∞ i=-∞ of the solution c(x; t m+1 ) at time t m+1 . Note that the cell average of reconstruction (3.4) is equal to c m i on the cell [x i-1=2 ; x i+1=2 ] for any choice of m i . Since Steps 2 and 3 are also conservative, the methods with slope limiters are conservative. Secondly, the choice of m i = 0 in Eq. (3.4) recovers Godunov's method. It is well known that Godunov's method generates solutions with excessive numerical di usion. More accurate reconstructions, such as Eq. (3.4), could be used to reduce the numerical di usion and to improve the accuracy of the numerical solutions.

In the context of the linear advection PDE (2.3), the solution of Step (ii) is simply c(x; t m+1 ) = ĉ(x-V t; t m ). Computing the cell average of c(x; t m+1 ) in Step (iii) leads to the following expression:

c m+1 i = c m i -Cr( c m i -c m i-1 ) - x 2 Cr(1 -Cr)( m i -m i-1 ): (3.5) A natural choice of m i = ( c m i+1 -c m i )=
x in Eq. (3.5) leads to the Lax-Wendro method. Thus, it is possible to obtain second-order accuracy by this approach. Secondly, the slope-limiter methods could generate oscillatory solutions (since the Lax-Wendro method could do so), if the slope limiters m i are not chosen properly. Geometrically, the oscillations are due to a poor choice of slopes, which leads to a piecewise-linear reconstruction ĉ(x; t m ) with much larger total variation than the given data { c m i } i=+∞ i=-∞ [START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. Hence, because of their importance, extensive research has been conducted on how to choose the slope m i in Eq. (3.4) to ensure the resulting methods to be total variation diminishing (TVD). These methods include the monotonic upstream-centered scheme for conservation laws (MUSCL) of van Leer and Minmod methods among others [START_REF] Colella | A direct Eulerian MUSCL scheme for gas dynamics[END_REF][START_REF] Goodman | A geometric approach to high-resolution TVD schemes[END_REF][START_REF] Van Leer | Towards the ultimate conservative di erence scheme II, Monotonicity and conservation combined in a second order scheme[END_REF][START_REF] Van Leer ; Engquist-Osher | On the relation between the upwind-di erencing schemes of Godunov[END_REF]. The simplest choice of the slope is probably the minmod slope deÿned by

m i = minmod c m i+1 -c m i x ; c m i -c m i-1
x with minmod(a; b) = 1 2 (sgn(a) + sgn(b))min(|a|; |b|).

In concluding this part, we notice the connection between the ux-and slope-limiter methods. Using formulation (3.3), we see that the numerical ux for the slope-limiter method (3.5) is

F(c m ; i) = Vc m i + x 2 V (1 -Cr) m i ;
which is of the same form as the ux-limiter method (3.3) if the slope-limiter m i is related to the ux-limiter (c m ; i) by m i = [(c m i+1 -c m i )= x] (c m ; i).

Essentially nonoscillatory (ENO) schemes and weighted essentially nonoscillatory (WENO) schemes

Traditional ÿnite di erence methods are based on ÿxed stencil interpolations of discrete data using polynomials. The resulting scheme is linear for linear PDEs. However, ÿxed stencil interpolation of second-or higher-order accuracy is necessarily oscillatory across a discontinuity; this is why the Lax-Wendro scheme (2.8) and the Beam-Warming scheme (2.9) introduce oscillations across shock discontinuities (see Remark 3). One common approach to eliminate or reduce spurious oscillations near discontinuities is to add a numerical di usion as in the SDFEM presented earlier. The numerical di usion should be tuned so that it is large enough near discontinuities but is small enough elsewhere to maintain high-order accuracy. One disadvantage of this approach is that it is hard to determine an appropriate amount of numerical di usion that introduces just enough dissipation without causing excessive smearing. Another approach is to apply ( ux or slope) limiters to eliminate the oscillations. By carefully designing such limiters (e.g., reducing the slope of a linear interpolant or using a linear rather than a quadratic interpolant near shock discontinuities), the TVD property could be achieved for some numerical schemes for nonlinear scalar conservation laws in one space dimension. Unfortunately, Osher and Chakravarthy proved that TVD methods must degenerate to ÿrst-order accuracy at local maximum or minimum points [START_REF] Osherm | High resolution schemes and the entropy condition[END_REF].

The ENO and WENO schemes are high-order accurate ÿnite di erence=volume schemes designed for nonlinear hyperbolic conservation laws with piecewise smooth solutions containing discontinuities [START_REF] Harten | Uniformly high-order accurate essentially nonoscillatory schemes, III[END_REF][START_REF] Harten | Uniformly high-order accurate nonoscillatory schemes, I[END_REF][START_REF] Jiang | E cient implementation of weighted ENO schemes[END_REF][START_REF] Liu | Weighted essentially nonoscillatory schemes[END_REF]. By delicately deÿning a nonlinear adaptive procedure to automatically choose the locally smooth stencil, the ENO and WENO schemes avoid crossing discontinuities in the interpolation procedure and thus generate uniformly high-order accurate, yet essentially nonoscillatory solutions. These schemes have been quite successful in applications, especially for problems containing both shock discontinuities and complicated smooth solution structures [START_REF] Shu | Essentially nonoscillatory (ENO) and weighted essentially nonoscillatory (WENO) schemes for hyperbolic conservation laws[END_REF].

The discontinuous Galerkin (DG) method

The original discontinuous Galerkin ÿnite element method was introduced by Reed and Hill for solving a linear neutron transport equation [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF], in which the method can be carried out element by element when the elements are suitably ordered according to the characteristic directions. Lesaint and Raviart [START_REF] Lesaint | On a ÿnite element method for solving the neutron transport equation[END_REF] carried out the ÿrst analysis for this method and proved a convergence rate of ( x) k for general triangular partitions and ( x) k+1 for Cartesian grids. Johnson and Pitkar anta [START_REF] Johnson | An analysis of discontinuous Galerkin methods for a scalar hyperbolic equation[END_REF] obtained an improved estimate of ( x) k+1=2 for general triangulations, which is conÿrmed to be optimal by Peterson [START_REF] Peterson | A note on the convergence of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF]. Chavent and Salzano [START_REF] Chavent | A ÿnite element method for the 1D water ooding problem with gravity[END_REF] constructed an explicit DG method for Eq. (2.11), in which piecewise linear FEM is used in space and an explicit Euler approximation is used in time. Unfortunately, the scheme is stable only if the Courant number Cr =O( √ x). Chavent and Cockburn [START_REF] Chavent | The local projection P 0 ; P 1 -discontinuous Galerkin ÿnite element method for scalar conservation laws[END_REF] modiÿed the scheme by introducing a slope limiter, and proved the scheme to be total variation bounded (TVB) when Cr 6 1 2 . However, the slope limiter introduced compromises the accuracy of the approximation in smooth regions. Cockburn and Shu [START_REF] Cockburn | The Runge-Kutta local projection P 1 -discontinuous Galerkin ÿnite element method for scalar conservation laws[END_REF] introduced the ÿrst Runge-Kutta DG (RKDG) method, which uses an explicit TVD second-order Runge-Kutta discretization and modiÿes the slope limiter to maintain the formal accuracy of the scheme at the extrema. The same authors then extended this approach to construct higher-order RKDG methods [START_REF] Cockburn | The Runge-Kutta local projection P 1 -discontinuous Galerkin ÿnite element method for scalar conservation laws II: general framework[END_REF], to multidimensional scalar conservation laws [START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin ÿnite element method for conservation laws IV: the multidimensional case[END_REF][START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin ÿnite element method for conservation laws V: multidimensional systems[END_REF], and to multidimensional systems. We refer interested readers to the survey article [START_REF] Cockburn | Finite element methods for conservation laws, this volume[END_REF] in this volume for detailed discussions on the DG methods.

Characteristic methods

Because of the hyperbolic nature of advective transport, characteristic methods have been investigated extensively for the solution of advection-di usion PDEs. In a characteristic (or Lagrangian) method, the transport of the uid is referred to a Lagrangian coordinate system that moves with the uid velocity. One tracks the movement of a uid particle and the coordinate system follows the movement of the uid. The time derivative along the characteristics of the advection-di usion PDE (3.1) is expressed as

Dc Dt = @c @t + u • c: (3.6)
Consequently, the advection-di usion PDE (3.1) is rewritten as the following parabolic di usionreaction PDE in a Lagrangian system:

Dc Dt -• (Dc) = cq (3.7)
and the advection has seemingly disappeared. In other words, in a Lagrangian coordinate system (that moves with the ow) one would only see the e ect of the di usion, reaction, and the the right-hand side source terms but not the e ect of the advection or moving steep fronts. Hence, the solutions of the advection-di usion PDEs are much smoother along the characteristics than they are in the time direction. This explains why characteristic methods usually allow large time steps to be used in a numerical simulation while still maintaining its stability and accuracy. Unfortunately, Eq. (3.7) is written in a Lagrangian coordinate system, which is constantly moving in time. Consequently, the corresponding characteristic or Lagrangian methods often raise extra and nontrivial analytical, numerical, and implementational di culties, which require very careful treatment. In contrast, Eq.

(1.2) or (3.1) is written in an Eulerian system which is ÿxed in space. Hence, Eulerian methods are relatively easy to formulate and to implement.

Classical characteristic or Eulerian-Lagrangian methods

The classical Eulerian-Lagrangian method is a ÿnite di erence method based on the forward tracking of particles in cells. In this method, the spatial domain is divided into a collection of elements or cells and a number of particles are placed within each cell. Then the governing PDE is used to determine the movement of the particles from cell to cell. In this algorithm, the solution is determined by the number of particles within a cell at any given time. Related works can be found in [START_REF] Farmer | A moving point method for arbitrary Peclet number multi-dimensional convection-di usion equations[END_REF][START_REF] Garder | Numerical calculations of multidimensional miscible displacement by the method of characteristics[END_REF][START_REF] Thomaidis | An explicit ÿnite di erence scheme based on the modiÿed method of characteristics for solving di usion-convection problems in one space dimension[END_REF]. In these methods, the di usion occurs at the time step t m and the solution is advected forward in time to the time step t m+1 , leading to the following scheme for Eq. (2.1):

c m+1 i -c m i t -D c m i+1 -2c m i + c m i-1 ( x) 2 = 0:
Here c m+1 i = c( xi ; t m+1 ) with xi = x i + V t. Because the advected nodes xi need not be nodes at time t m+1 , they are irregular, in general.

Neuman developed an Eulerian-Lagrangian ÿnite element method using a combination of forward and backward tracking algorithms [START_REF] Neuman | An Eulerian-Lagrangian numerical scheme for the dispersion-convection equation using conjugate space-time grids[END_REF][START_REF] Neuman | Adaptive Eulerian-Lagrangian ÿnite element method for advection-dispersion equation[END_REF]. Near a steep front, a forward tracking algorithm is used to move a cloud of particles from time t m to new positions at time t m+1 according to the advection, as done by Garder et al. [START_REF] Garder | Numerical calculations of multidimensional miscible displacement by the method of characteristics[END_REF]. An implicit scheme is then used to treat the di usion at time t m+1 . Away from a front, a backward tracking algorithm is used, in which one ÿnds a point that ends up at position x at time t m+1 .

Eulerian methods carry out the temporal discretization in the time direction, so they cannot accurately simulate all of the wave interactions that take place if the information propagates more than one cell per time step (i.e., if the CFL condition is violated), either for the reason of stability (for explicit methods) or for the reason of accuracy (for implicit methods). By using characteristic tracking, characteristic methods follow the movement of information or particles as well as their interactions. However, forward tracked characteristic methods often distort the evolving grids severely and greatly complicate the solution procedures, especially for multi-dimensional problems.

The modiÿed method of characteristics (MMOC)

In this part we brie y review the MMOC, which was proposed by Douglas and Russell for solving advection-di usion PDEs in a nonconservative form [START_REF] Douglas | Numerical methods for convection-dominated di usion problems based on combining the method of characteristics with ÿnite element or ÿnite di erence procedures[END_REF] and can be viewed as a representative of the Eulerian-Lagrangian methods developed during the same time period [START_REF] Benque | Quelques di culties des modeles numeriques en hydraulique[END_REF][START_REF] Pironneau | On the transport-di usion algorithm and its application to the Navier-Stokes equations[END_REF][START_REF] Pudykiewicz | Some properties and comparative performance of the semi-Lagrangian method of robert in the solution of the advection-di usion equation[END_REF]. Using the Lagrangian form (3.7), we can combine the ÿrst two terms on the left-hand side of (3.1) to form one term through characteristic tracking (3.6) (see, e.g. [START_REF] Douglas | Numerical methods for convection-dominated di usion problems based on combining the method of characteristics with ÿnite element or ÿnite di erence procedures[END_REF])

Dc(x; t m+1 ) Dt ≈ (x) c(x; t m+1 ) -c(x * ; t m ) t (3.8)
with x * = x -u(x; t m+1 ) t= (x). Substituting (3.8) for the ÿrst two terms on the left-hand side of Eq. (3.1) and integrating the resulting equation against any ÿnite element test functions w(x), one obtains the following MMOC scheme [START_REF] Douglas | Numerical methods for convection-dominated di usion problems based on combining the method of characteristics with ÿnite element or ÿnite di erence procedures[END_REF][START_REF] Ewing | Simulation of miscible displacement using mixed methods and a modiÿed method of characteristics[END_REF] for Eq. (3.1):

(x) c(x; t m+1 ) -c(x * ; t m ) t w(x) dx + w(x) • Dc(x; t m+1 ) dx = cq(x; t m+1 )w(x) dx: (3.9)
Eq. (3.9) follows the ow by tracking the characteristics backward from a point x in a ÿxed grid at the time step t m+1 to a point x * at time t m . Hence, the MMOC avoids the grid distortion problems present in forward tracking methods. Moreover, MMOC symmetrizes and stabilizes the transport PDEs, greatly reducing temporal errors; therefore MMOC allows for large time steps in a simulation without loss of accuracy and eliminates the excessive numerical dispersion and grid orientation e ects present in many Eulerian methods [START_REF] Douglas | Simulation of the transient behavior of a one-dimensional semiconductor device[END_REF][START_REF] Ewing | The Mathematics of Reservoir Simulation[END_REF][START_REF] Russell | Finite element and ÿnite di erence methods for continuous ows in porous media[END_REF]. However, the MMOC and the characteristic methods presented earlier have the following drawbacks:

Remark 4. In the context of the MMOC and other characteristic methods using a backtracking algorithm, the (x)c(x * ; t m )w(x) dx term in Eq. (3.9) is deÿned on the domain at time t m+1 . In this term, the test functions w(x) are standard FEM basis functions on at time t m+1 , but the value of c(x * ; t m ) has to be evaluated by a backtracking method where x * = r(t m ; x; t m+1 ) is the point at the foot corresponding to x at the head [START_REF] Douglas | Numerical methods for convection-dominated di usion problems based on combining the method of characteristics with ÿnite element or ÿnite di erence procedures[END_REF][START_REF] Ewing | Simulation of miscible displacement using mixed methods and a modiÿed method of characteristics[END_REF]. For multidimensional problems, the evaluation of this term with a backtracking algorithm requires signiÿcant e ort, due to the need to deÿne the geometry at time t m that requires mapping of points along the boundary of the element and subsequent interpolation and mapping onto the ÿxed spatial grid at the previous time t m [START_REF] Binning | A ÿnite volume Eulerian-Lagrangian localized adjoint method for solution of the contaminant transport equations in two-dimensional multi-phase ow systems[END_REF][START_REF] Morton | Stability of the Lagrangian-Galerkin method with nonexact integration[END_REF]. This procedure introduces a mass balance error and leads to schemes that fail to conserve mass [START_REF] Celia | A new numerical approach for the advective-di usive transport equation[END_REF][START_REF] Morton | Stability of the Lagrangian-Galerkin method with nonexact integration[END_REF][START_REF] Wang | Eulerian-Lagrangian localized methods for convection-di usion equations and their convergence analysis[END_REF]. Moreover, in these methods it is not clear how to treat ux boundary conditions in a mass-conservative manner without compromising the accuracy, when the characteristics track to the boundary of the domain [START_REF] Celia | A new numerical approach for the advective-di usive transport equation[END_REF][START_REF] Russell | Eulerian-Lagrangian localized adjoint methods with variable coe cients in multiple dimensions[END_REF][START_REF] Wang | An ELLAM Scheme for advectiondi usion equations in two dimensions[END_REF][START_REF] Wang | A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations[END_REF][START_REF] Wang | Eulerian-Lagrangian localized methods for convection-di usion equations and their convergence analysis[END_REF].

The modiÿed method of characteristics with adjusted advection (MMOCAA)

Recently, Douglas et al. proposed an MMOCAA scheme to correct the mass error of the MMOC by perturbing the foot of the characteristics slightly [START_REF] Douglas | On the numerical simulation of water ooding of heterogeneous petroleum reservoirs[END_REF][START_REF] Douglas | The modiÿed method of characteristics with adjusted advection[END_REF]. For Eq. (3.1) with a no-ow or periodic boundary condition, the summation of Eq. (3.9) for all the test functions (that add exactly to one) yields the following equation: (x)c(x; t m+1 ) dx -(x)c(x * ; t m ) dx = t cq(x; t m+1 ) dx:

Recall that the term on the right-hand side of this equation is obtained by an Euler approximation to the temporal integral in this term. On the other hand, integrating the original PDE (1.2) on the domain × [t m ; t m+1 ] yields the following equation:

(x)c(x; t m+1 ) dx - (x)c(x; t m ) dx = t m+1
t m cq dx dt:

Therefore, to maintain mass balance, we must have

(x)c(x; t m ) dx ≡ Q m = Q m * ≡ (x)c(x * ; t m ) dx:
For some ÿxed constant Ä ¿ 0, we deÿne

x * + = x - u(x; t m+1 ) (x) t + Ä u(x; t m+1 ) (x) ( t) 2 ; x * -= x - u(x; t m+1 ) (x) t -Ä u(x; t m+1 ) (x) ( t) 2 :
We also deÿne

c # (x * ; t m ) = max{c(x * + ; t m ); c(x * -; t m )}; if Q m * 6Q m ; min{c(x * + ; t m ); c(x * -; t m )}; if Q m * ¿ Q m : Because c(x; t m+1
) is unknown in the evaluation of Q m , an extrapolation of 2c(x; t m ) -c(x; t m-1 ) is used. We set

Q m # = (x)c # (x * ; t m ) dx: If Q m # = Q m * , we let Ä c(x * ; t m ) = c(x * ; t m
). In this case, the mass is not conserved. Otherwise, ÿnd

 m such that Q m =  m Q m * + (1 - m )Q m # and let Ä c(x * ; t m ) =  m c(x * ; t m ) + (1 -Â)c # (x * ; t m ).
In latter case, one has

(x) Ä c(x * ; t m ) dx = Q m :
Hence, mass is conserved globally. In the MMOCAA procedure one replaces c(x * ; t m ) in (3.8) and (3.9) by Ä c(x * ; t m ).

The Eulerian-Lagrangian localized adjoint method (ELLAM)

The ELLAM formalism was introduced by Celia et al. for the solution of one-dimensional advection-di usion PDEs [START_REF] Celia | An Eulerian-Lagrangian localized adjoint method for the advection-di usion equation[END_REF][START_REF] Herrera | Eulerian-Lagrangian localized adjoint methods: the theoretical framework[END_REF]. It provides a general characteristic solution procedure for advectiondominated PDEs, and it presents a consistent framework for treating general boundary conditions and maintaining mass conservation. The ELLAM formulation directly applies to Eq. (1.2) in a conservative form. Multiplying Eq. (1.2) with space-time test functions w that vanish outside × (t m ; t m+1 ] and are discontinuous in time at time t m , and integrating the resulting equation over the space-time domain ×(t m ; t m+1 ], we obtain a space-time weak formulation for Eq. (1.2) with a no ow boundary condition (x)c(x; t m+1 )w(x; t m+1 ) dx +

t m+1 t m w • (Dc) dx dt - t m+1 t m c ( w t + u • w) dx dt = (x)c(x; t m )w(x; t m + ) dx + t m+1
t m cqw dx dt; (3.10) where w(x; t m + ) = lim t→t m + w(x; t) takes into account that w(x; t) is discontinuous in time at time t m . Motivated by the localized adjoint method, the ELLAM formalism chooses the test functions from the solution space of the homogeneous adjoint equation of Eq. (1.2) (e.g. see [START_REF] Celia | An Eulerian-Lagrangian localized adjoint method for the advection-di usion equation[END_REF][START_REF] Herrera | Eulerian-Lagrangian localized adjoint methods: the theoretical framework[END_REF])

-(x) @w @t -u • w -• (Dw) = 0: (3.11)
Because the solution space for Eq. (3.11) is inÿnite dimensional and only a ÿnite number of test functions should be used, an operator splitting technique is applied to Eq. (3.11) to deÿne the test functions.

(i) In the ÿrst splitting, the two terms involving spatial derivatives are grouped together, leading to the following system of equations:

-(x) @w @t = 0; -u • w -• (Dw) = 0:
This splitting leads to a class of optimal test function methods involving upstream weighting in space [START_REF] Barrett | Approximate symmetrization and Petrov-Galerkin methods for di usion-convection problems[END_REF][START_REF] Celia | A new numerical approach for the advective-di usive transport equation[END_REF][START_REF] Christie | Finite element methods for second order di erential equations with signiÿcant ÿrst derivatives[END_REF], which yield solutions with signiÿcant temporal errors and numerical di usion. (ii) In the second splitting, the terms involving ÿrst-order derivatives are grouped together, leading to the following system of equations:

-(x) @w @t -u • w = 0; -• (Dw) = 0:
(3.12)

The ÿrst equation in (3.12) implies that the test functions should be constant along the characteristics deÿned by dr dt = u(r; Â) (r) ;

which re ects the hyperbolic nature of Eq. (1.2) and assures Lagrangian treatment of advection.

The second equation in (3.13) is an elliptic PDE, so standard FEM approximations would be a natural choice for the spatial conÿguration of the test functions.

Using splitting (3.12), we deÿne the test functions to be standard FEM basis functions on the spatial domain at time t m+1 and extend them by a constant into the space-time strip × [t m ; t m+1 ] along the characteristics deÿned by (3.13). Incorporating these test functions into the reference equation (3.10), we obtain an ELLAM scheme as follows:

(x)c(x; t m+1 ) dx + t (w [START_REF] Russell | Eulerian-Lagrangian localized adjoint methods for advection-dominated problems[END_REF][START_REF] Wang | An ELLAM Scheme for advectiondi usion equations in two dimensions[END_REF][START_REF] Wang | A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations[END_REF]. Second, it is proved that the ELLAM scheme conserves mass [START_REF] Celia | An Eulerian-Lagrangian localized adjoint method for the advection-di usion equation[END_REF][START_REF] Russell | Eulerian-Lagrangian localized adjoint methods with variable coe cients in multiple dimensions[END_REF]. Third, in contrast to the MMOC and many other characteristic methods that treat general boundary conditions in an ad hoc manner, the ELLAM formulation can treat any combinations of boundary conditions and provides a systematic way to calculate the boundary conditions accurately [START_REF] Celia | Eulerian-Lagrangian localized adjoint methods for contaminant transport simulations[END_REF][START_REF] Celia | An Eulerian-Lagrangian localized adjoint method for the advection-di usion equation[END_REF][START_REF] Wang | An ELLAM Scheme for advectiondi usion equations in two dimensions[END_REF][START_REF] Wang | A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations[END_REF]. Thus, the ELLAM formulation overcomes the drawbacks of many previous characteristic methods while maintaining their numerical advantages.

Remark 6. Most integrals in the ELLAM scheme (3.14) are standard in FEMs and can be evaluated in a straightforward manner. The only exception is the (x)c(x; t m )w(x; t m + ) dx term on the right-hand side of Eq. (3.14). This term corresponds to the (x)c(x * ; t m )w(x) dx term in the MMOC scheme (3.9). As discussed in Remark 4, in the MMOC and other characteristic methods using a backtracking algorithm the evaluation of the (x)c(x * ; t m )w(x) dx term requires signiÿcant e ort and introduces mass balance error [START_REF] Binning | A ÿnite volume Eulerian-Lagrangian localized adjoint method for solution of the contaminant transport equations in two-dimensional multi-phase ow systems[END_REF][START_REF] Morton | Stability of the Lagrangian-Galerkin method with nonexact integration[END_REF]. In the ELLAM scheme (3.14), the term (x)c(x; t m )w(x; t m + ) dx is evaluated by a forward tracking algorithm that was proposed by Russell and Trujillo [START_REF] Russell | Eulerian-Lagrangian localized adjoint methods with variable coe cients in multiple dimensions[END_REF]. In this approach, an integration quadrature would be enforced on at time t m with respect to a ÿxed spatial grid on which c(x; t m ) is deÿned. The di cult evaluation of w(x; t m + ) = lim t→t m + w(x; t m ) = w( x; t m+1 ) is carried out by a forward tracking algorithm from x at time t m to x = r(t m+1 ; x; t m ) at time t m+1 . Because this forward tracking is used only in the evaluation of the right-hand side of (3.14), it has no e ect on the solution grid or the data structure of the discrete system. Therefore, the forward tracking algorithm used here does not su er from the complication of distorted grids, which complicates many classical forward tracking algorithms.

In the past few years, Wang et al. developed ELLAM schemes for multidimensional advectiondi usion PDEs [START_REF] Wang | Eulerian-Lagrangian localized adjoint methods: analyses, implementations, and applications[END_REF][START_REF] Wang | An ELLAM Scheme for advectiondi usion equations in two dimensions[END_REF][START_REF] Wang | An ELLAM scheme for advection-di usion equations in multi-dimensions[END_REF]; Ewing and Wang [START_REF] Ewing | Eulerian-Lagrangian localized adjoint methods for linear advection or advection-reaction equations and their convergence analysis[END_REF] and Wang et al. [START_REF] Wang | A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations[END_REF] also developed ELLAM schemes for multidimensional advection-reaction PDEs; Celia and Ferrand [START_REF] Celia | A comparison of ELLAM formulations for simulation of reactive transport in groundwater[END_REF] and Healy and Russell [START_REF] Healy | A ÿnite-volume Eulerian-Lagrangian localized adjoint method for solution of the advection-dispersion equation[END_REF][START_REF] Healy | Solution of the advection-dispersion equation in two dimensions by a ÿnite-volume Eulerian-Lagrangian localized adjoint method[END_REF] developed ELLAM schemes in a ÿnite-volume setting. Dahle et al. developed ELLAM for two-phase ow [START_REF] Dahle | Eulerian-Lagrangian localized adjoint methods for a nonlinear convection-di usion equation[END_REF][START_REF] Ewing | Simulation of multiphase ows in porous media[END_REF]. The computational experiments carried out in [START_REF] Wang | An ELLAM Scheme for advectiondi usion equations in two dimensions[END_REF][START_REF] Wang | A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations[END_REF] showed that the ELLAM schemes often outperform many widely used and well-received numerical methods in the context of linear advection-di usion or advection-reaction PDEs. In addition, Binning and Celia developed a backtracking ÿnite-volume ELLAM scheme for unsaturated ow [START_REF] Binning | A ÿnite volume Eulerian-Lagrangian localized adjoint method for solution of the contaminant transport equations in two-dimensional multi-phase ow systems[END_REF], Wang et al. developed an ELLAM-MFEM solution technique for porous medium ows with point sources and sinks [START_REF] Wang | An ELLAM-MFEM solution technique for compressible uid ows in porous media with point sources and sinks[END_REF]. These works illustrate the strength of the ELLAM schemes in solving the coupled systems of advection-di usion PDEs. From a viewpoint of analysis, ELLAM methods introduce further di culties and complexities to the already complicated analyses of characteristic methods. We refer readers to the works of Wang et al. for the convergence analysis and optimal-order error estimates for the ELLAM schemes for advection-di usion or advection-di usion-reaction PDEs [START_REF] Wang | A family of ELLAM schemes for advection-di usion-reaction equations and their convergence analyses[END_REF][START_REF] Wang | An optimal-order error estimate for an ELLAM scheme for two-dimensional linear advection-di usion equations[END_REF][START_REF] Wang | Optimal-order convergence rates for Eulerian-Lagrangian localized adjoint methods for reactive transport and contamination in groundwater[END_REF][START_REF] Wang | Eulerian-Lagrangian localized methods for convection-di usion equations and their convergence analysis[END_REF], and the corresponding analysis of Ewing and Wang for the ELLAM schemes for advection-reaction PDEs [44 -46].

The characteristic mixed ÿnite element method (CMFEM)

The CMFEM was presented by Arbogast et al. in [START_REF] Arbogast | A characteristic-mixed method for contaminant transport and miscible displacement[END_REF][START_REF] Arbogast | A characteristic-mixed ÿnite element method for advection-dominated transport problems[END_REF] and Yang in [START_REF] Yang | A characteristic-mixed method with dynamic ÿnite element space for convection-dominated di usion problems[END_REF], and can be viewed as a procedure of ELLAM type [START_REF] Arbogast | A characteristic-mixed ÿnite element method for advection-dominated transport problems[END_REF]. It is also based on the space-time weak formulation (3.10), but uses a mixed ÿnite element approach by introducing the di usive ux z = -D as a new variable. Let V h × W h be the lowest-order Raviart-Thomas spaces [START_REF] Raviart | A mixed ÿnite element method for second order elliptic problems[END_REF], and Ŵ h be the space of discontinuous piecewise-linear functions on the same partition. Then, the CMFEM scheme can be formulated as follows: ÿnd c(x; t m+1 ) ∈ W h and z(x; t m+1 ) ∈ V h such that (x) c(x; t m+1 ) -ĉ(x * ; t m ) t w(x) dx + • z(x; t m+1 )w(x) dx = ( c -c)q(x; t m+1 )w(x) dx; ∀w ∈ W h ; D -1 z(x; t m+1 ) dx -c(x; t m+1 ) • C dx = 0; ∀C ∈ V h ; where x * is deÿned in (3. 8) and cˆ(x;t m ) ∈ W ˆ h is a post-processing of c(x;t m ) and z(x;t m ) deÿned by (x) ĉ(x; t m )w(x) dx = c(x; t m )w(x) dx; ∀w ∈ W h ; ŵ(x) • (D ĉ)(x; t m ) dx = -z(x; t m ) ŵ(x) dx; ∀ ŵ ∈ Ŵ h :

It is well known that in the mixed method, the scalar variable c(x; t m+1 ) is of ÿrst order accuracy in space. This post-processing procedure is used to improve the accuracy to the order of O(( x) 3=2 ) [START_REF] Arbogast | A characteristic-mixed ÿnite element method for advection-dominated transport problems[END_REF].

Remark 7. Theoretically the CMFEM is locally mass conservative. The situation might not be so clear numerically due to the following reasons: (i) The post-processing procedure is anti-di usive and, hence, could yield ĉ with undershoot or overshoot. A slope limiter has been used in the implementation of CMFEM to overcome this problem [START_REF] Arbogast | A characteristic-mixed method for contaminant transport and miscible displacement[END_REF]. It is not clear how the local mass conservation is achieved in this case. (ii) The CMFEM inherently requires a backtracking procedure and thus has to exactly determine the backtracked image at the previous time step t m of each cell at the future time step t m+1 in order to conserve mass. Since the backtracked image of each cell typically has curved boundaries in general, it is not clear how to trace these cell boundaries exactly to conserve mass numerically. Finally, the theoretically proved error estimate for the CMFEM is obtained only for Eq. (1.2) with a periodic boundary condition and is of O(( x) 3=2 ) which is suboptimal by a factor O(( x) 1=2 ).

Characteristic methods for immiscible uid ows, operator splitting techniques

In the governing equation (1.5) for immiscible ows, the hyperbolic part is given by Eq. (2.11) with a typically S-shaped function of the unknown, while the unknown function is a decreasing function in space. Hence, Eq. (2.11) could develop a non-unique solution [START_REF] Buckley | Mechanism of uid displacement in sands[END_REF][START_REF] Finlayson | Numerical Methods for Problems with Moving Fronts[END_REF][START_REF] Lax | Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves[END_REF][START_REF] Leveque | Numerical Methods for Conservation Laws[END_REF]. Thus, characteristic methods do not apply directly. Espedal and Ewing [START_REF] Espedal | Characteristic Petrov-Galerkin sub-domain methods for two-phase immiscible ow[END_REF] presented an operator-splitting technique to overcome this di culty. The fractional ow function f(c) is split into an advective concave hull f(c) of f(c), which is linear in what would be the shock region of Eq. (2.11), and a residual anti-di usive part. The modiÿed advection PDE @c @t + @ f(c) @x = 0; x ∈ ; t ∈ [0; T ] yields the same entropy solution as the PDE (2.11), and thus deÿnes characteristic directions uniquely. The residual anti-di usive advection term is grouped with the di usion term in the governing PDE so that correct balance between nonlinear advection and di usion is obtained. Numerically, the PDE is solved by a quadratic Petrov-Galerkin FEM. This technique has been applied in numerical simulation for immiscible ow by Espedal, Ewing, and their collaborators [START_REF] Dahle | Characteristic adaptive sub-domain methods for reservoir ow problems[END_REF][START_REF] Espedal | Simulation techniques for multiphase and multicomponent ows[END_REF]. Subsequently, Ewing [START_REF] Ewing | Operator splitting and Eulerian-Lagrangian localized adjoint methods for multi-phase ow[END_REF] and Dahle et al. have applied the operator-splitting technique to develop an ELLAM scheme for nonlinear advection-di usion PDEs, which has shown very promising results.

  The ELLAM scheme(3.14) symmetrizes the transport PDE (1.2), and generates accurate numerical solutions without excessive numerical di usion or nonphysical oscillation even if coarse spatial grids and large time steps are used

		• (Dc))(x; t m+1 ) dx
	=	(x)c(x; t m )w(x; t m + ) dx + t ( cqw)(x; t m+1 ) dx:	(3.14)
	Remark 5.


Cyril Voyant
email: voyant@univ.corse.fr

Gilles Notton
email: notton@univ.corse.fr

Marie Laure Nivet
email: nivet@univ.corse.fr

Fabrice Motte
email: motte@univ.corse.fr

Alexis Fouilloy
email: fouilloy@univ.corse.fr

Christophe Paoli

Bounded global irradiation prediction based on multilayer perceptron and time series formalism

As global solar radiation forecasting is a very important challenge, several methods are devoted to this goal with different levels of accuracy and confidence. In this study we propose to better understand how the uncertainty is propagated in the context of global radiation time series forecasting using machine learning. Indeed we propose to decompose the error considering four kinds of uncertainties: the error due to the measurement, the variability of time series, the machine learning uncertainty and the error related to the horizon. All these components of the error allow to determinate a global uncertainty generating a prediction bands related to the prediction efficiency. We also have defined a reliability index which could be very interesting for the grid manager in order to estimate the validity of predictions. We have experimented this method on a multilayer perceptron which is a popular machine learning technique. We have shown that the global error and its components are essentials to quantify in order to estimate the reliability of the model outputs. The described method has been successfully applied to four meteorological stations in Mediterranean area.

Incertainty;time series;confidence bands;prediction interval

I.

INTRODUCTION

Solar radiation is one of the principal energy sources for physical, biological and chemical processes, occupying the most important role in many engineering applications [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF]. The process of converting sunlight to electricity without combustion allows to create power without pollution. The major problem of such energy source is its intermittence and its stochastic character which make difficult their management into an electrical network [START_REF] Lorenz | Benchmarking of different approaches to forecast solar irradiance[END_REF]; Thereby, the development of forecasting models is necessary to use ideally this technology. By considering their effectiveness, it will be possible for example to identify the most optimal locations for developing a solar power project or to maintain the grid stability and security of a power management system [START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF]. Thus the solar energy forecasting is a process used to predict the amount of solar energy available for various time horizons [START_REF] Voyant | 24-hours ahead global irradiation forecasting using Multi-Layer Perceptron[END_REF]. Several methods have been developed by experts around the world and the mathematical formalism of Times Series (TS) has been often used for the short term forecasting (among 6 hours ahead) [START_REF] Diagne | Review of solar irradiance forecasting methods and a proposition for small-scale insular grids[END_REF]. TS is a set of ordered numbers that measures some activities over time. It is the historical record of global horizontal irradiance with measurements taken at equally spaced intervals with a consistency in the activity and the method of measurement. Some of the best predictors found in literature are Autoregressive and moving average, Bayesian inferences, Markov chains, k-Nearest-Neighbors predictors, support vector machine, regression tree, or artificial neural network (ANN). All these approaches are related to the machine learning application. The most often used is the last presented method: the artificial neural network and particularly the multilayer perceptron (MLP) [START_REF] Voyant | Optimization of an artificial neural network dedicated to the multivariate forecasting of daily global radiation[END_REF]. In the present study, we focus on this prediction method, the goal being to detail the uncertainties related to the global radiation prediction. These uncertainties can be decomposed into several components that will be explained and developed.

II. DATA

In Corsica Island, the data used to build the models are GHI measured in the meteorological stations of Ajaccio (41°55'N, 8°44'E, 4m asl) and Bastia (42°42'N, 9°27'E, 10m asl). They are located near the Mediterranean Sea and nearby mountains (1000 m altitude at 40 km from the sites). The data representing the global horizontal solar radiation were measured on an hourly basis from 1998 to 1999 (exactly two years). The two last studied stations are Montpellier (43.6°N and 3.9°E, 2 m asl) and Marseille (43.4°N and 5.2°E, 5 m asl) concerning the years 2008 and 2009. All these stations are equipped with pyranometers (CM 11 from Kipp & Zonen). The choice of these particular places is explained by their closed geographical and orographical configurations. These stations are located near the Mediterranean Sea and mountains. This specific geographical configuration of the four French meteorological stations makes cloudness difficult to forecast. Mediterranean climate is characterized by hot summers with abundant sunshine and mild, dry and clear winters. Irradiance nighttime values are not being used, the first morning data forecast are operated with the day before evening data

III. PREDICTION METHODOLOGY AND ERROR DECOMPOSITION

We chose to develop error propagation in the GHI prediction for the most common used predictor: the MLP. The base of this model is the time series approach (TS). A TS x(t) can be defined by a linear or non-linear model called fn (see Equation 1 where t = n, n-1,…, p+1, p with n, the number of observations and p the number of parameters of the model ; n ≫ p; h is the horizon of prediction and ϵ_(t+h) the committed error) [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF].

x(t+h)=fn (x(t),x(t-1)….,x(t-p+1))+ϵ(t+h) Eq 1

To estimate the fn model, a stationarity hypothesis is often necessary. This condition usually implies a stable process [START_REF] Hornik | Multilayer feedforward networks are universal approximators[END_REF]. This notion is directly linked to the fact that whether certain feature such as mean or variance change over time or remain constant. Previous studies [START_REF] Voyant | Predictability of PV power grid performance on insular sites without weather stations : use of artificial neural networks[END_REF][START_REF] Paoli | Solar Radiation Forecasting Using Ad-Hoc Time Series Preprocessing and Neural Networks[END_REF][START_REF] Voyant | Multihorizon solar radiation forecasting for Mediterranean locations using time series models[END_REF][11] show that the use of clear sky index (CSI) allows to make stationary the time series and so to correctly use the MLP forecasting.

A. Stationary process

In previous studies [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF][START_REF] Paoli | Use of Exogenous Data to Improve Artificial Networks Dedicated to Daily Global Radiation Forecasting[END_REF], it was demonstrated that the clear sky index calculated with the simplified Solis model [START_REF] Mueller | Rethinking satellitebased solar irradiance modelling: The SOLIS clear-sky module[END_REF] is the most reliable for our locations. The Solis model generates a clear sky hourly irradiation (CS) expressed by Eq. ( 2), the use of this model requires fitting parameter (g), extraterrestrial radiation (I0), solar elevation (h) and total measured atmospheric optical depth ():

𝐶𝑆(𝑡) = 𝐼 0 (𝑡). 𝑒𝑥𝑝 ( -𝜏 𝑠𝑖𝑛 𝑔 (ℎ(𝑡))
) . 𝑠𝑖𝑛 (ℎ(𝑡))

Eq 2

The simplified "Solis clear sky" model is based on radiative transfer calculations and the Lambert-Beer relation [START_REF] Ineichen | A broadband simplified version of the Solis clear sky model[END_REF]. The expression of the atmospheric transmittance is valid with polychromatic radiations, however when dealing with global radiation, the Lambert-Beer relation is only an approximation because of the back scattered effects. According to [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF] this model remains a good fitting function of the global horizontal radiation. The new computed time series (CSI) can be directly used with the MLP forecasting and is described be the equation 3:

CSI(t)=GHI(t)/CS(t) Eq 3

B. MLP prediction

Although a large range of different architectures of ANN is available [START_REF] Paoli | Solar Radiation Forecasting Using Ad-Hoc Time Series Preprocessing and Neural Networks[END_REF], MultiLayer Perceptron (MLP) remains the most popular [35]. In particular, feed-forward MLP networks with two layers (one hidden layer and one output layer) are often used for modeling and forecasting time series. Several studies [START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF] validated this approach based on ANN for the non-linear modeling of time series. To forecast the time series, a fixed number p of past values are set as inputs of the MLP, the output is the prediction of a future value [START_REF] Kalogirou | Artificial neural networks in renewable energy systems applications: a review[END_REF]. Considering the initial time series equation (Equation 1), this equation can be adapted to the non-linear case of one hidden layer MLP with b related to the biases, f and g to the activation function of the output and hidden layer, and to the weights. The number of hidden nodes (H) and the number of the input node (In) allow to detail this transformation. The number of layer 1 and 2 is given in superscript. (Equation 4):

𝐶𝑆𝐼 ̂(𝑡 + 1) = 𝑓(∑ 𝑦 𝑖 𝐻 𝑖=1 𝜔 𝑖 2 + 𝑏 2 ) with 𝑦 𝑖 = 𝑔(∑ 𝐶𝑆𝐼(𝑡 -𝑗 + 1) 𝐼𝑛 𝑗=1 𝜔 𝑖𝑗 1 + 𝑏 𝑖 1 ) Eq4
In the presented study, the MLP has been computed with the Matlab© software and its Neural Network toolbox. The characteristics chosen and related to previous work are the following: one hidden layer, the activation functions are the continuously and differentiable hyperbolic tangent (hidden) and linear (output), the Levenberg-Marquardt learning algorithm with a max fail parameter before stopping training equal to 5.

This algorithm is an approximation to the Newton's method. The prediction of the GHI is obtained using the equation:

(GHI) (t+1)=(CSI) (t+1).CS(t+1) Eq 5
To customize the input layer of the MLP we choose the use of the mutual information to determine In as described in [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF][START_REF] Voyant | Statistical parameters as a means to a priori assess the accuracy of solar forecasting models[END_REF][START_REF] Voyant | Numerical weather prediction or stochastic modeling: an objective criterion of choice for the global radiation forecasting[END_REF][START_REF] Parviz | Mutual Information Based Input Variable Selection Algorithm and Wavelet Neural Network for Time Series Prediction[END_REF][START_REF] Jiang | Mutual information algorithms[END_REF]. According the results obtained in these papers, we use H equal to In for all the experiments conducted in this study. Furthermore in order improve the learning of the MLP, it is a common practice to filter out the data removing night hours. Indeed we consider only periods between sunrise and sunset [41,42]. We have chosen to apply a selection criterion based on the solar zenith angle (SZA): solar radiation data for which the solar zenith angle is greater than 80° have been removed [START_REF] Lauret | A benchmarking of machine learning techniques for solar radiation forecasting in an insular context[END_REF]. This transformation is equivalent to a filtering related to the solar elevation angle lower than 20.

C. Error decomposition

In these section, we propose to decompose the error considering four kinds of uncertainties: the error due to the measurement, the error due to the variability of the time series, the error related to the machine learning uncertainty and the error related to the horizon.

In our assumption, all the previous  terms are independent random variables that are normally distributed (and therefore also jointly so), then their sum is also normally distributed and the global form of the standard deviation 𝜎 𝑡𝑜𝑡 (𝑡 + 1) becomes (in this equation the quick fluctuations are taken into account with 𝜎 𝑣𝑎𝑟 , but it is also possible to consider 𝜎 𝑖𝑛ℎ ):

𝜎 𝑡𝑜𝑡 (𝑡 + 1) = √( 𝜎 𝑚𝑒𝑎𝑠 ) 2 + (𝜎 𝑠𝑎𝑚𝑝 (𝑡 + 1)) 2 + (𝜎 𝑖𝑛𝑖 (𝑡 + 1)) 2 + (𝜎 𝑣𝑎𝑟 (𝑡 + 1)) 2 Eq 6
Considering that there is a persistence of the variability for a short horizon, 𝜎 𝑣𝑎𝑟 (𝑡) = 𝜎 𝑣𝑎𝑟 (𝑡 + 1)

thus:

𝜎 𝑡𝑜𝑡 (𝑡 + 1) = √( 𝜎 𝑚𝑒𝑎𝑠 ) 2 + (𝜎 𝑠𝑎𝑚𝑝 (𝑡 + 1))

2 + (𝜎 𝑖𝑛𝑖 (𝑡 + 1)) 2 + (𝜎 𝑣𝑎𝑟 (𝑡)) 2

Eq 7

With

𝜎 𝑚𝑒𝑎𝑠 = 𝜎(𝐶𝑆𝐼)(∑ (∑ 𝜔 𝑖 2 𝜔 𝑖𝑗 1 𝑁𝑐 𝑖=1 ) 2 𝑁𝑒 𝑗=1
)

1/2 , 𝜎 𝑠𝑎𝑚𝑝 and 𝜎 𝑖𝑛𝑖 are computed respectively with k-fold and 50 random initializations and 𝜎 𝑣𝑎𝑟 = 𝐺𝐻𝐼 ̂(𝑡 + 1). 𝑔(𝑉𝑜𝑙2(𝑡)). For an easier computing, it is also possible to use 𝜎 𝑖𝑛ℎ replacing 𝜎 𝑣𝑎𝑟 (𝑡) with a less robust result but not dependent on the instant of the prediction.

𝜎 𝑡𝑜𝑡 (𝑡 + 1) = √( 𝜎 𝑚𝑒𝑎𝑠 ) 2 + (𝜎 𝑠𝑎𝑚𝑝 (𝑡 + 1)) 2 + (𝜎 𝑖𝑛𝑖 (𝑡 + 1)) 2 + (𝜎 𝑖𝑛ℎ ) 2
Eq 8

It is possible to define a prediction band taking into account all the uncertainties (Eq 9).

𝐺𝐻𝐼 ̂(𝑡 + 1) = 𝐺𝐻𝐼 ̂𝑀𝐿𝑃 (𝑡 + 1) ± 𝜎 𝑡𝑜𝑡 (𝑡 + 1)

Eq 9

Such prediction intervals were often proposed in the literature [START_REF] David | Probabilistic forecasting of the solar irradiance with recursive ARMA and GARCH models[END_REF][START_REF] Join | Solar energy production: Short-term forecasting and risk management[END_REF][START_REF] Trapero | Calculation of solar irradiation prediction intervals combining volatility and kernel density estimates[END_REF]; they refer to machine learning method (𝜎 𝑀𝐿 (𝑡 + 1)) 2 = (𝜎 𝑠𝑎𝑚𝑝 (𝑡 + 1))

2 + (𝜎 𝑖𝑛𝑖 (𝑡 + 1)) 2 [START_REF] David | Probabilistic forecasting of the solar irradiance with recursive ARMA and GARCH models[END_REF][START_REF] Trapero | Calculation of solar irradiation prediction intervals combining volatility and kernel density estimates[END_REF] or to volatility and 𝜎 𝑣𝑎𝑟 (𝑡) [START_REF] Join | Solar energy production: Short-term forecasting and risk management[END_REF] but rarely both to the two kinds of uncertainty and never concerning 𝜎 𝑚𝑒𝑎𝑠 . Note that in the case of other machine learning method used the term 𝜎 𝑖𝑛𝑖 can be equal to zero (e.g. support vector regression, regression tree etc.). The ideal case would be to systematically propose a confidence interval of prediction related to the three sorts of uncertainty (with 𝜎 𝑇𝑆 = 𝜎 𝑣𝑎𝑟 (𝑡)𝑜𝑟𝜎 𝑖𝑛ℎ considering the desired reliability).

𝜎 𝑡𝑜𝑡 (𝑡 + 1) = √(𝜎 𝑚𝑒𝑎𝑠 ) 2 + (𝜎 𝑀𝐿 (𝑡 + 1)) 2 + (𝜎 𝑇𝑆 (𝑡)) 2 Eq 10 Now, considering the horizon of prediction, we define the new global uncertainty with the equation 30 with 𝜎 ℎ𝑜𝑟 = 𝐺𝐻𝐼 ̂(𝑡 + ℎ). 𝛼(ℎ).

𝜎 𝑡𝑜𝑡 (𝑡 + ℎ) = √(𝜎 𝑚𝑒𝑎𝑠 ) 2 + (𝜎 𝑀𝐿 (𝑡 + 1)) 2 + (𝜎 𝑇𝑆 ) 2 + (𝜎 ℎ𝑜𝑟 (ℎ)) 2 Eq 11 IV. RESULTS

The previous components allows to calculate the global uncertainty and to propose two prediction bands: UB for upper band and LB lower band [START_REF] Gneiting | Probabilistic Forecasting[END_REF]. Thus, the quality of the prediction can be defined by the triplet {𝐺𝐻𝐼 ̂(𝑡 + ℎ); 𝐿𝐵; 𝑈𝐵} [START_REF] Wilks | Statistical Methods in the Atmospheric Sciences An Introduction[END_REF]. We can also estimate the reliability of the prediction considering that the prediction is efficient when UB-LB is very lower than 𝐺𝐻𝐼 ̂(𝑡 + ℎ) and inefficient when UB-LB is equal or upper to 𝐺𝐻𝐼 ̂(𝑡 + ℎ) value. From this hypothesis, we can define the reliability  as t(UB(t+1)-LB(t+1))𝐺𝐻𝐼 ̂(𝑡 + ℎ) Lower is this parameter, more efficient is the prediction. We construct a reliability index between 0 and 1 considering that if (UB(t+1)-LB(t+1))𝐺𝐻𝐼 ̂(𝑡 + ℎ) > 1thent i.e. the prediction is not sureThe final prediction becomes: We can see that 𝜎 𝑚𝑒𝑎𝑠 is the parameter the less interesting for the bands construction and that it is necessary to consider the coupling of 𝜎 𝑖𝑛ℎ and 𝜎 𝑀𝐿 (related to 𝜎 𝑠𝑎𝑚𝑝 𝑎𝑛𝑑𝜎 𝑖𝑛𝑖 ) for a good prediction interval definition. For other sites the obtained curve are similar and no more information is observed. In the figure 6, the top curve (same prediction configuration that previously) compared the average prediction 𝐺𝐻𝐼 ̂(𝑡 + ℎ) (marks) versus the GHI measurement (line). The bottom curve shows the associated reliability index (t We see that when the variability is low (two first day from 3711 to 3726) the reliability is important (close to 70%) but when cloud occurs the value is much lower and can reach 0%..

V. CONCLUSIONS

In this paper we have shown that it is possible to compute a prediction band in the context of global radiation time series forecasting using machine learning. We have define for a popular machine learning technique, the multilayer perceptron, four kinds of uncertainties: the error due to the measurement, the variability of time series, the machine learning uncertainty (initialization and sampling) and the error related to the horizon. In literature, rarely both to the two kinds of uncertainty 𝜎 𝑖𝑛𝑖 and 𝜎 𝑀𝐿 are studied, and never 𝜎 𝑚𝑒𝑎𝑠 . We have also defined a reliability index which could be very interesting for the grid manager in order to estimate the validity of predictions. The described method has been successfully applied to five meteorological stations in Mediterranean area. We are sure that it is possible to generalize the approach to other sites and other machine learning tools. Thereby in future, we will try to apply the methodology to other time granularities and predictor as SVM, regression tree or random forest.

-

  𝐺𝐻𝐼 ̂(𝑡 + ℎ)= 〈𝐺𝐻𝐼 ̂(𝑡 + ℎ)〉, average of 50 simulations (50 training and initialization weights, 50 different training sets) LB= √( 𝜎 𝑚𝑒𝑎𝑠 ) 2 + (𝐺𝐻𝐼 ̂𝑚𝑖𝑛 (𝑡 + ℎ) -〈𝐺𝐻𝐼 ̂(𝑡 + ℎ)〉) 2 + (𝜎 𝑖𝑛ℎ ) 2 + (𝜎 ℎ𝑜𝑟 (ℎ)) 2 Eq 34 UB= √( 𝜎 𝑚𝑒𝑎𝑠 ) 2 + (𝐺𝐻𝐼 ̂𝑚𝑎𝑥 (𝑡 + ℎ) -〈𝐺𝐻𝐼 ̂(𝑡 + ℎ)〉) 2 + (𝜎 𝑖𝑛ℎ ) 2 + (𝜎 ℎ𝑜𝑟 (ℎ)) 2 Eq 35 With 𝜎 𝑚𝑒𝑎𝑠 (𝐺𝐻𝐼 ̂(𝑡 + 1)) = 1%. 𝐺𝐻𝐼 ̂(𝑡 + 1), 𝜎 𝑖𝑛ℎ = 𝐺𝐻𝐼 ̂(𝑡 + 1). 𝑛𝑅𝑀𝑆𝐸(𝐺𝐻𝐼 𝑡𝑟𝑒𝑛𝑑 (𝑡) -𝐺𝐻𝐼(𝑡)), 𝜎 ℎ𝑜𝑟 (ℎ) == 𝐺𝐻𝐼 ̂(𝑡 + ℎ). 𝛼(ℎ) and 𝐺𝐻𝐼 ̂𝑚𝑖𝑛/𝑚𝑎𝑥 (𝑡 + ℎ) are the min and max values of the 50 predictions generated with 50 simulations. The figure 5 shows for Ajaccio an example of the prediction bands, considering all the kind of uncertainty with horizon h=1 hour. Line represents measurement and dashed lines the upper and lower bands concerning each kind of uncertainties.

Figure 1 .

 1 Figure 1. Uncertainty in the GHI predictions for the horizon h=1 for Ajaccio case

Figure 2 .

 2 Figure 2. Comparison for the horizon h=1 for Ajaccio of GHI predictions (mark) and GHI measurement (line) on the top and associated reliability index in the bottom

ACKNOWLEDGMENT

This work is supported by European Union's Horizon 2020 research and innovation programme under grant agreement No 646529 through the TILOS project (Technology innovation for the Local Scale, Optimum integration of Battery Energy Storage). The data sets were provides by the CNRS UMR SPE 6134 Laboratory..