M4EKR, Multiphysics for ElectroKinetic Remediation of Polluted Soils

Electrokinetic soil remediation (EKR) is one of the most promising techniques for decontamination of low permeability soils, in which the most classical techniques have been found to be less efficient. However, its practical application on a real scale has been rather limited since the phenomena involved in these processes are very complex. For this reason, it is essential to use numerical models that allow us to know the main trends in the behaviour of soils and natural waters subjected to EKR processes. This study presents the numerical model M4EKR (Multiphysics for ElectroKinetic Remediation). The M4EKR module is a reactive transport model for partially saturated soils that allows reproducing the transport of species due to electroosmosis, electromigration, diffusion and advective flow. The model was completely implemented in COMSOL Multiphysics, a partial differential equation solver. The system of differential and algebraic equations to solve the chemical and transport problem was fully defined by the authors, and they were solved by the M4EKR module in COMSOL (monolithic approach). The scope of the model is illustrated by simulating an EKR process of a natural soil and porewater contaminated with a polar pesticide: 2,4-Dichlorophenoxyacetic acid. For simplicity, the M4EKR version used in this study does not include gas transport, it does not consider the deformability of the soil and it is assumed the processes occur under isothermal conditions.


Introduction
Global population growth has been accompanied by a significant geographical shift. By 2014, more than half of the world's population was already concentrated in cities, but by 2050 nearly two-thirds of the world's population will be living in urban areas [1]. This will mean the rapid growth of cities, especially in developing countries, but also in the most advanced countries. On many occasions, urban expansions will be made on old industrial spaces that in many cases present an important level of soil pollution. The remediation and restoration of these soils is undoubtedly a major challenge that human beings must face in the coming years.
Electrokinetic soil remediation is one of the most promising techniques for decontamination of low permeability soils, in which the most classical techniques [2] have been found to be less efficient. However, its practical application on a real scale has been rather limited [3] since the phenomena involved in these processes are very complex [4]. For this reason, it is essential to use numerical models that allow us to know the main trends in the behaviour of soils and natural waters subjected to EKR processes [5]. However, one of the most complex simulation tasks is the calculation of the multiple equilibria present in the geochemical systems of the soil and porewaters. For this reason, the use of reactive transport models coupled to electrokinetic models is essential [4,[6][7][8][9].
This study developed a model, M4EKR (Multiphysics for ElectroKinetic Remediation), that is applicable to electrokinetic remediation processes in natural unsaturated soils for the removal of ionic species present in porewater. An extended geochemical system with 36 species that are involved in 25 chemical equilibria resulting from the presence of the most common major ions in many natural groundwaters was studied. The model includes the transport of species due to electro-osmosis, electromigration, diffusion and advective flow in unsaturated conditions. This model can be used as a valuable tool to increase the knowledge of this treatment by a phenomenological understanding study. In this context, after describing the proposed formulation, the scope of the model is illustrated by simulating a synthetic EKR process of a natural soil and porewater contaminated with a polar pesticide: 2,4-Dichlorophenoxyacetic acid. The results show the descriptive capacity of the proposed model.

Conceptual model
The M4EKR module [5], which was used in this work, is a reactive transport model for partially saturated soils. For simplicity, the M4EKR version used in the present study does not include gas transport, and it does not consider the deformability of the soil. In addition, it is assumed that the processes occur under isothermal conditions (298.15 K).

Water mass balance
The equation used for the water mass balance is given by the expression where ∇ is the divergence operator, w l is the mass flow of water (kg m -2 s -1 ) and w m is the mass of water per unit total volume (kg m -3 ). The van Genuchten retention curve model [10] is used to calculate m w , and the soil porosity is considered constant [5]. The water mass flux ( w l ) is given by the expression  [11]. Electroosmotic permeability was defined as the product of the saturated electroosmotic permeability and a relative permeability function. A Brooks and Corey-type power function with an exponent of 3 was selected. The saturated electroosmotic permeability is assumed constant, independent of soil pH.

Reactive transport
The general geochemical system simulated in the M4EKR module considers a total of J components capable of generating a total of N chemical species during their combination. The total amounts of each component should be known through their corresponding mass balance equations, which have the following general form j j j where j m is the total mass of component j per unit total volume (mol m -3 ), j l is the total molar flux of component j (mol m -2 s -1 ), and j R is the rate of production or consumption of component j (mol m -3 s -1 ). However, it is not necessary to solve J partial differential equations but rather to find the total of J-2, as the total mass of a component of the system can be determined by fulfilling the electroneutrality condition, and the mass balance of the water species is solved by equation 1.
The total molar flux of component j is calculated as follows h eo em dif j j j j j = + + + l l l l l (4) where h j l , eo j l and em j l are the advective flows generated by the unsaturated hydraulic flow, the electroosmotic flow and the electromigration of ions, respectively, and dif j l is the Fickian diffusive-dispersive flow defined as The chemical speciation problem is solved using a classical stoichiometric approach by solving a system of mass-balance and mass-action equations [12]. The mathematical formulation of this approach has been widely described in the literature [13,14] and has been applied in M4EKR [5]. WATEQ Debye Hückel formulation has been used to estimate the activity coefficients.

Electric charge balance
The balance equation of the total electric charge, assuming electroneutrality throughout the system and that there is not charge accumulation capacity, is given by the following expression where i is the total current density (A m -2 ), which is calculated by applying Ohm's law [5,15] and using the empirical Rhoades formulation [16] for the apparent electrical conductivity of the soil and the formulation proposed by Appelo [17] for the electrical conductivity of porewaters.

Mass balance in electrolyte wells
The electrolyte wells are necessary to the EKR processes for adding the washing fluid and for removing the collected pollutants. In the present study, two wells in which the electrodes are placed are considered. In addition, since the geochemical model used here has a total of N species generated through the combination of J components [14] and it has been assumed that the electrolyte compartments behave as ideal continuous stirred tank reactors, J-2 balance equations are used to determine the component mass evolution (the same as in the soil). Each electrolyte compartment has the following ordinary differential equation for the mass balance of each component (10) where the * superscript is equal to A for the anode or C for the cathode, depending on the affected well, j R is the production or consumption flow of the component (mol s -1 ) according to the electrochemical reactions, j & M is the total mass of component j as expressed in moles and represents the input (superscript in) or output (superscript out) for the mass flows of component j (mol s -1 ). These flows, as well as j R and j & M , are defined in a previous study in which the complete formulation of the numerical model is described in detail [5]. It is interesting to note that the output of the pollutant from the system only occurs by overflow in the well where the cathode is located, and it moves towards the direction in which the electroosmotic water flow is moving.

Numerical model
The M4EKR model was completely implemented in the partial differential equation solver known as COMSOL Multiphysics [18]. This software is based on the application of the finite element method with Lagrange multipliers. The system of differential and algebraic equations to solve was fully defined by the authors. This approach is possible due to the versatility and adaptability of this type of program because of its multiphysics environment [19]. In addition, the automatic differentiation techniques [20] included in some programs, such as COMSOL Multiphysics, provide symbolic expressions for defining the iteration matrix [21], which improves the convergence behaviour of the model and allows the complexity of the problem to increase while the solution efficiency is maintained [22,23]. One of the differentiating features of the code is that it includes the mass-balance and mass-action algebraic equations of the stoichiometric approach. This approach is the most common one in geochemical codes [14,24,25], some of which, such as the PHREEQC program [13], are coupled to other modules to create reactive transport codes, which use common operatorsplitting procedures [26][27][28]. However, none of these procedures is used for speciation in the M4EKR module, but the system of equations in the geochemical system is solved by COMSOL Multiphysics (Monolithic approach). In this way, a large number of state variables must be solved in a coupled manner at each time step and in all cases. Finally, for more clarity regarding the analysis, a one-dimensional configuration is used.

Simulation
The electrokinetic remediation of a soil polluted with 2,4-D pesticide at an initial concentration of 20 mg/kg drysoil applying an electrical potential gradient of 1 V cm -1 has been studied. The mineralogical, textural and index properties of the simulated soil are described by López-Vizcaíno et al. [5], as well as the unsaturated hydraulic properties, which correspond to a natural soil taken from a clay quarry in Mora, Toledo, Spain. Figure 1 shows the dimensions of the experimental setup. The electrokinetic reactor is composed of a 15 × 15 × 3 cm central compartment that houses the soil to be treated and two side compartments with dimensions of 15 × 3 × 3 cm that are used as electrolyte wells, in which the electrodes are located. The electrodes have a prismatic geometry with an electrochemically active surface equal to the well/soil interface (S A =S C =45 cm 2 ). The electrolyte volume in the wells is assumed to be constant throughout the simulation process (V A =V C =135 cm 3 ).  Table 1 and 2 shows the components and species, respectively, that make up the geochemical model proposed.
The equilibrium reactions, their thermodynamic constants [29] and the physicochemical properties of each species (hard core diameter and binary diffusion coefficients in water) are also shown. The values of the binary diffusion coefficients in water of each species were not found in the literature, so they were estimated using the Pikal's expression [30].  Table 3 shows the chemical initial conditions of the problems simulated in this study. Soil is assumed to be saturated. The compositions of the electrolyte and porewater were defined from a local sample of groundwater with a significant carbonate concentration. Figure 1 shows the boundary conditions applied in the simulations.

Results and discussion
In order to assess the functionality of the M4EKR model, the EKR process proposed was simulated for a total test time of 3 days. Figure 2 shows the spatial distribution of a set of variables of interest for different times. The pH distribution follows the typical pattern observed in literature [31], including an acidic front and a basic front that move towards the centre of the domain from the anodic and cathodic wells, respectively (Figures 2a and   2b). Both are located in a central position (although closer to the cathode due to the higher ionic mobility of H + ), and the resulting front moves slowly to the position of the cathode well. The behaviours of representative ions (Na+ and Ca 2+ as cations and SO 4 2as anion) are similar to what is expected theoretically (Figures 2c to  2e). Thus, increased concentrations of the anionic and cationic components form near the anode and the cathode, respectively, which causes long-term increases in the ionic strengths with respect to the initial values in these zones (Figure 2f). Furthermore, a progressive desaturation of the soil sample takes place (Figures 2g  and 2h) as the electrokinetic remediation process progresses as was observed by Airoldi et al. [31]. This can lead to cracking and other geotechnical problems such as those described by Lopez-Vizcaíno et al. [32]. The spatial distribution of pH determines the speciation of some components such as carbonates and especially the pollutant substance: the 2,4-D pesticide. As can be seen in Figure 3, only acidic (and electrically neutral) species are found behind the acid front. This fact conditions how the pesticide or some other species move in the soil. For example, the pesticide in acid form can only be displaced by hydraulic and electroosmotic flow (mainly towards the cathode), but not by electromigration. Consequently, an accurate calculation of the pH front position is necessary to understand the phenomena that occur during an EKR process.
Numerical models, such as M4EKR, allow knowing the geochemical effects that control the electrokinetic soil remediation processes. Although the example presented is based on a polar-type pesticide, these effects are much more important in the case of heavy metals that can precipitate forming mineral phases difficult to extract from the soil by means of the most common remediation techniques.

Conclusion
There is a very wide range of physical and chemical phenomena in EKR processes that overlap and, depending on their magnitude, can produce very different results. The use of numerical models is essential for interpreting laboratory tests and for optimising the planning of EKR processes in real problems. This paper presents the M4EKR model, which is capable of simulating both the geochemical system and the electrokinetic effects that occur during an EKR process. The most important effects of the displacement of the pH front between anode and cathode have been identified. These are fundamentally associated with the properties of the pesticide studied. The proposed numerical model, therefore, allows the development of operation strategies for soil decontamination in the most efficient way possible, based on an understanding of the physical-chemical phenomena involved and their interaction.