Mechanisms underlying the freezing and melting behavior of water confined in silica nanopores as a function of pore size and pore filling

Abstract

Uncovering the mechanisms of freezing and melting behavior in nanoconfined fluids can unlock fundamental insights into the fate and transport of fluids in soils present in cold climates. From a scientific perspective, the structural and thermodynamic behavior of confined and interfacial water has sparked significant discussions, particularly regarding the characteristics of phase transitions and spatial heterogeneity as a function of temperature and pressure. Observations frequently report interfacial unfrozen liquid layers on hydrophilic surfaces, distorted ice crystals and suppressed freezing and melting points in confined water compared to bulk water. These effects are often attributed to the restricted molecular mobility and the influence of pore surfaces. However, the exact nature of these phase transitions and the specific characteristics of the layered arrangement remain uncertain. In this study, we present an approach to elucidate the layered structural arrangement and phase transitions of confined water through integrated thermal analysis and molecular dynamics (MD) simulations. By employing differential scanning calorimetry (DSC) on samples with precisely controlled pore fillings achieved through the vapor loading method, we reveal evidence of the formation of sequential water layers and estimate the individual layer widths. The observed reduction in the enthalpy of fusion during freezing relative to bulk water, alongside distinct heat flow peaks, indicates the potential occurrence of weak first-order transitions in the frozen layers. Experimental findings are further supported by classical MD simulations conducted on analogous systems confined in amorphous silica slit pores with widths ranging from 4 to 8 nm. The combined influence of cooling rate, surface hydrogen bonding, and non-bonding interactions between silica and confined water critically impact the development of interfacial ice polymorphs.

---

Greeshma Gadikota

Dr Greeshma Gadikota is an Associate Professor and Croll Sesquicentennial Fellow in Civil and Environmental Engineering at Cornell University. She directs the Sustainable Energy and Resource Recovery Group. Her postdoctoral training was at Princeton University and Columbia University. Her PhD in Chemical Engineering and MS degrees in Chemical Engineering and Operations Research are from Columbia University. Her BS in Chemical Engineering is from Michigan State University. Her scientific contributions are recognized by the DOE, NSF and ARO CAREER Awards, Sigma Xi Young Investigator Award, AICHE Sabic Award for Young Professionals, and ACS Women's Chemists Committee (WCC) Rising Star Award.

---

1. Introduction

Uncovering key properties of confined and interfacial water such as freezing and melting points has far-reaching implications across diverse fields, including environmental science,¹ biomedical research,² and chemical processes.³ The ability to architect materials with uniform pore size distributions^4,5 enables detailed investigation of the distinct phase behavior of fluids in confinement.^6–9 Water, when confined in nanopores, exhibits several unusual properties. The thermodynamic variable that is highly affected by the confined spaces is pressure due to the anisotropy offered by the confinement. At thermal equilibrium, the difference between the parallel and perpendicular components of the pressure tensor affects the interfacial tension of the solid–liquid interface, which in turn, influences the confined water structure.¹⁰ As a result of the solid–liquid interaction and pressure-induced distortion in the hydrogen bonded water molecules, the density of water is often observed to be reduced by 15–20% compared to that of bulk water in materials with pore diameters of 2.2–7.7 nm.¹¹ Another variable that significantly influences the structural and dynamical characteristics of confined water is temperature. For example, to deduce damage mechanisms in cement and concrete that consist of a large array of various sized nanopores as a function of temperature, Bonnaud and co-workers¹² investigated the water properties in the pores of calcium–silicate–hydrate (C–S–H) within a temperature range of 100 K to 575 K. A rigorous characterization of excess potential energy, excess entropy, and mean water density in their study revealed a liquid–liquid transition of the confined water within 180 K to 195 K followed by a glass transition at ∼170 K.¹² Differential behavior of confined water has also inspired researchers to examine the effect of confinement over several key biological processes associated with crucial biomolecules such as proteins that exhibit strong variability in folding pattern depending on the biological milieu, and the surrounding spatial, and thermodynamic constraints. To elucidate the effects of the latter two on the folding stability and initial disorder of a protein, Layek and co-workers¹³ reported the response of an albumin binding domain at lower temperatures within an armchair single-walled carbon nanotube confinement. Interestingly, the structural fluctuations in the protein that directly result during the unfolding episode, were suppressed due to the confinement along with a direct impact of thermal energy on the protein's folding propensity.¹³ Analogous investigations of the behavior of water in nanoscale confinement revealed several key phenomena underlying the differences in the melting and freezing points compared to bulk water. Zaragoza and co-workers¹⁴ assessed the diffusion and viscosity of water in the pores of carbon nanotube that led to the development of a confined Stokes–Einstein relation. At temperatures ≤263 K, the diffusion coefficient of confined water was found to resemble that of bulk water; while at high temperatures, the values were higher than their bulk counterparts. In fact, viscosities computed from diffusivities yielded by confined Stokes–Einstein relation differed significantly from that obtained using traditional Green–Kubo relationships. However, the authors observed the viscosity to be unaffected by the confinement geometry.¹⁴ These findings motivate the development of a rigorous research approach that links the phase transitions of confined water to the underlying molecular scale mechanisms.

Contrary to the relatively simpler phase transition of bulk water, phase transitions of water in nanopores is considered to be more complex with significant shifts in freezing and melting points due to the interfacial interactions with the pore surface.^15–17 Early theoretical models, such as the Gibbs–Thomson (G–T) equation attempted to describe the relationship between pore size and freezing and melting point depressions.^18,19 However, deviations observed between experimental data and theoretical predictions have led to modifications of the G–T equation, including the incorporation of an unfrozen layer at the pore surface.²⁰ Fitting phase transition onset points obtained from thermal or kinetic measurements with modified G–T equations that include a 0.3–0.8 nm unfrozen layer is widely practiced to characterize the unfrozen water layer in confined environments.^6,8,20,21

Advanced characterization techniques, such as differential scanning calorimetry (DSC), nuclear magnetic resonance (NMR), broadband dielectric spectroscopy (BDS), small- and wide-angle X-ray scattering (SAXS/WAXS), along with molecular dynamics (MD) simulation methods, have been utilized to reveal critical insights into the thermodynamic parameters, kinetic processes, structural organization, and molecular interactions characteristic of confined water.^{6–9,20,22,23} Growing experimental evidence and simulation studies support the existence of unfrozen interfacial water layer in hydrophilic confinements. Observations from NMR and BDS provide compelling evidence of unfrozen water persisting at low temperatures within nano-confinement, aligning with predictions made in theoretical models.^9,17,24 In addition to distinguishing between unfrozen interfacial layers and confined core water, more sophisticated multi-layered models have been proposed to elucidate the complex phase transition behaviors observed in the thermal analysis of confined water within pores larger than 2.2 nm. Various methodologies and simulation approaches suggest the presence of 2–4 structural layers, with dynamic shifts in these layers reported under different experimental conditions.^8,11,15 Among these studies, Xia and co-workers, through DSC and NMR studies, proposed a three-layer freezing model in materials with pore diameters of 4–15 nm to elucidate the fragile-to-strong dynamical crossover observed in supercooled water.⁶ The intricate nature of nano-confined water phase transitions challenges the classification of these transitions as either first-order or glass transitions, with emerging evidence indicating that confined water undergoes gradual, continuous transformations in ice structure.²⁵ For example, Erko and co-workers highlighted that the characteristics of confined water cannot be directly inferred from bulk water behavior, with phase transitions showing unique intermediate states.²⁶

Understanding the layer-by-layer or core-to-surface freezing mechanisms is another crucial factor for elucidating the structural and dynamic evolution of confined water during phase transitions. One promising approach involves capturing thermal snapshots of confined water samples with precisely controlled pore fillings, which could reveal discrete phase transition peaks corresponding to each layer, thus providing valuable insights into the sequential nature of these transitions. Additionally, MD simulations offer a molecular-scale perspective on the evolution of structures of water clusters as a function of the most crucial thermodynamic variables such as temperature and pressure, and the role of solid–liquid interaction in influencing the fate of the confined liquid, complementing experimental observations.^7,15,22In silico studies also provide detailed information on the coordination behavior of water molecules in confinement that is unattainable directly from experiments. These studies also offer insight into the behavior of confined water as a function of the cooling rate, which has direct consequences on the resultant polymorphs.²⁷ In this study, we propose to leverage vapor-loading techniques for pore filling to achieve accurate control over sample composition²⁸ enabling robust analyses of phase transitions in confined water. We segregate the regions of structurally different water layers, thoroughly investigate their coordination environment, and quantify the contribution of the solid surface in modulating the structure of water through non-bonding interactions. Employing MD simulations in tandem with DSC measurements (Fig. 1) enhances our understanding of the dynamic structural evolution of confined water and informs hypotheses about layer-specific phase transitions in nano-confinement.

Fig. 1 Snapshots of water-filled amorphous silica slit pores used in molecular dynamics simulations, with widths of (a) 4 nm, (b) 6 nm, and (c) 8 nm. Silicon, oxygen, and hydrogen atoms on the silica surface are shown as yellow, red, and white vdW spheres, respectively. Water molecules are depicted as blue lines.

Therefore, the three research questions this work addresses are: (1) How does pore filling and pore size influence phase transition behavior of nano-confined water and layer-by-layer freezing and melting behavior? (2) How does the coordination environment of water vary in the layered structures of nano-confined water? (3) What is the energetic basis underlying the phase transition behavior of nano-confined water? Addressing these questions is crucial for unpacking complex phase transition behavior of water in confinement on freezing and melting. To summarize, this study offers several new insights that are not discussed in prior studies. First, the layer-by-layer freezing behavior as a function of pore filling and pore size are determined using differential scanning calorimetry (DSC) measurements. Second, the molecular scale basis underlying the observed freezing and melting behavior in nanoscale confinement from the perspective of the coordination environment, hydrogen bonding, electrostatic and van der Waals interactions, are discussed.

2. Materials and methods

2.1. Confined water sample preparation via vapor loading and characterization using differential scanning calorimetry (DSC)

SBA-15 samples with pore diameters of 4, 6, 8, and 12 nm obtained from Sigma-Aldrich are used for loading water into these materials. The pore size distributions of SBA-15 are determined using an Anton Paar ASIQ C-MP model BET surface analyzer via nitrogen adsorption–desorption isotherms. The pore size distributions are calculated using the nonlocal density functional theory (NLDFT) method²⁹ (see Table S1). The vapor loading method adapted from Agrawal and co-workers²⁸ is used to fill the SBA-15 nanopores. Deionized water is prepared using a Milli-Q Advanced A10 model. To achieve vapor loading, 10 mL of DI water is added to a 35 mL sealed container, in which an open vial containing 0.1 g of SBA-15 is placed. The sealed setups are subsequently immersed in a water bath maintained at 308 K using a heating plate for samples with pore diameters of 4, 6, and 8 nm, while samples with a pore diameter of 12 nm are kept at ambient temperature. By monitoring the change in sample weight (m_sample) over time (Fig. 2), the variation in pore filling with loading time at a fixed temperature is obtained. The pore filling level (ϕ) is calculated using eqn (1).

Fig. 2 Variation of normalized pore filling level (volume of loaded water divided by the pore volume of SBA-15) during the vapor loading process. The water volume was calculated from weight changes using the bulk density. The inset shows a magnified view of the first 70 hours.

m _water represents the weight of the confined liquid at each time instant, ρ_bulk is the bulk liquid density, and V_SBA-15 represents the measured pore volume. ϕ is calculated by normalizing the confined liquid volume at each time point to the fully filled pore capacity (V_SBA-15). Density of the confined liquids is assumed to be equivalent to the bulk liquid density (ρ_bulk), as indicated by previous research.⁷ The freezing and melting phase transitions of the confined liquids are analyzed using differential scanning calorimetry (DSC) on a TA Instruments DSC-25 model with a ramp rate of 5 K min⁻¹ from room temperature down to 103 K, followed by equilibration at 103 K, and subsequent heating to room temperature at the same rate.

2.2. Analysis of DSC data

The onset temperatures of the DSC curves during freezing and melting are defined as the freezing point and melting point, respectively.²⁰ The sample-specific freezing and melting enthalpies (ΔH_process) are obtained by integrating the areas under the respective curves.³⁰ The phase transition enthalpy of confined water is obtained by dividing the sample-specific phase transition enthalpy by the weight proportion of confined water in the sample (α_water), as shown in eqn (2) and (3) below.

As confined water consists of both freezable and unfrozen water, samples that do not exhibit phase transition peaks (above 103 K) are considered to contain only unfrozen interfacial water. The corresponding pore filling levels of these samples represent the volume of unfrozen liquid. The fraction of freezable water (α_freezable) is determined using eqn (4).

In the sample group with the highest pore filling level (ϕ) containing only interfacial water, the specific pore filling level (ϕ_int) is used as a reference to calculate the fraction of freezable water (α_freezable) in samples with higher pore filling levels (ϕ > ϕ_int).

The normalized phase transition enthalpy is calculated using eqn (5) by comparing the phase transition enthalpy of confined freezable water with that of bulk water to quantify confinement-induced enthalpy suppression.

is the freezing or melting enthalpy of confined water. ΔH_bulk is the enthalpy of fusion of bulk water. ΔH_normalized is obtained by normalizing to the enthalpy of bulk water (ΔH_bulk) to investigate the order of transition in freezable water, thereby elucidating the nature of the phase transition (eqn (3)). The enthalpy of fusion for water (ΔH_bulk) used is 334 J g⁻¹.³¹ The unfrozen water, shown to have no contribution to the enthalpy change, is excluded from the calculation by introducing the term ‘fraction of freezable water’ (α_freezable)³² defined in eqn (4). To estimate the layer widths observed during the freezing process, we simplified the formation of confined water as a radially expanding ring from the pore surface with increasing pore filling. The relationship between the pore filling level (ϕ), pore diameter (D) (in nm), and the corresponding ring width (d) (in nm) is shown in eqn (6).

2.3. Simulation methodology: construction of slit pores, force field parameters, and molecular dynamics workflow

Slit pore models of amorphous silica are built to replicate the amorphous SBA-15 nanopores used in the experiments. Amorphous silica membranes are built starting from crystalline silica in the following steps. First, β-cristobalite structure of SiO₂ is obtained from Crystallography Open Database (COD ID: 1010944) with unit cell parameters as a = b = c = 7.12 ± 0.01 Å and α = β = γ = 90°.³³ The unit cell is then replicated in space by 8 × 6 × 6 to obtain a larger membrane of crystalline silica with lengths along x, y, and z dimensions as approximately 5.7, 4.3, and 4.3 nm respectively. This crystalline silica membrane is then made amorphous following the methodology provided by Stallons and co-workers.³⁴ First, the membrane structure is relaxed using the steepest descent minimization algorithm for 50 000 steps and then melted at 6000 K during the first 3000 NVT MD steps. This is followed by a 7 ps NVE run. A more realistic appearance of the modeled silica membrane is obtained by adding hydroxyl functionalities and H atoms to the valency unsaturated Si and O atoms respectively on the amorphous silica membrane, replicating the intrinsic hydroxyl groups observed on SBA-15 surfaces.^35,36 This addition is confined within 5 Å from the top of the silica surface.³⁷ The bonding and non-bonding potential in surface hydroxylated silica are modeled using ClayFF.³⁸ Therefore, the net charge obtained after hydroxyl addition is neutralized by distributing the excess charge among the Si atoms furthest from the exposed silica surface.³⁷ We call this membrane as SIL_hpl which is further duplicated in space to create slit pores of widths 4, 6, and 8 nm (see Fig. 1 and Fig. S1). It is important to note that in silico construction and simulation of a slit pore of width 12 nm is challenging in terms of system size and therefore not modeled in the present study, rather the phase transition behavior of water in such a large confinement is evaluated solely from the experiments.

In line with previous literature reports, the modeled slit pores are filled with OPLS compatible TIP5P water assuming water density same as its bulk density at 25 °C.^7,39–41 Each of the water filled slit pores is then energy minimized for 50 000 steps using steepest-descent algorithm followed by 50 ps of NVE run to balance the shear forces and 50 ns of NVT run to equilibrate the system temperature. Three starting structures are randomly sampled from the last 10 ns of the 50 ns long NVT run to conduct simulation triplicates using widely different random number seeds. Statistics are accumulated for all the systems from an additional 79 ns NVT run with simulated annealing. The annealing temperatures are selected based on the end set points of each phase transition observed in the freezing and melting profiles of neat water (Fig. 3 and Fig. S2, S3). A single sequence of annealing points is used with the control points at 298 K, 232 K, 225 K, 198 K, 103 K, 217 K, 260 K, and 298 K with a ramp rate of 10 K ns⁻¹. The chosen ramp rate is a trade-off between the computation time and the ramp rate-independent relaxation behavior of supercooled water.^42–44 In all the NVE and NVT simulations, a timestep of 1 fs is used to solve Newton's equations of motion using the leap-frog integrator. All bonds containing H are constrained. Neighbor searching for non-bonding interactions is completed using a Verlet cut-off scheme and the neighbor list is updated every 10 MD steps. Short-range electrostatics and vdW interactions are cut off at 12 Å with long-range dispersion correction for energy and pressure. Long-range electrostatic interactions are calculated using the Particle Mesh Ewald (PME) method with a grid spacing of 0.16 nm along with cubic interpolation.⁴⁵ Temperature coupling is achieved by velocity rescaling with a stochastic term using a coupling time constant of 0.1 ps and a reference temperature of 298 K. Periodic boundary conditions are applied along three dimensions of the simulation cells. GROMACS 2024⁴⁵ and VMD⁴⁶ are used as the simulation software and to visualize the trajectories, respectively. VMD, MATLAB (R2023b), and gmxtools⁴⁷ are used to analyze the trajectories and extract key molecular scale insights regarding the coordination environment.

Fig. 3 Heat flow data from differential scanning calorimetry (DSC) measurements during a freezing–melting cycle of water confined in SBA-15 materials with pore diameters of 4 nm. The black curve corresponds to pore fillings below 33%, where only interfacial water is present with no phase transition peaks. The red and green curves represent the formation of the water shell and core layers, respectively. During melting, all profiles exhibit a single peak, indicating simultaneous melting of the water layers. Schematics depict the interfacial, shell, and core layers in coral red, orange, and brown, respectively. Vertical dashed lines at 233 K, 239 K, and 245 K indicate the freezing points of the shell and core layers and the melting point of freezable water.

3. Results and discussions

3.1. The influence of confinement size and pore filling on freezing behavior of nanoconfined water

Multiple phase transition peaks during the freezing process and a peak during the melting are observed in water confined in pore sizes ranging from 4–12 nm (Fig. 3, 4 and Fig. S2–S4). The precise pore filling control achieved via vapor loading enables detailed observation of each phase transition. For confined water samples that did not exhibit phase transition peaks, the confined, unfrozen water is interpreted as the interfacial layer as noted when pore filling is below 31%, 26%, 22%, and 13% for materials with pore diameters of 4 nm, 6 nm, 8 nm, and 12 nm, respectively.

Fig. 4 Heat flow data from differential scanning calorimetry (DSC) measurements during a freezing–melting cycle of water confined in SBA-15 materials with pore diameters of 4, 6, 8, and 12 nm corresponding to their respective highest pore filling. During melting, all profiles exhibit a single peak, indicating simultaneous melting of the water layers. Vertical dotted lines in the freezing profiles indicate the freezing points of the shell and core layers and in the melting profiles, they correspond to the melting points of freezable water.

With increasing pore filling levels, two distinct phase transition peaks emerge sequentially. The first peak is indicative of the formation of a “shell” layer adjacent to the interface and the second peak corresponds to the “core” layer within the pore interior. This distinction is based on the sequential appearance of these peaks in the DSC profile, with the shell forming closer to the solid–liquid interface and the core occupying the pore center. The presence of multiple freezing peaks supports the hypothesis that water in hydrophilic confinements forms a layered structure during the freezing process.^6,48 Conversely, the single melting peak suggests that each layer transitions as a unified phase, further corroborating the layered structure hypothesis by providing evidence against the likelihood of multiple-phase transitions within a single layer. We attributed the layered structure to the solid–liquid interfacial interactions. Consequently, water molecules experience different degrees of influence from the interfacial forces, leading to distinct phase transitions for each layer.⁴⁹

For confinements with varying sizes, our observations suggest that the formation and structural characteristics of confined water layers in 4–8 nm pores are predominantly governed by the surface effect across all layers, whereas in the 12 nm pore, only the interfacial layer is primarily influenced by the surface, while the shell and core layers are influenced by the pore size. Based on the freezing points of individual layers (Table 1), confined water across all pore sizes exhibits a core-to-surface freezing pattern, consistent with findings reported in previous studies.^9,15 For confinement with specific pore filling values (ϕ_int), the absence of any phase transition peaks as shown by the black curves in Fig. 3 and Fig. S2–S4 implies the existence of unfrozen water. This unfrozen water is hypothesized to form as an interfacial layer on the hydrophilic confinement surface, predominantly influenced by hydrogen bonding between water molecules and hydroxyl groups on the silica surface.^50,51 Notably, ϕ_int decreases with increasing pore size. By calculating the thickness of this interfacial layer based on the fraction of unfrozen water via solving eqn (4), we estimated its width (d) to be approximately 0.35–0.43 nm (Table 2), aligning closely with previously reported values.^20,52,53 Since the spatial distribution of confined water density remains uncertain, both bulk water density assumption and a confined water density model proposed by Etzler et al. for comparative analysis are used.⁵⁴

Table 1 Freezing and melting points of confined water across different pore diameters, pore filling levels, and layered structures. For confined water in SBA-15 materials with pore diameters of 4–8 nm, a three-layered water structure forms and stabilizes as pore filling increases. In contrast, for confined water in SBA-15 materials with pore diameter of 12 nm, phase transition peaks indicate a three layered structure when pore filling exceeds 65%. However, the noticeable shift of Peak II, corresponding to the water core, suggests the formation of a core with bulk-like properties as the number of water molecules increases and is marked with asterisks. When pore filling exceeds 83%, Peak I disappears and a new Peak II emerges, indicating the merging of the original water shell and water core. “NA” indicates the non-availability of freezing/melting points for systems containing only the interfacial water layer and thus no phase transition within the temperature range explored. Similarly, “—” indicates the absence of the corresponding layer (shell/core) at the given pore filling level

Freezing pattern	Pore diameter (nm)	Pore filling (%)	Freezing point (K)				Melting point (K)
			Interfacial	Shell Peak I	Core Peak II	Core* Peak II*
One water layer	4	<33	NA	—	—	—	NA
	6	<26	NA	—	—	—	NA
	8	<22	NA	—	—	—	NA
	12	<13	NA	—	—	—	NA
Two water layers	4	33–54	NA	232.0 ± 3.4	—	—	231.4 ± 1.3
	6	26–39	NA	231.6 ± 2.4	—	—	238.7 ± 8.1
	8	22–30	NA	232.4 ± 1.2	—	—	242.4 ± 6.7
	12	13–53	NA	237.6 ± 0.4	—	—	259.5 ± 1.9
Three water layers	4	>54	NA	233.6 ± 1.4	239.9 ± 0.6	—	245.0 ± 5.4
	6	>39	NA	234.0 ± 0.8	240.2 ± 0.3	—	253.1 ± 4.8
	8	>30	NA	233.6 ± 1.1	239.2 ± 0.3	—	256.4 ± 2.4
	12	53–65	NA	237.0 ± 0.3	243.1 ± 0.1	—	260.5 ± 0.9
Three water layers*	12	65–83	NA	237.3 ± 0.3	—	247.3 ± 0.6	259.2 ± 1.2
Two water layers*	12	>83	NA	—	—	250.5 ± 0.3	259.2 ± 0.7

---

Table 2 Widths of confined water layers at different confinement sizes obtained from experiments and simulations. The layer width of interfacial water and combined interfacial and shell layers in 4, 6, and 8 nm pores are calculated from ϕ_int in experiments using two estimation approaches: (1) assuming bulk water density, and (2) applying the confined water density from Etzler and co-workers (ref. 54). Simulation results for the combined interfacial layer and water shell are provided for comparison in (3)

Water layers		Layer width (nm)
		4 nm pore	6 nm pore	8 nm pore	Average
(1)	Interfacial water	0.35 ± 0.01	0.41 ± 0.02	0.41 ± 0.02	0.39 ± 0.04
	Combination of interfacial layer and the water shell	0.62 ± 0.03	0.64 ± 0.02	0.62 ± 0.02	0.62 ± 0.02
(2)	Interfacial water	0.38 ± 0.01	0.45 ± 0.03	0.45 ± 0.02	0.43 ± 0.04
	Combination of interfacial layer and the water shell	0.68 ± 0.03	0.70 ± 0.03	0.68 ± 0.02	0.69 ± 0.03
(3)	Combination of interfacial layer and the water shell	0.63 ± 0.06	0.67 ± 0.03	0.70 ± 0.06	0.67 ± 0.04

---

In the materials with pore diameters in the range of 4–8 nm, a pattern is observed wherein an initial freezing peak (Peak I) appeared with onset temperatures around 232–234 K, followed by a second peak (Peak II) with onset temperatures around 240 K (Fig. 3 and Fig. S2, S3). Peaks I and II correspond to the freezing of the water shell and water core, respectively.^54,55 As pore filling increased, the intensity of Peak II increases, while Peak I remains stable. However, when the pore diameter is 12 nm (Fig. S4), after the water shell formed at pore fillings above 13%, Peak II with onset temperatures around 243 K emerged not as a new peak, but by separating from the existing shell peak. When the pore filling exceeded 65%, Peak II shifted to a higher freezing temperature of approximately 247 K, suggesting an increase in the freezing point of the water core. At a filling level above 83%, Peak I disappeared, and the freezing point of the core continued to rise to 250 K. Similar dynamic evolution of layered freezing structures has been reported by Findenegg and co-workers and Schreiber and co-workers.^8,20 We interpreted the shift of Peak II in water confined in pores with diameters of 12 nm as an evidence that larger confinements allow confined water to behave similarly to bulk water, due to the reduced surface tension of water confined in pores with larger diameters.^11,56,57 As the core region expands within the 12 nm confinement, interactions between water molecules in adjacent layers can compete with solid–liquid interactions at the interface, leading to a transition from a three-layered structure to a two-layered structure. Conversely, in pore diameters of 4–8 nm, the three-layered structure remains stable across varying pore filling levels. Given the stable three-layered structure observed in the materials with pore diameters of 4–8 nm during freezing, we applied eqn (4) to estimate the combined layer width of the interfacial water and the water shell, which aligns with those obtained from simulations (Table 2). Similar to the observation of interfacial layer, we found that this combined width remained consistent across all confinement sizes, further supporting the surface-dominated structuring of the interfacial layer.^58–60

3.2. Melting points depression and phase transition enthalpy of water in confinement

The melting process of confined water is associated with a single-phase transition peak across all pore sizes and filling levels studied, as opposed to the multiple phase transition peaks observed during the freezing process. The widely reported melting point depression is evident in our results as well, with greater depression (ΔT_m) observed with decreasing pore size.^61,62 We applied 273.15 K as the melting point of bulk water (T_mb) under ambient pressure⁶³ for comparing with the melting point of confined water (T_m). For samples with saturated pores, a decrease in pore diameter from 12 nm to 4 nm led to an intensified melting point depression, decreasing from approximately −13 K to −28 K, underscoring the significant role of confinement size in phase transition behavior. Independently, within each pore of diameter 4, 6, and 8 nm, a notable melting point depression under reduced filling levels is noted. This depression, approximately 14 K lower than that of fully saturated pores, can be attributed to the formation of low-density liquid phases and the structuring of hydrogen bonds near the pore walls.²³ In contrast, for the 12 nm confinement, the melting point remained constant as pore filling increased, supporting the hypothesis that a 12 nm confinement has sufficient space to facilitate bulk-like water behavior.¹¹ The end set points pertaining to the freezing and melting transitions are provided in Table S2 and the corresponding peak positions in the DSC profiles are provided in Table S3. The classical Gibbs–Thomson (G–T) relation, which describes the depression of melting point in confinement,⁶⁴ is typically expressed as eqn (7).where γ_sl is the solid–liquid interfacial energy, ρ_s is the solid density, ΔH_f is the molar enthalpy of fusion, and D is the effective pore radius.²³ It is worth noting that the derivation of the classical Gibbs–Thomson equation (eqn (7)) relies on the assumption that the volumetric free energy difference between the solid and liquid phases can be approximated aswhere ΔH_f is the molar enthalpy of fusion and T_mb is the bulk melting point. This approximation is valid and allows a linear relation between the melting point depression and inverse pore radius.⁶⁵ To account for the presence of an unfrozen interfacial water layer, a modified form has been proposed by Jähnert et al.,²³ expressed aswhere captures material-specific properties including interfacial energy, enthalpy of fusion, and solid density, and d represents the average thickness of the interfacial layer. Applying the average estimated layer width of 0.39 nm as d, which is calculated based on ϕ_int, we found that the melting point depressions aligned well with predictions from the modified G–T equation (Fig. 5), further validating our estimation of layer widths.

ΔT_m = T_m − T_mb = C_GT/(D − 2d) (10)

Fig. 5 Fitting of melting point depression using the modified Gibbs–Thomson equation (eqn (9)) for water confined in SBA-15 materials with various pore sizes. The fitted curve incorporates a correction term for the interfacial layer width. Data points show experimental results; the fit supports a 0.4 nm interfacial layer. Error bars represent the standard deviation from multiple measurements.

It is also crucial to note at this point that confinement plays a dominant role over the surface constituents as the spatial heterogeneity in the phase transition of water has been observed in several other confinements such as carbon nanotubes and model Lennard-Jones pores.^66–68

By comparing the enthalpy change of confined water during freezing and melting relative to that of bulk water, we observed an overall reduction in enthalpy change within confined systems (Fig. 6). This reduction has been similarly reported for confined water in hydrophilic pores^19,20 and bound water on hydrophilic surfaces,³⁰ suggesting that phase transitions in confined water are unlikely to exhibit the characteristics of a purely first-order transition. However, the presence of intense peaks during these phase transitions also indicates that they are unlikely to resemble glass transitions. Therefore, these phase transitions may be more appropriately classified as weak first-order transitions.^69,70 Notably, as pore size increased, the normalized enthalpy changes also increased, indicating a gradual shift in the nature of phase transitions with varying pore sizes. In particular, for fully filled samples, the normalized enthalpy (relative to bulk water) ranged from approximately 0.27–0.32 for 4 nm pores to 0.33–0.38 for 6 nm pores, and 0.44–0.47 for 8 nm pores, reaching around 0.65 for 12 nm pores. Furthermore, differences in phase transition behavior are apparent across various layers as pore filling increased and the enthalpy change from core water freezing became dominant, the normalized enthalpy change rose significantly from values below 0.1. These results suggest that the phase transitions are comparatively weaker within the water shell layer than in the water core, consistent with the distinct hydrogen bonding structures observed in our simulations as discussed later.

Fig. 6 Normalized enthalpy of confined water during freezing and melting as a function of normalized pore filling. Solid lines represent freezing; dashed lines represent melting. Data are shown for 4, 6, 8, and 12 nm pores. The increasing trend in normalized enthalpy with pore size and filling level indicates a gradual transition toward bulk-like behavior.

3.3. Coordination environment of water and its layered structure in confinement

A general idea of the first coordination shell of water in the nanoconfinement can be obtained from the radial distribution function of the water oxygen atoms (g_OO(r)) in the pore as shown in Fig. 7. Only the first coordination shell is prominent for profiles at room temperature whereas the second coordination shell appears to be more notable at decreasing temperatures. This observation can be attributed to the preferred orientation of the hydrogen bonding network of water at those temperatures (see Section 3.6). The locations of the first and second coordination shells are 2.7 Å and ∼4.3 Å respectively that align well with prior literature reports.^7,71–74 A well-defined minimum between two maxima in the g(r) profiles denotes compact coordination layers of water and minimal probability of molecular exchange between the coordination shells. Moreover, variation of the coordination number (n_OO(r)) with radial distance exhibits interesting features as profiles at different temperatures intersect at a distance of ∼0.32 nm with a coordination number of ∼4, which corresponds to an isosbestic point. This distance of 0.32 nm in fact marks the boundary of the first coordination shell of water oxygen atoms in pores of all widths. This observation suggests that the first coordination shell of water is tetrahedral in these confinements irrespective of the pore width and temperature.^75–77

Fig. 7 Radial distribution function of water oxygen atoms (g_OO(r)) (top panel) and coordination number of the same atoms (n_OO(r)) (bottom panel) confined within the slit pores of width 4, 6, and 8 nm as a function of distance (r) at the control temperatures used in the simulations. The arrows in the top panel indicate the positions of the first and second coordination shells. Isosbestic point located within the 1st coordination shell for all the three systems indicates tetrahedral coordination of water independent of temperature and pore width.

Since the fluid layer in direct contact with the solid surface is known to have different properties and phase transition behavior compared to water in the center of the pore due to short-range interactions,^77,78 the number and mass densities of atoms have been used by multiple studies to distinguish between these two fluid layers.^79–83 To distinguish between the water layers, the variation of mass densities of water along the major axis (Z) of the simulation cell that is normal to the solid surface is shown in Fig. 8. The width of the first peak in the density profile is considered to be equivalent to the width of the combined interface and shell layer of water as observed in the experiments (Fig. 3, 4, Fig. S2, S3 and Table 2). Magnitudes of the average width of combined interface and shell layers computed from simulations are provided in Table 2 (see entry 3) which align well with the experimental findings (see entry 1 and 2 in Table 2). The width is measured from the mean position of the surface hydroxyl groups.⁷⁹ It is important to note that in the presence of TIP5P water used in the simulations, the nature of the water layer within 1 nm of the hydroxylated silica surface exhibits mass density profiles with variable characteristics as a function of temperature. Multiple sub-layers are observed in this region as temperature decreases, forming one or two water layers. This is in line with previous studies^21,79–83 and can originate from the variation in hydrogen bonding pattern of water at lower temperatures (see Section 3.5). Moreover, although the mass density of water is seen to converge to its bulk density at 25 °C in regions beyond 1 nm from the solid surface (Fig. 8), it exhibits a slightly higher density in the combined interface and shell region.^60,81 Interestingly, femtosecond acoustic pulse experiments by Mante and co-workers found that the local density of water is up to five times higher than that of bulk water within 1 nm from the solid surface due to the strong solid–liquid interactions and the roughness of the solid surface.⁶⁰

Fig. 8 Mass density of water (ρ_water,Z) as a function of the distance along z axis of the simulation cell for slit pore of width 4 nm at the control temperatures as used in the simulations. The slit pore being symmetric, the arrows indicate the width of the (interface + shell) layer of confined water in contact with the amorphous silica surfaces. The first and fourth dash-dotted lines in each profile represent the z index of the geometric center of the hydroxyl groups on the silica surface exposed to the water medium that essentially lies on the plane of the hydroxyl groups, whereas the second and third dash-dotted lines indicate the boundaries of the (interface + shell) layer formed by confined water. The average width of the (interface + shell) layer of water in different slit pores is provided in Table 2(3). The profiles for 6 and 8 nm pores are provided in Fig. S5.

3.4. Energetic basis underlying silica–water interactions

To determine the energetic basis underlying the freezing behavior of water in nanoscale confinement, the electrostatic and van der Waals (vdW) energies were determined. Fig. 9 represents the variation of the electrostatic and van der Waals (vdW) energies between silica and water as well as between water molecules themselves in confinements of different sizes normalized by the total number of water molecules in the pores. A clear trend of electrostatic stabilization and vdW destabilization with decreasing temperature is prominent in these profiles. The electrostatic interaction between silica surface and water is most stable in the 4 nm pore at all temperatures which decreases as the pore size is increased. This is quite similar to the findings of Zou and co-workers where the ratio of the total surface potential energy experienced by liquid atoms to their kinetic energy diminishes as the pore size increases following a power law relationship at a specific temperature.⁷⁸ On the other hand, the vdW interaction between silica and water is slightly more repulsive in pores of smaller sizes at all temperatures studied. For interactions between water molecules themselves, the electrostatic interaction attains moderately higher stability as the pore size is increased. However, the magnitudes of electrostatic stability remain nearly invariant in 6 and 8 nm pores. Besides, the vdW interactions between the water molecules are nearly similar in all the slit pores investigated. Density functional theory-based water absorption study by Chen and co-workers using defective α-quartz (0001) surface reveals that there exists a difference of one order of magnitude between water–silica interaction energy and that between hydrogen bonded water molecules where the latter dominates.⁸⁴ This phenomenon also appears to be consistent for water confined within the amorphous silica slit pores examined in this study. Notably, cusps characterized by moderate slope changes observed in the profiles at 103 K indicate a weak first-order phase transition, consistent with experimental observations.^84,85

Fig. 9 Variation of the non-bonding interaction energies (electrostatic (E_{electrostatic}) and van der Waals (E_vdW)) in silica–water and water–water combinations normalized by the number of water molecules in each system as a function of the control temperatures used in the simulations of neat water systems of pore width 4, 6, and 8 nm. The standard deviations are obtained from triplicates of each simulation.

3.5. Hydrogen bonding as a determinant of structure in confined water layers

Hydrogen bonding network, being the most dominant form of dynamic interaction in water, controls nearly all of the physical properties of water such as its spatial structure, density, surface tension, dielectric constant, and colligative properties. For example, in bulk water, each water molecule forms up to four hydrogen bonds with neighboring water molecules that gives rise to a density of 0.9950 g cm⁻³ at 298 K. When cooled to 277 K (4 °C), water achieves its maximum density, and if cooled further, it becomes less dense by expanding its spatial structure.⁸⁶ However, a significant variation in intermolecular hydrogen bonding is noticed in confined water compared to the bulk water. Fig. 10(a) represents the variation of the number of hydrogen bonds per water molecule for combinations such as silica–water and water–water with temperature for all the slit pores. A distance cut-off of 3.5 Å and a donor–H–acceptor angle cut-off of 30° is considered for evaluating the hydrogen bonds.⁸⁷ The hydroxyl groups embedded in the amorphous silica surface along with the non-hydroxyl oxygens are considered in this regard. It is important to note that the hydroxyl O (O₁) and H (H₁) attached to the silica surface have partial charges of −0.95e and 0.425e respectively whereas the O (O₂) and H (H₂) in TIP5P water have partial charges of 0e and 0.241e respectively. The non-hydroxyl oxygen (O₃) in silica has a partial charge of −1.05e. For a straightforward comparison, the H-bonding interaction that is primarily an electrostatic interaction is equivalent to energies (Coulomb's interaction) of ∼−91 kJ mol⁻¹ for O₁–H₂, ∼−100 kJ mol⁻¹ for O₃–H₂, and 0 kJ mol⁻¹ for O₂–H₁ interactions, respectively, at a distance of separation same as the hydrogen bond cut-off. In line with the strong electrostatic attraction between O₁ and H₂, a layer of strongly bound water is evident in ¹H solid-state NMR studies by Grünberg and Vyalikh and co-workers.^88,89 However, the number of O₃–H₂ hydrogen bonds is far less than that of O₁–H₂ and O₂–H₁ in our case because of the nearly complete hydroxylation of the silica surface. Thompson and co-workers found that the number of hydrogen bonds per water molecule decreases to ∼2.2 in hydrophilic porous silica confinement compared to ∼3.6 in bulk.⁹⁰ It is evident from Fig. 10(a) that the number of hydrogen bonds per water molecule is slightly more between water molecules in the core compared to those present in the interface and shell region. Besides, both of them are higher than the silica–water hydrogen bonds. However, in all the combinations, the number of hydrogen bonds increases as temperature is reduced. The maximum number of hydrogen bonds per water molecule seen in the case of silica–water is 0.29 whereas they are 1.10 and 1.45 in the case of water–water combination in (interface + shell) and core region, respectively. A clear trend of increase in the number of hydrogen bonds between silica and water is noticed at all temperatures with decreasing pore sizes. However, this behavior is similar in the 6 nm and 8 nm pores near the region of the lowest temperature. This observation is in contrast to the water–water hydrogen bonds, which decrease with reduction in the pore diameter while water confined in 6 and 8 nm pores behave fairly similarly to each other. This result is in good agreement with the equilibrium MD simulations reported by Zhang and co-workers at 298.15 K where ∼0.8 hydrogen bonds per water molecule is observed to form between water molecules at a distance of ∼3.7 Å from the pore surface of a 4.75 Å cylindrical silica nanopore.⁹¹

Fig. 10 (a) Variation of the average number of hydrogen bonds per water molecule with respect to different temperatures as used in the simulations for 4, 6, and 8 nm slit pores between ‘silica and water’ and water molecules themselves in (interface + shell) region as well as the core region of the pores. The standard deviations are obtained from simulation triplicates. A distance cut-off of 3.5 Å and donor–H–acceptor angle cut-off of 30° are used for calculating the hydrogen bonds. All Oxygens and hydrogens present on the silica surface are considered in this regard; (b) snapshots representing the spatial orientation of the water molecules present within 4 Å from a small representative section of the silica surface in 4, 6, and 8 nm slit pores at 103 K representing majority of the water molecules acquiring H-down orientation with respect to silica. The silica surface is represented as vdW spheres of Si (yellow), O (red), and H (white) whereas water molecules are shown in CPK with O and H in red and white respectively; (c) different shapes of hydrogen bonded water network as seen in the (interface + shell) and core region of 4, 6, and 8 nm slit pores at temperatures 298 K and 103 K. The structures are viewed along the xy plane where they vary from simple tetrahedral network to regular/distorted pentagon, regular/distorted hexagon and higher order geometries. However, the majority of the water molecules display amorphous hydrogen bonded network in all the cases.

3.6. Structural distortion of ice in nanopores

The solid–liquid phase diagram of water exhibits multiple stable and metastable phases of ice. Liquid water shares a phase boundary with normal or hexagonal ice I_h, and ices III, IV, VI, and VII. Ice I_h is stable within a wide temperature range of −200 °C (73 K) to 0 °C (273 K) and a short range of pressure. On the other hand, ice VIII is stable within the same temperature range, but a wide range of pressure (13 to 25 kbar). Therefore, it is evident that different phases of ice can be observed in nanoconfinement depending on the pore pressure.^78,86

It is widely known that pressure exhibits wide variability with respect to that of the bulk water in the confined environments. It has been seen that competing effects of solid–liquid and liquid–liquid interactions determine the pressure in the nanopore and consequently, the liquid layer in direct contact with the surface has a significant influence on the pore pressure. Depending on the pore size, the pressure can be either positive or negative. A surface influence number as developed by Zou and co-workers predicts that the effect of solid–liquid interactions originating from the surface diminishes following a power law as the pore size is increased at a particular temperature. The change in pressure in the confinement can distort the H–O–H angle in water giving rise to distortion in the ice lattice that can include both orientational and ionic defects.^78,86

Strong electrostatic attraction between the water dipoles and the silica surface is known to force the water molecules in hydrogen down orientation⁹² as evidenced by the predominance of hydrogen-down water molecules in the simulation snapshots (Fig. 10(b)). Moreover, the widely spread hydrogen bond network present in water gives rise to water clusters of different geometries ranging from simple tetrahedral network of water molecules, regular, or distorted pentagons and hexagons to higher order clusters consisting of seven or eight water molecules.⁹³Fig. 10(c) exhibits some of these clusters as observed during the simulations. Thorough inspection of the water layers at different temperatures confirms the fact that amorphous solid water (ASW), a glassy state of water is formed during the simulation. Prior literature reports suggest that a cooling rate of the order of 10⁶ K s⁻¹ is sufficiently rapid to form the non-crystalline solid water (i.e. glassy water) by completely avoiding the freezing phenomenon.^27,94 The fast kinetics does not allow the water molecules enough time to arrange themselves in a crystalline state at the supercooling temperatures and consequently, the glassy state is reached. Low density amorphous (LDA), high density amorphous (HDA), and very high density amorphous (VHDA) ice are the three major forms of glassy water and their formation and interconversion depend on the experimental or simulation protocols. The cooling rate used in our simulations is of the order of 10¹⁰ K s⁻¹ which is fast enough to allow the glass transition of water leading to the formation of ASW in combination with effects from pore pressure and surface–water interactions. Fig. 11 shows the ASW structures observed in the (interface + shell) and core region of the slit pores at the lowest temperature sampled. Coupled with the density profiles and number of hydrogen bonds shown in Fig. 8 and 10(a) respectively, the ASW structures reveal that HDA water is found in the (interface + shell) region whereas the core is populated by the LDA water. This in turn, also suggests that the total pressure in the core is relatively less than that in the interface and shell.

Fig. 11 Representative hydrogen bonding network of amorphous solid water (ASW) in (interface + shell) and core regions of the 4, 6, and 8 nm slit pores. The water clusters are shown in the xy and xz plane respectively for (interface + shell) and core. Water molecules are represented in CPK with O and H in red and white repectively. Hydrogen bonds are represented as red dahsed lines evaluated with a distance cut-off of 3.5 Å and donor–H–acceptor angle cut-off of 30°. ASW present in the combined interface and shell is of type HDA and that present in the core is of the type LDA as evident from the higher number of hydrogen bonds and lower density of water in the core compared to those in (interface + shell) region.

4. Conclusion

In this study, a comprehensive methodology for determining the mechanisms underlying the freezing and melting behavior of water confined in silica nanopores with pore diameters in the range of 4 nm to 12 nm, is established. Layered structures of water are observed on freezing in silica nanopores and these structures and the corresponding freezing points are determined using differential scanning calorimetry (DSC) measurements. Multiple layers of water with varying freezing points are noted when the pore filling is typically greater than 45%. Specifically, a stable three-layer structure comprising an interfacial layer, a shell, and a core was identified for water confined in pores with diameters 4 and 8 nm. From the pore filling levels associated with each peak formation, we estimated the interfacial layer width at approximately 0.39 nm and the combined interfacial and shell layers at around 0.62 nm. For the 12 nm pores, a dynamic structural adaptation is observed, where the water core merged with the shell to form a unified layer at higher pore filling levels. The melting point depression depends on pore size across all samples, which we described using the modified Gibbs–Thomson equation with an estimated unfrozen layer of 0.39 nm based on the DSC data.

In the case of water confined in silica pores with diameters of 4–8 nm, the additional melting point depression correlated with pore filling levels, likely resulting from enhanced solid–liquid interactions in the smaller, hydrophilic confinements compared to the 12 nm case. Enthalpy analyses indicated a trend of increasing enthalpy with both pore size and filling level, with lower enthalpy changes consistently noted for water near the pore surface and in smaller nanopores. Molecular dynamics (MD) simulations further elucidated these interactions, using amorphous silica slit pores with diameters ranging from 4 to 8 nm to examine the coordination behavior of water and its layered structure during freezing. The primary density peak, with a width relatively stable across pore sizes and temperatures, aligned with the combined interfacial and shell layers observed experimentally. Electrostatic interactions between the silica surface and confined water molecules predominantly stabilized the water structure within smaller pores, while van der Waals interactions slightly destabilized the system, especially at lower temperatures and in narrower confinements. The non-bonding interactions between the silica surface and interfacial water significantly influence the hydrogen bonded network of water, serving as the primary factor behind the formation of the unfrozen interfacial water layer. Distorted ice polymorphs in a matrix of amorphous solid water are also identified, with their formation driven by the cooling rate, pore diameter, and silica–water interactions. Based on observations from experimental and simulation data, the phase transition of water in SBA-15 with pore diameters of 4–8 nm is characterized as weak first-order transitions that varied depending on the pore diameter and the number of water layers within the pores. The insights obtained on phase transition and energetics would guide further research on various geological processes, energy storage, and controlled drug delivery.

Conflicts of interest

There are no conflicts of interest to declare for this work.

Data availability

The data supporting this article have been included in the main manuscript and as part of the SI. See DOI: https://doi.org/10.1039/d5cp02430g

Acknowledgements

This work is supported by the Army Research Office (Proposal Number: 79587-BB-ECP). The authors gratefully acknowledge the use of the shared supercomputing resources, G-2 Clusters, at Cornell Ann S. Bowers College of Computing and Information Science. The contributions of A. W. and S. M. P. as undergraduate students are supported by the Army Outreach Education Program.

References

1. P. Bampoulis, K. Sotthewes, E. Dollekamp and B. Poelsema, Surf. Sci. Rep., 2018, 73, 233–264.
2. M. Tanaka, S. Morita and T. Hayashi, Colloids Surf., B, 2021, 198, 111449.
3. M. Jaugstetter, N. Blanc, M. Kratz and K. Tschulik, Chem. Soc. Rev., 2022, 51, 2491–2543.
4. S. A. Bagshaw, E. Prouzet and T. J. Pinnavaia, Science, 1995, 269, 1242–1244.
5. D. Zhao, J. Feng, Q. Huo, N. Melosh, G. H. Fredrickson, B. F. Chmelka and G. D. Stucky, Science, 1998, 279, 548–552.
6. Y. Xia, H. Cho, M. Deo, S. H. Risbud, M. H. Bartl and S. Sen, Sci. Rep., 2020, 10, 5327.
7. S. Mohammed, H. Asgar, C. J. Benmore and G. Gadikota, Phys. Chem. Chem. Phys., 2021, 23, 12706–12717.
8. G. H. Findenegg, S. Jähnert, D. Akcakayiran and A. Schreiber, Chem. Phys. Chem., 2008, 9, 2651–2659.
9. C. Lederle, M. Sattig and M. Vogel, J. Phys. Chem. C, 2018, 122, 15427–15434.
10. S. Ruiz-Barragan, H. Forbert and D. Marx, Phys. Chem. Chem. Phys., 2023, 25, 28119–28129.
11. A. W. Knight, N. G. Kalugin, E. Coker and A. G. Ilgen, Sci. Rep., 2019, 9, 8246.
12. P. A. Bonnaud, H. Manzano, R. Miura, A. Suzuki, N. Miyamoto, N. Hatakeyama and A. Miyamoto, J. Phys. Chem. C, 2016, 120, 11465–11480.
13. S. Layek and N. Sengupta, ChemPhysChem, 2024, 25, e202400618.
14. A. Zaragoza, M. A. González, L. Joly, I. López-Montero, M. A. Canales, A. L. Benavides and C. Valeriani, Phys. Chem. Chem. Phys., 2019, 21, 13653–13667.
15. J. C. Johnston and V. Molinero, J. Am. Chem. Soc., 2012, 134, 6650–6659.
16. F. Mallamace, M. Broccio, C. Corsaro, A. Faraone, D. Majolino, V. Venuti, L. Liu, C.-Y. Mou and S.-H. Chen, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 424–428.
17. M. Sattig, S. Reutter, F. Fujara, M. Werner, G. Buntkowsky and M. Vogel, Phys. Chem. Chem. Phys., 2014, 16, 19229–19240.
18. E. W. Hansen, H. C. Gran and E. J. Sellevold, J. Phys. Chem. B, 1997, 101, 7027–7032.
19. C. Faivre, D. Bellet and G. Dolino, Eur. Phys. J. B, 1999, 7, 19–36.
20. A. Schreiber, I. Ketelsen and G. H. Findenegg, Phys. Chem. Chem. Phys., 2001, 3, 1185–1195.
21. K. Ishikiriyama and M. Todoki, Thermochim. Acta, 1995, 256, 213–226.
22. I. C. Bourg and C. I. Steefel, J. Phys. Chem. C, 2012, 116, 11556–11564.
23. S. Jähnert, F. Vaca Chávez, G. E. Schaumann, A. Schreiber, M. Schönhoff and G. H. Findenegg, Phys. Chem. Chem. Phys., 2008, 10, 6039.
24. M. Brodrecht, E. Klotz, C. Lederle, H. Breitzke, B. Stühn, M. Vogel and G. Buntkowsky, Z. Phys. Chem., 2018, 232, 1003–1015.
25. S. Han, M. Y. Choi, P. Kumar and H. E. Stanley, Nat. Phys., 2010, 6, 685–689.
26. M. Erko, G. H. Findenegg, N. Cade, A. G. Michette and O. Paris, Phys. Rev. B:Condens. Matter Mater. Phys., 2011, 84, 104205.
27. P. G. Debenedetti and H. E. Stanley, Phys. Today, 2003, 56, 40–46.
28. S. Agrawal, M. Elmehlawy and M. P. Hoepfner, J. Phys. Chem. C, 2021, 125, 11097–11106.
29. G. Kupgan, T. P. Liyana-Arachchi and C. M. Colina, Langmuir, 2017, 33, 11138–11145.
30. A. Tahmasebi, J. Yu, H. Su, Y. Han, J. Lucas, H. Zheng and T. Wall, Fuel, 2014, 135, 243–252.
31. H. Kumano, T. Asaoka, A. Saito and S. Okawa, Int. J. Refrig., 2007, 30, 267–273.
32. S. Li, L. C. Dickinson and P. Chinachoti, J. Agric. Food Chem., 1998, 46, 62–71.
33. R. Wyckoff, Z. Kristallogr. - Cryst. Mater., 1925, 62, 189–200.
34. J. M. Stallons and E. Iglesia, Chem. Eng. Sci., 2001, 56, 4205–4216.
35. C. Pirez, J.-C. Morin, J. C. Manayil, A. F. Lee and K. Wilson, Microporous Mesoporous Mater., 2018, 271, 196–202.
36. L. M. Yang, Y. J. Wang, G. S. Luo and Y. Y. Dai, Microporous Mesoporous Mater., 2005, 84, 275–282.
37. J. Geske and M. Vogel, Mol. Simul., 2017, 43, 13–18.
38. R. T. Cygan, J.-J. Liang and A. G. Kalinichev, J. Phys. Chem. B, 2004, 108, 1255–1266.
39. G. A. Kaminski, R. A. Friesner, J. Tirado-Rives and W. L. Jorgensen, J. Phys. Chem. B, 2001, 105, 6474–6487.
40. M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys., 2000, 112, 8910–8922.
41. P. Hirunsit and P. B. Balbuena, J. Phys. Chem. C, 2007, 111, 1709–1715.
42. N. Pingua and A. K. Gautam, Mater. Today Commun., 2024, 39, 108777.
43. Z. Wang, Y. Zhou, Y. Feng, J. Zhang, R. Yu and Z. Yu, J. Am. Ceram. Soc., 2024, jace.20047.
44. A. Eltareb, G. E. Lopez and N. Giovambattista, Commun. Chem., 2024, 7, 36.
45. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
46. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
47. ed. R. A. de Groot and J. Nadrchal, Physics Computing 92: Proceedings of the 4th International Conference, World Scientific, Singapore, 1993.
48. I. R. Piletic, D. E. Moilanen, D. B. Spry, N. E. Levinger and M. D. Fayer, J. Phys. Chem. A, 2006, 110, 4985–4999.
49. I. Brovchenko and A. Oleinikova, Interfacial and confined water, Elsevier, Amsterdam Oxford, 2008.
50. S. Varghese, S. K. Kannam, J. S. Hansen and S. P. Sathian, Langmuir, 2019, 35, 8159–8166.
51. I. Hanasaki and A. Nakatani, J. Chem. Phys., 2006, 124, 174714.
52. S. Kittaka, S. Takahara, H. Matsumoto, Y. Wada, T. J. Satoh and T. Yamaguchi, J. Chem. Phys., 2013, 138, 204714.
53. R. Schmidt, E. W. Hansen, M. Stoecker, D. Akporiaye and O. H. Ellestad, J. Am. Chem. Soc., 1995, 117, 4049–4056.
54. F. M. Etzler and D. M. Fagundus, J. Colloid Interface Sci., 1987, 115, 513–519.
55. M. D. Fayer and N. E. Levinger, Annu. Rev. Anal. Chem., 2010, 3, 89–107.
56. X. X. Liu, Q. Wang, X. F. Huang, S. H. Yang, C. X. Li, X. J. Niu, Q. F. Shi, G. Sun and K. Q. Lu, J. Phys. Chem. B, 2010, 114, 4145–4150.
57. D. Feng, K. Wu, S. Bakhshian, S. A. Hosseini, J. Li and X. Li, Langmuir, 2020, 36, 8764–8776.
58. T. Iiyama, K. Nishikawa, T. Otowa and K. Kaneko, J. Phys. Chem., 1995, 99, 10075–10076.
59. R. Horstmann, L. Hecht, S. Kloth and M. Vogel, Langmuir, 2022, 38, 6506–6522.
60. P.-A. Mante, C.-C. Chen, Y.-C. Wen, H.-Y. Chen, S.-C. Yang, Y.-R. Huang, I. Ju Chen, Y. W. Chen, V. Gusev, M. J. Chen, J. L. Kuo, J. K. Sheu and C. K. Sun, Sci. Rep., 2014, 4, 6249.
61. K. Morishige, J. Phys. Chem. C, 2018, 122, 5013–5019.
62. T. J. Rottreau, C. M. A. Parlett, A. F. Lee and R. Evans, Microporous Mesoporous Mater., 2018, 264, 265–271.
63. Y. Xiong, P. S. Shabane and A. V. Onufriev, ACS Omega, 2020, 5, 25087–25094.
64. C. L. Jackson and G. B. McKenna, J. Chem. Phys., 1990, 93, 9002–9011.
65. L. Makkonen, Langmuir, 2002, 18, 1445–1448.
66. A. Tinti, G. Camisasca and A. Giacomello, Philos. Trans. R. Soc., A, 2021, 379, 20200403.
67. Y. Liu, J. Jiang, Y. Pu, J. S. Francisco and X. C. Zeng, ACS Nano, 2023, 17, 6922–6931.
68. P. Pugliese, M. M. Conde, M. Rovere and P. Gallo, J. Phys. Chem. B, 2017, 121, 10371–10381.
69. B. Van Roie, J. Leys, K. Denolf, C. Glorieux, G. Pitsi and J. Thoen, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2005, 72, 041702.
70. P. Losada-Pérez, J. Leys, G. Cordoyiannis, C. Glorieux and J. Thoen, in Enthalpy and Internal Energy, ed. E. Wilhelm and T. Letcher, Royal Society of Chemistry, Cambridge, 2017, pp. 364–379.
71. A. Fouzri, R. Dorbez-Sridi and M. Oumezzine, J. Chem. Phys., 2002, 116, 791–797.
72. P. Smirnov, T. Yamaguchi, S. Kittaka, S. Takahara and Y. Kuroda, J. Phys. Chem. B, 2000, 104, 5498–5504.
73. K. Yoshida, T. Yamaguchi, S. Kittaka, M.-C. Bellissent-Funel and P. Fouquet, J. Chem. Phys., 2008, 129, 054702.
74. A. Fouzri, R. Dorbez-Sridi and M. Oumezzine, Eur. Phys. J.:Appl. Phys., 2003, 22, 21–28.
75. C. Benmore, L. C. Gallington and E. Soignard, Mol. Phys., 2019, 117, 2470–2476.
76. K. Koga, H. Tanaka and X. C. Zeng, Nature, 2000, 408, 564–567.
77. A. G. Kalinichev, J. Wang and R. J. Kirkpatrick, Cem. Concr. Res., 2007, 37, 337–347.
78. A. Zou and S. C. Maroo, Phys. Fluids, 2021, 33, 042007.
79. D. Argyris, T. Ho, D. R. Cole and A. Striolo, J. Phys. Chem. C, 2011, 115, 2038–2046.
80. D. Argyris, D. R. Cole and A. Striolo, J. Phys. Chem. C, 2009, 113, 19591–19600.
81. A. A. Milischuk and B. M. Ladanyi, J. Chem. Phys., 2011, 135, 174709.
82. M. Sulpizi, M.-P. Gaigeot and M. Sprik, J. Chem. Theory Comput., 2012, 8, 1037–1047.
83. C. Weinberger, F. Zysk, M. Hartmann, N. K. Kaliannan, W. Keil, T. D. Kühne and M. Tiemann, Adv. Mater. Interfaces, 2022, 9, 2200245.
84. Y.-W. Chen and H.-P. Cheng, J. Chem. Phys., 2011, 134, 114703.
85. V. F. Correa and F. J. Castro, Entropy, 2021, 24, 31.
86. The Physics and Physical Chemistry of Water, ed. F. Franks, Springer, New York, 2012.
87. S. Han, P. Kumar and H. E. Stanley, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2009, 79, 041202.
88. B. Grünberg, T. Emmler, E. Gedat, I. Shenderovich, G. H. Findenegg, H. Limbach and G. Buntkowsky, Chem. – Eur. J., 2004, 10, 5689–5696.
89. A. Vyalikh, Th Emmler, B. Grünberg, Y. Xu, I. Shenderovich, G. H. Findenegg, H.-H. Limbach and G. Buntkowsky, Z. Phys. Chem., 2007, 221, 155–168.
90. H. Thompson, A. K. Soper, M. A. Ricci, F. Bruni and N. T. Skipper, J. Phys. Chem. B, 2007, 111, 5610–5620.
91. Q. Zhang, K.-Y. Chan and N. Quirke, Mol. Simul., 2009, 35, 1215–1223.
92. P. A. Bonnaud, B. Coasne and R. J.-M. Pellenq, J. Phys.: Condens. Matter, 2010, 22, 284110.
93. C. G. Salzmann, P. G. Radaelli, B. Slater and J. L. Finney, Phys. Chem. Chem. Phys., 2011, 13, 18468.
94. P. Jenniskens and D. F. Blake, Science, 1994, 265, 753–756.

---