Bioremediation of Vanadium from Contaminated Water in Bioreactor Using Methylocystis hirsuta Bacterium: Comparisons with In Silico 2D and 3D Simulations

: The elimination of poisonous wastes (e.g., heavy metals) from polluted water remains challenging, both in industrialized societies and developing countries. To overcome this human health and environmental issue, biotechnology (e.g., biosorption, bioaccumulation) is being applied as an economic and eco-friendly option compared to physicochemical methods (e.g., adsorption, membrane ﬁltration, and coagulation–ﬂocculation). The development of the appropriate biotechnology process (i.e., bioremediation) requires more accurate information and details, which are possible to obtain through the design of a set of resources and various computer applications. In sustainable remediation, microorganisms are one of the feasible choices for modifying and remaking the natural condition. In this in silico study, the methanotroph Methylocystis hirsuta (M. hirsuta ) was used for the ﬁrst time to simulate the removal of vanadium (Vn) from contaminated water through two-dimensional (2D) and three-dimensional (3D) modeling using COMSOL 4.4 software. Rotating machinery-laminar ﬂow, transport of diluted species, and reaction engineering physics were also used. Independency analyses of the numerical network, concentration contour, velocity contour, concentration–time, and velocity–distance charts were also calculated. The data consistently showed that the removal of Vn increased with increasing velocity (which depends on time). Indeed, the amount of pollutant removal at 120 rpm, 160 rpm, and 200 rpm was maintained at 10%, 12%, and 12%, respectively. The simulation results showed excellent conformity (less than 20%) with previously reported laboratory results. This proposed model of bioremediation is thus a reliable and accurate solution for the removal of heavy metals (i.e., Vn and possibly others) from polluted areas (such as contaminated water).


Introduction
It is worth mentioning that industries producing heavy metals most often discharge contaminated wastewater with large amounts of metals such as cadmium (Cd), chrome (Cr), copper (Cu), nickel (Ni), arsenic (As), lead (Pb), vanadium (Vn), and zinc (Zn). De facto, these industries represent the most dangerous chemical-intensive industries [1][2][3]. Cleaning up the contaminated water and soil from heavy metals is challenging because they are nonbiodegradable [3,4]. Since many natural environments worldwide are contaminated by heavy metallic waste, the development of cost-effective technologies to remove this waste from polluted water and soil has become urgent and important [3,5,6].
The different physicochemical methods currently applied to remove heavy metals from wastewater, such as chemical precipitation, ion exchange, adsorption, membrane filtration, coagulation-flocculation, flotation, and electrochemical methods have been largely described [3,7,8]. Herein, we report the application of biological and computerbased methods. If the concentration of heavy metals is extremely low, physicochemical methods (e.g., remediation by solvent, adsorption, and absorption by some materials) are unsuitable or costly. In this situation, to remove heavy metals, biological methods such as biosorption and bioaccumulation are noteworthy options [9]. In sustainable remediation, using microorganisms and plants is one of the feasible choices for modifying and remaking the natural condition of the soil. Indeed, although high levels of heavy metals in the soil lead to disadvantageous modification of the microbial community, low concentrations are crucial for microorganism growth [10].
The calculated use of biological degradation methods (such as plants, microorganisms, and microbial processes or products) to remove, degrade, or decrease the concentration of environmental pollutants/contaminants from places where they have been released is called bioremediation [3,11,12]. Pollutant concentrations are decreased to acceptable levels in the areas owned by farmers and/or regulatory agencies [13][14][15][16]. Bioremediation is most often used around the world to create healthier soils, to decontaminate sediments, to clean up water sources (e.g., groundwater aquifers), and to improve air quality [16]. Usually, the control measures required to reach optimal process conditions are evaluated [17,18]. Some of the advantages and disadvantages of using bioremediation processes in comparison with other remediation technologies can be inserted. Therefore, unlike the disruptive excavationbased remediation processes, bioremediation is less intrusive and facilitates remediation of environmental impacts without damaging delicate ecosystems [16,[19][20][21]. Susceptibility of bioremediation to environmental changes, such as differences in available oxygen, pH, or temperature [22], must be considered. Interestingly, bioremediation is much less expensive for treating a hazardous waste compared to that of conventional treatment methods such as vacuuming, absorbing, burning, dispersing, or moving/transferring contaminants from one environmental medium to another [14,16,23]. Bioremediation often requires very long periods of time [24,25]. The products of biodegradation may be more persistent or toxic than the parent compound [26,27].
To remove Cu, Cd, and Pb, some bacteria, including Bacillus sp., Pseudomonas sp., Micrococcus sp., Thiobacillus ferroxidans, and leptospirillum ferroxidans, are used to oxidize iron and sulfur [28]. In a study to investigate biosorption characteristics of Cr (VI), Cu (II), and Fe (III) ions, the living biomass of Bacillus licheniformis was used [29]. Vn is used widely in steelmaking as an industrially vital element. A high concentration of Vn in the human body causes dizziness, vomiting, nausea, and kidney damages; its pollution in water is a global issue [30]. Nowadays, to remove sulfide from streams, chemical methods are switched to biological methods, including the use of sulfur oxidizing bacteria (SOB). Due to the presence of different types of bacteria and diversity in their capabilities, the abilities of seven well-known SOBs were studied [31].
Herein, the methanotroph Methylocystis hirsuta was used to remove Vn. M. hirsuta, isolated from a groundwater aquifer, is an anaerobic and Gram-negative bacterium with single or cluster cells of relatively small size (0.3-0.6 × 0.7-1 µm). M. hirsuta reproduces by normal cell division; budding division does not occur. Cells are not motile but possess a spiny surface layer composed of polysaccharides [23]. In this study, a one-dimensional (1D) reactive transport model that was time-dependent was used. Several biogeochemical processes were simulated during the trace-metal-contaminated bioremediation of aquifers. The numerical solution of equations was determined by finite difference approximation. Before designing uranium bioremediation, some key processes must be simulated to achieve the highest accuracy [32].
In this study, and to the best of our knowledge, the Vn reduction through 2D and 3D simulations using Comsol Multiphysics software was reported for the first time. The rotating machinery-laminar flow, transport of diluted species, and reaction engineering physics were used. The in silico results were compared with experimental data. Independency anal-yses of the numerical network, concentration contour, velocity contour, concentration-time chart, and velocity-distance chart were also calculated.

Material and Methods
In this simulation, the M. hirsuta performance on the removal of biological contaminants was simulated. Two geometries were designed in 2D and 3D using the Comsol Multiphysics 4.4 software. Several simplifying assumptions, boundary conditions, and meshing geometry were used in the system.

Simplifying assumptions
The assumptions for the simulation were the following: • The culture medium with microorganisms and the inside metal was considered a single-phase.

•
The vanadium concentration and microorganisms were spread homogeneously throughout the system. • The physical properties of the solution changing with time were ignored.

•
The system temperature was constant and set at 25 • C [33,34].

Model Geometry
The model comprised two parts (i.e., a vessel and an impeller with 4 blades with pitch angel of 45 • ) and was made by employing 2D and 3D COMSOL software. Figure 1 shows the 2D and 3D structures of the model. In a 2D structure, the vessel diameter and the impeller diameter were assumed at 0.2 m and 0.08 m, respectively. In 3D, the vessel diameter, height, and the impeller diameter were assumed at 0.2 m, 0.3 m, and 0.03 m, respectively. physics were used. The in silico results were compared with experimental data. Independency analyses of the numerical network, concentration contour, velocity contour, concentration-time chart, and velocity-distance chart were also calculated.

Material and Methods
In this simulation, the M. hirsuta performance on the removal of biological contaminants was simulated. Two geometries were designed in 2D and 3D using the Comsol Multiphysics 4.4 software. Several simplifying assumptions, boundary conditions, and meshing geometry were used in the system.

Simplifying assumptions
The assumptions for the simulation were the following: • The culture medium with microorganisms and the inside metal was considered a single-phase.

•
The vanadium concentration and microorganisms were spread homogeneously throughout the system. • The physical properties of the solution changing with time were ignored.

•
The system temperature was constant and set at 25 °C [33,34].

Model Geometry
The model comprised two parts (i.e., a vessel and an impeller with 4 blades with pitch angel of 45°) and was made by employing 2D and 3D COMSOL software. Figure 1 shows the 2D and 3D structures of the model. In a 2D structure, the vessel diameter and the impeller diameter were assumed at 0.2 m and 0.08 m, respectively. In 3D, the vessel diameter, height, and the impeller diameter were assumed at 0.2 m, 0.3 m, and 0.03 m, respectively. The whole model in 2D was meshed as a free triangle and had a finer size, and it was divided into 9630 elements. In contrast, the model in 3D was meshed as a free tetrahedron and was divided into 38,091 elements.
Governing Equation COMSOL has several types of physics for different simulations [35]. In the current simulation, three physics were used: (i) rotating machinery-laminar flow, (ii) transport of diluted species, and (iii) reaction engineering.
Rotating machinery-laminar flow Due to the Reynolds number, the flow is laminar, and because of the existence of the impeller in bioreactor rotating machinery, laminar flow physics was selected. The whole model in 2D was meshed as a free triangle and had a finer size, and it was divided into 9630 elements. In contrast, the model in 3D was meshed as a free tetrahedron and was divided into 38,091 elements.
Governing Equation COMSOL has several types of physics for different simulations [35]. In the current simulation, three physics were used: (i) rotating machinery-laminar flow, (ii) transport of diluted species, and (iii) reaction engineering.
Rotating machinery-laminar flow Due to the Reynolds number, the flow is laminar, and because of the existence of the impeller in bioreactor rotating machinery, laminar flow physics was selected.

Transport of diluted species
The transport of diluted species interface is used to compute the concentration field of a dilute solute in a solvent. The transport and reactions of the species dissolved in gas, liquid, or solid can be computed. The driving forces for transport can be diffusion by Fick's law, convection when coupled to fluid flow, and migration when coupled to an electric field.
Depending on the licensed products, modeling multiple species transport is possible. Fick's law (Equation (3)) is used in this physics [35].
where C, D, R, u, and t are the concentration of Vn (mol/m 3 ), the diffusivity coefficient (m 2 /s), reaction rate (mol/m 3 ·s), the velocity vector (m/s), and time (s), respectively. The main equation that describes the solute transfer from aqueous phase is the continuity equation. This equation is also called the convection and diffusion equation [34]. The continuity equation should be solved numerically to obtain the concentration's solute distribution. To solve the continuity equation, the velocity distribution is needed. Velocity distribution is obtained by solving the momentum equation, i.e., Navier-Stokes equations. Therefore, the momentum and continuity equations should be coupled and solved simultaneously to obtain the concentration distribution of solute [34].

Reaction engineering
The reactor type in which the reaction occurs is very important [35]. The design of an industrial chemical reactor must satisfy the requirements, which are the chemical factors, the mass transfer factors, the heat transfer factors, the safety factors, and economic factors [36]. The reactor is a batch, and the reaction is of the first order and is irreversible by governing the Equation (4): where c, t, r, and ν i are concentration (mol/m 3 ), time (s), rate (mol/m 3 ·s), and the Stoichiometric coefficients, respectively [35].

Independency Analysis of Numerical Network
To analyze the numerical network, the geometry was meshed into three types. Table 1 shows the number of mesh elements in various mesh types. The final concentration of Vn (in the last time step) was checked. There was no significant difference in the final concentration of Vn, and independency of the numerical network was established. Finally, the fine mesh was selected.     Figure 2 shows the concentration contour of Vn at 120, 160, and 200 rpm. As can be seen from Figure 2, the Vn concentration gradient is significant at the regions near the center of the bioreactor. From the center to the side, the concentration of Velocity (rpm)

Concentration Contour
Vn is removed during the reaction of Va→W in the presence of Methylocystis hir bacteria. To investigate the concentration of Vn, contours of Vn concentration were ca lated. After simulation, concentration contour was defined in 21 days. Figure 2 shows concentration contour of Vn at 120, 160, and 200 rpm. As can be seen from Figure 2, the Vn concentration gradient is significant at the gions near the center of the bioreactor. From the center to the side, the concentration  Figure 2 shows the concentration contour of Vn at 120, 160, and 200 rpm.  Table 2 shows the final concentration of vanadium in the last time step.  Figure 2 shows the concentration contour of Vn at 120, 160, and 200 rpm. As can be seen from Figure 2, the Vn concentration gradient is significant at the regions near the center of the bioreactor. From the center to the side, the concentration of Vn decreased. High mass transfer flux occurs in these regions so most elimination occurred near the wall. The more the velocity increased, the smaller the area was around the bioreactor center (which has the maximum concentration). Additionally, the best removal As can be seen from Figure 2, the Vn concentration gradient is significant at the regions near the center of the bioreactor. From the center to the side, the concentration of Vn decreased. High mass transfer flux occurs in these regions so most elimination occurred near the wall. The more the velocity increased, the smaller the area was around the bioreactor center (which has the maximum concentration). Additionally, the best removal rate happens within a shorter distance from the center. This result means that an increment in velocity makes the final removal rate happen in a shorter space.

Velocity Contour
The velocity field in the bioreactor was simulated by solving the Navier-Stokes's equations. Figure 3 shows the result of spatial velocity variation in the bioreactor. The right side of the figure shows the assigned speed for each color.
Sustainability 2022, 14, 8807 6 of 10 rate happens within a shorter distance from the center. This result means that an increment in velocity makes the final removal rate happen in a shorter space.

Velocity Contour
The velocity field in the bioreactor was simulated by solving the Navier-Stokes's equations. Figure 3 shows the result of spatial velocity variation in the bioreactor. The right side of the figure shows the assigned speed for each color.

Velocity-Distance Chart
The model's equations with the appropriate boundary conditions were solved numerically using COMSOL software. This software employs the finite element method (FEM) for numerical solutions of model equations [34].
Momentum simulation was performed at steady state; for this reason, the velocitydistance chart was examined (Figure 4). As represented in Figure 4, the velocity of the center is zero. By gradually increasing the distance from the center, the velocity chart increased with increasing slope, and by the

Velocity-Distance Chart
The model's equations with the appropriate boundary conditions were solved numerically using COMSOL software. This software employs the finite element method (FEM) for numerical solutions of model equations [34].
Momentum simulation was performed at steady state; for this reason, the velocitydistance chart was examined (Figure 4). ment in velocity makes the final removal rate happen in a shorter space.

Velocity Contour
The velocity field in the bioreactor was simulated by solving the Navier-Stokes's equations. Figure 3 shows the result of spatial velocity variation in the bioreactor. The right side of the figure shows the assigned speed for each color.

Velocity-Distance Chart
The model's equations with the appropriate boundary conditions were solved numerically using COMSOL software. This software employs the finite element method (FEM) for numerical solutions of model equations [34].
Momentum simulation was performed at steady state; for this reason, the velocitydistance chart was examined (Figure 4). As represented in Figure 4, the velocity of the center is zero. By gradually increasing the distance from the center, the velocity chart increased with increasing slope, and by the As represented in Figure 4, the velocity of the center is zero. By gradually increasing the distance from the center, the velocity chart increased with increasing slope, and by the end of the impeller, it reached its maximum amount. The velocity then decreased until it reached zero on the side. The velocity profile is almost parabolic in the last region.

Vanadium Concentration vs. Time
The concentration-time changes in Vn were obtained. The results at three velocities (i.e., 120 rpm, 160 rpm, and 200 rpm) are illustrated. At the beginning of the reaction (at t = 0), the concentration of Vn is the highest and equals 0.98 mol/m 3 . At 21 days after the initial 2D simulation, this rate was decreased to 86% at 120 rpm and to 81% at 160 rpm and 200 rpm.

Concentration Contour
The simulation in 3D was more precise than that in 2D. In comparison with the 2D simulation, the final removal rate was achieved at a shorter distance from the center, as represented in Figure 5.

Vanadium Concentration vs. Time
The concentration-time changes in Vn were obtained. The results at three velocities (i.e., 120 rpm, 160 rpm, and 200 rpm) are illustrated. At the beginning of the reaction (at t = 0), the concentration of Vn is the highest and equals 0.98 mol/m 3 . At 21 days after the initial 2D simulation, this rate was decreased to 86% at 120 rpm and to 81% at 160 rpm and 200 rpm.

Concentration Contour
The simulation in 3D was more precise than that in 2D. In comparison with the 2D simulation, the final removal rate was achieved at a shorter distance from the center, as represented in Figure 5.

Velocity Contour
The local variations in velocity resulting from the bioreactor at 120, 160, and 200 rpm are depicted in Figure 6.

Velocity Contour
The local variations in velocity resulting from the bioreactor at 120, 160, and 200 rpm are depicted in Figure 6.

Vanadium Concentration vs. Time
After the 21-day 3D simulation, the Vn concentration was reduced to 82% at 120 rpm, to 76% at 160 rpm, and to 74% at 200 rpm. Table 3 shows the comparison of the elimination rates in the 2D and 3D simulation cases. According to the table, in the 2D simulation, the Vn removal amounts at 160 rpm and 200 rpm were the same but the amount was lower at 120 rpm. In comparison, 200 rpm was the best parameter for Vn removal throughout the 3D simulation.

Vanadium Concentration vs. Time
After the 21-day 3D simulation, the Vn concentration was reduced to 82% at 120 rpm, to 76% at 160 rpm, and to 74% at 200 rpm. Table 3 shows the comparison of the elimination rates in the 2D and 3D simulation cases. According to the table, in the 2D simulation, the Vn removal amounts at 160 rpm and 200 rpm were the same but the amount was lower at 120 rpm. In comparison, 200 rpm was the best parameter for Vn removal throughout the 3D simulation. In this study, COMSOL was used to simulate a 1D equilibrium transport lead-phenol binary system. The saturated sandy soil was considered as an aquifer and granular dead anaerobic sludge (GDAS) [37]. In order to describe absorption data in a binary system, three isothermal models were applied: the combined Langmuir-Freundlich model, the Extended Langmuir model, and the Extended Freundlich model [37]. We observed that the root mean squared error between the numerical and experimental results did not exceed 0.04, which is in line with previously published data [37]. In another study, three similar isothermal models were applied.
Taken together, the simulation data presented in this work compared with experimental data [38][39][40], including the experimental data on Vn previously reported by our team [40], show consistency in the application of our model for the bioremediation of deleterious heavy metals, such as Vn, found in contaminated water.

Conclusions
Bioremediation is a technology where biological organisms are used to remove or neutralize an environmental pollutant via a metabolic process. Herein, the COMSOL 2D and 3D simulations along with Fick's law and Navier-Stokes equations were used, for the first time, to understand how the elimination process of Vn from contaminated water occurs in the presence of M. hirsuta. Concentration-time chart, velocity-distance chart, concentration contour, and velocity contour were investigated. At higher velocity and time increments, the connection time between the microorganism (M. hirsuta) and the metal (Vn) increased and led to an increase in bioremoval by the microorganism. The simulation data match quite well with previously reported experimental data. It confirms that the proposed model is reliable and accurate and can be applied as a sustainable bioremediation option for the removal of Vn pollutants and possibly others to be tested.