A fully Lagrangian computational model for the integration of mixing and biochemical reactions in anaerobic digestion

The impact of mixing on biochemical reactions is of high importance in anaerobic digestion (AD). In this paper, a novel 2D fully Lagrangian computational model for the integration of mixing and biochemical reactions in AD is developed and presented. The mixing-induced ﬂuid ﬂow is modeled by smoothed particle hydrodynamics (SPH). The computational domain is discretized by SPH particles, each of which carries the information of biologically active compounds and follows the ﬂow ﬁeld. In this natural way, advection is reproduced, which is the main advantage of SPH for this type of problems. A mathematical model that governs the biochemical reactions is integrated in time for each particle, which allows to spa- tially resolve the biological concentrations. Mass transfer interactions between particles are reproduced by the diffusion equation to directly link mixing to biochemical reactions. The total biogas production is obtained by integrating over all the particles. Both SPH and biochemical models are veriﬁed against existing data in the literature and the integrated model is then applied to a real world anaerobic digester. The application of a novel fully Lagrangian method to AD is a stepping stone to future possible developments. However, in the simulation of such problems, SPH is still uncompetitive if compared to other mainstream methods and industrial application of the model depends on the computational eﬃciency of future SPH solvers. © This is the CC BY ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Anaerobic digestion (AD) is a widely applied technology for advanced treatment of biodegradable materials. Due to the biochemical processes in AD, biogas is ultimately formed under anaerobic conditions from the biological residues, subsequently utilized to generate energy. Although the environmental and economic advantages of AD are undisputed, the net energy yield of AD is not overwhelming compared to other renewable energy sources and hence the AD energy consumption needs to be optimized [1] .
The performance of AD depends on several physical and biochemical parameters, mainly including biomass characteristics, feeding patterns, temperature, pH, hydraulic retention time (HRT) and mixing. The individual steps of the AD process have been sub-ject of considerable research activity in recent decades. Comprehensive mathematical models have evolved for the biochemical reaction kinetics [2][3][4] , while a large number of studies have been dedicated to investigate flow patterns and mixing. Early experimental studies on mixing methods in AD are limited to tracer investigations [see e.g. [5][6][7] ], where a chemical material is injected into the digester and its evolution is investigated over time. The rate and the pattern of this evolution reveals the mixing effects in a simplified manner. Non-invasive experimental methods like computer automated radioactive particle tracking (CARPT) and computer tomography (CT) have also been employed to characterize the flow pattern inside lab-scale digesters [8] . However, experimental techniques are costly, resource-intensive and not applicable for all full-scale digesters.
Thanks to modern high performance computers, computational fluid dynamics (CFD) methods have emerged as promising tools to study the fluid flow in bioenergy systems [9] . Vesvikar and Al-Dahhan [10] performed a steady state CFD simulation of a lab-scale gas-lift digester, resolving water, air and the lifting gas as separate media. CFD has been also used to simulate circular and egg-shaped digesters with regard to various real-life experimental cases to calculate velocity fields, particle trajectories, dead zones and turbulence impacts [11] . Lab-scale digesters have been also studied by Bridgeman [12] using CFD, through which new findings for the velocity distribution and the biogas yield are achieved. Full-scale digesters are also studied using CFD [13,14] . Due to the importance of turbulent mixing in AD, Coughtrie et al. [15] compared different turbulence models for a bench scale gas-lift digester. They simplified the effects of the gas phase by an empirical velocity formulation. CFD has been also used for optimization of wastewater digestion by Hurtado et al. [16] . They studied the sludge flow behavior under steady and non-steady conditions and determined an intermittent operating regime with a high mixing quality. Zhang et al. [17] used CFD to investigate the mixing mode and power consumption in an unbaffled stirred tank bioreactor. Recently, CFD has been also employed to assess mixing quality in full-scale biogas-mixed anaerobic digestion [18] .
Despite the remarkable improvements in design and optimization of AD offered by the aforementioned studies on either physical or biochemical compartments, the intrinsic connection between the flow patterns (mixing) and the biochemical reactions is still poorly studied [9] . As an energy-intensive factor of AD, mixing should be optimized in order to maximize the net energy yield [19] . The direct impact of mixing on the biogas production has to be known to achieve an optimized process. Therefore, a fully integrated model is required to study the relation between mixing and the biochemical reactions. However, to our knowledge, only Wu and Chen [20] and Wu [21] have addressed this topic with a computational model. The former propose a model, which calculates the biomass resident time for each computational cell by solving a scalar transport equation based on a fully converged hydrodynamics solution and expresses the methane yield rate as a linear function of temperature. The latter similarly calculates first the resident time and accordingly solves the detailed biokinetic mathematical equations for the determined resident time to ultimately compute the total methane yield. Wu [21] validated the computationally predicted methane yield against experimental data. However, the proposed model therein is constrained by the fixed Eulerian grid points, which restricts the extension of the model to encompass the multi-phase nature and inherent heterogeneities of the biomass fed into digesters.
Since meshless Lagrangian computational methods show more accurate field-scale predictions of reactive transport processes [22] , in the present study smoothed particle hydrodynamics (SPH), as the most advanced member of this category, is considered as the ideal choice for being coupled to the biological model. SPH has already been employed to introduce integrated models for wastewater treatment by coupling hydrodynamics and biokinetics [23] . However, to our best knowledge, this is the first-ever application of a fully Lagrangian method to couple mixing and biochemical reactions in AD. There exist Euler-Lagrange CFD models in the literature, which have benefited from a Lagrangian approach to simulate the multi-phase medium of gas mixing in AD [18,19] . Very recently, Dapelo et al. [24] also used the Euler-Lagrange approach in a Lattice Boltzmann framework, which considerably accelerates the computations.
In this paper we propose a direct link between SPH as the hydrodynamics model and the versatile AD models [2,3] governing the biochemical reaction. The concentrations of biological compounds are assigned to SPH particles, which follow the flow field. In the presented 2D model, the biochemical reactions are solved for every single particle such that each particle behaves as an individual bioreactor. In order to realize the local mass transfer interactions, the diffusion equation is discretized and solved over the particles by SPH, thereby the spatial interactions between differ-ent compounds are taken into account. Finally, the overall biogas production is obtained by integrating over all the particles.
Another decisive factor in the performance analysis of AD is sedimentation. As such, the overachieving goal must be developing a model capable of mimicking multiple physical phases. The meshless nature of SPH is a proven advantage in the simulation of multi-phase flows [25,26] . In contrast to mesh-based Eulerian computational methods, like finite volume method (FVM), SPH needs no specific algorithm to identify the multi-phase interface, since different phases are naturally distinguished by particles. However, in multi-phase problems characterized by large density differences, which are not present in our study, the numerical instabilities caused by abrupt density discontinuities across the interface need to be overcome by special interface treatments [27][28][29] . The mentioned beneficial features of SPH make it possible to develop the approach further towards a fully integrated numerical model of the AD process including sedimentation effects. This is an asset to SPH, which is not found in previous numerical simulations of AD by mesh-based methods [21] . It is worth noting that the dissipative viscosity treatments in traditional SPH [30] impact turbulence quantities that may be important in mixing over the time scales of AD. SPH is also regarded to be well suited for fee-surface gravity dominated flows. However, the aforementioned assets still nominate SPH as an attractive method to study the AD process. Overall the advantages of using SPH in this type of processes can be summarized here as 1. Advection is reproduced in a natural way. 2. Reactive transport processes are predicted more accurately in comparison to mesh-bashed methods [22] . 3. The Lagrangian nature allows to naturally couple mixing to biochemical reactions. 4. Carrying the biological information on SPH particles spatially resolves the biological concentrations. 5. SPH allows to develop models capable of mimicking multiple physical phases (viz. flocculant and sedimenting matter) by naturally realizing multi-phase interfaces.
Despite the opening above-mentioned peculiarities offer with regard to future developments, industrial application of the proposed model depends on computational efficiency of the future SPH solvers. It is generally observed that 3D simulations of AD are of importance due to turbulence and rheological effects, which are costly prohibitive within SPH with the lengthy AD processes in mind. Nevertheless, extending the present model to non-Newtonian flows with the consideration of heat transfer in 2D is our ongoing work for more in-depth studies. We showed a first attempt in [31] .
The paper is organized as follows: Section 2 describes the basic principles of SPH together with the mathematical models for both the fluid flow and the biological processes. The employed numerical scheme to solve the time-dependent models is then summarized in Section 3 . Section 4 includes the validation of the mathematical models, mesh convergence studies and the SPH simulations of a real life anaerobic digester in comparison with FVM-based CFD simulations. Finally, the concluding remarks of the present study are highlighted in Section 5 .

Fundamentals of SPH
SPH is a fully Lagrangian method which describes the fluid flow by a finite number of particles. The particles carry the information regarding mass and other physical parameters and follow the flow field, while they are mathematically considered as interpolation points. The main idea of SPH is to interpolate physical quantities over the neighboring particles of a particle of interest using a weighting function W . According to this principle, any arbitrary function f ( r ) at position r in space is approximated by an integral formulation [32] where angle brackets denote the approximation of f ( r ) obtained from using the weighting function instead of the Dirac delta function, δ(r − r ) , h is the smoothing length, is the support domain of W , and δV is the differential volume of the particles. The weighting function, which in SPH literature is called as smoothing kernel, has to satisfy the normalization condition, W (r − r , h ) δV = 1 and the symmetry condition. In addition, it should approach the Dirac delta function as the support domain tends to zero, lim h → 0 W (r − r , h ) = δ(r − r ) . Representing the computational domain by particles and considering a set of particles inside the support domain of the particle i , one can write the discrete approximation of Eq. (1) where for the sake of simplicity W (r i − r j , h ) is substituted by W ij .
There are various smoothing kernels in literature and throughout the present work the Wendland kernel [33] is used for all numerical calculations. A smoothing length of h = 2 dx is used, where dx is the initial particle spacing.

Governing equations
The continuity equation relates the change of density to the advection of the continuum and is discretized in SPH form as where t, ρ and m represent the time, density and mass, respectively and v i j = v i − v j is the relative velocity between particle i and j . The momentum conservation equation in a Lagrangian framework has the following form (5) and is discretized using the formulation presented by Adami et al. [34] as where r i j = r i − r j , r i j = r i j and ∂W The coefficients h ij , c ij and ρ ij are the averages of the smoothing length, speed of sound and density between the two particles. Variables V i and η i represent the volume and dynamic viscosity of particle i , the parameter = 0 . 01 is included to ensure a non-zero denominator and the constant α is chosen such that the solution is not affected by the artificial viscosity.

Advection-diffusion equation
Diffusion is the only method of transporting material from the biomass particles into the bulk medium in a laminar flow [35] . In order to resolve the mass transfer interactions at a particle level, in the present work the advection-diffusion equation is solved over the particles. The original form of this equation includes Fickian diffusion and the advection term dC dt where C is the concentration of a compound, D is the diffusion coefficient, and v is the velocity field. There are several formulations in the SPH literature to discretize Eq. (7) [see e.g. [32,36,37] ]. The accuracy analysis of Aristodemo et al. [37] shows that their formulation leads to more accurate results. For this reason, here we use the formulation proposed in [37] to discretize Eq. (7) in the SPH form It should be noted that since the velocity fields in the herein considered problems are not strong, we can neglect the advection term of the Eq. (8) . The concentrations are naturally advected with the particles due to the Lagrangian nature of the SPH method. This characteristic is another advantage of the SPH method for the simulation of processes dealing with reactiondiffusion phenomena like AD.

Anaerobic digestion model
The anaerobic digestion model No. 1 (ADM1) [4] has been reported as the most comprehensive mathematical model for the AD process. However, as also mentioned in [38] and [21] , because of the large number of constants and coefficients in this model it is nearly impossible to validate the model according to any of the available experimental data sets. On the other hand, ADM1 is based on a completely stirred reactor with a constant input and output stream for the whole reactor. The constant input and output rates together with the total volume of the reactor are included in all the biochemical mass balance equations, including the methane production rate, which is not compatible with our particle-based description of the system. In our approach, the input and output streams are emulated by the SPH particles, where at each time instant a number of particles are allowed to leave/enter the reactor such that their occupied volume is equal to the desired volumetric flow rate. Moreover, in order to couple the biochemical reactions to the physical processes we need to solve the biochemical mass balance equations locally for each particle and reproduce the mass transfer interactions between particles by the diffusion equation. With the above considerations in mind, we use the model proposed in [2] and [3] , but adopting the total input/output flow rates and the total volume.
The employed model herein includes four processes: hydrolysis of particulate substrate, acidogenesis (consumption of soluble substrates by acidogenic bacteria), acetogenesis (consumption of volatile fatty acids and formation of acetate by propionate and butyrate degrading acetogenic bacteria), and finally methanogesis (consumption of acetate and methane generation by methanogenic bacteria). These biological reactions are modeled as a system of ordinary differential equations (ODEs). Since the rates of change of different active biological compounds differ significantly, this system of ODEs is denoted stiff in mathematical terms. It is thus proposed to reduce the stiffness by approximating the rate of change of the cation of hydrogen ( H + , a very quickly changing concentration) as an algebraic non-linear equation [2] .
In the present work, as illustrated in Fig. 1 , SPH particles describe the biomass flow while carrying the information regarding the biological concentrations, C i . For each particle of interest, i , the hydrodynamic equations along with the diffusion equation are discretized using the smoothing kernel and solved over the neighboring particles which lie within the support domain, W ij (described in Section 2 ). The radius of the support domain is κh where h is the smoothing length and κ is a constant, which for the Wendland kernel is equal to 2. At each time step, for each particle, the system of ODEs is solved using an Euler forward integration and the non-linear algebraic equation for H + is iteratively solved using the Newton-Raphson method, such that each particle is performing as an individual reactor. Having determined the concentration of H + , the pH value for each particle is calculated as pH = −log 10 (H + ) .
In order to take the effects of hydrodynamics into account and achieve a coupled model, the local mass transfer interactions between each particle and its neighboring particles are considered by solving the diffusion equation.
The aforementioned asset of SPH to simulate the multi-phase medium of AD is also better comprehensible in Fig 1 . In contrast to mesh-based Eulerian methods, the Lagrangian framework of SPH lets the particles follow the flow field without the constraint of fixed computational grid points. Each SPH particle can mimic either a soluble, an inert or any non-degradable material (e.g. metal, glass, etc.). This is an advantageous feature of SPH in computational studies of AD, which makes it better suited to include the effects of sedimentation, non-degradable materials, etc.

Numerical scheme
The gpuSPHASE software [39] is used here as the computational framework in which the time integration is based on the velocity-Verlet scheme. In order to calculate the local pH value for each particle, at the beginning of each time step, the non-linear algebraic equation is solved as a function of the local concentrations of the biological compounds. In addition, to take the biochemical reactions into account the system of ODEs are solved for each particle at each time step. With the newly predicted concentrations from the biochemical reactions at the beginning of each time step, the predictor-corrector scheme is used for the diffusion equation [36] . This form of time integration introduces a hybrid time-scheme with the use of Verlet and predictor-corrector algorithms, together with an iterative solver to calculate pH values. Simplifying the solver of the system of ODEs by an explicit Euler integration, the time stepping scheme takes the following form where f (pH) is the non-linear equation of pH, C i,bio is the calculated concentration according to the biochemical reactions, f bio is the general form of the ODEs and f C is the rate of change of concentrations in time, obtained from Eq. (8) . Since the continuity equation Eq. (15) and the force calculations Eq. (17) are computationally expensive, they are performed only once per step [34] .
To maintain the stability of the numerical algorithm, the time step size is restricted by several conditions including the Courant-Friedrichs-Levy (CFL) condition [32] , the conditions imposed by viscous and gravitational forces, and the diffusion condition , where 0 < C CFL < 1 is the CFL number, which in the present work is set equal to 0.25 and c is the speed of sound. In the weakly compressible SPH (WCSPH) scheme employed herein, the fluid is assumed to be weakly compressible and an equation of state is used to calculate pressure at each time step. The speed of sound c is a decisive parameter in the equation of state and is suitably chosen to limit the density fluctuations to 1 %, thereby ensuring the incompressibility constraint. gpuSPHASE utilizes graphics processing units (GPU) in order to accelerate the calculations. All the simulations presented herein are performed in parallel using a Nvidia GeForce GTX TITAN X (Pascal architecture) graphics card installed on a workstation with 16GB RAM main memory and an Intel(R) Xeon(R) CPU E5-1620 v3 CPU with 8 cores.

Validation of the diffusion equation
It is due to the nature of large scale anaerobic digester installations that a detailed, 3D and dynamic monitoring scheme is not feasible. As real world data for full model validation do not exist, the novel submodules of the presented Lagrangian-biokinetic model are validated individually in the following.
The gpuSPHASE framework has already been validated against several experimental, analytical and numerical fluid flow benchmark cases [39] and thus the implementation of SPH code is assumed to be tested and correct. In order to validate the novel diffusion term, the analytical benchmark case for a 1D diffusion is chosen which is previously used by several SPH studies [see e.g. [36,37,40] ]. The schematic of the problem is depicted in Fig. 2 . In this case, a 0.4 × 1 m rectangle is filled with water ( ρ = 10 0 0 kg . m −3 ) and a finite horizontal band of pollutant ( ρ = 10 0 0 kg . m −3 ) is located at the middle of the rectangle confined to the region z 1 ≤ z ≤ z 2 . The initial concentration of the pollutant is equal to C 0 = 1 kg . m −3 in the band and zero elsewhere. The diffusion coefficient is equal to 10 −4 m 2 . s −1 .
For the above test case, an analytical solution is proposed in [41] as follows  The SPH numerical predictions of the concentration at t = 1 . 0 s is shown in Fig. 3 . The evolution of the concentration is well predicted and the bell shaped distribution of the concentration is in agreement with the analytical solution. This simulation was also performed with the drift-kick-drift [42] and the Verlet [43] time stepping schemes. The best results were obtained using our hybrid scheme.
This problem is further investigated with an exponential initial distribution of the pollutant all over the rectangle height. The initial distribution of the concentration has the following form, which reaches its maximum at the middle of the rectangle ( z 0 = 0 . 5 m ) and the analytical solution for this problem proposed in [41] is as follows where t 0 = 1 s , c 0 = 1 kg . s 1 / 2 . m −3 and D = 10 −4 m 2 . s −1 . Fig. 4 illustrates the comparison of the SPH numerical predictions of the concentration distribution against the analytical solution for the exponential initial distribution at t = 1 s . As for the first case with a constant initial distribution, this case was also performed with the drift-kick-drift [42] and the Verlet [43] time stepping schemes. The best results were also in this case obtained using our hybrid scheme. The convergence of the concentration distributions with increasing resolution are also shown in Figs. 3 and 4 .
In order to quantify the convergence rate of the SPH method, the L 2 norm of error [44] for the 1D diffusion test case with two different initial distributions is plotted in Fig. 5 . Four particle resolutions are used, namely, dx = 0.005, 0.01, 0.02 and 0.03 m. A near optimum rate of convergence is demonstrated in both cases with constant and exponential initial distribution. Evidently, the accuracy of the SPH scheme improves with increasing particle resolutions and the scheme converges faster with constant initial concentration distribution.

Validation of the AD model
As mentioned in Section 2.4 , the AD model used in the present work consists of a system of ODEs and a non-linear algebraic  equation. At each time step the system of ODEs is solved using the forward Euler integration and the non-linear equation for the concentration of H + is iteratively solved using the Newton-Raphson method.
In order to validate the AD model implementation, the experimental data reported by Mackie and Bryant [45] are used. A labscale digester with a volume of 0.003 m 3 in a mesophilic temperature environment (40 °C) is used in [45] . The digester has different feeding patterns according to the HRT. It has been fed in a semicontinuous way twice a day for the shortest HRT (5 day) and once a day at the longer HRTs (13, 10 and 9 days). The initial concentrations are computed using the mathematical expressions proposed in [21] and are summarized in Table 1 where X A , X AP , X AB , X M , C Is , C S , C Ac , C Pr , C But , C Am , C CO 2 and C CH 4 denote the concentration of acidogenic bacteria, propionate degrading acetogenic bacteria, bu- Fig. 6. Comparison of the computationally predicted methane yield against the experimental results [45] . tyric degrading acetogenic bacteria, methanogenic bacteria, insoluble substrate, soluble substrate, acetate, propionate, butyrate, ammonium, carbon dioxide, and methane, respectively. Fig. 6 depicts the computationally predicted methane production compared with the experimental results of [45] for different organic loading rates (OLR). These methane yield values are obtained after a 24 h simulation of the AD model in a continuously stirred-tank reactor (CSTR). As it is shown, a fairly good agreement is achieved between the experimental results and the numerical predictions according to the mathematical model.

Geometry layout
In order to compare the velocity profiles between SPH and AN-SYS Fluent in Section 4.4 , and also as a case study in Section 4.6 , in this paper we investigate an egg-shaped anaerobic digester of a real-life waste water treatment plant (WWTP). The schematic configuration of the digester is shown in Fig. 7 . The WWTP of the Achental-Inntal-Zillertal wastewater management association (AIZ) is located in the community of Strass (Austria). The AIZ treatment plant is designed for 167,0 0 0 population equivalents. The AD process is realized by two identical egg-shaped digesters with a volume of 2500 m 3 each. The digesters are operated at a temperature of 35 °C. The mixing is realized by a mechanical draft tube, which is vertically located at the middle of the digester. The mixing effects are enhanced with the recirculation of the biomass. The flow rate in the recirculation stream is 105 m 3 . day −1 . Further details about the technical and operating data of the digester are found in [46] . In the subsequent analysis, the actual operating conditions of the AIZ digester are adopted.

Comparison of the velocity profiles between SPH and 2D ANSYS
While the SPH code implementation is already validated elsewhere (see Section 4.1 ), this section verifies the resulting velocity field in a real world anaerobic digester. Detailed 3D CFD Table 1 The initial concentrations (kg.m −3 ) of the biologically active compounds for solving the system of ODEs. analysis of the mixing in the real anaerobic digester described in Section 4.3 (also considered as case study in Section 4.6 ), has been conducted by applying FVM-based simulation in [46] . However, the comparison between 2D and 3D models is difficult for a number of reasons. At first, the egg-shaped form of the digester induces a higher dissipation of the velocity of the fluid in a 3D simulation, resulting in steeper velocity gradients and thus a slower average flow velocity in the body of the digester. Second, the third dimension is required to model axial rotation produced by an impeller. Third, the 2D plane obtained by cutting through an asymmetric 3D model can vary substantially according to the chosen cutting position, as is the case in the digester because of the inlet and outlet pipes. Thus, the numerical results obtained from the simulation of a custom-made 2D model using the commercial CFD package ANSYS Fluent are used to verify the hydrodynamic features of the Lagrangian model proposed herein. The vertical velocity profile predicted by the time dynamic SPHbased model is compared against the FVM steady state solution shown in Fig. 8 . Single-frame snapshots of the SPH simulation showed high fluctuations in the particlewise velocity values, meaning that a snapshot of the velocity is not representative. We attribute this to the small time step size between each SPH step. Therefore, after letting the simulation run for 120 seconds, the average velocity of 10 0 0 frames at an interval of 0.015 seconds was taken in order to smooth out these fluctuations. The resulting velocity profiles are plotted for a horizontal section of the eggshaped digester 6 m above the digester bottom with a total solid concentration (TS) of 12.1%. It can be observed that the velocity values predicted by the SPH method are in good agreement with the ANSYS Fluent results. These differences can be attributed to the different nature of the simulations, but the existing discrepancies do not influence the overall biomass flow. The ANSYS Fluent results are obtained from a fully converged steady-state solution, while the SPH results are presented for a short time interval. For this reason, the ANSYS Fluent simulation yields a smoother velocity profile.
Contours of the velocity fields are compared in Fig. 9 . A good agreement is observed between the two simulations, especially in the top and bottom areas. The contours of the fluid flow in the middle region of the digester are less pronounced in the SPH simulation, which can be attributed to the fluctuations and the subsequent averaging described above. The recirculation stream between the inlet and the outlet (see Fig. 7 for details of the SPH simulation) is implemented using a periodic boundary condition, which is a well-established treatment in SPH to emulate infinite flows [47] . In the SPH simulations presented herein, particles leave the digester from the outlet, with their velocity being controlled by a backward force field. When they enter the digester from the inlet again, their physical and biochemical properties are set equal to those of newly added biomass.

Mesh convergence study and GCI of the 2D ANSYS model
We generate 2D quadrilateral meshes inside the digester using ANSYS meshing application (version 19.0) with a 4 core Intel Xeon E5-1620 v3 CPU @ 3.50 GHz and 16GB of RAM. By setting the maximum allowed face size of a cell we can control how coarse the resulting grid should be. The iterative solution algorithm converges given an error bound for the residuals and runs in parallel using an NVIDIA GeForce GTX TITAN X. For a description of the three different meshes we used for the refinement study as well as the corresponding computation times see Table 2 .
Given the meshes in Table 2 , we compute the grid convergence index (GCI) according to the procedure presented in [48] . The parameter we are interested in is the pumping rate measured inside the draft tube of the digester denoted by φ j for mesh j ∈ {1, 2, 3}, respectively.     (23) yields an order of convergence of about p = 5 . 020 . The approximate relative error between high-and medium-resolution mesh is given by Together with a safety margin of 25% we compute the GCI via In summary, our results show that the value of the GCI for the pumping flow rate measured inside the draft tube is very low, thus confirming a robust model setup and a well chosen mesh resolution ( Table 3 ).

Full-scale anaerobic digester
To demonstrate the capability of the proposed model for the simulation of full-scale reactors, here we consider the real-life anaerobic digester described in Section 4.3 . In our numerical model, the fluid flow is induced using a force field inside the draft tube to accelerate the biomass to the desired velocity, which is known from the actual operation of the AIZ treatment plant. The no-slip boundary condition is applied on all solid walls of the digester and the draft tube. A periodic boundary condition is applied at the inlet and outlet of the recirculation flow. In order to control the recirculation flow rate a backward force field is applied at the outlet to damp down the high velocity field caused by the hydrostatic pressure at the bottom of the digester. The initial concentration of the biologically active compounds are summarized in Table 4 . These concentrations are obtained after a 4 h simulation of the AD model in a CSTR. Concentrations at the inlet are the same as the initial conditions summarized in Table 1 . Both the fluid phase and the solid walls of the egg-shaped anaerobic digester are discretized by the SPH particles. The particles are initially placed on a regular grid of points with an initial distance of 0.05 m. This configuration leads to a number of 103,868 particles for the whole computational domain (described in Section 2.1 ).
The number of created particles for a 2D SPH simulation of a full-scale digester is noteworthy. With the well known heavy computational demand of Lagrangian methods [49,50] in mind, the simplification of the AD process to two dimensions is nearly inevitable, although this affects turbulence quantities and mixing becomes highly directionally biased in 2D. With the extension of the simulations to three dimensions the number of particles drastically increases, which renders the computational cost prohibitive for real world applications. Therefore, the simulations presented herein are simplified to two dimensions and for this reason, the gpuSPHASE framework which has been developed and optimized Table 4 The initial concentrations (kg.m −3 ) of the biologically active compounds for solving the system of ODEs for the AIZ anaerobic digester. for the simulation of 2D long-running hydraulic phenomena is employed [39] .
Simulations of the AIZ digester were run with six different mixing intensities within the draft tube. We started the numerical simulations with a velocity of 1.5 m . s −1 within the draft tube, which is the value actually used in the AIZ WWTP. In order to investigate the effect of hydrodynamics on the biochemical AD reactions, five more simulations have been carried out with different values for the mixing velocity within the draft tube, namely 0.5, 1.0, 2.0, 2.5 and 3.0 m . s −1 . The mixing velocity at the recirculation flow is chosen equal to 0.95 m . s −1 for all the cases, as actually used in the AIZ WWTP. The methane production values are plotted in Fig. 10 for the six hydrodynamically different configurations. The methane production amounts are calculated by integrating over all the SPH particles; each of which acting as an individual bioreactor and having the mass transfer interactions with the neighboring particles by the diffusion equation. Evidently, the methane production is subject to the mixing velocity, with a defined optimum at about 2 m . s −1 .
A very important advantage of the proposed integrated model is that the direct link between hydrodynamics and the biochemical AD reactions makes it possible to optimize the mixing intensity of the process, thereby enhance the net energy yield. As shown in Fig. 10 , the optimum mixing intensity of the considered egg-shaped digester occurs at a mixing velocity of 2 . 0 m . s −1 within the draft tube. It can be observed that the methane production increases until the optimum point ( 2 . 0 m . s −1 ) and decreases afterwards. The methane yield drops to an almost constant value ( 0 . 631 kg . m −3 ) and does not increase with higher mixing velocities (from 2 . 0 m . s −1 on). The methane yield fall after the optimum mixing intensity is mainly attributed to the excessive mixing impact on the biomass, which as demonstrated in [51] causes the particles leave the digester as residual-undegraded biomass and consequently leads to a lower amount of produced methane. Fig. 11 (a) shows the distribution of methane concentration after 10 0 0 s of physical simulation time with a mixing velocity of 1 . 5 m . s −1 . This spatial illustration makes it possible to observe that the concentration of methane increases with time and reaches to its maximum ( 5 . 171 × 10 −4 kg . m −3 ) at points far from the new biomass entrance. Furthermore, the non-mixed dead zones of the digester are seen at the regions close to the solid walls which are far from the mechanical draft tube as well as the inlet/outlet areas. These observations are insightful to improve and to optimize the mixing methods in AD. The same parameter of the same digester is depicted in Fig. 11 (b) after a longer simulation time (at t = 3600 s ). Evidently, the methane concentration has increased  The distribution of methane concentration after 10 0 0 s of physical simulation time with a mixing intensity of 3 . 0 m . s −1 is shown in Fig. 12 . As it is observed, with the higher mixing velocity the biomass is better mixed within the digester. There are smaller nonmixed areas close to the walls and the diffusion effects are more homogeneous over the whole volume. Furthermore, the higher mixing intensity leads to higher concentration peaks for methane and consequently to higher methane production rates (see Fig. 10 ). The proposed integrated model is able to predict the effects of mixing intensities, which is the most energy consuming part of the process. Being complementary to the quantitative analyses (e.g. in Fig. 10 ), these qualitative analyses lead to optimal design of the AD process.
Another important feature of the AD process is the methane production potential of the biomass through time and space, which is also influenced by hydrodynamics of the process. A thorough investigation of the effect of mixing on the local methane production rates is necessary to improve the net energy yield. Fig. 13 illustrates the pointwise methane production rates throughout the whole volume of the digester. As it is shown, methane production is in its highest rate where the undegraded biomass is entering the digester. This rate decreases with time as the biomass particles move downward and is in its lowest rate at the bottom of the digester. In addition, at regions close to the solid walls this rate is lower than other regions. It is worth noting that optimization studies (see e.g. Fig. 10 ) can be carried out also for the local methane production rates according to different hydrodynamic circumstances, which for the sake of brevity are not reported here.
The pH value is a crucial parameter in the AD process. The acidity of the digester environment directly affects growth and de-  cay of the anaerobic bacteria, and consequently the overall performance of the process. The pH value distribution is depicted in Fig. 14 . In a similar pattern as the methane production rate, the pH value decreases as the particles settle down. As expected, the methane production rate changes faster than the pH value and for this reason, pH values show more uniform distribution. The rate of change of pH values is higher in the regions far from the solid walls.

Conclusions
Whilst many studies investigated digester mixing or biochemical processes of AD, the important link between hydrodynamics and biokinetics has gone largely unaddressed. In this paper a novel model has been developed based on SPH to create a direct link between mixing and biochemical processes. The computational results show that the integrated model makes it possible to directly study the impact of mixing on biogas production and opens future possible developments in AD studies.
It is shown that the biological concentrations can be spatially resolved due to the particle-based disrcretization of the computational domain. The peculiarities of SPH in dealing with multi-phase problems facilitate the description of the biokinetics.
Both SPH and biological models are validated against existing analytical or experimental results, showing reasonable agreements. Through application of the model to a real world egg-shaped anaerobic digester, it has been shown that the model could successfully replicate both mixing and biochemical processes. However, industrial application of the proposed model depends on computational efficiency of the future SPH solvers. Concerning the next developments, this work is a stepping stone to develop a fully integrated model, which couples the mixing and biochemical processes of AD and encompasses the sedimentation effects, non-Newtonian flows and heat transfer. In the long term, it is envisaged that the presented Lagrangian assets of SPH can improve the computational investigations of AD.