A multi-scale perspective of gas transport through soap-film membranes

Abstract

Soap films represent unique aqueous systems, whose physical properties can be tuned by acting on their nanoscale structure. Here, we specifically focus on transport properties through membranes realized in the form of soap films. While diffusion phenomena in the water core and surfactant monolayers are described using a continuum model, molecular dynamics is used to compute the static and dynamical properties of water, gases and the surfactant in the monolayers which is hexaethylene glycol monododecyl ether (C₁₂E₆). The obtained atomistic details are then incorporated into a drift-diffusion model for consistently extracting a boundary condition for the above continuum model describing transport phenomena at a larger scale. Numerical predictions are validated against experimental data from both properly designed experiments and the literature. Finally, the developed model is used to estimate the characteristic time for disparate gas mixing when initially separated by soap film membranes.

---

Design, System, Application

Here we target the accurate molecular design of thin water layers (whose thickness can be in the range from hundreds down to a few nanometers) realized in the form of soap films, for accurately controlling gas permeation. This is conducted by setting up a comprehensive multi-scale model to elucidate the role of all main design parameters (including the choice of self-assembled surfactant molecules, their concentration and film thickness). Possible future applications that can be envisioned for this work are in the energy sector. In fact, this work is part of a broader multi-disciplinary research project (http://sofiaproject.eu/), recently funded by the European Commission (grant agreement No. 828838), which aims at designing flexible (soap-film based) compartments for oxygen-fuel (mainly CO) generation (by photoreactions) and separation starting from coalescent bubbles initially filled with CO₂. In order to accomplish the application above, soap-film membranes should be able to keep separated disparate gases (e.g. O₂–CO) generated into separated compartments by photoreactions for dozens of seconds or longer preventing substantial oxygen-fuel mixing.

1 Introduction

Thin layers of liquid water within a gaseous environment can be generated in the form of soap films, with a range of thickness from hundreds down to a few nanometers (i.e. Newton black films).^1–3 Stabilization of soap films is due to the presence of surfactant molecules that self-assemble at the gas–liquid interface, thus lowering surface tension.⁴ The presence of these intricate assemblies of surfactant molecules makes the transport of species across soap films a non-trivial phenomenon.

Gas permeation through soap films has been extensively studied since the seminal work by Princen and co-workers^5,6 for a number of common gaseous species. Moreover, due to their high selectivity, soap films and bubbles have been also investigated as possible gas mixture separation membranes.^7,8 Despite the several studies reported in the literature, it is fair to say that most of the theoretical models that have been proposed to describe gas transport across soap films still require the evaluation of a number of empirical parameters.⁹ In this work, we propose a comprehensive multi-scale approach that is capable of consistently coupling molecular dynamics simulation and a continuum model in order to predict soap film permeability, which is a critical parameter affecting gas permeation across soap-film based membranes. In particular, the molecular simulations provide crucial input such as the interfacial thickness and the molecular interactions that both the water molecules and gas molecules (CO₂) make with the surfactant. This manuscript is organized in sections as follows. In section 2 the atomistic model is described and validated against experimental data reported in section 3. The continuum model is described in section 4 and validated in section 5, while conclusions are drawn in section 6.

2 Molecular model of the surfactant–water system

Fig. 1 shows a snapshot of the surfactant–water system that is considered in our atomistic simulations. The surfactant used for these simulations is hexaethylene glycol monododecyl ether (C₁₂E₆) which has been the subject of several previous experimental and modelling studies.^10–16 C₁₂E₆ is a surfactant consisting of a hydrophobic tail made up of 12 alkyl groups and a hydrophilic head which comprises 6 ethylene oxide (EO) groups and 1 hydroxyl (OH) group. The inset of Fig. 1 illustrates an enlarged picture of the surfactant with the head (H) and tail (T) groups.

Fig. 1 The left panel shows a snapshot taken from our simulations of the two surfactant–water interfaces. The right panel shows a blow-up of the C₁₂E₆ molecule that is used in our simulations. The labels H and T correspond to the hydrophilic head and hydrophobic tail, respectively, of the molecule.

Classical molecular dynamics simulations of C₁₂E₆ were performed using the GROMACS package.¹⁷ The force-field for the surfactant was adapted from a previous study by Striolo and co-workers. For details on the model, the reader is referred to the original work.¹⁸ Briefly, a united atom approach is used for the alkyl groups, while the EO and OH groups are treated explicitly. A combination of the TRaPPE-UA¹⁹ and OPLS^20,21 force-fields was used for the various bonded and non-bonded interaction potentials.

The system shown in Fig. 1 consists of 96 surfactant molecules, 48 on each surface, and 4055 water molecules. The water molecules were treated using the SPC/E²² potential. The choice of SPC/E was dictated by the fact that it has been already successfully implemented for reproducing the co-existence of liquid H₂O and gaseous CO₂.²³ In addition, the SPC potential accurately reproduces thermodynamic properties such as the second virial coefficient, enthalpy of vaporization and the surface tension of the air–water interface.²⁴ Periodic boundary conditions are applied in the x, y and z directions with a vacuum buffer of 215 Å between the two surfactant layers. The electrostatics were treated with the particle mesh Ewald (PME) method²⁵ using a real-space cutoff of 10 Å. The cutoff for the Lennard-Jones interactions was set to 10 Å. We studied the behavior of the surfactant close to the critical micelle concentration (CMC), which corresponds to 52 Ų per surfactant (approximately 48 surfactant molecules on each surface). The simulations were equilibrated first within the NVT ensemble for 11 ns using the Nosé–Hoover thermostat^26,27 at 298 Kelvin and using a time constant of 2.0 ps. NVT simulations were used to extract the surface tension which will be reported later in the paper. The total simulation time for these runs was 60 ns.

In order to assess the validity of our model, we performed some constant surface-tension simulations, which can be favourably compared with the experimental data reported below. Using the recently developed algorithms in ref. 28, we performed simulations constraining the surface tension under a surface normal pressure of 1 bar. In the latter simulations, the interfacial area fluctuates and allows us to assess the validity of our model.

2.1 Surface tension

From our simulations on the surfactant–water interface, we first determined the surface tension of the C₁₂E₆–water interface by examining the difference between the normal and tangential pressures defined as follows:where P_xx, P_yy and P_zz are the diagonal components of the pressure tensor along the x, y and z directions. L_z is the length of the simulation cell containing two interfaces.

The surface tension extracted from our NVT simulations reaches up to 37.8 ± 0.3 mN m⁻¹, which is in remarkable agreement with the experimental data reported in section 3. We also performed a set of simulations at a constrained surface tension of 41 mN m⁻¹ (experimental value).²⁹ Interestingly, upon monitoring the interfacial surface area, we predict a surface coverage value of 54.4 Ų, which is again consistent with our chosen critical micelle concentration of 52 Ų in the NVT runs. The above results give us confidence in the accuracy of our molecular model, allowing us to make some important inferences on the structure and dynamics of the surfactant–water interface.

2.2 Surfactant and water density profiles

In the top panel of Fig. 2, the probability density profiles associated with the water and surfactant are shown. For the surfactant, we used the center of mass of the molecules to determine the densities. The bottom panel of Fig. 2 shows a zoom-in of densities, where we separate out the contributions coming from the hydrophobic (shown in green) and hydrophilic (shown in blue) parts of the surfactant. The density profiles reveal some important details regarding the structure of the surfactant and water in the interfacial region. Firstly, we observe that the water density covers a range spanning from the bulk to the vacuum of approximately ∼20 Å. In order to quantify the width of the interfacial region, the water density profile was fitted to the following functional form as done in previous studies (Fig. 3):³⁰

Fig. 2 Mass density profile. The red curve corresponds to the water density. The surfactant density was separated into the hydrophobic tail and hydrophilic head. From the figure, the permeation of the water onto the hydrophobic region is shown. Also, the width of the tail/head layer of the surfactant is computed and presented in Table 1 using eqn (3).

Fig. 3 The water density profile was fitted to the function presented in eqn (2). With it, an average Gibbs dividing interface is defined at . Also the hydration layer is estimated from the fitted d.

Table 1 The thickness reported for the water corresponds to the “10–90” thickness parameter which is the length scale over which the density changes from 0.9ρ_L to 0.1ρ_L assuming ρ_V = 0 as discussed in the main text. For the total surfactant, head and tail contributions, the thickness reported is twice that of the full width at half maximum obtained from the Gaussian fits

	Thickness (Å)
	Our MDS	Previous MDS	Experimental
Water	17.8
All	10.15	13.4 (ref. 16)	13.25 (ref. 10 and 12)
Head	12.97	7.6 (ref. 16)	9.5 (ref. 10 and 12)
Tail	6.90	10.5 (ref. 16)	9.75 (ref. 10 and 12)

---

In eqn (2), ρ_L and ρ_V are the densities corresponding to the liquid and the gas phases, respectively, z₀ is the position of the Gibbs dividing surface (where ), and d is the thickness parameter for the interface. As in previous studies,^30,31 we computed the “10–90” thickness parameter representing the length-scale over which the density decays from 90 to 10% of the bulk value. This number reported in Table 1 is 2.1972d yielding a value of 17.8 Å.

For the surfactant and its individual head and tail distributions, the densities were fitted to the following Gaussian distribution:

The thickness associated with the surfactant part and the individual head and tail contributions that are shown in Table 1 corresponds to the full-width at half maximum (FWHM). The above results point to some relevant structural features of this interface. Firstly, the boundary between the hydrophobic and hydrophilic regions of the surfactant is rather fuzzy, which leads to a significant degree of penetration of water into the surfactant and, in particular, the hydrophobic regions. This is seen very clearly in the overlap of the water and hydrophobic density profiles in Fig. 2. Secondly, regarding the thickness of the surfactant and water interfaces, our results are in reasonable agreement with the experiments as well as previous simulations.¹⁶ Nevertheless, previous experiments using neutron reflectivity showed a certain uncertainty on the actual thickness of the hydrophobic and hydrophilic regions as well as the extent of water penetration because of the interfacial fluctuations and the roughness of the interface.^10,12

2.3 Surfactant orientation and flexibility

In order to examine the orientation of the surfactant molecules at the surface, we determined the angle θ between the vector formed by the two ends of the molecule and the surface normal reported in the top panel of Fig. 4. The distribution of θ peaks at approximately 31.6° which shows that the surfactant molecules are tilted at the interface. It should be stressed, however, that the distribution is quite broad indicating a rather heterogeneous environment.

Fig. 4 a) The tilt angle distribution, b) the end to end distance and c) the autocorrelation function associated with the end to end distance.

The middle panel of Fig. 4 shows the distribution of the end-to-end distance of the surfactant, which again confirms that the surfactant molecules are quite flexible and characterized by a wide range of different conformations. To assess the timescales associated with these fluctuations, we computed the equilibrium time correlation function associated with the changes in the conformation as probed by the end-to-end distance. The bottom panel of Fig. 4 shows the relaxation dynamics of all the surfactant molecules which features an exponential decay occurring on the timescale of 1.7 ns.

2.4 CO₂ concentration in the surfactant monolayer

We performed simulations of a single CO₂ molecule that was added to the above surfactant–water system. It is worth noticing here that, under usual operating conditions (i.e. low pressure), the concentration of CO₂ in the simulated small system corresponds to one single CO₂ molecule. A total of 5 independent trajectories were run each of 10 ns long within the NVT ensemble. The model used for the force field of the CO₂ was obtained from the literature.²⁴ Our simulations show that while the CO₂ molecule undergoes free diffusion in the gas phase, when it approaches the surfactant it appears to be trapped within the hydrophobic pockets of C₁₂E₆ as seen in Fig. 5, which shows a zoom in plot together with the CO₂ densities. The CO₂ densities are shown in solid black lines for regions that are adequately sampled while the dashed black lines are regions where we did not sample adequately enough the transitions from the hydrophobic part to the hydrophilic part and then into the water.

Fig. 5 Density profile of C₁₂E₆ (head with blue lines and tail in green) and water (red) with the probability to find CO₂ (dark line).

Due to the limited statistics involving the CO₂ molecule diffusing through the surfactant, we used the CO₂ density profiles to infer the local potential that the CO₂ experiences in the hydrophobic part of the surfactant. In order to obtain the potential, we computed the potential of the mean force associated with the density: −k_BT log(p(r)) which is shown in Fig. 6. The dashed lines correspond to the regions of the potential that are not sampled from the molecular dynamics simulations. We used this potential to infer the minimum activation barrier associated with the transfer of the CO₂ across the surfactant monolayer into the liquid. As described below, this potential is a crucial input for a continuum model to be fitted into eqn (12) and (13).

Fig. 6 Underlying potential of CO₂ moving along the z axis obtained as described in the text.

3 Surface tension and molecular area measurements

Dynamic surface tension values were measured with a drop tensiometer (Tracker, Teclis) by analyzing the axial symmetric shape (Laplace profile) of a rising air bubble in aqueous solution of C₁₂E₆ at different (bulk) concentrations. The air bubble is generated in a thermostated cuvette containing the water phase; the bubble of about 10 μl is delivered from a syringe controlled using a motor drive as shown in Fig. 7.

Fig. 7 The pendant bubble tensiometer (Tracker). (A) Motor to drive the syringe. (B) Camera and optics. (C) Motor driven syringe. (D) Light source. (E) Electronic box. (F) Software.

All measurements were performed at a controlled temperature of 25 °C. C₁₂E₆ with 98% purity from Sigma Aldrich was used without modification and dissolved in distilled water.

Results are reported in Fig. 8, where the surface tension versus time due to adsorption of the surfactant at the air/water interface at several concentrations of C₁₂E₆ is reported. For each concentration, the measurements are performed until the equilibrium surface tension is reached. The surfactant concentration is increased until its critical micellar concentration (CMC) is reached; after that point, the equilibrium surface tension is no longer concentration dependent.

Fig. 8 Dynamic surface tensions of clean surface adsorption in C₁₂E₆ aqueous solution at different concentrations.

Fig. 9 Equilibrium surface tension from Fig. 8 as a function of bulk concentration and compared with the results obtained using the Gibbs adsorption eqn (4).

The equilibrium surface tension versus the concentration of C₁₂E₆ (logarithmic scale) is plotted in Fig. 9. The value of the CMC obtained for C₁₂E₆ is 0.034 g l⁻¹ (0.075 mM) with a minimum surface tension of 31 mN m⁻¹.

The data in Fig. 9 are compared to the Gibbs adsorption equation. For non-ionic systems, before the CMC, the amount of surfactant adsorbed at the interface at equilibrium (Γ_eq) is related to the surfactant (bulk) concentration C_b by:

The best fit was obtained with the surface concentration at CMC of 0.0032 mmol m⁻² leading to a molecular area at CMC of 52 Ų.

4 Permeation of gaseous species across soap films

Among the various theoretical models for mass transport across a soap film, a widely accepted description in the literature is the one proposed by Princen and co-workers,^5,6 who further elaborated this in previous studies.^32,33 In this model, permeation of gas molecules through soap films depends on three main resistances due to the liquid water core inside the film and the two surfactant monolayers at the gas–liquid interface. Assuming Fick's diffusion, these authors suggested a linear proportionality between the molar flux of permeating gaseous species and the corresponding concentration difference across the soap film, namely:with n being the number of moles, k the soap film permeability in m s⁻¹, A the surface area of the monolayer in m² and (C₁ − C₂) the concentration difference across the film in mol m⁻³.

Clearly, the permeability depends on the thickness of the liquid water core. However, k also includes a non-trivial dependence on the surfactants and the involved gas, and it can be more specifically linked to other relevant quantities as follows:

with k_ML1 and k_ML2 being the monolayer permeability of the two surfactant monolayers in m s⁻¹ forming the soap film, D_w the diffusion coefficient of the gas in water in m² s⁻¹, H the Henry's law constant (dimensionless) and h the water core thickness in m. For the sake of generality, it is worth stressing that a lack of symmetry in the surfactant monolayers can be reflected in eqn (6) as a disparity in the monolayer permeability values. However, for most of the practical cases the soap film is symmetric, hence k_ML1 = k_ML2.

The monolayer permeability k_ML1 depends on the molecular structure of surfactant molecule (e.g. polarity of the hydrophilic head, length of the hydrophobic tail), on permeating gas, on temperature and surfactant concentration.⁹ In the case of thick films the permeability is controlled by the liquid core yielding:

On the other hand for very thin films (Newton black films) the water core resistance is negligible, so that:⁵

4.1 Monolayer permeability

Focusing specifically on the gas–film interface, the molar rate of a given gas diffusing through a single surfactant monolayer can be estimated by the following equation:with C′ and C̄ being the actual concentration and the equilibrium concentration of the considered gaseous species in the film liquid core, respectively.⁶Eqn (9) predicts a net gas flux from the gas to the liquid phase, which only stops when the following equilibrium condition (Henry's law) holds:with C denoting the actual concentration in the gas phase.

Clearly, in order to have predictive results from eqn (9), there is a need to link the monolayer permeability to more fundamental quantities. In this respect, a rather popular approach is as follows:⁹

where D_ML denotes an effective diffusion coefficient within the surfactant monolayer whose thickness is δ.

However, it is fair to say that the physical process of gas permeation through surfactant monolayers is not yet completely understood. In fact, as testified in the review by Farajzadeh and co-workers, disparate theoretical models have been suggested so far in the literature, thus also introducing a number of empirical parameters that are often hard to estimate.⁹

It is interesting to notice that, even when no surfactant molecules are present, classical atomistic simulations suggest that the transport from the gas to the liquid phase occurs within a non-flat potential landscape. In this respect, authors in ref. 34 and 35 report a free energy barrier at the gas–water interface. Hence, based on the above discussion in section 2, we expect an even more complicated picture when also surfactant molecules are present and specifically suggest investigating the perturbation in the free-energy landscape to be linked to eqn (11). Towards this end, here we propose a comprehensive and multi-scale perspective, where a molecular model is utilized for extracting relevant information to be used at a larger scale in a continuum framework.

In particular, throughout the surfactant monolayer, gas molecules are subject to an underlying energy potential. Hence, gas transport can be effectively described by means of the following Smoluchowski model:³⁶

where C_ML is intended to be the gas concentration in the monolayer, D_w the diffusion coefficient of the gas in water in m² s⁻¹, R the universal gas constant in J mol⁻¹ K⁻¹ and T the temperature in K. The function U(x) denotes the underlying potential experienced by CO₂ molecules at a generic position x in the monolayer.

Assuming for simplicity a one-dimensional (1D) stationary process along the coordinate z orthogonal to the monolayer, eqn (12) and (13) yield:

We propose to estimate the effective diffusion coefficient in the monolayer, D_ML, by imposing flux conservation between eqn (15) and the following (effective) Fickian description:Hence, the monolayer permeability can be calculated by eqn (11), where δ represents the width of the potential barrier (along the z coordinate).

5 Model results and experimental validation

Firstly, we set-up and validated a (COMSOL³⁷ based) continuum model by using the experimental data from Princen and co-workers,⁵ who accurately characterized the gas diffusion across the soap film of a shrinking bubble (diminishing bubble method). In their experimental setup, a bubble containing a specific gas (i.e. N₂, O₂, H₂, CO₂, He, Ne, Ar and air) was formed inside a closed cell filled with the same gas using a 4% solution of CTAB (hexadecyltrimethylammonium bromide). Once the bubble film had reached the equilibrium thickness, they measured the time variation of the following dimensionless parameter β₀, defined as:with d being the water density, g the acceleration of gravity, γ = 33 mN m⁻¹ the surface tension and r the radius of a sphere of the same volume as the bubble. From these data, we can estimate the gas flow from the shrinking bubble to the outer atmosphere subject to the following Young–Laplace pressure difference:with X being the curvature radius of the bubble and γ the surface tension.

In our continuum model, we focused on the gas transport across the soap film whose thickness is in the order of tens to hundreds of nanometers. A 2D axial-symmetric geometry was chosen as shown in Fig. 10: the lower chamber corresponds to a small section of the bubble's inner gas volume, while the upper chamber represents a small section of the outer gas volume. Fickian diffusion was assumed both in the gas phase and in the soap film. We considered the presence of the surfactants at the air–water interface by implementing eqn (9) in our model. Furthermore, we imposed a fixed pressure difference between the two chambers expressed in terms of concentration difference by the perfect gas law:

Finally, we fix the water core thickness at 5 nm as estimated in ref. 5.

Fig. 10 Left-hand side: A sketch of the shrinking bubble in the Princen³⁸ experimental setup. Middle: Continuum model results for oxygen diffusion through the soap film. Right-hand side: O₂ concentration distribution across the soap film.

By varying the ΔC, we calculate the molar flux reduction caused by the bubble shrinking with time. In Fig. 11, our simulation results were validated against the experimental data from ref. 5 by using both k_ML = 0.188 m s⁻¹ and k_ML = 0.242 m s⁻¹ as suggested in ref. 5.

Fig. 11 Dependence of the molar flux of oxygen on the permeability for a film thickness of 5 nm against experimental data from ref. 5.

5.1 Monolayer permeability, k_ML

In this section, we coupled the above atomistic description with a continuum model implementing the Smoluchowski diffusion equation. In particular, we solved eqn (12) and (13) in the COMSOL environment under 1D stationary conditions. Firstly, the free-energy landscape from the molecular dynamics simulations reported in Fig. 6 was fitted by using both a harmonic and a higher polynomial model, U₁ and U₂, as shown in Fig. 12. The main motivation for choosing a harmonic potential can be found in the literature, where gas permeation in the surfactant monolayer is described as a phenomenon mainly dictated by the activation energy (i.e. energy barrier height).^39–41 In addition to this, in our work, we decided also to apply a more accurate fitting, U₂.

Fig. 12 Harmonic (U₁ = a₁x² + a₂ × x + a₃) and polynomial (U₂ = b₁x⁷ + b₂x⁶ + b₃x⁵ + b₄x⁴ + b₅x³ + b₆x² + b₇x + b₈) fitting corresponding to the free energy profile from MD simulations where E_a is the activation energy. Polynomial coefficients: a₁ = 4.890, a₂ = −9.292, a₃ = 4.4135, b₁ = 7.583, b₂ = −38.670, b₃ = 63.370, b₄ = −24.720, b₅ = −21.140, b₆ = 12.960, b₇ = −3.877, and b₈ = 4.4135. x is in nm and U in kJ.

The simulation domain length was chosen on the basis of the potential barrier width, 1.9 nm, as shown in Fig. 12. The diffusion coefficient of CO₂ in water is assumed to be 1.67 × 10⁻⁹ m² s⁻¹ at T = 20 °C.⁴² We imposed Dirichlet boundary conditions at the ends of the simulation domain.

We determined the gas flux as a function of the concentration difference in the case of flat and non-flat potentials. We notice that, as shown in Fig. 13, the presence of asymmetric potentials may lead to a non-vanishing flux (J_S,U2) even at ΔC = 0. In contrast, for symmetric potentials, J_S,U1 = 0 is found at ΔC = 0.

Fig. 13 Molar flux depending on the concentration difference across the monolayer, diffusion under: a flat potential (blue triangles), a harmonic potential U₁ (green circles), and a polynomial potential U₂ (red squares).

Finally, we estimated the effective diffusion coefficient in the C₁₂E₆ surfactant monolayer by measuring the slope of curves reported in Fig. 13. Here, the diffusion coefficient of a gas within water is denoted as D_w, whereas D_ML is used as the effective diffusion coefficient in the presence of an underlying potential. By using potential U₁ and δ = 1.9 nm, we found a value for of 0.478 × 10⁻⁹ m² s⁻¹, with the corresponding monolayer permeability being = 0.25 m s⁻¹.

On the other hand, was found to be 0.541 × 10⁻⁹ m² s⁻¹ when implementing potential U₂, which led to a slightly higher monolayer permeability of = 0.28 m s⁻¹. Thus, the shape of the potential is found to have a limited influence on the resulting monolayer permeability (as compared to its height).

As already discussed, due to the possible poor sampling issues in the computation of the energy potential landscape in Fig. 6, we decided to perform additional computations assuming an underlying harmonic potential with a perturbed height and width as compared to the one reported in Fig. 12. More specifically, we assumed a conservative ±10% variation to both the activation energy E_a reported in Fig. 12 and the monolayer thickness δ. The above analysis provided a range of values for the monolayer permeability from 0.12 up to 0.36 m s⁻¹. To the best of our knowledge, values of the monolayer permeability of C₁₂E₆ have not been reported yet for CO₂ and under the conditions reported in this work.

However, our results appear to be consistent with similar studies in the literature. For instance, Krustev and co-workers^43,44 reported the k_ML from experimental measurements of the permeability for air through SDS (sodium dodecyl sulfate) + electrolyte solutions obtaining values which ranged from 0.063 to 0.278 m s⁻¹ depending on the type of electrolyte and its concentration.

5.2 Sensitivity analysis, film thickness and monolayer permeability

So far we have characterized the CO₂ behaviour in the surfactant monolayer. Here, we decided to carry out a wider sensitivity analysis on the monolayer permeability and on the water core thickness to extend our research and investigate the influence of different gases and surfactants on the transport phenomena across soap films.

The same simulation setup validated previously was adopted imposing a concentration difference of ΔC = Δp/RT = 0.0184 mol m⁻³ corresponding to a pressure difference of 45 Pa at the boundaries. We let both the monolayer permeability and the water core thickness vary from 0.05 up to 0.5 m s⁻¹ and from 5 nm up to 500 nm, respectively. We highlight that the film thickness has a major impact on gas diffusion as compared to the monolayer permeability. For instance, for h = 300 nm, a monolayer permeability variation from 0.05 to 0.5 m s⁻¹ reduces the flux by about 18%. Nevertheless, as the water core becomes thinner the effect of the monolayer permeability becomes more significant: e.g. for a film of about 50 nm an increment of the monolayer permeability from 0.05 to 0.5 m s⁻¹ halves the gas flux and this behaviour is even emphasized for thinner films as shown in Fig. 14.

Fig. 14 Influence of the gas monolayer permeability and water core thickness on the dimensionless molar flux. The diffusive flux was normalized by the maximum flux, equal to 9.303 × 10⁻⁵ mol m⁻² s⁻¹, obtained with the maximum value of k_ML = 0.5 m s⁻¹ and for the minimum value of h = 5 × 10⁻⁹ m.

5.3 Estimate of the cross-flow characteristic time through soap-film membranes

In this section, we apply the above models to investigate the characteristic time for gas mixing in two coalescent bubbles with different gas compositions. In particular, we consider the case of a double-bubble configuration which consists of two coalescent spherical shells of soap films (with an equal size and radius of 2 cm) separated by a flat soap film (surface area of 9.425 cm²).^45,46 A lower bound for the time needed before 10% of the initial number of moles of gas i inside one of the two bubbles flows to the other was evaluated as:where n_i = p_iV/RT is the initial number of moles of gas i in one of the bubbles and J_i is the maximum flux of gas i across the film from a bubble to the other. We assumed a 9.7% O₂, 88% CO₂ and 2.3% water vapor atmosphere for the first bubble and a 19.5% CO, 78.2% CO₂ and 2.3% water vapor atmosphere for the second one. The self-diffusion coefficients of gases are reported along with the Henry's law constants and the diffusion coefficients of gases in water in Table 2. We used the k_ML for CO₂ calculated before and equal to 0.25 m s⁻¹. Princen and co-workers⁵ observed that for a given surfactant, the monolayer permeability depends slightly on the gas nature; in contrast the same coefficient is strongly affected by the surfactant.⁴⁴ Thus, since no specific literature data were found, we safely decided to set the monolayer permeability of oxygen and carbon monoxide in C₆E₁₂ to 0.278 m s⁻¹, which is close to the value reported in ref. 44.

Table 2 Diffusion coefficients and Henry's law constants

Gas	D self	D w	H	Source
	m^2 s^−1	m^2 s^−1	—	Source
CO2	1.13 × 10^−5	1.67 × 10^−9	0.9520	Lide ref. 42
O2	2.09 × 10^−5	2.01 × 10^−9	0.0336	Lide ref. 42
CO	1.90 × 10^−5	2.03 × 10^−9	0.0256	Lide ref. 42

---

In order to have a conservative estimate, we considered a fixed molar flux constantly driven by the (initial) maximum gas concentration (ΔC(t = 0 s)). Hence, we calculated the time spent to reduce the initial number of moles of the i-th gas inside one of the two bubbles by 10% according to eqn (20). We set the water core thickness to 300 nm and, by applying an uncertainty to the monolayer permeability of ±50%, t_O2 was found to vary from 13.7 s up to 14.6 s, whereas t_CO was found to vary from 17.9 s up to 19.0 s. It is worth stressing once again that, due to the assumption on the fixed gas concentration, the above mixing times are to be considered as lower bounds.

6 Conclusions

Thin layers of water can be easily realized in the form of soap films. Interestingly, the surfactant nature and the thickness can be easily varied. As a representative example, a thin film pressure balance (TFPB) represents a convenient system where nanoscale properties can be precisely controlled and investigated in a relatively easy manner as shown in numerous studies.^47–52 In this work, we have focused on the gas transport across soap films by means of predictive multi-scale numerical models. In particular, we used a combination of molecular and continuum models to characterise and study the transport of gases through soap film membranes. Classical molecular dynamics simulations are used to study the structural and dynamical properties of the surfactant molecules at the air–water interface. The above atomistic model is validated against experimental data from purposely designed experiments in terms of surface tension and molecular area. Simulation results provide crucial insights into the thickness of the hydrophobic and hydrophilic parts of the surfactant as well as the extent of water penetration in the membrane. The molecular simulations of CO₂ at the interface reveal that the gas tends to get trapped within the hydrophobic part of the surfactant. These insights provide valuable input for studying the diffusion of CO₂ through a soap film using a continuum model. Specifically, the free-energy landscape is incorporated in a Smoluchowski model at a continuum level and utilized to estimate the surfactant monolayer permeability without resorting to empirical parameters. The above procedure is implemented (although it is not limited) to CO₂, and it can be easily extended to other gaseous species. Finally, a lower bound for a characteristic mixing time of two different gases initially separated into two coalescent bubbles is provided.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work has been funded by the European Commission in the framework of the project SoFiA – Soap Film based Artificial Photosynthesis – grant agreement 828838.

References

1. D. Exerowa and P. M. Kruglyakov, Foam and foam films: theory, experiment, application, Elsevier, 1997, vol. 5.
2. J. Clunie, J. Corkill and J. Goodman, Discuss. Faraday Soc., 1966, 42, 34–41.
3. G. Andreatta, L.-T. Lee, F. K. Lee and J.-J. Benattar, J. Phys. Chem. B, 2006, 110, 19537–19542.
4. I. J. Black and T. M. Herrington, J. Chem. Soc., Faraday Trans., 1995, 91, 4251–4260.
5. H. Princen and S. Mason, J. Colloid Sci., 1965, 20, 353–375.
6. H. Princen, J. Overbeek and S. Mason, J. Colloid Interface Sci., 1967, 24, 125–130.
7. R. L. Cook and R. W. Tock, Sep. Sci., 1974, 9, 185–193.
8. J. Reznickova, R. Petrychkovych, J. Vejrazka, K. Setnickova and P. Uchytil, Sep. Purif. Technol., 2016, 166, 26–33.
9. R. Farajzadeh, R. Krastev and P. L. Zitha, Adv. Colloid Interface Sci., 2008, 137, 27–44.
10. J. R. Lu, Z. X. Li, R. K. Thomas, E. J. Staples, I. Tucker and J. Penfold, J. Phys. Chem., 1993, 97, 8012–8020.
11. J. Lu, E. Lee, R. Thomas, J. Penfold and S. Flitsch, Langmuir, 1993, 9, 1352–1360.
12. J. R. Lu, T. J. Su, Z. X. Li, R. K. Thomas, E. J. Staples, I. Tucker and J. Penfold, J. Phys. Chem. B, 1997, 101, 10332–10339.
13. J. Lu, Z. Li, R. Thomas, B. Binks, D. Crichton, P. Fletcher, J. McNab and J. Penfold, J. Phys. Chem. B, 1998, 102, 5785–5793.
14. J. Chanda and S. Bandyopadhyay, J. Chem. Theory Comput., 2005, 1, 963–971.
15. R. Muruganathan, H.-J. Müller, H. Möhwald and R. Krastev, Langmuir, 2005, 21, 12222–12228.
16. L. Shi, N. R. Tummala and A. Striolo, Langmuir, 2010, 26, 5462–5474.
17. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1, 19–25.
18. L. Shi, N. R. Tummala and A. Striolo, Langmuir, 2010, 26, 5462–5474.
19. M. G. Martin and J. I. Siepmann, J. Phys. Chem. B, 1999, 103, 4508–4517.
20. J. M. Briggs, T. Matsui and W. L. Jorgensen, J. Comput. Chem., 1990, 11, 958–971.
21. W. L. Jorgensen, J. Phys. Chem., 1986, 90, 1276–1284.
22. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
23. T. Somasundaram, M. In het Panhuis, R. Lynden-Bell and C. Patterson, J. Chem. Phys., 1999, 111, 2190–2199.
24. H. Zhang and S. J. Singer, J. Phys. Chem. A, 2011, 115, 6285–6296.
25. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
26. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
27. W. G. Hoover, Phys. Rev. A: At., Mol., Opt. Phys., 1985, 31, 1695–1697.
28. K. M. Aoki, M. Yoneya and H. Yokoyama, J. Chem. Phys., 2006, 124, 064705.
29. E. Lucassen-Reynders, A. Cagna and J. Lucassen, Colloids Surf., A, 2001, 186, 63–72.
30. R. S. Taylor, L. X. Dang and B. C. Garrett, J. Phys. Chem., 1996, 100, 11720–11725.
31. A. M. Prpich, Y. Sheng, W. Wang, M. E. Biswas and P. Chen, PLoS One, 2009, 4, 1–6.
32. A. Brown, W. C. Thuman and J. McBain, J. Colloid Sci., 1953, 8, 491–507.
33. A. Brown, W. C. Thuman and J. McBain, J. Colloid Sci., 1953, 8, 508–519.
34. H. Zhang and S. J. Singer, J. Phys. Chem. A, 2011, 115, 6285–6296.
35. R. Vácha, P. Slavíček, M. Mucha, B. J. Finlayson-Pitts and P. Jungwirth, J. Phys. Chem. A, 2004, 108, 11573–11579.
36. P. Ansalone, M. Chinappi, L. Rondoni and F. Cecconi, J. Chem. Phys., 2015, 143, 154109.
37. COMSOL Multiphysics ®v. 5.5 www.comsol.com, COMSOL AB, Stockholm, Sweden.
38. H. Princen and S. Mason, J. Colloid Sci., 1965, 20, 156–172.
39. I. Langmuir and D. B. Langmuir, J. Phys. Chem., 1927, 31, 1719–1731.
40. I. Langmuir and V. J. Schaefer, J. Franklin Inst., 1943, 235, 119–162.
41. G. T. Barnes and V. K. La Mer, in Retardation of Evaporation by Monolayers, Elsevier, 1962, pp. 9–33.
42. D. R. Lide, CRC handbook of chemistry and physics, CRC press, 2004, vol. 85.
43. R. Krustev, D. Platikanov and M. Nedyalkov, Colloids Surf., A, 1993, 79, 129–136.
44. R. Krustev, D. Platikanov and M. Nedyalkov, Colloids Surf., A, 1997, 123–124, 383–390.
45. J. Hass and R. Schlafly, Ann. Math., 2000, 151, 459–515.
46. R. Urikovič, Comput. Graph. Forum, 2001, 67–76.
47. K. J. Mysels and M. N. Jones, Discuss. Faraday Soc., 1966, 42, 42–50.
48. B. M. Folmer and B. Kronberg, Langmuir, 2000, 16, 5987–5992.
49. N. Kristen and R. von Klitzing, Soft Matter, 2010, 6, 849–861.
50. H. Fauser and R. von Klitzing, Soft Matter, 2014, 10, 6903–6916.
51. A. Mamane, E. Chevallier, L. Olanier, F. Lequeux and C. Monteux, Soft Matter, 2017, 13, 1299–1305.
52. D. Langevin, Adv. Colloid Interface Sci., 2019, 102075.

---

Footnote

† These authors contributed equally to this work.

---