Aurore Loisy
email: aurore.loisy@bristol.ac.uk

Aurore Naso

Peter D M Spelt

The effective diffusivity of ordered and freely evolving bubbly suspensions

We investigate the dispersion of a passive scalar such as the concentration of a chemical species, or temperature, in homogeneous bubbly suspensions, by determining an effective diffusivity tensor. Defining the longitudinal and transverse components of this tensor with respect to the direction of averaged bubble rise velocity in a zero mixture velocity frame of reference, we focus on the convective contribution thereof, this being expected to be dominant in commonly encountered bubbly flows. We first extend the theory of Koch et al. (1989) (which is for dispersion in fixed beds of solid particles under Stokes flow) to account for weak inertial effects in the case of ordered suspensions. In the limits of low and of high Péclet number, including inertial effect of the flow does not affect the scaling of the effective diffusivity with respect to the Péclet number. These results are confirmed by direct numerical simulations performed in different flow regimes, for spherical or very deformed bubbles and from vanishingly small to moderate values of the Reynolds number. Scalar transport in arrays of freely rising bubbles is considered by us subsequently, using numerical simulations. In this case, the dispersion is found to be convectively enhanced at low Péclet number, like in ordered arrays. At high Péclet number, the Taylor dispersion scaling obtained for ordered configurations is replaced by the one characterizing a purely mechanical dispersion, like in random media, even if the level of disorder is very low.

Introduction

Bubble columns are commonly used in a broad range of technologies, notably in the chemical and biochemical industry. Simple bubble columns do not require active stirring and can therefore operate without interior moving parts. The large surface area between gas and liquid is useful for mass and species transfer, possibly involving chemical reactions (e.g., [START_REF] Deckwer | Bubble column reactors[END_REF]) such as in air-lift bioreactors, an example of which is in the treatment of wastewater. Bubble columns are also used for this reason in direct contact heat transfer (e.g., [START_REF] Hewitt | Process heat transfer[END_REF]. Besides offering a large surface area, rising bubbles agitate the liquid flow, which results in enhanced mixing that usually is desired, but this also poses a modelling difficulty. Similar mixing arises also in the diffusion through porous media in the presence of flow, which has been well studied previously, but mostly for fixed beds of particulates, often under creeping flow (e.g., [START_REF] Batchelor | Transport properties of two-phase materials with random structure[END_REF][START_REF] Koch | Dispersion in fixed beds[END_REF]. Mixing in bubble columns is complicated further by the fact that liquid velocity fluctuations are coupled with the dynamics of (deformable) bubbles, usually beyond creeping flow (e.g., [START_REF] Alméras | Mixing by bubble-induced turbulence[END_REF].

In the present study, we consider transport of a scalar (such as the concentration of a chemical species, or the temperature) through incompressible bubbly flows. Gradients of temperature and concentration may, in general, induce fluid motion and influence the velocity field through changes in density and viscosity, or through the interface rheology. If these effects are small, as assumed herein, temperature and solute concentration can be considered as passive scalars. Although the arbitrary choice was made in this study to use the terminology of the mass transfer problem, the results carry over to thermal applications (upon assuming that effects of viscous heating can be ignored).

Our present main interest is the formulation and closure of conservation equations and constitutive relations governing the dispersion of such a scalar in a bubbly suspension over scales (termed hereinafter the "macroscale") that are much larger than the bubble size (termed hereinafter the "microscale"). Under the assumption of macroscale homogeneity and stationarity, scalar dispersion in multiphase systems can be described by a macroscale version of Fick's (or Fourier's) law which relates the macroscale scalar flux to the macroscale scalar gradient through an effective diffusivity tensor (or effective conductivity tensor in thermal applications) [START_REF] Batchelor | Transport properties of two-phase materials with random structure[END_REF][START_REF] Koch | Dispersion in fixed beds[END_REF][START_REF] Koch | The symmetry properties of the effective diffusivity tensor in anisotropic porous media[END_REF]. This effective diffusivity is defined from an Eulerian perspective. Experimentally, scalar dispersion is usually investigated from a Lagrangian point of view. In the Lagrangian framework, the effective diffusivity is defined as the long-time limit of the time rate of change of a fluid tracer's mean-square displacement, that is, as a measure of spread about the mean position. [START_REF] Koch | The symmetry properties of the effective diffusivity tensor in anisotropic porous media[END_REF] demonstrated that the Lagrangian effective diffusivity is equivalent to the symmetric part of the Eulerian effective diffusivity, and that the antisymmetric part of the Eulerian effective diffusivity is associated with anisotropic microstructures.

Scalar dispersion in a suspension of particulates (bubbles, drops, or rigid particles) results from two processes of very different nature: the diffusion by Brownian motion of the molecules, and the convection by the fluid velocity disturbances induced by the particulate motion. The relative importance of these two processes is measured by the Péclet number Pe = U d b /D, where U is the characteristic velocity of the particulates relative to that of the system (defined in section 2), d b is the characteristic size of the particulates, and D is the diffusivity of the bulk. In the limit Pe = 0, the effective diffusivity is purely diffusive and depends only on the particulate-to-bulk diffusivity ratio, possible discontinuity of the scalar at the interface, particulate volume fraction, and suspension microstructure (i.e., the positions, shapes, and orientations of the inclusions). This particular situation is essentially relevant to heat and electricity conduction in composite materials. When Pe 1, the dominant contribution to the effective diffusivity is due to convective mixing. This last regime is that generally encountered in bubbly flows.

Recently [START_REF] Alméras | Mixing by bubble-induced turbulence[END_REF] investigated experimentally the dispersion of a lowdiffusive dye within a homogeneous swarm of high-Reynolds-number rising bubbles at Pe = O(10 6 ); herein we define the Reynolds number as Re = U d b /ν c where ν c is the kinematic viscosity of the liquid. They showed that scalar mixing primarily results from pseudo-turbulence, i.e., from the liquid agitation produced by bubble wake interactions, and can be modeled in a manner analogous to dispersion in shear-induced turbulence [START_REF] Taylor | Diffusion by continuous movements[END_REF]. Apart from the work of [START_REF] Alméras | Mixing by bubble-induced turbulence[END_REF], the only other experimental investigation of mixing in homogeneous bubbly flows reported in the literature is the preliminary study of [START_REF] Mareuge | Bubble induced dispersion of a passive scalar in bubbly flows[END_REF] which consists in a single data point. To the best of our knowledge, neither theoretical nor numerical investigations of scalar mixing in homogeneous bubbly flows have been reported thus far. Theoretical work is, however, available for other types of multiphase systems, and we shall review these now.

The determination of such an effective diffusivity, at the macroscale, necessitates consideration of the conditions at the microscale. One class of analytical work is devoted to the study of dilute systems with fixed random microstructure, for instance, as a model of a porous medium. In the absence of convection (Pe = 0), the analytical expression of the effective diffusivity is available in the dilute limit from analysis of the corresponding problem in conduction of heat or electricity through a dispersed medium (e.g., [START_REF] Maxwell | A treatise on electricity and magnetism[END_REF], [START_REF] Jeffrey | Conduction through a random suspension of spheres[END_REF]). The problem of scalar dispersion in the presence of a bulk convective motion (Pe > 0) has been analyzed by [START_REF] Koch | Dispersion in fixed beds[END_REF] for Stokes flow through a random bed of fixed solid spheres. Using the method of conditional averaging pursued earlier by [START_REF] Hinch | An averaged-equation approach to particle interactions in a fluid suspension[END_REF], they carried out an asymptotic analysis in low volume fraction of the effective diffusivity for all values of the Péclet number. Three mechanisms causing dispersion at high Péclet number were identified: mechanical dispersion resulting from the stochastic velocity field in the bulk, which is independent of Brownian diffusion and grows as U d b , holdup dispersion in stagnant and recirculating regions which is proportional to U 2 d 2 b /D, and boundary-layer dispersion which grows as U d b ln(U d b /D) near the solid particle surfaces.

Another class of analytical studies assumes a periodic microstructure. For the pure diffusion problem (Pe = 0), analytical solutions have been derived for a composite material consisting of regularly arranged spheres embedded in a homogeneous matrix [START_REF] Rayleigh | LVI. On the influence of obstacles arranged in rectangular order upon the properties of a medium[END_REF][START_REF] Sangani | The effective conductivity of a periodic array of spheres[END_REF], and the effect of anisotropy has been investigated by considering periodic arrangements of spheroidal inclusions [START_REF] Kushch | Conductivity of a periodic particle composite with transversely isotropic phases[END_REF][START_REF] Harfield | Conductivity calculation for a two-phase composite with spheroidal inclusions[END_REF]). In the presence of convection (Pe > 0), the general theory of dispersion developed by [START_REF] Brenner | Dispersion resulting from flow through spatially periodic porous media[END_REF] and [START_REF] Brenner | Dispersion resulting from flow through spatially periodic porous media II. Surface and intraparticle transport[END_REF] provides a consistent framework for determining the effective diffusivity in spatially periodic media. [START_REF] Koch | The effect of order on dispersion in porous media[END_REF] carried out explicit calculations for a periodic porous medium consisting of fixed solid particles arranged in a cubic lattice and embedded in a continuous phase under Stokes flow conditions. They showed that in ordered systems, the mechanical dispersion encountered in random media is absent, and that at high Péclet number, either Taylor dispersion, growing as U 2 d 2 b /D, or enhanced diffusion, which is proportional to D, is obtained depending on the direction of the mean flow relative to the lattice structure.

In bubbly flows, the spatial arrangement of the inclusions evolves in time, the microstructure of the suspension is unknown a priori, and Stokes flow is usually not applicable. For these reasons, prior analyses are, a priori, not applicable to bubbly suspensions. Nevertheless, we showed in prior work [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]) that the dynamics of freely evolving bubbly suspensions at moderate Reynolds number shares some common features with that of ordered arrays of bubbles. It is therefore of fundamental interest to investigate, contrast and compare the mixing properties of ordered and freely evolving bubbly suspensions in light of prior asymptotic analyses for ordered and random arrangements of rigid particles.

In this paper we investigate scalar dispersion, by determining the effective diffusivity in ordered and freely evolving bubbly suspensions, specifically, the contribution of bubbleinduced velocity disturbances thereof. The prior work outlined above has established that in the systems studied therein, the effective diffusivity can be much larger than that in each of the fluids involved, even if the diffusivity in the two media is the same and the scalar is continuous at the surface of particulates. In view of the already significant number of parameters involved, we shall therefore adopt this restriction here. Such a simplified approach will not provide an accurate description of real bubbly flows, but should shed some light on the fundamental mechanisms of mixing in these systems.

The paper is organised as follows. The theoretical framework and problem statement are provided in section 2. Our numerical approach to compute the effective diffusivity is presented, and followed by a description of the regimes and the range of parameter values that are investigated herein, in section 3. The first objective (in section 4) is to elucidate the role played by liquid inertia in ordered suspensions, using direct numerical simulation and analysis. The second objective (in section 5) is to investigate the effective diffusivity of freely evolving suspensions for a wide range of Péclet numbers, to compare it with that obtained for ordered suspensions, and to evaluate the effect of introducing additional degrees of freedom in the system. Finally, the main results and perspectives of this work are provided in section 6.

Problem statement

The local evolution of the passive scalar c in each fluid is governed by

∂ c ∂t + ∇ • q = 0 (2.1a)
where q is the flux of scalar given by

q = uc -D∇c (2.1b)
with u the fluid velocity and D the constant scalar diffusivity. We assume that the scalar and its gradient are continuous across the interface, and phase change is not considered in this study. Under these assumptions, no distinction between the phases is needed for the scalar transport, which is described by (2.1) in the entire system. We return to these restrictions in section 2.1 and in the Conclusions section; the objective here is to study this key basic reference problem.

In the context of heat transfer, (2.1) derives from the energy balance upon neglecting viscous heating, in this case c would represent the temperature, continuous at the interface, and D the thermal diffusivity as defined by Fourier's law, assumed to be equal in both gas and liquid. In the context of mass transfer, (2.1) describes the transport of a chemical species present at very low concentration c so that Fick's law describes the conservation of mass, neglecting any difference in molecular diffusivity D and solubility of the species in the two phases. While the assumption of equal molecular diffusivities is never satisfied in real systems, the assumption of a unit dimensionless Henry's constant is reasonably applicable to, e.g., carbon dioxide dispersion in the air-water system.

The fluid motion in the gas and liquid is governed by the incompressible Navier-Stokes equations, which are coupled at the interface by the appropriate jump conditions, namely the continuity of velocity and of tangential traction across the interface, and a jump in normal traction due to surface tension.

Macroscale description

The problem we are concerned with here is the modeling of scalar transport at a macroscale, that is, at the scale whereat the suspension may be seen as a homogeneous continuum, without distinction between the two phases. In order to obtain such a macroscopic description, we consider an ensemble of realizations of the suspension, these realizations having the same macroscopic conditions (e.g., fluid properties, gas volume fraction) but different microscopic configurations (e.g., bubble individual positions, shapes and velocities), and average over those realizations. In concrete terms, ensemble averaging would be realized by averaging over a large number of experiments run under identical macroscopic conditions. The ensemble-averaged transport equation is obtained from ensemble averaging the local transport equation (2.1). It reads

∂ c ∂t + ∇ • q = 0, (2.2)
where denotes the ensemble average operator, and where the ensemble-averaged flux is given by

q = u c -D∇ c + u c (2.3)
where the velocity fluctuations are defined by u = uu and the scalar fluctuations by c = c -c . Under the restrictions set out above, the average flux consists of three contributions: (i) u c is the advection of the average scalar field at the average system velocity; (ii) -D∇ c is the diffusion of the average scalar field directly by the average scalar gradient; (iii) u c corresponds to the advection of the scalar fluctuations by the velocity fluctuations induced by bubble motion. When the suspension is statistically homogeneous and in a statistically stationary state, the linearity in c of the local flux (2.1b) results, in the presence of an imposed constant average scalar gradient, in a macroscale constitutive relation of the form [START_REF] Koch | Dispersion in fixed beds[END_REF][START_REF] Koch | The symmetry properties of the effective diffusivity tensor in anisotropic porous media[END_REF]:

q = u c -D eff • ∇ c (2.4)
where D eff is a constant effective diffusivity tensor. Comparison of the effective diffusivity definition (2.4) with the average flux expression (2.3) yields the expression of the effective diffusivity. In order to reflect the contributions to the scalar flux identified above, it is customary to write the effective diffusivity as

D eff = DI + D conv (2.5)
where

D conv • ∇ c = -u c (2.6)
is the convective contribution arising from bubble-induced velocity fluctuations. For this model to be complete, one must find a closure relation for D conv only in terms of macroscopic quantities appearing in the problem statement. We recall here that further contributions to the average flux (2.3) and hence to the effective diffusivity (2.5) arise if the diffusivity in the fluids are not the same, or if the concentration is discontinuous at fluid/fluid interfaces (e.g., [START_REF] Batchelor | Thermal or electrical conduction through a granular material[END_REF]; [START_REF] Koch | Dispersion in fixed beds[END_REF]). We return to the significance of this in section 6 below.

Effective transport properties

To determine the effective diffusivity for (unbounded) homogeneous bubbly suspensions, we represent such flows by the periodic repetition of a cubic unit cell containing a finite number N b of freely moving bubbles of equal volume, building on our prior work on the dynamics of bubbles for this model system [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]). In the limit N b = 1, one obtains a simple cubic array of bubbles, which is of interest as a model of perfectly ordered suspensions. The opposite limit of large N b is of interest as a model of real suspensions, although convergence with the number of bubbles would have to be verified. We shall refer hereinafter to this setup with one bubble in the cell as an ordered array, and to that with more than one bubble in the unit cell as a free array.

The bubbles rise under the sole effect of buoyancy. Herein, an upward-pointing primary axis e 3 of the periodic arrangement is taken to be aligned with gravity (with the exception of the more general analysis presented in section 4.1). From symmetry arguments, and adopting a Cartesian coordinate system,

D conv =   D conv ⊥ D conv 12 D conv 13 D conv 12 D conv ⊥ D conv 13 D conv 31 D conv 31 D conv   (2.7)
where we have introduced the longitudinal and transverse components of the convective contribution to the effective diffusivity, denoted D conv and D conv ⊥ , respectively, and defined by

D conv = D conv 33 and D conv ⊥ = D conv 11 = D conv 22 .
(2.8)

Our first goal is to characterize the effects of liquid inertia (through Re) on the dependence of D conv on Pe for ordered suspensions (N b = 1), thereby extending prior work on dilute ordered arrays of rigid spheres in Stokes flow conditions [START_REF] Koch | The effect of order on dispersion in porous media[END_REF]. Our second goal is to evaluate the effect of introducing additional degrees of freedom in the system (through increasing N b ), and to investigate the dependence of D conv on Pe in freely evolving suspensions (sufficiently large N b ). As we found the offdiagonal components to be zero in all configurations that we investigated, only results for the longitudinal and the transverse components of D conv will be presented.

In dimensionless groups, we shall use as characteristic length scale the bubble size d b , which is defined, since bubbles are deformable, as the (equivalent) diameter of a sphere of the same volume. The characteristic velocity U is taken here as the bubble rise velocity in the frame of the suspension (the so-called drift velocity U = u d -u , where the first term is the volume average of velocity on the disperse phase only and the second one is the same average in the entire system). As already mentioned, a key dimensionless group appearing in the scalar transport problem is the Péclet number Pe = U d b /D which compares advective and diffusive transport. Our main objective is to elucidate the effect of the value of Pe on the effective diffusivity using analytical and numerical methods.

The effective diffusivity necessarily also depends on the gas volume fraction φ = (N b πd 3 b )/(6h 3 ) (h is the linear size of the unit cell); the analytical and computational methods used here pose some restrictions on the range of φ values that can be studied herein, we postpone discussion of that to the pertinent sections below. We also consider the effects of the number of bubbles in the periodic cell, N b , which affects the order in the suspension: N b = 1 corresponds to a cubic array, whereas more bubbles results in a different microstructure (the latter term encompasses all the information about the statistical distribution of the bubble positions, shapes, orientations, etc.). Since scalar transport is coupled to momentum transport, the bubble Reynolds number Re = U d b /ν c may also play a significant role that will be investigated here as well. The ranges of φ, N b and Re studied here are summarized in table 1.

As the bubbly flows we consider are buoyancy-driven, a difficulty arises from the fact that U is a priori unknown, and depends in a complex manner on N b , φ, the density and viscosity ratios between both phases, the Archimedes (or Galileo) number Ar = ρ c |ρ d -ρ c |gd 3 b /µ c , and the Bond (or Eötvös) number Bo = |ρ d -ρ c |gd 2 b /γ, where the subscripts d and c refer to the disperse (gas) and continuous (liquid) phases, respectively, g is the magnitude of the gravitational acceleration, ρ denotes density, µ is the dynamic viscosity, and γ is the surface tension. In most bubbly flows of practical relevance, the gas-to-liquid density and viscosity ratios are vanishingly small. Their precise values are not important from a physical point of view as long as they are small enough; in the simulations, the gas-to-liquid density and viscosity ratios were set to ρ d /ρ c = 10 -3 and µ d /µ c = 10 -2 , respectively. The dependence of U on (Ar , Bo, φ, N b ) has been addressed in [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF] and is not further discussed here. In the present study, we shall therefore assume that U is known.

Methodology

For convenience of numerical implementation, we reorganise the problem formulation by introducing the decomposition

c = c + c (3.1)
where c is the imposed constant linear scalar field

c = ∇ c • x. (3.2)
The advantage of this decomposition is that the disturbance field c is then spatially periodic. The governing equation for this disturbance field is

∂ c ∂t + ∇ • (uc) -∇ • (D∇c) = -u • ∇ c (3.3)
which is the equation we integrate numerically. The convective contribution to the effective diffusivity is then calculated from

D conv • ∇ c = -u c (3.4)
which can be shown to be equivalent to (2.6). In this expression, is defined as an ensemble average operator, as above. For statistically homogeneous and stationary systems, as considered here, it is inferred from ergodicity that ensemble averaging is identical to volume and time averaging. As a consequence, D conv is computed from (3.4) with the ensemble average being replaced in practice by a volume average combined with a time average over an appropriate time period.

Numerical method

Thus, the components of D conv are obtained from direct numerical simulations (DNS) by imposing a constant linear scalar field c and determining the resulting periodic disturbance scalar field. Two distinct simulations are required to fully determine the five independent components of D conv : in one simulation, ∇c = e 3 , which yields D conv 13 and D conv , in the other simulation, ∇c = e 1 , which yields D conv ⊥ , D conv 12 , and D conv 31 . The off-diagonal components of D conv were found to be zero (up to computer accuracy for ordered arrays, and statistical uncertainty for free arrays) for all the sets of parameters we considered, and therefore will not be shown.

The numerical methods employed to solve the two-phase flow have been described in detail in [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]. In short, we employ a standard projection method [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF]) to integrate the incompressible Navier-Stokes equations, a level-set method (e.g., [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF])) to capture the moving gas-liquid interface, and surface tension is accounted for using the continuum surface force model [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF].

Our algorithm proceeds iteratively through the following steps: (i) The position of the interface is first advanced in time according to the modified level-set method of [START_REF] Sabelnikov | Modified level set equation and its numerical assessment[END_REF] using a third-order total-variationdiminishing (TVD) Runge-Kutta scheme. The level-set function is then reinitialized using the procedure of [START_REF] Russo | A remark on computing distance functions[END_REF], and a correction is finally applied to enforce volume conservation.

(ii) The scalar transport equation (3.3) is advanced by using a mixed Crank-Nicolson/third-order Adams-Bashforth time-stepping scheme.

(iii) The time integration of the incompressible Navier-Stokes equations is then carried out using a mixed Crank-Nicolson/third-order Adams-Bashforth scheme and consists in the combination of a predictor step, where a temporary velocity field is estimated by ignoring the effect of pressure, and of a corrector step, where the velocity field is corrected by the pressure gradient term computed from the divergence-free condition. Spatial discretization relies on a mixed finite difference/finite volume approach on a fixed, staggered, Cartesian grid. Second-order centered schemes are generally employed, except for advective terms which are discretized using fifth-order weighted-essentiallynonoscillatory (WENO) schemes.

Results of numerical tests are presented in the Appendix.

Parametric study

Four different flow regimes, as defined by the set (Ar , Bo), are considered here. These are described in table 1, and have been studied in [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF] (the same case code names are used). In case S0, the bubbles are spherical and the Reynolds number is vanishingly small, which approaches Stokes flow conditions. In case S1, the bubbles are (nearly) spherical and Re 1. In case C, the bubbles are skirted, and Re ≈ 8. In case E1, the bubbles are ellipsoidal, and Re ≈ 30 -40.

Ordered arrays of bubbles in these four flow regimes have been considered for the smallest volume fraction numerically accessible (value provided in table 1). After a transient regime, all ordered suspensions considered here are in a strictly steady state (for the flow and the scalar) during which the results presented in section 4 were obtained. Simulations of scalar transport in free arrays have been performed for 2 N b 12 in case E1 at φ = 2.4 %. In these conditions, coalescence is indeed absent (it does occur at larger φ), whereas simulations at lower φ for free arrays are excessively expensive for the method and facilities used. In this regime, the system is in an unsteady but statistically stationary state (for the flow and the scalar), during which the statistics presented in section 5 have been measured. For each of these configurations (Ar , Bo, φ, N b ), the drift velocity (and thereby the Reynolds number) is known from [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]. This allowed us to impose the Péclet number a priori.

The numerical simulation results for ordered arrays are compared with the results of analysis at small (but possibly finite) Reynolds number and small volume fraction.

Ordered suspensions

We examine in this section the dispersion of a passive scalar in ordered suspensions of deformable bubbles. Our main objective here is to elucidate the effects of inertia on dispersion, using theoretical analysis and numerical simulation.

Asymptotic analysis

We first determine analytically the convective contribution to the effective diffusivity of ordered suspensions of spherical fluid particulates (bubbles or drops). The Reynolds number of the particulates is assumed to be small so that the Navier-Stokes equations can be approximated by the Oseen equations.

General solution

An ordered array of particulates translating at a drift velocity U is equivalent to an ordered array of fixed particulates immersed in a viscous fluid moving with an average system velocity u = -U . The centers of the particulates are located on the nodes of a simple cubic lattice:

r n = h (n 1 e 1 + n 2 e 2 + n 3 e 3 ) n 1 , n 2 , n 3 = 0, ±1, ±2, . . . (4.1)
where h is the lattice spacing and e i are the unit vectors aligned with the primitive axes of the cubic lattice. In the dilute limit (d b /h 1), the action of these particulates on the fluid can be represented by point forces -f . The convective contribution to the effective diffusivity arising from the far field has been derived by [START_REF] Koch | The effect of order on dispersion in porous media[END_REF] for an ordered array of rigid spheres in the Stokes flow regime. In what follows we extend their result to the case of spherical fluid particulates at small but finite Re.

When Pe 1, the convective contribution to the effective diffusivity arising from the far field can be approximated by [START_REF] Koch | The effect of order on dispersion in porous media[END_REF]):

D conv D = k =0 k 2 û (k)û (-k) (2π) 2 k 4 D 2 + (U • k) 2 , (4.2)
where the summation is over all vectors k in the reciprocal lattice

k = 1 h (n 1 e 1 + n 2 e 2 + n 3 e 3 ) (4.3)
and where û is the three-dimensional Fourier transform of the velocity disturbance u = uu . In Oseen flow past an ordered array of point particulates, û is given by

û (k) = f • (kk/k 2 -I) (2πk) 2 h 3 µ c + i2πh 3 ρ c U • k k = 0, (4.4)
where f is the hydrodynamic force exerted by the ambient fluid on a particulate. In the dilute limit, f can be approximated by the Oseen drag exerted on a single spherical fluid particulate:

f = F f 0,Stokes (4.5)
where f 0,Stokes is the Stokes drag on that particulate [START_REF] Hadamard | Mouvement permanent lent d'une sphere liquide et visqueuse dans un liquide visqueux[END_REF][START_REF] Rybczynski | Uber die fortschreitende Bewegung einer flussigen Kugel in einem zahen Medium[END_REF]):

f 0,Stokes = -2πµ * µ c d b U , with µ * = µ c + 3µ d /2 µ c + µ d , (4.6)
and where F accounts for the finite-Re correction to the Stokes drag [START_REF] Brenner | The resistance to a particle of arbitrary shape in translational motion at small Reynolds numbers[END_REF]:

F = 1 + 1 8 µ * Re. (4.7)
The convective contribution to the effective diffusivity of a dilute ordered array of fluid particulates in Oseen-flow conditions is therefore: where C is the dimensionless tensor: The tensor C only depends on Pe h , Re h , and on the orientation of U relative to the reciprocal lattice (which structure is, for cubic arrays, identical to that of the direct lattice). As highlighted by [START_REF] Koch | The effect of order on dispersion in porous media[END_REF], the asymptotic behavior of C, and hence of D conv , depends on whether there exists any k such that U • k = 0, that is, on whether there exists any separation vector r n in the real space which is perpendicular to U . The asymptotic behavior of D conv , where denotes the tensorial Frobenius norm, is provided in table 2. The results show that the dependence of D conv on Pe in the limits Pe h 1 and Pe h 1 is, qualitatively, not affected by (weak) inertial effects.

D conv D = µ * 2 (2π) 2 d 2 b h 2 F 2 C, (4.8a) regime D conv /(DF 2 d 2 b /h 2 ) Pe h = U h/D Re h = ρcU h/µc if ∃rn | U ⊥ rn if rn | U ⊥ rn Pe h 1 Re h 1 Pe 2 h Pe 2 h Re h 1 Pe 2 h Pe 2 h /Re 2 h Pe h 1 Re h 1 Pe 2 h 1 Re h 1 Pe 2 h 1/Re 2 h Table 2.
C = k * =0 U * • k * k * k * 2 -I 2 k * 2 (2π) 2 k * 4 Pe 2 h + (U * • k * ) 2 1 + Re 2 h (U * • k * ) 2 (2π) 2 k * 4 (4.8b) with U * = U /U , k * = kh, Re h = ρ c U h/µ c ,

Application to ordered arrays rising vertically

Let us now come back to our original problem of an ordered array of particulates rising under the effect of buoyancy. The gravitational acceleration is oriented along a primary axis of the array, g = -ge 3 , and although this is not the only possible solution (see, e.g., [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]), we restrict the analysis to the simplest case of bubbles rising vertically. In this case the hydrodynamic force exerted by the fluid on a particulate is parallel to the drift velocity, and, since this force balances the buoyancy force at steady state, F is related to U through

F = U 0,Stokes U (4.9)
where U 0,Stokes is the terminal velocity of an isolated spherical fluid particulate in Stokes flow:

U 0,Stokes = 1 12 |ρ c -ρ d |gd 2 b µ * µ c , with µ * = µ c + 3µ d /2 µ c + µ d . (4.10)
Note that F can also be expressed in terms of commonly employed dimensionless groups: , and F is given by (4.9)).

F = 1 12µ * Ar 2 Re . ( 4
Note that in (a), D conv is compensated by Pe 2 .

In the "sedimentation" problem considered here, F is generally not known (as U is generally not known): it is a non-trivial function of the flow regime and volume fraction which reduces to (4.7) when φ → 0 and when Oseen-flow approximation is applicable.

The longitudinal and transverse components of the convective contribution, D conv and D conv ⊥ respectively, have been calculated from (4.8) for d b /h = 10 -6 as a function of Pe h for various Re < 1. This very low value of d b /h is required to allow Pe h 1 while satisfying the condition Pe = Pe h d b /h 1 under which the analytical solution has been derived. The results, shown in figure 1, indicate that the asymptotic dependences of D conv and D conv ⊥ on Pe are independent of Re. The sole effect of inertia is to modify the proportionality constants (by a substantial amount for the transverse component though).

In the limit of low Pe h (say, Pe h < 10 1 ), both the transverse and the longitudinal components of D conv exhibit a quadratic dependence on the Péclet number (D conv ⊥, ∝ DPe 2 ). In this regime, diffusion is much faster than convection. As the scalar is advected by velocity disturbances, it rapidly spreads out owing to diffusion, and convective dispersion (measured through D conv ) is influenced by both mechanisms. This regime corresponds to the "convectively enhanced dispersion" regime in [START_REF] Koch | The effect of order on dispersion in porous media[END_REF].

In the limit of high Pe h (say, Pe h > 10 3 ), the transverse component of D conv is independent of the Péclet number (D conv ⊥ ∝ D) whereas its longitudinal component grows quadratically with the Péclet number (D conv ∝ DPe 2 ). In this regime, convection dominates, but owing to the spatial periodicity of the flow, convective dispersion is obtained only if molecular diffusion across streamlines is considered [START_REF] Koch | The effect of order on dispersion in porous media[END_REF]). This regime is termed "Taylor dispersion" owing to the formal analogy, pointed out by [START_REF] Brenner | Dispersion resulting from flow through spatially periodic porous media[END_REF], with one-dimensional shear-induced Taylor dispersion in a capillary tube.

We emphasize that the expression (4.8) has been derived from the approximation (4.2), the validity of which is established only for Pe 1 (which is, in practice, of limited use). Using symmetry arguments, [START_REF] Koch | The effect of order on dispersion in porous media[END_REF] (section 4.2 therein) showed that in the limit Pe 1, Taylor dispersion is obtained if the average flow is perpendicular to a set of planes of both translational and reflectional symmetry, such as Stokes flows parallel to the primary axis of an ordered array of spheres. Taylor dispersion is then easily understood by remarking that, owing to the symmetries of the flow, a fluid tracer particle entering the unit cell at one point, say x, exits the cell at the equivalent point in the next cell, that is, ,⊥ are those suggested by the asymptotic analysis (identical to those used in figure 1), and F is given by (4.9). The lines are drawn to guide the eyes. Note that in (a), D conv is compensated by Pe 2 .

x + he 3 , so that dispersion can only occur if diffusion across streamlines is present [START_REF] Koch | The effect of order on dispersion in porous media[END_REF]). In the presence of inertial effects, the reflectional symmetry is lost, hence this argument does not hold. [START_REF] Koch | The effect of order on dispersion in porous media[END_REF] also demonstrated that, for Stokes flow, the solution for Pe h/d b is identical, at lowest order, to that obtained for Pe 1 (section 4.3 therein, note that their Pe corresponds to Pe h in our notations). Such a demonstration for Oseen flow will not be attempted here. Instead, the range Pe 1 will be explored using direct numerical simulations.

Numerical results

The above analysis provides explicit expressions of D conv and D conv ⊥ . These are valid for spherical bubbles rising at Re < 1 (strictly speaking, at a Reynolds number sufficiently small to assume Oseen flow, in terms of Archimedes and Bond numbers this regime would be reached for Bo < 1 and Ar 1), and in the limits φ → 0 and Pe 1. We shall now determine using numerical simulations whether these restrictions can be relaxed, and if so, to which extent.

We examine the case of suspensions at low (but not vanishing) volume fraction in order to approach the dilute limit assumption, and to focus on the sole effect of inertia. The longitudinal and transverse components of the convective contribution to the effective diffusivity have been computed for h/d b = 6.4, which corresponds to a gas volume fraction of φ = 0.2 % (the smallest volume fraction accessible with the method and facilities used), for each of the four flow regimes listed in table 1, and Pe has been varied from 10 -1 to 10 3 . The results are shown in figure 2 as a function of Pe h , the Péclet number based on the lattice spacing, which is the parameter governing the transition between the two asymptotic limits (see section 4.1 and figure 1). The different colors, symbols and line styles depict the different flow regimes (the lines are drawn to guide the eyes). Qualitatively, figure 2 (D conv ,⊥ ∝ DPe 2 ). At high Péclet number (Pe h 10 3 ), Taylor dispersion is the dominant process, with very efficient mixing in the flow direction (D conv ∝ DPe 2 ) and negligible mixing in the transverse one (D conv ⊥ ∝ D). Inertial effects and bubble deformation only affect the proportionality constants, rather weakly for D conv but substantially for D conv ⊥ . To allow a quantitative comparison between the DNS and the analysis, we present in figure 3

the ratio of D conv ,⊥ to D conv,anal ,⊥
where D conv,anal ,⊥ is given by (4.8) with F computed directly from its definition (4.9). As the range of validity of the analysis is defined in terms of Pe (Pe 1, with Pe the Péclet number based on the bubble diameter), the data are presented here as a function of Pe rather than Pe h . For the longitudinal component, the numerical solution does not deviate by more than 5 % from the theoretical prediction, as can be seen from figure 3(a). The fact that the low-Pe, Oseen-flow analysis yields accurate predictions for D conv at Pe = 10 3 and Re = O(10) is not surprising, as the behavior of D conv /(DF 2 Pe 2 ) is rather insensitive to both the flow regime and the Péclet number (as shown in Fig. 1(a) and 2(a), this quantity does not vary more than 15% for the cases studied). We conclude that, at small volume fraction, D conv can be predicted within ±5 % from (4.8) at any Péclet number up to 10 3 and any Reynolds number up to 40, even when the bubbles are strongly deformed. For the transverse component, the asymptotic analysis underpredicts the value of D conv ⊥ at high Péclet number, even for Re 1. As a consequence, D conv ⊥ cannot be accurately estimated from our analytical solution when the assumptions underlying its derivation are not satisfied. It must be kept in mind though that this component varies much more than the longitudinal one between the regimes of small and large Péclet numbers, and is much more sensitive to the flow regime (Re, shape), which means that its value is more difficult to predict. In all, it is worth stressing that the asymptotic analysis yields the correct qualitative behavior and order of magnitude for D conv ⊥ at least up to Pe = 10 3 and Re ≈ 10, even for strongly deformed bubbles. Finally, we emphasize that we found D conv /D conv ⊥ 10 2 , so the most important component of the effective diffusivity tensor is the longitudinal one, except in situations where there is no longitudinal component of the gradient of the scalar on the macroscale.

To illustrate the dispersion regimes at low and high Péclet number, we present in figure 4 D conv and D conv ⊥ , respectively. In each of these figures, the field of c is represented for each flow regime in a vertical symmetry plane passing through the center of a bubble for Pe = 10 -1 (left) and Pe = 10 3 (right), and the Reynolds number increases from top to bottom. The field of c associated with D conv , shown in figure 4, exhibits similar features at low and high Pe. In contrast, the field of c associated with D conv ⊥ , represented in figure 5, is qualitatively different in these two limits. This illustrates qualitatively why the regimes at low and high Pe are similar for D conv (D conv ∝ Pe 2 ), whereas the scaling laws identified for D conv ⊥ are different in both limits (see figure 2). In addition, the Reynolds number and the bubble shape affect the fore-and-aft symmetry and the details of c , but not its essential features, which results in quantitative but not qualitative effects on D conv and D conv ⊥ .

Freely evolving suspensions

We examine in this section scalar mixing in freely evolving suspensions as represented by the periodic repetition of a unit cell containing several independent bubbles ("free arrays"). Our objective here is threefold: (i) to investigate the effective diffusivity of freely evolving suspensions at small and high Péclet numbers, (ii) to compare and contrast these results with those obtained in ordered systems, and (iii) to evaluate the effect of the system size (number of bubbles in a unit cell, N b ).

For that purpose, we considered a single flow regime (ellipsoidal bubbles at Re = O(10), corresponding to case E1 in table 1) at intermediate volume fraction (φ = 2.4 %) and explored the effect of varying the number of free bubbles N b on the dependence of D conv on the Péclet number. Due to the multiplicity of simulations involved and to their duration (typically several months on 64 cores), only a few different values of N b belonging to a rather limited range have been considered (namely N b = {2, 3, 5, 8, 12} in the simulations for the determination of D conv , and N b = {2, 8} in those for D conv ⊥ ). For the same reason, investigations of the effects of volume fraction and flow regime could not be undertaken.

The longitudinal and transverse components of D conv are plotted in figure 6 as a function of the Péclet number for various values of N b . Note that a very wide range of Péclet numbers is considered. Convergence of D conv with the system size is very fast: the values of D conv are essentially independent of the number of free bubbles for 2 N b 12 at all Péclet numbers. This suggests that D conv is independent of the system size N b , although this would need to be confirmed by considering larger values of this parameter. Our data for D conv ⊥ suggest that convergence with N b is slower for this quantity, especially at high Péclet number, although conclusions can hardly be drawn on this point due to the few values of N b considered.

We first examine the dependence of D conv on the Péclet number in free arrays of bubbles (N b > 1). At small Pe, D conv ,⊥ ∝ DPe 2 , whereas at high Pe, D conv ,⊥ ∝ DPe = U d b . Note that the scaling at high Pe is expected from a simple dimensional analysis in a convection-dominated regime where diffusion plays no role. This regime corresponds to the "mechanical dispersion" regime in [START_REF] Koch | Dispersion in fixed beds[END_REF]. The different dispersion regimes at low and high Pe can also be identified from the features of the scalar fluctuation field c . Instantaneous snapshots of c associated with D conv and D conv ⊥ are shown in figure 7 and figure 8, respectively, for an array of 8 free bubbles at Pe = 10 -1 (left) and at Pe = 10 6 (right). For a given component, the isocontours of c follow markedly different patterns at low and high Pe.

We now compare these results with those obtained for ordered arrays (black crosses in figure 6) and discuss the effect of the microstructure. At small Pe, D conv and D conv quadratically with Pe in both free and ordered arrays. This scaling was also obtained by [START_REF] Koch | Dispersion in fixed beds[END_REF] for low-Pe dispersion in porous media with random microstructure (albeit in the Stokes flow limit). Since in the low-Pe regime, diffusion by the random motion of molecules is much faster than convection by the flow, the microstructure has only a quantitative incidence on D conv , and dispersion is qualitatively identical in ordered and freely evolving suspensions. Note that similar features in the spatial distribution of c can be identified in ordered and free arrays at low Pe (see tubular structures in the left side of figures 4 and 7 for D conv , and quadrupolar ones in the left side of figures 5 and 8 for D conv ⊥ ). We however emphasize that precise quantitative agreement between the results for one and for many bubbles at low Péclet number in figure 6 is not expected, as the flows and the microstructures in the two systems are different [START_REF] Bunner | Dynamics of homogeneous bubbly flows Part 1. Rise velocity and microstructure of the bubbles[END_REF][START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]. At high Pe, the Taylor dispersion scaling obtained in ordered arrays is replaced, in both directions, by a scaling similar to the one characterizing mechanical dispersion, as soon as the relative motion between bubbles is allowed. In this regime, the transverse dispersion is indeed governed by mechanical dispersion. Irrespective of the value of Pe, any Taylor dispersion in the vertical direction is limited by transverse diffusion or dispersion, the latter becoming more significant at large Pe. This results in a scaling similar to that of mechanical dispersion in the longitudinal direction as well, such that a distinction between these two mechanisms (pure mechanical dispersion, or Taylor dispersion limited by transverse mechanical one) cannot be made. Incidentally, mechanical dispersion is also obtained at high Pe in random media in Stokes flow conditions [START_REF] Koch | Dispersion in fixed beds[END_REF]. Although the microstructure of the present bubbly suspensions has not been evaluated quantitatively, visual inspection and prior results on their dynamics [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF] showed that it is not random, but rather characterized by a certain "organization". Despite the fact that freely evolving suspensions resemble ordered ones with respect to their dynamics, scalar dispersion is extremely sensitive to the presence of disorder, and is fundamentally different in perfectly ordered and weakly disordered suspensions at high Péclet number. It does not, however, seem to be sensitive to the degree of disorder, as suggested by the fact that the same scalings with Pe are obtained for random porous media and weakly disordered suspensions. We stress that this last statement is purely speculative, and would require a quantitative study of the effect of the microstructure to be confirmed.

We finally attempt a comparison of our results with the experimental data of [START_REF] Alméras | Mixing by bubble-induced turbulence[END_REF], who measured the effective diffusivity of a homogeneous swarm of high-Reynolds-number rising bubbles at Pe ≈ 1.75×10 6 for gas volume fractions ranging from 1 % to 13 %. It is important to stress that in these experiments, Re ≈ 700, whereas in the simulations, Re ≈ 30, so the comparison is only indicative. Interpolation (by eye) of their data at φ ≈ 2.4 % (figure 10 in their paper) yields D eff /D = 1×10 6 and D eff ⊥ /D = 5×10 5 . These experimental values are represented by purple stars in figure 6. Note that at such high Péclet number, the dominant contribution to D eff is due to D conv , so it seems reasonable to assume that these are equivalent. The order of magnitude of D eff /D is comparable in the experiment and in the simulation, whereas D eff ⊥ /D is much higher in the experiment. This difference can be explained from the different properties of the numerical and experimental flows considered: partition coefficient (the dye concentration in the gas is presumably zero in the experiments from [START_REF] Alméras | Mixing by bubble-induced turbulence[END_REF]), diffusivity ratio, and bubble-induced liquid agitation in the horizontal direction. In our simulations of free arrays at moderate Re, the bubbles were indeed observed to rise along nearly straight vertical lines, and the anisotropy ratio characterizing the liquid velocity variance, 2 u 3 u 3 / u 1 u 1 + u 2 u 2 , is approximately 8 (for N b = 8), whereas in the experiment at high Re, the bubble motion is fully three-dimensional, and the anisotropy ratio is approximately 2. Finally, as only one value of the Péclet number was considered in the experiments of [START_REF] Alméras | Mixing by bubble-induced turbulence[END_REF], no comparison of their data with our results can be offered regarding the dependence of the effective diffusivity on the Péclet number.

Conclusions

In this study we investigated scalar dispersion in homogeneous bubbly suspensions as described by an effective diffusivity tensor. The longitudinal and transverse components of the convective contribution to the effective diffusivity, denoted D conv and D conv ⊥ , respectively, have been computed for bubbly suspensions in various flow regimes. This convective contribution is that associated with bubble-induced agitation, and is the dominant contribution to the effective diffusivity in commonly encountered bubbly flows.

The dispersion theory of [START_REF] Koch | The effect of order on dispersion in porous media[END_REF] indicates that convective mixing mechanisms in ordered suspensions in Stokes-flow conditions differ at low and high Péclet numbers. According to this theory, when the bulk flow is aligned with a primary axis of a simple cubic lattice of spheres, convectively enhanced dispersion is expected at low Péclet number, whereas Taylor dispersion should dominate at high Péclet number. In the present study, we have extended this theory to account for weak inertial effects, and we have shown that these two dispersion regimes are qualitatively unchanged in the presence of (weak) inertia. This result has been confirmed by direct numerical simulations for values of the Reynolds number ranging from vanishingly small to moderate. In all investigated cases, D conv was found to be significantly larger than D conv ⊥ , and theoretical predictions have been shown to yield the correct qualitative behaviour and order of magnitude of both D conv and D conv ⊥ in a variety of flow regimes (spherical to strongly deformed bubbles with Reynolds numbers up to 10) at small volume fraction.

Direct numerical simulations of scalar transport in freely evolving bubbly suspensions, as represented by free arrays of bubbles, have been carried out for a wide range of Péclet numbers, and the effect of introducing additional degrees of freedom in the system has been evaluated. At low Péclet number, dispersion in free arrays is convectively enhanced, as in ordered ones. At high Péclet number, in freely evolving suspensions wherein at least two bubbles are present in a unit cell, the longitudinal component of the effective diusivity exhibits a scaling that is similar to that characterizing mechanical dispersion. This suggests that the limiting role of molecular diffusion to Taylor dispersion is taken over by mechanical dispersion, or that mechanical dispersion itself dominates. Besides, the effective diffusivity seems to be weakly sensitive to the number of bubbles present in a unit cell. This last assertion requires more thorough investigations to be confirmed, but is encouraging regarding the possibility of computing the effective diffusivity of homogeneous bubbly flows from direct numerical simulations of systems of relatively small size. This would allow in particular a thorough investigation of the roles played by the volume fraction and the flow regime, which could not be undertaken as part of the present study.

The results presented in this paper are restricted to bubbles having the same diffusivity as that of the surrounding liquid, and to scalar fields that are continuous across the interface, and therefore cannot be straightforwardly compared to those obtained in real bubbly flows. A jump in the scalar field, which represents the difference in solubilities given by Henry's law in the context of chemical species transport, as well as a difference in diffusivities, would introduce a diffusive contribution to the effective diffusivity tensor (2.5) in addition to the convective one considered in this study. The present results show the convective contribution at large Péclet numbers and modest volume fraction to be substantially larger than the diffusive contribution from nonequal diffusivities or solubilities [START_REF] Maxwell | A treatise on electricity and magnetism[END_REF][START_REF] Jeffrey | Conduction through a random suspension of spheres[END_REF][START_REF] Koch | Dispersion in fixed beds[END_REF]. A difference in diffusivities or solubilities would however also have some indirect effect on the convective contribution, which magnitude should be investigated in the future.

Besides the effective diffusivity, another quantity of practical importance is the rate of interfacial scalar transport in the presence of an average scalar gradient between the disperse phase and the bulk. Heat and mass exchanges across phase boundaries are traditionally expressed as dimensionless transfer coefficients called the Nusselt and the Sherwood numbers, respectively. Their functional dependences on suspension properties, in particular the volume fraction, have been the subject of analytical [START_REF] Acrivos | Heat transfer to a slowly moving fluid from a dilute fixed bed of heated spheres[END_REF], numerical [START_REF] Aboulhasanzadeh | Effect of bubble interactions on mass transfer in bubbly flow[END_REF], and experimental [START_REF] Colombet | Experimental study of mass transfer in a dense bubble swarm[END_REF][START_REF] Colombet | Dynamics and mass transfer of rising bubbles in a homogenous swarm at large gas volume fraction[END_REF] studies. Formally, the Nusselt and the Sherwood numbers are closure coefficients for the conditionally averaged scalar transport equation, where the conditional average is defined as an ensemble average over the subset of realizations wherein a particulate is present at a given position. Less formally, the Nusselt and Sherwood numbers are related to a "mesoscale" description of scalar transfer between the two phases, whereas the effective diffusivity is associated with a "macroscale" description of scalar transport through a two-phase mixture seen as a continuum. They correspond to different closure problems, and one cannot be inferred from the other. Nevertheless, the present work will be primarily important for mass transfer processes in bubbly flows that are liquid-phase controlled. This is because then the mixture concentration distribution is key, whereas if it is gas-phase controlled, the concentration in the liquid will be almost uniform and one is primarily concerned by the circumstances inside each bubble. This work benefited from the financial support of the French research agency (grant ANR-12-BS09-0011), and was performed using the HPC resources provided by GENCI-CINES and GENCI-IDRIS (grant x20162b6893), PSMN ( École Normale Supérieure de Lyon), P2CHPD (Université Claude Bernard Lyon 1) and PMCS2I ( École Centrale de Lyon).

Appendix: Spatial convergence tests

We present the results of some spatial convergence tests of the algorithm solving the scalar transport equation. The results of similar tests for the algorithm solving the flow are shown in [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]; [START_REF] Loisy | Direct numerical simulation of bubbly flows: coupling with scalar transport and turbulence[END_REF].

The effect of the grid spacing on D conv and D conv ⊥ has been assessed for case E1 at Pe = 10 3 for one value of the volume fraction (φ = 2.4 %), in both ordered and free configurations. For ordered arrays, three different resolutions were tested, namely d b /∆x = {20, 40, 60} with ∆x the grid spacing. The results are shown in figure 9. The error in the values of D conv and D conv ⊥ arising from spatial discretization is less than 1 % when a resolution of 40 grid cells per bubble diameter is used. This resolution is the same as that used for the simulation of the corresponding bubbly flow in [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]. In practice, we used for each configuration the same resolution as that selected for the simulation of the corresponding ordered bubbly suspensions (see [START_REF] Loisy | Buoyancy-driven bubbly flows: ordered and free rise at small and intermediate volume fraction[END_REF]), namely 60 grid cells per diameter for case C and 40 grid cells per for the other cases.

For free arrays, due to the computational cost of the simulations, only two different resolutions were tested, namely 20 and 30 grid cells per bubble diameter, for an array of 8 bubbles. Simulations at higher resolution were too expensive to be continued over sufficiently long times to allow a quantitative estimate of the uncertainty. However the values of D conv and D conv ⊥ obtained with the finer grid, depicted by filled red squares in figure 6, are nearly indistinguishable from those obtained with the coarser grid. A resolution of 20 grid cells per diameter was therefore concluded to be sufficient for free arrays in view of the present purposes.

For a given case, the same resolution was used for all Péclet numbers. Note that when the gas diffusivity differs from that of the liquid (a situation not considered here but frequently encountered in practice), finer resolutions may be required, as thin scalar boundary layers around the bubbles would then need to be resolved.

  Asymptotic order of D conv depending on Pe h , Re h , and on the orientation of the mean flow relative to the real lattice, based on the solution (4.8), derived for an ordered array of point particulates in Oseen flow conditions (F is the Oseen drag divided by the Stokes drag).

  and Pe h = U h/D. The solution given by[START_REF] Koch | The effect of order on dispersion in porous media[END_REF] (equation (4.5) therein) for rigid spheres and Stokes flow is recovered in the limit Re → 0 and µ d /µ c → ∞.

Figure 1 .

 1 Figure 1. Longitudinal (a) and transverse (b) components of D conv as a function of the Péclet number based on the lattice spacing (Pe h = U h/D) for ordered arrays of point particulates at various small but finite Reynolds numbers (U = U e3, d b /h = 10 -6, and F is given by (4.9)). Note that in (a), D conv is compensated by Pe 2 .

Figure 2 .

 2 Figure 2. Longitudinal (a) and transverse (b) components of D conv as a function of the Péclet number based on the lattice spacing (Pe h = U h/D) for ordered arrays in various flow regimes at small volume fraction (φ = 0.2 %). The normalizations of D conv,⊥ are those suggested by the asymptotic analysis (identical to those used in figure1), and F is given by (4.9). The lines are drawn to guide the eyes. Note that in (a), D conv is compensated by Pe 2 .

Figure 3 .

 3 Figure 3. Numerical solution D conv divided by the analytical solution D conv,anal as a function of the Péclet number based on the bubble diameter (Pe = U d b /D) for ordered arrays in various flow regimes at small volume fraction (φ = 0.2 %): longitudinal (a) and transverse (b) components. D conv,anal is given by (4.8).

Figure 4 .

 4 Figure 4. Scalar fluctuation field c associated with D conv , shown in a vertical symmetry plane passing through the center of a bubble, for ordered arrays in various flow regimes at Pe = 10 -1 (left) and Pe = 10 3 (right). The imposed scalar field c increases linearly within the cell from bottom to top (φ = 0.2 %, the entire cell is shown, and gravity is pointing downward).

Figure 6 .

 6 Figure 6. Longitudinal (a) and transverse (b) components of D conv as a function of the Péclet number for various numbers of free bubbles N b in the unit cell (N b = 1 corresponds to an ordered array). Symbols other than purple stars: DNS (Re ≈ 30, φ = 2.4 %); purple stars: experimental data of Alméras et al. (2015) (Re ≈ 700, φ ≈ 2.4 %). A spatial resolution of d b /∆x = 20 was used for N b > 1, the effect of increasing resolution to d b /∆x = 30 is illustrated by the filled red squares for N b = 8 and Pe ≈ 10 3 (d b is the bubble volume-equivalent diameter and ∆x is the grid spacing).

Figure 7 .

 7 Figure 7. Instantaneous scalar fluctuation field c associated with D conv for a free array of 8 bubbles, at Pe = 10 -1 (left) and Pe = 10 6 (right). The gradient of c is vertical (the entire cell is shown, and gravity is pointing downward).

Figure 8 .

 8 Figure 8. Instantaneous scalar fluctuation field c associated with D conv ⊥ for a free array of 8 bubbles, at Pe = 10 -1 (left) and Pe = 10 6 (right). The gradient of c is horizontal (the entire cell is shown, and gravity is pointing downward).

Figure 9 .

 9 Figure 9. Spatial convergence for an ordered array of bubbles in case E1 at Pe = 10 3 : relative error in D conv and D conv ⊥ as a function of the grid spacing ∆x (d b is the bubble volume-equivalent diameter; D conv ∆x=0 is extrapolated assuming D conv = D conv ∆x=0 -k∆x n , where the values of the three parameters D conv ∆x=0 , k and n are fitted from numerical data).

Table 1 .

 1 Simulated flow configurations: Bo and Ar define the flow regime, N b is the number of free bubbles in the unit cell, φ is the gas volume fraction. The resulting bubble Reynolds number (Re) and shape are also provided.

	case Bo Ar N b	φ	Re	bubble shape
	S0 0.38 0.15 1	0.002 0.00164 spherical
	S1 0.38 5.03 1	0.002 1.72	spherical
	C	243 15.2 1	0.002 9.44	skirted
	E1 2.0 29.9 1	0.002 39.9	ellipsoidal
	E1 2.0 29.9 [1, 12] 0.024 ≈ 30	ellipsoidal