# Modelling of biological methanation by CFD and 1D model in an industrial bubble column

Vincent Ngu, Jérôme Morchain, Benjamin Casale, Arnaud Cockx

► To cite this version:

Vincent Ngu, Jérôme Morchain, Benjamin Casale, Arnaud Cockx. Modelling of biological methanation by CFD and 1D model in an industrial bubble column. Chemical Engineering Science, 2025, 316, pp.121940. 10.1016/j.ces.2025.121940. hal-05371717

HAL Id: hal-05371717

https://insa-toulouse.hal.science/hal-05371717v1

Submitted on 18 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modelling of biological methanation by CFD and 1D model in an industrial bubble column*

Vincent Ngu ID, Jérôme Morchain *ID, Benjamin Casale, Arnaud Cockx* ID

Toulouse Biotechnology Institute, Université de Toulouse, CNRS, INRAE, INSA, Toulouse, France

## ARTICLE INFO

### Keywords:

- Bubble column
- Multiscale simulations
- Multiphase reactors
- Gas-liquid mass transfer
- Hydrogen
- CFD

## ABSTRACT

In this paper, a Euler-Euler CFD model based on the two-fluid model is developed for the industrial-scale biological methanation process. The model resolves the hydrodynamics of the bubbly flow, along with the effects of the multispecies mass transfer (H₂, CO₂ and CH₄) and bioreaction kinetics. Three-dimensional simulations were performed with the commercial CFD code ANSYS Fluent® version 2021 R1 and compared with one-dimensional simulations solved with a previously validated in-house 1D code Matlab®. Both modelling tools implemented closure relations in a consistent manner, ensuring proper coupling in terms of interfacial force interactions, mass transfer, and bioreaction uptake rate. Following validation on a pilot methanation bubble column of 1 m high, the results focus on a 10 m high industrial bubble column where the hydrogen gas-liquid transfer rate is a limiting factor of the methanation process.

At this scale, hydrostatic pressure effects and longer gas residence times strongly intensify mass transfer and reduce gas hold-up and bubble size, leading to inverse density gradients and instability of the bubbly flow. Due to the substantial gas consumption, complex hydrodynamics prevail in large reactive bubble columns and significant concentration heterogeneities appear. The application of unsteady 3D two-phase flow simulations provides a more comprehensive description of these phenomena than reduced 1D models. Nevertheless, the 1D model remains a valuable tool for predicting overall performance, while its parameters can be adjusted to reflect hydrodynamic instabilities and to facilitate rapid investigation of operating scenarios.

## 1. Introduction

Bubble columns are multiphase contactors in which coupled multi-physical phenomena occur simultaneously at different spatial and temporal scales. Depending on the problem and the scale to be addressed, numerical simulations of bubble columns can be performed using different modelling approaches. Among these, multiphase flow simulations (using Computational Fluid Dynamic tools) is now routinely used but it remains a challenging task when a strong coupling exists between hydrodynamics, mass transfer and reaction. In contrast, 1D model are easier to handle and remain of valuable interest to perform sensitivity analysis on operating conditions and/or (bio) reaction kinetics (Siebler et al., 2020; Ngu et al., 2023).

Successful scale-up of biological methanation reactors requires a deep understanding of the mass transfer efficiency, especially that of hydrogen, to avoid poor methanation performance and high operating costs. The modelling of large-scale bubble column bioreactors requires a

detailed mechanistic description of the bubbly flow and closures on the gas-liquid interfacial mass transfer as well as the biological kinetics. In particular, the stoichiometry of that reaction is such that the gas fraction is highly dependent on the bioreaction rate with major consequences on two-phase flow. This work focuses on this two-way coupling through unsteady Euler-Euler CFD simulations of an existing large-scale biological methanation plant operated by Electrochaea GmbH, Denmark. This unit was previously studied by Ngu et al. (2022) using a 1D spatiotemporal model, which allowed a transient description of the route towards steady state for that biological process. Nevertheless, CFD models provide more accurate spatially resolved local interactions between phases which is valuable for studying large-scale multiphase flow, mass transfer, and bioreaction coupled problems. Despite the higher computational cost, CFD models are reliable conceptual reactor design tools for industrial-scale bubble column bioreactors, as demonstrated, for examples, in recent studies by Puiman et al. (2022), Siebler et al. (2019 and 2020), Nadal-Rey et al. (2022) and Mast et al. (2024).

CFD simulations of biological methanation with the validation of the

<footnote>* This article is part of a special issue entitled: 'GLS-16' published in Chemical Engineering Science.</footnote>

<footnote>* Corresponding authors.
E-mail addresses: morchain@insa-toulouse.fr (J. Morchain), cockx@insa-toulouse.fr (A. Cockx).</footnote>
Nomenclature

Roman symbols

<table><thead><tr><th>Symbol</th><th>Description (Unit)</th></tr></thead><tbody><tr><td>a</td><td>Volumetric interfacial area (m<sup>2</sup>·m<sup>-3</sup>)</td></tr><tr><td>C<sup>*</sup></td><td>Solubility/equilibrium concentration (mg·L<sup>-1</sup>)</td></tr><tr><td>C<sub>G</sub></td><td>Gas concentration (mg·L<sup>-1</sup>)</td></tr><tr><td>C<sub>L</sub></td><td>Liquid concentration (mg·L<sup>-1</sup>)</td></tr><tr><td>D<sub>i</sub></td><td>Species diffusivity (m<sup>2</sup>·s<sup>-1</sup>)</td></tr><tr><td>d<sub>b</sub></td><td>Sauter Bubble diameter (m)</td></tr><tr><td>D<sub>G</sub></td><td>Gaseous axial dispersion coefficient (m<sup>2</sup>·s<sup>-1</sup>)</td></tr><tr><td>D<sub>L</sub></td><td>Liquid axial dispersion coefficient (m<sup>2</sup>·s<sup>-1</sup>)</td></tr><tr><td>D<sub>r</sub></td><td>Reactor diameter (m)</td></tr><tr><td>H</td><td>Reactor Height (m)</td></tr><tr><td>j<sub>G</sub></td><td>Superficial gas velocity (m·s<sup>-1</sup>)</td></tr><tr><td>j<sub>L</sub></td><td>Superficial liquid velocity (m·s<sup>-1</sup>)</td></tr><tr><td>k<sub>L</sub></td><td>Liquid side mass transfer coefficient (m·s<sup>-1</sup>)</td></tr><tr><td>k<sub>L</sub>a</td><td>Volumetric mass transfer coefficient (s<sup>-1</sup>)</td></tr><tr><td>M</td><td>Molar mass (g·mol<sup>-1</sup>)</td></tr><tr><td>ṁ<sub>G</sub></td><td>Interfacial mass transfer rate (mg·L<sup>-1</sup>·h<sup>-1</sup>)</td></tr><tr><td>P</td><td>Pressure (Pa)</td></tr><tr><td>S<sub>CB</sub></td><td>Coalescence and breakage in the bubble population (m<sup>-1</sup>·s<sup>-1</sup>)</td></tr></tbody></table>

<table><tr><td>t</td><td>Time (s)</td></tr><tr><td>u<sub>rel</sub></td><td>Relative velocity between phases (m·s<sup>-1</sup>)</td></tr><tr><td>u</td><td>Local phase velocity (m·s<sup>-1</sup>)</td></tr><tr><td>X</td><td>Biomass concentration (mg·L<sup>-1</sup>)</td></tr><tr><td>y</td><td>Gas molar fraction (–)</td></tr><tr><td>z</td><td>Axial coordinate (m)</td></tr></table>

Greek symbols

<table><tr><td>a</td><td>Local phase holdup (–)</td></tr><tr><td>ε</td><td>Cross-section averaged phase holdup (–)</td></tr><tr><td>ρ</td><td>Fluid density (kg·m<sup>-3</sup>)</td></tr><tr><td>φ<sub>i</sub></td><td>Specific species uptake rate (g<sub>i</sub>·g<sub>X</sub><sup>-1</sup>·h<sup>-1</sup>)</td></tr></table>

Subscripts

<table><tr><td>G</td><td>Gaseous phase</td></tr><tr><td>i</td><td>Species</td></tr><tr><td>L</td><td>Liquid phase</td></tr></table>

Abbreviations

<table><tr><td>CFD</td><td>Computational Fluid Dynamics</td></tr><tr><td>IAC</td><td>Interfacial Area Concentration</td></tr><tr><td>RANS</td><td>Reynolds Averaged Navier Stokes</td></tr><tr><td>RNG</td><td>Renormalization Group</td></tr><tr><td>UDF</td><td>User Defined Function</td></tr></table>

bioreaction remain scarce in the literature. To the best of our knowledge, the only work with the full validation of the CFD model for biological methanation focuses on a laboratory-scale trickle bed reactor (Markthaler et al., 2022; Markthaler et al., 2020). Furthermore, very few large-scale biological methanation plants are reported in the literature, which makes CFD validation more risky in the absence of local measurements (Rafrafi et al., 2020). However, some open literature data (“Electrochaea - Carbon and Energy storage,” 2017; Jensen et al., 2021; Lardon et al., 2018; Rusmanis et al., 2019; Sveinbjörnsson and Münster, 2017) provide sufficient information for global validation. As an industrial example, it is mentioned that the Electrochaea biological methanation plant in Avedore, Denmark is the first and largest full-scale demonstration plant in the world (Rusmanis et al., 2019). The aim of this paper is to perform CFD simulations of the Electrochaea’s large-scale biological methanation reactor and to compare the results with 1D model and experimental data from the literature. Similar to Siebler et al. (2020), but in a different type of bioreactor, this paper “deals with the trade-off between sophisticated spatio-temporally resolved large-scale simulations using computationally intensive Euler-Euler approaches (CFD) and coarse-grained 1D models enabling fast performances evaluations”.

Both modeling approaches consider two-phase flow, multispecies mass transfer (H₂, CO₂ and CH₄) and their coupling with the bioreaction kinetics. Where relevant, comparisons of the CFD results with those of the 1D model are presented and discussed to shed light on the unexpectedly complex hydrodynamic behaviour of that bioreactor.

In this paper, after presenting the industrial context of biological methanation, the CFD model setup and interfacial force interactions, mass transfer and biological uptake rate implementations are described. The two-phase 1D model uses the same closure relations (parts 1.1 and 1.2). After a preliminary validation of the implementation (part 2.1), transient 3D CFD simulations of the industrial bubble column are presented and time-averaged for the nominal methane production (parts 2.2 and 2.3). The comparison with 3D and 1D simulations is then detailed (part 2.4). Finally, the 1D model is used to show the sensitivity of the methane purity to the gas and liquid dispersion.

## 2. Experimental and numerical methods

### 2.1. Biological methanation at industrial scale

The industrial bubble column considered is designed to treat up to 50 Nm³/h of CO₂ from a nearby anaerobic digester with a typical biogas mixture of 60 % CH₄ and 40 % CO₂. With the CO₂ from the biogas, the bioreactor is expected to produce 125 Nm³/h of CH₄ with a purity between 95 and 97 % molar fraction in order to feed the natural gas network according to the standard required in European countries (Thema et al., 2019). This requires the addition of 200 Nm³/h of H₂ produced from the electrolysis of water using surplus renewable electricity. The composition of the gas at the inlet matches then the following stoichiometry expected in the biological reaction:

$$4 H_2 + CO_2 \rightarrow CH_4 + 2 H_2O$$

The industrial bubble column operates under thermophilic conditions at a temperature of 65 °C and a pressure of 8 bar.

The dimensions of the large-scale bubble column, in Fig. 1, were estimated from publicly available data and images (“Electrochaea -

<figure><img src="image_2.png" /><figcaption>Fig. 1. Industrial site with biomethanation column in Denmark (Credit: www.electrochaea.com).</figcaption></figure>
Carbon and Energy storage,” 2017; Lardon et al., 2018). The initial liquid height is assumed to be 10 m but, to account for the increase in water level due to gassing, the total height of the column is set to 11 m. The internal diameter is set at 0.72 m according to Rusmanis et al. (2019).

## 2.2. CFD mesh and model setup

The 3D geometry was constructed with a gas sparger homogeneously placed on the reactor base. The mesh was generated using a sweeping technique throughout the column with a meshing strategy similar to that of Siebler et al. (2019). A mesh consisting of 141 000 hexahedral cells with an average size of 5 cm, good orthogonal quality and skewness was constructed, as shown in Fig. 2.

The transient multiphase simulation was performed in the Eulerian framework using ANSYS Fluent 2021R1. The CFD model setup of the large bubble column was largely based on the work of Siebler et al. (2019). The Reynolds-Averaged Navier-Stokes (RANS) equations are solved with dispersed RNG $k-\epsilon$ turbulence model. Both phases were set as mixtures to investigate the concentration of each species H₂, CO₂, CH₄.

The interfacial forces modelling is identical to that used in the large-scale bioreactor simulation work of Puiman et al. (2022) and Siebler et al. (2019) which considered the drag, lift, wall lubrication, turbulent dispersion, and turbulence interaction forces (Table 1). The inlet bubble diameter at inlet is assumed to be 2 mm as in the work with the 1D model of Ngu et al. (2022). It was found that small bubbles (≤ 2 mm) are necessary to ensure high H₂ mass transfer flux and to achieve high CH₄ purity at the industrial scale of 10 m height and under 8 bars pressure (Table 2).

A superficial gas velocity of 0.031 m.s⁻¹, corresponding to the nominal industrial gas flow rate, was set at the inlet sparger with the velocity magnitude equal to the bubble terminal velocity calculated using the Tomiyama et al. (1998) drag model (0.286 m.s⁻¹). A pressure outlet with gas backflow was imposed at the top of the column, with the headspace pressure set at 8 bar. All walls were set to a no-slip condition for the liquid and a free-slip condition for the gas. The inlet gas species mole fraction ratio of H₂:CO₂:CH₄ was 61.5 %:15.3 %:23.2 % (Table 2). The ideal gas law was applied to the multi-species gas phase.

The transient resolution strategy is detailed as follows: first, the unsteady gas-liquid flow is solved without considering mass transfer and bioreaction to obtain a first indication of the flow structure; followed by the species transport equations. Lastly, mass transfer and bioreaction are activated, and all equations were solved simultaneously. A total of 500 s of physical time is required to achieve a pseudo-stationary regime. Finally, data sampling at 1 Hz was performed during an additional 100 s period in order to generate time-averaged quantities that could be compared to experiments and 1D model. All simulations were performed with a time step of 0.01 s and a maximum number of 10 iterations per time step.

<figure><img src="image_3.png" /><figcaption>Fig. 2. (a) Axial and (b) radial section of the mesh.</figcaption></figure>

<caption>**Table 1**</caption>

Summary of interfacial force correlations used in the CFD model for industrial case.

<table><caption>**Table 2**</caption><thead><tr><td>Force</td><td>Correlation</td></tr></thead><tbody><tr><td>Drag</td><td>(Tomiyama et al., 1998)</td></tr><tr><td>Lift</td><td>(Tomiyama et al., 2002)</td></tr><tr><td>Wall</td><td>(Antal et al., 1991)</td></tr><tr><td>Turbulent dispersion</td><td>(Burns et al., 2004)</td></tr><tr><td>Turbulent interaction</td><td>(Sato and Sekoguchi, 1975)</td></tr></tbody></table>

Physical properties and parameters used in the CFD model.

<table><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>Reactor diameter</td><td>0.72</td><td>m</td></tr><tr><td>Reactor height</td><td>11</td><td>m</td></tr><tr><td>Initial liquid height</td><td>10</td><td>m</td></tr><tr><td>Temperature</td><td>65</td><td>°C</td></tr><tr><td>Headspace pressure</td><td>8</td><td>Bar</td></tr><tr><td>Inlet superficial gas velocity</td><td>0.031</td><td>m.s<sup>-1</sup></td></tr><tr><td>Inlet molar fraction H<sub>2</sub>: CO<sub>2</sub>: CH<sub>4</sub></td><td>61.5: 15.3: 23.2</td><td>%</td></tr></tbody></table>

In line with the outcome of the 1D model established in Ngu et al. (2022), the bubble size variation is connected to (i) the hydrostatic pressure changes and (ii) the huge mass loss in the gas phase (4 mol of H₂ and 1 mol of CO₂ reduced to 1 mol of CH₄). To account for these effects, the interfacial area transport equation was activated in ANSYS Fluent. It is an additional scalar quantity called the interfacial area concentration (IAC) for which equation (1) is solved.

$$ \frac{\partial(\rho_G a)}{\partial t} + \nabla \cdot (\rho_G \overline{u_G} a) = \frac{1}{3} \frac{D \rho_G}{D t} a + \frac{2}{3} \frac{\dot{m}_G}{\alpha_G} a + \rho_G (S_{CB}) \quad (1) $$

In this equation, *a* is the volumetric interfacial area, $\dot{m}_G$ the interfacial mass transfer rate and $S_{CB}$ a term related to coalescence and breakage of bubbles. With this model, the change in bubble size due to volume contraction or expansion, mass transfer, bubble coalescence and breakup is translated into changes in the interfacial area. Similar to the spatio-temporal 1D model, coalescence and break-up phenomena are assumed to be balanced, so that the bubble size variation is essentially controlled by mass transfer and pressure effects. The third term in the left-hand side of (1) was therefore discarded.

Two user-defined functions (UDFs) for mass transfer and bioreaction uptake rate have been developed. The species transfer rate $\dot{m}_i$ is given by equation (2), which is a function of the species dependent mass transfer coefficient $k_{L,i}$, the interfacial area *a* and the driving force ($C_i^* - C_{L,i}$). The mass transfer coefficient is modelled with the Higbie (1935) model given as:

$$ \dot{m}_i = k_{L,i} a (C_i^* - C_{L,i}) \quad (2) $$

$$ k_{L,i} = 2 \sqrt{\frac{u_{rel} D_i}{\pi d_b}} \quad (3) $$

$$ d_b = \frac{6\alpha_G}{a} \quad (4) $$

With $d_b$ is the Sauter mean diameter calculated from the resolved interfacial area transport equation. The diffusion coefficients $D_i$ for the three species H₂, CO₂, and CH₄ at 65 °C are 1.2 × 10⁻⁸, 5.4 × 10⁻⁹, and 3.95 × 10⁻⁸ m².s⁻¹, respectively. The dimensionless Henry's constant at 65 °C was used to calculate the solubility $C_i^*$ with the corresponding dimensionless Henry's constant for the species H₂, CO₂, and CH₄ being 56.27, 2.79, and 54, respectively (Sander, 2015).

For biological methanation, it is known that the H₂ gas-liquid mass transfer rate is the limiting factor of the overall biological methanation process (Jensen et al., 2021; Lecker et al., 2017; Rafrafi et al., 2020). It is
then assumed that there is no enhancement of the physical mass transfer by the biological reaction as confirmed in the review by Volger et al. (2024) on the effect of broth composition on mass transfer in bubble columns. It has been proposed that the steady state is characterised by a balance between the specific hydrogen uptake rate (mass transfer rate divided by the biomass concentration) and the utilisation rate for maintenance (Ngu et al., 2022). From this point, constant rates of carbon dioxide consumption and methane production are achieved, both being proportional to the H₂ mass transfer rate. Since the hydrogen mass transfer rate is the limiting factor for the metabolism, the specific hydrogen uptake rate at steady state is equal to the maximum hydrogen mass transfer rate per cell mass:

$$ \varphi_{H2,max} = \frac{k_L a(C_{H_2}^* - 0)}{a_L X} \qquad (5) $$

The specific carbon dioxide and methane uptake rates are expressed stoichiometrically as

$$ \varphi_{CO_2} = \frac{\varphi_{H2,max}}{4} \cdot \frac{M_{CO_2}}{M_{H_2}} \quad \text{and} \quad \varphi_{CH_4} = \frac{\varphi_{H2,max}}{4} \cdot \frac{M_{CH_4}}{M_{H_2}} \qquad (6,7) $$

Since the aim here is to predict the bioreactor performance at steady state and to compare with the literature value, a constant biomass concentration, X = 35 g.L⁻¹ is considered in line with the results obtained with the 1D model for the industrial case (Ngu et al., 2022). As the growth rate of methanogens is usually very low, it is not possible and therefore not attempted to simulate the transient growth using the CFD. This approach is similar to the work of Mast et al. (2024) where the biomass concentration is estimated by a one compartment model to define the operating parameters used in the CFD at several quasi-steady state fermentation times.

# 3. Results and discussion

## 3.1. Preliminary test case on a small methanation bubble column

Prior to large-scale biological methanation CFD simulations, we performed the CFD workflow in a 1 m high reactive bubble column. The test case was a gas-liquid co-current bubble column (gas and liquid

superficial gas velocities of 5 and 1 mm.s⁻¹, respectively) in a 2D geometry. Fig. 3(a) shows the mesh of the simplified test case, which is presented in 2D with symmetry boundary conditions on the vertical sides to mimic a 1D system. Here, the validation case consists of the hydrogen and carbon dioxide mass transfer in co-current bubbly flow. In Fig. 3(b), the CFD results are compared with those of the 1D model of Ngu (2022) in order to validate the reliability of the CFD setup (in particular, user subroutine implementation). Note that the range has been chosen to emphasize the differences between the two approaches while the relative error is below 3 % in the worst case (gas retention). Fig. 3(b) and (d) show the hydrogen concentration evolution in the gas and liquid phases, respectively. In Fig. 3(b), a 10 % increase in the hydrogen gas concentration (from 0.035 to 0.038) is observed in the first half of the column due to the large loss of carbon dioxide in the gas phase. The liquid hydrogen concentration further increases due to mass transfer, without reaching its equilibrium concentration because of an insufficient contact time with the liquid phase. Fig. 3(c) and (e) show that the bubble diameter and gas holdup decreased in the lower half of the column due to the high mass transfer flux, but increased slightly in the upper half of the column due to the volume expansion caused by the decreasing hydrostatic pressure. Thanks to the Interfacial Area Concentration (IAC) equation, it is shown that when mass transfer dominates in Eq. (1), both the gas holdup and the bubble diameter are reduced, since the shrinkage of the bubble size follows the cube root of the gas volume (Bauer and Eigenberger 1999; Ngu et al. 2022). On the contrary, the pressure evolution tends to increase the gas holdup along the height of the bubble column, which is observed more in the second half of the bubble column. At this stage, it can be concluded that the two-way coupling between multispecies transfer and two-phase hydrodynamics is validated by the good agreement with the in-house 1D Matlab® equivalent model. The same setup could then be used on an industrial scale in a 3D cylindrical bubble column.

## 3.2. Pseudo-stationary hydrodynamics at industrial scale

The unsteady gas-liquid flow is first solved without considering mass transfer and bioreaction, but with passive species transport equations to get a first indication of the flow structure. Fig. 4 shows the gas-liquid flow and the species contours after 150 s. The hydrodynamic results

<figure><img src="image_4.png" /><figcaption>Fig. 3. (a) Rectangular 2D mesh for the test case. Comparison of CFD (red circles) and 1D (continuous line) predictions for co-current mass transfer in a 1-m high bubble column. (b) Hydrogen gas concentration (c) bubble diameter, (d) hydrogen liquid concentration (e) gas holdup.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Contours of (a) axial gas velocity component (b) axial liquid velocity component (c) gas holdup (d) gas species molar fraction.</figcaption></figure>

show that the bubbles rise rather uniformly due to a homogeneous gas injection covering the entire column cross-section, with a bubble terminal velocity of about 0.3 m/s according to the Tomiyama et al. (1998) drag model. The gas holdup increases slightly and linearly with decreasing hydrostatic pressure (9 > P > 8 bars), in accordance with the ideal gas law. The total gas holdup was 11.7 %, which is in good agreement with the values obtained from the 1D model (11.6 % in Ngu et al., 2022). Following the work of Puiman et al. (2022), who also compared the averaged gas holdup with the global correlation of Zuber and Findlay (1965) and Heijnen and Van't Riet (1984), it is found that the correlations give 10.9 % and 12.4 % respectively, which is close to the CFD predicted value. The fluid flow is uniform with very weak recirculation, resulting a homogeneous bubbly regime where the gas fraction evolution is only due to hydrostatic pressure (Fig. 4). The gas species are also homogeneously distributed in the column.

The multispecies mass transfer and bioreaction kinetics were then activated to interact with the hydrodynamics for 500 s of simulation time.

### 3.3. Hydrodynamics coupled to transfer and bioreaction

In the presence of bioreaction and mass transfer, the hydrodynamic behaviour is drastically different from the previous case. In this section, the results of the large-scale biological methanation are discussed. First, the predictions of the hydrodynamic variables such as bubble diameter, gas holdup, gas velocity and interfacial area are presented, as shown in Fig. 5. In each subfigure, an instantaneous contour map is presented on the left and a time-averaged contour map over 100 s is shown on the right.

Due to huge transfer and consumption of H₂ and CO₂, the gas holdup varies significantly from the inlet to the outlet, ranging from about more than 10 % at the bottom to 5 % at the top, as shown in Fig. 5(c). This almost two-fold reduction in gas holdup creates an axial gradient of gas holdup, resulting in a heavier two-fluid mixture at the top of the column compared to the bottom. This density gradient is responsible for the large unsteadiness of the two-phase flow. As a consequence, large-scale mixing is favoured by meso-scale structures due to local gas accumulation, resulting in axial and radial dispersions.

The bubble diameter varies also spatially and temporally due to these density gradient effects. The time-averaged data show that the mean bubble size decreases axially from the bottom to the top of the column, as shown in Fig. 5(a), from 2 mm to about 1.5 mm. This is less than the

<figure><img src="image_6.png" /><figcaption>Fig. 5. Contour plots of instantaneous and time-averaged profiles on a vertical plane and selected cross-sections (every meter) for (a) bubble diameter (b) axial gas velocity component (c) gas holdup and (d) volumetric interfacial area.</figcaption></figure>

decrease in gas hold-up, as the shrinkage of the bubble size follows the cubic root of the gas volume (Bauer and Eigenberger 1999; Ngu et al. 2022).

The gradient profiles of gas hold-up and bubble size result in a gradient of the interfacial area (Equation (4)) with maximum values also located at the bottom in Fig. 5(d). The combination of higher interfacial area and hydrostatic pressure is then expected to result in a high mass transfer rate in this region. At the top of the column, the interfacial area for mass transfer is lower, but still meets the biological demand enabling biomethane production. The spatially dispersed axial gas velocity in Fig. 5(b) reflects the influence of the different bubble sizes in Fig. 5(a) on the gas-liquid hydrodynamics, as different bubble sizes rise at a different velocity, but overall, the influence is that of the unsteady axial liquid velocity. It can also be observed that the time-averaged cross-section of the gas holdup and the interfacial area were mostly homogeneous from the mid-height up to the top of the column.

Complementing the previous contour maps, Fig. 6 shows both the 10 s time step instantaneous and 100 s time-averaged radial profiles of the liquid axial velocity and gas holdup for z/H of 0.25 (top row), 0.5
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Radial profiles of (a, c, and e) liquid axial velocity and (b, d, and f) gas holdup for z/H of 0.25 (top row), 0.5 (middle row) and 0.75 (bottom row). Red circle lines are instantaneous profiles, Black triangle lines are time-averaged profiles.</figcaption></figure>

(middle row), and 0.75 (bottom row). For the liquid axial velocity, it fluctuates approximately between ± 0.30 m.s⁻¹ at z/H of 0.25, between ± 0.20 m.s⁻¹ at z/H of 0.5, and between ± 0.10 m.s⁻¹ at z/H of 0.75. At each height, all instantaneous radial profiles of liquid velocity are asymmetric, reflecting the unsteady bubbly flow in the bubble column. At z/H of 0.25, the liquid mainly ascends on the right side and descends on the left side but at some instant, inversion can also occur, as shown in Fig. 6 (a). As z/H increases, the liquid velocity is found to decrease. This shows that the two-phase flow becomes more stable with weaker liquid recirculation in the column. At the top of the column, the axial liquid velocity becomes negligible. In addition, as z/H increases, the time-averaged radial profiles (black triangular lines) become more uniform with less radial gradients. Meanwhile, the gas holdup varies approximately between 0.06 and 0.07 at z/H of 0.25, between 0.05 and 0.06 at z/H of 0.5, and between 0.04 and 0.05 at z/H of 0.75. For the gas holdup, the radial profiles are rather uniform. At z/H = 0.5 and 0.75, three radial profiles show a slightly higher values than other profiles, reflecting the time-dependent behaviour of the column. This is also observed in the time profile of the gas holdup in Fig. 7 (c).

Fig. 7 shows the time profiles for different variables at three different z-locations. The time profiles were recorded for 100 s at the axis position (x = 0, y = 0) for three different heights. The Sauter mean diameter and the gas holdup show a slight increase in value at less than 20 s, but decrease immediately afterwards. The gas holdup remains rather stable after 40 s, indicating that the pseudo-steady state has been reached. At these measuring points, the axial liquid velocity fluctuates between ± 0.05 m.s⁻¹. The gas species molar fraction of the substrates (hydrogen and carbon dioxide) decreases with z, while the methane gas molar fraction increases, reflecting the successful conversion of the substrates

into the gas product. For carbon dioxide, the substrate is mostly depleted after the mid-height.

Fig. 8 shows the molar fraction surface plots for hydrogen and methane. Carbon dioxide is not plotted as it is completely transferred in the first meter of the column due to its high solubility and the high gas-liquid mass transfer rate. Fig. 8(a) shows that the hydrogen molar fraction decreases from about 60 % at the inlet to less than 10 % at the outlet. Successful biogas upgrading was achieved as the methane molar fraction increased from 25 % at the inlet to 95 % at the outlet, as shown in Fig. 8(b) It is close to the purity reported at full industrial scale (between 95 and 97 %). This proves that the chosen bioreactor design as a “slim and high bubble column” is suitable to produce high quality methane by promoting H₂ gas-liquid mass transfer. High bubble columns are advantageous for mass transfer due to the increased gas residence time and gas solubility (hydrostatic pressure), which has been previously discussed in the literature (Deckwer et al., 1978; Ngu et al., 2022; Siebler et al., 2019). Furthermore, the CFD results showed that the radial gradients are negligible compared to the axial gradients for each variable. This is an indication that a 1D approach may be relevant for comparison.

### 3.4. CFD comparison with the 1D comprehensive model

In this part, some time-averaged CFD results are compared with the previously published 1D results, as shown in Fig. 9. Globally, the CFD model follows the same trend as the 1D model. The hydrostatic pressure decreases linearly from bottom to top, as shown in Fig. 9 (a). Along with the decreasing hydrostatic pressure and gas-liquid transfer, the gas holdup (Fig. 9 (c)) and the volumetric gas flow rate (Fig. 9 (d)) decrease
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Temporal profiles of (a) Sauter mean diameter, (b) hydrogen gas molar fraction, (c) gas holdup, (d) carbon dioxide gas molar fraction, (e) liquid axial velocity, and (f) methane gas molar fraction for a centred point at $z/H = 0.25$ (red points), 0.5 (blue points) and 0.75 (black points).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Contour plots of instantaneous and time-averaged profiles for (a) hydrogen gaseous molar fraction (b) methane gaseous molar fraction.</figcaption></figure>

sharply from the bottom to the top of the column. As a result of the decreasing gas holdup, the gas flow rate was also reduced by almost a factor of two. According to the stoichiometry of the reaction, the flow

rate should be reduced fivefold but since the hydrostatic pressure varies by a factor of two in a 10 m column, the theoretical flow rate reduction is 2.5, which is consistent with the twofold reduction. The differences between the CFD and 1D models in terms of bubble diameter in Fig. 9(b) may appear large, but this is due to the scale chosen for the x-axis, as the volume-averaged relative error in bubble size is only 1.25 %.

Fig. 10 shows the axial gas and liquid concentration profiles for the hydrogen (Figures (a) and (b)) and for the methane (Figures (c) and (d)) obtained with the 3D and 1D approaches. It can be seen that the hydrogen gas concentration decreases throughout the reactor as all the transferred hydrogen is converted to methane. In the liquid, the concentrations of hydrogen and carbon dioxide (not shown) reach zero, indicating the total consumption of the dissolved substrates by the bioreaction.

For the methane concentrations, the differences between the CFD model and the 1D model were the greatest in the lower part of the column. This could be due to the unsteadiness of the two-phase flow in the large bubble column with the presence of bioreaction. Due to the nature of the bioreaction, a large loss of gas was expected, contributing to the gas holdup axial gradients. The spatial evolution leads to large scale unsteady liquid recirculation and bubble plumes due to density gradients, as shown in Fig. 11 (a) at the same time as in Fig. 5. Fig. 11(b) (c) and (d) illustrate the unsteady, oscillating behaviour of the bubble column. Near the inlet, the liquid recirculates downwards near the right wall and pushes the sparged gas to the left side, resulting in a higher local liquid velocity on the left side, thus forming a circulation loop. Axial liquid velocities between –0.2 and + 0.2 m.s⁻¹ can be observed locally with asymmetric liquid recirculation profiles. The gas holdup
<figure><img src="image_10.png" /><figcaption>Fig. 9. Comparison between CFD and 1D axial profiles. (a) Hydrostatic pressure (b) Bubble diameter (c) Gas holdup (d) Gas flowrate.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Comparison between CFD and 1D axial profiles. (a) Hydrogen gas concentration (b) Hydrogen liquid concentration (c) Methane gas concentration (d) Methane liquid concentration (e) Gas molar fraction.</figcaption></figure>

radial profiles were fairly uniform with a slightly higher radial gradient observed at the mid column height ($z/H = 0.5$).

Overall, the CFD predictions are in good agreement with the 1D model results above mid-height. This is expected as the radial gradients after this height were relatively weak, as shown in the radial contour plots in Fig. 5 and Fig. 6. It can also be observed in Fig. 5 that the axial evolution of the gas holdup was weaker at the half top of the column, resulting in a relative stability of the two-phase flow. The slope of the gas

holdup at the half bottom of the column is steeper, reflecting a higher axial gradient ($\left. \frac{\partial \epsilon_G}{\partial z} \right|_{6-10m} \approx 0.2\%/m < \left. \frac{\partial \epsilon_G}{\partial z} \right|_{0-5m} \approx 1\%/m$).

This justifies the choice of the 1D model to simulate the industrial bubble column of biomethanation as it is sufficient to obtain satisfactory results. At the outlet, the CFD model predicts a CH₄ purity of 95 % and a productivity of 126.5 Nm3/h which is in good agreement with the 1D model prediction of 95 % and 123 Nm3/h (Ngu et al., 2022) and the
<figure><img src="image_12.png" /><figcaption>Fig. 11. (a) Contour of liquid axial velocity superposed with the arrows indicating the velocity vectors. (b) Zoom of the liquid velocity vector near the inlet circled with red contours in (a). Radial profiles of (c) liquid axial velocity and (d) gas holdup for z/H = 0.25, 0.5, and 0.75.</figcaption></figure>

reported value of 90–95 % and 125 Nm3/h of CH₄ (“Electrochaea - Carbon and Energy storage,” 2017; Sveinbjörnsson and Münster, 2017). The similarity in the results essentially confirms that the implementation of the different models and the selection of the closure law are correct and physically consistent. The unsteadiness of the flow simulated in the 3D CFD study is modelled as a dispersion term in the 1D model. The usual correlations for estimating this dispersion coefficient are found to be relevant in the configuration simulated here.

In the publication of Ngu et al. (2022), the axial dispersion coefficients for the gas and liquid phases used in the 1D model were estimated from the correlation of Wachi and Nojima (1990) and Deckwer et al. (1974), respectively, but kept constant for the sake of simplicity at values equal to D<sub>G</sub> = 0.05 m<sup>2</sup>/s and D<sub>L</sub> = 0.1 m<sup>2</sup>/s, corresponding to reference values based on the inlet gas velocity. Nonetheless, it could be adjusted as the bioreactor is found to be particularly sensitive to density gradients effects.

Fig. 12 shows that methane purity is more sensitive to the dispersion in the gas phase than in the liquid phase. In this particular case of biological methanation, the substrate consumption induces very low concentrations in the liquid phase and then the dispersion in this phase is not really important because we are close to a perfectly mixed bubble column with regard to biokinetics. On the contrary, the produced

methane flux and its purity in the gas are directly related to the hydrogen flux $\phi_{CH_4} = \frac{q_{H_2max}}{4} \cdot \frac{M_{CH_4}}{M_{H_2}}$ (Ngu et al. 2023). In this large column, the transfer from the bubbles is more sensitive to the gas dispersion, usually close to a plug flow behaviour but here perturbed by a huge dispersion due to the density gradient of the equivalent mixture phase.

## 4. Conclusion

Fully resolved CFD simulations of biological methanation in a large-scale bubble column were performed and they were compared with open literature data for the nominal industrial gas flow rate. To the best of our knowledge, this was the first CFD attempt on a large-scale ex situ biological methanation reactor.

The bioreactor performance predicted by the CFD model was also validated by comparison with that of the 1D model and the results were satisfactory. It was found that when the radial profiles in the second half of the column were fairly uniform, the CFD model gave a very close value to the 1D model. However, due to the nature of the reaction (5 mol of gas substrates being reduced to 1 mol of gas product), large gas holdup depletions are expected in the axial direction. These density gradients lead to unsteadiness in the bubbly flow, which is more

<figure><img src="image_13.png" /><figcaption>Fig. 12. Comparison between CFD and 1D axial profiles for the methane concentration (a) Sensitivity to gas (left) and (b) liquid dispersion coefficients (right).</figcaption></figure>
pronounced in the lower half of the column where carbon dioxide is completely transferred in the first meter.

CFD modelling becomes crucial and advantageous when one is interested in quantifying the local phenomena that are difficult to model using a coarse 1D approach. In this case, the CFD model is the ideal tool to characterise the instabilities and the effects of spatio-temporal heterogeneities on the bioreaction. This study highlights that the CFD model becomes useful in modelling bioreactors when the main substrates are gas fed and the stoichiometry of the reaction leads to gas phase depletion, especially in the case for highly soluble gas such as carbon dioxide or with low soluble gas such as hydrogen but high gas residence time, i.e., important bubble column height. In the future, a finer mesh and further work on bubble induced turbulence may be required to better describe the bubble plume dynamics and recirculation flow structure. Despite the higher computational cost, detailed CFD simulations provide valuable insights into the fluid flow thanks to a higher spatial resolution that could not be achieved with the 1D model, especially the oscillating behaviour of the bubble plume near the gas sparger. Nevertheless, the 1D model is a useful tool for exploring a wide range of parameters prior to performing computationally constrained CFD calculations. In particular, the 1D model shows that the purity of the methane production is very sensitive to the dispersion of the gas phase in an industrial-scale bubble column. The conceptual design of the bioreactor should then be based on the synergistic application of the high-resolution CFD model and the practical 1D model.

## CRediT authorship contribution statement

**Vincent Ngu:** Writing – original draft, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Jérôme Morchain:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. **Benjamin Casale:** Writing – review & editing, Validation, Software, Investigation. **Arnaud Cockx:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Funding acquisition, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

This work was supported by a government grant managed by the Agence Nationale de la Recherche (ANR) under the France 2030 PEPR BEST program with the reference ANR-22-PEBB-0008. The CFD work was performed using the HPC resources of the CALMIP supercomputing centre under the 2022 grant (P22025).

## Data availability

No data was used for the research described in the article.

## References

Antal, S.P., Lahey, R.T., Flaherty, J.E., 1991. Analysis of phase distribution in fully developed laminar bubbly two-phase flow. Int. J. Multiph. Flow 17 (5), 635-652.

Bauer, M., Eigenberger, G., 1999. A concept for multi-scale modeling of bubble columns and loop reactors. Chem. Eng. Sci. 54, 5109.

Burns, A.D., Frank, T., Hamill, I., Shi, J.-M.M., 2004. The Favre Averaged drag model for turbulent dispersion in Eulerian multi-phase flows. In: Fifth International Conference on Multiphase Flow, ICMF-2004, No. 392, pp. 1–17.

Deckwer, W.-D., Burckhart, R., Zoll, G., 1974. Mixing and mass transfer in tall bubble columns. Chem. Eng. Sci. 29 (11), 2177-2188. https://doi.org/10.1016/0009-2509 (74)80025-4.

Deckwer, W.-D., Adler, I., Zaidi, A., 1978. A comprehensive study on co2-interphase mass transfer in vertical cocurrent and countercurrent gas-liquid flow. Can. J. Chem. Eng. 56, 43-55. https://doi.org/10.1002/cje.5450560107.

Electrochaee - Carbon and Energy storage, 2017.

Heijnen, J.J., Van'T Riet, K., 1984. Mass transfer, mixing and heat transfer phenomena in low viscosity bubble column reactors. Chem. Eng. J. 28, B21-B42. https://doi.org/10.1016/0300-9467(84)85025-X.

Higbie, R., 1935. The rate of absorption of a pure gas into a still liquid during short periods of exposure. Trans AICHE 31, 365-389.

Jensen, M.B., Ottosen, L.D.M., Kofoed, M.V.W., 2021. H₂ gas-liquid mass transfer: a key element in biological Power-to-Gas methanation. Renew. Sustain. Energy Rev. 147, 111209. https://doi.org/10.1016/j.rser.2021.111209.

Lardon, L., Thorberg, D., Krogaard, L., 2018. Biogaz valorization and efficient energy management [WWW Document]. URL http://powerstep.eu/system/files/generated/files/resource/d3-2-technical-and-economic-analysis-of-biological-methanation deliverable.pdf.

Lecker, B., Ill, L., Lemmer, A., Oechsner, H., 2017. Biological hydrogen methanation – a review. Bioresour. Technol. 245, 1220-1228. https://doi.org/10.1016/j.biortech.2017.08.176.

Markthaler, S., Plankenbühler, T., Miederer, J., Kolb, S., Herkendell, K., Karl, J., 2022. Combined two-model CFD simulation of trickle bed reactors with head-sump extension: case study on hydrodynamics and biological methanation. Ind. Eng. Chem. Res. 61, 4134-4152. https://doi.org/10.1021/acs.iecr.1c04262.

Markthaler, S., Plankenbühler, T., Weidlich, T., Neubert, M., Karl, J., 2020. Numerical simulation of trickle bed reactors for biological methanation. Chem. Eng. Sci. 226, 115847. https://doi.org/10.1016/j.ces.2020.115847.

Mast, Y., Ghaderi, A., Takors, R., 2024. Real case study of 600 m³ bubble column fermentations: spatially resolved simulations unveil optimization potentials for l-phenylalanine production with *Escherichia coli*. Biotechnol. Bioeng. https://doi.org/10.1002/bit.28869, 122, 2, (265-286).

Nadal-Rey, G., McClure, D.D., Kavanagh, J.M., Cassells, B., Cornelissen, S., Fletcher, D.F., Gernaey, K.V., 2022. Computational fluid dynamics modelling of hydrodynamics, mixing and oxygen transfer in industrial bioreactors with Newtonian broths. Biochem. Eng. J. 177, 108265. https://doi.org/10.1016/j.bej.2021.108265.

Ngu, V., Morchain, J., Cockx, A., 2022. Spatio-temporal 1D gas-liquid model for biological methanation in lab scale and industrial bubble column. Chem. Eng. Sci. 251, 117478. https://doi.org/10.1016/j.ces.2022.117478.

Ngu, V., Fletcher, D.F., Kavanagh, J.M., Rafrafi, Y., Dumas, C., Morchain, J., Cockx, A., 2023. H2 mass transfer – a key factor for efficient biological methanation: Comparison between pilot-scale experimental data, 1D and CFD models. Chem. Eng. Sci. 268, 118382.

Puiman, L., Abrahamson, B., van der Lans, R.G.J.M., Haringa, C., Noorman, H.J., Picioreanu, C., 2022. Alleviating mass transfer limitations in industrial external-loop syngas-to-ethanol fermentation. Chem. Eng. Sci. 259, 117770. https://doi.org/10.1016/j.ces.2022.117770.

Rafrafi, Y., Laguillaumie, L., Dumas, C., 2020. Biological Methanation of H2 and CO2 with Mixed Cultures: Current Advances, Hurdles and Challenges. Waste Biomass Valorization. Doi: 10.1007/s12649-020-01283-z.

Rusmanis, D., O'Shea, R., Wall, D.M., Murphy, J.D., 2019. Biological hydrogen methanation systems – an overview of design and efficiency. Bioengineered 10, 604-634. https://doi.org/10.1080/21655979.2019.1684607.

Sander, R., 2015. Compilation of Henry's law constants (version 4.0) for water as solvent. Atmos. Chem. Phys. 15, 4399-4981. https://doi.org/10.5194/acp-15-4399-2015.

Sato, Y., Sekoguchi, K., 1975. Liquid velocity distribution in two-phase bubble flow. Int. J. Multiph. Flow 2 (1), 97-195.

Siebler, F., Lapin, A., Hermann, M., Takors, R., 2019. The impact of CO gradients on C. ljungdahlii in a 125m³ bubble column: Mass transfer, circulation time and lifeline analysis. Chem. Eng. Sci. 207, 410-423. https://doi.org/10.1016/j.ces.2019.06.018.

Siebler, F., Lapin, A., Takors, R., 2020. Synergistically applying 1-D modeling and CFD for designing industrial scale bubble column syngas bioreactors. Eng. Life Sci. 20, 239-252. https://doi.org/10.1002/elsc.201900132.

Sveinbjørnsson, D., Münster, E., 2017. WP1: Gas conditioning and grid operation: Upgrading of Biogas to Biomethane with the Addition of Hydrogen from Electrolysis. PlanEnergi, Skørping.

Thema, M., Bauer, F., Sterner, M., 2019. Power-to-gas: electrolysis and methanation status review. Renew. Sustain. Energy Rev. 112, 775-787. https://doi.org/10.1016/j.rser.2019.06.030.

Tomiyama, A., Kataoka, I., Zun, I., Sakaguchi, T., 1998. Drag coefficients of single bubbles under normal and micro gravity conditions. JSME Int J., Ser. B 41, 472-479. https://doi.org/10.1299/jsmeb.41.472.

Tomiyama, A., Tamai, H., Zun, I., Hosokawa, S., 2002. Transverse migration of single bubbles in simple shear flows. Chem. Eng. Sci. 57, 1849-1858. https://doi.org/10.1016/S0009-2509(02)00085-4.

Volger, R., Puiman, L., Haringa, C., 2024. Bubbles and Broth: A review on the impact of broth composition on bubble column bioreactor hydrodynamics. Biochem. Eng. J. 201, 109124. https://doi.org/10.1016/j.bej.2023.109124.

Wachi, S., Nojima, Y., 1990. Gas-phase dispersion in bubble columns. Chem. Eng. Sci. 45 (4), 901-905. https://doi.org/10.1016/0009-2509(90)85012-3.

Zuber, N., Findlay, J.A., 1965. Average volumetric concentration in two-phase flow systems. J. Heat Transf. 87, 453-468. https://doi.org/10.1115/1.3689137.