Impact of edetate sodium on calcium carbonate scaling: Insights from molecular dynamics simulations

Abstract

In order to mitigate calcium carbonate (CaCO₃) scaling in oil and gas wells, antiscaling agents are commonly used. This research employed the molecular dynamics method to simulate the crystallization process of a supersaturated CaCO₃ solution, with particular emphasis on examining the influence of the antiscalant edetate sodium (EDTA). Analysis of the spatial distribution of Ca²⁺ in proximity to CO₃²⁻ revealed that EDTA effectively impedes the interaction between Ca²⁺ and CO₃²⁻, resulting in a reduction in the mean square displacement and diffusion coefficient of CO₃²⁻. The examination of intermolecular interaction energies indicates that the binding energy between Ca²⁺ and CO₃²⁻ is approximately −516 kcal mol⁻¹, while the binding energy between Ca²⁺ and EDTA is estimated at −718 kcal mol⁻¹. These findings suggest that EDTA exhibits a higher affinity for binding with Ca²⁺ compared to CO₃²⁻, consequently hindering the formation of CaCO₃.

Keywords

energy decomposition molecular simulation preventing scaling scale inhibitor

Schematic diagram of EDTA molecular complexation with calcium ions

Introduction

In Moxi Gas Field in Sichuan, the Longwangmiao formation gas wells exhibit a high salinity water production capacity, and multiple scaling ions can be detected in the formation water.¹ The occurrence of scale blockage in oil pipes, primarily in the form of calcium carbonate (CaCO₃) scale, frequently leads to a notable reduction in gas well productivity,^2,3 and the application of hydrochloric acid is the most common strategy for removing carbonate deposits.^4,5 However, implementation of this strategy results in increased levels of carbon dioxide (CO₂) and hydrogen sulfide (H₂S) in the gas, along with significant corrosion to metal pipes.⁶ In addition, frequent blockage of wellbores results in costly removal procedures, as well as production interruptions. To control blockages caused by the CaCO₃ scaling, the addition of scale inhibitors to aqueous solutions is often used.^7,8

Phosphate is a kind of high-performance scale inhibitor.⁹ However, due to its potential environmental harm, research is increasingly concentrated on the development of environmentally friendly scale inhibitors. In the work by Chen et al.,¹⁰ a novel and highly efficient polyaspartic acid/graphene oxide grafted copolymer (PASP/GO) was prepared. They evaluated the scale inhibition performance of PASP/GO and studied the mechanism by using the molecular dynamics (MD) simulation method. According to the results, PASP/GO exhibited excellent inhibition performance and dispersion capacity, and PASP/GO adsorbs scale particles by chelation, crystal distortion, and dispersion. In the work by Zuo et al.,¹¹ they studied the performance of several common scale inhibitors by using MD methods and demonstrated that 2-phosphanobutane-1,2,4-tricarboxylic acid (PBTCA) has the best scale inhibition performance.

In the work by Zeng et al.,¹² MD simulation was conducted to investigate the interaction between various configurations of poly(aspartic acid) (PASP) and different surfaces of CaCO₃ crystals in aqueous solution. The findings indicate that PASP effectively inhibits the growth of CaCO₃ scale, with the binding energies primarily arising from coulombic interactions, including ionic bonds. In the work by Zhang et al.,¹³ experimental and MD simulation methodologies were both used to meticulously investigate the scale inhibitory properties of carboxymethyl dextran (CMD) on the nucleation and crystallization kinetics of CaCO₃. The findings revealed significant interactions between Ca²⁺ in solution and the carboxyl groups of CMD. As the concentration of CMD increased, the likelihood and strength of adsorption of Ca²⁺ and CO₃²⁻ onto the solution’s surface decreased.

Due to environmental concerns, phosphorus-based scale inhibitors are not included in our study. Edetate sodium (EDTA) serves as a scale inhibitor whose functional groups can adsorb Ca²⁺ in industrial circulating water.^2,14,15 To investigate the fundamental principles underlying the function of EDTA as a scale inhibitor, this study examines the ion aggregation in CaCO₃ solution through MD simulation with and without EDTA; additionally, weak interactions between ions were also evaluated. The findings suggest that the interaction energy between EDTA and Ca²⁺ surpasses that of CO₃²⁻ and Ca²⁺, indicating a higher likelihood of EDTA binding with Ca²⁺ and preventing the precipitation of CaCO₃.

Computational details

In this study, GROMACS software version 2022.5^16,17 was employed to simulate the aggregation of Ca²⁺ and CO₃²⁻ ions with and without the scale inhibitor EDTA, the molecular structure formula is shown in Figure 1. The General AMBER Force Field (GAFF) was used for the characterization of small molecules,¹⁸ with atomic charges being assigned as the advanced restrained electrostatic potential (RESP2) charges.¹⁹ The ion parameters were obtained from existing literature sources,^20,21 and OPC3 water model²² was used to describe water. Packmol software²³ was used to build the simulation box with initial size of 10 nm × 10 nm × 10 nm. Within the control group, the concentration of the scale inhibitor EDTA is approximately 3%, the number of molecules and atoms present in the two systems are listed in Table 1. During the simulation, the temperature was set at 320 K and the pressure at 10 MPa in order to mimic the conditions in the wellbore. Following the minimization of system energy utilizing the steepest descent algorithm, leapfrog algorithm was used to perform 100 ns MD simulation under the canonical ensemble with constant pressure and temperature (NPT) ensemble; Parrinello-Rahman algorithm²⁴ with pressure coupling constant time of 2.5 ps and the V-rescale thermostat algorithm²⁵ with coupling constant time of 2 ps were both used to control the pressure and temperature within the simulation system. The Coulomb interaction calculation method was particle-mesh Ewald (PME) algorithm²⁶ and the van der Waals interaction algorithm was cut-off; the cutoff were both 1 nm. Visual molecular dynamics (VMD 1.9.3) software²⁷ was used for visual observation.

Figure 1.

Molecular structural formula of EDTA (structural and ball-stick model).

Table 1.

The number of molecules and atoms in the simulation box.

System	Ca²⁺	CO₃²⁻	EDTA	Na⁺	H₂O
A, without scale inhibitor	100	100	0	0	33,509
B, with 3% scale inhibitor EDTA	100	100	50	200	32,781

Independent gradient model based on Hirshfeld partition (IGMH) analysis²⁸ steps: first, select all atoms/molecules within 3.5 nm radius of a specified atom/molecule in VMD, followed by the utilization of ORCA 5.0.4^29–31 for the acquisition of wavefunction data at the B3LPY D3 ma-def2-TZVP level. Subsequently, Multiwfn³² was used for IGMH analysis based on the obtained wavefunction information.

Results and discussion

Figure 2(a) and (c) depict the initial state prior to MD simulation with and without the scale inhibitor EDTA, respectively, showing all molecules/atoms in a disordered arrangement. In contrast, Figure 2(b) and (d) illustrate the equilibrium state, where molecules/atoms aggregated together due to electrostatic attraction between anions and cations. The trajectory videos can be accessed at the following URLs https://youtu.be/IUDTrv2TSDk and https://youtu.be/dsZkiOZBqAM. There is a notable disparity in the distribution of molecules/atoms before and after simulation, with the aggregation of molecules/atoms after equilibrium showing an area of interest for further investigation.

Figure 2.

Screenshots before and after MD simulation (: Ca²⁺; :CO₃²⁺; : EDTA; : Na⁺).

The radial distribution function (RDF) is capable of characterizing the spatial arrangement of Ca²⁺ and CO₃²⁺. It indicates the position of a given particle $α$ , the distance of particle $β$ from particle $α$ (the probability of distribution around $α$ ), that can represent the order degree of the two components

g_{α β} (r) = \frac{V_{s}}{N_{α} N_{β}} \sum_{i = 1}^{N_{α}} \frac{n_{i β} (r)}{4 π r^{2} Δ r}

(1)

where $g_{α β} (r)$ is the average probability of particle $β$ appearing within the distance of $r$ and $r + Δ r$ away from the $i$ th particle $α$ . $N_{α}$ and $N_{α}$ are the number of $α$ and $β$ particles, respectively. $V_{s}$ is the volume of the simulation box; $r$ is the distance between two kinds of particles; $n_{i β} (r)$ is the number of $β$ particles within the distance of $r$ and $r + Δ r$ away from the $i$ th particle $α$ .

Figure 3 shows the RDF analysis results. The peak occurrence of Ca²⁺ ions at a distance of 0.23 nm from CO₃²⁻ ions suggests the formation of the first coordination layer. A minor peak is observed at approximately 0.38 nm, which is believed to be associated with the second coordination layer. From the perspective of peak intensity, the second coordination layer can be ignored. This observation suggests that CO₃²⁻ and Ca²⁺ ions exhibit a tendency to form pairs due to the opposing charges attraction, but not all anions are capable of pairing with cations. In contrast, the presence of EDTA in the solution leads to a significant decrease in the RDF value between Ca²⁺ and CO₃²⁻, suggesting that the interaction between Ca²⁺ and CO₃²⁻ is impeded, thereby inhibiting the binding of Ca²⁺ and CO₃²⁻. This phenomenon also elucidates the mechanism by which scale inhibitors function.

Figure 3.

RDF of Ca²⁺ and CO₃²⁻ with and without EDTA.

Based on Figure 3, it is evident that Ca²⁺ and CO₃²⁻ exhibit the greatest likelihood of being situated at a distance of 0.23 nm. To determine the precise quantity in the MD process, we quantified the number of CO₃²⁻ ions in close proximity to Ca²⁺ ions, utilizing a distance threshold of 0.25 nm, marginally exceeding the 0.23 nm distance observed in Figure 3.

According to Figure 4, in the system without scale inhibitor, Ca²⁺ ions and CO₃²⁻ ions exhibit rapid aggregation as soon as the simulation begins, with an average of 6.5 CO₃²⁻ ions in close proximity to each Ca²⁺ ion when reaching equilibrium. This observation aligns with the characteristic structure of CaCO₃ crystals.³³ Conversely, in the system with scale inhibitor EDTA, the proximity of CO₃²⁻ ions to Ca²⁺ ions is significantly reduced, with only about 3 CO₃²⁻ ions in close proximity to each Ca²⁺ ion when reaching equilibrium, suggesting that the inhibitor effectively hinders the formation of CaCO₃ crystals.

Figure 4.

Number of Ca²⁺ near 0.25 nm of CO₃²⁻.

In the course of MD simulations, particles within the system exhibit continuous movement from their initial positions, and their positions are different at each moment. The quantification of this displacement is commonly referred to as the mean square displacement (MSD). The MSD is able to assess the extent to which the position of a particle diverges from its initial state over a defined duration. In this study, we used MSD to quantify the movement of Ca²⁺ and CO₃²⁻. The degree of deviation of the positions of these molecules from their initial reference positions is determined over time through the application of equation (1), where $N$ is the number of target molecules and $r_{i} (t)$ is the position of molecule i at time $t$

M S D = \frac{1}{N} \sum_{i = 1}^{N} {[r_{i} (t) - r_{i} (0)]}^{2}

(2)

Figure 5 demonstrates that the MSD curves of Ca²⁺ and CO₃²⁻ ions without scale inhibitor exhibit close proximity, suggesting a consistent pairing of these ions; this is caused by the attraction of charges, and the diffusion coefficient is very close, about 0.3769 cm² s⁻¹. On the other hand, the presence of scale inhibitors in the solution results in a notable decrease in the diffusion coefficient of Ca²⁺ and CO₃²⁻ ions, and it is about 0.1745 cm² s⁻¹. This reduction in diffusion coefficient indicates that EDTA impedes the mobility of these ions, thereby hindering the formation and growth of CaCO₃ crystals.

Figure 5.

MSD of Ca²⁺ and CO₃²⁻ in the two systems.

The weak interactions between molecules/atoms are also a very worthwhile aspect to study. In 2022, Lu introduced a simple and intuitive method, known as IGMH,²⁸ to study the weak interactions between molecules/atoms. Following the simulation of equilibrium, our research primarily focuses on exploring the interactions between ions. Of particular interest are the interactions between Ca²⁺ and CO₃²⁻ ions, as well as those between EDTA and Ca²⁺ ions.

As presented in Figure 6, there are two kinds of weak interactions in the IGMH analysis graphs, appearing as blue and green isosurfaces. There are several very clear blue isosurface between Ca²⁺ and CO₃²⁻ ions; unquestionably, they are attributed to the attraction of opposite charges. The green isosurfaces between CO₃²⁻ indicate the existence of van der Waals interactions. It should be noted that a crystal-like configuration is formed between Ca²⁺ and CO₃²⁻; please refer to the following URL for observation: https://youtu.be/IgpQ9f_bAIM. On the other hand, the isosurface between EDTA and Ca²⁺ ions is also blue, indicating that it is also formed by opposite charge attraction. It should be noted that both carboxylate groups in EDTA bind to the same Ca²⁺.

Figure 6.

IGMH analysis of weak intermolecular interactions.

From the colors of the isosurface in Figure 6, it seems that the attraction between CO₃²⁻ and Ca²⁺ is stronger than that between EDTA and Ca²⁺. This contradicts the theory of EDTA as a scale inhibitor. Quantitative assessment of interaction strength serves as a key indicator for evaluating binding stability between molecules or ions. The method of energy decomposition is used to provide a quantitative characterization of the various sources of intermolecular interactions. To investigate the interaction energy between ions, the PSI4 code³⁴ was used for high-order symmetry-adapted perturbation theory (SAPT) analyses. Because Ca belongs to the fourth period element, and due to our limited computing resources, we cannot use extremely advanced base sets aug-cc-pVTZ, the def2-TZVPP basis set was used, with the keywords SAPT2+(3) δ MP2/def2-TZVPP. The SAPT input files of PSI4 were generated with the help of Multiwfn program.³²

From the data presented in Table 2, it is evident that the binding energy between ions is negative, suggesting that the attraction between these species occurs spontaneously without the need for an external driving force. Further analysis through energy decomposition reveals that electrostatic interactions predominantly contribute to this phenomenon, underscoring the spontaneous nature of the mutual attraction between anions and cations. It can be inferred that the interaction energy between EDTA and Ca²⁺ ions is greater than that between CO₃²⁻ and Ca²⁺ ions, indicating a higher propensity for Ca²⁺ ions to bind with EDTA and thereby prevent the precipitation of CaCO₃. This fundamental principle of utilizing EDTA as a scale inhibitor is also applicable in the research and development of novel scale inhibitors.

Table 2.

Analysis of interaction energy between ions (kcal mol⁻¹).

System	Total energy	Electrostatics	Exchange	Induction	Dispersion
Ca²⁺···CO₃²⁻	–516.82	–502.12	58.17	–68.27	–4.60
Ca²⁺···EDTA	–718.04	–665.89	42.71	–89.76	–5.09

Conclusion

Based on MD simulations, this study investigates how EDTA works as a scale inhibitor. Through the application of the IGMH method, this study illustrated the weak interactions between CO₃²⁻ ions, EDTA, and Ca²⁺ ions, revealing that mutual attraction is predominantly driven by charge interactions. In the presence of scale inhibitors EDTA, a significant reduction in the number of CO₃²⁻ ions in close proximity to Ca²⁺ ions was observed. The energy decomposition method was used to quantify ion interaction strength. The findings indicated that the interaction energy between EDTA and Ca²⁺ ions approached −718 kcal mol⁻¹, whereas the interaction energy between CO₃²⁻ and Ca²⁺ ions was approximately −516 kcal mol⁻¹. Evidently, EDTA exhibits a greater affinity for Ca²⁺ ions, impeding the association between Ca²⁺ ions and CO₃²⁻ ions, thereby inhibiting the formation of CaCO₃ crystals.

Footnotes

Authors’ Note

X.J. and B.Z. wrote the main manuscript text, and L.Z. and F.F. prepared Figures 1 – and Y.L. polished the language. All authors reviewed the manuscript.

Availability of data and materials

All the data that support the findings of this study are available on request from the corresponding author.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by Sichuan Provincial Science and Technology Plan Project (2023YFG0092).

Consent for publication

The authors declare that they have no known competing financial interests or personal relationships that may affect the work reported in this paper.

ORCID iD

Xianwu Jing

References

Wei

Jinli

Junwei

, et al. Origin of upper triassic formation water in middle Sichuan Basin and its natural gas significance. Petrol Explor Develop 2009; 36(4): 428–435.

Kamal

Hussein

Mahmoud

, et al. Oilfield scale formation and chemical removal: a review. J Petrol Sci Eng 2018; 171: 127–139.

Al Rawahi

Shaik

Rao

. Studies on scale deposition in oil industries & their control. Int J Innovat Res Sci Technol 2017; 3(12): 152–167.

Loureiro

Martins

Gonçalves

, et al. Large-scale pipe flow experiments for the evaluation of nonchemical solutions for calcium carbonate scaling inhibition and control. SPE J 2023; 28(1): 201–214.

Wang

Chen

. Inhibition of calcium carbonate scale under severe conditions. Petrol Sci Technol 2022; 42: 2234–2249.

Cirimello

Otegui

Carfi

, et al. Failure and integrity analysis of casings used for oil well drilling. Eng Fail Anal 2017; 75: 1–14.

Jafar Mazumder

. A review of green scale inhibitors: process, types, mechanism and properties. Coatings 2020; 10(10): 928.

Mpelwa

Tang

. State of the art of synthetic threshold scale inhibitors for mineral scaling in the petroleum industry: a review. Petrol Sci 2019; 16: 830–849.

Housse

Hadfi

Karmal

, et al. Toxicity profile, phytochemical composition, and anti-scaling properties of the aqueous extract of Ocimum basilicum L. leaves as novel green and cost-effective inhibitor: experimental, MC/SAA and DFT approach. Waste Biomass Valoriz 2023; 14(11): 3553–3573.

10.

Chen

Liang

. Synthesis of polyaspartic acid/graphene oxide grafted copolymer and evaluation of scale inhibition and dispersion performance. Diamond Rel Mater 2020; 108: 107949.

11.

Zuo

Yang

Zhang

, et al. Effect of scale inhibitors on the structure and morphology of CaCO₃ crystal electrochemically deposited on TA1 alloy. J Colloid Interf Sci 2020; 562: 558–566.

12.

Zeng

Wang

Gong

. Molecular dynamics simulation of the interaction between polyaspartic acid and calcium carbonate. Mol Simul 2013; 39(3): 169–175.

13.

Zhang

Zhao

, et al. Inhibition of calcium carbonate nucleation and crystallization by carboxymethyl dextran: experiments and molecular dynamics simulations. Desalination 2024; 573: 117205.

14.

Gopi

Subramanian

. Polymorphism in CaCO₃—effect of temperature under the influence of EDTA (di sodium salt). Desalination 2012; 297: 38–47.

15.

Nirmala Devi

Sanjiv Raj

Subramanian

. Synergistic effects of magnesium and EDTA on polymorphism and morphology of CaCO₃ and its influence on scale. J Cryst Growth 2021; 564: 126108.

16.

Abraham

Murtola

Schulz

, et al. GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX 2015; 1–2: 19–25.

17.

Van Der Spoel

Lindahl

Hess

, et al. GROMACS: fast, flexible, and free. J Comput Chem 2005; 26(16): 1701–1718.

18.

Wang

Wolf

Caldwell

, et al. Development and testing of a general amber force field. J Comput Chem 2004; 25(9): 1157–1174.

19.

Schauperl

Nerenberg

Jang

, et al. Non-bonded force field model with advanced restrained electrostatic potential charges (RESP2). Commun Chem 2020; 3(1): 44.

20.

Sengupta

Song

, et al. Parameterization of monovalent ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB water models. J Chem Inform Model 2021; 61(2): 869–880.

21.

Song

L F

, et al. Systematic parametrization of divalent metal ions for the OPC3, OPC, TIP3P-FB, and TIP4P-FB water models. J Chem Theor Comput 2020; 16(7): 4429–4442.

22.

Izadi

Onufriev

. Accuracy limit of rigid 3-point water models. J Chem Phys 2016; 145(7): 074501.

23.

Martínez

Andrade

Birgin

, et al. PACKMOL: a package for building initial configurations for molecular dynamics simulations. J Comput Chem 2009; 30(13): 2157–2164.

24.

Carretero-González

Kevrekidis

, et al. A Parrinello–Rahman approach to vortex lattices. Phys Lett A 2005; 341(1–4): 128–134.

25.

Bussi

Donadio

Parrinello

. Canonical sampling through velocity rescaling. J Chem Phys 2007; 126(1): 14101.

26.

Darden

York

Pedersen

. Particle mesh Ewald: an N log (N) method for Ewald sums in large systems. J Chem Phys 1993; 98(12): 10089–10092.

27.

Humphrey

Dalke

Schulten

. VMD: visual molecular dynamics. J Mol Graph 1996; 14(1): 33–38.

28.

Chen

. Independent gradient model based on Hirshfeld partition: a new method for visual study of interactions in chemical systems. J Comput Chem 2022; 43(8): 539–555.

29.

Neese

. The ORCA program system. Wiley Interdiscip Rev: Comput Mol Sci 2012; 2(1): 73–78.

30.

Neese

Wennmohs

Becker

, et al. The ORCA quantum chemistry program package. J Chem Phys 2020; 152(22): 224108.

31.

Neese

. Software update: the ORCA program system—version 5.0. WIREs Comput Mol Sci 2022; 12(5): e1606.

32.

Chen

. Multiwfn: as multifunctional wavefunction analyzer. J Comput Chem 2012; 33(5): 580–592.

33.

De Villiers

. Crystal structures of aragonite, strontianite, and witherite. Am Mineral: J Earth Planet Mater 1971; 56(5–6): 758–767.

34.

Parrish

Burns

Smith

, et al. Psi 4 1.1: an open-source electronic structure program emphasizing automation, advanced libraries, and interoperability. J Chem Theor Comput 2017; 13(7): 3185–3197.