Modelling of Bubbly Flow Using CFD-PBM Solver in OpenFOAM: Study of Local Population Balance Models and Extended Quadrature Method of Moments Applications

: In order to optimize and design new bubbly ﬂow reactors, it is necessary to predict the bubble behavior and properties with respect to the time and location. In gas-liquid ﬂows, it is easily observed that the bubble sizes may vary widely. The bubble size distribution is relatively sharply deﬁned, and bubble rises are uniform in homogeneous ﬂow; however bubbles aggregate, and large bubbles are formed rapidly in heterogeneous ﬂow. To assist in the analysis of these systems, the volume, size and other properties of dispersed bubbles can be described mathematically by distribution functions. Therefore, a mathematical modeling tool called the Population Balance Model (PBM) is required to predict the distribution functions of the bubble motion and the variation of their properties. In the present paper, two rectangular bubble columns and a water electrolysis reactor are modeled using the open-source Computational Fluid Dynamic (CFD) package OpenFOAM. Furthermore, the Method of Classes (CM) and Quadrature-based Moments Method (QBMM) are described, implemented and compared using the developed CFD-PBM solver. These PBM tools are applied in two bubbly ﬂow cases: bubble columns (using a Eulerian-Eulerian two-phase approach to predict the ﬂow) and a water electrolysis reactor (using a single-phase approach to predict the ﬂow). The numerical results are compared with measured data available in the scientiﬁc literature. It is observed that the Extended Quadrature Method of Moments (EQMOM) leads to a slight improvement in the prediction of experimental measurements and provides a continuous reconstruction of the Number Density Function (NDF), which is helpful in the modeling of gas evolution electrodes in the water electrolysis reactor.


Introduction
In recent years, there has been a growing interest in two-phase flows as they can be observed in various industries such as petroleum, mining, chemical and biotechnology [1][2][3][4]. In a two-phase flow, the dispersion of the particles (bubbles) in the continuous phase plays an important role in the process. One of the main challenges is to use the population balance model to predict the evolution of the number of bubbles and of their size in bubbly flow. To give an example, in order to optimize and enhance the efficiency of electrochemical reactors, the evolved bubbles should be under control, since the presence of bubbles has favorable and unfavorable effects, such as a resistive film for the electric current flow and ions supply heading to the electrode to participate in the electrode reaction.
Bubbles can also play a role as a turbulence promoter and inducer of convection over the electrode surface [5][6][7][8]. The accurate description of bubble size evolution is then beneficial to formulate predictive models for design optimization.
In the literature concerning the simulation of gas-liquid systems, several authors have used a monodisperse model only accounting for a mean bubble size (Deen et al. [9]; Holzinger [10]; Friberg [11], Morud and Hjertager [12]; Ranade and Deshpande [13]; Ranade [14]; Schwarz and Turner [15]; Schwarz [16]). While these models reduce the computational cost of the numerical simulation, they are unable to predict the evolution of the bubble size distribution, limiting their applicability and reliability for industrial design purposes [17]. One of the main problems in modeling of multiphase dispersion is the prediction of bubble size distribution. Global polydisperse models, in which a spatially-varying field of bubble sizes is considered, can globally account for the polydispersity of the bubbles, in addition to the phase continuity equation. As global polydispersity models yield a non-constant bubble size, the disperse phase is global polydisperse with a locally monodispersed size distribution (Kerdouss et al. [18] and Ishii et al. [19]). However, by applying this approach, the local probability distribution of the bubble size is not considered. Detailed models using the locally polydisperse approach give more information on the secondary phase behavior (Dhanasekharan et al. [20]; Venneker et al. [21]). The most used methods of locally polydisperse models are the Method of Classes (CM) (Balakin et al. [22]; Bannari et al. [1]; Becker et al. [23]; Kumar and Ramkrishna [24]; Kumar and Ramkrishna [25]; Puel et al. [26]), Quadrature Method of Moments (QMOM) (McGraw [27]; Marchisio et al. [28]; Marchisio et al. [29]; Marchisio et al. [30]; Sanyal et al. [31]) and Direct Quadrature Method of Moments (DQMOM) (Silva and Lage [32]; Selma et al. [33]; Marchisio and Fox [34]). The method of classes, while intuitive and accurate, is computationally intensive due to the large number of classes required to finely discretize the Number Density Function (NDF) with a large number of classes. Compared with CM, the QMOM can consider a wide range of bubble sizes with a reduced number of equations for the moments of the NDF. However, in some evaporation and combustion problems (Fox et al. [35]; Yuan et al. [36]), the value of the NDF for null internal coordinates needs to be known, which is not the case if the QMOM method is used. DQMOM solves the equations for weights and abscissae directly. Shortcomings related to the conservation of moments affect the DQMOM approach since weights and abscissas are not conserved quantities (Yuan et al. [36]). In order to overcome these limitations, Yuan et al. [36] introduced the Extended Quadrature Method of Moments (EQMOM), which enables the shape of NDF to be reconstructed from a moment set using continuous kernel density functions instead of Dirac delta functions.
QMOM and EQMOM have been recently compared in the study of liquid-liquid dispersion in a stirred tank [33,37]. Li et al. [37] performed a comparison between QMOM and EQMOM in turbulent liquid-liquid dispersion, and they observed that these two methods provided similar predictions. They managed to reconstruct the droplet size distribution by using EQMOM. CM and DQMOM have been compared in bubbly flow by Selma et al. [33], as well. Their study was carried out based on two test cases, one involving a bubble column and the other a stirred tank. They reported that it is more computationally efficient to use DQMOM compared to CM even with a relatively low number of classes. A summary of some works concerning the coupling of E-E with Population Balance Model (PBM) and the comparison among PBMs are provided in Tables 1 and 2, respectively.  In the present work, we simulate a gas-liquid flow in two rectangular bubble columns using CM, QMOM, DQMOM and EQMOM. The numerical results obtained in these two cases are used to compare the four solution methods for the PBE. The implementation of EQMOM provided by OpenQBMM [41] is coupled to a two-fluid solver in OpenFOAM to describe the bubble evolution in bubble columns. The EQMOM approach, coupled with a single-phase CFD solver in OpenFOAM, is then used to describe the evolution of the bubble phase in an electrochemical cell. This simplification, in the case of the electrochemical system, is possible because the bubble sizes of interest in this system are of the order of a micron, significantly smaller and having a lower Stokes number than in the case of the bubble column. In this case, then, a one-way coupling approach is sufficient.
The remainder of this article is organized as follows. First, the modeling approach is explained in Section 2. The numerical solution method used in the present work is summarized in Section 3. The numerical results obtained for the bubble columns and the electrochemical cell are presented in Section 4, comparing them to the experiments for validation purposes. Conclusions are drawn in Section 5.

Numerical Model
The governing equations of single-phase and two-phase system (two-fluid model) used in this model are shown in Table 3. Table 3. Governing equations of the CFD-PBM model.

Population Balance Modeling
The evolution of the size distribution of a particle population, which may consist of solid particles, bubbles or droplets, is described studying the changes in space and time of the Number Density Function (NDF) n(ζ; x, t). Here, ζ indicates the internal coordinate representing the size of the discrete element of the disperse phase, x is the position vector in physical space and t is time.
Assuming the velocity of the disperse phase is known, the evolution of the NDF is defined by the population balance equation (Marchisio et al. [28], Marchisio and Fox [49], D. Ramkrishna [50]): where U G is the velocity of the disperse phase, Γ is the diffusivity and G(ζ) the continuous rate of change in the space of internal-coordinate. The first term of Equation (1) represents accumulation; the second term describes convection; and the third diffusion in physical space. The source terms B br (ζ; x, t), B br (ζ; x, t), D ag (ζ; x, t) and D ag (ζ; x, t) are birth rate due to breakage, birth rate due to coalescence, death rate due to breakage and death rate due to coalescence, respectively. Finally, N(ζ; x, t) is the rate of change of the NDF due to nucleation. In the present study, the NDF has the form of being length-based and G(ζ) is defined as growth rate. Bubbles may nucleate when the liquid is supersaturated with gas. When the dissolved gas concentration reaches a critical value, bubbles nucleate. This critical value might be theoretically obtained from Classical Nucleation Theory (CNT) (Abraham [51]). The nucleation theory provides information about the generation of nuclei (the formation of cluster of molecules after reaching some critical size) per unit time and volume of the liquid as a function of the local parameters. The bubbles can grow in size if the liquid is saturated with dissolved gas. The remaining gas is transported into the liquid on a molecular level and gives rise to supersaturation that causes an increases of the growth of bubbles that move through liquid. In this case, the growth term is included in the PBE as a source term. In some systems such as a gas-evolving electrode, bubble growth and nucleation terms might be treated as boundary conditions, since the sources of nucleation are usually irregularities of the electrode surfaces, and once a nucleus exists, the bubble growth occurs on the electrode surface in a concentration boundary layer. Gas bubbles develop at nucleation sites on the electrode surface, grow in size until they reach a critical break-off diameter and detach into the electrolyte afterwards (Nierhaus [3], Tomasoni [52]). In other words, the bubble formation happens on the electrode. Consequently, it can be described by means of a boundary condition (an ingoing flux boundary).
Growth, nucleation and diffusivity were not taken into account in the example applications presented in this work, since the focus is on bubble coalescence and breakage. Therefore, the final form of the evolution equation of the NDF is: The breakage and coalescence source terms are modeled as (Marchisio et al. [28], Marchisio and Fox [49]): Here, β(ζ, ζ ) is the coalescence rate between bubbles of size ζ and ζ ; a(ζ) is the break-up frequency of a bubble with size ζ; b(ζ|ζ ) represents daughter distribution function generated from the breakup of a bubble of size ζ .

Class Methods
The class method solves the bubble' number density, directly [50]. In CM, the continuous size range of bubbles can be realized through the discretization of the bubbles size distribution into a number of classes of discrete sizes. For each class, the equation of the number density of bubbles is solved and coalescence and breakup rates are transformed into birth and death rates. The population balance equation for the i-th bubble class is represented as: where n i is the number of the bubbles from group i per unit volume. B br and B agr are the birth rates caused by breakup and coalescence, respectively, and D agr and D br the corresponding death rates, from coalescence and breakage, respectively. The relationship between the volume fraction and the number density is: where ν i is the volume of a bubble of class i.
where α G is the volume fraction of the dispersed phase.
To solve the population balance equation using the scalars f i , Equation (7) is changed to the following equation: Here, f i is the bubble volume fraction of group of size i. As Equation (10) shows, all bubbles in a computational cell move with the same velocity U G . This approach is called the MUlti SIzeGroup (MUSIG) [53]. It is worth observing that assuming all the bubbles in a computational cell move with the same velocity is a limitation of the approach, which may be acceptable for narrow bubble size distributions, but not in general. A class of Quadrature-based Moments Methods (QBMM) that is not affected by this limitation was proposed by Yuan et al. [54], and its adoption will be the topic of future work.
The method of classes calculates the Sauter mean diameter d 32 as follows: The Sauter mean diameter represents the average bubble size, and it is defined as the diameter of a sphere that has the same volume to surface area ratio as the bubble under consideration. The Sauter diameter is often used in problems where the active surface area is the relevant parameter [49].
Ramkrishna [50] proposed an approach named fixed pivot in order to discretize the source terms in the PBE. The approach assumes that the population of bubbles is distributed on pivotal grid points x i with x i+1 = sx i and s > 1, where i refers to the class i with i < n. Assuming spherical bubbles, In the CM technique, the bubble size is divided into n = 2r + 1 classes, where n is odd in order to have symmetrical divisions. This gives the following relation: Bannari et al. [1] evaluated the the values of r and s for different classes assuming d max = 5 mm and d mean = 10 mm (Table 4). Breakage and aggregation may create bubbles with volume ν such that x i < ν < x i+1 . This bubble must be split by assigning respectively fraction γ i and γ i+1 to x i and x i+1 . The following limitations preserve the number balance and mass balance.
Ramkrishna [50] has also reported the birth in class i due to coalescence in this way: where θ is a test function expressed as: and: The death rates D i,agr in class i due to coalescence is defined as follows: while the birth rate B i br in class i due to breakup as: and the death rate D i,br in class i due to break-up: where: The mentioned integrals are solved by the Gaussian quadrature integration as follows: P n is a Legendre polynomial and can be formulated using the recurrence relation: W j is the weighting function of the orthogonal polynomials as shown in Table 5. If j equals five, adequate accuracy is achieved. Boundary and initial conditions to solve CM equations are shown in Table 6. It is assumed that the bubble size at the inlet is monodispersed and spatially uniform, equal to the mean diameter (d 2r ). The first-order upwind scheme is used to discretize the advection term of the f i equations.
Evolution equations for the moments of the NDF read: By solving Equation (26) for a set of at least four moments, the Sauter mean diameter d 32 = m 3 /m 2 can be calculated.

Quadrature Method of Moments
In the QMOM technique, the unknown NDF is represented by a weighted summation of Dirac delta distributions δ(L − L α ): where W i are non-negative weights of each kernel density function, L i are the corresponding quadrature abscissae and N is the number of kernel density functions to approximate the NDF. Source terms in the moment Equation (26) become: where the N is the number of weights w α , and the corresponding abscissae L α are determined from the first 2N integer moments of the NDF. β αβ is the aggregation kernel for the bubbles of size L α and L β ; a α is the breakage kernel for the bubble size of L α ; and b (k) α represents daughter bubble distribution function.
Boundary conditions and initial conditions to solve the moment equations are shown in Table 7. Gauss upwind is utilized as the divergence scheme. If the bubble size, which can be represented by a constant value or a distribution function, is known, the moments will be easily calculated. The moments for constant bubble size are found by: where ζ is bubble size and α G is gas phase fraction. DQMOM is based on the direct solution of the transport equations for weights and abscissas of the quadrature approximation (Fan et al. [55]). The transport equations for weights and abscissas are written as: where a i and b i are found by the solution of the following linear system in terms of the unknown a i and b i (Bove [56]): The source terms S (k) are taken to be the as same as for the QMOM method. Boundary and initial conditions to solve for weights and abscissae are shown in Table 8. The first-order upwind scheme is used to discretize the advection term in both Equations (33) and (34). The linear system obtained from Equation (35) is solved using the Gauss-Seidel technique. The EQMOM approach approximates the unknown NDF with a weighted sum of smooth, non-negative kernel density functions δ σ (L, L α ) [36,57]: In this work, the log-normal kernel density is used [58]: Source terms are then closed in terms of the primary and secondary quadrature found with the EQMOM procedure, leading to: where the N primary weights W α , the corresponding primary abscissae L α , together with the parameter σ are determined from the first 2N + 1 integer moments of the NDF. The 2N α quantities W αγ , called secondary weights and abscissae, respectively, are computed using the standard Gaussian quadrature formulae for known orthogonal polynomials to the kernel NDF [36,57]. β α 1 γ 1 α 2 γ 2 is the aggregation kernel for the bubbles of size L α 1 γ 1 and L α 2 γ 2 ; a αγ is the breakage kernel for the bubbles size of L αγ ; and b αγ represents the daughter distribution function. Boundary conditions and initial conditions used to solve the moment equations are the same as those used for the QMOM approach (Table 7).

Closure Models for Coalescence and Breakage
There are many breakage and coalescence kernels available for bubbly flow, but they are essentially written in a similar form with some minors differences in the model constants or assumptions [38]. The discussion of aggregation and breakage kernels is beyond the scope of this paper. Therefore, it was decided to utilize the ones (Luo and Svendsen [46] and Hagesather et al. [47]) that have been widely adopted for the bubble columns (Bannari et al. [1], Kerdouss et al. [18], Kerdouss et al. [59], Kerdouss et al. [60] and Selma et al. [33]). The equations of coalescence and breakage kernels are shown in Table 3.

Numerical Solution
The models described in the present work were solved using the OpenFOAM library [61]. Two-fluid model and single-phase flow libraries in OpenFOAM (twoPhaseEulerFOAM and pimpleFoam) were customized in order to couple the PBE approaches used in this work to the two -phase flow model, to simulate bubbly flows in bubble columns, and to the single-phase flow model (dilute system), to simulate the electrochemical cell. Equations are solved iteratively, relying on the semi-implicit PIMPLE (merged PISO-SIMPLE) approach provided by OpenFOAM, which is a combination of the PISO (Pressure Implicit with Split Operator) and SIMPLE (Semi Implicit Method for Pressure Linked Equations) procedures.

Test Case 1: Pseudo-2D Bubble Column
The first test case for the coupled Population Balance Model and Two Fluid Model (PBM-TFM) approach is a simple geometry (a rectangular bubble column [1]). The same test case was used in [33] to compare CM and DQMOM. Consequently, only the simulations with EQMOM were necessary to allow the comparison.
The dimensions and boundary conditions used to perform the simulations are shown in Figure 1, while Table 9 summarizes the models used in the simulations for this test case. Tables 10-12 demonstrates the boundary and initial conditions applied in CM, QMOM, EQMOM and DQMOM, respectively.

Inlet Wall
Outlet Inlet value Neumann  Table 11. Boundary and initial conditions for m i equations used in the QMOM and EQMOM methods (Test Cases I and II).

Boundary Conditions Initial Condition
Inlet Wall Outlet Inlet value Neumann Neumann Table 12. Boundary and initial conditions for W i and L i equations used in the DQMOM method [33].

Inlet Wall
Outlet Neumann Neumann Figure 2 depicts a comparison between experimental measurements [62] and axial liquid velocity provided by the CFD-PBM solver on a line along y = 37 cm. The comparison confirms that the solver using EQMOM works properly for the bubble column. Figure 2a shows that the variation of secondary nodes does not have a significant impact. However, the increase of primary nodes from two to three provides better accuracy, as was expected. Figure 2b demonstrates the differences among population balance models. It is observed that EQMOM using three nodes represents a slight improvement in comparison with CM (25 classes) and DQMOM (three nodes). It is noteworthy that the CM (25 classes) computationally is very expensive [33]. By contrast, DQMOM and EQMOM have less computational demand.  Figure 3 demonstrates the radial bubble segregation based on the Sauter mean diameter contour. In fact, the bubble size decreases heading toward the wall due to breakage phenomena. However, the high gas phase fraction and coalescence events at the center of column lead to the creation of larger bubble sizes. As Figure 3 shows, Sauter mean diameter distributions predicted by the applied PBMs (DQMOM, EQMOM, CM) are qualitatively confirmed by comparison with experimental observation.

Experiment
DQMOM [33] EQMOM CM EQMOM is capable of providing a smooth reconstruction of the NDF for each arbitrary zone or cell in the computational domain. While this reconstruction is not unique for the set of transported moments, examining the reconstructed NDF can provide a better insight into the behavior of the dispersed phase. In order to plot NDF, weights and abscissae calculated by EQMOM and extracted from CFD-PBM solver are utilized to approximate Equation (37). Figure 4 indicates the shape of NDF in the liquid phase (water zone) at t = 60 s. It can be observed that the mean bubble diameter ranges from 1 mm to 6 mm, and thus, the monodispersity is not a suitable approximation. As can be observed from Figure 4a, the NDF is represented by a Dirac delta function for N = 2, while Figure 4b exhibits a continuous NDF for N = 3. This means the increase of the nodes in EQMOM configuration applied in the bubble column leads to the creation of a spread in bubble sizes, even locally.

Test Case 2: 3D Bubble Column
For the second case, the coupled PBM-TFM solver was tested in a square bubble column (Deen [64]). The models used in the simulations are summarized in Table 13. This configuration is based on the work of Holzinger [10] in which the same test case was studied using a monodisperse bubble size distribution. In Test Case 2, QMOM and EQMOM simulations were performed.
The column has a square cross-section with W (width) = D(depth) = 0.15 m. The sparger is in the form of a square, the are of which is A in = 0.03 × 0.03 m. The dimensions and boundary conditions are shown in Figure 5. The domain is discretized into 15 × 15 × 60 control volumes, a total of 13,500 cells. Figure 6 shows the comparison between QMOM and EQMOM with the experimental measurements reported by Deen [64]. Three nodes were used for QMOM, and three primary nodes were used for EQMOM. The EQMOM approach yielded a minor improvement in the gas velocity profile, while no change was observed through the liquid velocity profile. Likewise, the color maps of the Sauter mean diameter are similar for both methods, as shown in Figure 7. In spite of the similarity, the computational cost of EQMOM is higher than the QMOM method, as Madadi-Kandjani and Passalacqua [58] reported for a zero-dimension case, as well as the current simulation process.

Test Case 3: Water Electrolysis Reactor
The bubbly flow in a water electrolysis reactor consisting of a gas evolution electrode was chosen as the third case. This decision is based on the simplicity of the case and on the availability of experimental data for the bubble size distribution. It was decided to investigate the Inverted Rotating Disk Electrode (IRDE) proposed by Van Parys et al. [65], which is composed of three electrodes as follows: When a potential difference is applied between electrodes, hydrogen is produced at the cathode and oxygen at the anode in the form of bubbles. In this case, we put aside the influence of the anode, because the experimental bubble distribution is only available over cathode and hydrogen bubbles. Consequently, the impact of oxygen bubbles is not taken into account in the simulation. The specifications of the IRDEreactor are presented in Figure 8. In this system, the rotating cathode is considered as an inlet boundary with imposed velocity, volume fraction and bubble size distribution. The angular component of the liquid velocity is set according to the rotational speed of working electrode (ω = 100 and 250 rpm). According to the applied rotational speeds, the calculated Reynolds number is less than the critical ones [66]. Therefore, the flow regime is laminar, which allows neglecting source terms in the population balance model because of = 0. In the present study, gas hold-up is extremely low. The bubbles follow the bulk flow and are affected by the continuous phase (water), but not vice versa. Hence, the flow field was calculated with a single-phase solver. The flow field obtained from the single-phase approach was then imported in the population balance solver, in order to advect the bubble size distribution imposed at the inlet and study how bubbles distribute in the IRDE. The settings, boundary and initial conditions used in the simulation are summarized in Tables 14 and 15.  A distribution of bubble sizes is observed at the electrodes of the IRDE reactor. For this reason, the continuous distribution function reported by Nierhaus et al. [66] was used in the simulation and imposed at the electrode surface, which is treated as an inlet boundary for the gas phase (Figure 9a).
A grid with 34,300 hexahedra yielded sufficiently accurate results for the proposed problem and was selected for the IRDE study.
Nierhaus et al. [66] reported the experimental bubble size distributions in an optical window (W1) for two different rotating velocities, 100 rpm and 250 rpm. Figure 8 illustrates how the W1 volume has been configured in the IRDE reactor. W1 was located above the electrode to enable tracking bubbles in the rising plume.  √ ω z ν ) as a function of the dimensionless height (γ = z δ ), where δ is the displacement thickness of the fluid boundary layer (δ = ν ω z ). The comparison shows that the numerical results match with the analytical solution [68] in the region close to the electrode. The confirmed flow field applied in population balance calculations for its NDF consists of the accumulation and convection term (physical space). Table 15. Boundary and initial condition for m i equations used in Test Case 3.

Boundary Conditions Initial Condition
Inlet Wall Outlet Inlet value Neumann Neumann To validate the CFD-PBM solver in IRDE, an analysis is applied to investigate the bubble size distribution in volume W1. The comparison of the bubble size distribution between the experimental study and current CFD-PBM using EQMOM (three nodes) is presented in Figure 10. The results thus obtained are compatible with experimental measurements (Figure 10a,b). The fair agreement confirms the assumptions in the PBM model, particularly for 150 rpm. In fact, since the bubble size is so small (low Stokes number), most of the effect is due to advection, and no segregation occurs.
In the case of 150 rpm, simulation data better match experimental data, while there is a disagreement for a higher rotational speed. The prime cause of the discrepancy might be the result of the neglect of size change effects in the PBM model.

Conclusions
In this paper, an analysis was performed for comparison among local Population Balance Models (PBM), CM, QMOM, DQMOM and EQMOM, based on three cases in the presence of bubbly flow. The purpose of the paper is to validate the CFD-PBM solver, which was applied for two different bubble columns and a water electrolysis reactor to predict the bubble size distribution. The originality of the OpenFOAM solver lies on the fact that it employs the novel method of PBM, which can accept the continuous bubble size distribution as a boundary condition. Moreover, the solver is able to export the distribution function for a specified region in an arbitrary time based on the EQMOM method. It was observed that the CFD-PBM using EQMOM provides a reasonable prediction, as well as CM consisting of 25 classes, but requires less computational demand compared with CM (Table 16). From the research that has been carried out, it is possible to conclude that QMOM and EQMOM have similar predictions. EQMOM is computationally more expensive than QMOM, although it is able to obtain a continuous NDF of the model. In order to simulate bubbly flow in the bubble column, it is proposed to use at least three nodes in the EQMOM technique. This minimum value is required to acquire a continuous NDF function. It is evident that the experimental and numerical NDF can be compared in the water electrolysis reactor (IRDE). Agreement is achieved using EQMOM in the IRDE reactor. The results proved that the most dominant term is advection in the PBM model.