Phase dependence of effective surface diffusivity in surfactant monolayers dilatated far from equilibrium

Abstract

Modeling the flow of insoluble surfactant monolayers undergoing rapid expansion and compression remains a long-standing challenge in hydrodynamics. In a recent letter, compressible flow of a DPPC monolayer in the coexisting liquid-expanded/liquid-condensed phase was accurately described with Newtonian and Fickian modeling by using large effective surface diffusivity, many orders of magnitude larger than the monolayer's equilibrium diffusivity, along with finite surface dilatational viscosity. The large effective diffusivity was attributed to monolayer phase coexistence. Here, the applicability of large effective diffusivity is investigated through flow measurements of DPPC and vitamin K1 monolayers over a wide range of concentrations, from gas phase to collapse. Effective surface diffusivity and surface viscosity were determined by fitting numerical simulations, via Navier–Stokes bulk flow coupled to a Boussinesq–Scriven interface with monolayer advection–diffusion, to spatio-temporal surface velocity measurements. Results reveal that large effective surface diffusivity is not limited to coexisting phases, producing good fits for all monolayers outside gaseous phases. The failure of this theoretical framework when applied to gaseous monolayers is likely due to their extreme compressibility. The extreme diffusivity in non-gaseous phases motivates future exploration of non-Newtonian and non-Fickian interfacial models.

---

1 Introduction

Insoluble surfactants (surface active agents) energetically favor the formation of monomolecular films at air–liquid interfaces, which in turn determine the hydrodynamic response of the system.^1–4 The presence of surfactant monolayers affects the transport of mass, momentum, and energy, with the monolayer acting as a boundary condition coupled to the flow of the bulk fluid.^5–7 Such effects are important in time-varying flow phenomena, including thin film stability,^8,9 coalescence/breakup,^10,11 wetting,^12,13 spreading,^14,15 and material properties.^16–18 Such phenomena are significant in numerous fluid systems that are far from equilibrium, ranging from our lung's liquid lining¹⁹ and the ocean's surface²⁰ to pharmaceutical production²¹ and inkjet²² printing. Respiration in the lung is a small-scale example that demonstrates the importance of surfactants in the dynamics of air–liquid interfaces. During respiration, the alveolar air–liquid interface, where gas exchange occurs, rapidly expands and contracts, while surfactants allow its functioning through modulating surface tension and surface viscosity.^19,23,24 A large-scale example of surfactants influencing flow is the significantly altered breaking behavior of meter-scale gravity waves in the presence of surfactants.²⁵

When surfactant monolayers are sheared but not expanded or compressed, i.e., not dilatated, their flow response is well described by either Newtonian or non-Newtonian constitutive relations involving only surface shear viscosity.^26–31 Numerical models have successfully predicted the flow of many different monolayers, both Newtonian and non-Newtonian, in systems driven at the interface or driven from the bulk.^5,32–35 Complications arise when monolayers are dilatated due to the simultaneous effects of surface tension gradients (Marangoni stress), surface diffusivity, surface shear viscosity, and surface dilatational viscosity.^1,36–43 These complications have produced anomalous results, including negative values of surface dilatational viscosity.^36,44,45 While equilibrium surface tension has been investigated for numerous insoluble monolayers,^46–48 inconsistencies such as negative surface dilatational viscosity remain, which may be rooted in inaccurate modeling of surface diffusivity under dynamic conditions. Enhancement of effective surface diffusivity above measured equilibrium values^49,50 and theoretical predictions⁵¹ has been observed in systems driven even slightly from equilibrium. Increases of 2 to 3 orders of magnitude have been reported in a steady-state thermo-Marangoni flow,⁵⁰ and 4 to 5 orders of magnitude in a relatively slowly oscillating Langmuir trough⁵² (oscillation period ∼30 s). Furthermore, in a recent study, surface diffusivity 7 orders of magnitude larger than equilibrium values was required to fit numerical simulations to surface velocity measurements of a DPPC monolayer driven far from equilibrium (oscillation period of ∼1 s with large amplitude).⁵³ The large magnitude of this effective diffusivity was attributed to dynamics of the coexisting liquid-expanded/liquid-condensed phase.

In the present study, effective surface diffusivity is investigated for surfactant monolayers hydrodynamically dilatated far from equilibrium using experimental measurements and numerical simulations. Experiments were performed over a wide range of monolayer concentrations, from gaseous phases to collapse, for two insoluble surfactants: 1,2-dipalmitoylphosphatidylcholine (DPPC) and phytomenadione (vitamin K1). The phospholipid DPPC was selected for its physiological relevance as the main constituent of lung surfactant, rich interfacial phase behavior, and extensive use in hydrodynamic research.^{23,24,54–56} The non-ionic surfactant vitamin K1 was selected for its multifaceted phase behavior and prevalence as a model monolayer.^36,46,47 This study examines the compressible flow of insoluble surfactant monolayers by focusing on interfaces that are flat and fixed in space. These restrictions circumvent complications for mechanistic modeling of air–liquid interfaces resulting from interface creation,⁵⁷ as occurs in Langmuir troughs, interfacial curvature,^57,58 as occurs in oscillating drops and bubbles, and surfactant transport to and from the bulk liquid,^59,60 as occurs in soluble surfactants. These constraints are achieved in the oscillatory-driven channel presented in Fig. 1, where spatio-temporal forcing of surfactant monolayers is a result of viscous traction by bulk flow with large inertia via a sinusoidally-driven floor. The rapid oscillation periods (∼1 s) and large amplitudes in this study force monolayers further away from equilibrium than most earlier work,^49,50,52 and thus require fluid inertia be accounted for using the full Navier–Stokes equations in numerical simulations, including nonlinear terms. Coupling between the bulk and interface is achieved using the Boussinesq–Scriven surface model along with Fickian diffusion of the monolayer. While other interfacial models exist, including non-Newtonian^{33,34,44,61–64} and non-Fickian^64–66 models, these often introduce additional parameters and inconsistencies between models.⁶⁷ This investigation seeks to identify a parameter space for surface diffusivity and surface viscosity that is self-consistent within a Newtonian interfacial model, with Fickian diffusion of the monolayer, coupled to Navier–Stokes bulk flow.

Fig. 1 Schematic of the inertia-dominated flow in the oscillatory-driven channel showing the x–y coordinate system, along with depth H, width W (= 2H), and a DPPC monolayer at the air–water interface. The floor is driven sinusoidally with maximum speed U and frequency ω (rad s⁻¹). Computed bulk flow (blue – minimum, red – maximum), along with representative streamlines, depicts the flow at the instant the floor is at maximum leftward travel, compressing the monolayer to the right.

2 Materials and methods

2.1 Materials

DPPC was sourced from a 25 mg mL⁻¹ stock solution in chloroform (Avanti Research 850355C-1g) and vitamin K1 from an undiluted liquid (V3501-1G, Sigma-Aldrich). All monolayer-forming surfactants were prepared into 0.1 mg mL⁻¹ working solutions in ACS-grade chloroform, which was used as a spreading agent.⁴⁶ Deionized (DI) water, with 18+ MΩ cm resistivity and less than 5 ppb dissolved organics, was used as the bulk fluid in all experiments. The materials comprising the oscillatory-driven channel were Delrin walls and a glass floor. Seal between the channel walls and floor was achieved using silicone gaskets. None of the materials used in the oscillatory-driven channel resulted in any measurable contamination of the air–liquid interface.

2.2 Experimental methods

All the experiments were conducted at 22.5 ± 0.5 °C. Measurements of surface tension equations of state were performed via slow compression of each monolayer in a custom-built Langmuir trough made of Teflon. Wilhelmy plate measurements of surface tension were obtained using an electrobalance (NIMA type PS4) and a plate comprised of clean, freshly-cut filter paper. Compressions were performed after a 10-min wait with an initial area of 161 cm² quasistatically compressed down to 21 cm² at the slow rate of 0.2 cm² s⁻¹, ensuring the system remained near equilibrium. These compressions were performed in triplicate, averaged, and subsequently fit using bimodal 5-parameter logistic functions that produced smooth functions that could be input into numerical simulations. Specifically, the form of these sigmoidal fits was: , where a_n is the response at zero concentration, b_n is the slope factor, c_n is the inflection point, d_n is the response at infinite concentration and e_n is an asymmetry factor. The measured equation of state for DPPC is presented in Fig. 2a. Interfacial compressibility, C_s = −(∂c/∂σ)/c, was calculated based on the fit to the equation of state, and is presented for DPPC monolayers in Fig. 2b.^68,69

Fig. 2 (a) DPPC surface tension equation of state σ(c), measured quasistatically via slow compression in a Langmuir trough, with an accompanying logistic function fit† used in numerical simulations. Experimental values are the average of 3 replicate measurements, with the gray error shade region indicating one standard deviation. Labels indicate midpoints of monolayer phases: gas (G), liquid expanded (LE), coexisting liquid expanded/condensed (LE/LC), liquid condensed (LC), solid condensed (SC), and point of collapse (C). Area per molecule (A_m) is provided below the concentration (c) axis for reference. (b) DPPC monolayer compressibility C_s calculated based on the fit to the equation of state.

The flow device was the oscillatory-driven channel, of depth H (= 0.5 cm) and width W (= 1 cm) with a sinusoidally-driven floor, as shown in Fig. 1. The span of the channel (7H in the z-direction) is sufficiently large so that flow in the midspan region is spanwise invariant, allowing for comparison between experiments and two-dimensional simulations.⁷⁰ Spatio-temporally resolved surface velocity field measurements were made across the full width of the channel. The sinusoidal oscillations of the floor were at angular frequency ω with displacement amplitude δ and the corresponding maximum speed U = δω. Experiments were performed for two flow conditions corresponding to combinations of δ and ω of (1.5 cm, 2.5 rad s⁻¹) and (1.5 cm, 7.4 rad s⁻¹). The dimensionless flow amplitude is given by the Reynolds number Re = ρUH/μ, where ρ is the bulk fluid density taken to be 0.998 g cm⁻³ for water, and μ is the dynamic viscosity of water taken to be 0.942 cP. The dimensionless frequency is Ω = ωH/U. The two driving flow conditions in the experiments correspond to Re = 200 and Re = 600, both at Ω = 1/3. In this range of flow parameters, the air–liquid interface remains essentially flat. The flatness of the interface was verified using a laser surface-reflectivity probe, consisting of a laser diode and a high-speed camera. At these forcing amplitudes and frequencies, the flow in the oscillatory-driven channel is known to be hydrodynamically stable.^70–73

Flow measurements were made for mean monolayer concentrations c₀ ranging from 0.1 to 2.2 mg m⁻² in increments of 0.1 mg m⁻² for DPPC, and 0.2 to 1.6 mg m⁻² in increments of 0.2 mg m⁻² for vitamin K1. All monolayers were spread onto the interface using the standard technique of dissolution in a water-immiscible solvent with high vapor pressure.⁴⁶ Subsequently, the floor was oscillated for 10 minutes before the start of data acquisition to allow for complete evaporation of the spreading solvent and avoid startup transients. This time corresponded to 240 or 710 oscillation cycles of the floor depending on the driving flow condition. After the initialization period, hollow-ceramic microsphere tracer particles (3M, ZEELAN Industries Inc., Z-Light, W-1800, sieved to between 75 and 150 µm) were sparsely scattered on the interface. Similar tracer-particle concentrations have been shown to exhibit no measurable effect on fluid flow by comparison between velocimetry via Brewster angle microscopy (with no tracer particles) and particle tracking velocimetry (PTV) measurement techniques.³¹ Immediately upon the addition of tracer particles, interfacial flow was imaged from above for 10 seconds, corresponding to ∼4 or 12 periods of oscillation for the two different flows. The resulting image sequences underwent PTV analysis and phase averaging over all captured oscillations to produce surface velocity distributions u^s(x, t) spanning the width of the channel for a full period of oscillation. These full-width surface velocity measurements were further averaged over three replicates, each consisting of an entirely distinct freshly-deposited monolayer. Between each experiment, the channel was washed sequentially with DI water, methanol, acetone, and then DI water again. Channel cleanliness was subsequently checked for expected performance with DI water and no monolayer.

2.3 Computational methods

Numerical simulations of the coupled bulk and interfacial flows were performed using COMSOL Multiphysics® version 6.1. Flow simulations were performed with a gas–liquid interface that remained flat and fixed in space. The 2D bulk fluid domain mesh consisted of 100 × 100 elements, spatially varying linearly in size from smallest in the corners to largest at the center, with a largest-to-smallest element ratio of 10. The 1D surface domain consisted of 250 equally-sized mesh elements. Verification of the COMSOL model was performed by comparing results for DPPC at c₀ = 1.5 mg m⁻² against a high-fidelity finite difference Fortran code.⁵³

The 2D bulk flow is modeled using the incompressible Navier–Stokes equations:

∇·u⃑ = 0 (1)

andwhere u⃑(x, y, t) = (u, v) is the velocity vector and g⃑ is the gravitational acceleration vector (0, −9.81 m s⁻²). A no-slip boundary condition (u = 0, v = 0) was prescribed for the two solid side walls (x = ±W/2). The floor of the channel (y = 0) was prescribed a tangentially-moving wall boundary condition (u = Usin(ωt), v = 0). The interfacial stress balance, coupling a flat interface to the bulk, was provided by the Boussinesq–Scriven interfacial stress balance:^1,38where u^s(x, t) is surface velocity, κ^s is surface dilatational viscosity, μ^s is surface shear viscosity, and σ(c) is surface tension as a function of monolayer concentration, c(x, t). The surface tension gradient, ∂σ/∂x, can be calculated using ∂σ/∂x = (∂σ/∂c) (∂c/∂x), where ∂σ/∂c is described by the equation of state. The second factor, ∂c/∂x, is computed using the 1D advection–diffusion equation for transport within the monolayer:^1,38,74where D^s is surface diffusivity.

The formulation of Navier–Stokes, Boussinesq–Scriven, and interfacial advection–diffusion forms a flow model with three interfacial material parameters: κ^s, μ^s, and D^s, along with the equation of state σ(c). The 1D interfacial flow of the oscillatory-driven channel (as expressed in eqn (3)) does not distinguish between the effects of κ^s and μ^s, resulting in an effective summed surface viscosity of κ^s + μ^s, used in this model as a single parameter. Physically, both surface viscosities must be positive, necessitating that κ^s + μ^s be greater than μ^s, which is well established.^28,30,31 D^s measured in systems at or near equilibrium as a function of concentration is termed equilibrium surface diffusivity D^s_e(c) and was obtained based on descriptions of D^s_e(c) as a function of monolayer phases.^49,50,75 Values of D^s_e(c) range from 10⁻² cm² s⁻¹ in gas phase and 10⁻⁶ cm² s⁻¹ in liquid phase to 10⁻¹⁰ cm² s⁻¹ in condensed phase. For simulations that did not use D^s_e(c), D^s was allowed to vary as necessary for the computed surface velocity to match the measurements. Computations utilized the Boussinesq approximation, where κ^s + μ^s and D^s only varied with c₀, and did not change with x, t. The Boussinesq approximation is necessary since parameter dependence on (x, t) is unknown and difficult to implement in practice.⁷⁶

Three interfacial models were evaluated through comparison to experimental measurements: model-i – an inviscid interface with equilibrium diffusivity (κ^s + μ^s = 0, D^s = D^s_e), model-ii – a viscous interface with equilibrium diffusivity (κ^s + μ^s > 0, D^s = D^s_e), and model-iii – a viscous interface with large effective diffusivity (κ^s + μ^s > 0, D^s ≫ D^s_e). Model-i required no fitting to experimental data as κ^s + μ^s (= 0) and D^s (= D^s_e(c₀)) were predetermined. For model-ii and model-iii, the computed spatio-temporal u^s(x, t) profiles were fit to experimental measurements by least squares optimization using κ^s + μ^s and D^s as the two fitting parameters (only κ^s + μ^s in model-ii and both κ^s + μ^s and D^s in model-iii). To assess the performance of each model, goodness of fit between measured and computed u^s(x, t) profiles was evaluated using the coefficient of determination R². Determination of the optimal R² for the two models that required fitting was first investigated with a sensitivity sweep in the κ^s + μ^s, D^s parameter space. The result of this coarse sweep was then used to provide initial values for the final fitting algorithm.

3 Results

This study evaluated insoluble monolayer behavior as a function of mean concentration c₀, first through general monolayer response using maximum surface velocity, then through spatio-temporal surface velocity measurements compared to numerical models. Fig. 3 shows the maximum surface velocity that was measured at the channel center, u^s_max|_x = 0 (normalized by floor speed U) as a function of c₀ for DPPC. u^s_max|_x = 0 data are essentially a coarse-grained measure of monolayer response at different c₀ to a given oscillatory driving flow. As c₀ increases, u^s_max|_x = 0 decreases monotonically from slightly more than 5% of the maximum floor velocity U at c₀ = 0 (clean air–water interface) to less than 1% at c₀ = 1.0 mg m⁻². After this initial monotonic decline with increasing c₀, u^s_max|_x = 0 plateaus for c₀ between 1.0 and 2.0 mg m⁻², followed by a further gradual decrease. Changes in u^s_max|_x = 0 trends coincide with monolayer phase, as depicted in the surface tension equation of state (Fig. 2a). Moreover, the qualitative resemblance between u^s_max|_x = 0 and monolayer compressibility C_s (Fig. 2b) demonstrates the same trend in the response regime. This phase-dependent response of u^s_max|_x = 0 was observed for both driving flow conditions and surfactant types, as well as at other channel locations such as the quarter widths (not shown for brevity). While the measured trends in u^s_max|_x = 0 demonstrate the dependence of hydrodynamic response on monolayer phase, analysis of full-field time-resolved u^s(x, t) measurements is necessary to evaluate the influence of κ^s + μ^s and D^s.

Fig. 3 Maximum surface velocity u^s_max (measured at the channel midpoint x = 0 and normalized by maximum floor speed U) versus mean DPPC monolayer concentration c₀ for the Re = 200 driving flow condition, with ±1 standard deviation error bars from three measurement replicates.

Spatio-temporal measurements of surface velocity u^s(x, t) made at the channel mid-span (z = 0) are used to evaluate the three interfacial models. Fig. 4 presents measurements of u^s profiles along with numerical solutions from all three models, each at quarter-period time points for the Re = 200 flow condition. LE (Fig. 4a) and SC (Fig. 4b) phases of DPPC (non-coexisting phases) are shown to examine the behavior of DPPC in regions on either side of recently presented⁵³ coexisting monolayer data. For all non-gaseous phases examined, model-i and model-ii perform equally poorly (the green model-ii curve falling atop the blue model-i) with best-fit u^s from these models falling far below measured magnitudes due to suppression of flow by Marangoni stresses computed based on equilibrium interfacial diffusivity D^s_e. Conversely, model-iii matches measured u^s magnitudes and shapes well by allowing effective interfacial diffusivity D^s to be much greater than D^s_e, significantly reducing Marangoni stresses. As quantified below, model-iii, using large effective D^s, is capable of fitting measured u^s in all non-gaseous monolayer phases for both surfactants under both driving flow conditions. No model fit measured u^s profiles in gaseous phases, where monolayer compressibility is extremely large and Marangoni stresses are minimal: model-i produces overly high u^s unrestrained by viscous stresses, while model-ii and model-iii over-attenuate u^s by viscous stresses.

Fig. 4 Normalized surface velocity u^s/U across channel width x/W, taken at quarter-period time points for DPPC at Re = 200. Curves show experimental measurements (black with ±1 standard deviation gray error shading from three replicates) and numerical solutions of the three interfacial models (color). (a) Shows an LE phase monolayer (c₀ = 1.1 mg m⁻²) and (b) shows an SC phase monolayer (c₀ = 2.2 mg m⁻²).

Model-iii, with finite κ^s + μ^s and large effective D^s, is the only model capable of capturing experimentally measured u^s(x, t) profiles in non-gaseous phases. The hydrodynamic effect of κ^s + μ^s and D^s is investigated through a sensitivity analysis that entails sweeps in the κ^s + μ^s, D^s parameter space over six orders of magnitude (10⁻³ < D^s < 10³ cm² s⁻¹ and 10⁻⁶ < κ^s + μ^s < 10⁰ g s⁻¹). Model-iii accuracy is quantified using R², with sensitivity results presented in Fig. 5 for gaseous (G) and non-gaseous (LE/LC) phase DPPC monolayers. The gas phase monolayer (c₀ = 0.5 mg m⁻²) is insensitive to D^s but sensitive to κ^s + μ^s, with a best matching value of ∼10⁻¹ g s⁻¹, see Fig. 5a. However, for gaseous phases, no combination of κ^s + μ^s and D^s provides an adequate fit to the measured u^s(x, t). Conversely, the non-gas phase monolayer presented (c₀ = 1.5 mg m⁻²) is significantly sensitive to D^s, indicated by the narrow horizontal region of large R² at D^s ≈ 10⁰ cm² s⁻¹, see Fig. 5b. Furthermore, κ^s + μ^s did not affect the solution until after the upper bound of κ^s + μ^s ≈ 10⁻¹ g s⁻¹, where fits become poor due to viscous quenching of interfacial flow.

Fig. 5 Coefficient of determination R² between experimentally measured surface velocity u^s(x, t) and computations (model-iii), presenting DPPC's sensitivity to surface viscosity (κ^s + μ^s) and surface diffusivity (D^s). (a) Shows a G phase (c₀ = 0.5 mg m⁻²) insensitivity to D^s and sensitivity to κ^s + μ^s, (b) shows an LE/LC phase (c₀ = 1.5 mg m⁻²) sensitivity to D^s with a κ^s + μ^s upper bound.

Nonlinear regression was performed following sensitivity analysis with κ^s + μ^s and D^s as fitting parameters. Goodness of fit for these regressions was again quantified using R², presented for DPPC monolayers at both Re in Fig. 6. As expected from earlier results, model-i does not fit experimentally measured u^s(x, t), with an average R² of 0 for both surfactants, both Re, and all c₀. Model-i produces u^s magnitudes much higher than measurements in gaseous phases due to the lack of interfacial viscous stresses, and much lower values in non-gaseous phases. Model-ii, though an improvement over model-i, still fits measurements poorly in gaseous phases with average R² = 0.42, under-predicting u^s magnitudes due to over-attenuation by viscous stresses. In non-gaseous phases, model-ii coincides with model-i and does not fit measured u^s, giving an average R² of 0, due to flow quenching by both interfacial viscous and Marangoni stresses. Model-iii fits similarly poorly to model-ii in gaseous phases, with an average R² of 0.45, with interfacial viscous stresses again producing low u^s unassisted by the large D^s that produce little effect by lowering the already minimal Marangoni stresses in gas-phase monolayers. Conversely, in non-gaseous phases, model-iii fits measured u^s magnitudes and shapes well with an average R² of 0.92, using nonzero κ^s + μ^s and large effective D^s to balance the effects of interfacial viscous and Marangoni stresses. The slight reduction of R² at c₀ ≥ 2.1 mg m⁻² in the figure for Re = 200 (see Fig. 6a) is likely due to measurement uncertainty associated with low flow magnitudes at such high c₀, as observed in Fig. 3.

Fig. 6 Coefficient of determination R² between experimentally measured surface velocity u^s(x, t) and the three interfacial models for DPPC monolayers at (a) Re = 200 and (b) Re = 600.

Fig. 7 presents model-iii fit results of effective D^s as a function of c₀ for DPPC at both Re, alongside literature-derived^50,75 equilibrium diffusivity D^s_e. In poorly-fit gaseous phases, solutions are insensitive to D^s due to the minimal contribution from Marangoni stresses. In well-fit non-gaseous phases, D^s increases from 6 to 10 orders of magnitude higher than D^s_e with increasing c₀, the largest changes occurring at c₀ corresponding to phase changes in the equation of state (see Fig. 2a). Additionally, κ^s + μ^s remains positive, well-fit solutions accepting the range μ^s < κ^s + μ^s ≈ 10⁻¹ g s⁻¹.^28,30‡

Fig. 7 Effective surface diffusivity D^s versus mean monolayer concentration c₀ obtained from fitting model-iii to DPPC measurements at each Re (gold and cyan), alongside literature estimates^49,50 of equilibrium diffusivity D^s_e (black).

An experimental and numerical analysis identical to that for DPPC was also performed for vitamin K1 monolayers. This protocol spanned c₀ of vitamin K1, ranging from the gas phase to just before monolayer collapse; see Fig. 8a. One distinction between monolayers of DPPC and vitamin K1 is that the latter exhibits two gaseous phases: G and G/L coexisting. D^s values as a function of c₀ produced from fitting model-iii to u^s profiles once again show large D^s in the well-fit non-gaseous phases. This large D^s represents an effective parameter required to produce accurate fits of the spatio-temporal surface velocity u^s(x, t) measurements using Newtonian and Fickian theory. Overall, vitamin K1 displayed trends similar to DPPC, with model-iii producing accurate fits in non-gaseous phases with large effective D^s presented in Fig. 8c, with a positive κ^s + μ^s. An observation distinct to vitamin K1 is that the stronger Re = 600 flow transitioned to large-diffusive behavior slightly earlier, at c₀ = 0.6 mg m⁻², than the Re = 200 flow case (though this early-transitioned case achieved only an R² of 0.62). For both surfactants at both Re, well-fit non-gaseous regions (at and above c₀ = 1.0 mg m⁻² for DPPC and c₀ = 0.8 mg m⁻² for vitamin K1) required similarly large diffusivities with D^s ≈ 10¹ cm² s⁻¹ (mean R² = 0.91).

Fig. 8 (a) Vitamin K1 surface tension equation of state σ(c), measured quasistatically via slow compression in a Langmuir trough, with an accompanying logistic fit§ used in numerical simulations. Experimental values are the average of three replicate measurements, with the gray error shade indicating one standard deviation. Labels indicate midpoints of monolayer phases: gas (G), coexisting gas/liquid (G/L), liquid (L), and collapse (C). (b) Coefficient of determination R² between experimentally measured surface velocity u^s(x, t) and model-iii for vitamin K1 monolayers at Re = 200 and Re = 600. (c) Effective surface diffusivity D^s versus mean monolayer concentration c₀ obtained from fitting model-iii to vitamin K1 measurements at each Re (gold and cyan), alongside literature estimates^49,50 of equilibrium diffusivity D^s_e (black). Area per molecule (A_m) is provided below the c axis for reference.

4 Discussion and conclusions

This study assessed the ability of hydrodynamic theory to capture phase-dependent behavior in insoluble monolayers of DPPC and vitamin K1 when forced far from equilibrium. The Boussinesq–Scriven surface model (Newtonian), with advection–diffusion (Fickian) for transport within the monolayer, was coupled to Navier–Stokes bulk flow. Surface diffusivity (D^s) and combined surface viscosity (κ^s + μ^s) were both studied through the performance of three interfacial models as compared to measured spatio-temporal surface velocity u^s(x, t), although the flow response was found to be insensitive to surface viscosity. Using equilibrium diffusivity D^s_e, a common practice for monolayer hydrodynamic models, neither inviscid (model-i) nor viscous (model-ii) interfaces could capture monolayer behavior. When D^s and κ^s + μ^s were allowed to vary freely as fitting parameters (model-iii), two key findings resulted: no fit was possible in gaseous phases, and non-gaseous phases could be fit by using extremely large effective D^s. Failure to fit in gaseous phases was likely due to the effect of extremely high compressibility (shown for DPPC in Fig. 2b), with surfactant concentration c(x, t) experiencing extreme variations due to the applied dynamic forcing. While there is a large variation in the computed c(x, t), this variation occurs in a monolayer phase that exhibits nearly no variation in surface tension, as shown in Fig. 2a and 8a, making Marangoni stresses weak. This extreme compressibility produces very large computed spatio-temporal variation in c, such that D^s and κ^s + μ^s, which are functions of c, likely violate the assumption of the Boussinesq approximation, where parameters do not vary in x, t. This observation suggests that equations of state for viscosity and diffusivity may be essential for capturing monolayer behavior in highly compressible regimes, circumventing the need for the Boussinesq approximation.⁷⁶

In non-gaseous phases, computations were able to provide a good fit to measured u^s(x, t) by using extremely large effective D^s of ∼10¹ cm² s⁻¹, displayed for DPPC in Fig. 7 and for vitamin K1 in Fig. 8c. These large D^s values represent an effective parameter, one that lessens surface tension gradients, to produce the balance required for fitting. In interfacial advection–diffusion, increasing D^s (via fitting) produces the necessary effect in the system, preventing large gradients in c through diffusion, by rapid homogenization. Consequently, within the Boussinesq–Scriven interfacial stress model, the computed Marangoni stress due to surface tension gradients was reduced by the large effective D^s, with the interface more rapidly responding in stress to applied strain. Such a near-instantaneous response in stress to strain is characteristic of linear elasticity. This implies that non-Newtonian and/or non-Fickian behavior, where surface viscous stress or surface diffusivity are functions of strain rate and/or strain, are a plausible explanation for such large effective D^s. Additionally, while the model was insensitive to κ_s + μ_s, good fits were possible for all κ^s + μ^s < 10⁻¹ g s⁻¹. These well-fit solutions accepted the range μ^s < κ^s + μ^s ≈ 10⁻¹ g s⁻¹.^28,30

The findings presented here add to the body of evidence showing that surfactant monolayers rapidly forced far from equilibrium require D^s seven to ten orders of magnitude above equilibrium values. Large effective D^s may also be implicated in the reported anomalies in κ^s when D^s_e is used in modeling, indicating that large parameter values are compensating for inadequate modeling and highlighting the limitations of Newtonian and Fickian interfacial models. The phase-dependent trends in D^s and the importance of compressibility demonstrated in this work motivate investigation into complex interfacial models.

Author contributions

Tyler J. Mucci: writing – original draft, writing – review & editing, methodology, formal analysis, data curation, software, investigation, validation, project administration. Joe A. Adam: writing – review & editing, project administration, supervision. Amir H. Hirsa: writing – review & editing, funding acquisition, project administration, supervision.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this study's findings and the codes used for simulation and analysis are available from the author upon reasonable request.

Acknowledgements

This work was supported by NSF grants 2204081 and 2323020, and NASA grant 80NSSC25K7630.

References

1. D. Edwards, H. Brenner and D. Wasan, Interfacial Transport Processes & Rheology, Butterworths-Heinemann, Stoenheam, MA, USA, 1991.
2. A. Adamson and A. Gast, Physical Chemistry of Surfaces, Wiley, New York, NY, USA, 1997.
3. R. J. E. Farn, Chemistry and Technology of Surfactants, Blackwell Publishing, Oxford, UK, 2006.
4. D. Myers, Surfactant Science and Technology, John Wiley & Sons, Hoboken, NJ, USA, 2020.
5. H. Manikantan and T. M. Squires, J. Fluid Mech., 2020, 892, P1.
6. J. Yalim, J. M. Lopez, S. R. Griffin, J. A. Adam, K. D. Brown, P. M. McMackin and A. H. Hirsa, Physica D, 2024, 463, 134150.
7. J. M. Lopez and A. H. Hirsa, Eur. Phys. J.:Spec. Top., 2026, 1–10.
8. E. Hermans, M. Saad Bhamla, P. Kao, G. G. Fuller and J. Vermant, Soft Matter, 2015, 11, 8048–8057.
9. M. S. Bhamla, C. Chai, M. A. Àlvarez-Valenzuela, J. Tajuelo and G. G. Fuller, PLoS One, 2017, 12, e0175753.
10. S. Narayan, A. E. Metaxas, R. Bachnak, T. Neumiller and C. S. Dutcher, Curr. Opin. Colloid Interface Sci., 2020, 50, 101385.
11. C. R. Anthony, H. Wee, V. Garg, S. S. Thete, P. M. Kamat, B. W. Wagoner, E. D. Wilkes, P. K. Notz, A. U. Chen, R. Suryo, K. Sambath, J. C. Panditaratne, Y.-C. Liao and O. A. Basaran, Annu. Rev. Fluid Mech., 2023, 55, 707–747.
12. K. Y. Chan and A. Borhan, J. Colloid Interface Sci., 2005, 287, 233–248.
13. A. H. Kumar, H. Wee, N. K. Dhanwani and O. A. Basaran, J. Fluid Mech., 2026, 1030, R1.
14. O. E. Jensen, J. Fluid Mech., 1995, 293, 349–378.
15. T. Bickel and F. Detcheverry, Phys. Rev. E, 2022, 106, 045107.
16. A. S. Luviano, J. Campos-Terán, D. Langevin, R. Castillo and G. Espinosa, Langmuir, 2019, 35, 16734–16744.
17. H. Manikantan and T. M. Squires, Phys. Rev. Fluids, 2017, 2, 023301.
18. S. He, M. Afshang, M. Caggioni, S. Lindberg and K. M. Schultz, Langmuir, 2023, 39, 12346–12356.
19. P. Mikolka, et al., Lung, 2024, 202, 299–315.
20. S. J. Pogorzelski and A. D. Kogut, J. Sea Res., 2003, 49, 347–356.
21. A. T. Florence and J. M. Gillan, Pestic. Sci., 1975, 6, 429–439.
22. M. Rump, C. Diddens, U. Sen, M. Versluis, D. Lohse and T. Segers, Phys. Rev. Appl., 2025, 23, 024076.
23. Y. Y. Zuo and F. Possmayer, J. Appl. Physiol., 2007, 102, 1733–1734.
24. H. Zhang, Y. E. Wang, Q. Fan and Y. Y. Zuo, Langmuir, 2011, 27, 8351–8358.
25. X. Liu and J. H. Duncan, Nature, 2003, 421, 520–523.
26. H. E. Gaub and H. M. McConnell, J. Phys. Chem., 1986, 90, 6830–6832.
27. J. T. Petkov, K. D. Danov, N. D. Denkov, R. Aust and F. Durst, Langmuir, 1996, 12, 2650–2653.
28. A. Hirsa, J. Lopez and R. Miraghaie, J. Fluid Mech., 2002, 470, 135–149.
29. P. Erni, P. Fischer and E. J. Windhab, 3rd International Symposium on Food Rheology and Structure, 2003, 411–415.
30. A. Sadoughi, J. M. Lopez and A. H. Hirsa, Phys. Fluids, 2013, 25, 032107.
31. A. Raghunandan, A. H. Hirsa, P. T. Underhill and J. M. Lopez, Phys. Rev. Lett., 2018, 121, 164502.
32. A. J. James and J. Lowengrub, J. Comput. Phys., 2004, 201, 685–722.
33. L. M. C. Sagis, Rev. Mod. Phys., 2011, 83, 1367–1403.
34. A. H. Hirsa and J. M. Lopez, Fluids, 2021, 6, 198.
35. N. O. Jaensson, P. D. Anderson and J. Vermant, J. Non-Newtonian Fluid Mech., 2021, 290, 104507.
36. J. M. Lopez, M. J. Vogel and A. H. Hirsa, Phys. Rev. E, 2004, 70, 056308.
37. P. Erni, P. Fischer and E. J. Windhab, Langmuir, 2005, 21, 10555–10563.
38. J. Slattery, L. Sagis and E. Oh, Interfacial Transport Phenomena, Springer, New York, NY, USA, 2007.
39. M. Pourali, M. Kröger, J. Vermant, P. D. Anderson and N. O. Jaensson, Phys. Fluids, 2021, 33, 062103.
40. E. Guzmán, A. Maestro, C. Carbone, F. Ortega and R. G. Rubio, Fluids, 2022, 7, 335.
41. C. O. Ciutara, S. Barman, S. Iasella, B. Huang and J. A. Zasadzinski, J. Colloid Interface Sci., 2023, 629, 125–135.
42. A. García-Gutiérrez, F. Parra, J. A. Mayo and C. Rubio, Theor. Comput. Fluid Dyn., 2025, 39, 44.
43. W. A. A. Sudjarwo and J. L. Toca-Herrera, Langmuir, 2025, 41, 16128–16138.
44. J. Giermanska-Kahn, F. Monroy and D. Langevin, Phys. Rev. E:Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1999, 60, 7163–7173.
45. V. E. Cuenca, H. Martinelli, M. D. L. A. Ramirez, H. A. Ritacco, P. Andreozzi and S. E. Moya, Langmuir, 2019, 35, 14300–14309.
46. G. L. Gaines, Insoluble Monolayers at Liquid-gas Interfaces, Wiley-Interscience, 1966.
47. A. B. Serfis and R. Katzenberger, Colloids Surf., A, 1998, 138, 91–95.
48. S.-M. Argyri, L. Evenas and R. Bordes, J. Colloid Interface Sci., 2023, 640, 637–646.
49. R. Peters and K. Beck, Proc. Natl. Acad. Sci. U. S. A., 1983, 80, 7183–7187.
50. A. Shmyrov and A. Mizev, Langmuir, 2019, 35, 14180–14187.
51. G. Homsy and E. Meiburg, J. Fluid Mech., 1984, 139, 443–459.
52. L. R. Arriaga, I. López-Montero, J. Ignés-Mullol and F. Monroy, J. Phys. Chem. B, 2010, 114, 4509–4520.
53. T. J. Mucci, B. L. Liu, J. A. Adam, A. H. Hirsa, J. Yalim and J. M. Lopez, Phys. Rev. E, 2025, 111, L013501.
54. M. Agassandian and R. K. Mallampalli, Biochim. Biophys. Acta., 2013, 1831, 612–625.
55. A. Stachowicz-Kusnierz, T. Seidler, E. Rogalska, J. Korchowiec and B. Korchowiec, Chemosphere, 2020, 240, 124850.
56. E. Guzmán, Coatings, 2022, 12, 277.
57. L. Scriven, Chem. Eng. Sci., 1960, 12, 98–108.
58. P. Griffiths, D. Xu and L. Davoust, J. Non-Newtonian Fluid Mech., 2025, 346, 105498.
59. C.-H. Chang and E. I. Franses, Colloids Surf., A, 1995, 100, 1–45.
60. Q. Zhao, W. Ren and Z. Zhang, Comput. Meth. Appl. Mech. Eng., 2021, 385, 114033.
61. H. Maru and D. Wasan, Chem. Eng. Sci., 1979, 34, 1295–1307.
62. D. Sharpe and J. C. Earnshaw, J. Chem. Phys., 1997, 107, 7493–7501.
63. F. Monroy, J. G. Kahn and D. Langevin, Colloids Surf., A, 1998, 143, 251–260.
64. L. M. Sagis, Adv. Colloid Interface Sci., 2014, 206, 328–343.
65. D. Edwards, J. Polym. Sci. B, 1996, 34, 981–997.
66. J. Ferreira, M. Grassi, E. Gudino and P. de Oliveira, Appl. Math. Model., 2015, 39, 194–204.
67. T. Mucci, PhD thesis, Rensselaer Polytechnic Institute, 2025.
68. D. Vollhardt and V. Fainerman, Adv. Colloid Interface Sci., 2006, 127, 83–97.
69. D. I. Krylov and V. V. Arslanov, Colloids Surf., A, 2025, 717, 136771.
70. M. J. Vogel, A. H. Hirsa and J. M. Lopez, J. Fluid Mech., 2003, 478, 197–226.
71. H. M. Blackburn and J. Lopez, J. Fluid Mech., 2003, 497, 289–317.
72. F. Marques, J. M. Lopez and H. M. Blackburn, Phys. D, 2004, 189, 247–276.
73. J. J. Leung, A. H. Hirsa, H. M. Blackburn, F. Marques and J. M. Lopez, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2005, 71, 026305.
74. J. M. Lopez and A. H. Hirsa, J. Colloid Interface Sci., 2001, 242, 1–5.
75. S. Peace, R. Richards and N. Williams, Langmuir, 1998, 14, 667–678.
76. J. Lopez and A. Hirsa, J. Colloid Interface Sci., 1998, 206, 231–239.

---

Footnotes

† Parameters for DPPC: (a₁, a₂) = (67.31, 4.80) mN m⁻¹, (b₁, b₂) = (9.06, 22.14), (c₁, c₂) = (2.63, 1.13) mg m⁻², (d₁, d₂) = (0, 0) mN m⁻¹, (e₁, e₂) = (1.889, 0.578) mN m⁻¹.

‡ Well-fit solutions in non-gaseous phases insensitive to κ^s + μ^s used the initial value of 0.0048 g s⁻¹.

§ Parameters for vitamin K1: (a₁, a₂) = (10.39, 61.76) mN m⁻¹, (b₁, b₂) = (12.95, 74.86), (c₁, c₂) = (1.76, 0.88) mg m⁻², (d₁, d₂) = (2.50, 51.30) mN m⁻¹, (e₁, e₂) = (235.81, 0.041).

---