Lattice-Boltzmann coupled models for advection–diffusion flow on a wide range of Péclet numbers

Traditional Lattice-Boltzmann modelling of advection–diffusion flow is affected by numerical instability if the advective term becomes dominant over the diffusive ( i.e. , high-Péclet flow). To overcome the problem, two 3D one-way coupled models are proposed. In a traditional model, a Lattice-Boltzmann Navier–Stokes solver is coupled to a Lattice-Boltzmann advection–diffusion model. In a novel model, the Lattice-Boltzmann Navier– Stokes solver is coupled to an explicit finite-difference algorithm for advection–diffusion. The finite-difference algorithm also includes a novel approach to mitigate the numerical diffusivity connected with the upwind differentiation scheme. The models are validated using two non-trivial benchmarks, which includes discontinuous initial conditions and the case Pe g → ∞ for the first time, where Pe g is the grid Péclet number. The evaluation of Pe g alongside Pe is discussed. Accuracy, stability and the order of convergence are assessed for a wide range of Péclet numbers. Recommendations are then given as to which model to select depending on the value Pe g —in particular, it is shown that the coupled finite-difference/Lattice-Boltzmann provide stable solutions in the case Pe → ∞ , Pe g → ∞ .


Introduction
Advection-diffusion occurs in a wide range of fluid-flow problems, including transport of chemicals, natural convection heat transfer when thermal expansion can be ignored, assessment of mixing through a passive scalar tracer, and many more.
Over the last decades, numerical modelling -in particular, computational fluid dynamics (CFD) -has been extensively used as a cheap but robust alternative to often lengthy and expensive experiments [1]. The lattice-Boltzmann is a particular type of explicit finite-difference method, with a tuneable diffusivity parameter [2]; this characteristics provides an advantage over traditional explicit finite-difference methods, as it allows the use of much larger timesteps than what would otherwise required to maintain stability. More significantly, it represents a valid alternative to implicit/segregated finite-volume CFD methods traditionally used to solve the Navier-Stokes equations because of a number of advantages in terms of numerical efficiency and parallelizability, viz.: (i) full explicitness, with no internal loop being required-therefore, every timestep is updated through a limited and well-defined number of floating-point operations; (ii) separation between non-linear and non-local part of the algorithm with the latter * Correspondence to: B2.20, Chesham Building, Bradford BD7 1DP, United Kingdom. E-mail address: ddapelo@bradford.ac.uk (D. Dapelo).
being usually limited to first neighbour-access, thereby allowing large parallel runs with no significant efficiency loss due to non-scalable inter-processor communication; and (iii) structure simplicity, allowing relatively simple extensions towards a wide variety of phenomena [2][3][4]. [5] showed that a Lattice-Boltzmann model for laboratory-scale gas mixing on a non-Newtonian fluid performed around 180 times faster than a finite-volume analogue [6], and ran on ten times more processors with no apparent loss of efficiency.
Whilst the Lattice-Boltzmann method has been applied to 1D and 2D advection-diffusion problems [7][8][9], the models demonstrate numerical instability in advection-dominant flow (high Péclet number) as the relaxation time approaches the critical value of 1/2. Proposed solutions valid for the 1D, 2D and 3D [10][11][12] mitigate the problem, but still exhibit instability for Pe ⪆ 10 4 . A more recent fractionalstep method [13] requires doubling the number of timesteps and the number of lattice points across the coordinate directions; furthermore, it performs calculations requiring access to second-neighbour lattice sites, thereby compromising parallelizability. Therefore, there is clear a need for an advection-diffusion Lattice-Boltzmann method with the https://doi.org/ 10 Numerical diffusion term connected to the upwind differentiation scheme for advection, m −3 Unit vector along the positive direction of the coordinate Navier-Stokes one-particle density function, kg m −3 Navier-Stokes discretized one-particle density function relative to the th lattice direction, kg m −3

(eq)
Navier-Stokes equilibrium one-particle density function, kg m −3 characteristics of high Péclet number stability, and requiring access to first-neighbouring cells only.
Within this article, two Lattice-Boltzmann models are proposed to reproduce a diffusive scalar field concentration advected by a hydrodynamic velocity field. In both the models, hydrodynamics and (eq) Navier-Stokes discretized equilibrium one-particle density function relative to the th lattice direction, kg m −3 Advection-diffusion discretized one-particle density function relative to the th lattice direction, kg m −3

(eq)
Advection-diffusion discretized equilibrium one- Millions of lattice updates per processor per second advection-diffusion are simulated in a coupled way. The two models are composed of a Lattice-Boltzmann submodel to solve hydrodynamics, and a second submodel for advection-diffusion. The coupling occurs one-way, viz.: the solution to the hydrodynamics subproblem influences the advection-diffusion solution, but the converse does not occur. In the ''traditional'' model, the advection-diffusion submodel is defined as a Lattice-Boltzmann model, whilst in the ''novel'' model, an explicit finite-difference submodel is used. The rationale for introducing an explicit finite-difference submodel for advection-diffusion is based on the fact that finite-difference is the most obvious choice of model when the computational domain is composed of a cubic lattice (as is the case of the numerical work presented here), and, for the very specific case of an explicit method with access to first neighbours only, it maintains the same level of parallelizability and efficiency as the standard Lattice-Boltzmann method. Scaling tests are performed in order to verify this claim. In both models, at each timestep, the Lattice-Boltzmann submodel solves the Navier-Stokes equations to update the hydrodynamics; the coupling is performed by passing the obtained velocity field to the advection-diffusion submodel; and finally, the latter solves the advection-diffusion equation and updates the scalar field concentration.
In the literature cited above, the assessments are performed in terms of the Péclet number, or other numbers tracing the balance between advective and a given form of diffusive (chemical, or thermal) transport. This approach provides a good assessment of the physics underlying the phenomena being analysed, but does not necessarily provide a reliable assessment of numerical stability. [9] defines the grid Reynolds number Re g from lattice quantities, and identifies that a simulation is stable if Re g ≳ 10, heuristically meaning that ''The lattice should always be sufficiently fine to resolve local vortices''. [14,15] introduce the grid Péclet number Pe g as the analogue to Re g in advection-diffusion models. This approach implies a separation between physical interpretation of a process (i.e., the balance between advective and diffusive transport), which is expressed through ℜ, Pe, . . . ; and stability assessment of a numerical simulation (i.e., whether or not the lattice is fine enough to resolve the local vortices), which is expressed through Re g , Pe g , . . . . In parallel with the assessment of Pe g , the evaluation of the numerical stability is also performed in [9,14,15] through methodologies involving the mutual balance between Re g or Pe g , the lattice relaxation time (or in the advection-diffusion) and the characteristic velocity or the courant number Co ≡ ∕ (where is the lattice spacing and the timestep). The above-mentioned work does not explicitly consider the case of zero diffusivity, i.e.Pe g → ∞ ; in fact, it can be implied from [9] and [14] that the simulations are intrinsically unstable in this specific case. However, an explicit assessment of the case Pe g → ∞ is justified by its importance in applications, e.g. in assessing mixing quality through scalar tracer concentration, or wherever diffusive processes are negligible altogether. For this reason, the work described within this article is focussed on the assessment of model stability in terms of Pe g alone, and the case Pe, Pe g → ∞ is explicitly taken into consideration for the first time within the research concerning the Lattice-Boltzmann method; the balance between Pe g , Co and relaxation time is considered only for the purpose of verifying the coherence of the simulations. In order to provide the most possible robust assessment, two benchmarks are implemented in a non-trivial way as a Gaussian concentration packet or a sinusoidal wave spreading and being advected by a constant velocity field, with a tracer concentration presenting a sharp discontinuity as initial condition; within both periodical and non-periodical computational domains. Previous work [16] shows that Lattice-Boltzmann advection-diffusion observes second-order convergence on simulations starting from continuous initial conditions, and first-order if the initial condition presents a discontinuity. This reduction of the order of convergence is explained by the fact that, in the discontinuous case, the eigenvalues of the Lattice-Boltzmann evolution matrix are close to −1, thereby causing the onset of oscillations. However, [16] is limited in analysis to 1D simulations: within this work, the analysis is extended to the 3D case, thereby providing a complete assessment of the convergence for the model and filling the knowledge gap in the literature.
This article is structured as follows. In Section 2, the two models are outlined. In Section 3, the Navier-Stokes Lattice-Boltzmann model and the advection-diffusion Lattice-Boltzmann are described. Then, the explicit finite-difference model is described in Section 4. Two benchmarks to assess performance and stability are described in Section 5. The results are reported in Section 6, and recommendations on which model to use depending on the value of Pe g are given in Section 7. Finally, conclusions are drawn in Section 8.

General outline of the model
Both the proposed coupled and the standard models follow the general algorithm described below.

Lattice-Boltzmann modelling
The Lattice-Boltzmann model is a mesoscopic model insofar as it aims at finding a trajectory ( , , ) in the phase space The ''one-particle density function'' ( , , ) describes the probability of finding one particle of fluid within the elemental cube ( , + ) with a velocity comprised within the cubic interval ( , + ) at the time . The observable fields (viz., density, velocity and shear stress) correspond to the particle density function's zeroth, first and second moments respectively [9]: The model solves the Boltzmann equation: Algorithm 1: General algorithm 1 Initialize hydrodynamics (density field and velocity field ) and advection-diffusion characteristics (scalar concentration field ); 2 for timestep: = 0 to max by do 3 Hydrodynamics update: a standard Lattice-Boltzmann submodel solves the weakly-compressible Navier-Stokes equations: (with being the pressure and the dynamic viscosity) and updates the velocity field ;

4
Pass the updated velocity field to the submodel for advection-diffusion; 5 Advection-diffusion update: an explicit finite-difference submodel (or an advection-diffusion Lattice-Boltzmann submodel for the comparison model) solves the advection-diffusion equation: (with being the diffusivity) and updates the concentration field ; Eq. (7) describes continuity for in the phase space, plus a source-sink term arising from inter-particle collisions. In the diluted gas hypothesis, only binary particle collisions contribute to the collision operator . Under the isotropy assumption, it is possible to model the potentially extremely complicated operator [18] as an average, relaxing towards the equilibrium particle density function (eq) : where is the relaxation time, and (eq) is the Maxwell equilibrium distribution [9]: s is the speed of sound, and density and velocity are evaluated through Eqs. (4) and (5). Space is discretized as a cubic lattice of lattice size , time is subdivided into timesteps of interval , and the three-dimensional velocity space is discretized as a finite set of vectors pointing to the zeroth, first, second and third neighbours of a generic lattice site, with moduli respectively √ 0, √ 1, √ 2 and √ 3 times ∕ . A given discretization is conventionally identified by a tag D Q , with the dimensionality of the problem and the dimension of the discrete velocity space. ( , , ) is redefined over the discretized lattice and timesteps, and rewritten as a set ( , ), each indicating the probability of finding a particle with the corresponding discretized velocity . The integrals in Eqs. (4), (5) and (6) are redefined as summations over the velocity set: The error arising from the discretization of the velocity set is cancelled by writing the Maxwell equilibrium function (Eq. (9)) in terms of orthonormal Hermite polynomials up to the second order: where the weights and the speed of sound are defined in a standard way for the specific D Q lattice, with the latter usually being defined as: The Boltzmann Eq. (7), considering the BGK assumption (Eq. (8)), becomes the Lattice-Boltzmann Equation: ( The updating process for each lattice site at a given timestep comprises two steps. First: a local, non-linear collision: Then: a linear, non-local streaming : In this way, the adiabatic dynamics with a compressibility error of Ma 2 with Ma being the Mach number, are recovered. In the limit Ma ≪ 1, the Lattice-Boltzmann Eq. (15) reproduces the Navier-Stokes (Eqs. (1)) and (2) in the incompressible limit [9], with pressure and kinematic viscosity defined as:

Lattice-Boltzmann model for the advection-diffusion equations
A comparison between the advection-diffusion (Eq. (3)) and the Navier-Stokes equation for the momentum (Eq. (2)) shows that the latter can be regarded as an advection-diffusion equation for the momentum density vector in place of the scalar concentration . This analogy allows the export of the same model developed in Section 3 for the Navier-Stokes case, to the advection-diffusion [9]. Alternatively, it is possible to derive the Lattice-Boltzmann advection-diffusion model directly from the macroscopic equation (3) following the approach of [19]. Either way, a discretized one-particle density function is conventionally defined as ( , ), and the concentration is obtained from its zeroth moment: The only difference from the Navier-Stokes case is that higher-order moments of do not possess physical meaning because only one single scalar conservation law (viz., mass) is approximated in the advection diffusion equation (Eq. (3)), and therefore are not evaluated. In fact, the physical advection velocity appearing in Eq. (3) is imposed externally (viz., copied from the hydrodynamics model). The rest of the model is left unchanged: i.e.the definition of speed of sound as a numerical constant (Eq. (14)) is maintained unaltered, the Lattice-Boltzmann equation is defined as in Eq. (15): (with an advection-diffusion relaxation time ), the equilibrium function follows Eq. (13) with in place of : and the diffusivity takes the place of the kinematic viscosity in Eq. (18): The same lattice of the Lattice-Boltzmann Navier-Stokes model is used in order to avoid interpolation, and with the same timestep to simplify the time iteration. Consequently, the values of the speed of found for the advection-diffusion and the Navier-Stokes models are the same.

Stability constraint of the lattice-Boltzmann model for advectiondiffusion
The Péclet number is defined as: where and are respectively the length and velocity scales of a given problem. As the Reynolds number is the ratio of inertial to viscous forces in the Navier-Stokes case, Pe can be interpreted as the balance between advective and diffusive transport. The advection-diffusion relaxation time approaches the stability-critical value of 1∕2 when Pe grows, exactly as in the Navier-Stokes case if ℜ grows [13].
Following the analogy between Navier-Stokes and advectiondiffusion, a first criterion to assess advection-diffusion Lattice-Boltzmann can be outlined through defining the grid Péclet number following [9]: The stability criterion reads: Within this work, the dimensionless version * of the quantity is computed from its units of measure { } as follows: Considering that, within this work, we have = with being the number of cells comprised within , it follows that Pe g can be equivalently defined as: A second criterion for advection-diffusion Lattice-Boltzmann stability comes from the same analogy. [20] studied the neutral stability Co- * curves, and [9] provided a simplified mathematical expression to defined the stability domain of the Lattice-Boltzmann method. This expression can be extended to the advection-diffusion as follows: As the focus of this work is to compare the stability of the Lattice-Boltzmann and the explicit finite-difference advection-diffusion models in terms of Pe g and in the specific case Pe g → ∞, the criterion of Eq. (25) is mainly employed. The robustness of the subsequent conclusions is then assessed through the criterion of Eq. (28).

Finite-difference model for the advection-diffusion equations
An explicit finite-difference discretization [21] of the macroscopic advection-diffusion equation (3 is performed on the same lattice of the Lattice-Boltzmann Navier-Stokes model in order to avoid interpolation, and with the same timestep to simplify the time iteration: The ''new'' superscript indicates that the value of the concentration field is referred to the timestep being currently resolved; conversely, its absence indicates that the value of the concentration field refers to the previous timestep. This formulation allows an explicit evaluation of new without timestep-internal loops, as in the Lattice-Boltzmann case. The differentiation schemes [ ⋅ ] sch are chosen as: the standard second-order central scheme for the diffusion term; and either standard first-order upwind, or standard second-order central schemes for the advection term. In three dimensions, the central diffusion scheme reads: The central advection scheme reads: The upwind: with: The orders of convergence of the schemes above (Eqs. (30)-(32)) can be verified by substituting Taylor expansions of ( ±̂) around ( ) in the respective equations.

Stability constraints of the finite-difference model for advectiondiffusion
Substituting the appropriate differentiation schemes (Eqs. (30)-(32)), it is possible to rewrite Eq. (29) in the following general form: As a necessary condition for the algorithm to be stable, all the coefficients in the equation above should be non-negative. This requirement, in turn, poses constraints on the values of velocity, diffusivity, lattice spacing and timestep, as follows.
A first constraint comes from the requirement that 0 should be non-negative. This is related to the explicit nature of the model, and consists of a superior limit to the timestep value. Within this work, it is formulated as follows: Whilst this constraint usually poses serious limits to explicit methods, it does not prevent a successful application of the coupled model because in Lattice-Boltzmann modelling, the timestep value is usually kept low by the assumption ∕ ∼ Ma ≪ 1. Conversely, once and are chosen following the balance between accuracy, stability and numerical expense usually performed in Lattice-Boltzmann modelling [9], Eq. (35) poses a superior limit to the diffusivity.
A second constraint comes from the requirement that − should be non-negative in the case of the central advection scheme: This constraint poses an inferior limit to the ratio of diffusivity over reference velocity. A comparison between Eqs. (27) and (36) shows that: Therefore, the constraint of Eq. (36) is equivalent to a constraint on Pe g : As in Section 3.2, the interest is directed towards the stability criterion of Eq. (38) in order to compare the stability of the Lattice-Boltzmann and the explicit finite-difference advection-diffusion models in terms of Pe g and in the specific case Pe g → ∞. The criterion of Eq. (35) is used as a background verification of the applicability of the explicit finite-difference method.

Negative-diffusivity correction
Although Eq. (36) poses a strict stability constraint only if the central advection scheme is adopted, bespoke numerical results in the upwind advection case show the presence of numerical diffusivity along the direction of the flow, proportional to . This is due to the fact that the first truncated term in Eq. (32), corresponding to second-ordered Taylor expansions of ( ±̂) around ( ), takes the form of: A comparison of this term to the diffusion differentiation scheme (Eq. (30)) shows that this term can be interpreted as an artificial anisotropic diffusivity contribution. In order to mitigate this effect, a novel approach, consisting of multiplying the artificial diffusivity (Eq. (39)) by a tuning factor and subtracting it from the diffusivity scheme, is proposed. The secondorder central differentiation scheme for diffusivity is thus redefined as follows when employed in combination with the upwind scheme for advection:

Numerical work
A periodical and a non-periodical benchmark are defined. Over the periodical benchmark, simulations are performed through only the Lattice-Boltzmann advection-diffusion model, with both smooth and sharp initial conditions for the concentration field . Over the nonperiodical benchmark, simulations are performed through both the Lattice-Boltzmann and the finite-difference advection-diffusion models, but only with a sharp initial condition. As detailed below, the smoothstart simulations reproduce a sinusoid, and the sharp-start one a series of Gaussian packets.

Periodical benchmark
A simplified version of the general algorithm (Section 2) is appliedspecifically, the hydrodynamics is not solved, and a constant velocity field = (2.5, 2.5, 2.5) m∕s is passed to the advection-diffusion model. A time-averaged discrete relative 2 -error norm is calculated via: where: is the numerical outcome of the simulations, while an is the bespoke analytical solution.

Non-periodical benchmark
The full general algorithm is applied as described in Section 2. The computational domain (Fig. 1) consists of a square-section pipe with = = 1 m and ≫ being set depending on the specific run. The flow velocity is a constant with magnitude of 0, 0.1 or 2 m∕s, directed towards the direction and is once again set depending on the specific run. The boundary conditions for the Lattice-Boltzmann Navier-Stokes solver are: inlet and outlet are defined as in the figure; free-slip on the walls. For the advection-diffusion solvers, the boundary conditions on the wall are defined as: bounce-back for the Lattice-Boltzmann, no-penetration for the finite-difference. The potential error arising from setting bounce-back and no-penetration conditions to inlet and outlet is avoided by choosing ≫ , and by stopping the simulations before the Gaussian packet's 3 tail reached the outlet.
The Navier-Stokes lattice relaxation time is set to = 0.8. The timestep is determined under diffusive scaling through the prescription = 2 ⋅ 1 s∕m 2 . Unless otherwise specified, it is set = 11. It follows that Co amounts to 0, 0.009 and 0.182 for corresponding respectively to 0, 0.1 and 2.0 m∕s.
Simulation results and errors are evaluated at each timestep separately, in terms of the value of along the dashed line in Fig. 1  (indicated as below). Only the last timestep is considered for further analysis. The error is computed as the 2 distance between numerical and analytical values of over the 2 norm of the analytical solution, evaluated along the same dashed line: (43)

Smooth initial condition
The initial condition is set as: It can be demonstrated [22] that the benchmark possesses the following analytical solution: An example of solution is provided in Fig. 2.

Non-periodical
A spreading Gaussian density packet advecting along is modelled. The initial condition for is defined as: The initial position 0 is either 0 or ∕4 depending on the particular run. It can be demonstrated that the benchmark possesses the following analytical solution: and the initial condition for the analytical solution is the Dirac's delta: A proof of the mutual coherence of Eqs. (46), (47) and (48) is given in the Supplementary Material.

Periodical
A superposition of spreading Gaussian density packets advecting along is modelled. The initial condition for is defined as: It can be demonstrated [23] that the benchmark possesses the following analytical solution: and the initial condition for the analytical solution is the Dirac's comb: The proof of the mutual coherence of Eqs.

Results
A series of numerical simulations reproducing the smooth-start benchmark was performed on a Lenovo X380 laptop equipped with a 4-core, 8-thread Intel Core i5-8250U CPU. Another series, reproducing the sharp-start, was performed on a Dell XPS 13 9600 laptop equipped with a 4-core, 8-thread Intel Core i7-1065G7. A third series of simulations for testing purposes was computed on a maximum of 10 nodes with two deca-core Intel Xeon E5-2660 v3 each. Finally, the scaling-up tests were performed on three 40-core Lenovo ThinkSystem SR65 units.

Smooth start: Lattice-Boltzmann
The convergence behaviour of the Lattice-Boltzmann model for advection-diffusion on the smooth start case is reported in Fig. 3. The second-order convergence is evident.

Sharp start on the periodical benchmark: Lattice-Boltzmann
The convergence behaviour of the Lattice-Boltzmann model for advection-diffusion on the smooth start case is reported in Fig. 4. A deterioration of the convergence from second to first order can be observed, as noted by [16]. However, such deterioration occurs only for Pe g ≲ 100-that is, a relaxed version of the stability criterion of Eq. (25). In other words, the violation of the above-mentioned criterion leads to results that still retain stability, but at the price of a lower order of convergence.
Concerning the second stability criterion (Eq. (28)), it can be verified that the simulations with Pe = 10 000 and Pe = 100 000 fail the second line of Eq. (28), despite the fact that the simulations with Pe = 10 000 and Pe g ≲ still display second-order convergence.

Sharp start on the non-periodical lattice: Lattice-Boltzmann
The convergence behaviour of the Lattice-Boltzmann model for advection-diffusion with Pe g = 0 on the sharp start case is reported in Fig. 5. The value of diffusivity is set to = 0.01 m 2 ∕s. In this case, the first-order convergence is evident, as in [16].
In Fig. 6(a), the value of along the coordinate is plotted for = 0 and for values of comprised between 0.001 and 1 m 2 ∕s. The error is shown to be below 4%. However, Fig. 6(b) evidences the presence of oscillations for = 0 m 2 ∕s.
The results of tests performed on > 0, corresponding to Co = 0.00909, are shown in Fig. 7. The corresponding values of , Pe, Pe g and * are reported in Table 1. In all the cases, the computational time was between 20 and 25 CPUs. Fig. 7(a) shows the onset of oscillations for Pe ⪆ 400, which results in an increase in the error (Fig. 7(b)). From this value of Pe onwards, the simulations are also found to fail the stability criterion of Eq. (28). Oscillations dominate the figure for Pe → ∞ and eliminate the model's predictive power.

Sharp start: finite-difference
The convergence behaviour of the explicit finite-difference model is reported in Figs. 8 (for diffusive-only systems) and 9 (for > 0). In the case > 0, the same values of , Pe and Pe g of Section 6.3 and Table 1 are used. Nd = 50 in all the runs, while NU is infinite in the = 0 runs (Fig. 8) and spans from 2 to 11.6 in the > 0 case (Fig. 9). Therefore, the stability constraints reported in Eqs. (35) and (36) are met. Second-order convergence is evident for both the central and the upwind advection schemes in the zero advective velocity case (Fig. 8). This is to be expected, as the only scheme in place is the secondorder diffusion scheme when advection is not present. In the case of   non-zero advection (Fig. 9), the central scheme displays second-order convergence, while the upwind converges to the first order for ≥ 20. Considering numerical simulations were run for > 0, as expected, divergent behaviour is observed almost immediately for Nd ≲ 0.5 regardless of the advection scheme. Consequently, only simulations with Nd > 50 are reported in the rest of this article. Fig. 10 shows the outcome of simulations where the central advection scheme is used. In the cases with = 0.1 m∕s, the computational time was between 20 and 25 CPUs. Oscillations appear at Pe ≃ 100 ( Fig. 10(a)), and eliminate the model's predictive power for higher values of Pe (Fig. 10(b)). A comparison between Figs. 7(b) and 10(c) shows that the error associated to the simulations is larger than the one experienced in the Lattice-Boltzmann case.
The simulations performed with the upwind scheme are reported in Fig. 11. In the cases with = 0.1 m∕s, the computational time was between 20 and 25 CPUs. Contrary to the upwind case (Fig. 10) and the Lattice-Boltzmann model (Fig. 7), the oscillations are suppressed and the instabilities removed-crucially, the model's predictive power is maintained also in the case Pe → ∞. As predicted, the numerical diffusivity is relevant: a comparison between Fig. 11(c), Figs. 7(b) and 10(c) shows that the error associated to the upwind finite-difference simulations up to Pe = 1 000 is comparable (albeit slightly smaller) to the central finite-difference, and larger than the one experienced in the Lattice-Boltzmann case.

Sharp start: finite-difference. negative-diffusivity correction
The convergence behaviour of the explicit finite-difference method with upwind advection scheme and negative-diffusivity correction is reported in Fig. 12. The model with the negative-diffusivity correction ( Fig. 12(b)) is shown to display the same order of convergence of the pure upwind ( Fig. 12(a)), but is able to display convergent behaviour at coarser meshes.
In Fig. 13, the action of the variation of the negative-diffusivity parameter over the results of the simulations is displayed. In all cases, the computational time was between 25 and 35 CPUs with peaks up to 67 CPUs. The results show that the numerical diffusivity cannot be D. Dapelo et al.  removed; however, it is shown that the introduction of the negativediffusivity term allows a significant reduction of the error of up to 50%. Fig. 13(f) shows that, for Pe ⪆ 1 000, values of larger than 0.8 result in the appearance of oscillations; it is therefore necessary to strike a balance between precision (viz., suppression of the numerical diffusion) and stability (that is, preventing oscillations). This balance is achieved by minimizing the value of the error as a function of and is achieved for = 0.8145 ± 0.0025 for Pe = 10, = 0.890 ± 0.001 for Pe = 100 and = 0.769 ± 0.001 for Pe = 1000.

Upscaling tests for the non-periodical case
In Fig. 14   processor per second). The speed-up of a given lattice at a number of processors 1 with respect to the initial number of processors 0 is given as: That gives, for instance, the value = 4.012 for 120 processors against a theoretical speed-up of th ≡ 1 ∕ 0 = 6 in the finite-difference model for the = 317 3 lattice. It can be seen that both the model display similar upscaling characteristics. This proves that the models behave similarly in terms of parallelizability.

Discussion
The Lattice-Boltzmann advection-diffusion simulations with smooth initial conditions (Section 6.1) display the expected second-order convergence. The Lattice-Boltzmann advection-diffusion simulations with non-smooth initial conditions (Section 6.3), are found to retain secondorder convergence if they respect relaxed versions of the stability conditions proposed in the literature; otherwise, they demonstrate a deterioration of convergence from second to first order. Specifically, the Pe g criterion (Eq. (27)) is found to be determinant, as transition from second to first order of convergence is observed to occur for Pe g ≳ 100 and in the case Pe g = 0 (which makes it reasonable to suppose that the same happens in a sufficiently small neighbour of Pe g = 0). Conversely, the stability criterion based on the simplified version of [20]'s observation (Eq. (28)) is found to be either irrelevant or a duplication of the Pe g criterion for the numerical setup taken in consideration within this work.
Lattice-Boltzmann and explicit finite-difference models are shown to be equivalent in terms of computational expense. The upwind finitedifference with the negative-diffusivity correction constitutes an exception, but the increase in computational expense is within 50%. Peaks of computational expense should be considered carefully, as they are likely to be generated by sudden increases of resource requirements in the laptop's operating system. Therefore, the comparison between Lattice-Boltzmann and finite-difference models is limited to numerical precision and predictive power.
A comparison of the error over Pe and a conversion of the latter into Pe g for the Lattice-Boltzmann (Fig. 7(b)), the central (Fig. 10(c)) and the upwind finite-difference models (Fig. 11(c)) in the context of sharp initial conditions shows that: (i) for Pe g < 10, the finite-difference model with the central scheme for advection produces the smallest amount of error; (ii) for 10 < Pe g < 100, the Lattice-Boltzmann is superior; and (iii) the finite-difference model with the upwind scheme for advection is never competitive. These statements are justified by the following considerations: (i) the central finite-difference performs better than the other methods for Pe g < 10 because it displays secondorder convergence whereas the upwind is only first-order convergent, and Lattice-Boltzmann's convergence deteriorates to the first order for Pe g → 0; (ii) having the same (viz., second) order of convergence, the central finite-difference performs worse than the Lattice-Boltzmann for 10 < Pe g < 100 because of the onset of oscillations (see Fig. 10(a)), and this is predictable by comparing Eq. (25) to Eq. (38); (iii) the upwind finite-difference model is inferior to the Lattice-Boltzmann because of its numerical diffusivity (see Fig. 11(a)).
A comparison between Figs. 7(a), 10(a) and 11(a) shows that the Lattice-Boltzmann and the central finite-difference models rapidly lose their predictive power when Pe g > 100 because of the oscillations, whereas the upwind finite-difference remains stable up to Pe g → ∞.
Finally, Fig. 13 shows that the introduction of the negative-diffusivity correction brings about a significant reduction of the numerical diffusivity error. However, oscillations are re-introduced if the weight of the correction is too high, i.e.: → 1. Therefore, the optimum value of must be chosen in order to reduce the error as much as possible, whilst maintaining system stability through oscillation removal. The proposed method of numerically finding the minimum of the error in function of is shown to successfully provide this balance. However, the magnitude of the oscillations -and therefore the optimum value of -depends on the values of Pe; this is exemplified by the fact that three different optima were found for three different values of Pe. Consequently, the operation of finding the optimum value of by finding the minimum value of the error should be repeated for each application of the method, and each set of data if possible.

Conclusions
The Lattice-Boltzmann advection-diffusion model is shown to possess advantages in terms of parallelizability comparable to pure Lattice-Boltzmann, and to exhibit second order of convergence for smooth initial conditions. For non-smooth initial conditions, the order of convergence is shown to primarily depend on Pe g . As this work is primarily concerned on assessing the dependence of Lattice-Boltzmann and explicit finite-difference models for advection-diffusion from Pe and Pe g with the explicit case Pe, Pe g → ∞ in mind, the dependence of Lattice-Boltzmann advection-diffusion's order of convergence from Co and * has not been thoroughfully assessed. Further work is required to address this specific gap and will also (i) investigate Lattice-Boltzmann's order of convergence for Pe g → 0, and (ii) find the critical value for Pe g in the smooth-initial condition case.
Based on the considerations reported in Section 7 on the mutual comparison between central finite-difference, upwind finite-difference and Lattice-Boltzmann advection-diffusion models, the following recommendations are given. a. Pe g < 10: use explicit finite-difference with the central differentiation scheme for advection. b. 10 < Pe g < 100: use Lattice-Boltzmann. c. 100 < Pe g < ∞: use explicit finite-difference with the upwind differentiation scheme for advection and the negative-diffusivity correction. Tune the parameter in order to strike the optimum balance between stability and precision (viz., find the minimum of the error as a single-variable function of ).
Further work will also include the application of the explicit finitedifference models to the continuous boundary condition benchmark, and the extension of the comparative work to models implementing multiple-relaxation-time and two-relaxation-time Lattice-Boltzmann for advection-diffusion.