Interfacial and bulk properties of concentrated solutions of ammonium nitrate †

We conducted molecular dynamics (MD) simulations to calculate the density and surface tension of concentrated ammonium nitrate (AN) solutions up to the solubility limit of ammonium nitrate in water, by combining the SPC/E, SPCE/F and TIP4P/2005 water models with OPLS model for ammonium and nitrate ions. This is the first time that the properties of concentrated solutions of nitrates, especially AN, have been studied by molecular dynamics. We eﬀectively account for the polarisation eﬀects by the electronic continuum correction (ECC), practically realised via rescaling of the ionic charges. We found that, the full-charge force field MD simulations overestimate the experimental results, as the ions experience repulsion from the interface and prefer to remain in the subsurface layer and the bulk solution. In contrast, reducing the ionic charges results in the behaviour that fits well with the experimental data. The nitrate anions display a greater propensity for the interface than the ammonium cations. We accurately predict both the density and the rise in the surface tension of concentrated solutions of AN, recommending TIP4P/2005 for water and the scaled-charge OPLS model (OPLS/ECC) for the ions in the solutions. We observe that, the adsorption of anions to the interface accompanies their depletion in the subsurface layer, which is preferentially occupied by cations, resulting in the formation of the electric double layer. We demonstrate the ion deficiency for up to 3 Å below the surface and establish the requirement to include the polarisability eﬀects in the OPLS model for AN. While these results confirmed the findings of the previous studies for dilute solutions, they are new in the solubility limit. Concentrated solutions exhibit a strong eﬀect of the abundance of solute on the coordination numbers of ions and on the degree of ion pairing. Surprisingly, ion pairing decreases significantly at the interface compared with the bulk. The present study identifies OPLS/ECC, along with TIP4P/2005, to yield accurate predictions of physical properties of concentrated AN, with precision required for industrial applications, such as a formulation of emulsion and fuel-oil explosives that now predominate the civilian use of AN. An application of this model will allow one to predict the surface properties of supersaturated solutions of AN which fall outside the capability of the present laboratory experiments but are important industrially.


Introduction
Ammonium nitrate (NH 4 NO 3 , AN) finds applications as an oxidiser in rocket propellant engines, gas generator systems (e.g., automotive inflator systems, automobile airbag systems and heavy-lift launchers) and pyrotechniques, due to its low cost, chemical stability, moderate flame temperature and halogen-free exhaust gases. The chemical serves as the main ingredient of fertilisers, herbicides, insecticides, AN emulsion and fuel-oil (ANFO) explosives, the most frequently employed commercial explosives. The physical stability and detonation of AN emulsions and ANFO highly depend on emulsion formulation and the wettability of AN prills, with the latter affected by several factors including surface tension, density, viscosity and polarisability of the solvent. 1 Therefore, determination of density and surface tension of aqueous solutions of AN, especially concentrated solutions, constitutes the essential requirement for engineering applications of this material. Besides the thermodynamic models suggested by several investigators, 2-5 molecular simulations have contributed significantly to the present understanding of the surface tension of ionic solutions at the microscopic level, in particular, by describing the distribution of ions near liquid-vapour interfaces. [6][7][8][9][10][11][12] Normally, the presence of ions at the interface, which attract each other through electrostatic forces, and interact with the solvent through the solute-solvent attraction, increases the surface tension, as there is more resistance to stretch the surface. This makes the formation of new interfaces more energy intensive during the emulsification. To the best of our knowledge, solutions of ammonium nitrate have been largely neglected as a research subject of molecular simulations, despite the practical importance of concentrated solutions of AN. Moreover, a study of the surface activity of germane ions offers insights into heterogeneous aerosol activities and their potential impact on the atmospheric and tropospheric chemistry. 13 Researchers studied the effect of size, charge and the polarisability of the ions on their location in the solution. [14][15][16][17][18][19][20] Vrbka et al. 21 reported that large polarisable ions such as heavier halides exhibit an affinity for water-vapour interface while small non-polarisable ions (e.g. alkali cations and fluoride) repel from the interface. Based on the results from the polarisable potential models, Dang 16 reported that anions display increased affinity, in comparison with cations, for adsorbing to the interface. The larger the anions, the stronger binding to the liquid-vapour interface. Jungwirth and Tobias 17 emphasised the effect of size and polarisability of ions at the interface and in the interior of a water slab, and showed that, for all concentrations of aqueous solutions of sodium chloride, the chloride anions occupy a significant portion of slab interfaces, while sodium cations locate themselves toward the interior of the slab.
A controversy exists in the literature about the inclination of nitrate anions at the liquid-vapour interface. Some researchers 22,23 suggested no affinity of nitrate cations for the interface, whereas others [24][25][26] argued for their weak propensity to partition to the interface. Tian et al. 27 studied the liquidvapour interfaces of several salt solutions using phase-sensitive sum-frequency vibrational spectroscopy. They revealed the presence of an electric double layer formed by cations and anions induced by the relative attraction to the surface of different ions, in decreasing order of I À , NO 3 À , NH 4 + , Cl À , K + , Na + , and SO 4 2À . Hua et al. 28 also studied the distributions of interfacial ions of some nitrate salt solutions (LiNO 3 , NaNO 3 , NH 4 NO 3 , and Mg(NO 3 ) 2 ) using the same experimental technique. They confirmed the appearance of an electric double-layer structure with the greater abundance of NO 3 À ions at the surface as compared with their counterions. From their measurements, the induced net electric field in AN solution was smaller compared to the other nitrate solutions which suggested the least charge separation distributions between the NH 4 + and NO 3 À among the studied cations. Therefore, NH 4 + tended to position themselves closer to NO 3 À .
It is well documented that, in systems of aqueous electrolyte solutions, the electronic polarisation effects play a critical role. For instance, molecular simulations, using non-polarisable force fields, introduce inaccuracies in describing the shortrange ion-water and ion-ion interactions, resulting in artificial ion clustering. [29][30][31] Although polarisable models describe more truthfully the interactions in these systems, they suffer from expensive computations and difficult parametrisation. For these reasons, previous researchers attempted to incorporate polarisation effects into non-polarisable force fields. Recently, Leontyev et al. [32][33][34] suggested an electronic continuum model that incorporates the electronic polarisation into a non-polarisable force field. In this method, all ion charges are scaled with the inverse square-root of the high-frequency dielectric permittivity of water (i.e., with (e el ) À1/2 E 0.75).
The ab initio molecular dynamics (AIMD) simulations provide an alternative approach to account for the electronic polarisation effects. 35 In an AIMD calculation, Kohn-Sham density functional theory (DFT) framework serves to compute the electronic structures of the ground states of the molecules. The polarisation of electron clouds can result from the response of the electron densities to the surrounding electric field. Forces obtained directly from the calculations of electronic structures, performed ''on the fly'', propagate the dynamics of the nuclei. AIMD simulations have been extensively applied to gain understanding of molecular environments in water 36 and in aqueous solutions. [37][38][39] A similar ''charge scaling'' commonly arises in simulations of ionic liquids using non-polarisable force fields. [40][41][42] We refer to this methodology throughout this study as the electronic continuum correction (ECC) model. Applying the ECC model for simulations of liquid water interacting with proteins and solutes in the bulk yielded the same results as those obtained from models that explicitly introduce the electronic polarisation. 32,43,44 The importance of the ECC model manifests itself for both monovalent and multivalent ions, as they are capable of polarising their surrounding water molecules. 30,[45][46][47][48][49] Mason et al. 46 and Pegado et al. 50 showed that, predictions of the properties of aqueous salt solutions in the bulk improved significantly using the electronic continuum correction model compared with non-polarisable simulations. However, the applicability of the ECC model has not been fully examined for liquid-liquid and liquid-vapour interfaces.
Vazdar et al. 51 investigated the affinity of halide ions to the water-oil and water-vapour interfaces using the Smith-Dang/ ECC model to account for the polarisation effects. They reported interfacial affinities of anions to the water-oil interfaces consistent with the results of experiments and previous explicitly polarisable models. However, their model overestimated the anionic surface affinities for water-vapour interfaces due to electronic discontinuity upon moving from liquid to the vapour phase.
Neyt et al. 18 reported a non-monotonic increase of the surface tension with rising concentration of sodium chloride using the OPLS/ECC and Reif/ECC models, where OPLS stands for optimised potential for liquid simulations. They emphasised a requirement for the consistent re-parametrisation of the non-polarisable force fields for the salt and the scaling of the charges. Contrary to Neyt's et al. unsuccessful application of the ECC model, Otten et al. 52 showed an enhanced surface affinity of iodide to a water-vapour interface once they reduced iodide charge, similarly to the prediction from the ECC model, as compared with the results of non-polarisable simulations using full charges on the ions. Breton and Joly, 53 in their paper on solutions of sodium chloride, have recently reported the importance of using long-range solvers for Lennard-Jones interactions in non-uniform systems, such as those featuring an interface, to predict the dispersion forces and applying the ionic charge rescaling factor to demonstrate good agreement between computed and measured surface tension and interfacial structure.
We were intrigued by the reports in the literature on the effect of rescaling of ionic charges on the ability of nonpolarisable force field to yield improved values of macroscopic properties, and were curious about the applicability of this approach to predict the bulk and surface properties of a wide range of solutions of ammonium nitrate, starting with dilute solutions up to highly concentrated systems; the latter close to the solubility limit of ammonium nitrate in water (65 wt% equivalent to 23 mol kg À1 at 298. 15 K 54 ). For this reason, we computed the bulk density, structure of the solvation shells and the surface tension, with and without rescaling of the charges. In the sections that follow, we complement these results with the discussion that compares the distribution and the number of molecules in the first solvation shell of the ions in the bulk and at the liquid-vapour interfaces to discover excellent agreement between our modelling results and the experimental measurements of Tian et al. 27 and Hua et al. 28 The same physical behaviours appear both in the model and in the experiments, such as the formation of the ionic double layer induced by the relative propensity of ions to the surface (NO 3 À 4 NH 4 + ), underpinning the observed accord.
This paper is organised as follows: in Section 2, we provide details of our simulations and the analytical methods. Sections 3.1 and 3.2 discuss the results for the bulk solutions of AN, with Section 3.3 focusing on predicting the surface tension of mixtures of AN with water. This leads us to pondering, in Sections 3.4 and 3.5, the distribution of ions and the hydration structure in the interior and at surface of the solution slabs. We summarise the major findings of this investigation in Section 4. As far as we are aware, this is the first MD investigation of concentrated solutions of nitrates, and the first MD study that examines the molecular properties of ammonium nitrate, a technologically important group of chemicals. In the chapters that follow, we will identify a MD model capable of predicting both bulk and surface properties of concentrated solutions of AN with accuracy required for practical applications, such as predictions of the surface tension in making AN emulsion explosives.

Potential models
In the simulations, we consider three non-polarisable models for water, SPC/E, 55 SPCE/F 56 and TIP4P/2005. 57 The SPC/E, also known as the extended simple point charge model, represents rigid water molecules with three-point masses, whereas the SPCE/F serves as a flexible version of the SPC/E water model and TIP4P/2005 considers a rigid four points; two hydrogens, oxygen and one Lennard-Jones interaction site located at the oxygen atom. Table S1 of ESI † lists molecular parameters of the three models. Both SPCE/F and TIP4P/2005 reproduce very well the surface tension of water at 298.15 K. [58][59][60][61][62][63] In this article, we compare the results from the SPCE/F and TIP4P/2005 models with the simple and commonly used SPC/E model to gain a clear overview on the performance and transferability of ion parameters to the three water models. 64,65 We represent the ammonium and nitrate ions (NH 4 + and NO 3 À ) by a molecular force field compatible with the nonpolarisable optimised potentials for liquid simulations all-atom (OPLS-AA), abbreviated in this contribution as OPLS, except for Table 1. 66 This model comprises harmonic bonds and angles, torsion energy defined by cosine series, and non-bonded interactions given by electrostatic partial charges (Coulomb) and 12-6 Lennard-Jones atomic sites. Eqn (1) presents the functional form of the total potential energy (U) where U INTRA and U INTER denote the intramolecular and intermolecular energy contributions, respectively. The parameter U INTRA vanishes for rigid molecules, however, for non-rigid molecules/ions, it is expressed as where the parameters are the force constant (k), equilibrium bond value (r 0 ), equilibrium angle value (y 0 ) and Fourier coefficients (V), with their units provided in Table 1.
The intermolecular interactions are composed of repulsion, dispersion and electrostatic contributions that are represented by the Lennard-Jones (LJ) and coulombic (ELEC) potentials, respectively: where s and e denote the Lennard-Jones radii and potential well depths, q a the charge of atom a and e 0 the vacuum permittivity that amounts to 8.854187 Â 10 À12 F m À1 . Table 1 assembles the OPLS molecular parameters of NH 4 + and NO 3 À ions. We selected these parameters because of their previous application to compute surface tension of aqueous ionic solutions. 10 To maintain compatibility with the OPLS force field, we apply the geometric mixing rules for the Lennard-Jones coefficients, i.e., s ij = (s ii s jj ) 1/2 and e ij = (e ii e jj ) 1/2 . In the current study, we have considered the effects of electronic polarisation of the ions using non-polarisable force fields for both ammonium and nitrate ions, through the ECC model. 32 In this approach, to account for the electronic polarisability, the fixed charges of the ions are scaled by the inverse square-root of the electronic part of the dielectric constant, (e el ) À1/2 , of the surrounding medium. Since the parameter e el has a value of 1.78 for water, the atomic charges are rescaled by the factor of (1.78) À1/2 . The last column of Table 1 lists the rescaled charges of ions deployed in this study.

Molecular dynamics
The LAMMPS package 70 served to perform the classical molecular dynamic simulations. We approximated the long-range forces due to Lennard-Jones and coulombic interactions with accuracy of 1 Â 10 À6 , with respect to untruncated values, by a particle-particle particle-mesh (pppm) method, implemented by Ismail et al. 63,71 in the package. This approach removes the need for the analytical tail correction in the computation of surface tension. Periodic boundary conditions were applied in the three directions, and the Verlet algorithm integrated the equation of motion with a time step of 1 and 0.5 fs for the rigid and non-rigid water models, respectively. In the rigid models of water, the SHAKE technique 72 constrained the bond lengths and bond angles.
The Packmol software 73 facilitated the random insertion of the ions and water molecules into a box with initial dimensions of L x L y L z (L x = L y = 40 Å). All simulations started from low density configurations with the total number of water molecules fixed to 2133, and the number of NH 4 + and NO 3 À ions increased from 25 to 720 to reflect the target range of molality (0.65 to 18.7 mol kg À1 ). Table 2 summarises the number of molecules and the concentrations of the simulated solutions. Before equilibrating the MD simulations, we applied the conjugate gradient method to minimise the energy and to remove strain from the initial configurations. To calculate the density of the solutions, the systems were subjected to 4 ns equilibration runs for the isobaricisothermal ensemble (NPT) at p = 0.1 MPa and T = 298.15 K, with Nosé-Hoover thermostat and barostat with relaxation time constants, respectively, of 0.1 and 1 ps for the rigid models of water and 0.05 and 0.5 ps for the flexible water. As the length and width of the simulation box remain constant, the changes in the box volume, due to the variation in density at different concentrations, vary the height of the box along the z direction from 41.27 to 75.45 Å. Once the bulk density of each system reached a constant, indicating that the equilibrium has been attained, we computed the bulk density by averaging the values obtained in the last 2 ns. We also identified the bulk solvation structure of ion assemblies in the system. For calculating the surface tension, we again randomly inserted specified number of molecules in a simulation box   (1) 120.00 and then extended the simulation box to three times its initial size along the z direction, to create two interfaces between the AN solution and the vacuum. 74 Fig. 1 depicts a typical molecular configuration of 0.65 mol kg À1 AN solution with 2133 water molecules and 25 pairs of NH 4 + and NO 3 À ions.
We placed two reflecting walls just outside of the interfaces to retain the molecules in the slab and the simulations were started with an NPT ensemble at p = 0.1 MPa and T = 298.15 K to equilibrate the system for 1 ns. Then, we switched to the NVT ensemble at 298.15 K with no reflecting walls with the Nosé-Hoover thermostat relaxation time of 0.1 ps and 0.05 ps for the rigid and flexible water models, respectively. The NVT runs involved 12 Å cut-off distance and 1 and 0.5 fs time steps for the models of rigid and flexible water, respectively. The systems were equilibrated for 1 ns followed by a production run in the range of 15 ns. We finally compared the results for density and surface tension with the available experimental values for solutions of AN.
To validate the computational methodology, we calculated the density and the surface tension of pure (i.e., neat) water within the temperature range of 298.15 to 600 K, and then juxtaposed the results of our calculations against the experimental measurements. Fig. S1 (ESI †) presents the temperature dependence of the surface tension and liquid-vapour coexistence curve of water obtained from the SPC/E, SPCE/F and TIP4P/2005 models and from experimental measurements. 75 [58][59][60][61]63 Table S2 of ESI † summarises all calculated values within a temperature range of 298.15-600 K in addition to the surface thickness in which the density of water increases from 10% to 90% of its bulk value, with this layer denoted as the 10-90 surface thickness (d t ). We described the calculation of d t further in this section.
Moreover, the simulation of pure water from the TIP4P/2005 model, employing a system of 4000 water molecules (L x = L y = 49 Å) or a larger system of 6500 molecules (L x = L y = 58 Å), resulted in essentially no differences in the calculated values of surface tensions, 67.42 AE 0.06 and 67.25 AE 0.07, respectively. Therefore, we expect that, our calculations for the system containing 2133 water molecules with dimensions of 40 Å along the x and y directions will not suffer from the artefacts observed for small systems. Furthermore, a comparison of the pure water-water radial distribution function (see below) with those for solutions of AN reveal the perturbation of the hydration shells by the salt.
We calculated the surface tension of the solutions by evaluating the normal (p N = p zz ), and tangential (p T = 0.5(p xx + p yy )) components of the Kirkwood-Buff (KB) stress tensor. The surface tension follows from the expression below with L z denoting the length of the simulation cell in the longest direction and p ii the ii th component of stress tensor. The h. . .i refers to time average and the factor 0.5 accounts for the presence of two symmetric interfaces.
We computed the radial distribution function (RDF, g(r)) in the bulk and at the interface of solutions from the recorded trajectories of the molecules, saved during the production runs. In addition, we calculated the coordination number, as function of distance from the central atom, from the RDF formula where r 0 denotes the lower limit of integration where g(r) disappears. For the first coordination number, r 1 stands for the first minimum in the RDF profile. Finally, the number density, r, for NO 3 À and NH 4 + is listed in the last column of Table 2, and for water molecules it can be easily calculated from the information provided in the same table.

Interfacial region detection
There are various methods of different suitability and quality for detection of the interfacial region. [77][78][79][80] In this study, we considered the Gibbs dividing surface (GDS) as a convenient point of reference to analyse the perturbations in molecular distributions at the liquid-vapour surface in detail. The location of the GDS can be obtained by fitting the total density profile, r(z), to a hyperbolic tangent function commonly used 81 of the following form where r liq and r vap define the bulk liquid and vapour density, z 0 is the z-position of the Gibbs dividing surface (GDS) parallel to the interface where the surface excess for the total density is zero and is a parameter related to the interface thickness in which the total density varies from 10% to 90% of its bulk value by the relation d t 10-90 = 2.1972d t . In the current study, r liq and r vap were both estimated from the MD computations. We have presented the mass and charge density profiles relative to the GDS, set at z = 0 along the longitudinal direction. We defined the interface (À2.5 Å o z o 4 Å), subsurface (À7 Å o z o À2.5 Å) and bulk slabs (z o À7 Å) as illustrated in Fig. 2, with the interface region comprising the GDS.  We observed that, the predictions of the OPLS model, without the electronic continuum correction for ions, for the three water models, agree quite well with the experimental values for concentrations up to 2 mol kg À1 . However, for higher molality, these models overestimate the density of the solutions.

Bulk density of AN solutions
The results of the OPLS/ECC model, with the selected force fields of water, demonstrate excellent agreement with the experimental data. 82,83 The maximum percentage error decreases from 8.3%, 9.5% and 9.9% to 1.7%, 2.7% and 1.5% for the SPC/E + OPLS, SPCE/F + OPLS and TIP4P/2005 + OPLS models, without and with the ECC, respectively. These models comprise different ion-water interactions but the same ion-ion interactions. The ion-ion interactions become increasingly crucial as the concentration of the salt increases. 84 Neyt et al. 18 also reported a reduction in the calculated density values with the ECC version of the Reif and Hünenberger 85 and OPLS 86 models, compared with the original models for the aqueous solutions of sodium chloride.

Structure of the bulk solvation shell
We studied the local structure of the bulk AN solutions by investigating the ion-water, ion-ion and water-water radial distribution functions (RDFs) for the OPLS/ECC force fields, and present the RDFs results for different molality of bulk AN solution and the TIP4P/2005 model at 298.15 K. Functions corresponding to the SPC/E and SPCE/F models have not been plotted since they do not differ significantly from those for the TIP4P/2005. Fig. 4 illustrates the calculated ion-water and pure water-water radial distribution functions obtained in a separate simulation (TIP4P/2005), with the nitrogen atoms in the nitrate and ammonium shown as N n , N a , respectively, and oxygen and hydrogen atoms of water presented as O w , H w .
The first maxima of g N n -O w and g N n -H w in the TIP4P/2005 + OPLS/ECC model in our simulation stand at 3.48 and 2.77 Å, in that order. These results indicate the distance of the hydration shell around the nitrate ions. The doublet character of the first maximum of g N n -O w , with a shoulder at 3.81 Å, corresponding to the second water shell, is also distinguishable in the studies of Dang et al. 87 and Thomas et al., 8 who applied a polarisable force field for nitrate ion in aqueous solution. The two peaks in the g N n -H w suggest that, one of the hydrogens in a water molecule points toward the nitrate ion and the other away from it (see Fig. 4). As the concentration of AN increases, the number of water molecules in the first hydration shell decreases, and the first minimum in g N n -O w becomes more pronounced, suggesting the formation of crystal-like domains in the solutions. In addition, g N a -O w and g N a -H w display the first peaks at 2.87 and 3.35 Å, correspondingly. Evidently, oxygen atoms in water molecules in the hydration shells of ammonium ions point toward the ammonium and both hydrogen atoms in a water molecule position themselves away from NH 4 + .
Comparing with g O w -O w of pure water presented in the left panels of Fig. 4, the hydration shell of g N n -O w is less prominent and is not as well-defined as the one of the g N a -O w . In addition, nitrate ions introduce a longer-range perturbation and the second hydration shell emerges as a shoulder of the first shell. Table 3 contrasts the distance of the first bulk hydration peaks for different pairs obtained from the TIP4P/2005 + OPLS/ ECC model employed in this study with their corresponding values published in the literature. Good agreement exists between our values and those reported by other researchers.  Therefore, we remark that the models for nitrate and ammonium ions implemented in this study accurately reproduce the structure of the hydration shells in the bulk that surround the solvated ions. The increased formation of ion pairs in concentrated solutions of AN represents the reason of slight drop in the intensity of the first hydration shell for the higher ion concentrations. In other words, the hydration numbers decrease with the increased formation of ion pairs at higher concentrations of AN. Table S3 (ESI †) lists the calculated coordination number (CN) for the first minimum in g(r) of Fig. 4. The force field used here accurately reproduces the structural properties of hydrated ammonium and nitrate ions. Similarly, to the previous studies, the original non-polarisable force field (OPLS) yields a more structured hydration shell than the OPLS/ECC model. 94 This is evident in sharpening of the peaks in RDF for non-ECC force fields. Fig. S3 of the ESI † presents a comparison between the calculated RDF from the OPLS and OPLS/ECC force-field models.
In Fig. 5, we present the ion-ion pair correlation functions (N n -N a , N n -N n and N a -N a ). The N n -N a RDF shows a pronounced peak at around 3.50 Å. The maximum of the peak decreases from 3.53 Å to 3.48 Å as molality increases from 0.65 to 18.7 mol kg À1 . These correlation functions indicate the presence of ion pairs in the solutions. Ions are paired if the distance between them is less than the cut-off radius of 4.87 Å. The cut-off radius is often selected as the first minimum of the pair correlation function between ions that carry the opposite charges. 95 The probability of ion-pair formation increases with increasing the concentration of the salt. The presence of two peaks in RDFs of N n -N n and N a -N a reveals the formation of the  Table 2.  88 3.65 2.75 Dang et al. 87 3.4 2.4 Banerjee et al. 89 3.5 2.6 Thomas et al. 8 3.5 2.6 Vchirawongkwin et al. 90 3 crystal-like domains, with increasing salt concentration in the solution. In addition, the coordination number increases with concentration due to formation of more ion pairs. The RDFs for the N n -N n and N a -N a pairs, at 18.7 mol kg À1 , are in excellent agreement with the RDF of melted pure AN at temperatures of 550 K in that Velardez et al. 96 study that investigated the structural changes of solid-state of AN using molecular dynamics. These authors reported the maximum of the first peak in both ion pairs (N n -N n and N a -N a RDFs) to appear around 4.6 Å. The RDFs calculated in our study exhibit the maxima of the first and second peaks for N n -N n pairs at 4.73 and 6.87 Å and for N a -N a at 4.64 and 7.03 Å, respectively. Wahab and Mahiuddin 82 concluded that, the free hydrated ions of NH 4 + and NO 3 À exist up to the molality of 9.1 mol kg À1 , whereas the solvent-separated ion pairs form in the concentration range of 9.1-12.0 mol kg À1 . Finally, beyond 12.0 mol kg À1 , a transition to the solvent-shared ion-pairs occurs as a result of a decrease in the number of available water molecules. Fig. 6 demonstrates lack of a well-defined second hydration shell at high ion concentrations of 13.6 and 18.7 mol kg À1 . This occurs because of the presence of ions in the interstitial positions, concurring with a similar behaviour observed for highly concentrated solutions of NaCl and KCl. 97 The function related to the low salt concentration are comparable with the one of the g O w -O w of pure water.

Surface tension of AN solutions
We computed the surface tension of AN solutions over the studied range of molality (0.65-18.7 mol kg À1 ). Fig. 7 compares the surface tension of AN solutions with the OPLS and the OPLS/ECC models for the ions, and the three models of water (SPC/E, SPCE/F, and TIP4P/2005). The experimental results of the surface tension of AN solutions at 298.15 K were obtained from the interpolation between the surface tension values of AN   Fig. 7.
The error between the computed and the experimental measurements arise due to underestimation of the surface tension of pure water by TIP4P/2005 compared with the experimental measurements. The surface tension calculated with the SPC/E, SPCE/F and TIP4P/2005 models underpredicted the measurements by 15%, 4.5% and 6.6%, in that order, in comparison to the experimental figure of 71.99 AE 0.05 mN m À1 at 298.15 K. 76 In a summary evaluation, only the TIP4P/2005 model in conjunction with the OPLS/ECC description of the AN ions provide precise predictions of the surface tension of AN solutions with the accuracy needed for industrial applications over the entire range of concentrations. The simulated results account well for the experimental trend in the measurement of surface tension with increasing concentration of the solute. Table 4 lists the computed surface tensions of AN solutions for different salt molality at 298.15 K for the OPLS/ECC models and different models of water compared with experimental values. Industrial applications, such as AN emulsion explosives, involve supersaturated solutions of AN, which properties are often outside the reach of the present experimental techniques. As far as we can ascertain, this finding is new and of great practical relevance.

Density profile
To quantify the physical behaviour of ions and their occurrence at the surface or in the interior of the slab, we plotted in Fig. 8 the water, ion and total density along the z direction for three concentrations under study (2.5, 8.6 and 18.7 mol kg À1 ) relative to the GDS; all results represent the use of the TIP4P/2005 water model. The GDS is positioned at z = 0 along the direction perpendicular to the open surface. The solid and dash lines reflect the simulation results for the OPLS and OPLS/ECC models, respectively, with charge density from the latter reporting to the secondary y-axis. To calculate the quantities displayed in Fig. 8, we have divided the simulation box into the vertical slabs of 0.2 Å in thickness.
We found that, the approach of ions to the interface and their separation differ significantly as predicted from the OPLS and the OPLS/ECC models. For the former, ions repel from the interface and form a depletion layer with a thickness of approximately 5 Å and the separation between layers of ions is rather small. While for the latter, the ions display pronounced segregation, with nitrate anions preferring to place themselves at the interface. This leads to the depletion of NO 3 À in the subsurface, resulting in the formation of a layer preferentially occupied by the cations. Fig. 8 illustrates the structure of this ionic double layer, arising for all concentrations of AN considered in this study, with weak affinities of anions and cations for the surface and subsurface, respectively. The predictions from the OPLS/ECC model reproduce well the results of the polarisable force fields for the ions. 100 We demonstrate the ion deficiency up to 3 Å below the surface of the solution and that polarisability effects  Although the previous investigation concluded that the surface excess of one type of ions may not necessarily induce the emergence of an electric double layer, we did not observe this behaviour in our investigation. Table 5 lists the values of density of the liquid and vapour phases, together with thickness of the interface, as computed from the hyperbolic tangent function (see eqn (10)) for all studied concentrations of AN. We observed a monotonic increase in the width of the interface with the rising concentration of solutions of AN. For the low-concentration solution (o2.5 mol kg À1 ), ions maintain the same number of water molecules in the first hydration shell in the interior of the slab. However, the present calculations establish that, ions at the interface hydrate with  fewer water molecules than those in the interior. This behaviour stems from the absence of water molecules in the vacuum regions above or below the slab. For the concentrated solutions (e.g., 8.6 and 18.7 mol kg À1 ), more water molecules exist in the first hydration shell that surrounds the ions. This behaviour arises in the region between À8 Å o z o À2 Å, compared with the interior due to the lower abundance of ions and a higher number density of water molecules in this region. The declining coordination number toward the surface reflects the decreasing number of water molecules present closer to the vacuum region. In addition, Fig. 9 reveals higher probability of detecting unpaired ions at the interface compared with the interior region of the slabs. This is a surprising finding considering the behaviour of germane systems revealed in literature. Minofar et al. 24 reported that, the ion pairing increases for Mg(OAc) 2 (aq) at the water/air interface compared to the bulk of the solution. According to the authors, the dominant reason for this behaviour is the intense attraction of acetate anions to the interface and, consequently, their pairing with magnesium cations present preferentially in the subsurface. This intense attraction induces powerful partitioning of ions away from the bulk into the interface that the authors described as the surfactant effect, as the phenomenon is associated with decreasing surface tension as function of the nominal concentration of the salt. More importantly, from the perspective of our results, Minofar et al. also concluded that, compared with Mg(OAc) 2 (aq), Mg(NO 3 ) 2 (aq) behaves differently and displays less affinity for ion pairing (no detectable pairs at 1 M concentration), but still the number of ion pairs increases at the interface compared with the interior, for the solution concentration of 2 M, as illustrated in Fig. 7 and 8 of ref. 24. In contrast to these findings, our simulations indicate that, while increased salt concentration induces higher ion pairing in the bulk, there is always decrease of ion pairing at the interface in comparison to the bulk. For instance, for the molality of 18.7 mol kg À1 , on average, every nitrate anion pairs with B4.3 and B1.3 ammonium cation in the interior and at the interface of the slab, respectively (see Fig. 9). The variation of ion pairing at the interface with respect to the bulk is dependent on the affinity of ions to the surface and the degree of cations and anions segregation in the interfacial layer. Moreover, the ion deficiency at the interface coincides with higher population of water molecules in the subsurface and in part of the interface layer, which prompts declining of ion pairing in that region (see Fig. 8). Fig. 9 demonstrates that although NO 3 À has the affinity for the surface, this region is not completely depleted of the NH 4 + cations. The results from the OPLS force field illustrate the same trend; albeit, for different average coordination numbers as discussed earlier. Also, the OPLS/ECC model yields the higher affinity of anions to the surface as compared with the OPLS force field. Lack of a plateau in the interior section of concentrated solution could be attributed to the nonhomogenous distributions of ions in bulk (see Fig. 8). From Fig. 9, one infers that larger ions (nitrate) show the greatest changes in the coordination number upon transfer into the vapour phase due to the elevated number of water molecules in the first solvation shell. Bauer et al. 101 explained the reason for this comportment by weaker attraction of water to larger ions. We also found that, the general features of the RDFs remain similar from the interior and the interface region (see Fig. S5 of the ESI, † where we present the functions for the salt molality of 18.7 mol kg À1 ). This results in the propensity of the ions to maintain the hydration shell structure and molecular interactions even at the interface. Snapshots of hydration shells (defined as o5.1 Å) of the nitrate ion for the AN molality of 18.7 mol kg À1 in the bulk and at the interface illustrate the water depletion zone at the interface (Fig. S6, ESI †).

Conclusions
The present contribution computed the bulk and surface properties of solutions of ammonium nitrate (AN), over a broad range of concentrations using a non-polarisable and implicitly polarisable models for ions (OPLS and OPLS/ECC, in that order) and three non-polarisable models for water molecules (SPC/E, SPCE/F and TIP4P/2005). The OPLS model overestimates the density and the surface tension of AN solutions and is unable to reproduce the experimental values. However, a simple and computationally efficient correction to non-polarisable simulations demonstrate that, the amended model works well for the bulk and liquid-vapour interface of AN solutions.
Three water models in conjunction with OPLS/ECC provide good predictions of the bulk properties (density and solvation structure) of solutions of ammonium nitrate. However, only TIP4P/2005 model yields correct trends in the surface tension and provides accurate estimates of the interfacial properties of the solutions in comparison to the experimental measurements. Density profiles reveal the repulsion of ions from the surface, as calculated from the OPLS model, while the OPLS/ECC model, similarly to polarisable models, induces the migration of anions toward the surface, and deficiency of negative charge in the sublayer of the surface. This behaviour generates an electric double layer that explains the increasing surface tension as the concentration of ammonium nitrate increases in the solution. This finding is of significant practical importance as it affords the calculation of the interfacial properties of supersaturated solutions of AN, with measurement of such properties presently outside the reach of the current experimental techniques. Yet, the aqueous (discontinuous) phase of AN emulsion explosives consist of supersaturated solutions of AN and knowledge of such properties is critical to improve the stability of emulsion explosives.
We evaluated the degree of ion pairing as function of concentration of ammonium nitrate. The hydration of ions by water molecules decreases, as the number density of nitrate and ammonium ions increases in the solution. While negligible at lower concentrations, ion pairing becomes significant as the concentration rises. However, we discovered that the ion pairing decreases at the interface in comparison to the bulk.

Conflicts of interest
There are no conflicts to declare.