Solvent Exchange Leading to Nanobubble Nucleation: A Molecular Dynamics Study.

The solvent exchange procedure has become the most-used protocol to produce surface nanobubbles, while the molecular mechanisms behind the solvent exchange are far from being fully understood. In this paper, we build a simple model and use molecular dynamics simulations to investigate the dynamic characteristics of solvent exchange for producing nanobubbles. We find that at the first stage of solvent exchange, there exists an interface between interchanging solvents of different gas solubility. This interface moves toward the substrate gradually as the exchange process proceeds. Our simulations reveal directed diffusion of gas molecules against the gas concentration gradient, driven by the solubility gradient of the liquid composition across the moving solvent-solvent interface. It is this directed diffusion that causes gas retention and produces a local gas oversaturation much higher near the substrate than far from it. At the second stage of solvent exchange, the high local gas oversaturation leads to bubble nucleation either on the solid surface or in the bulk solution, which is found to depend on the substrate hydrophobicity and the degree of local gas oversaturation. Our findings suggest that solvent exchange could be developed into a standard procedure to produce oversaturation and used to a variety of nucleation applications other than generating nanobubbles.

Introduction

Numerous studies confirmed the existence of interfacial nanobubbles[1−7] that preferentially nucleate on hydrophobic solid surfaces immersed in solutions with dissolved gas.[8−13] Their stability has been traced back to contact line pinning originating from the chemical and physical heterogeneity of the substrate.[2,14−18] This view is confirmed in numerical simulations.[19−23] In the pinned state, the gas pressure in the solution and the Laplace pressure in the bubble can stably balance.[2,17,19,21] Because nanobubbles have a variety of potential applications, such as boundary slip in fluids, flotation of minerals, bimolecular adsorption, and immersion lithography,[24−28] their formation and unusual properties have drawn much attention for intensive investigations.

For the formation of nanobubbles, the most-used protocol is the solvent exchange procedure.[11,29−32] As the first example, Lou et al.[1] used ethanol to clean the solid surface, and then water was injected to exchange the ethanol, after which nanobubbles were clearly observed on the substrate. The standard protocol of the solvent exchange process normally includes three stages. A hydrophobic substrate is first contacted with water which is then replaced by ethanol. At this stage, no nanobubbles can be observed. But after the ethanol is replaced by water, nanoscale bubbles covering the substrate surface can be found. Besides ethanol, other organic solvents,[11] such as methanol and 2-propanol, can also be used in the solvent exchange process. Similarly, other exchange processes such as exchanging cold water against warm water[33] and ethanol solution against salt solution[34] have been successfully used to produce nanobubbles. Moreover, the solvent exchange process can also be applied to produce interfacial nanodroplets[2,35−39] at an interface between a solid and an immiscible liquid.

It has been assumed that during the exchange process, when the ethanol is replaced gradually by water, gas molecules cannot diffuse into the atmosphere but stay in the water. As air has a higher solubility in ethanol than in water, the exchange process thus first leads to gas oversaturation and then to nanobubble nucleation. Though the solvent exchange process has become the most used protocols to produce nanobubbles in experiments, the exact mechanism to produce the gas saturation and then the nanobubbles is yet unknown. Furthermore, this method cannot precisely control nanobubble formation, because a variety of factors such as the saturation level of gas, the exchange rate, or the liquid shear and flow boundary conditions cannot be perfectly controlled in experiments. Therefore, it is hard to experimentally study how a solvent exchange process dynamically nucleates nanobubbles. In this work, alternatively, we use molecular dynamic (MD) simulations to study the dynamical characteristics of solvent exchange and the corresponding mechanisms of generating nanobubbles and to identify the regions in parameter space where surface nanobubbles are formed with the solvent exchange process. The advantage of the MD simulations is that the effect of various factors can be precisely controlled and separated. The obvious disadvantage is that the length and time scales of simulations remain very limited.

Simulation Model and Method

To investigate the molecular mechanisms of the solvent exchange process for nanobubble formation, MD simulations were performed by using LAMMPS,[40] an open source program for massively parallel simulations. For our simulations, we establish a model system to simulate the solvent exchange process, as shown in Figure a. First, a quasi-two-dimensional simulation box with a size of 22.4 × 2.24 × H nm3 was built with the height H of the simulation box fluctuating at a given pressure. Periodic boundary conditions were used in the x and y directions, while in the z direction two solid substrates that consist of frozen solid molecules on a FCC lattice with a lattice parameter of 5.606 Å and the (100) surface were used to restrain the fluid. The bottom substrate was fixed during the simulations, and a square pore with a width of 13.64 nm and a depth of 3.36 nm was introduced on the substrate to pin the contact line of the generated nanobubbles.

Figure 1

(a) Early (t = 1 ns) and (b) late (t = 600 ns) representation of the simulation box of 22.4 × 2.24 × H nm3. The green particles represent solvent II having a poor gas solubility (e.g., water), the gray ones represent solvent I having a high gas solubility (e.g., ethanol), the red ones represent the gas molecules, the blue ones represent the solid particle of the top substrate, and the orange ones represent the solid particles of the bottom substrate. The shaded area shows the source region that controls the gas concentration in the reservoir of solvent II. Regions marked by 1, 2, 3, and 4 represent bulk liquid near the source region, and bulk liquid far from both the source region and the substrate, near-substrate region, and inside-pore region, respectively. (c) The calculated Young contact angle for a droplet of solvent II as a function of the interaction strength ϵSL between molecules of solvent II and the bottom substrate at T = 81.2 K.

A minimal model for the solvent exchange between a good solvent (solvent I) and a bad solvent (solvent II) should catch the following essential factors: The solvent exchange itself, diffusion of dissolved gas, and contact with an infinite gas reservoir in solvent II. To mimic the solvent exchange process in contact with the bulk solution of solvent II with a dissolved gas, a source region representing the infinite reservoir was included in the simulation box to control the gas concentration far from the substrates (see the shaded region in Figure a). Solvent exchange was then allowed by molecular diffussion. In practice, after a given time interval of regular MD simulations, the identity exchange between liquid and gas molecules in the source region was periodically carried out to keep the gas concentration in the reservoir of solvent II fixed. Once the gas concentration in the reservoir became smaller than the target value as a result of more gas molecules leaving this reservoir, a number of solvent II molecules in the region were randomly chosen and their identify was changed to that of gas molecules. Similarly, if more gas molecules from outside diffused into the region than left, a number of randomly chosen gas molecules in the region would be changed into molecules of solvent II. At the same time interval, the solvent exchange in the region was also conducted in a similar way via interchanging the particle identities between solvent I with a relatively high gas solubility(xsolubilityI = 0.14) and solvent II of smaller gas solubility xsolubilityII = 0.00031, maintaining the ratio of solvent II to solvent I. Note that in this work we assume that there is no solvent I molecules in the bulk solution of solvent II. In this way, an explicit reservoir with fixed gas concentration and solvent composition, which here represents the bulk solvent II with dissolved gas, was included in our model to provide a concentration gradient as in real processes. The interfacial region between the source region and the solid surface is the diffusion zone, through which solvents and gas are diffusively transported under the constraint of the chemical potential gradient imposed by the source region. Thus, with conditions close to those in a real solvent exchange experiments, in which gradients in the chemical potential drive the flow, this model enables us to simulate the phenomenon of solvent exchange and diffusion. Rough estimate from a typical solvent exchange process[38] shows a flow velocity of ∼10–6 nm/ns–1 at 30 nm above the substrates. Because we only investigate the dynamic mechanism of solvent exchange for nucleating nanobubbles near the solid surface where the fluid velocity is negligible, we did not consider the effect of shear force in this work, as for gaseous liquids without solvent exchange.[41,42] Because we only considered a region of ∼30 nm above the bottom substrate, here we changed xsourcegas from 0.001 to 0.005, a little larger than xsolubilityII, leading to an oversaturation of ζ > 0 after the solvent exchange process but much smaller than xsolubilityI.

For the intermolecular interactions, a truncated Lennard–Jones (LJ) 12-6 potential was employed with a cutoff distance of 1.1 nm (see Table for LJ parameters). For solvent II, in particular, its interaction with the bottom substrate was varied to model different substrate hydrophobicity, and the given interaction parameter here would produce a contact angle for a sitting droplet (θ) within the range from 31° to 130° (see Figure b). Although reduced units were used in our simulations, all variables were reported here in their actual physical units. To convert reduced units to their real units, both mass m and LJ parameters were chosen as those of the argon atom.

Table 1

Parameters for Lennard–Jones Interaction between Different Molecules

molecules	σ (nm)	ϵ (meV)
solvent I–solvent I	0.341	10.30
solvent I–solvent II	0.341	10.30
solvent I–solid (bottom)	0.341	8.24
solvent I–solid (top)	0.341	10.30
solvent I–gas	0.341	8.24
solvent II–solvent II	0.341	10.30
solvent II–solid (bottom)	0.341	4.12 to 8.24
solvent II–solid (top)	0.341	10.30
solvent II–gas	0.341	3.09
gas–solid (bottom)	0.341	1.89
gas–solid (top)	0.341	1.72
gas–gas	0.341	3.43

We carried out isothermal, isostress (NPT) ensemble MD simulations with a fixed number of fluid molecules N = 35520, T = 81.2 K, and P = 5 atm. An external force along the z direction was exerted on the smooth top substrate to maintain the given pressure. The integration of equations of motion was the classical velocity Verlet algorithm with a time step of 5 fs. The fluid temperature was controlled by the Nosé–Hoover method with a time constant of 0.5 ps.[43] Note that we also performed MD simulations at another pressure of P = 1 atm, and because the same solvent exchange mechanism was obtained, here we only reported the results for P = 5 atm.

Results and Discussion

Directed Diffusion of Gas Molecules against Its Concentration Gradient for Generating Local Oversaturation

To understand the dynamics of solvent exchange and how it leads to nanobubble nucleation, we fixed the gas concentration in the source region xsourcegas = 0.004 and varied θ from 131° to 31°. We recorded the time evolution of the local gas density at different locations that include regions 1 to 4 (see their locations denoted in Figure a) as well as several typical snapshots (see Figure ). The figure reveals that the whole process of solvent exchange can be divided into two stages. During the first stage, the gas concentration in the bulk solution increased with time and with the distance from the source region. Using θ = 130° as an example, we can find from Figure d that although near the source region the gas concentration Cgas (in units of mol/L) seems to remain unchanged for the whole process, Cgas in the bulk liquid near the substrate or inside the pore (see regions 2, 3, and 4 in Figure d) increases continuously within 0–200 ns. Unexpectedly and remarkably, the first stage finally produced a local gas concentration in those regions much higher than Csourcegas. The first stage ended when the first small bubble appeared at the hydrophobic substrate at ∼250 ns (Figure a), after which Cgas in the bulk solution began to decrease.

Figure 2

Nanobubbles are nucleated as a result of the temporal–spatial evolution of solvent composition and gas concentration. (a–c) Typical snapshots at different simulation time and (d–f) time evolution of the local gas density in different regions (see regions 1, 2, 3, and 4 marked in Figure a) during the solvent exchange processes at xsourcegas = 0.004 and different values of θ: (a, d) θ = 130°; (b, e) 91°; (c, f) 31°. The color code of the snapshots is the same as denoted in Figure . Note that although (a, d) and (c, f) lead to basically the same final bubble morphologies, their nucleation pathways are totally different because of different substrate hydrophobicity: One nucleated from the substrate and the other from the bulk solution.

The second stage features the growth of the nucleated nanobubble toward its stable state, and meanwhile the gas density in the bulk solution decreases to its value in the control zone, Csourcegas (see Figure a,d). The decrease of Cgas is because when the bubble formed, surrounding gas molecules tend to diffuse into the bubble, which thus reduces the chemical potential and the density of gas molecules in the surrounding environment. Note that at the second stage, the gas concentration in regions 2, 3, and 4 may increase until they reach a very high value. Inspection of the corresponding snapshots shows that the increase is induced by the first nucleus located in those regions (region 3 in Figure d, region 4 in Figure e, regions 2 and 3 in Figure f) or the gas enrichment in the corner of hydrophobic pore (region 4 in Figure d). This is obviously different from the increase of gas density at the first stage, which is due to the gas oversaturation. Note that the produced nanobubbles in Figure have not yet reached their equilibrium shape that is solely determined by the degree of gas oversaturation and independent of the hydrophobicity of the surface.[17]

To interpret why a much higher Cgas than Csourcegas is produced for locations far from the source region at the first stage, we show typical density profiles for solvent I and gas molecules in Figure . The figure indicates that during the first stage there exists an interface (either sharp or a little blurred) between solvent I and solvent II, and the interface moves toward the bottom substrate gradually as the exchange process proceeds (see Figure ). When the moving interface either reached the bottom substrate or a bubble was nucleated, the first stage ended.

Figure 3

Evolution of solvent composition and local gas concentration during solvent exchange processes. Density profiles for solvent I (a, c) and gas molecules (b, d) for several typical snapshots during solvent exchange processes of xsourcegas = 0.004. In the pictures we changed the substrate hydrophobicity and set θ = 91° in (a, b) and θ = 31° in (c, d).

The most surprising observation at the first stage is the existence of a directed diffusion of gas molecules toward the substrate against its concentration gradient, which causes gas retention. This contradicts naive expectations concerning the gas diffusion. However, the existence of a moving interface between solvent I and solvent II (Figure ) indicates that the gas diffusion against its density gradient is in fact due to the increase of gas solubility toward the region of solvent I at the interface. In other words, during the exchange process, the liquid near the moving solvent–solvent interface has different solvent compositions, and thus the solubility of the gas is solvent-composition-dependent. Therefore, there actually can be a flux of diffused gas from lower concentration to higher concentration, because the latter corresponds to a region with the solvent composition having a high solubility of gas molecules. In other words, the diffusion of gas molecules tends to point to the solvent region that has a higher ratio of the solvent with a higher solubility, although the gradient of gas concentration is against this direction. The forced diffusion prevents the gas from being washed away and causes the enrichment of gas molecules in the solvent I-rich region. Thus, a local oversaturation degree of gas molecules appears near the interface, and it increases with time and as the interface moves toward the bottom substrate. The nucleation of nanobubbles occurs once the local gas oversaturation near the substrate is sufficiently high.

Nucleation Pathway Depends on Local Gas Oversaturation and Solid Hydrophobicity

Next we performed extensive simulations to investigate how the nucleation pathway changes with gas concentrations in the source region xsourcegas and the hydrophobicity (expressed through the contact angle θ) of the bottom substrate. Here we considered different concentrations of gas molecules in the source region, including xsourcegas = 0.001, 0.002, 0.003, 0.004, and 0.005, which are larger than the gas solubility in solvent II (0.00031) but much smaller than that for solvent I (0.14). We also considered the effect of substrate hydrophobicity via changing the interaction between solid and solvent II that produces the Young contact angle θ from 131° to 31°. The obtained nucleation pathways from the extensive simulations are summarized in Figure . The figure clearly indicates that both Csourcegas and θ play key roles in nucleating nanobubbles.

Figure 4

Phase diagram for pathways of nucleating nanobubbles as a function of xsourcegas and substrate hydrophobicity θ. In this figure, bulk nucleation indicates formation of nanobubbles that is initiated in the bulk solution, while interface nucleation indicates nucleation that occurs at the substrate.

In general, interface nucleation of nanobubbles, for which the nucleus is initially formed on the top of the bottom substrate or inside the substrate pore, can be observed at high xsourcegas and large θ. This is because for substrates with a large θ, the hydrophobic nature of the solid surface would result in the enhancement of gas molecules near the substrate and thus promote the surface nucleation. The same effect occurs for the case with high xsourcegas: the high xsourcegas produces a higher local gas oversaturation near the substrate. In contrast, for hydrophilic substrates and small dissolved gas content (e.g., xsourcegas < 0.003), there is no bubble formation observed, and instead a Wenzel state was found, as the solvent II would finally wet the pore. This observation is in agreement with experimental observation that hydrophilic surfaces and degassing inhibit nanobubble formation.[2,44−46] When xsourcegas is sufficiently high while the bottom substrate is hydrophilic, however, the nucleation event appears in the bulk liquid (denoted as bulk nucleation in Figure ). In this case, the hydrophilic substrate is unfavorable for the nucleation of nanobubbles on the surface, and instead the bubble is nucleated in the bulk due to the sufficiently high gas oversaturation generated by the solvent exchange. Note that the formed nanobubble in the bulk is thermodynamically unstable. Nucleation in the bulk also indicates that substrate roughness or chemical heterogeneity is not vital for nucleating nanobubbles, but it is essential for nanobubble stability.[14−22]

Figure a, as an example, shows a typical time evolution of the simulation box height at different gas concentrations in the source region, for θ = 91°. For xsourcegas = 0.002, there is no nucleation event observed, and this is confirmed by the fact that the simulation box height oscillated around its initial value. But when xsourcegas = 0.003, 0.004, and 0.005, the simulation box height increased significantly after the occurrence of the nucleation event at the substrate, after which the newly formed nanobubbles further grew. We can also find that the nucleation time decreases with increasing xsourcegas (see inset of Figure a). This is because higher xsourcegas produces higher local oversaturation that promotes nucleation. For xsourcegas = 0.003 and 0.004, nucleation of nanobubbles took place inside the pore. While for xsourcegas = 0.005, even two nuclei with one on the outer surface of the substrate and the other inside the pore were initially found and then coalesced subsequently (Figure b).

Figure 5

(a) Time evolution of the simulation box height at different concentrations of gas in the source region xsourcegas, with the inset showing the period of time for the occurrence of the nucleation event: the nucleation time tbubble ∼ 300 ns for xsourcegas = 0.003, 250 ns for xsourcegas = 0.004, and 230 ns for xsourcegas = 0.005. (b) Several typical snapshots show the coalescence of small bubbles that simultaneously appeared at different places of the bottom substrate in the case with a high xsourcegas = 0.005. In this figure we set θ = 91°.

Different pathways for nanobubble nucleation for different values of θ at xsourcegas = 0.004 are shown in Figure . Figure a–c show corresponding snapshots for nanobubbles nucleated from different pathways, that is, nucleation on the substrate but out of the pore (θ = 130° in Figure a), nucleation inside the pore (θ = 91° in Figure b), and nucleation in the bulk solution (θ = 31° in Figure c). From Figure a we can find that the nucleation took place both on the outer surface of the substrate and inside the pore nearly at the same time, but the smaller nucleus finally disappeared as the bigger one grew. In this case, it is the high hydrophobicity of the substrate (θ = 130°) that induces the nucleation. For θ = 91°, nucleation always started from the pore after a longer waiting time, and the nucleus fluctuated strongly until its size exceeded the critical size and then grew spontaneously to fill the pore, followed by growth into a nanobubble without significant energy barrier. For the case of θ = 31°, the nucleus was first generated from the bulk liquid and increased its size for a long time, showing random Brownian movement (Figure c). Ultimately, it contacted the solid interface as a result of directed diffusion of gas molecules.

Conclusions

While the solvent exchange procedure has become the most-used protocol to produce surface nanobubbles, the molecular mechanism behind the solvent exchange are far from being understood. Our MD simulations show a two-stage mechanism for nucleating nanobubbles via solvent exchange. During the first stage an interface between two interchanging solvents is found. The interface moves toward the substrate gradually as the exchange process proceeds. We find that there exists directed diffusion of gas molecules against gas concentration gradient driven by the solubility gradient of liquid composition across the moving solvent–solvent interface. The forced diffusion against the gas density gradient prevents the gas molecules from washing away and produces a locally very high gas oversaturation near substrates. As a result, the locally high gas oversaturation nucleates nanobubbles via different pathways of forming nanobubbles initially either on the solid surfaces or in the bulk solution, depending on the substrate hydrophobicity and the degree of local gas oversaturation; see the phase space in Figure . The findings of our work suggest that solvent exchange could be developed into a standard procedure to produce oversaturation and used to a variety of nucleation applications other than generating nanobubbles.