Hydrodynamics of an Electrochemical Membrane Bioreactor

An electrochemical membrane bioreactor (EMBR) has recently been developed for energy recovery and wastewater treatment. The hydrodynamics of the EMBR would significantly affect the mass transfers and reaction kinetics, exerting a pronounced effect on reactor performance. However, only scarce information is available to date. In this study, the hydrodynamic characteristics of the EMBR were investigated through various approaches. Tracer tests were adopted to generate residence time distribution curves at various hydraulic residence times, and three hydraulic models were developed to simulate the results of tracer studies. In addition, the detailed flow patterns of the EMBR were acquired from a computational fluid dynamics (CFD) simulation. Compared to the tank-in-series and axial dispersion ones, the Martin model could describe hydraulic performance of the EBMR better. CFD simulation results clearly indicated the existence of a preferential or circuitous flow in the EMBR. Moreover, the possible locations of dead zones in the EMBR were visualized through the CFD simulation. Based on these results, the relationship between the reactor performance and the hydrodynamics of EMBR was further elucidated relative to the current generation. The results of this study would benefit the design, operation and optimization of the EMBR for simultaneous energy recovery and wastewater treatment.


Results and Discussion
Hydraulic performance at various HRTs. As shown in Fig. 1, two or three peaks are observed from each RTD curve at four different HRTs; therefore, a short circuiting stream or preferential flow existed in the EMBR because water leaked from different locations of graphite felt in the cathode and produced different channels of water from the inlet to the outlet in the EMBR 19 . Moreover, a long tail appeared in each RTD curve, especially at HRT values of 3.12 and 7.02 h, to imply that stagnant or dead zones were present in the EMBR and the release of the tracer Li + was slow with flow stream in these regions 20 . In the EMBR, most of the anode chamber was taken up by the graphite felt, generating numerous micro Comparison of the hydraulic models. Both TIS and AD models are widely used to simulate the hydrodynamics of the non-ideal flow reactors 21 , while Martin model could give better description on the preferential flow or short circuiting in a reactor 19,22 . Therefore, these three hydraulic models were selected to describe the hydrodynamic properties of the EMBR for comparison in this study.
As shown in Fig. 1, the Martin model generated better simulation results from the experimental data in each RTD curve compared to the others. In other words, neither the TIS nor AD models could characterize the hydrodynamics of the EMBR effectively. Moreover, the simulation reveals that three parallel sub-flows were involved in the Martin model for each HRT. Therefore, the water could flow primarily through the EMBR in three different channels, and each channel could be characterized by the TIS model. As shown in Table 1, the number (N) of CSTR in each strand of the Martin model was not below 3.4 at four HRTs. As discussed above, the degree of back-mixing in the reactor is insignificant if the value of N exceeds 3.01. Therefore, the flow pattern of each strand in the EMBR was similar to a plugged flow.
The dead volume of the EMBR at various HRTs estimated with both the Martin model and an RTD analysis are summarized in Table 1. The values of dead volume are almost identical at each HRT for the two methods, further implying that the Martin model could accurately simulate the hydrodynamic performance of the EMBR. In addition, Table 1 also shows that the dead zone became larger at higher or lower HRT in the EMBR. The liquid had a larger velocity at a shorter HRT; therefore, a vortex flow could form at the corner, possibly generated the higher values for the dead zones in the EMBR. However, a dead corner or stagnant zone would emerge in the reactor if the HRT was too long, possibly generating a larger dead zone in the EMBR.
The porous graphite felt was padded into the anode as the electrode material and dispersed unevenly inside, which resulted in a heterogeneous mixing in the EMBR. Although TIS and AD models are pervasive, however, both of them were developed based on the uniform dispersion degree in a reactor and thus wasn't able to well simulate the heterogeneous mixing in the EMBR. Moreover, both TIS and AD models gave little information about dead volume and short circuiting of the EMBR. On the contrary, Martin model could describe the heterogeneous mixing in the reactors by adopting a number of separate strands, and therefore shown a better simulation performance on the hydrodynamics of the EMBR compared to TIS and AD ones.
Flow pattern visualization through CFD simulation. In order to conduct CFD simulation, the permeability of both graphite felt and non-woven cloth in the EMBR was initially estimated through the porous media model. A porous media model does nothing but add a momentum sink to the governing momentum equations, generating pressure gradients in the porous media 23 . The momentum sink term in the model for simple homogeneous porous media can be described as: where S i is the momentum source, and K and C 2 are the permeability and inertial resistance, respectively. In our study, the inertial resistance of both the graphite felt and non-woven cloth could be treated as zero because water flowed slowly within the entire EMBR. In addition, all of the fibers in the graphite felt and non-woven cloth were treated as though they were randomly located in planes running parallel to gravity. Therefore, the permeability of the graphite felt and non-woven cloth could be determined using Eq. (2) according to a previous report: 24 where ε is the porosity, d f is the fiber diameter, and ε p and α are constants that depend on the arrangement of the fibers. Values of 0.11 and 0.52 were suggested for ε p and α for parallel permeability K P , while 0.11 and 0.785 were used for normal permeability K N
Based on the above estimated values, a transient CFD simulation was performed over 70 s at each HRT, and the results show that the fraction of the liquid volume was always 100% in the entire EMBR. Specifically, no air bubbles could enter into the reactor during the simulation process, and therefore the velocity magnitude contours and vector fields of the mixture only the covered liquid phase in the EMBR region. When using an HRT of 3.12 h as an example, the variations in the liquid contour of velocity magnitude in the EMBR is shown in Fig 2a. The similar contours at 40 and 70 s imply that the fluid flow in the EMBR region had already reached a steady state by 40 s. Additionally, the flow velocity in the graphite felt zone was much lower than that in the non-graphite felt zones because the graphite felt imposed flow resistance on the water stream, causing a momentum loss in the water stream (Fig. 2a). Vesvikar and Al-Dahhan reported that the region where the flow velocity is less than five percent of the maximum among simulation zones could be considered the dead zone 25 . Consequently, the dead zones of the EMBR are located primarily at the bottom and upper outer regions of the reactor, as shown in Fig. 2b. Based on the velocity vector field at 40 s in Fig. 2c, the liquid flowed regularly with less back mixing or mixing with others in the anode without graphite felt, indicating that the flow pattern was closed to generate a plug flow in this area. However, the main flow of the stream would disperse to numerous branches after it entered the graphite felt zone of the anode and then discharge from different locations in the graphite felt in the cathode. This behavior would allow the generation of various channels of liquid from inlet to outlet in the EMBR, indicating the possibility of existing preferential or circuitous flow in the reactor.

Relationship between current generation and hydrodynamics. The simulated DO distributions
in the cathode of the EMBR at various HRTs are shown in Fig. S1. The DO concentration should have decreased when increasing the distance from the outside surface of the cathode at each HRT. However, a high flow velocity of liquid at a short HRT produced a low DO concentration in the cathode because the DO diffused in a direction opposing the flow of water in the EMBR.
The existence of dead zones will decrease the effective volume in bioreactors, affecting the extent of the biochemical reactions. However, the oxygen transfer process is critical for air-cathode MFCs because it can determine the rate of the oxygen reduction reaction. Fig. 3 shows the variations in the current density in the EMBR with the effective volume of the reactor and average DO concentration in the cathode and multiplier. A positive relationship was found only between the current density and multiplier of two factors rather than a single one in the EMBR. Therefore, the power production in the EMBR was simultaneously determined through the effective volume of the reactor and the DO concentration in the cathode, which were both strongly affected by the hydrodynamics of the reactor, as discussed above. The DO had a lowest concentration because it had the highest flow velocity at HRT = 3.12 h (Fig. S1), therefore, the electricity generation in the EMBR was seriously limited, resulting in a minimum current density of 4.7 A/m 3 . Longer HRTs would produce higher DO concentrations in the cathode and may be beneficial for electricity production, but the dead volume of the reactor would increase and thus reduce the capability for electricity production in the EMBR (Table1). At HRT = 17.73 h, the ratio of the dead zones (29.3%) was significantly higher compared to the other HRTs, generating a lower current density of 9.2 A/m 3 in the EMBR. These results demonstrate that the hydrodynamic characteristics of the reactor are critical for the performance of the EMBR. An optimized HRT could reduce the dead volume of the reactor while also providing an appropriate DO concentration in the cathode, thereby promoting electricity production in the EMBR.
In conclusion, the revealed hydrodynamic characteristics of the EMBR would help us know the degree of mixing as well as the value and location of dead zones, which will be conducive to improvement of the EMBR structure. Moreover, after an analysis on the relationship between the revealed hydrodynamic properties and the performance of the EMBR, it has been confirmed that HRT had a significant impact on the current generation, and thus it is imperative to choose an optimized HRT for the operation of the EMBR.

Methods
EMBR structure and operation. The schematic diagram of EMBR is shown in Fig. S2a. Non-woven cloth with 75% porosity and 50-mm pores supported by a perforated polyvinyl chloride (PVC) tube was placed between cathode and tubular anode chambers. The non-woven cloth could separate anodic and cathodic electrodes thus avoiding short circuit. Graphite felt with 93.5% porosity and 150-μ m pores (Sanye Carbon Co., China) was served as the electrode both in anode and cathode. In the anode chamber (height 20 cm, diameter 4.5 cm), the graphite felt with a total volume area of 197.6 mL was rolled into a swirl so that the water from the top of the EMBR could flow more uniform, and the working volume of the anode chamber was 247.1 mL. The cathode of the EMBR was composed of a graphite felt with a projected surface area of 294 cm 2 that wrapped around the non-woven cloth. Besides, electrodes were connected to the circuit through titanium wires across an external resistor (100 Ω), the reactor was operated in a continuous-flow mode at around 25 o C 9,26 . The electrons and protons were generated when consuming substrate by the electrochemically active bacteria at the anode, while at the cathode electrons from the anode was reacted with oxygen in the air and protons diffused from anode.
The electrochemical reactor usually consists of anode and cathode as well as the separator membrane. In the EMBR, the membrane was replaced by the non-woven cloth, which could be used as not only a separator but also a filter for treated water from anode. Besides, the bacteria was inoculated both in the anode and cathode of the EMBR to consume organic and inorganic matters. Additionally, both anodic and cathodic pHs of the EMBR were kept at a natural value because of bacteria inside, and moreover the operational HRT would not be too short due to the biofilm formation in both compartments.
Tracer test. Tracer test is often used to produce RTD curves in different systems. Lithium-ion has been widely used as the tracer not only due to its simple measurement procedure, low detection limit and low background, but also because it does not volatilize, precipitate or absorb in the reactor. Based on lithium-ion tracer experiments, Capela et al. obtained information about mixing degree and non-effective volume in a full-scale anaerobic contact reactor 27 , while Wang et al. found that dosing of ferrous salts at two different points would bring two distinct effects on the reactor in a pilot scale membrane reactor 28 . Besides, lithium tracer has also been applied to provide RTD curves in an electrochemical system 29 .
In this study, a 2 g/mL aqueous Li 2 SO 4 solution was used as the tracer in the experiments, and the pulse input method was adopted through the use of deionized water as the carrying fluid. One milliliter of the Li 2 SO 4 solution was injected into the flowing water of the EMBR over 0.5 s, producing an average Li + concentration of 10 mg/L in the entire reactor. As summarized in Table S1, four different hydraulic retention times (HRTs) were investigated by controlling the flow rate of the EMBR with a peristaltic pump (Longer Co., China). The samples were collected from the EMBR effluent during injection before stopping after fourfold HRTs 30 . The Li + concentration in the sample was measured through flame emission spectroscopy (Vario 6, Analytik Jena AG, Germany) above 670.8 nm according to the Standard Methods 31 .

RTD analysis. The exit age-distribution function E(t) can be expressed as follows:
where c(t) is the tracer concentration at the outlet of the reactor 32 . The mean residence time t m is calculated as follows: when comparing the mixing performance of the differently sized flow systems, all of the parameters should be normalized versus normalized time θ: where τ is the hydraulic residence time, V R is reactor volume, and Q is volume flow rate. Therefore, the normalized RTD function E(θ) can be presented as follows: where M is the total amount of tracer. The dimensionless variance 2 σ θ , which characterizes the degree of back-mixing degree in the reactors, is defined as follows: N is a criterion that predicts the flow pattern in a reactor 32 . The flow pattern of the reactor is closer to a complete mixing flow if N remains below 3.01; in all other cases, the pattern is more similar to a plug flow 21 . The exit age-distribution function E(θ) of the TIS model is given as follows: Where D is the axial dispersion coefficient, and D/uL is defined as the Peclet number (Pe), which represents the ratios of mass transport attributed to advection and dispersion, respectively 32 . The Peclet number is also a criterion that predicts the degree of back-mixing in a reactor, and the relationship between 2 σ θ and Pe can be described as follows: Pe Pe e 2 2 1 1 12 Martin model. This model can evaluate the short circuiting stream or preferential flow in reactors while determining the dead volume accurately 19 . A schematic diagram of this model for the EMBR is shown in Fig. 4c, consisting of several parallel strands; each strand is considered as the number of CSTRs in series. Therefore, the following equation was obtained based on the mass balance:   where t peak and θ peak are the time and dimensionless time of the maximal tracer concentration at each strand outlet, respectively. In this study, the Microsoft EXCEL built-in solver ® was used to determine the values of various parameters in each model to reach the best correlation between the experimental data and the simulation results.
CFD simulation. A 2-D axisymmetric unsteady CFD simulation was used to investigate the hydrodynamics of the EMBR through commercially available CFD software called FLUENT. The Gambit pretreatment software was chosen to divide grids in the simulated region, which was treated as a 2-D plane because an axisymmetric space was adopted for the simulation in this study, as shown in Fig. S2b. The detailed boundary conditions of this simulation are provided in the Supporting Information (SI). The Volume of Fluid (VOF) model, which is an Euler-Euler approach, was applied to address the multiple phases: gas and water 33 . In the VOF model, two or more phases share a single set of conservation equations, and the volume fraction of each phase is tracked in each computational cell throughout the domain. The conservation equations for mass and momentum in the VOF formulation are defined as follows, respectively: The density (ρ) and viscosity (μ) in a two-phase system can be estimated using Eq. (21) and Eq. (22):