2.1. Poro-elastocapillary substrate

The deformation of the polymer matrix is characterised by a mapping from a reference configuration prior to deformation (i.e. a dry gel without any solvent) to the current state. By standard notation, this mapping is expressed as

(2.1) \begin{align} \boldsymbol{x}=\boldsymbol{\psi }(\boldsymbol{X}, t). \end{align}

Here, $\boldsymbol{X}$ represents the position of a material point in the reference domain, which is mapped to its current position $\boldsymbol{x}$ by the deformation map $\boldsymbol{\psi }$ . We assume that the geometry is invariant along the contact line, allowing us to treat the problem as two-dimensional (plane strain elasticity). The polymer matrix can be swollen by a solvent, so the mapping $\boldsymbol{\psi }$ does not preserve volume. Defining the deformation gradient tensor as $\boldsymbol F =\partial \boldsymbol{x} / \partial \boldsymbol{X}$ , the relative change in volume is given by

(2.2) \begin{align} J=\mathrm{det}(\boldsymbol F). \end{align}

Using the dry network without any solvent as the reference configuration, $J$ naturally defines the ‘swelling ratio’. Next, it is assumed that both the polymer molecules that make up the matrix and the solvent molecules are incompressible. This assumption implies that the volumetric change of the network directly determines the local number density of the solvent molecules. Defining $C(\boldsymbol X,t)$ as the number of solvent molecules per unit volume in the reference configuration, the molecular incompressibility is ensured by the constraint (Hong et al. Reference Hong, Zhao, Zhou and Suo2008)

(2.3) \begin{align} J = 1+v C, \end{align}

where $v$ is the volume per solvent molecule. The system is thus described by two fields $\boldsymbol{\psi }(\boldsymbol X,t)$ and $C(\boldsymbol X,t)$ , which are constrained by (2.3).

The free energy of the polymer matrix is described using a hyperelastic formulation based on a strain-energy density, $W_{\textit{el}}(\boldsymbol{F})$ , per unit volume in the reference configuration. In the present study, we adopt a common neo-Hookean model based on Gaussian chain elasticity, which for plane-strain conditions gives the energy density

(2.4) \begin{align} W_{\textit{el}}(\boldsymbol{F})=\frac {1}{2} G[\boldsymbol{F}: \boldsymbol{F}-2-2 \log (\operatorname {det}(\boldsymbol{F}))], \end{align}

where $G$ is the shear modulus of the network. The network elasticity is of entropic origin, which can be appreciated from the connection $G=NkT$ , where $N$ represents the number of chains per unit volume (Boyce & Arruda Reference Boyce and Arruda2000; Marckmann & Verron Reference Marckmann and Verron2006; Hong et al. Reference Hong, Zhao, Zhou and Suo2008; Kim, Yin & Suo Reference Kim, Yin and Suo2022), and $kT$ is the thermal energy. The free energy for the mixing of the solvent and polymer can be described using the Flory–Huggins theory. The corresponding energy density with respect to the volume in the reference configuration takes the form (Hong et al. Reference Hong, Zhao, Zhou and Suo2008)

(2.5) \begin{align} W_{\textit{FH}}(C)=-\frac {k T}{v}\left [v C \log \left (1+\frac {1}{v C}\right )+\frac {\chi }{1+v C}\right ]\!. \end{align}

The first term represents the entropy of mixing, while the second term is the enthalpy of mixing, the strength of which is set by the interaction parameter $\chi$ . The free-energy density for the swollen poroelastic network comprises the sum $W_{\textit{el}}(\boldsymbol{F})+W_{\textit{FH}}(C)$ .

Next, we address the capillary contributions. These contributions come in two forms. First, we endow the surface of the substrate with a surface energy $\gamma _s$ . Adopting the convention that is common in capillarity, this surface free energy is measured per unit area in the current configuration. In principle, this energy can be a function of the degree of swelling at the interface, which would give rise to the Shuttleworth effect (Pandey et al. Reference Pandey, Andreotti, Karpitschka, Van Zwieten, van Brummelen and Snoeijer2020; Henkel et al. Reference Henkel, Essink, Hoang, van Zwieten, van Brummelen, Thiele and Snoeijer2022); for simplicity it is assumed that $\gamma _s$ is independent of deformation. The second capillary effect is due to the liquid–vapour surface tension $\gamma _{\textit{LV}}$ . This can be represented by a localised force (Andreotti & Snoeijer Reference Andreotti and Snoeijer2020), pulling along the liquid–vapour interface, and leads to the formation of the elasto-capillary ridge, as indicated in figure 1(b). Accordingly, we introduce a perfectly localised external traction, characterised by the associated work functional $\mathcal{W}(\boldsymbol{\psi })=\gamma _{\textit{LV}} \boldsymbol{t}_{\textit{LV}} \boldsymbol{\cdot }\boldsymbol{\psi }(\boldsymbol{X}_{{cl}})$ , where $ \boldsymbol{X}_{{cl}}$ denotes the position of the contact line on the solid surface in the reference configuration, and $\boldsymbol t_{\textit{LV}}$ is the pulling direction.

In summary, the functionals describing the free energy of the poro-elastocapillary substrate, and the work performed by the liquid–vapour surface tension, per unit length along the contact line, are given by

(2.6a) \begin{align} \mathcal{E}[\boldsymbol{\psi },C,\varPi ] &= \int {}{\mathrm{d}}^2X\,\big (W_{\textit{el}}+W_{\textit{FH}} + \varPi (1+vC-J)\big ) +\int {}\mathrm{d}s\,\gamma _s, \end{align}

(2.6b) \begin{align} \mathcal{W}[\boldsymbol{\psi }]&=\gamma _{\textit{LV}} \boldsymbol{t}_{\textit{LV}} \boldsymbol{\cdot }\boldsymbol{\psi }(\boldsymbol{X}_{{cl}}), \end{align}

respectively, where we incorporated the incompressibility constraint (2.3) into the free-energy function via the Lagrange multiplier $\varPi$ . This Lagrange multiplier turns out to be the osmotic pressure. The integration measure $\mathrm{d}s$ represents the arc length in the deformed configuration. We remark that realistic experimental systems can have more intricate energies than assumed here (i.e. beyond neo-Hookean, Flory–Huggins and constant surface tension). However, this does not change the form (2.6a ), and can in principle be included in the proposed modelling framework. In the Conclusion we return to this point.

In the model that we consider, the poroelastic substrate is contiguous to a liquid reservoir. The absorption of solvent molecules into and deformation of the poroelastic substrate also affect the energy of this liquid reservoir. We assume that during its interaction with the substrate, the liquid reservoir remains at constant pressure, $p_0$ , and chemical potential, $\mu _0$ . Accordingly, the deformation of the substrate produces additional energy in the amount of $p_0V$ via work on the reservoir, and the absorption of solvent molecules consumes energy in the amount $\mu _0N_c$ , where $V$ and $N_c$ denote the volume of the substrate and the number of solvent molecules it contains, respectively. Noting that $C=(J-1)/v$ by (2.3), the energetic contributions from the interaction of the substrate with the liquid reservoir can be expressed in terms of the deformation and concentration as

(2.7) \begin{align} \begin{aligned} p_0V-\mu _0N_c &= p_0\int {}{\mathrm{d}} ^2X\,J-\mu _0\int {}{\mathrm{d}} ^2X\,C \\ &= \tilde {p}_0\int {}{\mathrm{d}} ^2X\,J+\operatorname {const}, \end{aligned} \end{align}

with $\tilde {p}_0=p_0-\mu _0/v$ . The identities in (2.7) convey that, for an incompressible system, an exchange of molecules cannot be distinguished from an exchange of volume. Hence, without loss of generality, we will in the remainder set the chemical potential of the reservoir to $\mu _0=0$ . For an incompressible system, the value of $\tilde{p}_0$ is without any consequence as well, as it only sets the gauge for the stress. By setting $p_0=0$ , it is thus understood that we measure the pressure with respect to that of the reservoir. Hence, in the remainder it will suffice to ignore the energetic contributions of the reservoir, under the convention that the gauge of the reservoir is chosen with $\mu _0=0$ and $p_0=0$ .

2.2. Dynamical equations and dissipation

We now consider the evolution equations that model the dynamics of wetting ridges in the presence of poroelasticity and viscoelasticity. In modelling the dynamics of poro-viscoelastic materials, it is generally appropriate to neglect inertial effects. Therefore, we consider the momentum balance in quasi-stationary form, and the solvent evolution in time-dependent form.

We present the conservation equation for the solvent in the current configuration. To this end, we denote by $\boldsymbol{v}(\boldsymbol{x},t)=\partial _t\boldsymbol{\psi }(\boldsymbol{X},t)$ the material velocity of the substrate, and by $c(\boldsymbol{x},t)=C(\boldsymbol{X},t)/J(\boldsymbol{X},t)$ the solvent concentration in the current configuration, under the deformation in (2.1). Conservation of the solvent is then described by a conservation equation of the form

(2.8) \begin{align} \partial _tc+\boldsymbol{\nabla }{}\boldsymbol{\cdot }(c\boldsymbol{v})+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{q}_c=0 ,\end{align}

where $\boldsymbol{\nabla}$ represents the gradient operator in the current configuration and $\boldsymbol{q}_c$ denotes the diffusive flux. To provide a constitutive relation for the diffusive flux, we introduce the chemical potential of the solvent in the reference configuration according to

(2.9) \begin{align} M:=\frac {\delta \mathcal{E}}{\delta {}C} = W'_{\textit{FH}}+\varPi {}v ,\end{align}

where $(\boldsymbol{\cdot })'$ denotes differentiation. The chemical potential in the current configuration is then given by $\mu (\boldsymbol{x},t)=M(\boldsymbol{X},t)$ . We define the diffusive flux using Fick’s law (corresponding to type I diffusion for very soft polymer networks; Hong et al. Reference Hong, Zhao, Zhou and Suo2008)

(2.10) \begin{align} \boldsymbol{q}_c=-D\boldsymbol{\nabla }\mu ,\end{align}

with a constant mobility parameter $D\gt 0$ . The mobility governs the diffusive transport of solvent, and will in the remainder be referred to as a diffusion coefficient. More sophisticated models of the diffusive process can be realised by replacing $D$ with a symmetric-positive-definite tensor $\boldsymbol{D}$ , possibly dependent on $\boldsymbol{\psi }$ and $C$ (for example see Bertrand et al. Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016).

Neglecting inertial effects in the motion of the substrate, conservation of linear momentum is expressed by the static-equilibrium relation

(2.11) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }=0, \end{align}

with $\boldsymbol{\sigma }$ the Cauchy stress tensor. We consider constitutive relations for $\boldsymbol{\sigma }=\boldsymbol{\sigma }_{\textit{el}} + \boldsymbol{\sigma }_{{v}}$ , that comprise an elastic part $\boldsymbol{\sigma }_{\textit{el}}$ and (potentially) a viscous part $\boldsymbol{\sigma }_{{v}}$ . For the considered neo-Hookean hyperelastic material model, the elastic Cauchy stress tensor follows from the free energy as

(2.12) \begin{align} \boldsymbol{\sigma }_{\textit{el}}(\boldsymbol{x},t)=\big [J^{-1}\boldsymbol{P}\boldsymbol{\cdot }\boldsymbol{F}^{T}-\varPi \boldsymbol{I}\big ](\boldsymbol{X},t) ,\end{align}

with $\boldsymbol{P}:=\mathrm{d}W_{\textit{el}}(\boldsymbol{F})/\mathrm{d}\boldsymbol{F}$ as the first Piola–Kirchhoff stress tensor corresponding to the elastic part. Note that we have integrated the pressure part $\varPi \boldsymbol{I}$ originating from the constraint (2.3) into the elastic Cauchy stress. The same Lagrange multiplier $\varPi$ appears in the chemical potential (2.9), and lies at the origin of the poroelastic coupling. For the viscous part, we adopt a Newtonian stress–strain-rate relation

(2.13) \begin{align} \boldsymbol{\sigma }_{{v}}(\boldsymbol{x},t) = \eta \big (\boldsymbol{\nabla }\boldsymbol{v} + (\boldsymbol{\nabla }\boldsymbol{v})^T\big ) ,\end{align}

with $\eta \gt 0$ as the viscosity. The combined elastic and viscous stress tensors provide a Kelvin–Voigt visco-elastic model.

At the interface between the poroelastic substrate and the liquid reservoir, the traction exerted by the elastic and viscous stresses are balanced by the surface tension. This dynamic interface condition is encoded as

(2.14) \begin{align} \boldsymbol{\sigma }\boldsymbol{\cdot }\boldsymbol{n} =\gamma _s\kappa \boldsymbol{n} ,\end{align}

where $\boldsymbol{n}$ denotes the exterior unit normal vector and $\kappa$ denotes the curvature of the interface, under the convention that curvature is negative if the osculating circle is located in the substrate. Condition (2.14) holds everywhere on the interface, except at the contact line, where the interface generally exhibits a kink and, accordingly, the curvature is undefined. A separate, pointwise condition must be specified at the contact line. We impose that the external traction is balanced by the surface-tension tractions acting on both sides of the contact line, according to

(2.15) \begin{align} \gamma _{\textit{LV}}\boldsymbol{t}_{\textit{LV}}-\gamma _s(\boldsymbol{t}_-+\boldsymbol{t}_+)=0 ,\end{align}

where $\boldsymbol{t}_{\pm }$ are the unit conormal vectors, normal to the contact line, and tangential and external to the smooth parts of the interface adjacent to the contact line. It is noteworthy that (2.15) in fact corresponds to the Neumann law. The conditions (2.14) and (2.15) are naturally satisfied in the equilibrium state (Pandey et al. Reference Pandey, Andreotti, Karpitschka, Van Zwieten, van Brummelen and Snoeijer2020), and imposing them dynamically ensures there is no dissipation associated with capillarity. Likewise, the adjacency of the liquid reservoir implies that the chemical potential at the interface vanishes at equilibrium, and we choose to impose $\mu =0$ dynamically to avoid any source of dissipation at the interface. This boundary condition implies that the substrate is assumed to be in contact with a solvent reservoir, allowing solvent exchange with the external phase. In terms of figure 1(a) this means that in our simulations both the droplet and the surrounding phase can swell the substrate.

The part of the substrate boundary complementary to the interface comprises the bottom and lateral boundaries. We assume that the poroelastic substrate is subjected to a preswelling treatment, in such a manner that the preswelling leads to a homogeneous isotropic expansion of the dry substrate with swelling ratio $\lambda$ . Subsequently, the bottom boundary is regarded as fixed in the preswollen configuration and, accordingly, we impose $\boldsymbol{\psi }(\boldsymbol{X},t)=\lambda \boldsymbol{X}$ . On the lateral boundaries, we impose free slip in the preswollen configuration, i.e. $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\psi }(\boldsymbol{X},t)=\boldsymbol{n}\boldsymbol{\cdot }\lambda \boldsymbol{X}$ and $\boldsymbol{n}\times (\boldsymbol{\sigma }\boldsymbol{\cdot }\boldsymbol{n})\times \boldsymbol{n}=0$ with $\times$ the cross-product. We regard the bottom and lateral boundaries as closed and therefore set $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\mu =0$ .

Our interest pertains to the evolution of the wetting ridge due to a point load at the interface associated with a liquid–vapour interface. We hence regard a preswollen initial configuration corresponding to the equilibrium of the poroelastic substrate in contact with the liquid reservoir, in the absence of the point load. The initial deformation then coincides with that of the preswollen configuration, i.e. $\boldsymbol{\psi }(\boldsymbol{X},t=0)=\lambda \boldsymbol{X}$ , and the initial concentration is uniform and such that the chemical potential (relative to the liquid reservoir) vanishes.

Thermodynamic consistency, in the sense of dissipation of the formulation without capillarity and without viscoelasticity, was discussed in Hong et al. (Reference Hong, Zhao, Zhou and Suo2008). The effect of $\boldsymbol \sigma _{ {v}}$ can be included in the energy-balance equation as well. Following the steps in Pandey et al. (Reference Pandey, Andreotti, Karpitschka, Van Zwieten, van Brummelen and Snoeijer2020) and making use of the dynamical equations (2.8) and (2.11) the free-energy functional evolves according to

(2.16) \begin{align} \mathrm{d}_t(\mathcal{E}-\mathcal{W}) = - \int {\mathrm{d}}^2x\, \big (2\eta |\boldsymbol{\nabla }\boldsymbol{v}|^2+D|\boldsymbol{\nabla }\mu |^2\big ) + \operatorname {bnd} ,\end{align}

where $|\boldsymbol{\nabla }\boldsymbol{v}|^2=\boldsymbol{\nabla }\boldsymbol{v}:\boldsymbol{\nabla }\boldsymbol{v}$ and $\operatorname {bnd}$ denotes terms on the part of the boundary of the substrate complementary to the interface; see Appendix B. Equation (2.16) conveys that the free energy is non-increasing, up to work supplied by liquid–vapour surface tension and, possibly, terms on the complementary part of the substrate boundary. Under the aforementioned boundary conditions, however, $\operatorname {bnd}$ vanishes so that dissipation only occurs in the bulk of the substrate. Moreover, because equilibrium configurations are characterised by vanishing dissipation, (2.16) conveys that $\boldsymbol{v}$ and $\mu$ are uniform in equilibrium. Subject to the aforementioned boundary conditions on the auxiliary part, $\boldsymbol{v}=0$ and $\mu =0$ in equilibrium.

2.3. Weak formulation

To further elucidate the poro-visco-elastocapillary model, and to provide a basis for the finite-element approximation in § 2.5, we regard the weak formulation of (2.8)–(2.15) constrained by (2.3).

A weak formulation of (2.11) subject to the incompressibility constraint (2.3), including the viscous part (2.13), with interface conditions (2.14) and Neumann law (2.15), can be conveniently expressed in terms of derivatives of the energy and work functionals as

(2.17) \begin{align} (\boldsymbol{\psi }(\boldsymbol{\cdot },t),\varPi (\boldsymbol{\cdot },t)):\\ \nonumber \partial _{\boldsymbol{\psi }}\mathcal{E}[\boldsymbol{\psi },C,\varPi ](\boldsymbol{W}) +\int {\mathrm{d}}^2x\,\boldsymbol{\nabla }\boldsymbol{w}:\boldsymbol{\sigma }_{{v}} \\ \nonumber +\partial _{\varPi }\mathcal{E}[\boldsymbol{\psi },C,\varPi ](V)=\partial _{\boldsymbol{\psi }}\mathcal{W}[\boldsymbol{\psi }](\boldsymbol{W}) \\ \nonumber \quad \quad \forall (\boldsymbol{W},V), \end{align}

where $\partial _{\boldsymbol{\psi }}\mathcal{E}[\boldsymbol{\psi },C,\varPi ](\boldsymbol{W})$ represents the Gateaux derivative of the free-energy functional with respect to its deformation argument in the direction $\boldsymbol{W}$ , in particular,

(2.18) \begin{align} \partial _{\boldsymbol{\psi }}\mathcal{E}[\boldsymbol{\psi },C,\varPi ](\boldsymbol{W}) =\frac {{\mathrm{d}}}{\mathrm{d}s}\mathcal{E}[\boldsymbol{\psi }+s\boldsymbol{W},C,\varPi ]\Big |_{s=0} .\end{align}

This notational convention for the Gateaux derivative carries over to other instantiations of $\partial _{(\boldsymbol{\cdot })}$ . The functions $(\boldsymbol{W},V)$ act as test functions in the reference configuration, while $\boldsymbol{w}(\boldsymbol{x})=\boldsymbol{W}(\boldsymbol{X})$ is the representation of $\boldsymbol{W}$ in the current configuration. The trial function $(\boldsymbol{\psi },\varPi )$ and test function $(\boldsymbol{W},V)$ are supposed to belong to suitable admissibility classes, and $\boldsymbol{\psi }$ and $\boldsymbol{W}$ are subject to conditions at the complementary boundary in accordance with imposed essential boundary conditions. By means of functional derivatives analogous to those in Appendix B, one can show that (2.17) is indeed equivalent to (2.11)–(2.15).

The weak formulation of (2.8) is most conveniently expressed in a form that represents the transport part in the reference configuration and the diffusive part in the current configuration

(2.19) \begin{align} C(\boldsymbol{\cdot },t):\quad \frac {{\mathrm{d}}}{\mathrm{d}t} \int {\mathrm{d}}^2X\,{\textit{CY}} + \int {\mathrm{d}}^2x\,D\boldsymbol{\nabla }{}\mu \boldsymbol{\cdot }\boldsymbol{\nabla }{}y=0 \quad \forall {}Y ,\end{align}

with $Y$ as the test function in the reference configuration and $y(\boldsymbol{x})=Y(\boldsymbol{X})$ its representation in the current configuration; see Appendix C. The definition of the chemical potential in (2.9) is expressed in weak form as

(2.20) \begin{align} M(\boldsymbol{\cdot },t):\quad \int {\mathrm{d}}^2X\,{\textit{MZ}} -\partial _C\mathcal{E}[\boldsymbol{\psi },C,\varPi ](Z) =0 \quad \forall {}Z ,\end{align}

where $Z$ acts as test function. The weak formulation of the poro-visco-elastocapillary problem consists of the aggregation of (2.17)–(2.20), equipped with suitable initial conditions for $\boldsymbol{\psi }$ and $C$ .

2.4. Scales and dimensionless numbers

The free-energy formulation provides the natural scales for energies and lengths. We start by comparing the Flory–Huggins energy density, scaling with $kT/v$ , to the elastic energy density, scaling with $G=NkT$ . Its ratio defines a first dimensionless parameter, $Nv$ . Combined with the Flory–Huggins parameter $\chi$ , this determines the natural degree of swelling of the substrate when fully immersed. These parameters are the natural ingredients of the Flory–Rehner theory (Flory & Rehner Reference Flory and Rehner1943).

Then, the ratio of surface energy to bulk elastic energy provides a length scale $\gamma _s/G$ . This defines the so-called elasto-capillary length, which gives a typical measure of elastic deformation. Another length scale of the system is the layer thickness of the preswollen substrate, which we denote as $H$ . The ratio $\gamma _s/(GH)$ geometrically describes the importance of the deformation. A final dimensionless parameter derived from the free energy is obtained from the ratio of surface tensions $\gamma _{\textit{LV}}/\gamma _s$ . This parameter is known to determine the solid angle of the solid, according to Neumann’s law (Pandey et al. Reference Pandey, Andreotti, Karpitschka, Van Zwieten, van Brummelen and Snoeijer2020; Flapper et al. Reference Flapper, Pandey, Essink, Van Brummelen, Karpitschka and Snoeijer2023).

In summary, there are 4 dimensionless numbers governing the ‘statics’ of the problem. For tractability, we keep them at fixed typical values, namely

(2.21) \begin{align} Nv=1, \quad \chi =0, \quad \frac {\gamma _s}{\textit{GH}}=1, \quad \frac {\gamma _{\textit{LV}}}{\gamma _s}=1. \end{align}

Changing the parameter values will lead only to quantitative differences, in particular the degree of preswelling ( $Nv$ , $\chi$ ), the relative size of the elastic deformation ( $\gamma _s/GH$ ) and the opening angle of the ridge ( $\gamma _{\textit{LV}}/\gamma _s$ ). However, our main interest here is to explore the dynamical behaviour of the ridge, which is not expected to change qualitatively with these parameters.

The dynamical model is defined by (2.8) and (2.11), each of which introduces a natural time scale. The equation for solvent transport relates $c$ to $\mu$ , which have typical scales $1/v$ and $kT \sim Gv$ . To estimate the diffusion time, one can use either $\gamma _s/G$ (the ridge size) or $H$ (the substrate thickness) as the length scale. We opt for the former to make the equations dimensionless, which introduces the diffusion time scale inside the ridge

(2.22) \begin{align} \tau _{\textit{D}} = \frac {(\gamma _s/G)^2}{D v^2 G} = \frac {\gamma _s^2}{D v^2 G^3}. \end{align}

It should be noticed that the choice of $\gamma _s/G$ does not imply that it is the only relevant scale. In particular, for thin substrates the thickness $H$ can strongly influence the deformation, especially under fixed-bottom conditions.

In the presence of viscous effects per (2.13), the momentum equation contains a visco-elastic relaxation time scale $\tau _{ {\textit{VE}}} = \eta /G$ , introducing a fifth dimensionless number

(2.23) \begin{align} \tilde {\tau }=\frac {\tau _{\textit{VE}}}{\tau _D} = \frac {\tau _{\textit{VE}} D v^2 G^3}{\gamma _s^2}. \end{align}

This parameter governs the cross-over from poroelastic to viscoelastic response.

In the remainder of the paper we use tildes to denote dimensionless variables. We scale lengths with $\gamma _s/G$ , energies with $vG$ , while time is scaled with $\tau _D$ . Hence, results will be presented in terms of the following dimensionless variables:

(2.24) \begin{align} \tilde {\boldsymbol{x}}=\frac {\boldsymbol{x}}{\gamma _s/G}, \quad \tilde {t} = \frac {t}{\tau _D}, \quad \tilde {\mu }=\frac {\mu }{vG}, \quad \tilde {c}=v c. \end{align}

Note that the solvent concentration is made dimensionless with molecular volume (rather than $\gamma _s/G$ ), so that $\tilde{c}$ directly gives the liquid volume fraction.

2.5. Numerical methods

We approximate the poro-visco-elastocapillary problem (2.17)–(2.20) by means of the finite-element method in the spatial dependence, and the implicit Euler method in the temporal dependence. The approximation method has been implemented in the open-source software framework Nutils (van Zwieten, van Zwieten & Hoitinga Reference van Zwieten, van Zwieten and Hoitinga2022).

To provide adequate resolution at the tip of the wetting ridge, we apply a mesh with a priori local refinements near the contact line. Starting with a coarse baseline mesh, we successively subdivide the elements contained in a semi-disc with increasingly small radius, centred at the contact line; see figure 1(c). The discrete approximation is obtained by replacing the ambient spaces for the trial functions $\boldsymbol{\psi },\varPi ,C,\mu$ and the test functions $\boldsymbol{W},V,Y,Z$ by finite-element spaces subordinate to the locally refined mesh. In particular, we apply continuous piece-wise quadratic approximation spaces for $\boldsymbol{\psi },C,\mu$ and $\boldsymbol{W},Y,Z$ and continuous piece-wise linear approximation spaces for $\varPi$ and $V$ . This implies that for the deformation/pressure pair $(\boldsymbol{\psi },\varPi )$ we apply so-called Taylor–Hood elements. The functional derivatives that appear in the weak formulation (2.17)–(2.20) are obtained from implementations of the free-energy and work functionals, by means of the automatic differentiation functionality in Nutils.

The nonlinear algebraic problem that emerges in each time step is solved by means of a Newton–Raphson method with line search. The linear tangent problems that occur in each iteration of the Newton–Raphson method are solved by a direct solver. We validate the numerical results by comparing the interface profile at $ \tilde {t} = 0^+$ and $\tilde {t} = \infty$ with the equilibrium results presented in Pandey et al. (Reference Pandey, Andreotti, Karpitschka, Van Zwieten, van Brummelen and Snoeijer2020) and Flapper et al. (Reference Flapper, Pandey, Essink, Van Brummelen, Karpitschka and Snoeijer2023). Further details are provided in Appendix A, and the code is made available at Zheng et al. (Reference Zheng, Chan, van Brummelen and Snoeijer2025).

3.1. Ridge growth

To illustrate the evolution of the free surface of the poroelastic gel, figure 2(a) plots the vertical displacement of the free surface, $\tilde {u}_z(\tilde {x},\tilde {t})$ , versus the horizontal position in the deformed configuration, $\tilde {x}$ , at various instants, $\tilde {t}$ . At $\tilde t = 0$ , the polymeric gel surface is flat, corresponding to the homogeneously preswollen initial state. Upon application of the capillary force at the contact line, the gel exhibits an instantaneous response, as evidenced by the profile at $\tilde {t} = 10^{-8}$ . We refer to this very early time limit as $\tilde {t} = 0^+$ . This instantaneous response is essentially incompressible because the solvent within the gel did not yet have any time to diffuse, preventing any noticeable change in swelling immediately after the capillary force is applied; see also Zhao et al. (Reference Zhao, Lequeux, Narita, Roche, Limat and Dervaux2018b ). Note that an instantaneous response would in reality be prevented by inertia, but this involves extremely small time scales. (Inertia acts on a time scale that can be estimated from the elastocapillary length and the speed of elastic waves, as $ (\gamma _s / G ) / \sqrt { G / \rho } = \gamma _s \rho ^{1/2} / G^{3 / 2}$ , typically of the order of $10^{-5} \mathrm{\,s}$ for soft gels. Formally, this inertial transient precedes the viscoelastic and poroelastic regimes of interest, which is of the order of seconds or larger.) Subsequently, as time progresses, the gel absorbs solvent from the surroundings and swells, resulting in an increase in its volume, as indicated by the curve at $\tilde {t} = 10^{-1}$ and the profile at $\tilde {t}=\infty$ , in figure 2(a). The substrate eventually reaches an equilibrium state as $\tilde {t}\rightarrow \infty$ .

To further characterise the dynamics of the ridge growth, we show the height of the elastocapillary ridge, $ \tilde {h}(\tilde {t}):=\tilde {u}_z(0,\tilde {t})$ , as a function of time in figure 2(b). We indeed observe an instantaneous displacement at $\tilde {t}=0^+$ , which is followed by a gradual increase until equilibrium is reached. The inset of figure 2(b) presents the same data, but instead plots $\tilde {h}(\tilde {t}) - \tilde {h}(0^+)$ versus $\tilde {t}$ on a logarithmic scale. This graph shows that after the instantaneous incompressible response, the growth of the ridge follows a power law with an exponent of 1/2. In dimensional terms, this growth implies that $h(t) - h(0^+)\sim (D v^2 G t)^{1/2}$ , which is scale free in the sense that it does not involve the elasto-capillary length $\gamma _s/G$ nor the system size $H$ . This scale-free behaviour reflects that the growth is initially due to diffusive imbibition of the solvent into the substrate. After this scaling regime, the equilibration is observed around $\tilde t \sim 1$ , suggesting that the diffusion time scale $\tau _D$ , is indeed the appropriate poroelastic time scale for the late-time process, where the ridge size $\gamma _s/G$ begins to play a role.

3.2. Chemical potential

Solvent transport plays a fundamental role in the dynamics of the wetting-ridge formation. The chemical potential $\mu$ per (2.9) effectively controls the magnitude and direction of the solvent flux, because the flux is proportional to $-\boldsymbol{\nabla }\mu$ . Hence, the spatio-temporal evolution of $\mu$ provides meaningful insight into the underlying solvent-transport mechanisms.

Figure 3.

Spatio-temporal evolution of the chemical potential for a poroelastic substrate ( $\tilde{\tau} =0$ ). (a) The rescaled chemical potential field at $\tilde {t} = 2 \times 10^{-4}$ ; (b) A zoomed-in view near the contact line position from (a); (c) Measurement of $\tilde {\mu }$ from the tip along the $z$ -direction at different times. The circles indicating the local minima define the boundary layer thickness $\delta$ . The dashed line indicates the prediction (3.1) for the instantaneous incompressible response; (d) Boundary layer thickness $\delta$ as a function of time $\tilde {t}$ .

Figure 3(a) presents the rescaled chemical potential $\tilde{\mu}$ in the porous substrate, at $ \tilde {t} = 2 \times 10^{-4}$ . The figure conveys that the chemical potential displays a localised minimum just below the contact line, represented in the figure by the blue region. This minimum in the chemical potential can equivalently be interpreted as a region of low ‘pore pressure’, and originates from the tensile elastic stress generated by the deformed polymer network. The solvent thus tends to diffuse towards this blue region. A magnified view near the contact line is given in figure 3(b). This zoom shows that $\tilde{\mu} =0$ at the free surface, reflecting the boundary condition that the surface is in equilibrium with the reservoir. The blue region is separated from the free surface by a thin boundary layer, that gives rise to a strong local gradient of $\tilde{\mu}$ . Hence, the solvent flux will be strongest near the interface, so that the predominant flux is from the reservoir into the ridge.

Next, we focus on the temporal evolution of the chemical potential. Figure 3(c) shows $\tilde{\mu}$ along the vertical line below the contact line ( $\tilde{x}$ = 0), at different times. The horizontal axis $\tilde h - \tilde z$ represents the distance to the contact line. The graphs confirm the appearance of a minimum of $\tilde{\mu}$ , which shifts and fades over time. At very long times $\tilde{\mu} \to 0$ everywhere, corresponding to the global equilibration. Interestingly, the data for early times ( $\tilde {t} = 10^{-8}$ , shown in black) seem to follow a logarithmic dependence $\tilde{\mu} \sim \ln (\tilde h - \tilde z)$ . This logarithmic behaviour can be understood in terms of the instantaneous incompressible response expected at $t=0^+$ . At this early-time limit, the concentration field did not yet have time to evolve, and the chemical potential, as defined in (2.9), follows the behaviour of $\varPi$ . For an incompressible corner, $\varPi$ plays the role of the pressure that is known to exhibit a logarithmic dependence on the distance to the contact line (Pandey et al. Reference Pandey, Andreotti, Karpitschka, Van Zwieten, van Brummelen and Snoeijer2020). Translated in terms of the scaled chemical potential, this takes the form

(3.1) \begin{align} \tilde {\mu } \simeq \left (\frac {\pi }{\theta _{\text{S}}}-\frac {\theta _{\text{S}}}{\pi }\right ) \ln ( \tilde r), \end{align}

where $\tilde{r} $ is the distance to the contact line and $\theta _{\text{S}}$ is the corner angle ( $\theta _{\text{S}}=2\pi /3$ when $\gamma _{\textit{LV}}/\gamma _s=1$ ). This prediction is superimposed as the grey dashed line in figure 3(c), and indeed perfectly captures the observed behaviour in the vicinity of the contact line at $\tilde {t}=10^{-8}$ . Note that the logarithmic divergence of $\tilde{\mu}$ is not specific to the neo-Hookean solid; a different constitutive relation only modifies the prefactor in (3.1).

The picture that emerges is that the instantaneous response dictates that the chemical potential diverges logarithmically, according to (3.1), in the vicinity of the contact line. This behaviour, however, is incompatible with the boundary condition $\mu =0$ that is imposed at the free surface. This conundrum is resolved by the emergence of a thin boundary layer near the contact line. Indeed, the blue and red curves in figure 3(c) closely follow the near-instantaneous response (black line) away from the contact line. At shorter distances, the logarithmic trend is not followed, and the data tend to $\tilde{\mu} =0$ at the contact line. We estimate the size of the boundary layer $\delta$ by tracing the minimum of the chemical potential, indicated by the circles. Figure 3(d) reports the boundary layer thickness over time, which again follows a diffusive $1/2$ scaling. Hence, we associate the dynamics of the boundary layer with the scale-free diffusive solvent transport.

3.3. Swelling dynamics

From the previous results we thus conclude that the incompressible response gives an excellent representation of the behaviour at early times, except for a narrow boundary layer $\delta \sim t^{1/2}$ . Swelling does occur in this boundary layer, which is characterised by a strong flux $\sim 1/\delta \sim 1/t^{1/2}$ . Now we turn to the question of how this translates to the spatio-temporal dynamics of the swelling inside the substrate.

Figure 4.

Swelling dynamics for a poroelastic substrate ( $\tilde{\tau} =0$ ). (a) Measurement of $J$ from the tip along the $z$ -direction at different times. The preswelling at $t=0$ corresponds to an initial $J_0=1.582$ . (b) Solvent volume fraction field $\tilde{c}$ at the final equilibrium state. Near the tip the solvent fraction $\tilde c \to 1$ , while at large distance one recovers the value due to preswelling $\tilde c \approx 0.37$ .

Figure 4(a) shows the swelling ratio $J$ as a function of the vertical distance $\tilde h - \tilde z$ , measured directly below the contact line. The initial value corresponds to the preswollen substrate with $J_0=1.582$ . In line with the evolution of the boundary layer for $\tilde{\mu}$ , we observe that swelling is initiated close to the contact line, and gradually invades the entire domain. The long-time equilibrium profile is given by the yellow line. The strongest swelling occurs at the tip of the wetting ridge, and exhibits a power-law singularity. Indeed, previous work on the equilibrium state (Flapper et al. Reference Flapper, Pandey, Essink, Van Brummelen, Karpitschka and Snoeijer2023) has analytically demonstrated that, for the neo-Hookean solid, the swelling ratio $J \sim \tilde r^\beta$ , where the exponent $\beta = 2(1 - \pi /\theta _{{S}})$ ; in this expression $\theta _{\text{S}}$ once again is the opening angle of the tip of the solid. The numerical results closely follow this analytical prediction, also transiently in the vicinity of the contact line.

The divergence of $J$ implies that the tip region is almost purely liquid. Specifically, the scaled concentration $\tilde c = 1 - 1/J$ directly gives the liquid volume fraction. The volume fraction at equilibrium is shown in figure 4(b), and indeed approaches unity near the tip. In the present formulation the liquid cannot phase separate between the matrix and the liquid reservoir, as observed in some experiments (Jensen et al. Reference Jensen, Sarfati, Style, Boltyanskiy, Chakrabarti, Chaudhury and Dufresne2015; Cai et al. Reference Cai, Skabeev, Morozova and Pham2021; Hauer et al. Reference Hauer, Cai, Skabeev, Vollmer and Pham2023). It is clear, however, that the strong aspiration of liquid in elastic corners provides a natural mechanism for this phenomenon.

4.1. Ridge growth

Figure 5 depicts a typical visco-poroelastic evolution of the free surface for the case where $\tilde{\tau} =10^{-3}$ . This value of $\tilde{\tau}$ implies that the network relaxation is faster compared with the diffusive transport of the solvent. Nonetheless, we observe that the growth of the ridge is fundamentally altered with respect to the purely poroelastic case: the instantaneous incompressible response is now replaced by a gradual viscoelastic growth of the wetting ridge, occurring on the time scale $\tilde t \lt \tilde{\tau}$ . The ridge emerges gradually in this regime, both in terms of its height and its width. It is noteworthy that the solid angle of the ridge forms instantaneously, in accordance with the Neumann law (2.15), which is implicit in (2.17) as a natural (weakly imposed) condition; see also Karpitschka et al. (Reference Karpitschka, Das, van Gorcum, Perrin, Andreotti and Snoeijer2015) and Chan (Reference Chan2022). We note that, after $\tilde t \approx \tilde{\tau} =10^{-3}$ (dark blue), the profiles closely follow those in figure 2(a). The limiting transport mechanism during this late-time dynamics is the slow poroelastic solvent flow, not the viscoelastic relaxation.

Figure 6.

Growth of the ridge for visco-poroelastic substrates. (a) The deformation $\tilde {h}$ of the tip as a function of time $\tilde {t}$ for various viscoelastic relaxation times $\tilde {\tau }$ . (b) Same data on a double logarithmic scale, showing the linear ridge growth when viscoelasticity is included.

Figure 7.

Swelling ratio (a) and chemical potential (b), for a visco-poroelastic substrate. Profiles along the $z$ -direction from the tip. The solid lines correspond to a purely poroelastic material ( $ \tilde{\tau} =0$ ). The dashed lines correspond to a visco-poroelastic material with $\tilde {\tau } = 10^{-3}$ .

To further quantify the ridge dynamics, we report in figure 6(a) the height of the ridge $\tilde {h}$ as a function of time $\tilde {t}$ in a semi-log plot, for various viscoelastic relaxation times $\tilde {\tau }$ . When $\tilde {\tau }= 0$ (represented by the black line), we recover the purely poroelastic case that we previously observed to exhibit an instantaneous response. Conversely, when $\tilde {\tau }$ is non-zero, one no longer observes this instantaneous deformation. Instead, the tip height grows gradually until joining the purely poroelastic curve approximately when $\tilde t \sim \tilde{\tau}$ . These results confirm the cross-over, from viscoelastic at early times to poroelastic at later times. Figure 6(b) shows the same data on a doubly logarithmic plot. We find a linear initial growth in the viscoelastic regime. The scaling follows $\tilde h \sim \tilde t/\tilde{\tau}$ , which in dimensional units amounts to $h \sim \gamma _s t/(G\tau _{\textit{VE}})$ . Bearing in mind that the network viscosity $\eta = G \tau _{\textit{VE}}$ , the growth of the ridge is dictated by the capillary velocity $\gamma _s/\eta$ . Hence, the linear ridge growth here arises from the Newtonian rheology employed in the viscoelastic term of the model. Note that if the network instead followed a nonlinear relaxation (e.g. a power-law rheology), the early-time scaling of the ridge height would be modified accordingly, and take over the rheological exponent.

4.2. Swelling dynamics

Next, we turn to the transport of solvent in the presence of viscoelasticity. Figure 7(a) shows the swelling ratio below the contact line, at different times. The dash-dotted lines represent the results for the case where $\tilde {\tau } = 10^{-3}$ , which was previously considered in figure 5. To clearly illustrate how viscoelasticity changes the transport dynamics, we also report the purely poroelastic results for $\tilde {\tau } = 0$ , as solid lines. For times $\tilde t \gt \tilde{\tau}$ , the solid lines and dash-dotted lines essentially overlap. This implies that the late-time swelling dynamics is not affected by viscoelasticity, but is purely poroelastic in nature. By contrast, large differences are observed for the early time dynamics, $\tilde t \lt \tilde{\tau}$ . There is still a tendency towards swelling near the tip, even at very short times, but this swelling is retarded compared with the purely poroelastic case. The suppression of solvent transport at early times can be attributed to the vanishing of logarithmic divergence in the chemical potential, as shown in figure 7(b). In other words, the viscoelastic stresses inhibit the development of strong ‘pore pressure’. Only when $\tilde t \geq \tilde{\tau}$ do the curves for $\tilde{\tau} =0$ and $\tilde{\tau} \neq 0$ start to overlap, indicating the regime where transport is purely poroelastic.

Figure 8.

Regime map illustrating the interplay between the viscoelastic and poroelastic dynamics. For a given $\tilde{\tau}$ , the initial response is viscoelastic followed by a poroelastic regime. In the latter regime, the response can initially be incompressible depending on the dimensionless distance to the contact line ( $\tilde d = dG/\gamma _s)$ .

To summarise the dynamical regimes, figure 8 presents a schematic diagram in terms of the normalised time $\tilde {t}=t / \tau _D$ and the ratio $\tilde {\tau }=\tau _{\textit{VE}} / \tau _D$ . The temporal evolution of the simulation crosses the different regimes in the diagram along a horizontal line, at a given value of $\tilde{\tau}$ . The gel initially behaves as an incompressible viscoelastic solid, since no solvent influx occurs in regions far from the tip. As time progresses and solvent diffusion penetrates the bulk, the substrate gradually enters the poroelastic regime, leading to swelling at larger scales. At a distance $d$ from the contact line, the onset of swelling is governed by the diffusive time $\tilde t \sim (dG/\gamma _s )^2$ . The cross-over line $\tilde t = \tilde {\tau }$ delineates the transition between viscoelastic- and poroelastic-dominated behaviour.