Droplet absorption and spreading into thin layers of polymer hydrogels

Abstract

From biological tissues to microactuators and absorption of solvents into layers of paint, macroscopically non-porous materials with the capacity to swell when in contact with a solvent are ubiquitous. In these systems, owing to strong solid–fluid interactions, chemically driven flows can yield large geometric changes. We study experimentally and theoretically the canonical problem of the swelling of a thin hydrogel layer by a single water drop. Using a bespoke experimental set-up, we observe fast absorption leading to a radially spreading axisymmetric blister. We use a fully three-dimensional linear poroelastic framework with nonlinear kinematic equations to obtain governing equations, which we then reduce with thin-layer scalings to a one-dimensional nonlinear diffusion equation for the evolution of the blister geometry. In the limits of large and small deformations, the evolution of the blister characteristic height and radius are self-similar, following power laws in time. Our experimental measurements show that the evolution of the blister is broadly within the theoretical predictions in the large and small deformation regimes. In the general intermediate deformation regime, the measurements are well captured by our reduced one-dimensional diffusion model, which does not require the sophisticated and computationally expensive numerical approaches necessary for the original two-dimensional nonlinear coupled transport problem. By adapting modelling techniques from the fluid dynamics of thin porous elastic layers to a polymer swelling problem, our modelling framework extends the range of polymer swelling problems that can be treated with semi-analytical methods. Moreover, our detailed experimental data can serve as a test case for future nonlinear poroelastic frameworks of swelling polymer materials.

1. Introduction

High-swelling polymer systems, such as hydrogels (a class of hydrophilic cross-linked polymers that form three-dimensional networks of macroscopic extent in water), have received considerable attention recently (Vervoort 2006; Laftah, Hashim & Ibrahim 2011). Large volumetric stresses and strains, induced by swelling in response to changing environmental conditions, lead to a fascinating class of problems, where chemically driven transport within the polymer induces large changes in the geometry of the hydrogel (e.g. Guvendiren, Burdick & Yang 2010; Etzold, Linden & Worster 2021). These systems largely behave as hyper-elastic incompressible solids on short time scales (Yoon et al. 2010). On longer time scales, hydrogels can respond to changing external conditions (stress, temperature, humidity) by reversible swelling, shrinking or deforming. Since the mechanical properties of hydrogels can be adjusted over a broad range by changing the polymer chemistry and cross-linker density (Ozcelik 2016; Tortorella et al. 2021), they have found widespread applications.

The ability of external stimuli to trigger swelling that drives large deformations has been harnessed to drive microfluidic pumps (Richter et al. 2009; Kwon et al. 2011) and has led to their use as micro-actuators for optical, flow control, sensor and microrobotics applications (Porter et al. 2007; Ionov 2014). Understanding the interaction between the swelling dynamics and the microscopic driving forces is crucial for the design of these systems.

Hydrogels can absorb large volumes of water (up to thousands of times their dry weight). This, and the fact that hydrogels have similar physico-chemical characteristics to many tissues (Drury & Mooney 2003; Hoffman 2012), makes hydrogels an excellent substitute for tissues in the laboratory (e.g. when studying decompression sickness; Walsh et al. 2017; Zhang, Etzold & Lefauve 2021). Indeed, for over sixty years, since the pioneering work of Wichterle & Lim (1960) who developed the first contact lens using a polyhydroxyethylmethacrylate bio-compatible hydrogel, hydrogels have been extensively used across many medical applications (Hoffman 2012). In tissue engineering, they are used as a matrix framework for the repair and regeneration of a wide range of tissues and organs. Lee & Mooney (2001) give a general review of all the applications, while Tang et al. (2020) and Madhusudanan, Raju & Shankarappa (2020) present applications in the spinal column and brain, respectively. Furthermore, for wound dressing, hydrogels maintain a moist healing environment whilst allowing gaseous exchange, thereby promoting wound healing (Stubbe et al. 2021; Tortorella et al. 2021; Deng et al. 2022).

These medical systems, characterised by complex interactions between biological substances and man-made hydrogels, require a thorough understanding of hydrogel swelling behaviour. An example of this interplay lies in the work of Cont et al. (2020), who explore experimentally how Vibrio cholerae biofilms deform both thin hydrogel layers and epithelial cell mono-layers attached to the surface of a soft extracellular matrix. The interactions between swelling and mechanical stresses cause the layers to buckle into the biofilm and often break up, thus allowing the biofilm to compromise the physiology of its host.

A range of fully three-dimensional poroelastic theories for hydrogel swelling have been proposed. These theories derive a constitutive equation for the stress tensor from assumptions about the strain energy density function. While linear theories assume this density function to be quadratic in the strain (Yoon et al. 2010), nonlinear theories (e.g. Hong et al. 2008; Doi 2009; Chester & Anand 2010) construct it from thermodynamic models describing polymer deformation and polymer–solvent mixing (Cai & Suo 2012). These models are equivalent in the limit of small deformation (Bouklas & Huang 2012). Nevertheless, most experimental and modelling studies have only considered one-dimensional or quasi-one-dimensional geometries such as the swelling of fibres (Van de Velde, Protière & Duprat 2021; Van de Velde et al. 2022), spheres (Tanaka & Fillmore 1979; Engelsberg & Barros 2013; Bertrand et al. 2016; Butler & Montenegro-Johnson 2022) or flat sheets bound to surfaces (Tanaka & Fillmore 1979). In other cases, finite-element methods are used to simulate more complicated problems such as the swelling caused by a succession of falling droplets (Phadnis et al. 2018). However, considerable challenges remain, especially for high-swelling systems producing large swelling gradients (Yu, Malakpoor & Huyghe 2020).

In this article, inspired by the flat geometries found in cell layers (Cont et al. 2020), wound dressings (Tortorella et al. 2021), hydrogel water harvesters (Li et al. 2018) and layers of paint (Varady et al. 2016), we investigate in detail the temporal dynamics of the absorption into and subsequent spreading of a water drop within a thin absorbent polymer layer of hydrogel. In particular, an axisymmetric bulge forms at the hydrogel surface, which we will henceforth denote as a blister.

In § 2 we describe experiments tracking the absorption of a single droplet of water by a thin layer of hydrogel. Particular care is taken to define the characteristic properties of the blister that emerges on the surface of the hydrogel, such as the characteristic radius of the blister (see § 2.4). Error bars are computed through a careful uncertainty analysis (see Appendix A for further details). The results of these experiments are described in § 2.5. The water first absorbs into the hydrogel, forming a blister. A reversible transient buckling/crumpling instability then forms on the surface of the blister. The blister slowly spreads outwards, eventually approaching what appears qualitatively to be a self-similar shape for this nonlinear diffusion process.

The requisite next step is to develop a theoretical model to probe the underlying physical processes. In § 3, we analyse this problem using a poroelastic framework that combines nonlinear conservation of volume, Darcy's law with a permeability dependent on polymer volume fraction and Biot's linearised equation for the stress tensor. Motivated by the apparent self-similar profile observed in the last phase of the experiments, we develop a mathematical model for the swelling hydrogel layer. A poroelastic framework is used to obtain a nonlinear diffusion equation that describes the resulting spreading dynamics of the system. Regarding the physical processes at the heart of our model, we consider a system in which the dynamics is governed by the interplay between the motion of both polymer and water and the balance of elastic stresses through the permeable polymer. Classical shallow-layer scalings reduce this axisymmetric three-dimensional problem to a one-dimensional nonlinear diffusion model. This model, which admits analytic similarity solutions in both the small (§ 4.1) and large deformation limits (§ 4.2), is studied numerically in § 4.3, using the finite-element software package FEniCS (Logg, Mardal & Wells 2012; Alnæs et al. 2015). The regimes of the parameter space in which the analytic solutions agree well with the numerical solutions are explored and discussed.

We demonstrate the effectiveness of our two-parameter one-dimensional nonlinear diffusion model by fitting numerical solutions of the model to the experimental results. The numerical predictions agree very well with the experimental data over the entire range of parameters considered. Our model can quantitatively predict the evolution of all observable quantities over hours of observation time. Finally, conclusions and possible next steps are detailed in § 6.

2. Experimental

Here, we describe a series of experiments that considered the absorption of single droplets of water with volumes between $5$ and $100\,{\mathrm {\mu }}{\rm l}$ by a thin layer of hydrogel, leading to the formation and then subsequent outward spreading of a blister. This configuration is analogous to the injection of a finite volume of fluid into a thin poroelastic layer, a well-studied canonical fluid dynamics problem (e.g. Hewitt, Neufeld & Balmforth 2015).

2.1. Materials and methods

Illustrated by two photographs (figure 1a,b) and a schematic (figure 1c), our experiments utilised commercially available medical-grade hydrogel pads (Hydrogel Nipple Pads, Medela, Switzerland), which were manufactured as approximately $8\times 8\,{\rm cm}^{2}$ sheets. Unfortunately, their proprietary nature precluded us from obtaining details on the synthesis, in particular the polymer/monomer fraction at synthesis. Upon delivery, the hydrogel pads appear as a rubbery sheet with a slightly tacky surface. After prolonged exposure to the atmosphere (ca. 25 %–40 % relative humidity, ca. $22\,^{\circ }{\rm C}$), their initial thickness $a$ was determined with a micrometer as $a=1.13\,{\rm mm}\pm 0.01\,{\rm mm}$. One side of a sheet was affixed to a thin plastic sheet (impermeable to water) of approximate thickness $0.04\,{\rm mm}$. The other side was initially protected by an easily removable plastic film. The hydrogel sheets incorporated a thin non-woven gauze during their manufacture. We assume that these do not affect swelling due to the lateral constraint provided by the bottom plastic sheet, which leads to predominantly upwards swelling. In our experiments, these larger sheets were cut into coupons approximately $25\times 25\,{\rm mm}^{2}$. The coupons were then glued via the thin impermeable plastic sheet onto $76\times 26\,{\rm mm}^{2}$ microscope slides (Menzel, Germany) with UV-cured Norland Optical Adhesive NOA 68 optical glue (Thorlabs, USA). The gluing of the coupon to the microscope slide ensured that the gauze did not affect the swelling of the sheets.

Figure 1. Laboratory experiments investigating the absorption and subsequent spreading of water into a hydrogel sheet. (a) Photograph of the preliminary experimental set-up without oil in a cell sealed from the atmosphere with an evaporation barrier. (b) Photograph of the full experimental set-up with an oil bath. In both cases, the side view camera is positioned outside the region of the apparatus captured by the photograph, pointing in the direction given by the red arrow with label camera. (c) Schematic of the experimental apparatus.

Utilising Karl-Fischer titration, a standard volumetric titration method to determine trace amounts of water in a sample, the water content of a typical hydrogel pad that had equilibrated with the laboratory air was measured to be 22 %. Polymer swelling is driven by entropic and pairwise interactions between solvent and polymer molecules and charged groups if these are present (Flory 1953). We tested for the presence of ionic groups by swelling tests in solvents that are less able to support the dissociation of the counter-ion. We found that isopropanol and ethanol did not swell the hydrogel perceptibly and, when placed in saturated aqueous sodium chloride solutions, the swelling was drastically reduced. We take these observations as evidence for the presence of charged groups within the hydrogel.

2.2. Suppression of vapour transport

In preliminary experiments, water droplets were placed onto mounted hydrogel that was otherwise exposed to air but enclosed into a sealed cell to prevent evaporation of the droplet into the broader atmosphere (see figure 1a). The edges of the hydrogel sheet were observed to swell perceptibly with this swelling beginning immediately after droplet placement. This was attributed to evaporation from the regions in the vicinity of the droplet, namely the most strongly swollen regions, leading to higher humidity and lateral transport of humidity in the sealed cell and thus reabsorption by the drier regions of the hydrogel at the edge of the sheet. This was confirmed by experiments in which water droplets were placed next to the hydrogel on separate microscope slides within the sealed cell. Swelling was observed that could only have come from vapour-phase transport.

As shown in figure 1(b,c), to eliminate vapour transport experimentally, the hydrogel was immersed in a water-immiscible oil. We chose HySpin AWS 32 (CASTROL) for its low solubility with water. We remark that the presence of the oil can modify the interfacial energies of the system compared with a droplet in ambient air, and thus wetting and contact angle properties. However, this should only influence the early stages of the experiments, when the water drop is in the process of being absorbed in the hydrogel. Since our study focuses on the later stages of the experiment, once the water has been fully absorbed into the hydrogel, the impact of using oil instead of air on the late-time transport dynamics in the hydrogel should be small.

It is important to recognise that even apparently immiscible liquids are sparingly soluble into each other. For the hydraulic oil, the water content was found utilising Karl-Fisher titration to be 6700 ppm after equilibration with the laboratory air. We also conducted a titration with oil that had been equilibrated with a small amount of liquid water, indicating an estimate for the saturated water content of the oil to be up to 1.1 %. The later value is likely to be an overestimate due to the possibility of an emulsion forming during the shipping of the oil–water mixture to the place of analysis.

To ensure that the oil and the hydrogel were in equilibrium during an experiment, they were both exposed to the same atmosphere before the experiment for at least two days. The oil was not changed between experiments and thus was always exposed to the atmosphere. Furthermore, the hydrogel sheets were allowed to equilibrate for at least 1–3 h within the oil before droplet placement. The oil was not stirred and disturbances were kept to a minimum to avoid introducing air bubbles in the system. Equilibration of the hydrogel with the oil was verified through careful tracking of the interface. The interface was stable, namely the oil did not noticeably absorb into the hydrogel. In particular, the swelling of the edges observed when the oil was not present did not occur any more.

2.3. Imaging system and data analysis

A monochrome camera (SP-5000M-CXP2, Jai, Denmark), with a variable magnification telecentric lens (TEV0305, VS Technology Corporation, Japan, set at 0.4$\times$ magnification) aligned with the plane of the hydrogel sheet, recorded the evolution of our samples. Background lighting was provided by a white LED light sheet (MiniSun Light Pad 17031), masked with black tape and placed approximately 4 m away from the hydrogel to ensure approximately collimated illumination.

DigiFlow was used for camera control (Dalziel 2017). Typical images are shown in figure 2 and were analysed using a Canny edge detector from the Python library Scikit Image Library (van der Walt et al. 2014). Suitable parameters for the edge detector (radius of Gaussian filter, upper and lower gradient thresholds for edge detection) were manually set for each dataset using an interactive interface. To mitigate edge effects arising from the finite size of the hydrogel sheet, we truncated each time series when the swelling from the blister approached the edges. For each dataset, we extracted the position of the upper edge of the unswollen hydrogel before the droplet was placed (see figure 2a) and the upper edge of hydrogel and droplet once the droplet was placed (see figure 2c–i). For analysis, the blister height (vertical displacement) was defined as

(2.1)\begin{equation} h_a = h(x,t)-a, \end{equation}

recalling that $h$ is the height of the hydrogel sheet, $a$ is the initial undeformed hydrogel sheet thickness and $0\leq x< x_{max}$, where $x_{max}$ is the right-hand boundary of the cropped image and the origin is the centre of the blister. This was then computed for each column of pixels. SciPy's Savitzky–Golay filter (window length $100\,{\rm pixel}$ (corresponds to 1.26 mm), polynomial order 1) was used to smooth the raw $h_a$ data (Virtanen et al. 2020). The total height of the hydrogel $h(x,t)$ was thus determined by adding the initial thickness of the hydrogel sheet (determined by a micrometer). This procedure was adopted since the lower edge of the hydrogel was not visible in the images.

Figure 2. Sequence of photographs showing the absorption of a $100\,{\mathrm {\mu }}{\rm l}$ droplet into the hydrogel sheet.(a–d) States with liquid water remaining. (e, f) A surface instability forms. ( f–h) Transition towards the long-time spreading regime which is shown in (h–l). The visible hydrogel sheet is marked blue in (a). From (c) onward a red line marks the hydrogel–oil interface.

2.4. Derived quantities

For each experimentally obtained profile $h(x,t)$, we found the centre position of the blister $x = r_0$ by fitting a Gaussian curve and taking $r_0$ to be the maximum of this curve. A Gaussian curve was selected instead of, e.g. a quadratic profile, since this is the expected form in the limit of linear diffusion. We define the corresponding characteristic radius of this blister profile $R(t)$ as the radius at which the blister height $h_a(r+r_0,t)$ has decreased from the maximum height of the blister $\max (h_a(r+r_0,t))$ to half of this value, thus,

(2.2) \begin{equation} R = \{r : h_a(r+r_0) = \tfrac{1}{2}(\underset{x}{\max}(h_a(x,t))) \}.\end{equation}

Due to the symmetry of the experimental profile about $r = r_0$, this leads to two values, $R_+$ and $R_-$, with $R_{-} \leq 0 \leq R_{+}$, which we chose to average (i.e. $R=(R_+ + |R_-|)/2$). Note that this approach implicitly assumes axisymmetry of the blister since we calculate $R$ from a single cross-section. We also compute the volume of the absorbed water using

(2.3)\begin{equation} V = {\rm \pi}\int_0^{x_{max}} (x-r_0) h_a(x)\, {{\rm d}\kern0.7pt x},\end{equation}

where the integral is evaluated over the full horizontal region of interest of the image (from $x = 0$ to $x = x_{max}$). As detailed in Appendix A, the uncertainty is dominated by the degree of alignment between the optical axis of the camera and the surface of the hydrogel sheet. This impacts most significantly our estimate of the radius $R$ where the error is

(2.4)\begin{equation} \delta R \approx \frac{\delta h_a}{\left(\dfrac{\partial h_a}{\partial x}\right)} {\bigg{\vert}_{x=r_0+R}}. \end{equation}

Here, $\delta h_a=(+22.7,-12.5)\, \mathrm {\mu }{\rm m}$ is the uncertainty in $h_a$ due to the pixel resolution and misalignment, respectively.

This also gives rise to an uncertainty in the volume, estimated as

(2.5)\begin{equation} \delta V \approx {\rm \pi}\int_{x_{mn}}^{x_{mx}} (x-r_0) \delta h_a(x) \,\mathrm{d}\kern0.06em x, \end{equation}

where the integration bounds $[x_{mn}, x_{mx}]$ are defined such that

(2.6)\begin{equation} \left\{x_{mn}\leq x\leq x_{mx}: h_a(x-r_0)>\delta z \right\}. \end{equation}

2.5. Observations

A montage of typical experimental images for a $100\,{\mathrm {\mu }}{\rm l}$ droplet is shown in figure 2. Refraction at the hydrogel–oil interface causes the swollen region to appear dark since the telecentric lens only collects the nearly collimated light directly from the light sheet. When released, the water droplet descended through the hydraulic oil to reach the hydrogel surface within a few seconds. Upon impact, the droplet retained its sphericity for approximately 30 s before relaxing to a sessile drop state (figure 2a). The subsequent hydrogel swelling dynamics had three distinct phases.

In the first phase, which lasted approximately 10 min once the water had reached a sessile drop state on the hydrogel surface, the surface water and swollen hydrogel coexisted as the sessile droplet transformed into a swollen blister (see figure 2a–c). This transition was particularly apparent at the droplet edge, which became progressively steeper (figure 2c) before decreasing again.

In the second phase, once all the surface water had been absorbed, a surface instability appeared, forming a blister with a crinkled surface (figure 2e–g). This buckling/crumpling instability arises in hydrogels that have strong swelling gradients or are laterally confined (Onuki 1988; Bouklas, Landis & Huang 2015) where the buckling relieves stress caused by mechanical connection with less swollen layers. Supplementary experiments in air (not shown here), where the water was removed at various times during the first phase, indicated that this instability with wavelength increasing with time was present once a thin swollen layer had been formed after droplet placement. Experimental images illustrating this instability are presented in Appendix B. We therefore conclude that the visible silhouette given in figure 2(a–d) consists of a swollen layer with buckling instability covered by water.

In the third phase, which lasted until the experiments were terminated once the blister front had reached the edge of the hydrogel or the resolution limit of the camera, the buckling/crumpling instability vanished. During this period, the blister assumed an apparent self-similar shape that resembled a nonlinear diffusion process.

Under the assumption that the blister spreads axisymmetrically, figure 3 plots raw data for the temporal evolution of the blister radius $R$ for droplet volumes 5, 10, 25, 30, 50 and 100 ${\mathrm {\mu }}{\rm l}$ (darker colours denote larger droplets). The data contain a full experiment from the moment that the droplet makes contact with the hydrogel surface until the point that meaningful observations were no longer possible. Initially, $R$ sharply increases as the water droplet, having made contact, spreads over the surface of the hydrogel (see figure 2a). Then $R$ decreases slightly, as the morphology of the blister changes from being relatively flat capped (caused by the horizontal redistribution of water above the droplet as it is absorbed at a rate that is nearly constant over the entire droplet radius) towards a smoother more self-similar shape (see the transition from figure 2b–h). Finally, we reach the self-similar regime where the radius of the swollen region increases monotonically (see figure 2i–l).

Figure 3. Raw data with corresponding error bars showing for a range of experiments how the swollen region radius varies as a function of time. Darker blue curves denote experiments with larger water droplets (see text for drop volumes). Error bars are generated through an uncertainty analysis that is described in § 2.4 and Appendix A.

Since conservation of water volume is assumed in § 3, we probed the experimental profile data for conservation of volume, as detailed in Appendix C.1. Overall the experimental uncertainties (Appendix A) were too large to reach definite conclusions. Furthermore, particularly for larger droplets, an error in the apparent volume could have arisen from the fact that the initial blister profile was not completely axisymmetric, resulting in spreading that was not quite radial. However, the main physical process which probably led to water loss was diffusion into the oil layer. Water loss from diffusion into the oil layer was estimated in Appendix C.2, enabling us to conclude that no significant water loss occurred during our experiments.

3. Poroelastic model

We now develop a mathematical model for the swelling hydrogel layer. A poroelastic framework is used to obtain a nonlinear diffusion equation that describes the resulting spreading dynamics of the system. We consider a system in which the dynamics is governed by the interplay between the motion of both polymer and water and the balance of elastic stresses through the polymer network.

3.1. Problem description

Figure 4 sketches the surface of a hydrogel layer of undeformed thickness $a$ that is fixed to a rigid horizontal boundary at $z = 0$. The gel is perturbed through the addition of a spherical water drop of diameter $d$, forming an axisymmetric blister of radius $R(t)$ and height profile $h(r,t)$ with $0\leq r<\infty$ and $0\leq t<\infty$. For an arbitrary function $f = f(z)$, we define the vertically averaging operator $\langle {\cdot } \rangle$ where $\langle f \rangle = ( \int ^{h}_{0} f \,{\rm d}z ) / h$. We examine an idealised hydrogel composition, namely a solution of polymer (volume fraction $\phi _p$, where $\phi _p = \phi _{0p}$ (cf. § 2.1 for experimental estimate of $\phi _{0p}$) at the start of the experiment before the water drop is added) and water (modelled as a Newtonian fluid with dynamic viscosity $\mu _f$ and volume fraction $1-\phi _p$). We denote the pore-averaged velocity and Cauchy stress tensor of the solid and liquid phases by $\{ \boldsymbol {u_p} = (u_p, w_p) , \boldsymbol {\sigma _{p}} \}$ and $\{ \boldsymbol {u_f} = (u_f, w_f) , \boldsymbol {\sigma _{f}} \approx -p\boldsymbol {I} \}$, respectively, where $p$ is the pore pressure and $\boldsymbol {I}$ the identity tensor.

Figure 4. Schematic of a hydrogel sheet. An axisymmetric blister of thickness $h(r,t)$ and characteristic radius $R(t)$ spreads out on the surface of a hydrogel sheet of undeformed thickness $a$, fixed to a rigid horizontal boundary at $z = 0$. The red dashed line highlights how $R(t)$ is calculated from experimental or numerical profile data using (2.2) or (4.33), respectively. Inset: the hydrogel is a solution of polymer (volume fraction $\phi _p$) and water (volume fraction $\phi _f = 1 - \phi _p$). The pore-averaged velocities of the solid and fluid phases are denoted by $\boldsymbol {u_p} = (u_p, w_p)$ and $\boldsymbol {u_f} = (u_f, w_f)$ respectively.

3.2. Governing equations

Conserving mass in both the solid and fluid phases gives

(3.1a)$$\begin{gather} \frac{\partial \phi_p}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\phi_p \boldsymbol{u_p}) = 0, \end{gather}$$

(3.1b)$$\begin{gather}-\frac{\partial \phi_p}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} ((1 - \phi_p)\boldsymbol{u_f}) = 0. \end{gather}$$

Defining the Terzaghi effective stress tensor as $\boldsymbol {\sigma } = \phi _p(\boldsymbol {\sigma _p} - \boldsymbol {\sigma _f})$ (Wang 2001), a momentum balance yields

(3.2) \begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot}(\phi \boldsymbol{\sigma}_p + (1 - \phi) \boldsymbol{\sigma}_f ) = \boldsymbol{\nabla} \boldsymbol{\cdot}(\phi \boldsymbol{\sigma}_p - (1 - \phi)p\boldsymbol{I}) = 0 \quad \Longrightarrow \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\sigma} = \boldsymbol{\nabla}p. \end{equation}

In general, $\boldsymbol {\sigma }$ obeys an elastic constitutive law of the form

(3.3)\begin{equation} \boldsymbol{\sigma} = \boldsymbol{\sigma} \left( \boldsymbol{\nabla} \boldsymbol{\xi} \right), \end{equation}

where $\boldsymbol {\xi } = (\xi, \zeta )$ is the two-dimensional deformation vector of the medium away from a reference state. The reference state that we chose was the hydrogel just before the droplet has been placed. The deformation vector is related to the velocity of the polymer phase through the material derivative

(3.4)\begin{equation} \boldsymbol{u_p} = \left(\frac{\partial}{\partial t} + \boldsymbol{u_p} \boldsymbol{\cdot} \boldsymbol{\nabla}\right)\boldsymbol{\xi}. \end{equation}

This is a common approach for many problems in biological physics (from lubrication of joints, Jensen et al. 1994, to cell cytoplasm, Charrras, Mitchison & Mahadevan 2009) and geophysics (from hydrology subsidence and pumping problems, Gibson, Schiffman & Pu 1970; Hewitt et al. 2015, to industrial filtration, Barry, Mercer & Zoppou 1997). Here, we consider the simplest case, namely that $\boldsymbol {\sigma }$ obeys a linear elastic constitutive equation of the form

(3.5)\begin{equation} \boldsymbol{\sigma}(\boldsymbol{\nabla} \boldsymbol{\xi}) = \left(K - \frac{2G}{3}\right)(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\xi})\boldsymbol{I} + G(\boldsymbol{\nabla} \boldsymbol{\xi} + \boldsymbol{\nabla} \boldsymbol{\xi}^{\rm T}), \end{equation}

where $K$ and $G$ are the osmotic and shear moduli, respectively (Doi 2009; Etzold et al. 2021).

As shown by both Doi (2009) and Bouklas & Huang (2012), the linear elastic equation (3.5) can be obtained by linearisation of nonlinear constitutive equations similar to those used by Engelsberg & Barros (2013), Bertrand et al. (2016) and others. It is of the same functional form as the constitutive equation of linear elasticity, and thus contains the same terms. However, it represents different microscopic physics that give rise to the same macroscopic behaviour of poroelasticity. Therefore, we refer to $K$, in line with Doi (2009) and Etzold et al. (2021), as the osmotic modulus of the hydrogel with dimensions of pressure, and $G$ as the shear modulus of the hydrogel. Similar to conventional poroelasticity, it is important to realise that $K$ is not the bulk modulus of the solid. Instead, it describes the (elastic) resistance of the porous matrix against volumetric changes (Biot 1941). In a hydrogel, $K$ represents the combination of osmotic (solution) effects and a contribution from the elasticity of the polymer matrix (Doi 2009). It is therefore appropriate to see (3.5) as a simplified version of physically more encompassing nonlinear elastic models. Note that the form of (3.5) implies that, when the system is in its reference state, $\boldsymbol {\sigma }=0$. While this is not always physically appropriate, it is immaterial for the current analysis in this paper; we will return to this point in § 6.

We close the system of equations (3.1a)–(3.5) by invoking Darcy's law for flow within the gel

(3.6)\begin{equation} (1 - \phi_p)(\boldsymbol{u_p} - \boldsymbol{u_f}) = \frac{\kappa}{\mu_f} \boldsymbol{\nabla}p, \end{equation}

where $\kappa = \kappa (\phi _p)$ is the effective hydrogel permeability with characteristic permeability scale $\kappa _0$.

To make progress, we make a further simplifying assumption by defining $\kappa ^{*} = \kappa / \kappa _0$ as an explicit dimensionless function of the polymer volume fraction $\phi _p$. Since the true structure of these materials at the pore scale (microscale) is complex, this can only be an approximation. In the literature on porous media, many different models for the saturation dependence of $\kappa ^{*}$ have been used. In a geophysical context, for mathematical simplicity, Hewitt et al. (2015) took the permeability of a shallow deformable porous layer as $\kappa ^{*} = 1$. Modelling bacterial biofilms, Seminara et al. (2012) set the permeability of a growing biofilm as $\kappa ^{*} = (1 - \phi _p)^2$ while Fortune, Oliveira & Goldstein (2022) constructed the system in such a way to not have to explicitly define $\kappa ^{*}$. When considering a spherical hydrogel, Bertrand et al. (2016) proposed $\kappa ^{*}=(1-\phi _p)/\phi _p^\beta$. For the present study, in the interest of generality, we follow Tokita & Tanaka (1991) and Etzold et al. (2021) by defining

(3.7)\begin{equation} \kappa^{*} = \left(\frac{\phi_{0p}}{\phi_p}\right)^{\beta}, \end{equation}

where $\beta$ is an empirical parameter chosen by matching to experimental data, recalling that $\phi _{0p}$ is defined as the polymer volume fraction at the start of an experiment before the water drop has been added, i.e. at ambient laboratory conditions. For a purely cuboidal network swelled isotropically, previous researchers have proposed $\beta = 2/3$ (Etzold et al. 2021). Al-Amodi & Hill (2022) argued that $\beta =2/3$ could also be obtained for a polyelectrolyte network by considering an affine transformation, a less restrictive condition than that imposed by Etzold et al. (2021). In contrast, experiments have measured values ranging between $\beta =1.5$ (Tokita & Tanaka 1991) and 1.85 (Grattoni et al. 2001).

3.3. Dimensionless shallow-layer scalings

We scale radial and vertical lengths with the initial radius $R_0 = R(t = 0)$ and height ${H_0 = h(r = 0, t = 0)}$ of the blister, respectively (i.e. $\{ r, R \} \sim R_0$ and $\{ z, \, h \} \sim H_0$). Since the characteristic time scale for the system is the poroelastic time scale for pressure-driven fluid flow in the horizontal direction $t_0$, we scale $t \sim t_0$ where

(3.8)\begin{equation} t_0 = \frac{\mu_f R_0^2}{(K+4G/3)\kappa_0}, \end{equation}

and the denominator has been chosen for algebraic convenience in what follows. From (3.1), we have $u_f \sim U_0 = R_0 / t_0$ and $\{ w_f, w_p \} \sim H_0 / t_0$. Since horizontal and vertical elastic stresses in the hydrogel are of the same magnitude, $u_p \sim w_p \sim H_0 / t_0$. Since $\boldsymbol {u_p}$ is defined as the material derivative of $\boldsymbol {\xi }$, we have $\{ \zeta, \xi \} \sim H_0$. By definition, $\kappa \sim \kappa _0$ while (3.6) implies that $\{ p, \tilde {p} \} \sim P_0 = K + 4G/3$. Finally, from volume conservation, $d \sim ( 12 R_0^2 H_0 )^{1/3}$.

3.4. Non-dimensional governing equations

We now exploit a lubrication approximation, noting that we have chosen a suitable starting time $t = 0$ so that the characteristic radial length scale $R_0$ is much greater than the characteristic vertical length scale $H_0$. We non-dimensionalise the governing equations anisotropically using the scalings given above, denoting the dimensionless form of a dimensional variable $f$ by $f^{*}$ and setting

(3.9a,b)\begin{equation} \mathcal{A} = \frac{a}{H_{0}}, \quad \mathcal{V} = \frac{d^3}{12 R_0^2 H_0}. \end{equation}

Note that, since the blister has positive curvature, $H_0 > a$ and thus $0 < \mathcal {A} < 1$. Keeping only leading-order terms in the aspect ratio $\epsilon$ (where $\epsilon = H_0/R_0\ll 1$), the system of (3.1a)–(3.6) becomes

(3.10a)$$\begin{gather} \frac{\partial \phi_p}{\partial t^{*}} + \frac{\partial}{\partial z^{*}}(\phi_p w_p^{*}) = {O}\left(\epsilon \right), \end{gather}$$

(3.10b)$$\begin{gather}-\frac{\partial \phi_p}{\partial t^{*}} + \frac{\partial}{\partial z^{*}}( (1 - \phi_p) w_f^{*}) + \frac{1}{r^{*}}\frac{\partial}{\partial r^{*}}(r^{*}(1 - \phi_p)u_f^{*}) = 0, \end{gather}$$

(3.10c)$$\begin{gather}w_f^{*} - w_p^{*} ={-} \frac{\kappa^{*}}{\epsilon^2(1 - \phi_p)} \frac{\partial P^{*}}{\partial z^{*}}, \end{gather}$$

(3.10d)$$\begin{gather}u_f^{*} ={-} \frac{\kappa^{*}}{(1 - \phi_p)} \frac{\partial P^{*}}{\partial r^{*}} + {O}\left( \epsilon \right), \end{gather}$$

(3.10e)$$\begin{gather}\frac{\partial P^{*}}{\partial r^{*}} = \frac{1}{\epsilon}\left( \frac{G}{K + 4G/3} \right)\frac{\partial^2 \xi^{*}}{\partial {z^{*}}^2} + \left( \frac{K + G/3}{K + 4G/3} \right)\frac{\partial^2 \zeta^{*}}{\partial r^{*} \partial z^{*}} + {O}\left( \epsilon \right), \end{gather}$$

(3.10f)$$\begin{gather}\frac{\partial P^{*}}{\partial z^{*}} = \frac{\partial^2 \zeta^{*}}{\partial {z^{*}}^2} + \epsilon \left( \frac{K + G/3}{K + 4G/3} \right)\frac{1}{r^{*}}\frac{\partial}{\partial r^{*}}\left( r^{*} \frac{\partial \xi^{*}}{\partial z^{*}} \right) + {O}( \epsilon^2 ), \end{gather}$$

where

(3.10g)$$\begin{gather} u_p^{*} = \frac{\partial \xi^{*}}{\partial t^{*}} + w_p^{*} \frac{\partial \xi^{*}}{\partial z^{*}} + {O}\left( \epsilon \right), \end{gather}$$

(3.10h)$$\begin{gather}w_p^{*} = \left(1 - \frac{\partial \zeta^{*}}{\partial z^{*}}\right)^{{-}1} \frac{\partial \zeta^{*}}{\partial t^{*}} + {O}\left( \epsilon \right). \end{gather}$$

3.5. Boundary conditions

The bottom of the hydrogel layer was affixed to a thin rigid plastic sheet which was in turn glued to a microscope slide. We model this by the boundary conditions

(3.11a)\begin{equation} w_f^{*} = 0 \text{ and } \xi^{*} = \zeta^{*} = 0 \quad \Longrightarrow \quad u_p^{*} = w_p^{*} = 0\quad \text{at } z^{*} = 0. \end{equation}

The kinematic boundary conditions for the free surface at $z^{*} = h^{*}$ are

(3.11b)$$\begin{gather} \left. w_p^{*}\right|_{h^{*}} = \frac{\partial h^{*}}{\partial t^{*}} + \left.\epsilon u_p^{*}\right|_{h^{*}} \frac{\partial h^{*}}{\partial r^{*}}, \end{gather}$$

(3.11c)$$\begin{gather}\left. w_f^{*}\right|_{h^{*}} = \frac{\partial h^{*}}{\partial t^{*}} +\left. u_f^{*}\right|_{h^{*}} \frac{\partial h^{*}}{\partial r^{*}}. \end{gather}$$

Vertically integrating (3.10a) and then applying (3.11b) gives the polymer volume conservation condition

(3.11d)\begin{equation} \frac{\partial}{\partial t^{*}}\left( \int^{h^{*}}_{0} \phi_p \,{\rm d}z^{*} \right) = {O}(\epsilon). \end{equation}

The stress-free boundary at $z^{*} = h^{*}$ requires

(3.11e)$$\begin{gather} \left.(\sigma^{*})_{r^{*} z^{*}} \right|_{h^{*}} = 0 \quad \Longrightarrow \quad \left. \frac{\partial \xi^{*}}{\partial z^{*}}\right|_{h^{*}} = 0 + {O}(\epsilon), \end{gather}$$

(3.11f)$$\begin{gather}\left. (P^{*} - (\sigma^{*})_{z^{*}z^{*}})\right|_{{\left(h^{*}\right)}^-} = \left. P^{*}\right|_{{\left(h^{*}\right)}^-} - \left.\frac{\partial \zeta^{*}}{\partial z^{*}}\right|_{ {\left(h^{*}\right)}^-} + {O}\left( \epsilon \right) =\left. P^{*}\right|_{ {\left(h^{*}\right)}^+ } = P^{*}_0, \end{gather}$$

where $\left. P^{*}\right |_{ {(h^{*})}^+} = P_{0}^{*}$ is the external fluid pressure, which we assume to be uniform and constant. Global water conservation requires

(3.11g)\begin{equation} \int^{\infty}_0 r^{*} ({h^{*}}-\mathcal{A}) \, {\rm d}r^{*} = \mathcal{V}. \end{equation}

Note that this condition is independent of the initial polymer volume fraction $\phi _{0p}$. Finally, at large radial distance in the hydrogel, the polymer is undeformed and in mechanical equilibrium such that

(3.11h)\begin{equation} h^{*} \rightarrow \mathcal{A}, \quad \phi_p \rightarrow \phi_{0p}, \quad \xi \rightarrow 0, \quad \zeta \rightarrow 0 \quad \text{as } r^{*} \rightarrow \infty. \end{equation}

4. Perturbation expansion in the aspect ratio

Considering the terms at ${O}(1/\epsilon )$, (3.10e) simplifies to give

(4.1)\begin{equation} \frac{\partial^2 \xi^{*}}{\partial {z^{*}}^2} = 0 + {O}\left( \epsilon \right), \end{equation}

with boundary conditions

(4.2)\begin{equation} \left.\xi^{*} \right|_{0} = \left.\frac{\partial \xi^{*}}{\partial z^{*}} \right|_{h^{*}} = 0 \Longrightarrow \xi^{*} = 0 + {O}\left( \epsilon \right), \end{equation}

(i.e. $\xi^*$ and thus $u_p^*$ are not leading-order terms). Similarly, to find $P^{*}$ we expand out (3.10c) and then apply (3.11f), yielding

(4.3)\begin{equation} \frac{\partial P^{*}}{\partial z^{*}} = 0 + {O}(\epsilon) \quad \Longrightarrow \quad P^{*} = P^{*}_0 +\left. \frac{\partial \zeta^{*}}{\partial z^{*}}\right|_{h^{*}} + {O}(\epsilon). \end{equation}

At ${O}(\epsilon ^{0})$ and using (4.3), we find (3.10 f) simplifies to give

(4.4)\begin{equation} \frac{\partial^2 \zeta^{*}}{\partial {z^{*}}^2} = 0 + {O}\left( \epsilon \right). \end{equation}

Using both (3.10h) and (4.2), the solid-phase kinematic boundary condition given in (3.11b) simplifies to

(4.5)\begin{equation} \frac{\partial h^{*}}{\partial t^{*}} \left( 1 - \left.\frac{\partial \zeta^{*}}{\partial z^{*}}\right|_{h^{*}} \right) = \left.\frac{\partial \zeta^{*}}{\partial t^{*}}\right|_{h^{*}} + {O}(\epsilon). \end{equation}

Hence, integrating (4.4) twice with respect to $z^{*}$ and applying (3.11a) gives

(4.6)\begin{equation} \zeta^{*} = C_{\zeta} z^{*} + {O}(\epsilon), \end{equation}

where the constant of integration $C_{\zeta } = C_{\zeta }(t^{*})$ is independent of $z^{*}$. Substituting this into (4.5), integrating with respect to $t^{*}$ and applying (3.11h) gives

(4.7)\begin{equation} \zeta^{*} = z^{*}\left( 1 - \frac{\mathcal{A}}{h^{*}} \right) + {O}\left( \epsilon \right). \end{equation}

Furthermore, substituting this into (3.10h) and (4.3) yields, respectively,

(4.8)$$\begin{gather} w_p^{*} = \frac{z^{*}}{h^{*}}\frac{\partial h^{*}}{\partial t^{*}} + {O}\left( \epsilon \right), \end{gather}$$

(4.9)$$\begin{gather}P^{*} = P^{*}_0 + 1 - \frac{\mathcal{A}}{h^{*}} + {O}(\epsilon). \end{gather}$$

As a sanity check, note that from (4.7) we recover $\zeta ^{*}|_{z^*=h^{*}} = h^{*} - \mathcal {A}$.

Using (4.8), the polymer volume fraction mass conservation equation (3.10a) simplifies to give

(4.10)\begin{equation} \frac{\partial}{\partial t^{*}}( h^{*}\phi_p ) + z^{*} \frac{\partial h^{*}}{\partial t^{*}} \frac{\partial \phi_p}{\partial z^{*}} = 0 + {O}(\epsilon). \end{equation}

Using the method of separation of variables, we find that this equation admits solutions of the form $\phi _p = f(r^*) {[h^*]}^{-v-1} {z^{*}}^{\nu }$ for any constant $v$ and function $f$. Therefore, a general separable solution is of the form

(4.11)\begin{equation} \phi_p = \phi_{0p}\int \left[ \left( \frac{\mathcal{A}}{h^{*}} \right)^{1 + \nu} {z^{*}}^{\nu} \right] \varPhi \, {\rm d}\nu + {O}(\epsilon), \end{equation}

where $\varPhi = \varPhi (r^{*},\nu )$ is independent of $t^{*}$ and $z^{*}$ and $\nu$ is the eigenvalue spectrum associated with the separable solution. For mathematical simplicity, we will assume that when $t^{*} = 0$, $\phi _{0p}$ is independent of $z^{*}$. Hence, the dominant mode in (4.11) is the $\nu = 0$ mode giving

(4.12)\begin{equation} \phi_{p} = \frac{\mathcal{A} \phi_{0p}}{h^{*}} + {O}\left( \epsilon \right). \end{equation}

Finally, integrating (3.10b) in the $z^{*}$ direction, using both (3.10d) and the boundary conditions given in (3.11b) and (3.11c), yields the continuity equation for the deformation of the sheet $\overline {h^{*}} = h^{*} - \mathcal {A}$, expressible as

(4.13)\begin{equation} \frac{\partial \overline{h^{*}}}{\partial t^{*}} = \frac{1}{r^{*}}\frac{\partial}{\partial r^{*}}\left( r^{*} \int^{\overline{h^{*}} + \mathcal{A}}_0 \kappa^{*} \frac{\partial P^{*}}{\partial r^{*}} \,{\rm d}z^{*} \right) = \frac{\mathcal{A}}{r^{*}}\frac{\partial}{\partial r^{*}}\left( \frac{r^{*} \kappa^{*}}{\overline{h^{*}}+\mathcal{A}}\frac{\partial \overline{h^{*}}}{\partial r^{*}} \right), \end{equation}

with $\mathcal {A}\leq \overline {h^{*}} \leq 1$ for $0\leq r^{*}<\infty$, where we have used (4.12) to write $\kappa ^{*}(\phi )$ in terms of $\overline {h^{*}}$ as $\kappa ^{*} = (1 + \overline {h^{*}}/\mathcal {A})^{\beta }$.

Motivated by the apparent self-similar profile of the experimental data (figure 2), we explore the existence of self-similar solutions to (4.13).

4.1. Self-similar solutions: small deformation limit

In the small deformation limit when $\max (\overline {h^{*}}) \ll 1$, the deformation height profile $\overline {h^{*}} = \overline {h^{*}}_s$ satisfies

(4.14a)\begin{equation} \frac{\partial \overline{h^{*}}_s}{\partial t^{*}} = \nabla^2 \overline{h^{*}}_s, \end{equation}

with volume conservation condition

(4.14b)\begin{equation} \int^{\infty}_0 r^{*} \overline{h^{*}}_s \,{\rm d}r^{*} = \mathcal{V}. \end{equation}

Trying a self-similar solution of the form $f_s(\eta )/{ ( R^{*}_s )^2}$, where $\eta = r^{*}/R^{*}_s$ and $R^{*}_s = R^{*}_s(t^{*})$, (4.14a) simplifies to

(4.15)\begin{equation} -R^{*}_s \frac{\partial R^{*}_s}{\partial t^{*}}\frac{\partial}{\partial \eta}(\eta^2 f_s) = \frac{\partial}{\partial \eta} \left(\eta \frac{\partial f_s}{\partial \eta}\right). \end{equation}

Using separation of variables and the initial condition $R^{*}_s(t = 0) = 1$, the $t^{*}$-dependent terms become

(4.16)\begin{equation} \frac{\partial}{\partial t^{*}}((R^{*}_s)^2 ) = \lambda \Longrightarrow R^{*}_s = (1 + \lambda t^{*})^{1/2}, \end{equation}

where $\lambda$ is a constant. Similarly, using regularity of $\overline {h^{*}}_s$ at $\eta = 0$ and the initial condition $\overline {h^{*}}_s(r^{*} = 0, t^{*} = 0) = 1-\mathcal {A}$, the $\eta$-dependent terms become

(4.17)\begin{align} -\frac{\lambda}{2}\frac{\partial}{\partial \eta}(\eta^2 f_s) &= \frac{\partial}{\partial \eta}\left( \eta \frac{\partial f_s}{\partial \eta} \right) \quad \Longrightarrow \quad \frac{\partial f_s}{\partial \eta} ={-}\frac{\lambda \eta}{2}f_s \nonumber\\ \quad \Longrightarrow\quad f_s &= (1-\mathcal{A}) \exp{\left( -\frac{\lambda \eta^2}{4} \right)}. \end{align}

Applying the volume conservation condition in (4.14b) sets $\lambda = 2(1-\mathcal {A})/\mathcal {V}$, and thus we recover the self-similar solution

(4.18)\begin{equation} h^{*}_{s} - \mathcal{A} = \frac{1}{\left(R_s^{*}\right)^2} \exp{\left(-\frac{(1-\mathcal{A})}{2\mathcal{V}}\left(\frac{r^{*}}{R_s^{*}}\right)^2\right)}, \end{equation}

where

(4.19)\begin{equation} R^{*}_s = \left( 1 + \frac{2(1-\mathcal{A}) t^{*}}{\mathcal{V}} \right)^{1/2}. \end{equation}

Note that this is independent of the permeability exponent $\beta$.

4.2. Self-similar solutions: large deformation limit

We also consider the case of large deformation, for which we pose

(4.20)\begin{equation} \overline{h^{*}} \gg \mathcal{A}. \end{equation}

Physically, it is clear that this cannot be true everywhere, since, as expressed in (3.11h), the height profile must approach the undeformed value $\mathcal {A}$ for large $r^{*}$. However, if $\max (\overline {h^{*}}) \gg \mathcal {A}$, it is reasonable to assume that (4.20) is valid for most of the domain. Combining (4.20) with (4.13) yields

(4.21a)\begin{equation} \frac{\partial \overline{h_l^{*}}}{\partial t^{*}} = {\mathcal{A}}^{1 - \beta}\frac{1}{r^{*}}\frac{\partial}{\partial r^{*}}\left( r^{*} {(\overline{h_l^{*}})}^{\beta - 1} \frac{\partial \overline{h_l^{*}}}{\partial r^{*}} \right), \end{equation}

with volume conservation condition

(4.21b)\begin{equation} \int^{\infty}_0 r^{*} \overline{h_l^{*}} \,{\rm d}r^{*} = \mathcal{V},\end{equation}

whose solutions are expected to describe the behaviour of the solutions to (4.13) well enough for strongly swollen systems.

Trying a self-similar solution of the form $f_l(\eta )/{(R^{*}_l)^2}$, where $\eta = r^{*}/R^{*}_l$ and $R^{*}_l = R^{*}_l(t^{*})$, (4.21a) simplifies to

(4.22)\begin{equation} -(R^{*}_l)^{2\beta - 1} \frac{\partial R^{*}_l}{\partial t^{*}}\frac{\partial}{\partial \eta}(\eta^2 f_l) = {\mathcal{A}}^{1 - \beta}\frac{\partial}{\partial \eta} \left(\eta f_l^{\beta - 1}\frac{\partial f_l}{\partial \eta}\right). \end{equation}

Using separation of variables and the initial condition $R^{*}_l(t = 0) = 1$, the $t^{*}$-dependent terms become

(4.23)\begin{equation} \frac{\partial}{\partial t^{*}}((R^{*}_l)^{2\beta}) = \lambda \Longrightarrow R^{*}_l = (1 + \lambda t^{*})^{1/2\beta}, \end{equation}

where $\lambda$ is a constant. Similarly, using regularity of $\overline {h_l^{*}}$ at $\eta = 0$ and the initial condition $\overline {h_l^{*}}(r^{*} = 0, t^{*} = 0) = 1-\mathcal {A}$, the $\eta$-dependent terms become

(4.24) \begin{align} -\frac{\lambda {\mathcal{A}}^{\beta - 1}}{2\beta}\frac{\partial}{\partial \eta}( \eta^2 f_l ) &= \frac{\partial}{\partial \eta}\left( \eta f_l^{\beta - 1} \frac{\partial f_l}{\partial \eta} \right) \quad \Longrightarrow \quad \frac{\partial}{\partial \eta}(\,f_l^{\beta - 1} ) ={-}\frac{\lambda \eta (\beta - 1) \mathcal{A}^{\beta - 1}}{2\beta}, \nonumber\\ \Longrightarrow \quad f_l &= (1 - \mathcal{A})( 1 - C \eta^2 ) ^{1/(\beta - 1)}, \end{align}

where $C$ satisfies

(4.25)\begin{equation} C = \frac{\lambda (\beta - 1) \mathcal{A}^{\beta - 1}}{4\beta{(1 - \mathcal{A})^{\beta - 1}}}. \end{equation}

The value of $\lambda$ is now determined from conservation of volume (4.21b). If $\beta >1$, solutions are only physically meaningful between $0< r^{*}< R_{0l}^{*}$, where $R_{0l}^{*} = R_l^{*}/\sqrt {C}$ denotes the first root of (4.24), since (4.24) is either complex (for non-integer $1/(\beta -1)$) or diverges (for integer $1/(\beta -1)$) for $r^{*}>R_{0l}^{*}$.

Thus, the volume conservation condition (4.21b) becomes

(4.26)\begin{equation} \int^{R^{*}_{0l}}_0 r^{*} \overline{h_l^{*}} \,{\rm d}r^{*} = \mathcal{V}, \end{equation}

and we find

(4.27)\begin{equation} \lambda_{\beta>1}= \frac{2(1 - \mathcal{A})^{{\beta}}}{\mathcal{V}\mathcal{A}^{\beta - 1}}. \end{equation}

If $\beta <1$ and $C<0$, (4.26) is integrable for $0< r^{*}<\infty$, whilst the integral diverges if $C>0$. Thus we integrate from $0< r^{*}<\infty$ and find that $\lambda$ similarly satisfies

(4.28)\begin{equation} \lambda_{\beta<1}=\frac{2(1-\mathcal{A})^{\beta}}{\mathcal{V} \mathcal{A}^{\beta-1}}.\end{equation}

Note that the expressions for $\lambda _{\beta >1}$ and $\lambda _{\beta <1}$ are functionally identical, despite the different algebra required in the two cases. Thus, regardless of the value of $\beta$, the large deformation solution is

(4.29)\begin{equation} h^{*}_{l} - \mathcal{A} = \frac{(1 - \mathcal{A})}{\left(R_l^{*}\right)^2} \left( 1 - \frac{{(1 - \mathcal{A})}(\beta - 1)}{2\beta \mathcal{V}}\left( \frac{r^{*}}{R^{*}_l} \right)^2 \right)^{1/(\beta - 1)}, \end{equation}

with $0\leq r^{*}\leq \infty$ if $\beta <1$ and $0\leq r^{*}\leq R^{*}_l$ if $\beta >1$ and

(4.30)\begin{equation} R^{*}_l = \left( 1 + \frac{2 t^{*}(1 - \mathcal{A})^{{\beta}} {\mathcal{A}}^{1-\beta}}{\mathcal{V}} \right)^{1/(2\beta)}. \end{equation}

Note that, as for $\lambda$, $\overline {h^*}$ and $R^*_l$ are also functionally identical in the two cases ($\beta > 1$ and $0 < \beta < 1$) but have different support. These large deformation solutions are, unlike the small deformation case (4.18), functions of $\beta$. In the limit when $\beta = 1$, $\overline {h^{*}_l}$ collapses to the small deformation similarity solution $\overline {h^{*}_s}$ from § 4.1.

The fact that the assumption (4.20) used to obtain (4.21a) becomes unphysical for large $r^{*}$ manifests in the fact that for $\beta >1$ no real solutions exist for $r^{*}>R_l^{*}$ and that the asymptotic decay of (4.29) towards zero as $r^{*}\rightarrow \infty$ contradicts (4.20). In both cases the dynamics in the outer regions of the swollen region cannot be captured correctly. As stated above, as long as $\max (\overline {h^{*}}) \gg \mathcal {A}$ (4.20) is valid for most of the blister we expect (4.29) to describe the dynamics reasonably well. This is numerically confirmed below in § 4.3.

4.3. Numerical exploration

The similarity solutions are only valid in asymptotic limits. Furthermore, any solution in the large deformation regime will eventually move as time progresses towards the small deformation regime. Hence, we explore cases with large, intermediate and small $\mathcal {A}$. For general $\mathcal {A}$, the nonlinear diffusion equation given in (4.13) does not admit an analytic solution. Therefore, it is solved numerically using the finite-element software package FEniCS (Logg et al. 2012; Alnæs et al. 2015) on a domain ${\mathcal {Q}:0< r^{*}<100}$, bounded by no-flux boundary conditions. The time derivative is discretised using the backward Euler scheme

(4.31)\begin{equation} \frac{\partial h^{*}}{\partial t^{*}} \approx \frac{(h^{*} - h^{*}_{-1})}{\delta t^{*}}, \end{equation}

where $h^{*}_{{-1}}$ is the value of $h^{*}$ at the previous time step. This leads to the weak form of (4.13)

(4.32)\begin{equation} \int r^{*} Q h^{*}_{{-1}} \,{\rm d}{r^{*}} = \mathcal{A}^{1 - \beta} \delta t^{*} \int r^{*} (h^{*})^{\beta-1} \frac{\partial h^{*}}{\partial r^{*}} \frac{\partial Q}{\partial r^{*}}\, {\rm d}{r^{*}} + \int r^{*} h^{*} Q \,{\rm d}{r^{*}}, \end{equation}

where $Q$ is a test function.

The solutions were calculated with timestep $\delta t^{*} = 10^{-3}$ and constant grid spacing ${\delta r^{*} = 10^{-2}}$.

For a given blister height profile $h^{*} = h^{*}(r^{*},t^{*})$, to be consistent with the experiments, the effective radius $R^{*} = R^{*}(\mathcal {A},t^{*})$ is defined as the point where the deformation height of the sheet is half the maximum deformation height of the sheet, e.g.

(4.33) \begin{equation} R^{*} = \{ r^{*} : h^{*}(r^{*})-\mathcal{A} = \tfrac{1}{2}(\underset{\widetilde{r^{*}}}{\max}(h^{*}(\widetilde{r^{*}}))-\mathcal{A} ) \}.\end{equation}

Note that this is the non-dimensional version of (2.2). Figure 5 illustrates the time evolution of numerical solutions of (4.13) for $R^{*}$ for cases with $\beta =3/2>1$ (shown in figure 5a) and $\beta =2/3<1$ (shown in figure 5b), both showing a range of values of $\mathcal {A}$ (darker blue curves denote larger $\mathcal {A}$). Dotted lines with circles denote numerical solutions of (4.13) evolved from an initial condition of the form given in (4.29). Solid lines denote the corresponding large deformation similarity solutions (given in (4.30)). If the initial deformation is large (small $\mathcal {A}$) and $t^{*}$ is small, the large deformation similarity solution agrees well with the numerics in both cases, where $R^{*} \sim (t^{*})^{1/3}$ for $\beta =3/2$ and $R^{*} \sim (t^{*})^{3/4}$ for $\beta =2/3$. Since the maximum height of the blister $\max (h^{*})$ decreases as time passes, the system eventually reaches the small deformation region with corresponding similarity solution for $R^{*}$ given in (4.18), with $R^{*} \sim ( t^{*})^{1/2}$. Consequently, the similarity solutions (4.29) and (4.30) are useful to describe the behaviour of the system if $\max (\overline {h^{*}})\gg \mathcal {A}$.

Figure 5. Swelling dynamics of thin hydrogel sheets according to the proposed poroelastic model for(a) $\beta =3/2>1$ and (b) $\beta = 2/3<1$; for both figures $\mathcal {V}=1-\mathcal {A}$ and $\mathcal {A} \in [0.001,0.01,0.1,0.5,0.9]$ (light blue to dark blue), namely darker shades of blue denote larger $\mathcal {A}$. Dotted lines with circles plot the scaled blister radius $R^{*} = R(t)/R(0)$ of numerical solutions of (4.13) as a function of time, starting from a self-similar initial height profile of the form given in (4.29). For comparison, solid lines show the radius predicted by the self-similar solution in the large deformation limit (see (4.30)).

5. Comparison between theory and experiment

5.1. Fitting procedure

The experiments reported here have $0.43 \leq \mathcal {A} \leq 0.68$ and so lie in the intermediate deformation regime. This, together with the need to determine virtual origins in both time and space for the experimental data, precludes a direct quantitative test of the similarity solutions for the small and large deformation regimes (which is expected to describe the late-time behaviour of the blister after an initial adjustment period) against our experimental data. Instead, we solve (4.13) numerically, using an experimentally observed initial condition. The initial conditions were taken to be the experimental height profile at $t = t_s$. We chose $t_s$ separately for each experiment to ensure that the model assumptions were justified, namely that the uneven surface caused by the buckling/crumpling instability had decayed sufficiently so that the assumption that the profile was axisymmetric was valid. We also took this as an indication that the vertical gradients in polymer volume fraction had decreased (see (4.12)). This also ensured that all water had been absorbed and a no-flux boundary condition would be present at the surface of the hydrogel.

Note that this means that, in the subsequent figures, $t^{*}=0$ corresponds to the experimental time $t_s$. Table 1 in Appendix D gives the associated experimental times $t_s$, the corresponding frame number and the derived values for each dataset for $\{R_0, H_0, \mathcal {A}, t_0 \}$.

Table 1. Experimental $t_s$ from which the initial condition for the numerical solution was obtained with corresponding values $\{ H_0, R_0, \mathcal {A}, t_0 \}$. The second column refers to the frame number in the corresponding attached raw dataset.

Since numerical solutions predict the evolution of the blister as a function of dimensionless time $t^{*}$ with a single fitted parameter $\beta$, the poroelastic time scale $t_0 = \mu _f R_0^2 / \kappa _0 (K + 4G/3)$ for the propagation of the swollen layer through the hydrogel during the initial absorption of the droplet has to be determined through fitting. Note that $t_0$ can be split into a part that is experiment dependent and a part that is experiment independent, namely $t_0 = \varOmega \, R_0^2$, where the experiment-dependent part $R_0^2$ is found through image analysis while the experiment-independent part $\varOmega = \mu _f / \kappa _0 (K + 4G/3)$ is the same unknown constant fitted across all the experimental data. Furthermore, note that we do not need to fit the effective elastic moduli $K$ and $G$ individually, decreasing the number of free parameters we need to fit from the experimental data and hence increasing the value of the fit. In particular, for a given value for $\beta$, we determined $\varOmega$ through a simultaneous fit to all the experimental datasets, minimising the normalised mean square error for the maximum observed height, as expressed by the objective function

(5.1)\begin{equation} \sum_i \left\{\underset{j}{\mathrm{avg}}\left( \left.\left[{ \max(h_{e}}(t_j))- \max{\left(h_\beta\left(\frac{t_j}{\varOmega R^2_0}\right)\right)} \right] \right/ \max(h_{e}(t_j)) \right)^2\right\}_i. \end{equation}

Here, $i$ corresponds to different experimental datasets, namely different initial droplet volumes. The index $j$ corresponds to different points in time $t_j$ with corresponding hydrogel height profile $h_e(t_j)$ for a given experiment. Finally, $h_{\beta _1}(t^{*}_1)$ denotes a height profile generated numerically at time $t^{*} = t^{*}_1$ for $\beta = \beta _1$. Note that we are minimising across all experimental datasets at once.

Both this minimisation and the subsequent linear interpolation to calculate $h^{*}_\beta$ at the experimental time steps were performed using standard NumPy and SciPy methods in Python (Harris et al. 2020; Virtanen et al. 2020). This fit was constructed with the aim to replicate the maximum (centre) height of the blister.

Figure 6 explores graphically how varying $\beta$ affects the temporal evolution of the numerical solutions for $\max {(h^{*})}$ and $R^{*}$. As a first sweep, we calculated the evolution of each profile for $2/3\leq \beta <3$ in increments of 1/3. From this, we could bound $\beta$ to the range $\beta \in [4/3,7/3]$ until we investigated $\beta =[25/12,26/12,27/12]$. The final values adopted were those obtained for the $\beta$ with the lowest value of the objective function in (5.1), namely $\beta = 2.25\pm 0.083$ and $\varOmega = 6.72 \times 10^{9}\,{\rm s} \,{\rm m}^{-2}$, which corresponds to poroelastic time scales $t_0 = \varOmega R_0^2 = 4.14-36.79\,{\rm h}$, as shown in table 1. This large variation in time scales arises since the initial radii vary by approximately a factor of three between datasets.

Figure 6. The experimental temporal evolution of (a) the scaled maximum blister height $\max (h^{*})$ and (b) the characteristic blister radius $R^{*}$ formed by a $50\,{\mathrm {\mu }}{\rm l}$ droplet (dots) superimposed onto a set of blue lines computed using the theoretical model (4.13), using the initial experimental height profile as the initial condition, where $\beta \in \{ 2/3,1,4/3,11/6,25/12,9/4,5/2,3 \}$. Curves in a darker shade of blue correspond to larger values of $\beta$. The dimensionless time scale $t^{*}$ was matched to real time $t$ using the relation $\varOmega = 6.72\times 10^{9}\,{\rm s}\,{\rm m}^{-2}$ obtained as best fit for $\beta =2.25$ (thicker blue line). Error bars are generated through an uncertainty analysis that is described in § 2.4 and Appendix A.

Whilst $\beta =2.25$ is considerably larger than the suggestion of $2/3$ made by Etzold et al. (2021), it is more comparable to many previous experimental studies which found $1.5\leq \beta \leq 1.85$ (Tokita & Tanaka 1991; Grattoni et al. 2001; Bertrand et al. 2016). This may be because the fit is not particularly sensitive to $\beta$ with the temporal evolution of $R^{*}$ and $h^{*}$ qualitatively maintaining the same key features. Quantitatively, when $\beta$ increases, $h^{*}$ decays faster while $R^{*}$ increases faster. For example, if we had chosen $\beta =1.83$, (approximately the highest value reported in the literature, Grattoni et al. 2001), $\varOmega$ would have decreased slightly to $5.82\times 10^{9}\,{\rm s}\,{\rm m}^{-2}$ (a change of around 15 %). This small change does not influence the characteristic trends that are observed as the system evolves over time. This insensitivity of the fit might be due to only modest changes of the observables over the observation period (factor 2). While we have not tried expressions for $\kappa ^{*}$ in terms of $\phi _p$ other than (3.7), at small $\phi _p$ the leading-order term in the Taylor expansion will be of the polynomial form we have chosen in our model (3.7).

5.2. Profile comparison

Figures 7 and 8 compare experimental data (dots) with numerical predictions (lines) generated using the model (4.13) and the fitted parameters ($\beta = 2.25$ and $\varOmega = 6.72 \times 10^{9}\,{\rm s}\,{\rm m}^{-2}$), taking initial conditions directly from the experimental data that they are fitted against (darker blue curves indicate smaller initial droplets). Figure 7 explores how the maximum scaled blister height $\max (h^{*})$ evolves temporally with $t^*$ for different droplet volumes. Agreement within the estimated uncertainty is found for all datasets for the full range of experimental data. The inset compares experimental height profiles of $h^{*}(r^{*})$ at a late time towards the end of each dataset. The red stars on the main plot indicate the corresponding points in time. Within the estimated uncertainty, the model and experiments agree well.

Figure 7. The temporal evolution of the maximum scaled height $\max (h^{*})$. Inset: comparison of experimental profiles (dots), at the times indicated by the red stars, with numerical profiles (solid lines). Error bars represent uncertainties. Curves in lighter shades of blue correspond to larger initial droplet volumes. Numerical solutions are generated using initial conditions taken directly from the experimental data that they are fitted against. Error bars are generated through an uncertainty analysis that is described in § 2.4 and Appendix A.

Figure 8. The temporal evolution of the scaled radius $R^{*}$ of the swelling region for a range of different initial droplet volumes, plotting experimental data with dots and numerical solutions with lines. Error bars represent uncertainties in the experimental data. Curves in lighter shades of blue correspond to larger initial droplet volumes, with the same volume/colour combinations as in figure 7. Numerical solutions are generated using initial conditions taken directly from the experimental data that they are fitted against once the crumpling instability has decayed. Red dots indicate the strongest deviations between numerics and experiment at late times. Error bars are generated through an uncertainty analysis that is described in § 2.4 and Appendix A.

Figure 8 compares the half-height radius $R^{*}$ of the blister obtained from experimental and numerical data. We find agreement, generally comfortably within the uncertainty range of the experiment for all observation times, noting that the experimental uncertainty associated with $R^{*}$ is much greater than the corresponding uncertainty associated with $\max {( h^{*} )}$ (see Appendix A for further details). However, both the $5\,{\mathrm {\mu }}{\rm l}$ drop (the darkest blue curve) and the $30\,{\mathrm {\mu }}{\rm l}$ drop (the third curve from the bottom) seem to indicate some systematic deviations at later times. Although these deviations could arise from a breakdown of our assumptions, given the agreement shown in figure 7, it is more likely that these arise from limitations in the experimental set-up. For example, the accurate determination of $R^{*}$ becomes more difficult with decreasing slope of the blister.

Through log–log plots (not shown here), we have explored how well the experimental data shown in figures 7 and 8 approach the small and large deformation regimes (4.18) and (4.29), respectively, which have distinct power law evolution at large time. It was found that little quantitative information could be gained since the experimental data mostly lie in an intermediate deformation region, and the duration of the experiments was not long enough to clearly observe any asymptotic power law dependence. Hence, the blister evolution in our experiments was still influenced by virtual origins in time, arising from the initial conditions.

Expanding on the inset of figure 7, figure 9 compares the spatial dependence of the experimental profiles (dots) with the numerical predictions (lines) for the $50\,{\mathrm {\mu }}{\rm l}$ dataset. The lightest blue curve corresponds to a time of 2 min after the first frame. Each subsequent darker blue line corresponds to twice the previous time. There is generally agreement between experiment and theoretical predictions. However, as time increases, the numerical solution appears to over-predict the position of the leading edge of the swollen region, although still remaining within experimental error.

Figure 9. Comparison between experimentally observed swelling and numerical solutions of the poroelastic model, plotting the spatial dependence of the scaled height $h^{*}$ for a blister formed by a $50\,{\mathrm {\mu }}{\rm l}$ droplet at different times during the experiment. Here, experimental data are denoted by circles while numerical solutions are denoted by lines. The lightest blue curve corresponds to $t - t_s = 2\,{\rm min}$ with each subsequent darker blue curve twice the previous time. Curves in darker shades of blue denote later times. Numerical solutions are generated using initial conditions taken directly from the experimental data that they are fitted against. Error bars are generated through an uncertainty analysis that is described in § 2.4 and Appendix A.

Discussed in more detail below in Appendix A, the imaging system cannot detect very small swelling (estimated to be less than $45\,{\mathrm {\mu }}{\rm m}$) due both to surface roughness and imperfect alignment between the hydrogel sheet and the camera. At later stages of the experiment, this could hide the tip of the swollen region from observation. This would also affect the computation of $R^{*}$ due to the shallow gradients but have less impact on $\max (h^{*})$, hence providing an explanation for why we see better agreement in figure 7 than in figure 8.

In our model, we have assumed that the elastic deformation gradients in the $z$ direction are negligible. This assumption is reasonable since the characteristic poroelastic time scale $t_0$ is much greater than the time scale on which these gradients decay. Indeed, using experimental observations, we can relate these gradients to the buckling/crumpling instability, which vanishes within approximately 30 min. Since we have found ${t_0 = 4.14 -36.79\,{\rm h}}$, this confirms empirically our assumption that the elastic deformation gradients in the $z$ direction are only important for the early-time swelling dynamics, not captured by our model.

6. Discussion and conclusions

In this paper, we have reported experimental results for the absorption and subsequent transport of a water droplet into a thin hydrogel sheet. Our experimental results are compared with predictions from a two-parameter model, which show good agreement. Our results fill multiple important gaps in the literature on hydrogel swelling. We present results on a previously unstudied canonical geometry, applying both experimental and theoretical methods to study it.

Our experimental data constitute a valuable test for the hydrogel models introduced in § 1 for simpler geometries (e.g. Hong et al. 2008; Chester & Anand 2010; Engelsberg & Barros 2013; Bertrand et al. 2016). We introduce an experimental method to set a near-stress-free (von-Neumann) boundary condition on the hydrogel surface for the fluid flow in the hydrogel using a water-immiscible oil. We then see a temporary shift towards a no-slip (Dirichlet) boundary condition after placing the drop, returning to a stress-free condition once there is no longer liquid water at the surface of the hydrogel. Hence, this enables the study of swelling problems with spatially and temporally varying boundary conditions.

From a modelling perspective, we chose a linear poroelastic framework with nonlinear kinematics (a common choice in the study of poroelasticity, e.g. MacMinn, Dufresne & Wettlaufer 2015) to model our observations. Furthermore, we showed how a thin-layer (lubrication) approximation can reduce the full axisymmetric three-dimensional swelling model to a single one-dimensional nonlinear diffusion equation for the evolution of the non-dimensional height of the layer $h^{*}$. This approach extends the geometries potentially amenable to analytical treatment beyond strictly one-dimensional geometries (Yoon et al. 2010; Engelsberg & Barros 2013; Bertrand et al. 2016). We hope this will inspire future work to be more ambitious than one-dimensional studies.

Our work demonstrated that a simple poroelastic modelling framework can, after fitting to the experimental data for the maximum blister height, predict the shape of the swollen blister over the observation time. This model therefore captures the key dynamics of the hydrogel system, as already demonstrated by Yoon et al. (2010) for a one-dimensional system.

Our model admits self-similar solutions of the one-dimensional time-dependent nonlinear diffusion transport equation. In the large deformation limit, the blister extent was captured by a $t^{1/(2\beta )}$ power law, with an exponent-dependent relationship between the permeability and polymer volume fraction via the parameter $\beta$. In the small deformation limit, the blister profile recovered the canonical diffusive Gaussian profile with the blister extent obeying a $t^{1/2}$ power law.

Our numerical exploration of the one-dimensional nonlinear diffusion equation showed that the blister evolution was broadly captured between the large and small deformation regimes. However, our experimental data showed the impact of virtual origins in time, which meant that self-similar large and small deformation regimes could only be reached in the limit of large time or large domain size. Hence, domain size constraints in our experimental set-up prevented us from capturing fully these two regimes. Therefore, we anticipate that for most finite-size domains, the blister lies in an intermediate deformation regime that is captured by our one-dimensional nonlinear diffusion equation, and whose solution does not require as sophisticated and computationally expensive numerical approaches as solving the full original transport problem.

Recent models for hydrogel swelling (e.g. Hong et al. 2008; Doi 2009; Chester & Anand 2010; Bertrand et al. 2016) use fully nonlinear kinematics. Such models require additional information to derive the constitutive equation for the stress tensor. For example, the Flory–Rehner and Flory–Huggins models (Treloar 1975; Cai & Suo 2012) are based on statistical thermodynamics. As demonstrated by Doi (2009) and Bouklas & Huang (2012), the linear poroelastic equation for the stress tensor that we have used can be obtained by linearisation of the above models. Therefore, the linearised model describes the same physics, albeit potentially with reduced accuracy over a larger range of strains.

Our choice of the linearised constitutive equation for the stress tensor was motivated by simplicity. The nonlinear models have their own caveats (see for example Chester & Anand 2010 on the limitations of the Gaussian-chain model for polymer chains used in many nonlinear hydrogel models). Furthermore, they introduce additional parameters. In our experimental system, independent determination of all these parameters is difficult. This is because we rely on a commercially available hydrogel that we cannot manufacture to the shapes necessary to conduct other measurements (e.g. with the device described in Li et al. 2012). Hence, these would become additional fit parameters which we felt prudent to avoid. More importantly, the linearised equations are mathematically simpler. This is desirable to make analytical progress and facilitate insight into the key macroscopic dynamics. Our choice is justified a posteriori by the agreement between experiment and model.

We recognise that our simplification has a potential price. Formally, we use a linear elastic constitutive equation for the stress tensor beyond the strict validity range of the linearisation of the kinematics, since our deformation gradients are of order one (Spencer 2004; Doi 2009). It is therefore possible that limitations of our model have been absorbed into fitted parameters, in particular the permeability exponent $\beta$. This could explain the difference we find here ($\beta \approx 2.25$), compared with the value of $\beta =2/3$ proposed by Etzold et al. (2021) and Al-Amodi & Hill (2022) based on theoretical geometric arguments. This difference warrants further investigation although it should be noted that our value is similar to other hydrogel experimental studies (e.g. 1.5 from Tokita & Tanaka 1991 and 1.85 from Grattoni et al. 2001).

A compelling idea, which could explain the range of $\beta$ found in our experiments and reported in the literature emerged during the review process for this article through communication with one of the anonymous reviewers. In our experiment, we consider a hydrogel that in its reference state has a relatively high polymer volume fraction. Therefore, in such a state, the polymer chains exist as collapsed coils. These will fill a considerable amount of the pore space with fluid flow being governed by the much smaller pores within the coil. As the polymer swells the chains stretch and the coils unravel. Fluid flow will become dominated by the larger pores between the polymer chains, whose subsequent size change (about the now swollen state) could be described by an affine transformation. Thus, $\beta =2/3$ might be a limiting case for lower polymer volume fractions, whilst drier gels see a larger change in permeability with swelling, leading to larger $\beta$. While its exploration is beyond the scope of this paper, we hope this will inspire future computer modelling studies.

Despite the simplifications discussed above, the agreement between the observed blister shapes and the predictions by the model shows that, at least from a mathematical point of view, our model includes all the key features needed to reproduce the data.

Another limitation of the current framework arises from the construction of the linearisation that leads to (3.2) and (3.5). We chose the unswollen hydrogel as the reference state, namely implying from (3.5) and (3.2) that the pressure $p$ is zero in this state. Since our model only has von-Neumann boundary conditions for $p$, its absolute value is not important. This is not appropriate when modelling situations with a Dirichlet boundary condition on the surface such as during the early stages of swelling, where this boundary condition for $p$ would model the water hydrogel boundary (e.g. Bertrand et al. 2016). In this case, the absolute value of $p$ becomes important since, for consistency, the pressure $p$ for pure water must be equal to ambient pressure. This is not captured by the forms of (3.2) and (3.5) used in this paper, but can likely be included by careful linearisation around an appropriate reference state (e.g. the fully swollen state such as in Etzold et al. 2021).

One of the key questions that arise from the above discussion is the need to understand under which circumstances the additional features of a nonlinear framework would be required to describe experimental data.

For example within the presented problem, the surface instability visible during the early stages of swelling indicates the presence of thin, highly swollen layers which might no longer be described by a linearised framework. A version of our model employing a nonlinear constitutive equation for the stress tensor would certainly contribute to answer these questions, with our dataset providing an experimental base to compare such a nonlinear model against. This would also pave the way towards studying the early stages of swelling. One should note, however, that the thin-layer approximation we have used is not valid for the finite amplitude manifestation of the instability nor when the length scale of the instability is comparable to or smaller than the thickness of the hydrogel.

An experimental challenge when exploring the above and related hydrogel swelling problems more generally is the material itself. Whilst hydrogel synthesis is in principle not too difficult, the reproducible synthesis of geometries similar to sheets adds considerable up-front technical challenges. Our choice to rely on commercially available material avoided these but introduces other limitations such as the gauze layer (§ 2.1). We therefore believe that the study of high-swelling hydrogels would benefit from a collection of ‘standard’ recipes leading to well-characterised materials accessible to a broad range of researchers. This would enable further studies, with a reduced need to determine material parameters by fitting.

In summary, our approach establishes experimental and theoretical connections between the polymer swelling literature and a large body of work considering spreading in thin porous layers (e.g. Nordbotten & Celia 2006; Phillips 2009; Rutqvist 2012; Hewitt et al. 2015). This will accelerate further progress for complex swelling problems in slender geometries, where the polymer interacts with its environment in a complex fashion. Examples are the absorption of droplets on polymer materials (Phadnis et al. 2018), mass transfer between a shear flow and swelling polymer surfaces (see Landel et al. 2016; Delavoipière et al. 2018, for non-swelling examples) and biofilm growth on tissue (Fortune et al. 2022).