Coupling Vortical Bulk Flows to the Air–Water Interface: From Putting Oil on Troubled Waters to Surfactants on Protein Solutions

: The air–water interface in ﬂowing systems remains a challenge to model, even in cases where the interface is essentially ﬂat. This is because even though each side is governed by the Navier–Stokes equations, the stress balance which provides the boundary conditions for the equations involves properties associated with surfactants that are inevitably present at the air–water interface. Aside from challenges in measuring interfacial properties, either intrinsic or ﬂow-dependent, the two-way coupling of bulk and interfacial ﬂows is non-trivial, even for very simple ﬂow geometries. Here, we present an overview of the physics associated with surfactant monolayers of ﬂowing liquid and describe how the monolayer affects the bulk ﬂow and how the monolayer is transported and deformed by the bulk ﬂow. The emphasis is primarily on cylindrical ﬂow geometries, and both Newtonian and non-Newtonian interfacial responses are considered. We consider interfacial ﬂows that are solenoidal as well as those where the surface velocity is not divergence free.


Introduction
In this paper, which is dedicated to the late Prof. William W. Willmarth, we provide an overview of the hydrodynamics associated with surfactants (surface-active materials), due to their omnipresence in both natural and manmade systems. Our emphasis is on theoretical modeling of vortical flows with significant inertia, inspired by the work of Willmarth in free-surface turbulence. Our perspective was also shaped by Willmarth, who advised the PhD work of the first author, and strongly believed in carefully conducted experiments to discover new phenomena and closely coupled with theory/numerics. In that spirit, here we present the theoretical foundation for two-way coupling of inertial bulk flow with surfactant-covered interfaces. The stress-strain constitutive relation for Newtonian interfaces is first discussed and how it may be generalized. This work is a culmination of more than 25 years of collaboration between an experimenter in free-surface flows and a theoretical hydrodynamist.
Interfacial hydrodynamics can drastically alter the transfer of mass, momentum, and energy across the surface of a liquid due to the coupling to the bulk hydrodynamics. The importance of interfacial hydrodynamics is evident everywhere we look in nature, ranging from small scale systems such as the alveoli in the lungs [1] to CO 2 transfer between the atmosphere and the oceans [2][3][4][5]. In manufacturing, interfacial hydrodynamics impacts everything from wet processing of micro-electronics to the production of polyurethane foam [6,7]. A recent development in the pharmaceutical manufacturing industry is the increased utilization of interfacial processing [8]. Yet, it is extraordinary how little is understood about interfacial hydrodynamics, how to predict it, and ultimately, how to control it.
Our understanding of interfacial hydrodynamics has come a long way from just knowing that oil can calm troubled waters, but we still have a long way to go. A major challenge of our time is manufacturing of pharmaceuticals, including monoclonal antibodies, at a global scale. For more than 50 years, pharmaceuticals have been produced using batch manufacturing which is a multi-step process involving synthesis, crystallization, blending, granulation and tableting. In between each stage, the products are tested, shipped and stored from one batch processing location to another. More recently, there has been a slow shift to continuous manufacturing, where the pharmaceuticals are produced simultaneously on a single process line, increasing both cost efficiency and product quality [9]. Among the challenges in designing these new processors for each manufacturing stage is how to mitigate adverse effects from hydrodynamic and interfacial stresses [10][11][12][13]. For example, in bioreactors where microbial factories such as E. coli are used to over-express therapeutic proteins, such as insulin via recombinant technologies, one must limit the range of shear stress imposed on the micro-organisms. After the protein is extracted, it must be purified and utilized at large concentrations-such as in injectable monoclonal antibodies. Proteins can denature and aggregate as a result of large shear and/or exposure to hydrophobic interfaces. Surfactant additives are utilized in various stages of pharmaceutical manufacturing processes in order to mitigate or at least minimize protein aggregation. The current state-of-the-art in pharmaceutical manufacturing resembles the aeronautics industry 30 years ago, which then relied on wind tunnel tests, as opposed to predictive computational fluid dynamics today.
Many of the difficulties in predicting transport processes result from inadequacies in our understanding of the coupled hydrodynamics between bulk flow, interfacial flow, and chemical species that contribute to the interfacial stress balance. Interfacial hydrodynamics entails many different facets. The problem of complex bulk flow with a highly curved free surface has been comprehensively treated in [14], where surface viscosity is not considered. Here, our interest is in complex bulk flows with viscous interfaces, such as large mixers employed during drug blending operations which are relevant in the pharmaceutical industry. In order to make this problem tractable, we begin by focusing on essentially flat interfaces. Thus, we need physically relevant constitutive relations for the interface as well as proper coupling between interfacial and bulk stresses.
The observations that liquid surfaces often have a non-Newtonian response and that their flow is not decoupled from the bulk flow motivate the need to solve the problem of interfacial flow coupled to bulk flow [19,25,30,35,36]. Such a capability is needed to determine material properties from experimental measurements of non-Newtonian interfaces, and finally to enable predictive modeling. Appropriate constitutive equations and intrinsic (material) properties are critical because one cannot use apparent properties measured at a given thermodynamic state (e.g. surface temperature and surfactant concentration) under one set of flow conditions to make predictions at the same thermodynamic state but under different flow conditions. Even if the intrinsic properties are determined under conditions where coupling to the bulk is minimal, it is important to understand how that interface couples to bulk flows in other conditions. The lack of such a capability is hindering both scientific progress in understanding nature and technological advances in health care and other industries significantly impacted by interfacial processes.
Recently, ref. [37] provided a review of interfacial rheology, where rheometric flows with low inertia are used to measure intrinsic and apparent properties of interfaces. These include various phospholipids and proteins. Two-way coupling of bulk and interfacial flow has been recently been described for cases where inertia is negligible [38]. In this paper, we focus on flows where inertia is important. Further, we restrict the discussion to insoluble surfactants. Soluble surfactants introduce additional phenomena, which are not addressed here, associated with their transport between the interface and the bulk [39][40][41][42]. Nevertheless, irrespective of whether the surfactants are soluble or insoluble, the interfacial velocity need not be divergence-free, even if the interface does not change shape or size, particularly when the bulk flow has non-negligible inertia [31,33]. These are central issues considered in this paper.

Interfacial Stress Balance
The coupling between the interfacial and bulk flows occurs because of a stress balance. The general stress balance is given in [43]. Ignoring stresses from the air, this balance relates the liquid stress to the interfacial stress. We begin by considering the Boussinesq-Scriven surface model [44], which is applicable to Newtonian interfaces: where u s is the surface velocity vector, σ is the surface tension, µ s is the surface shear viscosity, κ s is the surface dilatational viscosity, T s is the surface stress tensor, I s is surface projection tensor, and the surface rate of deformation tensor is The stress at the interface is Usually, µ s and κ s are taken as constants (often zero). Here, we treat these as functions of the surface concentration c, as is σ. The surface concentration c is determined by a surface advection-diffusion equation where d s is the surface diffusivity. For a flat interface, the normal stress balance gives the normal velocity at the interface w s = 0. The azimuthal and radial interfacial velocity components are determined by the tangential stress balances, given by the azimuthal and radial components of (3). These have the form where µ is the shear viscosity of the bulk fluid, and f 1 and f 2 are functions that depend on everything (σ, µ s , κ s , c, u s , v s ); detailed derivations of these are presented in [45]. In the past, we have considered various regimes in which these equations are greatly simplified [46]. These restrictions need to be relaxed in order to simulate complex bulk flows coupled to complex interfaces. In the following section, we outline how to generalize this formulation to non-Newtonian interfaces (including shear-rate dependency), in the axisymmetric regime for ease of exposition.
In purely shearing interfacial flows (with no interfacial dilation or compression), many simple monomolecular films exhibit a Newtonian response under varying flow conditions and different flow geometries, yielding a consistent measurement of their surface shear viscosity [18,47]. These interfaces also exhibit well-understood Marangoni stresses [17,48]. Measuring the properties of Newtonian interfaces is now a mature field, whereas that of complex (non-Newtonian) interfaces is not as well developed [49]. Importantly, interfaces with large macro-molecules, such as proteins or many surfactants at large surface pressures (surface pressure Π is the difference in surface tension between a clean interface and the interface covered by the monolayer), do not exhibit a Newtonian response. Furthermore, often the interface is not only sheared but is subjected to a combination of both shear and dilation, which introduces Marangoni effects.
It is widely understood that the effects of surface shear viscosity µ s are dominant at small scales. However, it is not widely appreciated that surface shear viscosity can dramatically alter the entire flow field even at large scales. Figure 1 shows experiments conducted with a length scale of approximately 10 cm, with Reynolds number ∼10 3 , where a columnar vortex in the bulk meets the air-water interface [50]. For an interface with an oleyl alcohol monolayer, with µ s ∼ 0, the flow looks identical to the case with a clean interface, while with a stearic acid monolayer on the interface, which has a non-zero µ s , the bulk flow is altered. The photographs were taken in the meridional (vertical) plane via laser-induced fluorescence. It is remarkable that a monomolecular film (thickness of about one nanometer) can dramatically affect the bulk flow, generating a toroidal vortex of diameter about 10 cm. The resultant bulk flow is practically indistinguishable from the case where a rigid solid lid is placed at the free surface immediately after the columnar vortex is generated [51].  Recently, a containerless biochemical reactor was introduced for the microgravity environment where containment is achieved by surface tension and shear is propagated primarily by the action of the surface shear viscosity [52][53][54][55]. Such a ring-sheared drop, shown schematically in Figure 2a, is driven by the rotation of a thin contact ring at some latitude in the top hemisphere along with a stationary ring in the bottom hemisphere. The ring-sheared drop was originally designed for the study of amyloid fibril formation aboard the International Space Station. Aside from mitigating settlement of aggregates in a microgravity environment, the ring-sheared drop minimizes complications associated with wall-nucleation during amyloidogenesis. Amyloid fibrils in flowing systems are widely studied since they are implicated in neurodegenerative diseases, including Parkinson's and Alzheimer's disease [56].
Numerical simulations of the flow in the ring-sheared drop, presented in Figure 2b,c, illustrate the effects of the surface shear viscosity on the vortex lines. For very small values of surface shear viscosity µ s (quoted non-dimensionally as a Boussinesq number Bq = µ s /(µR) = 10 −4 , where µ is the dynamic viscosity of the bulk drop liquid and R is the radius of the rings), the drop is essentially in solid-body rotation at the average angular velocity of the two rings; the vortex lines are essentially aligned with the rotation axis. On the other hand, for a large value of surface shear viscosity, corresponding to Bq = 10 2 , the vortex lines show significant bending near the interface. A key observation is that surface viscosity does not directly affect the incidence angle of vortex lines at the interface, but rather it leads to a decelaration of the azimuthal flow component near the interface resulting in axial bending of the vortex lines, which in turn drives a secondary meridional flow [52]. The results shown in the figure are for Reynolds numbers Re = ΩR 2 /ν = 10 2 and 10 3 , where ν is the drop kinematic viscosity, further emphasizing that the effects of surface shear viscosity are not only manifested at small scales or low inertia.

Modeling the Non-Newtonian Interfacial Response
The simplest linear constitutive equation corresponding to a Newtonian interface is the Boussinesq-Scriven surface model (1). In regimes were the flow is axisymmetric, then the interfacial flow is purely shearing, the effects of surface dilatational viscosity are absent, Marangoni stress keeps the radial velocity zero, and the surface shear viscosity µ s is the only surface property involved in the surface stress balance. In [36], the constitutive relation was generalized in order to model shear-thinning interfacial effects due to the presence of large macromolecules, such as proteins or densely-packed condensed-phase monolayers of phospholipids, by replacing µ s with µ s eff , which is the effective surface shear viscosity, a function of the magnitude of the local interfacial shear rateγ. The choice for the functional form of µ s eff was motivated by the multiple reports of shear-thinning in DPPC films at high surface packing (corresponding to large surface pressure Π 30 mN/m) and its Newtonian response at low packing (Π < 15 mN/m) [25,27,57], and inspired by the Sisko model in bulk (3D) rheology [58]: where µ s ∞ is the surface shear viscosity at very large shear rates,K is the consistency index, and n is the power-law index. µ s ∞ ,K, and n are material properties of the interface, and as such are only functions of the thermodynamic state of the material at the interface, i.e., the surface concentration c of the monolayer at a given temperature. DPPC is widely studied since it is the main constituent of lung surfactants, and also plays an important role in the bilayer that forms mammalian cell membranes.

Two-Way Coupled Flows
In order to gain a fundamental understanding of the two-way coupling between bulk and interfacial flows, it is desirable to keep the geometry as simple as possible, while allowing for a fully 3D unsteady bulk flow coupled to a nominally flat 2D unsteady and potentially non-Newtonian interface. Such a system is shown schematically in Figure 3, consisting of an open stationary cylinder filled to the brim with water, with a monolayer on the surface. The flow is driven by the constant rotation of the floor. When the floor is driven at a relatively low rotation rate Ω, the bulk flow is axisymmetric. Fast rotation of the floor leads to fluid inertia causing instabilities breaking the symmetry, and the coupling between complex bulk flow and interface can be studied. This open cylinder flow is an idealization of many chemical and biological reactors and mixers, with the impeller/rotor replaced by the rotating floor [59,60]. The simplification afforded by this generic flow system allows us to focus on coupled hydrodynamics between the bulk and the interface without complications associated with the impeller/rotor, such as cavitation and 3D intricacies of the blade design. This flow geometry has been attracting some broad attention [61][62][63][64][65].
In Figure 3, two different monolayers were used. One was a monomolecular film of vitamin K 1 , which is insoluble in water. The other was a single layer of micron-scale plastic particles (polystyrene) treated to be surface bound and scattered on the surface. Vitamin K 1 forms a well-behaved Langmuir film on water, with no noticeable hysteresis (surface tension vs. concentration measured in a Langmuir trough is identical during compressions and expansions), and no measurable surface viscosities. It only exhibits elasticity (Marangoni effect), which can be captured accurately with computations [17], or even simple analytic modeling [48]. The surface-bound polystyrene monolayer behaves as a fluid marker, with little or no elasticity until the particles are closely packed. The bulk flow in the open cylinder system (Figure 3) is governed by the rotation rate of the floor, Ω, the cylinder radius, R, and the kinematic viscosity of the bulk fluid (nominally water), ν. These combine to give the Reynolds number Re = ΩR 2 /ν. The cylinder height-to-radius aspect ratio, H/R, is also important-it determines how much of the shear introduced by the rotating floor reaches the surface. Then, of course, what material is on the surface also impacts the coupled flow, as was illustrated in a related problem in Figure 1. Regardless of whether it is clean or not, the deformability of the surface is governed by the Froude number Fr = Ω 2 R 2 /gH, where g is the gravitational acceleration. For the experiment shown in Figure 3, Fr ≈ 0.03, and the largest surface deformation was measured to be 0.006R and the largest slope of the surface was 0.008 [20]. For flows with larger aspect ratio H/R, Fr is even smaller and surface deformations will remain negligible. This flow with a nominally flat and clean interface has been studied experimentally and numerically [67], showing excellent agreement capturing the 3D instability. It is a bulk flow instability which is manifest at the interface. Figure 4 shows a direct comparison between the computed axial vorticity from 3D Navier-Stokes simulations and the experimentally measured axial vorticity using PIV, at mid-depth, together with isosurface plots of the computed perturbation velocity components. If the interface is not perfectly clean-which is generally the case with water no matter how pure the source or the container-then the flow can be drastically different. Surfactants generally reduce the surface tension. Gradients in the surfactant concentration result in surface tension gradients, i.e., Marangoni stresses. For a flat interface, there are no curvature effects and so the absolute value of surface tension is not dynamically important, but its gradients are and these must be balanced by viscous stresses in the bulk fluid evaluated at the interface. If the surfactant is viscous, it will also impart a viscous resistance on the interfacial flow, which must also be balanced by viscous stresses in the bulk fluid evaluated at the interface.
In the past we had focused on interfacial property measurements and so kept the flows as simple as possible. For example, to quantify the Marangoni stress or the surface shear viscosity, measurements and computations were done in an annular system, where the flow was driven fast enough to get a good signal-to-noise ratio, but not so fast as to break the axisymmetry [17,18]. We have also shown through experiments and computations that in oscillatory-driven axisymmetric flows, the phase lag between the oscillatory driving mechanism and interfacial velocity response can be fully accounted for by the coupling between interfacial and bulk flows, an effect that can be misinterpreted as a non-Newtonian interfacial response [29].
We have also studied non-Newtonian interfaces using axisymmetric flow to determine shear-thinning responses under steady conditions [36]. The non-Newtonian interface model (Section 3) was implemented and validated in the knife-edge flow geometry, consisting of a stationary open cylinder also filled with water to the brim, with a monolayer on the surface. A schematic is shown in Figure 5. The monolayer is sheared by the constant rotation of a circular knife edge of radius a and thickness , which just touches the interface. The model and experiments were able to access many regimes of non-Newtonian behavior. Interfacial azimuthal velocity profiles of DPPC at surface pressure Π = 40 dyne/cm over a range of Re and R/a in the knife-edge flow geometry are shown in Figure 6. These were measured using particle tracking velocimetry. We found good agreement between model predictions and the measurements over two decades of Re and for flows in different aspect ratio containers. Figure 6 shows that the measured interfacial velocity profiles (symbols) at all Re and R/a are clearly slower than the predicted Newtonian response (cyan curves). With the coupling to the bulk flow fully accounted for, this departure clearly shows the non-Newtonian nature of DPPC at this high Π. In contrast, the measured velocity profiles are in excellent agreement with predictions made by treating the interface as shear-thinning with a unique set of material properties [36]. Perhaps more importantly, the results show that the effective surface shear viscosity is modeled consistently over a range of six decades in shear rate variation, as shown in Figure 7. Those results also show that we were able to recover measurements by Hermans and Vermant [57] for the same monolayer but in a different flow geometry. In [36], the viscous response of the entire interfacial flow field was measured and matched against computed interfacial velocity using our non-Newtonian interfacial model to determine interfacial properties, in contrast to the typical practice of using a local estimate of the imposed shear. This allowed us to explore global/macroscopic effects of the imposed shear to quantify interfacial properties.    Figure 6. The black curve µ s eff computed from the non-Newtonian model [36], and the green diamond symbols are data from [57] taken in a double-wall ring apparatus with DPPC at Π = 40 mN/m. (Reproduced with permission from [36], published by American Physical Society 2018).

Marangoni Effects
In the previous section, interfacial flows were considered where the monolayer concentration remained uniform, and hence the Maragoni stress was zero, in order to isolate the effects of interfacial shear viscosity. In this section, we consider interfacial flows that drive surfactant concentration gradients and hence have non-zero Marangoni stress. Monolayers of vitamin K 1 are particularly useful to ellucidate these as they have essentially zero surface shear viscosity.
The deep-channel surface viscometer, schematically depicted in Figure 8, provides a sensitive means for determination of the surface shear viscosity when the radial component of the surface flow is negligible, i.e., Re = ΩR 2 o /ν is sufficiently small [68]. At large Re, where inertia is dominant, the channel can be utilized for the study of Marangoni elasticity, namely the balance between the surface tension gradient and liquid stress evaluated at the surface. A strong secondary flow is generated by the rotation of the floor, driving the fluid out radially in the rotating floor boundary layer. This radial boundary layer flow is turned up along the outer cylinder, which then turns inward at the surface. Measurements of the radially inward flow at the free surface in a relatively large annular channel, with outer radius R o ≈ 10 cm and inner radius R i ≈ 8 cm, are shown in Figure 9a. The initial (uniform) concentration of a vitamin K 1 monolayer was varied. At sufficiently small initial concentrations, the radial inward flow near the surface is strong enough to completely clear the outer region of the flow and compress the monolayer toward the inner cylinder. Numerical simulations using the equation-of-state, measured in a Langmuir trough, provide surprisingly good agreement with the experiments (Figure 9b), especially in light of the fact that the capillary number, Ca = µΩr 0 /σ, is very small for the experiments with water at that scale. The small Ca makes the results very sensitive to the details in the equation-of-state. The main discrepancy is that with an initial concentration c 0 0.8 mg/m 2 , the simulations predict that the monolayer cannot be completely cleared near the outer cylinder, whereas in the experiments that limit is c 0 ∼ 1 mg/m 2 .

Surface Dilatational Viscosity and Non-Equilibrium Monolayer State
Surface tension is a well-defined equilibrium quantity. It is commonly measured quasi-statically using a Langmuir trough, and is used successfully in static equilibrium problems. Shear viscosity has also been consistently measured by a variety of techniques, all of them implemented at steady state [18,47,68]. Oscillatory sheared monolayers can also be treated via equilibrium interfacial rheology provided that the bulk flow inertia is accounted for in the interfacial stress balance, and the monolayer concentration remains essentially uniform [29].
The situation with dilatational viscosity κ s is much more complicated; its measurement using techniques with different time scales for the surface strain varies by as much as a factor of 10 5 for a given surfactant [69]. In a compressing or dilating interfacial flow, the time scale for the flow process is usually different from the time scales associated with monolayer equilibration. Thus, resistance to compression or dilation as predicted using equilibrium surface tension gradients may not be reconciled with the resultant stress. In fact, there have been several reports of apparent negative surface dilatational viscosity measurements [33,[70][71][72][73][74][75][76][77][78]. Negative viscosity is physically inconsistent as it violates the second law of thermodynamics, and some have conjectured this anomaly to be due to monolayer phase behavior [33,75], but it remains a hotly debated issue [79][80][81].
Usually, the equation-of-state is measured by spreading a uniform monolayer on a Langmuir trough, allowing it to equilibriate for a long time and then subjecting it to a quasi-static compression and/or expansion and measuring the surface tension (using a Wilhelmy plate) as a function of monolayer concentration, giving the equation-of-state. However, if the monolayer is first forced, for example by shearing, then the resultant equation-of-state is found to be different because the monolayer phase morphology is different [82]. The point is that for continuously compressed and expanded monolayers, finding the equation-of-state for a freshly spread static monolayer is not what is needed. It is the thermodynamic properties (surface tension and surface viscosities) of the monolayer in its morphological state under oscillatory dilatational flow that are needed.

Conclusions
Interfacial systems that have focused on modeling and measuring their rheological properties [37] have been such that the bulk flow is generally in the inertialess Stokes flow regime. In those cases, the role of the bulk flow is merely to support the interfacial film.
The challenge is what happens when the bulk flows are hydrodynamically important, e.g., when inertia is not negligible. Although the bulk flow is generally still incompressible, the interfacial flow is not divergence-free. Sagis [83] discusses how to treat dilating interfaces, but the challenge remains how to couple these interfacial flows to hydrodynamically important and inertia-dominated bulk flows. When one is growing a micro-organism, such as E. coli, as part of vaccine manufacturing to inoculcate billions of people or producing huge quantities of monoclonal antibodies to treat hundreds of millions of people, large systems are utilized. Such biochemical reactors, with thousands of liters of flowing liquid, are undoubtedly inertia-dominated flows.