Absorption of surfactant-laden droplets into porous media: A numerical study

a r t i c l e i n f o


a b s t r a c t
Hypothesis: Droplets can absorb into permeable substrates due to capillarity. It is hypothesized that the contact line dynamics influence this process and that an unpinned contact line results in slower absorption than a pinned contact line, since the contact area between the droplet and the substrate will decrease over time for the former. Furthermore, it is expected that surfactants can be used to accelerate the absorption. Simulations: Lubrication theory is employed to model the droplet and Darcy's law is combined with the conservation law of mass to describe the absorption dynamics. For the surfactant transport, several convectiondiffusion-adsorption equations are solved. Findings: It is found that moving contact lines result in a parabola-shaped wetted area and a slower absorption and a deeper penetration depth than pinned contact lines. The evolution of the penetration depth was quantitatively validated by comparison with two experimental studies from literature. Surfactants were shown to accelerate the absorption process, but only if their adsorption kinetics are slow compared to the absorption. Otherwise, all surfactant adsorbs onto the pore walls before reaching the wetting front, resulting in the same absorption rate as without surfactants. This behavior agrees with both experimental and analytical literature.

Introduction
Porous media are encountered everywhere. The soil in which seeds are planted, the paper that is used for books and even the human skin are all permeable solids. Despite this widespread occurrence, the flow through porous media is not trivial to understand. Because of the thousands of tiny, twisting, interconnected channels it is complex to track a fluid flow in a porous medium, both numerically and experimentally. Nevertheless, research into this topic is valuable: understanding the flow through porous media allows one to develop or improve a broad range of technologies, such as inkjet printing [1][2][3], irrigation [4], oil recovery [5,6] and even medical treatments [7][8][9].
One of these uninvestigated factors is the dynamics of the contact line of an absorbing sessile droplet on top of the porous medium. In case of partial wetting, dynamic contact line behavior can be roughly categorized into two modes: pinning behavior, where the contact angle changes over time, while the contact radius remains constant, and slipping behavior, where the contact radius changes over time, while the contact angle remains constant [18]. Typical factors that promote slipping behavior over pinning behavior are smoother substrates [19], higher contact angles and increasing deviations from the equilibrium contact angle [20]. Apart from the two extreme possibilities, also a combination is often encountered, called stick-slip behavior, in which the contact line is intermittently pinned and slipping [21].
If the contact line is allowed to slip, the total contact area of the droplet will decrease over time during the absorption process. It is therefore hypothesized that droplets with moving contact lines absorb more slowly than droplets with pinned contact lines. Furthermore, it can be expected that for unpinned contact lines the wetted volume will have a different shape than for pinned contact lines. This is valuable information given that parameters as absorption rate and penetration depth are important in several technologies. For example, in case of inkjet printing slow absorption can cause smudging and shallow penetration may reduce the lifetime of the product.
Another factor in droplet absorption that has not been fully investigated, is the composition of the fluid. Especially the effects of surfactants on droplet absorption have received only limited attention in the experimental literature [22,23] and next to none in numerical articles. This is despite the fact that surfactants are known to be able to 'enhance' (i.e. accelerate) the absorption process, in case a liquid reservoir is absorbed into a porous medium (1D absorption), which was demonstrated both experimentally [5,6,[24][25][26] and theoretically [27]. It is therefore expected that surfactants can cause a similar enhancement to droplet absorption, which is of relevance given the wide usage of surfactants in droplet technology.
In this work the effects of contact line dynamics and surfactants on the absorption of droplets into porous media are considered. Results are obtained with an axisymmetric 3D numerical model that is based on a combination of lubrication theory and Darcy's law. By assuming a relatively small contact angle, an efficient yet accurate model can be developed to simulate droplet dynamics [28,29]. The surfactant transport is modelled by several convection-diffusion-adsorption equations [30,31].
The article has the following structure: first, in Section 2, the equations are introduced that describe the problem and the numerical procedure is outlined. The droplet model and the absorption model are presented and an explanation is given for the surfactant transport equations and the corresponding closure relations. After that, in Section 3, the numerical results are analyzed and validated by comparison with experimental literature. The differences in absorption of a droplet with a pinned and a moving contact line are discussed and the effect of surfactants is considered both qualitatively and quantitatively for various pore sizes. Lastly, in a concluding section, the results are summarized and an evaluation is given of the current state of the art and potential future research.

Mathematical model
In this section, the equations that describe the system are given and explained.

Drop evolution
A sessile droplet on a porous substrate is considered. The droplet initially has a height H at its center, a contact line radius R and a contact angle h. The fluid the droplet consists of is nonvolatile, incompressible and isothermal. Thus, no evaporation occurs and mass density q and dynamic viscosity l remain constant. Only cases are assumed for which H is smaller than the capillary length where c lg denotes the surface tension between the liquid-gas interface, q g the mass density of the gas and g the gravitational acceleration. This means that the effects of gravity can be neglected. The Reynolds number is much smaller than unity and a cylindrical coordinate system (r; /; z) is used to describe the system. The problem is considered axisymmetric, which implies that @ @/ ¼ 0 and it is assumed that there is no swirl (U / ¼ 0). The contact angle h is relatively small (about h < 40 ), meaning that lubrication theory (e ¼ H=R ¼ 1) can be applied, as shown by [28,29]. Given lubrication theory, an evolution equation for the droplet height profile hðr; tÞ can be derived [32]: Here, b is the Navier slip length, which is an indicator for the degree of slip that is allowed, p denotes the pressure and W p the volume flux into the substrate, caused by capillary suction. This volume flux is equal to the vertical fluid velocity w p at z ¼ 0, which will be derived later on. The boundary and initial conditions that define the droplet evolution are given by: @h @r @p @r hðR; tÞ ¼ 0; These equations denote the symmetry conditions at r ¼ 0 (for h; p and c lg ), the contact line position at r ¼ R (where h ¼ 0) and the initial drop profile h 0 , which is a spherical cap. Given that the droplet is considered under partial wetting conditions, there is an associated equilibrium shape (a spherical cap) with a corresponding equilibrium contact angle. This as opposed to complete wetting conditions, where there is no equilibrium associated and the droplet spreads out indefinitely. The equilibrium contact angle can be anywhere between the advancing and receding contact angle of the liquid, depending on the degree of hysteresis. Since in the considered cases the spherical cap shape is formed nearly instantaneously compared to the absorption process (t cap $ lR=c lg % 10 À5 s ( t abs % 1 À 10 À2 s) and the contact angle equilibrates nearly instantaneously (numerical experiments typically show t h % t abs =100), it can be assumed that the droplet is initially in quasi-equilibrium, i.e. the equilibrium shape as if the substrate is nonpermeable (also see [10,11]). The contact line can be either pinned, for which b ¼ 0 and R remains constant, or moving, for which b > 0 and a constitutive relation for R is required that describes the contact line velocity [33]. In this work, we use the relation: Here, k is a typical sensitivity of the contact line position to deviations of the contact angle h from the receding contact angle h rec or the advancing contact angle h adv and a is a power-law index, which can range from 1 to 3. This relation is extensively used and experimentally validated in literature [32][33][34][35][36][37][38][39]. In this work, unpinned contact lines are only considered for pure droplets, because surfactants have a strong tendency to pin the contact line [32,[40][41][42][43].
The pressure in the droplet is dominated by curvature effects of the surface. This means that p can be given by the Laplace pressure: Here, @ r denotes the partial derivative with respect to r. Note that p only depends on r and is independent of z. This is a consequence of the lubrication approximation, meaning that results are more accurate for lower values of e (see e.g. [28,29]).

Absorption model
The droplet is absorbed into the porous substrate as a result of capillary action. This flow can be modeled on a macroscopic level by applying Darcy's law, which is often used for the flow through porous media [2,14,44,45]. This model is chosen for its simplicity over other models, while still being sufficiently accurate [46] and able to deal with the boundary conditions that are involved (e.g. for an impermeable boundary Brinkman's extension would be required [47,48]).
Darcy's law gives a relation for the velocity field u ! p ¼ ðu p ; w p Þ: Here, j is the permeability of the substrate and p p the pressure in the wetted region. A good measure for j can be given by the Carman-Kozeny equation [49,50], which is a model for the flow through a packed bed of solid spheres with diameter d. The permeability is subsequently given by: with g the porosity of the porous medium.
Given mass conservation, it follows that the pressure field can be found by solving the Laplace equation: The boundary conditions which p p is subjected to are given by: Eq. (13) describes the pressure that the droplet exerts on the substrate, Eq. (14) is the no penetration condition at the substrate surface next to the droplet and Eq. (15) is the capillary suction the fluid in the wetted region experiences at its interface, defined at z ¼ Àh p . The corresponding capillary pressure p c can be estimated by considering the capillary action in a single, round channel: Here, h adv is the advancing contact angle. Note that the channel diameter used in the expression for p c is assumed to be equal to the sphere diameter that is used for estimating j. With an expression for the velocity field, the evolution equation for h p can be found, similar to Eq. (1): In the porous medium there are also symmetry conditions at there initially is a thin fluid film h Ã in the porous medium, just below the area covered by the droplet, which is required to remove the incompatibility of Eqs. (1) and (1) for h p ¼ 0.

Surfactant in the droplet
At the liquid-air interface of the droplet a surfactant concentration Cðr; tÞ is defined. The transport equation for Cðr; tÞ is given by [51]: Here, U t is the fluid velocity tangential to the liquid-air interface, D C is the surface diffusion coefficient and J CC is the sorptive transport between the interface and the bulk. The fraction p @ @r denotes the surface derivative. The right-hand side of Eq. (18) consists of several terms of which each corresponds to a specific transport contribution. The first term denotes the convective fluid transport tangential to the surface, the second term is the change rate in concentration as a result of changes in curvature, the third term is the diffusion rate and the fourth term corrects for the transformation of the surface coordinates to the cylindrical coordinates [31,32].
Similarly, a transport equation for the surfactant bulk concentration Cðr; tÞ can be given. This concentration is assumed to be independent of z, meaning that any vertical concentration gradient is considered insignificant. This is a generally used assumption in literature when the lubrication approximation is used [31,52,53]. The transport equation is defined in terms of wðr; tÞ ¼ Cðr; tÞhðr; tÞ, because in that case the dependent variable becomes independent of hðr; tÞ [30]: Here, D C denotes the bulk diffusion coefficient and C p ðr; z; tÞ is the bulk surfactant concentration in the pores. The first term on the right-hand side consists of three parts: a pressure-driven convection part, a part that accounts for the Marangoni effect and a diffusion part. The second term is the sorptive exchange with the liquid-air interface, including a factor that takes into account the slope of the interface ( @s @r ). The third term is the convective transport flux with the porous medium and the last term is the diffusive exchange with the porous medium.
The concentrations Cðr; tÞ and Cðr; tÞ are subject to the following boundary conditions and initial conditions: wðr; 0Þ ¼ C 0 hðr; 0Þ: ð23Þ The boundary conditions denote the symmetry condition and the no-flux condition at the contact line.
The surfactants at the interface Cðr; tÞ tend to reduce the liquidgas surface tension c lg , meaning that an equation of state is required to close the problem. In this work, the Szyszkowski equation of state is chosen for this purpose, which takes into account the repelling effect individual surfactant molecules have on each other [54]. This closure relation is typically valid for intermediate interfacial concentrations that are not too close to the maximum concentration C 1 . For lower concentrations the equation reduces to the linear, dilute equation of state. The Szyszkowski equation is given by: In this equation, c lg;0 denotes the liquid-gas surface tension for a surfactant-free liquid, R u is the universal gas constant and T the temperature.
The sorptive exchange depends on both the bulk concentration Cðr; tÞ and the interfacial concentration Cðr; tÞ and on the available space at the interface [31,53,55]. Thus, the following equation for J CC can be assumed: Here, k C a and k C d are the adsorption and desorption coefficients of the liquid-air interface respectively. In this equation it can be recognized that the positive adsorption term increases for higher bulk concentrations and lower interfacial concentrations, while the negative desorption term only depends on C [55].

Surfactant in the porous medium
In contrast to the surfactant bulk concentration Cðr; tÞ in the droplet, the bulk concentration in the porous medium C p ðr; z; tÞ is also assumed to depend on the axial coordinate z. The reason for this is that the vertical dimension of the wetted region typically is of the same order of magnitude as the horizontal dimension. Furthermore, the pressure gradient has a significant component in axial direction. The evolution of C p ðr; z; tÞ is thus to be described by an axisymmetric 3D convection-diffusion-adsorption equation: Here, J CS accounts for the sorption between the bulk and the walls of the pores. This effect is not considered in the droplet itself, because there the total liquid-solid interface is much smaller. Furthermore, adsorption onto the solid-liquid interface does not directly influence the flow behavior (it cannot cause Marangoni flow) and will not change the contact angle of the droplet, because only pinned cases are considered when surfactants are involved. The concentration C p ðr; z; tÞ is subject to the following boundary and initial conditions: C p ðr; 0; tÞw p ðr; 0; tÞ ¼ Cðr; tÞW p ðr; tÞ; ð28Þ Furthermore, surfactant cannot be transported outside the wetted region.
The sorption term J CS is given by: with k S a and k S d the solid-liquid adsorption and desorption coefficient respectively, Sðr; z; tÞ the amount of adsorbed surfactant per unit of volume and S 1 the maximum adsorbed surfactant concentration per unit of volume. The 4 d prefactor originates from the expression for the area of the channel walls A p for a given control volume V ¼ A p d=ð4gÞ (see [2] for a derivation). At t ¼ 0 no surfactant has adsorbed on the pore walls (Sðr; z; 0Þ ¼ 0).
As an equation of state for the surface tension of the solid-liquid interface c sl a variant of the Sheludko equation [56,57] is used: Adsorbing surfactants will therefore increase the magnitude of the capillary suction pressure. The volume-averaged value of Sðr; z; tÞ at the interface, i.e. S À int , is used for calculation of p c . This results in better stability, since p c becomes more uniform.
The surfactant concentration at the liquid-air interface in the porous medium is not taken into account, because it does not affect p c (as shown by [58,59]). Furthermore, due to the adsorption kinetics Ch ) C, meaning that adsorption onto the liquid-air interface has no significant influence on the bulk concentration. This type of surfactant is typically called a 'penetrant' and is used to influence liquid absorption.
The numerical procedure that is employed to solve the equations outlined in this section, is given in Section S.1 of the Supplementary Material.

Results
In this section the results of the simulations are presented and discussed. First, pure droplets with both pinned and moving contact lines are considered and quantitatively compared with experimental literature. After that, droplets with surfactants are examined and results are shown to be consistent with literature as well.

Pure droplets
Droplets with both a pinned contact line and a moving contact line can be absorbed by a porous substrate. This results in different shapes of the absorbed wetted region as can be seen in Fig. 1.
In the pinned case (Fig. 1a) the depth profile in the porous medium propagates with a flat front and the wetted region also tends to expand in radial direction. In the unpinned cases (Fig. 1b, c, d), however, the wetting front becomes increasingly parabolic in shape and the 'contact line' in the porous medium behaves as if it were pinned. This behavior seems to be independent of initial contact angle, although the final shape becomes flatter as the initial contact angle is decreased, because of the change in aspect ratio of the droplet.
Note that decreasing the initial contact angle can be done in two ways: changing the initial radius, while keeping the initial volume constant (Fig. 1c) and changing the initial volume, while keeping the initial radius constant (Fig. 1d). Regardless of the way in which the initial contact angle is reduced, the qualitative behavior remains the same: a relatively flatter wetted region.
An explanation for these differences and similarities can be found by analyzing the pressure in the porous medium, as for example given by Fig. 2.
The pressure right underneath the droplet is equal to the Laplace pressure in the droplet (~Oð1Þ Pa) and from there decreases towards the edge of the wetted region, where the pressure equals the capillary pressure. The vertical gradient towards the horizontal wetting front is more or less identical for all cases, while the horizontal gradient towards the side of the wetted region is different, because of the differences in contact radii. Since for the moving contact line cases, the droplet contact line moves away from the wetted region contact line, there is barely any radial expansion. This is not the situation for the pinned contact line case, because there the droplet contact line remains close to the wetted region contact line, giving it the opportunity to expand in radial direction as well. The same mechanism also explains the differences in shape of the wetted regions: as the contact area shrinks for the moving contact line cases, an increasing proportion of the wetted region stops expanding, because it is too far away from the droplet in radial direction.
This behavior can also be observed in the depth and volume evolution of the wetted region. As can be seen in Fig. 3 the penetration depth H p , defined as the deepest point of the wetted region, evolves nearly identical for both the pinned and unpinned cases.
However, the final penetration depth tends to be larger for the unpinned cases, because of the shape of the wetted region. Since for the unpinned cases the contact area shrinks over time, the result is a more pointed shape than for the pinned case (also see Fig. 1) and a relatively slower volume evolution, as seen in Fig. 3b. It also shows that decreasing the initial contact angle generally increases the absorption rate: for the constant initial volume because the total contact area becomes larger and for the constant initial radius because the total volume to be absorbed becomes smaller.
Note that the pinned contact line case seems to have a smaller final absorbed volume than the other moving contact line cases with the same initial drop volume. This is, however, a result of the simulations being cut off for different remaining drop volumes. The pinned contact line case tends to be less stable than the moving contact line cases, because its aspect ratio changes more severely. Therefore the simulations with pinned contact lines are stopped when the remaining drop volume is larger than the drop volume that remains when the simulations with moving contact lines are stopped. The input parameters that were used for Figs. 1-3 are given in Section S.2 of the Supplementary Material.
The implication of these findings is that the evolution of the penetration depth is mostly independent of the shape of the liquid reservoir on the surface. However, the final penetration depth and the volume evolution do depend on this shape. These trends have not been noted before to the knowledge of the authors. A practical consequence of this is that the simulations can be compared to experiments that involve droplets, but also to experiments that use a container filled with fluid as a reservoir (1D absorption). This is illustrated in Fig. 4, where simulation results are compared to droplet experiments by Nees (2011) [2,60] and container experiments by Starov et al. (2004) [27].
The experiments by Nees (2011) involve a 10 lL droplet that absorbs into a porous substrate that is made by melting glass beads together. There is some spread in the experimental data, which is explained by the inhomogeneities in the porous medium. The droplet experiments agree with the predictions made by the numerical model, as can be seen in Fig. 4a.
The experiments by Starov et al. (2004) are performed by submerging a bound nitrocellulose membrane in a liquid-filled bath, resulting in a uniform, onedimensional evolution of the liquid front. They do this for both pure water and SDS solution. As can be seen in Fig. 4b, the droplet simulations agree well with the experiments. Note that the experimental results that are shown in Fig. 4b were actually performed with a small SDS concentration. Starov et al. (2004) show however that in this case the results with surfactant are the same as without surfactants. Since these were the only time-dependent results that were given in the article, they are used for this comparison with pure droplets and can safely be assumed appropriate. In This approximation follows both the simulations and experiments rather well, as can be seen in Fig. 4. As a consequence, it is possible to derive parameters (j, p c , g) from results by fitting a square root function to it. Subsequently, a single, unknown parameter (in this work jp c ) can be extracted from the fitted prefactor. This method will be applied in the next subsection.

Droplets with surfactants
When surfactants are added to the droplet, this can influence the absorption process. As surfactant molecules adsorb onto the pore walls, the capillary pressure has a tendency to increase (as implied by Eqs. (31) and (32)), which potentially accelerates the absorption. As for example illustrated in Fig. 5, adding surfactants results in a faster evolution of the penetration depth H p and absorbed volume V p . This is indeed confirmed by literature [27,61,62,63,64].
The magnitude of the acceleration depends on the adsorbed surfactant concentration at the liquid front and is mainly limited by surface energy effects, i.e. c sg and the range of c sl (which ranges from c sl;0 to c sl;1 ). High values of c sg and low values of c sl;0 À c sl;1 will typically diminish the potential for accelerated absorption by surfactants, because they decrease the relative effect of adsorbed surfactant on p c .
The adsorption kinetics of the surfactant also determine whether the absorption process accelerates. As illustrated in Fig. 6, a relatively fast adsorption causes all surfactant to be consumed before it reaches the wetting front.
In Fig. 6a there is only a significant surfactant concentration present in the top part of the wetted region, i.e. close to the droplet, which functions as a surfactant source. However, since the adsorption kinetics are relatively fast compared to the liquid absorption, only surfactant is adsorbed in the upper part of the wetted region (see Fig. 6b). As a result, there is no surfactant left to adsorb in the bottom part. The result is that the surfactant has no influence on the wetting properties of the liquid, because the capillary pressure is the same as for a pure liquid.
If the pore diameter is larger, however, the droplet absorbs faster, meaning that the adsorption and absorption processes have similar time scales. Furthermore, for larger pore diameters the specific surface area is smaller, hence a smaller amount of surfactant can be adsorbed. As a result, the surfactant can get closer to the wetting front before being adsorbed. This can be seen in Fig. 6c, where the surfactant concentration is nearly uniformly distributed over the wetted region, except close to the interface, where it is smaller. Naturally, a similar concentration distribution can be found for the adsorbed concentration S (see Fig. 6d). This causes the absorption process to accelerate, because the magnitude of the capillary pressure -and thus the pressure gradient in the wetted region -increases.
An effective way of actually quantifying the absorption process is through the value of jp c (the product of the permeability and the capillary pressure). As shown by Eq. (33), this product mainly determines the evolution velocity of H p and is not known As can be seen from Fig. 7, the value of jp c indeed increases with C 0 for larger pore diameters, while it remains constant for smaller pore diameters. As shown before in Fig. 6, smaller adsorption rates (small pore diameters) will cause all surfactant to adsorb before reaching the wetting front, while larger adsorption rates (large pore diameters) result in an increase of the magnitude of p c by surfactant. For the latter case the amount of surfactant adsorbing at the wetting front naturally increases as the surfactant concentration increases, which can be seen for d ¼ 3:0 lm in Fig. 7. This is useful data, because it shows that increasing the amount of surfactant in a solution does not necessarily accelerate the absorption. If the surfactant has adsorption kinetics that are simply too fast, it will never affect the absorption dynamics. Therefore, this could, for example, guide engineers to use a slower surfactant if an enhanced droplet absorption is desired.
Similar behavior was also found for 1D absorption by Starov et al. (2004) [27]. They performed experiments on the absorption   As a final note, jp c at C 0 ¼ 0 increases linearly with d, which is what would be expected from Eqs. (11) and (16).

Conclusion
In this work a numerical model was created for the absorption of surfactant-laden sessile droplets into porous media. The droplet was modeled using axisymmetric lubrication theory and the absorption process using mass conservation and Darcy's law. In the droplet, the surfactants were considered both at the liquid-air interface and in the bulk using convection-diffusion-adsorption equations and in the porous medium the surfactants were modeled in the bulk through another convection-diffusion-adsorption equation and at the pore walls through an adsorption equation.
The contact line of the droplet was considered for both a pinned case and an unpinned case with a slip model. It was shown that if the contact line is pinned, the droplet absorbs with a mostly flat wetting front, which also expands in radial direction. On the other hand, for a slipping contact line, the wetting front becomes increasingly parabolic in shape and the 'contact line' of the wetted region remains nearly pinned. This results in similar evolution of the penetration depth for both contact line models, while the volume of the wetted region evolves more slowly for the slip model. Therefore, the final penetration depth tends to be larger for an unpinned contact line. The evolution of the penetration depth was validated by quantitative comparison with experimental results from literature. It was also found that a reduced initial contact angle results in faster absorption and a flatter -although still parabolic -shape of the wetted region, which is a result of the relatively larger contact area and flatter droplet.
When surfactants are involved, it was shown that there is a potential for the absorption process to accelerate. Whether this acceleration occurs depends on the time scale of the adsorption kinetics of the surfactants compared to the time scale of absorption of the droplet. Small absorption rates (small pore sizes) result in surfactant adsorbing on the pore walls before reaching the wetting front. Thus, in that case the absorption process of the surfactantladen droplet will be indistinguishable from the absorption of a pure droplet. If the absorption process is fast (for large pore sizes), however, surfactant does adsorb at the wetting front, resulting in more suction. As a consequence, the absorption rate will increase with initial surfactant concentration if the pore size is large enough. These results are in agreement with 1D experimental and analytical results by Starov et al. (2004) [27].
These findings agree with the hypothesis that was made before, namely that a moving contact line slows down the absorption process and that surfactants can be used to accelerate the absorption. This acceleration only happens if the absorption -before taking into account surfactants -is fast enough.
Possible applications lie in the ability to control the absorption rate, final penetration depth and final shape of the wetted region. For example, if the wetted region requires a final shape that is mostly flat, it helps to add components that promote pinning, such as surfactants that reduce the contact angle [32,43]. Furthermore, if the absorption needs to be accelerated, it is important to choose a surfactant that does not adsorb too fast (typically larger surfactants) and that have a tendency to adsorb onto the liquid-solid interface. These type of methods for control can be employed in several technologies that involve porous media, such as inkjet printing, irrigation and medical treatmens that involve absorption through the skin. Previous numerical studies often employed 1D absorption models [27,[70][71][72] rather than 2D or 3D and although there are several articles on 2D or 3D droplet absorption models [10][11][12][13][14][15][16][17], none considered the effect of the contact line dynamics and/or involved surfactants. Furthermore, while a significant number of experimental studies on surfactant-enhanced liquid absorption has been carried out [5,6,[24][25][26], it tends to be rather difficult to actually visualize the flow, let alone to quantify properties as pressure and surfactant concentration. Therefore, numerical studies, like this one, are a valuable contribution to the field.
Further research opportunities lie in additional parametric analysis of the relevant surfactant properties and comparison with experiments. While extensive analysis has been done on the effect of several dimensionless numbers on the absorption of pure droplets (e.g. [2,14]), this has not been done for surfactants. This article has made a first step in that direction. Furthermore, experimental studies that systematically consider the effect of different surfactants on absorption are rare. This would be a requirement for validating potential numerical results.
Also, Molecular Dynamics (MD) simulations can help to increase the accuracy of our model. In the initial stages of the absorption process and for nanoscale pores the current hydrodynamic model tends to be less valid. For both these limits, MD simulations can offer corrections to the equations (e.g. see [65][66][67]). Furthermore, MD models can be employed to get better insight into the adsorption kinetics of surfactants and the factors that influence this process. This way, it becomes possible to compare our model quantitatively with experiments, because real-life surfactants can be implemented if their adsorption parameters are known [68,69].
An additional suggestion for improving the model is to define a local capillary pressure rather than one that is averaged over the wetting front. This would allow one to capture more details of the flow. The corresponding numerical instabilities may be solved by employing a different, 2D method to define the interface position, since in the current model this is done by a 1D height profile h p ðrÞ. Also, in many cases porous substrates swell as a result of the absorbed liquid (e.g. paper, see [70][71][72]), which can be modeled as well. Furthermore, the consideration of volatile fluids, i.e. simultaneous evaporation and imbibition, is of relevance and can be taken into account in future studies.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.