Influence of Electrodes Configuration on Hydraulic Characteristics of Constructed Wetland–Microbial Fuel Cell Systems Using Graphite Rods and Plates as Electrodes

<jats:p>Constructed wetland–microbial fuel cell coupling systems (CW–MFCs) have received significant academic interest in the last decade mainly due to the promotion of MFCs in relation to pollutants’ degradation in CWs. Firstly, we investigated the effect of hydraulic retention time (HRT) and electrode configuration on the flow field characteristics of CW–MFCs using graphite rods and plates as electrodes, as well as the optimization of electrode configuration using computational fluid dynamics (CFD) numerical simulation. The results showed that: (1) the apparent HRT was the most influential and decisive factor, with a contribution of over 90% for the average HRT of CW–MFCs; (2) anode spacing was the most influential factor for the hydraulic performance of CW–MFCs, with contributions of over 50% for water flow divergence and hydraulic efficiency (λ) and over 45% for effective volume ratio (e); (3) anode size was significant for e and λ, with a contribution of over 20%; (4) cathode position and cathode size had no statistically significant effect on the hydraulic performance of CW–MFCs. It was mainly through the blocking of water flows, flows around, compressing water flow channels and boundary layer separation that the MFC electrodes influenced the hydraulic characteristics of the flow field in CW–MFCs. Optimizing the flow field by optimizing the electrode configuration helped to facilitate electricity generation and pollutants’ removal in CW–MFCs. This study offers a scientific reference for improving the hydraulic performance of CW–MFCs, and it also provides a new research perspective for improving the wastewater treatment and electricity production performance of CW–MFCs.</jats:p>

section_content: Constructed wetlands (CWs) have found applications worldwide, owing to their characteristics of convenient operation, easy maintenance, significant ecological and environmental benefits, wastewater resource utilization, and low costs [1, 2]. Microbial fuel cells (MFCs) use microorganisms as biocatalysts to oxidize organic substrates and generate electricity [3, 4]. Compared with traditional chemical fuel cells, MFCs possess the advantages of high efficiency, zero pollution, wide fuel sources, mild operating conditions, and high biocompatibility [4]. Since wastewaters, especially domestic wastewater, contain a large amount of organic components and nutrients, many studies have combined MFCs with other wastewater treatment technologies to accelerate pollutants' degradation and realize resource recovery [5, 6]. The anode and cathode of a typical MFC are required to remain anaerobic and aerobic, respectively [4]. Since the redox gradient naturally formed in CWs could satisfy the requirement for MFC operations, many studies conducted in the last decade have incorporated MFCs into CWs to achieve simultaneous wastewater treatment improvement and bioelectricity production [7, 8]. 

The requirements of MFCs for electrode materials include high conductivity, noncorrosion, non-toxicity, good chemical stability, a high specific surface area, a low cost, and availability. Based on the above requirements and combined with the characteristics of CWs, the most commonly used electrodes in CW-MFC systems are carbon-based materials, such as carbon cloth, carbon felt, carbon fibers, carbon brushes, activated carbon, graphite rods, graphite plates, and graphite particles [8]. MFC electrodes are bound to influence the water flow's path and the hydraulic retention time (HRT) distribution of wastewater in CWs. Considering that the water flow inside CWs is an important factor influencing the wastewater treatment performance, it is necessary to reveal the influence mechanism of MFC electrode configurations on the hydraulic characteristics of CW-MFCs. 

Tracer experiment and numerical simulation are the two main methods employed to analyze the hydraulic characteristics of CWs. The tracer experiment possesses higher reliability, but requires the construction or modification of the CW system, which is more time-consuming and requires a higher workload [9]. The numerical simulation can obtain results quickly and economically on a computer by simulating the operation of devices under various conditions, especially complex conditions [10]. As long as the physical object is reasonably simplified and the error is controlled within a certain range during the simulation process, the numerical simulation can play an important role in the development, optimization, and performance prediction of devices. 

Computational fluid dynamics (CFD) simulation software presented by FLUENT is most often used to simulate the flow field characteristics of CWs. Fan et al. [10] employed CFD simulation to investigate the effect of wetland configuration (inlet location, constructed media, and protection layer) and operating conditions (inlet velocity and outlet pressure) on the hydraulic performance of a subsurface flow wetland (SSFW) and found that the hydraulic efficiency of SSFW was predominantly affected by the wetland configuration; in particular, (1) the SSFW with an inlet centrally located at the edge of the upper media had better hydraulic efficiency than the wetland did with an inlet at the top and bottom edge of the upper media, (2) the effect of the constructed media was complicated and must be carefully adjusted to increase the hydraulic efficiency, and (3) the higher resistance in the protection layer benefited the hydraulic efficiency of SSFW. Han et al. [11] employed CFD simulation to predict the flow characteristics and distribution of suspended solids (SS) in surface flow constructed wetlands (SFCWs), and thus determined the effectiveness of SFCWs in removing SS, and the results showed that the 3D CFD model performed reasonably well at predicting complex flow fields associated with complex wetland geometry. Rajabzadeh et al. [12] developed a robust CFD model, accounting for both spatial and temporal dynamics of a subsurface vertical flow wetland by combining fluid transport, solute transport, biokinetics, biofilm development, and biofilm detachment/sloughing using COMSOL Multiphysics TM, and this model can be applied to predict bio-clogging processes in a spatial manner similar to what would be realistically expected. Rengers et al. [13] employed CFD simulation to determine the empirical effects of geometric and flow parameters on the hydraulic performance of a horizontal SSFW and found that the length, the baffles, and the interaction between the length and baffles had significant statistical influence on the hydraulic efficiency. Noh et al. [14] designed the cross-section of CWs using CFD, and the simulation results showed that installing baffles improved the precipitation efficiency of particulates by decreasing the water flow velocity; in addition, the vertical baffle functioned better than the angled baffle did. Yang et al. [15] analyzed the effect of clogging on hydraulic behavior in a vertical-flow CW using CFD simulation combined with conservative tracer tests and found that the reduction of the pore volume caused by clogging was not the only reason for the decrease in the actual HRT. Maurer et al. [16] employed CFD simulation to distinguish velocity areas in CWs and investigate the influence of water flow on micropollutants' degradation and found that the number and concentration of micropollutants were independent of the water flow, but micropollutant metabolites were detected in higher proportions (both number and concentration) in lower flow rate areas. Wang et al. [17] employed CFD simulation to systematically evaluate the hydraulic performance of an up and down baffled SSFW with different baffle settings and substrate laying methods and found that (1) it was the position of baffles, not the number of baffles, that had significant influence on the hydraulic efficiency of the SSFW, and (2) by increasing the thickness of the middle substrate, the low flow rate phenomenon of the upper substrate and rapid outflow phenomenon of the lower substrate can be improved to a certain extent, thereby improving the hydraulic performance. 

However, no information is available yet regarding the use of a numerical simulation method to assess the hydraulic performance and electrode configuration optimization of CW-MFC systems. Due to such a lack of knowledge, the present study aimed to develop a method to simulate the flow field characteristics of CW-MFCs using CFD, and the specific objectives were to: (1) investigate the influence mechanism of MFC electrode configurations and HRT on the hydraulic characteristics of CW-MFCs by designing a multifactor orthogonal experiment; (2) select and explain the electrode configuration parameters mostly influencing the hydraulic performance of CW-MFCs; (3) identify the optimal electrode configuration to optimize the construction of CW-MFCs. It is expected that the results of this study will reveal the influence mechanism of MFC electrode configurations on the wastewater treatment performance of CW-MFCs from a hydrodynamics perspective.
section_name: Introduction

section_name: Materials and Methods

section_content: The numerical simulation and analysis of the porous medium flow field inside the CW-MFC system was carried out using FLUENT, and the porous medium model and particle trajectory model were adopted to investigate the HRT distribution of water, especially focusing on flow field distribution features near the MFC electrodes.
section_name: Numerical Simulation Method

section_content: The flow of water through the substrates of CWs can be regarded as the movement of water in pores of the porous medium, which can therefore be simulated using the porous media model in FLUENT. In essence, the porous media model adds a source term (S i ) that represents momentum consumption to the momentum equation, so as to simulate the flow resistance generated by the porous media. The formula is as follows:

where S i is the source term in the i-th momentum equation; D and C are given matrices;

2 ρv mag v j is the inertial loss item; µ is the dynamic viscosity of the fluid (Pa•s); v j is the velocity of a point in the flow field along a certain direction (m/s); ρ is the density of the fluid (kg/m 3 ). The negative source term is also called the "sink". The contribution of momentum sink to the momentum gradient of the porous media unit is to produce a pressure drop on the unit that is proportional to the fluid velocity (or velocity squared). 

For simple and uniform porous media, D and C can be defined as diagonal matrices, with 1/α and C 2 as diagonal units, respectively. The mathematical model is as follows:

where α represents the permeability of the porous media (m 2 ), and 1/α is the viscous drag coefficient (m -2 ); C 2 is inertial drag coefficient (m -1 ); d p is the average diameter of substrate particles (m); ε is the porosity of substrates.
section_name: Mathematical Model Porous Media Model

section_content: The particle trajectory model, also called the Euler-Lagrange model, treats the particles as a discrete phase and the fluid as a continuous phase. The motion of discrete particles and the continuous fluid is calculated in the Lagrangian coordinate system and Eulerian coordinate system, respectively. The macroscopic trajectory of particle swarms is obtained by counting the instantaneous position of a large number of particles, and the residence time distribution (RTD) of particle swarms in the particle trajectory model is approximated as the RTD of water flow inside the wetland. 

When a real particle moves in a liquid, the main forces acting on the particle include gravity, resistance, virtual mass force, pressure gradient force, Basset force, Saffman force, Magnus force, etc. As for a virtual particle in a liquid, only gravity, resistance, and virtual mass force are taken into consideration. According to the force analysis, the equation of particle motion is as follows:

where m p is the mass of particle (kg); v p is the velocity of particle (m/s); t is the time (s); F G is the gravity of particle (N); F R is the resistance (N); F V is the virtual mass force (N); F P is the pressure gradient force (N); F B is the Basset force (N); F X is the sum of other external forces (N). In this paper, the particles on the orbit were set to be shot randomly into the wetland model from the inlet surface. The number of particles was determined according to the selected particle concentration and area of the model inlet surface, which was a medium value of 1500 particles/m 2 and 0. 2 m 2, respectively, in this study, and thus the number of particles was a total of 300.
section_name: Particle Trajectory Model

section_content: The process of the numerical simulation method for the flow field of CW-MFCs is shown in Figure 1, and the main steps were carried out as follows:

Step one: Establish the geometric model of the CW-MFC system using ANSYS modelling software ICEM. Specifically, determine the geometric size of the simulated CW, the height and particle size distribution of the CW substrates, and the size and spacing of electrodes, and then input the abovementioned parameters in the software to build the geometric model of the CW-MFC system. The geometric model of the CW-MFC system in this study is illustrated in Figure 2A, and its main structural parameters are listed in Table 1. 

Step two: Analyze the geometric model of the CW-MFC system, identify the fluid part (porous media basin), the solid part (electrodes), and the boundary part (inlet cross-section, outlet cross-section, internal cross-section, and wall surface), and establish a meshed model of CW-MFC. Specifically, first define the incorporated MFC electrodes as solid; secondly, define the inlet cross-section, outlet cross-section, internal cross-section, and wall surface of the model and confirm the different flow area inside the CW substrates; finally, perform automatic volume recognition for the area through which the fluid flows by using the volume setting function of software ICEM and mesh each layer of fluid using an unstructured meshing method; then, output the mesh file. The 3D mesh model of the CW-MFC system in this study is illustrated in Figure 2B. The hexahedral mesh was adopted; the maximum mesh size was 0. 004, and the total number of calculated mesh units was 1,428,680. It should be emphasized that proper mesh encryption was required near the electrodes. 

Step three: Simulate the meshed model of the CW-MFC system using FLUENT and obtain the flow velocity field characteristics of CW-MFC. Specifically, first import the mesh file obtained in step two into FLUENT, replace the CW substrates with the porous media model, and set the corresponding porous media parameters in different flow areas according to the substrate distribution; secondly, set the porous media flow area as the laminar flow model, and set the outlet area as the turbulent flow model; thirdly, set the boundary conditions, including the wastewater inlet as the velocity inlet boundary condition, the outlet as the pressure outlet boundary condition, and the wall as an adiabatic and non-slip wall surface; finally, obtain data of velocity distribution by simulating changes in the flow field in the geometric model of CW-MFC under a certain inlet velocity. 

Step four: Obtain hydraulic characteristic parameters of the flow field of CW-MFC. Specifically, first, according to the velocity distribution data obtained in step three, inject orbital particles from the inlet cross-section on the basis of a stable flow field, and set their physical parameters to obtain the orbital particle movement data of the CW-MFC model under a certain inlet velocity; secondly, approximate the HRT distribution of particle swarms in the model as the HRT distribution of pollutant inside CW-MFC, and calculate the actual HRT inside CW-MFC; finally, compare the actual HRT with the set HRT to obtain the hydraulic efficiency value of CW-MFC. 

Step five: Adjust the size and spacing of electrodes and influent HRT, repeat steps one to four to obtain the hydraulic efficiency values of CW-MFC under different electrode configurations and inlet velocities, analyze the influence of electrode configurations and influent HRT on the hydraulic efficiency of CW-MFC, and improve the hydraulic efficiency of CW-MFC by continuously optimizing the electrode configuration. The following assumptions were made while modeling: (1) the porous media were considered to have a uniform dispersion structure and be isotropic; (2) the inlet velocity was evenly distributed, and the water flow was steady; (3) the influence of reaction exotherm and solid substrates on the temperature was not considered, that is, the entire flow was an isothermal process; (4) the adiabatic wall was adopted, and the influence of reaction and heat transfer was not considered. The upwind style with second-order precision was chosen for discrete interpolation using the pressure-based three-dimensional steady solver (FLUENT-2021). When the residuals of the variables in the continuity equation and momentum equation were both below 10 -5 and did not change with the calculation, the result was considered to be convergent. During the calculation, the compression factor needed to be continuously adjusted to make the results converge.
section_name: Process of Numerical Simulation Method

section_content: The hydraulic evaluation parameters calculated from the RTD included the apparent HRT (T n ), average HRT (T m ), variance (σ 2 ), water flow divergence (σ 2 0 ), effective volume ratio (e) and hydraulic efficiency (λ). These parameters are defined by the following equations [17] [18] [19]. 

where T n is the apparent hydraulic residence time (h), also known as the theoretical hydraulic residence time; V is the volume of the wetland (m 3 ); ε is the porosity of the substrate; Q is the volume flow rate (m 3 /h); T m is the average hydraulic residence time (h) and is obtained by calculating the average of the hydraulic retention times of all particles in the particle trajectory model; σ is the standard deviation from the average residence time, and σ 2 is the variance; σ 2 0 is the water flow divergence and is used to judge the water flow pattern in the wetland (σ 2 0 = 1 indicates fully mixed flow, and σ 2 0 = 0 indicates plug flow); e is the effective volume ratio, and when e < 1, it indicates that short-cut flow exists in the wetland, but when e > 1, it indicates that dead zones or stagnant water zones exist in the wetland; λ is the hydraulic efficiency, and the larger the value of λ is, the higher was the internal hydraulic efficiency of wetlands is. 

Additionally, the inlet flow velocity was calculated as follows,

where v is the inlet flow velocity (m/s); V is the wetland volume (m 3 ); T n is the apparent hydraulic residence time (h); A is the area of the inlet section (m 2 ); ε is the porosity of substrates.
section_name: Evaluation of Hydraulic Characteristics

section_content: The parameters of electrode configuration in CW-MFCs include electrode size, electrode spacing, and electrode position. Based on the requirement of MFCs for the redox gradient, in a CW-MFC system, the MFC anode is normally installed in the interior of CW substrates, and the MFC cathode is normally installed on the CW substrate surface or near the plant roots located in the upper layer of substrates. Hence, the anode size, anode spacing, cathode position, and cathode size were selected as the factors affecting the flow field of CW-MFCs. As a common and important factor, HRT was also taken into consideration. 

In order to investigate the influence of the above factors on the hydraulic characteristics of CW-MFCs more effectively with fewer experiments, An L 25 (5 5 ) orthogonal experiment was designed in this study assuming that no strong interaction effects existed among these five factors. Details of the factors and levels are provided in Table 2. Analysis of variance (ANOVA) was performed to illustrate the statistical significance of the impact of investigated factors on the hydraulic performance of CW-MFCs. All of the statistical analyses were conducted using SPSS and Origin software. iment was designed in this study assuming that no strong interaction effects existed among these five factors. Details of the factors and levels are provided in Table 2. Analysis of variance (ANOVA) was performed to illustrate the statistical significance of the impact of investigated factors on the hydraulic performance of CW-MFCs. All of the statistical analyses were conducted using SPSS and Origin software.
section_name: Design of Orthogonal Experiment

section_content: According to the numerical simulation results, five groups of representative electrode configuration combinations (A1B3C3D1E3, A2B3C4D2E1, A3B2C4D1E3, A4B3C1D2E2, and A5B3C2D1E5) were selected for tracer experiments to verify the accuracy of the numerical simulation. Lab-scale physical CW-MFC reactors were built according to the dimensions and substrate layers of the geometric CW-MFC model (Figure 2A ). NaCl was chosen as the tracer. At the beginning of the tracer experiment, NaCl solution was instantly added to the influent of the CW-MFC systems one time, so that the electrical conductivity (EC) rose to more than 10 times larger than the background value. Sampling was performed every 20 min at the top outlet to monitor the changes in conductivity until the EC returned to the background value. When the conductivity fluctuated violently, the sampling interval can be shortened appropriately. Each experiment was repeated three times, and the result with the highest tracer recovery rate was selected for the calculation of hydraulic efficiency. 

According to the fluid reactor theory, the NaCl concentration obtained by converting the effluent EC value was taken as the density distribution of residence time [20]. Therefore, the effluent EC value was standardized according to Equation ( 6 ) to obtain the RTD density (N(t)), and then we plotted the time t as the abscissa to obtain the RTD [21, 22].
section_name: Verification of the Accuracy of Numerical Simulation

section_content: According to the numerical simulation results, five groups of representative electrode configuration combinations (A1B3C3D1E3, A2B3C4D2E1, A3B2C4D1E3, A4B3C1D2E2, and A5B3C2D1E5) were selected for tracer experiments to verify the accuracy of the numerical simulation. Lab-scale physical CW-MFC reactors were built according to the dimensions and substrate layers of the geometric CW-MFC model (Figure 2A ). NaCl was chosen as the tracer. At the beginning of the tracer experiment, NaCl solution was instantly added to the influent of the CW-MFC systems one time, so that the electrical conductivity (EC) rose to more than 10 times larger than the background value. Sampling was performed every 20 min at the top outlet to monitor the changes in conductivity until the EC returned to the background value. When the conductivity fluctuated violently, the sampling interval can be shortened appropriately. Each experiment was repeated three times, and the result with the highest tracer recovery rate was selected for the calculation of hydraulic efficiency. 

According to the fluid reactor theory, the NaCl concentration obtained by converting the effluent EC value was taken as the density distribution of residence time [20]. Therefore, the effluent EC value was standardized according to Equation ( 6 ) to obtain the RTD density (N(t)), and then we plotted the time t as the abscissa to obtain the RTD [21, 22]. 

where t is the time after the tracer experiment started (h); N(t) is the standardized residence time distribution density (h -1 ); E(t) is the effluent electrical conductivity at time t (mS/m); E w is the influent background conductivity (mS/m); M NaCl is the molar mass of NaCl (g/mol; M NaCl = 58. 44 g/mol); λ Na + and λ Cl -are the molar conductivities of Na + and Cl -, respectively (S•m 2 /mol; λ Na + = 5. 01 × 10 -3 S•m 2 /mol and λ Cl -= 7. 63 × 10 -3 S•m 2 /mol); Q is the volume flow rate (m 3 /h); m is the total mass of the added tracer (g). 

The average HRT ™ and variance (σ 2 ) obtained from the RTD are calculated as follows,

Sustainability 2023, 15, 6397 9 of 21

The other four hydraulic evaluation parameters (T n, σ 2 0, e, and λ) are calculated using Equations ( 1 )-( 4 ). 

The tracer recovery rate (F(t)) is calculated as follows,
section_name: Verification of the Accuracy of Numerical Simulation

section_content: Pollutants' removal and electricity generation were investigated to evaluate the performance of CW-MFCs with optimized electrode configurations. Three lab-scale physical CW-MFC reactors were built according to the dimensions and substrate layers of the geometric CW-MFC model (Table 1 ). Canna indica was selected as the wetland plant. The anode (graphite rods) and cathode (graphite plates) were connected across a variable external resistor (0-9999. 9 Ω) with titanium wires, and the external resistors were set to 1000 Ω. The synthetic domestic wastewater was continuously dosed into the reactors from the bottom using a peristaltic pump. 

The influents were prepared using peptone, glucose, CH Sampling began when the output voltages of CW-MFCs stabilized. Water samples were collected every three days from the effluent and were analyzed immediately for chemical oxygen demand (COD), ammonia (NH 3 -N), total nitrogen (TN), and total phosphorus (TP) using a multi-parameter water quality tester (MI-200H, Tianjin ZKCP Technology Co., Ltd., Tianjin, China). The tests were repeated twice. A total of 15 water samples were taken for each CW-MFC reactor. The pollutant removal percent was calculated as follows:

where R is the pollutant removal percent (%); C i and C e are the mean influent and effluent concentrations (mg/L), respectively. The parameters of electricity generation performance include output voltage, current and power density. The voltage was collected using a multi-channel data logger (Model CT-4008-5v10mA-164, Shenzhen Neware Electronics Co., Ltd., Shenzhen, China) and recorded by a computer at intervals of 30 min. The current was calculated by Ohm's law. The power density was calculated as follows:

where P d is the power density (mW/m 3 ), P is the power (mW), V is the volume of the wetland (m 3 ), U is the output voltage (V), and R ex is the external resistance (Ω).
section_name: Evaluation of Pollutants' Removal and Electricity Generation

section_name: Results and Discussion

section_content: Table 4 lists the design matrix of the L 25 (5 5 ) orthogonal array and numerical simulation data of the hydraulic characteristics. Figure 3 illustrates the hydraulic characteristics of CW-MFCs obtained from the tracer experiments, so as to verify the accuracy of the numer-ical simulation. The effective volume ratio (e) and hydraulic efficiency (λ) of CWs reported by the authors of previous studies were 0. 49-0. 87 and 0. 301-0. 775, respectively, which fluctuated violently under different wetland configurations and operating conditions [10, 17]. The values of e and λ obtained for the CW-MFC systems in this study were 0. 469-0. 763 and 0. 350-0. 737, respectively, which were basically within the scope of previous studies. Using the data in Table 4 and Figure 3, we calculated that under the same HRT and electrodes configuration, the average HRT (T m ), effective volume ratio (e), and hydraulic efficiency (λ) of the CW-MFC obtained from the numerical simulation were lower than those obtained from the tracer experiment by 4. 2-10. 5%, 4. 2-10. 5%, and -0. 89-7. 6%, respectively, and the corresponding averages were 7. 7%, 7. 7%, and -2. 4%, respectively. Moreover, the values of water flow divergence (σ 2 0 ) obtained from both the numerical simulation (0. 006-0. 434) and tracer experiment (0. 007-0. 291) were relatively low, which indicated basically the same water flow pattern in the CW-MFCs. All the values of e obtained from both the numerical simulation (0. 469-0. 763) and tracer experiment (0. 625-0. 806) were less than 1. This indicated that a short-cut flow existed in the CW-MFC system, and the water flowed out of the system in a short time through fast channels, which suggest that the theoretical HRT could not be reached. Overall, the errors of the hydraulic characteristic parameters obtained by the numerical simulation were within the acceptable range. ANOVA was performed to evaluate the significance and contribution of the investigated factors on the hydraulic performance of CW-MFCs, and the results are listed in Table 5. The apparent HRT (A) had an extremely significant and decisive influence on the average HRT, with a PC value of 92. 0%, which was totally in accordance with expectations. The anode size (B) was extremely significant for the hydraulic efficiency (PC = 24. 44%) and significant for the effective volume ratio (PC = 22. 32%). Anode spacing (C) was found to be extremely significant for the hydraulic efficiency (PC = 60. 0%) and significant for the water flow divergence (PC = 52. 58%) and effective volume ratio (PC = 45. 54%). Additionally, it can be seen from Table 5 that no significant interaction effects existed among the investigated five factors. A, apparent hydraulic residence time; B, anode size; C, anode spacing; D, cathode position; E, cathode size; T m, average hydraulic residence time; σ 2, variance in hydraulic residence time distribution; σ 2 0, water flow divergence; e, effective volume ratio; λ, hydraulic efficiency; SS, sum of squares; DOF, total degree of freedom; MS, mean square; PC, percent contribution, indicating the contribution of the factor to the system performance (PC = SS/total SS); Sig., significance; *, significant influence (0. 01 < p < 0. 05); **, extremely significant influence (p ≤ 0. 01). 

Overall, according to the variance analysis, the sequences and degrees of the influence of the tested factors were A** > C ≈ B > D > E for the average HRT (T m ), A ≈ C > E > B > D for the variance in HRT distribution (σ 2 ), C * > A > E > B > D for the water flow divergence (σ 2 0 ), C* > B* > D > E > A for the effective volume ratio (e), and C** > B** > E > D > A for the hydraulic efficiency (λ) (* denotes significant influence (0. 01 < p < 0. 05) and ** denotes extremely significant influence (p ≤ 0. 01)). It can be concluded that (1) the HRT and electrode configurations barely influenced HRT distribution divergence in CW-MFCs, (2) the size and spacing of anodes were significant and important for the hydraulic performance of CW-MFCs, and (3) the cathode position and size had no statistically significant effect on the hydraulic performance of CW-MFCs.
section_name: Variance Analysis of the Orthogonal Experiment Results

section_content: The streamline diagram of the flow field inside the CW-MFCs obtained from the numerical simulation is provided in However, the flow fields of all the CW-MFCs shared a common characteristic: the flow velocity inside the CW-MFCs kept increasing along the water flow path from the bottom to higher up and reached the maximum near the top of the graphite rod anodes, then slowed down and gradually tended to be steady over the top of the anodes. When the water flowed upward, the water-carrying section area decreased due to the blocking of the electrodes, the water body tended to shrink, and the pressure decreased, and the flow velocity thus increased sharply. After the water flow passed through the electrodes, the water-carrying section area suddenly expanded, resulting in the boundary layer separation of the water flow above the electrodes. After separation, the shear layer spread rapidly and continuously exchanged momentum with the surrounding water body through convection and diffusion; therefore, a smooth water flow and uniform flow field distribution were achieved somewhere above the electrodes. Additionally, the substrate distribution was the secondary reason for the continuous increase in the flow velocity during the upward flow of water, since the particle size and porosity of the three layers of substrates decreased from the bottom to the top; thus, the water-carrying section area also decreased layer by layer. 

Tendency analysis was carried out to determine the optimal level of each factor, and thus the optimal operating conditions for the hydraulic performance of CW-MFCs, and the tendency chart is illustrated in Figure 5. 4 for each impact factor at level i (i = 1, 2, 3, 4, and 5) (A, apparent hydraulic residence time; B, anode size; C, anode spacing; D, cathode position; E, cathode size).
section_name: Influencing Mechanism Analysis and Electrode Configuration Optimization

section_content: As shown in Figure 5, the average HRT (Tm) grew almost linearly with the increase in apparent HRT (Tn), which visually demonstrates the decisive influence of the apparent HRT on the average HRT. The variance in HRT distribution (σ 2 ) showed a general upward trend as the apparent HRT increased; however, the water flow divergence (σ 2 0 ), effective 4 for each impact factor at level i (i = 1, 2, 3, 4, and 5) (A, apparent hydraulic residence time; B, anode size; C, anode spacing; D, cathode position; E, cathode size).
section_name: Apparent Hydraulic Residence Time

section_content: As shown in Figure 5, the average HRT (T m ) grew almost linearly with the increase in apparent HRT (T n ), which visually demonstrates the decisive influence of the apparent HRT on the average HRT. The variance in HRT distribution (σ 2 ) showed a general upward trend as the apparent HRT increased; however, the water flow divergence (σ 2 0 ), effective volume ratio (e) and hydraulic efficiency (λ) did not show clear changes with the apparent HRT. Moreover, lower values of σ 2 0 (0. 071-0. 196) indicated that the internal water flow in CW-MFCs was closer to push flow under the investigated levels of the apparent HRT (0. 25-1. 5 d), and the values of e (0. 612-0. 658) and λ (0. 535-0. 575) varied a little with the apparent HRT. 

The apparent HRT determined the influent flow velocity. The longer the apparent HRT is, the lower the influent flow velocity is, and the greater the disturbance and blocking of the electrodes on the water is. Hence, as the apparent HRT increased, the uniformity of the flow field distribution reduced to some extent. Despite the change in the hydraulic characteristic parameters the apparent HRT, variance analysis showed that the apparent HRT in the range from 0. 25 d to 1. 5 d has no statistically significant effect on the HRT distribution divergence, water flow pattern, or hydraulic efficiency of CW-MFCs. However, HRT is the most influential factor for the removal of pollutants in wetlands, since HRT determines the time required for pollutants to be adsorbed and degraded in wetlands. Therefore, the longer HRT of 1. 5 d was selected as the optimal level to ensure the best performance by CW-MFCs for treating pollutants.
section_name: Apparent Hydraulic Residence Time

section_content: As shown in Figure 5, when the anode size increased from B1 (0. 02 m) to B2 (0. 04 m), the average HRT (T m ) slightly increased from 0. 543 to 0. 556, respectively; then, it kept decreasing to 0. 450 as the anode size increased to B4 (0. 08 m) and increased to 0. 514 at B5 (0. 10 m). The changing trend of variance (σ 2 ) with the anode size was similar to that of the average HRT (T m ); however, the variance was much greater. The water flow divergence (σ 2 0 ) kept increasing from 0. 069 to 0. 213 as the anode size increased from B1 (0. 02 m) to B5 (0. 10 m), respectively. The effective volume ratio (e) kept decreasing from 0. 679 to 0. 593 as the anode size increased from B1 (0. 02 m) to B4 (0. 08 m), respectively, and then increased to 0. 623 at B5 (0. 10 m). The increase in variance (σ 2 ) and divergence (σ 2 0 ) resulted in the decrease in hydraulic efficiency (λ); thus, the change law of hydraulic efficiency (λ) with the anode size was opposite to that of water flow divergence (σ 2 0 ). The anodes occupied part of the space in the CW-MFC, and thus, they reduced the water-carrying section area, which led to an increase in the flow velocity after the water flowed around the bottom of the anodes. Meanwhile, boundary layer separation occurred in the water flow above the anodes due to the sudden expansion of the water-carrying section (Figure 4 ). This led to irregular fluctuations in the uniformity of the CW-MFC flow field with the anode size. Anodes with larger sizes occupied more space and further compressed the water flow channel, which caused the water flow pattern to be more mixed flow from push flow (Figure 5c ), and generally exacerbated the short-cut flow in CW-MFCs (Figure 5d ); finally, the hydraulic efficiency of CW-MFCs decreased as the anode size increased. 

Since the anode size was extremely significant for the hydraulic efficiency (λ) and significant for the effective volume ratio (e), B1 (0. 02 m), at which the highest values of λ and e occurred, was chosen to be the optimal level.
section_name: Anode Size

section_content: As shown in Figure 5, the average HRT (T m ) slightly decreased from 0. 516 at C1 (0. 1 m) to 0. 464 at C3 (0. 3 m); then, it increased to 0. 568 at C4 (0. 4 m), and decreased to 0. 544 at C5 (0. 5 m). The changing trend of the effective volume ratio (e) with the anode spacing was similar to that of the average HRT (T m ), and the highest value of e (0. 701) was achieved at C4 (0. 4 m). Variance (σ 2 ) decreased from 0. 224 to 0. 025 as the anode spacing increased from C1 (0. 1 m) to C3 (0. 3 m), respectively; then, it kept slightly increasing to 0. 084 at C5 (0. 5 m). Water flow divergence (σ 2 0 ) showed a similar changing trend with the anode spacing to that of variance (σ 2 ). The change law of hydraulic efficiency (λ) with the anode spacing was generally opposite to those of variance (σ 2 ) and divergence (σ 2 0 ), and the highest value of λ (0. 664) was achieved at C4 (0. 4 m). 

Anode spacing determined the distance between anodes, as well as the distance between the anodes and the wall of CW-MFC. When the anode spacing increased, the velocity of the water flow between the anodes slowed down since the water-carrying section area expanded (see Figure 4g,m), which prolonged the average HRT, enhanced the uniformity and stability of the flow field, and thus improved the hydraulic efficiency. However, the increase in anode spacing reduced the distance between the anode and wetland wall, and further, the water-carrying section area between the anode and wall reduced, and thus the flow velocity sped up, which shortened the average HRT, weakened the uniformity and stability of the flow field, and thus decreased the hydraulic efficiency. This trade-off in the flow field between anodes and between the anodes and the wall determined the variation in the hydraulic characteristic parameters of the CW-MFC with different anode spacings. From the higher values of T m, e, and λ and lower values of σ 2 and σ 2 0 at C4 (0. 4 m) and C5 (0. 5 m), it can be inferred that relatively sufficient anode spacing was beneficial to improve the overall hydraulic performance of CW-MFCs. 

Considering that the anode spacing was extremely significant for the hydraulic efficiency (λ) and significant for the water flow divergence (σ 2 0 ) and effective volume ratio (e), C4 (0. 4 m), at which the highest values of λ and e and lowest value of σ 2 0 occurred, was chosen to be the optimal level.
section_name: Anode Spacing

section_content: As shown in Figure 5, the values of T m, σ 2, σ 2 0, e, and λ at D1 (substrate surface) were higher than those at D2 (near plant roots). It was easily understood that the cathode plate located on the substrate surface hardly influenced the flow field inside the substrates of the CW-MFC, while the cathode plate placed near plant roots affected the hydraulic performance of CW-MFCs because of the blocking of the water flow, the flow around them, and boundary layer separation by the cathodes. However, the values of T m, e, and λ at D2 were lower than those at D1 only by 7. 6%, 8. 1%, and 5. 4%, respectively, showing relatively small differences. This probably because that the cathode near the plant's roots was located in the third layer of substrates from the bottom to higher up, and the affected water flow space was relatively limited. 

Since the cathode position was not significant for any of the hydraulic characteristic parameters, both D1 (substrate surface) and D2 (near plant roots) were selected as the optimal levels; however, D1 was recommended.
section_name: Cathode Position

section_content: Unsurprisingly, the larger the cathode size is, the greater the disturbance of the cathode in water is. Hence, as the cathode size increased, the values of σ 2 and σ 2 0 generally increased, and the values of e and λ generally decreased despite the fluctuations. However, almost no changes were found in the average HRT (T m ) with the variation in cathode size. Further, as the cathode size increased, the stability and uniformity of the flow field in CW-MFCs reduced, the mixed flow and short-cut flow were exacerbated, and the hydraulic efficiency of CW-MFCs decreased. However, since the water flow affected by the cathodes was relatively limited, the cathode size had no significant effect on the characteristics of the flow field in CW-MFCs. Therefore, from the perspective of hydraulics, there was no specific optimal level for the cathode size. 

In summary, MFC electrodes influenced the water flow velocity distribution and its stability and uniformity mainly through the blocking of the water flow, the flow around them, compressing water flow channels, and boundary layer separation, thus influencing the hydraulic efficiency of CW-MFCs. The comprehensive optimal electrode configurations in terms of the hydraulic performance of CW-MFCs were the apparent hydraulic residence time = 1. 5 d, anode size = 0. 02 m, anode spacing = 0. 4 m, cathode position = substrate surface or near plant roots, and no specific optimal level for cathode size (A5B1C4D1E or A5B1C4D2E).
section_name: Cathode Size

section_content: Confirmation tests were carried out to investigate the performance of hydraulics, electricity generation, and pollutants' removal by CW-MFCs under the optimal electrode configuration (A5B1C4D1E2). Since the hydraulic performance of CW-MFCs under the electrode configuration of B2C4D1E3 was the best among the 25 conditions in Table 4, CW-MFC performance under A5B2C4D1E3 was also investigated in comparison. The results are provided in Figure 6 and Table 6. performance under A5B2C4D1E3 was also investigated in comparison. The results are provided in Figure 6 and Table 6. As shown in Figure 6, in CW-MFC under A5B1C4D1E2, the values of σ 2 and σ 2 0 were much lower than those in CW-MFC under A5B2C4D1E3, while the values of Tm, e, and λ were higher than those in CW-MFC under A5B2C4D1E3 by 6. 7%, 6. 6%, and 8. 1%, respectively. This conclusively proved that the hydraulic performance of CW-MFC under A5B1C4D1E2 was better than that under A5B2C4D1E3. Additionally, the internal flow velocity in CW-MFC under A5B1C4D1E2 was 2. 6 × 10 -5 m/s-2. 5 × 10 -4 m/s, demonstrating less variation than that in CW-MFC under A5B2C4D1E3 (3. 2 × 10 -5 m/s-3. 2 × 10 -4 m/s). Moreover, the more concentrated and uniform streamlines visually demonstrated that the CW-MFC under A5B1C4D1E2 possessed a more uniform and stable flow field when it As shown in Figure 6, in CW-MFC under A5B1C4D1E2, the values of σ 2 and σ 2 0 were much lower than those in CW-MFC under A5B2C4D1E3, while the values of T m, e, and λ were higher than those in CW-MFC under A5B2C4D1E3 by 6. 7%, 6. 6%, and 8. 1%, respectively. This conclusively proved that the hydraulic performance of CW-MFC under A5B1C4D1E2 was better than that under A5B2C4D1E3. Additionally, the internal flow velocity in CW-MFC under A5B1C4D1E2 was 2. 6 × 10 -5 m/s-2. 5 × 10 -4 m/s, demonstrating less variation than that in CW-MFC under A5B2C4D1E3 (3. 2 × 10 -5 m/s-3. 2 × 10 -4 m/s). Moreover, the more concentrated and uniform streamlines visually demonstrated that the CW-MFC under A5B1C4D1E2 possessed a more uniform and stable flow field when it was compared to that of CW-MFC under A5B2C4D1E3. 

Electricity generation and pollutants' removal in CW-MFCs were influenced by various factors, such as the organic loading, redox conditions, types and growth of plants, and arrangement and materials of electrodes [7, 8]. According to previous studies, the output power density of CW-MFCs varied between 0. 2 mW/m 3 and 19. 6 W/m 3, and it was below 25 mW/m 3 in most cases [8, 23, 24], and the pollutants' removal in CW-MFCs was 44. 5-99. 0% for COD, 65. 9-96. 7% for TP, 36. 2-97. 0% for NH 3 -N, and 49. 7-99. 0% for TN [8, 25, 26]. Hence, in this study, the power density (2. 1-3. 4 mW/m 3 ) and removal rate of COD (82. 6-88. 2%), TP (74. 5-76. 6%), NH 3 -N (60. 2-66. 9%), and TN (64. 4-68. 6%) of CW-MFCs were in line with normal values. As shown in Table 6, the average output voltage, current, and power density in CW-MFC under A5B1C4D1E2 were higher than those in CW-MFC under A5B2C4D1E3 by 27. 3%, 27. 3%, and 62. 2%, respectively. The average removal rates of COD, TP, NH 3 -N, and TN in CW-MFC under A5B1C4D1E2 were higher than those of CW-MFC under A5B2C4D1E3 by 5. 6%, 2. 1%, 6. 7%, and 4. 2%, respectively. It can be concluded that optimizing the electrode configuration enabled the CW-MFC to obtain a more uniform and stable HRT distribution, thereby enhancing the contact between wastewater and the substrate and microorganisms in CW-MFC; thus, the adsorption and degradation of pollutants, and as a result, electricity generation and pollutants' removal in CW-MFCs were improved. Moreover, the improvements achieved in electricity generation can in turn facilitate the pollutants removal, especially the removal of organics and nitrogen, which is consistent with the results of previous research carried out by Wang et al. [27]. 

Overall, in comprehensively considering the hydraulic performance, electricity generation and pollutants' removal, the selected optimal electrode configurations (A5B1C4D1E or A5B1C4D2E) are suitable for the CW-MFCs.
section_name: Performance of CW-MFCs under the Optimal Electrode Configuration

section_content: The hydraulic characteristic was an important factor influencing the performance of CW-MFCs, since the water flow and HRT distribution of wastewater determines the contact time between the pollutants and the substrates, electrodes, and biofilms. This study first verified the effects of MFC electrode configurations on the flow field characteristics of CW-MFCs by CFD numerical simulation and obtained optimized values of electrode configurations. This proves the feasibility of optimizing the hydraulic performance of CW-MFCs by adjusting electrode configurations through CFD numerical simulation, and thus avoids the heavy workload of experiments. However, the flow field in CW-MFCs can be influenced by various factors, such as the properties of substrates, plant roots, type and arrangement of MFC electrodes, and system scale. The obtained results in this study might only work under certain restrictions, but they still could offer scientific reference for improving the hydraulic performance of CW-MFCs, and moreover, this study provided a new research perspective for improving the wastewater treatment and electricity production performance of CW-MFCs. The future research will mainly focus on: (1) combined effects of multiple factors on the hydraulic performance of CW-MFCs, and (2) the cost control of MFC electrodes, so as to provide a theoretical basis for regulating the flow field of CW-MFCs, and further, promote the practical application of CW-MFCs.
section_name: Practical Applications and Future Research Prospects

section_content: The apparent HRT (A) was the most influential and decisive factor, with a contribution of over 90% for the average HRT of CW-MFCs. Anode spacing (C) was the most influential factor for the hydraulic performance of CW-MFCs, with contributions of over 50% for water flow divergence (σ 2 0 ) and hydraulic efficiency (λ) and over 45% for the effective
section_name: Conclusions

section_content: Funding: This research was funded by the Hebei Province's Funding Program on Introducing Overseas Students (Grant No. C20200316 ) and the Natural Science Foundation of Hebei Province of China (Grant No. E2020202027 ). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable.

section_name: Data Availability Statement:

section_content: The authors confirm that the data supporting the findings of this study are available within the article.
section_name: Data Availability Statement:

section_name: Conflicts of Interest:

section_content: Author Contributions: Conceptualization, X. W. and C. Z. ; methodology, X. W. and C. Z. ; formal analysis, Y. Z., J. L., M. X. and W. X. ; investigation, Y. Z., J. L., M. X. and W. X. ; writing-original draft preparation, X. W. and Y. Z. ; writing-review and editing, X. W. ; supervision, X. W. ; funding acquisition, X. W. All authors have read and agreed to the published version of the manuscript.

section_content: The authors declare no conflict of interest.
section_name: Conflicts of Interest: