Molecular dynamics simulation of ammonium ion removal by freezing concentration

Ammonium wastewater is a serious and common water pollutant that can have harmful effects on the environment. Freeze concentration, as an energy-efficient and environmentally friendly method, is used to treat ammonium wastewater by ice-water phase transition. The simulation results show that most of the ions are retained in the liquidphase, and it is reported for the first time that the probability of NH4 + (90%) remaining in the water is significantly higher than that of Cl− (67%). We have analyzed the influence of ions on ice/water structure from the perspective of structure and energy and explained the reason for the difference in the probability of NH4 + and Cl− remaining in the liquid phase. We find that the coordination number (CN) of NH4 + decreases from 6 to 4 when one NH4 + permeates the ice layer, indicating that the first hydration layer of ammonium ions underwent significant reorganization during this period. In contrast, a similar reduction in CN was not observed during the entry of Cl− into the ice layer. Moreover, the hydration energy shows that NH4 + prefers to stay in the liquid phase than in the ice phase because of the higher hydration energy difference compared with that of Cl−. The results of this work indicate that freeze concentration can efficiently remove NH4 + by ice-water phase transition, which greatly reduces the discharge of ammonium wastewater and pave the way for further study of the freezing process for wastewater treatment.


1. Introduction
Excessive release of nitrogen can have harmful consequences for the water environment [1], such as eutrophication [2,3], and toxic algal blooms [4,5]. Therefore, nitrogen-containing wastewater must be effectively treated before discharge. To remove nitrogen from wastewater, a lot of work has been done, such as ammonia stripping [6], breakpoint chlorination [7], ion exchange [8], electrodialysis [9], and biological nitrification [10], etc. The ammonia stripping method changes free ammonia from liquid phase to gas phase by steam stripping in high pH value, which is suitable for high-concentration ammonia nitrogen wastewater but consumes a lot of heat energy [11]. Meanwhile, the treatment efficiency declines as the ion concentration decreases [11]. At this point, chemical and biological methods are generally used for low-concentration wastewater below 1000 mg l −1 . For example, the breakpoint chlorination method converts ammonia nitrogen into nitrogen by adding chlorine and the ion exchange method uses the exchanger and adsorbent made of nanomaterials to remove ammonia nitrogen [12][13][14][15]. However, these methods either cause secondary pollution or incur additional costs. Biological nitrification oxidizes ammonia nitrogen by bacteria, but this microbial activity is susceptible to environmental influences and the reaction speed is slow [16]. In addition to the above methods, freeze concentration, as an environment-friendly, cost-effective, and efficient method, has been developed recently [17][18][19][20][21][22].
Freeze concentration refers to the process where the dissolved solutes are expelled from the ice phase during freezing [23]. As the impurities are rejected ahead of the ice crystals grow under freezing conditions, freeze concentration can be used for the treatment of wastewater within a relatively large concentration range [22]. Because the latent heat of freezing (334 kJ kg −1 ) is significantly smaller than that of evaporation (2257 kJ kg −1 ) [24], freeze concentration is an energy-efficient and environment-friendly process with respect to traditional methods [17-22, 25, 26]. Especially, there is a long and cold winter in most of northern China, which provides suitable operating conditions for this method. As a result, freeze concentration can be applied to the primary treatment of most wastewater and has led to extensive research in recent years. Wang R et al used the natural freezing method to remove ammonia nitrogen, and the removal rate reach about 80% [26]. Larsen T et al investigated the removal effect of the freeze-thaw on ammonia nitrogen and reported removal rates ranging from 80% to 90% [25]. Furthermore, Okada T et al found an imbalance in the distribution of anions and cations at the ice interface with the liquid phase in the frozen electrolyte by fluorometric measurements of pH [27], but the reason for the imbalance of ion distribution is unclear and requires observation at the molecular mechanism level, which is experimentally challenging. As an alternative, molecular dynamics (MD) simulations have been used to probe the molecular mechanism of ions removal during freezing.
MD simulations have made a lot of progress in understanding the motion of ions during freeze concentration [28][29][30][31][32]. Taking the crystallization of sodium chloride solution as a model, Gallo P et al showed that the presence of ions reduces the chemical potential of water, thereby lowering the freezing point of water [29]. Furthermore, Gillet P et al found that ions significantly affect the structure of the hydrogen bond network of water and consequently the ice crystallization process [31]. Gallo P et al indicated that one chloride ion can replace two oxygen sites in the ice crystal, while one sodium ion occupies the void space in the ice crystal, and the addition of sodium ion can cause local deformation of the ice structure [30]. Most of the previous MD simulation studies have focused on sodium and chloride ions, while ammonium ions have been less explored and require extensive research.
In this work, we perform molecular dynamics simulations to explore the possibility of removing NH 4 + by modeling the crystallization process of the NH 4 Cl solution. We find that NH 4 + and Cl − can be separated upon freezing and calculate the removal rate of NH 4 + for the first time. Then we explain the removal rate difference between NH 4 + and Cl − from structural and energetic aspects. The structural reason is presented by analyzing the radial distribution function, the distribution of angles in the first hydration shell, and the coordination number (CN) around the ions. When NH 4 + enters the ice phase from the liquid phase, the hydration shell changes from an octahedral-like structure to a tetrahedral-like structure, and the CN of NH 4 + decreases from 6 to 4. While for Cl − , the CN peak of Cl − is around 6 in both the ice and liquid phase, indicating a significant reorganization of the first hydration layer of NH 4 + . The hydration energy calculation further confirms that NH 4 + prefers to stay in the liquid phase because the ionic hydration energy of NH 4 + in ice crystals is higher than that of NH 4 + in aqueous solution. Furthermore, we compute the free energy barrier of both ions escaping from the ice to the liquid phase and demonstrate that it is easier for NH 4 + to overcome the energy barrier to return to the liquid phase compared with Cl − . These simulation results provide a feasible theoretical basis at the molecular level for removing NH 4 + from the NH 4 Cl solution through freeze concentration.

Simulation system
The simulated box consists of an NH 4 Cl solution and two layers of hexagonal ice, serving as a prismatic plane for triggering ice growth (as shown in figure 1). The oxygen atoms of the initial ice molecules are harmonically confined to their original positions with a force constant of 1000 kJ mol −1 nm −2 . The NH 4 Cl solution includes 8340 water molecules and 15 pairs of NH 4 + and Cl − , which make the concentration equal to 0.1 mol l −1 . At the initial time, the ions are randomly distributed in the solution.
To avoid unnecessary interactions between the ice slab and the solution at the boundary of the simulation box, a 2.5 nm vacuum region was created on the right side of the solution. Then, the dimensions of the whole system box are 6.35 × 5.17 × 11.0 nm 3 . Periodic boundary conditions are applied in three directions.

Molecular dynamics simulation
We performed all MD simulations using the GROMACS package [33] to investigate the removal of ions during ice crystal growth. The OPLS-AA force field [34] and TIP4P model [35] are used to describe the interactions between ions and water, respectively. The intermolecular van der Waals and Lennard-Jones interactions are truncated at 1 nm and the long-range electrostatic interactions are evaluated with the particle-mesh Ewald algorithm [36]. The leapfrog algorithm with a time step of 2 fs is applied to solve the equations of motion. The Nosé-Hoover thermostat is employed to control the temperature of the system. After the steepest descent energy minimization, the system was first equilibrated through a 100 ps NVT at 200 K. Then, a 600 ns NVT simulation was performed at 200 K as the production run.

Removal rate
The removal rate of ions during the freezing progress is: where C ice and C water represent the ion concentration in the ice and liquid phase, respectively. The ice and water molecules are identified by the CHILL + algorithm [37]. The ice/waterinterface and the water/vapor interface are defined as the 50 wt% of water along the z direction (red dashed line in figure 1).

Free energy
To understand the ion removal mechanism, the free energy of an ion moving from the ice/water interface into bulk water is computed as the potential of mean force (PMF). The simulation temperature is set to 237 K, 5 K above the freezing point of the TIP4P water model, to prevent the growth of ice. A series of NVT simulations for 65 umbrella windows are performed with harmonic-restraint potentials, ( ) where the force constant k = 1000 kJ mol −1 nm −2 , the equilibrium distance d 0 spans from 0.2 to 1.5 nm every 0.02 nm, and z is the absolute position of the center of mass of the solute in the direction perpendicular to the ice surface. We ran 20 ns for each window and collected the last 15 ns for the free energy calculation. The resulting z distributions are combined using the weighted histogram analysis method [38]. The PMFs as a function of the absolute z value are later replotted against the distance from the ice surface, where zero is defined at the averaged coordinate of oxygen of the restrained ice molecules. The average PMF at distances from 1.2 to 1.5 nm is chosen as the baseline (i.e., the zero of the free energy profiles).

3. Results and discussion
The simulation system and the crystallization progress of the NH 4 Cl solution are shown in figure 1 and Movie S1. With the continuous growth of ice and the reduction of water molecules in the ionic solution, the NH 4 + and Cl − are carried by the water molecules and slowly enter the ice layer, which is consistent with previous MD simulations [29]. Interestingly, there are more Cl − ions trapped in the ice phase than NH 4 + , which may be related to the hydration shell of the two ions. The ice/waterinterface and the water/vapor interface (two red dash lines in figure 1) are defined according to the 50 wt% of water along the z direction.
We quantitatively analyze the process of ion separation by calculating the ion removal rate. As shown in figure 2, The removal rate of NH 4 + is significantly higher than that of Cl − . During the first 200 ns of the simulation, only the ice-water phase transition occurs in the system, and the ions do not enter the ice phase, which means that both NH 4 + and Cl − are staying in the liquid phase. As the ion concentration increases due to the growth of ice, the NH 4 + and Cl − start to enter the ice phase after 200 ns. Once ions enter the newly grown ice phase, the ion concentration of the ice phase increases, and the ion removal rate decreases. Figure 2 shows that the removal rate of NH 4 + is always 10%-30% higher than that of Cl − during 200 ∼600 ns, which means that NH 4 + has a greater possibility to stay in the liquid phase than Cl − during the freeze concentration process. We find that in 0.23 mol l −1 NH 4 + solution, the removal rate of NH 4 + could be maintained at about 90% even at 600 ns. This indicates that freeze concentration can effectively treat ammonia wastewater.
Here, we investigate the reasons for the difference in separation efficiency of NH 4 + and -Cl from the structural point of view. We calculate the ion-oxygen radial distribution function, g ion-o (r), to analyze the structural properties of the hydration shells in both the ice and liquid phase ( figure 3). For the ions in the liquid phase, the NH 4 + -O correlation contains the first peak at r = 2.70 Å and for -Cl -O at r = 3.20 Å. This difference is due to the influence of the positive charge of the NH 4 + which attracts the oxygen, whereas in the case of -Cl , the negative charge attracts the hydrogen and results in larger -Cl -O distances. In addition, the positions of the two ions in the second peak are similar, which indicates that the influence of the two ions on the water structure is mainly within the first hydration shell. For the ions in the ice phase, the two ions have the same first peak position as that in the liquid phase, but both peaks are higher than that in the liquid phase. This is consistent with the fact that water molecules in the ice phase are more ordered compared to that in the liquid phase. Meanwhile, the position difference of the second peak of the two ions in the ice phase is larger than that in the liquid phase, which indicates that ions have a longer-range influence on the structure of the ice phase with respect to that of the water phase.  Since the influence of the two ions on the ice and water structure is mainly within the first hydration shell as defined by the first minimum of the radial distance of g ion-o (r), we further analyze the spatial geometry of the first hydration shell by comparing the O-ion-O angular distribution. As shown in figure 4, the angular distributions of both NH 4 + and Cl − in the liquid phase are broader than that in the ice phase, which is consistent with the fact that the distribution of water molecules around the ion in the liquid phase is disordered. For the ions in the ice phase, the main peak of the angular distribution of NH 4 + locates at about 105°, and the water molecules within the first hydration shell exhibit a tetrahedron-like structure ( figure 4). However, the main peak of the Cl-angular distribution locates at about 90°and 170°where the water molecules surrounding the ion exhibit a twisted octahedron-like structure ( figure 4). The reorientation of the first hydration shell of the ions indicates that the surrounding water molecules undergo significant structural changes during the phase transition.
In addition to the angular distribution analysis, the aforementioned tetrahedron-like and twisted octahedron-like hydration shells can be further characterized by the CN of the ion. As shown in figure 5, the CN distributions of NH 4 + (5.54) and Cl − (5.73) in the liquid phase are similar. Their CN distribution ranges from 3 to 8, with the highest peak at CN = 6. While for the ions in the ice phase, the peak of CN distribution of NH 4 + and Cl − moves to 4 and 6 respectively. The corresponding number of CN peaks in the ice phase is consistent with the tetrahedron-like and twisted octahedron-like hydration shells of NH 4 + and Cl − . Especially, the overall CN of Cl − is obviously larger than that of NH 4 + in the ice phase, with an average CN of 5.48 for Cl − and 4.19 for NH 4 + ,  respectively (see table 1). This indicates that Cl − attracts more ice molecules than NH 4 + within the first hydration shell and suggests a higher interaction between ice and Cl − compared with that between ice and NH 4 + .
Changes in water and ice structure lead to the variation of the ionic hydration energy. To account for the changes in the structure of water molecules around the ions, we compute the hydration energy of NH 4 + /Cl − in liquid and ice phases during the freezing process. The hydration energy of each ion is defined as the average potential energy of ions and water/ice molecules in the first hydration shell, as shown in table 2. The hydration energy of NH 4 + in the liquid phase is 17.07 kcal mol −1 lower than that in the ice phase, which indicates that NH 4 + tends to stay in solution as the ice grows. The hydration energy difference ΔE is 36.27 times that of the typical thermal energy K B T, which indicates that the hydration energy difference is not caused by the thermal motion of NH 4 + ions, but rather by the repulsion of NH 4 + by the ice crystals. However, the hydration energy of Cl − in the liquid phase is 0.62 kcal mol −1 higher than that of the ice phase, suggesting that as the ice grows, the Cl − ions prefer to stay in the ice. The hydration energy difference ΔE is 1.32 times the typical thermal energy K B T, which is much smaller than the hydration energy difference ΔE of NH 4 + . The hydration energy difference suggests that a large hydration energy barrier needs to be overcome when NH 4 + enters the ice phase, but not for Cl − .
Therefore, NH 4 + is more likely to stay in the liquid phase than Cl − .
Once trapped in the ice phase, the ability of ions to escape from the ice phase and return to the liquid phase also affects the separation efficiency of the freeze concentration. In order to investigate which ion is more likely to escape the ice phase, we compute the free energies of NH 4 + and Cl − , as shown in figure 6. It can be seen that the free energy distribution exhibits multiple local minima. The free energy of Cl − to ice/waterinterface has two minima with comparable stability, corresponding to Cl − in contact with ice and separated from ice by a monolayer of water. While the free energy of NH 4 + on the ice surface has two shallower minima, indicating NH 4 + is unstable on the ice/waterinterface. On the other hand, the free energy variation generates an energy barrier ΔG, which is the PMF difference between a local minimum and an adjacent maximum. ΔG of ions reflects the    Table 2. Average hydration energy E per ion for the interactions of NH 4 + and Cl − with water and ice. ΔE denotes the energy difference between ion-water and ion-ice interactions. probability of an ion during the diffusion from the ice to the liquid phase. The higher energy barrier for Cl − than NH 4 + means that it is easier for NH 4 + to escape from the ice to the liquid phase compared with Cl − .

4. Conclusions
To conclude, we have performed MD simulations to investigate the removal efficiency of NH 4 + and Cl − during the freezing of NH 4 Cl solutions. The simulation results show that the removal rate of NH 4 + (90%) is significantly higher than that of Cl − (65%) during the ice-water phase transition, which strongly proves that freeze concentration is another new solution to treat ammonium wastewater besides reverse osmosis treatment and biological treatment. The structural and energy analysis of the first hydration shell of NH 4 + and Cl − revealed that NH 4 + has a greater resistance to entering the ice layer than Cl − due to the reorganization of surrounding water molecules and the increased hydration energy. On the other hand, the lower free energy barrier also makes it easier for NH 4 + to return to the liquid phase compared with Cl − . These simulation results of this work provide insights into the microscopic mechanism of the removal of NH 4 + and Cl − , which may facilitate the development of preliminary treatment of ammonium wastewater by freeze concentration.