# Droplet absorption and spreading into thin layers of polymer hydrogels

Merlin A Etzold, George T Fortune, Julien R Landel, Stuart B Dalziel

<figure><img src="image_1.png" /></figure>

► To cite this version:

Merlin A Etzold, George T Fortune, Julien R Landel, Stuart B Dalziel. Droplet absorption and spreading into thin layers of polymer hydrogels. Journal of Fluid Mechanics, 2023, 974, pp.A7. 10.1017/jfm.2023.655. hal-04651763

HAL Id: hal-04651763

https://hal.science/hal-04651763v1

Submitted on 17 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Droplet absorption and spreading into thin layers of polymer hydrogels

Merlin A. Etzold<sup>1,2,†</sup>, George T. Fortune<sup>2</sup>, Julien R. Landel<sup>3</sup> and Stuart B. Dalziel<sup>2</sup>

(Received 1 December 2022; revised 17 June 2023; accepted 10 July 2023)

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.

**Key words:** porous media, lubrication theory, coupled diffusion and flow

† Email address for correspondence: metzold@dstl.gov.uk

<footnote n="1"><sup>1</sup>The Defence Science and Technology Laboratory, Porton Down, Salisbury, Wiltshire SP4 0JQ, UK</footnote>

<footnote n="2"><sup>2</sup>Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK</footnote>

<footnote n="3"><sup>3</sup>Department of Mathematics, University of Manchester, Alan Turing Building, Oxford Road, Manchester M13 9PL, UK</footnote>
# 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).
Droplet absorption and spreading into thin hydrogel layers

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 µ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 × 8 cm² 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°C), their initial thickness *a* was determined with a micrometer as $a = 1.13 \text{ mm} \pm 0.01 \text{ mm}$. One side of a sheet was affixed to a thin plastic sheet (impermeable to water) of approximate thickness 0.04 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 × 25 mm². The coupons were then glued via the thin impermeable plastic sheet onto 76 × 26 mm² 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.

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.
Droplet absorption and spreading into thin hydrogel layers

<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

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,
<figure><img src="image_5.png" /><figcaption>Figure 2. Sequence of photographs showing the absorption of a 100 µ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.</figcaption></figure>

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× 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

$$h_a = h(x, t) - a, \qquad (2.1)$$

recalling that $h$ is the height of the hydrogel sheet, $a$ is the initial undeformed hydrogel sheet thickness and $0 \le 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 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.
## 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,

$$R = \{r : h_a(r + r_0) = \frac{1}{2}(\max_x(h_a(x, t)))\}. \qquad (2.2)$$

Due to the symmetry of the experimental profile about $r = r_0$, this leads to two values, $R_+$ and $R_-$, with $R_- \le 0 \le 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

$$V = \pi \int_{0}^{x_{\max}} (x - r_{0})h_{a}(x) dx, \qquad (2.3)$$

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

$$\delta R \approx \left. \frac{\delta h_a}{\left( \frac{\partial h_a}{\partial x} \right)} \right|_{x=r_0+R} . \qquad (2.4)$$

Here, $\delta h_a = (+22.7, -12.5)$ µ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

$$\delta V \approx \pi \int_{x_{mn}}^{x_{mx}} (x - r_0) \delta h_a(x) dx, \qquad (2.5)$$

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

$$\{x_{mn} \le x \le x_{mx} : h_a(x - r_0) > \delta z\}. \qquad (2.6)$$

## 2.5. Observations

A montage of typical experimental images for a 100 µ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.
M.A. Etzold, G.T. Fortune, J.R. Landel and S.B. Dalziel

<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

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 µ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).
Droplet absorption and spreading into thin hydrogel layers

<figure><img src="image_7.png" /><figcaption>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 $\mathbf{u}_p = (u_p, w_p)$ and $\mathbf{u}_f = (u_f, w_f)$ respectively.</figcaption></figure>

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 \le r < \infty$ and $0 \le t < \infty$. For an arbitrary function $f = f(z)$, we define the vertically averaging operator $\langle \cdot \rangle$ where $\langle f \rangle = (\int_0^h f dz)/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 $\{\mathbf{u}_p = (u_p, w_p), \boldsymbol{\sigma}_p\}$ and $\{\mathbf{u}_f = (u_f, w_f), \boldsymbol{\sigma}_f \approx -p\mathbf{I}\}$, respectively, where $p$ is the pore pressure and $\mathbf{I}$ the identity tensor.
## 3.2. Governing equations

Conserving mass in both the solid and fluid phases gives

$$ \frac{\partial \phi_p}{\partial t} + \nabla \cdot (\phi_p \mathbf{u}_p) = 0, \qquad (3.1a) $$

$$ -\frac{\partial \phi_p}{\partial t} + \nabla \cdot ((1 - \phi_p)\mathbf{u}_f) = 0. \qquad (3.1b) $$

Defining the Terzaghi effective stress tensor as $\sigma = \phi_p(\sigma_p - \sigma_f)$ (Wang 2001), a momentum balance yields

$$ \nabla \cdot (\phi \sigma_p + (1-\phi)\sigma_f) = \nabla \cdot (\phi \sigma_p - (1-\phi)p\mathbf{I}) = 0 \quad \Rightarrow \quad \nabla \cdot \sigma = \nabla p. \qquad (3.2) $$

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

$$ \sigma = \sigma (\nabla \xi), \qquad (3.3) $$

where $\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

$$ \mathbf{u}_p = \left( \frac{\partial}{\partial t} + \mathbf{u}_p \cdot \nabla \right) \xi. \qquad (3.4) $$

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 $\sigma$ obeys a linear elastic constitutive equation of the form

$$ \sigma(\nabla\xi) = \left( K - \frac{2G}{3} \right) (\nabla \cdot \xi)I + G(\nabla\xi + \nabla\xi^T), \qquad (3.5) $$

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, $\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

$$ (1 - \phi_p)(\mathbf{u}_p - \mathbf{u}_f) = \frac{\kappa}{\mu_f} \nabla p, \qquad (3.6) $$

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

$$ \kappa^* = \left( \frac{\phi_{0p}}{\phi_p} \right)^{\beta}, \qquad (3.7) $$

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

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

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 $u_p$ is defined as the material derivative of $\xi$, we have $\{\xi, \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 (12R_0^2H_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

$$A = \frac{a}{H_0}, \quad V = \frac{d^3}{12R_0^2 H_0}. \qquad (3.9a,b)$$

Note that, since the blister has positive curvature, $H_0 > a$ and thus $0 < 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

$$\frac{\partial \phi_p}{\partial t^*} + \frac{\partial}{\partial z^*} (\phi_p w_p^*) = O(\epsilon), \qquad (3.10a)$$

$$-\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, \qquad (3.10b)$$

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

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

$$\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 \xi^*}{\partial r^* \partial z^*} + O(\epsilon), \qquad (3.10e)$$

$$\frac{\partial P^*}{\partial z^*} = \frac{\partial^2 \xi^*}{\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), \qquad (3.10f)$$

where

$$u_p^* = \frac{\partial \xi^*}{\partial t^*} + w_p^* \frac{\partial \xi^*}{\partial z^*} + O(\epsilon), \qquad (3.10g)$$

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

### 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

$$w_f^* = 0 \text{ and } \xi^* = \zeta^* = 0 \quad \Rightarrow \quad u_p^* = w_p^* = 0 \quad \text{at } z^* = 0. \qquad (3.11a)$$

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

$$w_p^*|_{h^*} = \frac{\partial h^*}{\partial t^*} + \epsilon u_p^*|_{h^*} \frac{\partial h^*}{\partial r^*}, \qquad (3.11b)$$

$$w_f^*|_{h^*} = \frac{\partial h^*}{\partial t^*} + u_f^*|_{h^*} \frac{\partial h^*}{\partial r^*}. \qquad (3.11c)$$

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

$$\frac{\partial}{\partial t^*} \left( \int_0^{h^*} \phi_p \, dz^* \right) = O(\epsilon). \qquad (3.11d)$$
The stress-free boundary at $z^* = h^*$ requires

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

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

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

$$ \int_{0}^{\infty} r^*(h^* - A) dr^* = V. \qquad (3.11g) $$

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

$$ h^* \to A, \quad \phi_p \to \phi_{0p}, \quad \xi \to 0, \quad \zeta \to 0 \quad \text{as} \quad r^* \to \infty. \qquad (3.11h) $$

### 4. Perturbation expansion in the aspect ratio

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

$$ \frac{\partial^2 \xi^*}{\partial z^{*2}} = 0 + O(\epsilon), \qquad (4.1) $$

with boundary conditions

$$ \xi^*|_0 = \left. \frac{\partial \xi^*}{\partial z^*} \right|_{h^*} = 0 \Rightarrow \xi^* = 0 + O(\epsilon), \qquad (4.2) $$

(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

$$ \frac{\partial P^*}{\partial z^*} = 0 + O(\epsilon) \quad \Rightarrow \quad P^* = P_0^* + \left. \frac{\partial \zeta^*}{\partial z^*} \right|_{h^*} + O(\epsilon). \qquad (4.3) $$

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

$$ \frac{\partial^2 \zeta^*}{\partial z^{*2}} = 0 + O(\epsilon). \qquad (4.4) $$

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

$$ \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). \qquad (4.5) $$

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

$$ \zeta^* = C_{\zeta} z^* + O(\epsilon), \qquad (4.6) $$

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

$$ \zeta^* = z^* \left( 1 - \frac{A}{h^*} \right) + O(\epsilon). \qquad (4.7) $$
Furthermore, substituting this into (3.10h) and (4.3) yields, respectively,

$$w_p^* = \frac{z^*}{h^*} \frac{\partial h^*}{\partial t^*} + O(\epsilon), \qquad (4.8)$$

$$P^* = P_0^* + 1 - \frac{A}{h^*} + O(\epsilon). \qquad (4.9)$$

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

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

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

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

$$\phi_p = \phi_{0p} \int \left[ \left( \frac{A}{h^*} \right)^{1+\nu} z^{*\nu} \right] \Phi \, d\nu + O(\epsilon), \qquad (4.11)$$

where $\Phi = \Phi(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

$$\phi_p = \frac{A\phi_{0p}}{h^*} + O(\epsilon). \qquad (4.12)$$

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 $\bar{h}^* = h^* - A$, expressible as

$$\frac{\partial \bar{h}^*}{\partial t^*} = \frac{1}{r^*} \frac{\partial}{\partial r^*} \left( r^* \int_0^{\bar{h}^*+A} \kappa^* \frac{\partial P^*}{\partial r^*} \mathrm{d}z^* \right) = \frac{A}{r^*} \frac{\partial}{\partial r^*} \left( \frac{r^* \kappa^*}{\bar{h}^*+A} \frac{\partial \bar{h}^*}{\partial r^*} \right), \qquad (4.13)$$

with $A \le \bar{h}^* \le 1$ for $0 \le r^* < \infty$, where we have used (4.12) to write $\kappa^*(\phi)$ in terms of $\bar{h}^*$ as $\kappa^* = (1 + \bar{h}^*/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(\bar{h}^*) \ll 1$, the deformation height profile $\bar{h}^* = \bar{h}^*_s$ satisfies

$$\frac{\partial \bar{h}^*_s}{\partial t^*} = \nabla^2 \bar{h}^*_s, \qquad (4.14a)$$

with volume conservation condition

$$\int_0^\infty r^* \bar{h}^*_s dr^* = V. \qquad (4.14b)$$
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

$$ -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). \qquad (4.15) $$

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

$$ \frac{\partial}{\partial t^*}((R_s^*)^2) = \lambda \implies R_s^* = (1 + \lambda t^*)^{1/2}, \qquad (4.16) $$

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

$$ \begin{aligned} -\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 \implies \quad \frac{\partial f_s}{\partial \eta} = -\frac{\lambda \eta}{2} f_s \\ \implies \quad f_s &= (1 - \mathcal{A}) \exp \left( -\frac{\lambda \eta^2}{4} \right). \end{aligned} \qquad (4.17) $$

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

$$ h_s^* - \mathcal{A} = \frac{1}{(R_s^*)^2} \exp \left( -\frac{(1-\mathcal{A})}{2\mathcal{V}} \left( \frac{r_s^*}{R_s^*} \right)^2 \right), \qquad (4.18) $$

where

$$ R_s^* = \left( 1 + \frac{2(1 - \mathcal{A})t^*}{\mathcal{V}} \right)^{1/2}. \qquad (4.19) $$

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

$$ \overline{h}^* \gg \mathcal{A}. \qquad (4.20) $$

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

$$ \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), \qquad (4.21a) $$

with volume conservation condition

$$ \int_0^\infty r^* \overline{h}_l^* dr^* = \mathcal{V}, \qquad (4.21b) $$

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

$$ -(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). \quad (4.22) $$

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

$$ \frac{\partial}{\partial t^*}((R_l^*)^{2\beta}) = \lambda \implies R_l^* = (1 + \lambda t^*)^{1/2\beta}, \quad (4.23) $$

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

$$ \begin{aligned} -\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) &\implies \frac{\partial}{\partial \eta} (f_l^{\beta-1}) &= -\frac{\lambda \eta (\beta-1) \mathcal{A}^{\beta-1}}{2\beta}, \\ &\implies f_l = (1 - \mathcal{A})(1 - C\eta^2)^{1/(\beta-1)}, \end{aligned} \quad (4.24) $$

where $C$ satisfies

$$ C = \frac{\lambda(\beta - 1)\mathcal{A}^{\beta-1}}{4\beta(1 - \mathcal{A})^{\beta-1}}. \quad (4.25) $$

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{\mathcal{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

$$ \int_0^{R_{0l}^*} r^* \overline{h_l^*} dr^* = \mathcal{V}, \quad (4.26) $$

and we find

$$ \lambda_{\beta>1} = \frac{2(1 - \mathcal{A})^{\beta}}{\mathcal{V}\mathcal{A}^{\beta-1}}. \quad (4.27) $$

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

$$ \lambda_{\beta<1} = \frac{2(1 - \mathcal{A})^{\beta}}{\mathcal{V}\mathcal{A}^{\beta-1}}. \quad (4.28) $$

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

$$ h_l^* - \mathcal{A} = \frac{(1 - \mathcal{A})}{(R_l^*)^2} \left( 1 - \frac{(1 - \mathcal{A})(\beta - 1)}{2\beta\mathcal{V}} \left( \frac{r^*}{R_l^*} \right)^2 \right)^{1/(\beta-1)}, \quad (4.29) $$

with $0 \le r^* \le \infty$ if $\beta < 1$ and $0 \le r^* \le R_l^*$ if $\beta > 1$ and

$$ R_l^* = \left( 1 + \frac{2t^*(1 - \mathcal{A})^\beta \mathcal{A}^{1-\beta}}{\mathcal{V}} \right)^{1/(2\beta)}. \quad (4.30) $$
Note that, as for $\lambda, \bar{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$, $\bar{h}_l^*$ collapses to the small deformation similarity solution $\bar{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^* \to \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(\bar{h}^*) \gg 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 $A$. For general $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 $Q: 0 < r^* < 100$, bounded by no-flux boundary conditions. The time derivative is discretised using the backward Euler scheme

$$ \frac{\partial h^*}{\partial t^*} \approx \frac{(h^* - h_{-1}^*)}{\delta t^*}, \qquad (4.31) $$

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

$$ \int r^* Q h_{-1}^* dr^* = A^{1-\beta} \delta t^* \int r^*(h^*)^{\beta-1} \frac{\partial h^*}{\partial r^*} \frac{\partial Q}{\partial r^*} dr^* + \int r^* h^* Q dr^*, \qquad (4.32) $$

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^*(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.

$$ R^* = \{r^*: h^*(r^*) - A = \frac{1}{2}(\max(h^*(\tilde{r}^*)) - A)\}. \qquad (4.33) $$

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 $A$ (darker blue curves denote larger $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 $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(\bar{h}^*) \gg A$.
<figure><img src="image_8.png" /></figure>

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 $V = 1 - A$ and $A \in [0.001, 0.01, 0.1, 0.5, 0.9]$ (light blue to dark blue), namely darker shades of blue denote larger $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 \le A \le 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, A, t_0$}.

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 = \Omega R_0^2$, where the experiment-dependent part $R_0^2$ is found through image analysis while the experiment-independent part $\Omega = \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 $\Omega$ 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

$$ \sum_i \left\{ \text{avg}_j \left( \left[ \max(h_e(t_j)) - \max\left(h_\beta\left(\frac{t_j}{\Omega R_0^2}\right)\right) \right] / \max(h_e(t_j)) \right)^2 \right\}_i \quad (5.1) $$

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 \le \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 $\Omega = 6.72 \times 10^9 \text{ s m}^{-2}$, which corresponds to poroelastic time scales $t_0 = \Omega R_0^2 = 4.14 - 36.79 \text{ 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.

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 \le \beta \le 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), $\Omega$ would have decreased slightly to $5.82 \times 10^9 \text{ s m}^{-2}$ (a change of around 15%). This
<figure><img src="image_9.png" /><figcaption>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 µ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 $\Omega = 6.72 \times 10^9 \text{ s 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.</figcaption></figure>

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 $\Omega = 6.72 \times$
$10^9 \text{ s 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 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 µl drop (the darkest blue curve) and the 30 µ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 µ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.

Discussed in more detail below in Appendix A, the imaging system cannot detect very small swelling (estimated to be less than 45 µ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$ 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.
Droplet absorption and spreading into thin hydrogel layers

<figure><img src="image_10.png" /></figure>

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 µ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$ 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.

## 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
*Droplet absorption and spreading into thin hydrogel layers*

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).

**Supplementary material.** Experimental blister shape profiles are available at https://doi.org/10.17863/CAM.99950.

**Acknowledgements.** The authors would like to thank the late Dr H. McEvoy (Dstl, Porton Down) and Professor M. Grae Worster (Cambridge) for proposing the problem. The authors also thank Dr S. Marriott (Dstl, Porton Down) for his help with the Karl-Fischer titration, M.G.W. for his comments on the manuscript and Mr P. Mitton and the late Mr C. Summerfield for their advice and help during the experimental development.

**Funding.** M.A.E., G.T.F., J.R.L. and S.B.D. acknowledge funding from Dstl under extramural research agreement DSTLX-1000138254. G.T.F. acknowledges a Doctoral Training Fellowship from the Engineering and Physical Sciences Research Council.

**Declaration of interests.** The authors report no conflict of interest.

**Author ORCIDs.**

● Merlin A. Etzold https://orcid.org/0000-0002-4914-4908;

● George T. Fortune https://orcid.org/0000-0003-0817-9271;

● Julien R. Landel https://orcid.org/0000-0003-3159-8749;

● Stuart B. Dalziel https://orcid.org/0000-0002-8487-2038.

**Author contributions.** M.A.E., S.B.D. and J.R.L. identified the problem's practical application and acquired the funding. M.A.E. developed the experiments with input from S.B.D. and J.R.L. G.T.F. developed the mathematical model with input from MAE as the experiments developed. G.T.F. and M.A.E. matched the model with experimental data. M.A.E. and G.T.F. prepared the first draft of the manuscript. All authors contributed on further revisions of the manuscript.

**Rights retention.** For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

**Appendix A. Uncertainty analysis**

We assume our imaging system to be accurate within one pixel. Hence, at a pixel pitch of 5 µm and a magnification of 0.4, this leads to imaging uncertainties $δr = δz = 12.5$ µm. The particular configuration of the side view imaging leads to further uncertainties on the determination of $h_a$ due both to surface roughness (imperfections in the manufacturing process) and alignment errors between the camera and the sample (illustrated in figure 10a). Figure 10(b) shows how a small misalignment of angle $δα$ causes some of the swelling to be hidden below the edge of the sample. This introduces an uncertainty in the position of the upper edge of the hydrogel $δa$. Assuming that the alignment between the camera and the sample is accurate within 0.05°, we find that, for 26 mm wide samples, the edge position of the unswollen hydrogel appears to be located $δa^+ = +22.7$ µm too high. We thus estimate the uncertainty of the blister profile as $δh_a = (+δa, -δz) = (+22.7, -12.5)$ µm.
Droplet absorption and spreading into thin hydrogel layers

<figure><img src="image_11.png" /></figure>

Figure 10. Illustration of the uncertainties inherent in the imaging. (a) Side view of the experimental configuration, showing how a misalignment (exaggerated for clarity) causes part of the droplet to be hidden from the camera by the edge of the hydrogel sheet. (b) Image (linear projection) of the droplet on the hydrogel sheet as seen by the camera. Dashed lines are not visible for the camera and give rise to an uncertainty in the position of the hydrogel surface $\delta a$, that creates a corresponding uncertainty $\delta R$ in the determination of $R$. The uncertainty due to finite pixel resolution in the radial and vertical directions are defined as $\delta r$ and $\delta z$, respectively.

Furthermore, figure 10(b) illustrates the potentially significant effect that $\delta h_a$ may have on the uncertainty in $R$. We estimate this from the slope of $h_a$ at $R$ as

$$ \delta R \approx \left. \frac{(\delta h)}{\left(\frac{\partial h_a}{\partial x}\right)} \right|_{x=r_0+R} , \qquad (A1) $$

where the derivative is approximated for each frame using a fourth-order central difference scheme. Utilising (A1), we see that if the slope of the blister surface decreases over time, the uncertainties grow over time. This makes the experimental determination of $R$ difficult at late stages when the blister becomes flat.

The uncertainty of the determination of $\delta h_a$ also affects the computation of $V$, (illustrated by the hatched region in figure 10(b)). We estimate the uncertainty $\delta V$ on the determination of the blister volume as

$$ \delta V \approx \pi \int_{x_{mn}}^{x_{mx}} (x - r_0) \delta h_a(x) \, dx, \qquad (A2) $$

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

$$ \{x_m \le x \le x_{mx} : h_a(x - r_0) > \delta z\}. \qquad (A3) $$

This approximates the hatched region given in figure 10(b). We did not attempt to extrapolate the exact edge position based on the slope since the slope of $h(x, t)$ becomes very small at the leading edge (indicated by both our theory and experimental data; see
<figure><img src="image_12.png" /><figcaption>Figure 11. Absorption of a dyed (0.01 % by weight methylene blue) droplet (100–200 µl) observed in plan view against a bright background. (a) Immediately after droplet placement the only visible structure is a pattern arising from gauze that was embedded within the hydrogel during the manufacturing process. (b–e) As the water droplet is absorbed into the hydrogel, a surface instability emerges, growing coarser as time progresses. (f) The instability has fully decayed as the droplet and hydrogel have dried out. Panels show (a) 0 s, (b) 13 s, (c) 87 s, (d) 493 s, (e) 1037 s, (f) 2083 s.</figcaption></figure>

figure 2 at late times). Similar to the error analysis for *R*, the above analysis shows that accurately determining *V* becomes progressively harder as the blister spreads, since $h_a$ decreases over time whilst the area covered by the blister, defined in (2.6), increases with time. We also note that the accurate determination of the volume relies on the assumption of cylindrical symmetry of the blister; we comment on this in Appendix C.

## Appendix B. Surface instability

Here, we briefly report additional footage to illustrate the development of the crumpling instability over time. We found that the instability becomes visible if a droplet of methylene blue solution (0.01 % by weight) is used rather than pure water. Figure 11 shows typical observations for the absorption of a single droplet with volume between 100 and 200 µl imaged in plan view against a diffusive white illumination using the same grey-scale camera and telecentric lens as the main experiment (cf. § 2.3).

Figure 11(a) shows the droplet just after making contact. The droplet appears dark against the overexposed background due to the presence of the dye. Within the droplet, the gauze layer embedded into the hydrogel is visible as random fibre pattern. Within seconds after droplet placement, shown in figure 11(b), a second regular pattern (dark lines) becomes visible. This pattern then coarsens over time (figure 11c,d). In figure 11(e), the edge of the droplet has become dry due to evaporation and absorption, while the pattern in the centre has coarsened even further.
Droplet absorption and spreading into thin hydrogel layers

<figure><img src="image_13.png" /></figure>

Figure 12. Raw data with corresponding error bars showing for a range of experiments how the apparent volume of the blister (found by direct integration) varies as a function of time. Darker blue curves denote experiments with larger water droplets and black dotted lines indicate the nominal droplet volume. Error bars are generated through an uncertainty analysis that is described in § 2.4 and Appendix A.

From the supplementary video, it can be deduced that the black lines at this point in the experiment are caused by solid dye that has precipitated in the creases of the instability. This movement of dye particulates in the creases indicates a complex flow pattern in the supernatant fluid at this stage as the instability vanishes (figure 11e, f).

Qualitatively, the pattern observed during this process corresponds to that observed in experiments using pure water (cf. § 2.5). We therefore propose that this sequence qualitatively represents the evolution of the instability during droplet absorption. At the same time, these experiment highlights the role specific chemical interactions play in hydrogel swelling, an area which certainly warrants further exploration.

# Appendix C. Conservation of water volume

## C.1. *Experimental test of conservation of volume*

Conservation of volume is an important assumption for our theoretical modelling. In this section we briefly investigate whether the experimental measurements are consistent with this assumption. In figure 12, we report the blister volume as a function of time, found by integrating the experimental height profile (2.3). Here, the nominal volume of the dataset is shown as a horizontal dotted line. Full lines map how the volume of the swollen region evolves for the different experiments. Corresponding error bars represent our uncertainty estimates (see (2.5) and Appendix A for further details).

Note that the integration (2.3) of the experimental measurement of $h_a$ at the start of all experiments overestimates the volume of the blister. As $h_a$ is effectively the silhouette of a projection of the blister, the localised surface depressions associated with the creasing of the surface hydrogel are hidden from the side-view camera and thus not fully accounted for.

For the smallest three droplets (namely those with volume 5, 10 and 25 µl), volume is largely conserved throughout the experiment. For the remaining larger droplets, we observe a continuous decrease in volume over the observed period. Our uncertainty
estimates capture the trend of this behaviour well. However, it does not accurately predict the magnitude of this decrease in volume for the two largest droplets.

However, as noted in the discussion of the error analysis given above in Appendix A, accurately determining the blister volume becomes progressively more difficult with increasing blister radius. Its accuracy relies strongly on the assumption of cylindrical symmetry. Given the instrumentation available at the time, the cylindrical symmetry of the droplet settling on the hydrogel could only be verified by eye and thus might not have been perfect in some cases. We also observed that the bounds of the leading edge of the droplet (see (2.6)) sometimes varied by up to 10%, indicating imperfect spherical symmetry. Hence, a plausible explanation to explain why smaller blisters seem to conserve their volume better than larger blisters is that smaller droplets experience stronger capillary forces and thus would settle more symmetrically on a surface compared with a big droplet.

In summary, conservation of volume for the smaller droplets and the success of our error analysis to explain the trends seen in figure 12 for larger droplets suggests that volume is nearly conserved. In particular, although we measure a significant volume loss, the majority of this loss can be attributed to measurement error rather than actual volume loss. However, for completeness, in Appendix C.2, we also estimated potential losses of water through the covering oil. From this analysis, we can conclude that it is unlikely that the swollen hydrogel lost a significant amount of water into the oil during the duration of an experiment.

## C.2. Analysis of water transport into or through the oil

Since there is a sparing solubility of water in the hydraulic oil, we investigated the possible water loss due to transport into the oil phase by making an order of magnitude estimation for this transport by considering the effect from diffusive mass flux in the z-direction. We assume that mass flux is due to diffusion only (no advective component) and that the oil is saturated with water at the interface to the swollen region of the hydrogel. We take the solubility limit as 1.1% water content (cf. § 2.2), which we believe to be an overestimation since the chemical potential of water in unsaturated hydrogel is less than of pure water (Etzold et al. 2021) and the oil is in equilibrium with the laboratory atmosphere far away from the hydrogel (estimated to be 0.67% water content by titration, see § 2.2). Hence, we model the difference in water volume concentration $\Delta c$ that drives water diffusion through the oil as

$$ \Delta c \approx \Delta\omega\rho_o, \qquad (C1) $$

where the difference in water mass fraction that drives diffusion $\Delta\omega = 0.43\%$, the density of mineral oil is estimated as $\rho_o = 833 \text{ kg m}^{-3}$ and we have assumed both isochoric mixing and that the water concentration in the oil is small. At a fixed radial distance $r = r'$ and in the absence of buoyancy-driven motion (we expect the water content to increase the density of the oil), the concentration field $c = c(r', z, t, t')$ for water in the oil can be described using the similarity solution for one-dimensional diffusion in the z-direction (neglecting in-plane diffusion within the $(r, \theta)$ plane) to give

$$ c = \Delta c \operatorname{erfc}(z/\sqrt{4D(t-t')}), \qquad (C2a) $$

with corresponding diffusive water mass flux at $z = 0$

$$ j = \frac{\Delta c \sqrt{D}}{\sqrt{\pi(t-t')}}, \qquad (C2b) $$
where erfc(·) is the complementary error function, D is the diffusion coefficient for water in the oil and t' is the time at which the swelling front reaches the point r = r' i.e. r' = R(t'). The total water volume lost for a swollen blister thus becomes

$$ \Delta V = \frac{1}{\rho_w} \int_0^t \left[ A_0 j(t') + \frac{\partial A}{\partial t}(t') \otimes j(t') \right] dt', \quad (C3) $$

where we have assumed that the blister has a circular base and surface area A(t) with A(0) = A₀. The convolution integral denoted by ⊗ accounts for the fact that diffusion starts only when the hydrogel starts to swell. For simplicity, if we assume the radius grows linearly with a growth rate χ, then (C3) predicts

$$ \begin{aligned} \Delta V &= \frac{2\rho_o}{\rho_w} \sqrt{\pi D \Delta\omega} \left[ R_0^2 \sqrt{t} + \int_0^t \int_0^t \chi \frac{(R_0 + \chi t'')}{\sqrt{t' - t''}} dt'' dt' \right] \\ &= \frac{2\rho_o}{\rho_w} \sqrt{\pi D \Delta\omega} \left[ R_0^2 \sqrt{t} + \frac{4}{15} \chi t^{3/2} (5R_0 + 2\chi t) \right]. \end{aligned} \quad (C4) $$

Utilising the experimental data for the largest blister (see figure 2 and table 1) with initial radius R₀ = 5 mm, approximating the radius of the hydrogel area covered by the blister as growing linearly with time yields an effective χ = 0.5 mm h⁻¹. This estimation was calculated based on imaging data (e.g. compare figures 2c and 2k) rather than using R evaluated from experimental measurements using (2.2) since for flat blisters our R underestimates the edge position. We found no data for the diffusion of water in HySpin AWS 32 and established that such data are not widely available for mineral oils. However, given its relatively high viscosity of μ₀ = 32 cSt at 40 °C, the water diffusion coefficient in oil is likely to be an order of magnitude less than the self-diffusion coefficient of bulk water, D = 3 × 10⁻⁹ m² s⁻¹. Furthermore, Hilder & van den Tempe (1971) reported water diffusion coefficients in groundnut oil (50 times more viscous than water) and in kerosene (twice as viscous as water) as 2.5 × 10⁻¹⁰ m² s⁻¹ and 8.5 × 10⁻¹⁰ m² s⁻¹ respectively. Hence, taking the diffusion constants in groundnut oil and in kerosene as lower and upper bounds respectively for the diffusion constant of water in our hydraulic oil, (C4) predicts diffusive volume losses of {0.3 – 0.6}, {1.2 – 2.2} and {2.7 – 5.0} μl after one, five and ten hours, respectively. As this is almost an order of magnitude too small to account for the apparent volume loss found in figure 12.

For completeness, a supplementary experiment was conducted where a 2 μl droplet of water was placed directly onto a microscope slide that was immersed into an oil bath with no hydrogel present. This droplet reached a sessile state after 25 h. Whilst the accuracy of our image analysis of this experiment was impacted from resolution issues in locating the diffuse edges of the droplet, we found that the water loss of the droplet was less than 0.1 μl over a time period of approximately 120 h. The surface area of the droplet remained approximately constant and was estimated to be 4.6 mm². Assuming a steady state mass loss from the droplet into the atmosphere through a 1 cm deep layer of oil, this analysis yields a diffusion coefficient of 1–2 × 10⁻¹⁰ m² s⁻¹. By approximating this axisymmetric three-dimensional diffusion problem using the equations above, assuming mass transport into an infinite oil bath yields an even lower diffusion coefficient of 7 × 10⁻¹¹ m² s⁻¹. We note that these droplets remained in the oil bath for approximately one month without vanishing.

In summary, it appears unlikely that, based on the water diffusion coefficient from the literature, the swollen hydrogel lost a significant amount of water into the oil during the duration of an experiment (corroborated by a supplementary experiment testing
<table><caption>Table 1. Experimental $t_s$ from which the initial condition for the numerical solution was obtained with corresponding values {$H_0, R_0, A, t_0$}. The second column refers to the frame number in the corresponding attached raw dataset.</caption><thead><tr><th>V (μl)</th><th>Frame</th><th>t<sub>s</sub> (min)</th><th>H<sub>0</sub> (mm)</th><th>R<sub>0</sub> (mm)</th><th>A</th><th>t<sub>0</sub> (h)</th></tr></thead><tbody><tr><td>100</td><td>180</td><td>82.0</td><td>2.61</td><td>4.44</td><td>0.43</td><td>36.79</td></tr><tr><td>50</td><td>220</td><td>39.5</td><td>2.44</td><td>3.27</td><td>0.46</td><td>19.96</td></tr><tr><td>30</td><td>782</td><td>21.0</td><td>2.28</td><td>2.78</td><td>0.49</td><td>14.42</td></tr><tr><td>25</td><td>235</td><td>26.0</td><td>2.18</td><td>2.61</td><td>0.51</td><td>12.71</td></tr><tr><td>10</td><td>109</td><td>20.0</td><td>1.82</td><td>1.94</td><td>0.62</td><td>7.02</td></tr><tr><td>5</td><td>358</td><td>9.5</td><td>1.66</td><td>1.49</td><td>0.68</td><td>4.14</td></tr></tbody></table>

the stability of a small water droplet in the same conditions). Also, the experimental observation that volume appears to be conserved for smaller droplets suggests that the extra volume loss apparent above in figure 12 for larger droplets beyond that attributed to measurement error is likely due to the initial blister profile not being completely axisymmetric, resulting in spreading that is not quite radial.

## Appendix D. Additional experimental information

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

## REFERENCES

* AL-AMODI, A. & HILL, R.J. 2022 Streaming potentials of hyaluronic acid hydrogel films. *Langmuir* **38** (44), 13370–13381.

* ALNÆS, M.S., BLECHTA, J., HAKE, J., JOHANSSON, A., KEHLET, B., LOGG, A., RICHARDSON, C., RING, J., ROGNES, M.E. & WELLS, G.N. 2015 The FEniCS project version 1.5. *Arch. Numer. Softw.* **3** (100). https://doi.org/10.11588/ans.2015.100.20553.

* BARRY, S.I., MERCER, G.N. & ZOPPOU, C. 1997 Deformation and fluid flow due to a source in a poro-elastic layer. *Appl. Math. Model.* **21** (11), 681–689.

* BERTRAND, T., PEIXINHO, J., MUKHOPADHYAY, S. & MACMINN, C.W. 2016 Dynamics of swelling and drying in a spherical gel. *Phys. Rev. Appl.* **6**, 064010.

* BIOT, M.A. 1941 General theory of three-dimensional consolidation. *J. Appl. Phys.* **12**, 155–164.

* BOUKLAS, N. & HUANG, R. 2012 Swelling kinetics of polymer gels: comparison of linear and nonlinear theories. *Soft Matt.* **8**, 8194–8203.

* BOUKLAS, N., LANDIS, C.M. & HUANG, R. 2015 A nonlinear, transient finite element method for coupled solvent diffusion and large deformation of hydrogels. *J. Mech. Phys. Solids* **79**, 21–43.

* BUTLER, M.D. & MONTENEGRO-JOHNSON, T.D. 2022 The swelling and shrinking of spherical thermo-responsive hydrogels. *J. Fluid Mech.* **947**, A11.

* CAI, S. & SUO, Z. 2012 Equations of state for ideal elastomeric gels. *Europhys. Lett.* **97** (3), 34009.

* CHARRAS, G.T., MITCHISON, T.J. & MAHADEVAN, L. 2009 Animal cell hydraulics. *J. Cell Sci.* **122** (18), 3233–3241.

* CHESTER, S.A. & ANAND, L. 2010 A coupled theory of fluid permeation and large deformations for elastomeric materials. *J. Mech. Phys. Solids* **58** (11), 1879–1906.

* CONT, A., ROSSY, T., AL-MAYYAH, Z. & PERSAT, A. 2020 Biofilms deform soft surfaces and disrupt epithelia. *eLife* **9**, e56533.

* DALZIEL, S.B. 2017 *DigiFlow User Guide*. Dalziel Research Partners.

* DELAVOPIÈRE, J., HEURTEFEU, B., TEISSEIRE, J., CHATEAUMINOIS, A., TRAN, Y., FERMIGIER, M. & VERNEUIL, E. 2018 Swelling dynamics of surface-attached hydrogel thin films in vapor flows. *Langmuir* **34** (50), 15238–15244.
Droplet absorption and spreading into thin hydrogel layers

DENG, P., YAO, L., CHEN, J., TANG, Z. & ZHOU, J. 2022 Chitosan-based hydrogels with injectable, self-healing and antibacterial properties for wound healing. *Carbohydr. Polym.* **276**, 118718.

DOI, M. 2009 Gel dynamics. *J. Phys. Soc. Japan* **78** (5), 052001.

DRURY, J.L. & MOONEY, D.J. 2003 Hydrogels for tissue engineering: scaffold design variables and applications. *Biomaterials* **24**, 4337-4351.

ENGELSBERG, M. & BARROS, W. JR. 2013 Free-evolution kinetics in a high-swelling polymeric hydrogel. *Phys. Rev. E* **88** (6), 062602.

ETZOLD, M.A., LINDEN, P.L. & WORSTER, M.G. 2021 Transpiration through hydrogel. *J. Fluid Mech.* **925**, A8.

FLORY, P.J. 1953 *Principles of Polymer Chemistry*. Cornell University Press.

FORTUNE, G.T., OLIVEIRA, N.M. & GOLDSTEIN, R.E. 2022 Biofilm growth under elastic confinement. *Phys. Rev. Lett.* **128** (17), 178102.

GIBSON, R.E., SCHIFFMAN, R.L. & PU, S.L. 1970 Plane strain and axially symmetric consolidation of a clay layer on a smooth impervious base. *Q. J. Mech. Appl. Maths* **23** (4), 505-520.

GRATTONI, C.A., AL-SHARJI, H.H., YANG, C., MUGGERIDGE, A.H. & ZIMMERMAN, R.W. 2001 Rheology and permeability of crosslinked polyacrylamide gel. *J. Colloid Interface Sci.* **240** (2), 601–607.

GUVENDIREN, M., BURDICK, J.A. & YANG, S. 2010 Solvent induced transition from wrinkles to creases in thin film gels with depth-wise crosslinking gradients. *Soft Matt.* **6**, 5795–5801.

HARRIS, C.R., et al. 2020 Array programming with NumPy. *Nature* **585** (7825), 357-362.

HEWITT, D.R., NEUFELD, J.A. & BALMFORTH, N.J. 2015 Shallow, gravity-driven flow in a poro-elastic layer. *J. Fluid Mech.* **778**, 335–360.

HILDER, M.H. & VAN DEN TEMPE, M. 1971 Diffusivity of water in groundnut oil and paraffi oil. *J. Appl. Chem. Biotechnol.* **21** (6), 176-178.

HOFFMAN, A.S. 2012 Hydrogels for biomedical applications. *Adv. Drug Deliv. Rev.* **64**, 18–23.

HONG, W., ZHAO, X., ZHOU, J. & SUO, Z. 2008 A theory of coupled diffusion and large deformation in polymeric gels. *J. Mech. Phys. Solids* **56** (5), 1779–1793.

IONOV, L. 2014 Hydrogel-based actuators: possibilities and limitations. *Mater. Today* **17** (10), 494-503.

JENSEN, O.E., GLUCKSBERG, M.R., SACHS, J.R. & GROTBERG, J.B. 1994 Weakly nonlinear deformation of a thin poroelastic layer with a free surface. *J. Appl. Mech.* **543**, 293-302.

KWON, G.H., JEONG, G.S., PARK, J.Y., MOON, J.H. & LEE, S.H. 2011 A low-energy-consumption electroactive valveless hydrogel micropump for long-term biomedical applications. *Lab on a Chip* **11** (17), 2910–2915.

LAFTAH, W.A., HASHIM, S. & IBRAHIM, A.N. 2011 Polymer hydrogels: a review. *Polym. Plast. Technol. Engng* **50** (14), 1475-1486.

LANDEL, J.R., THOMAS, A.L., MCEVOY, H. & DALZIEL, S.B. 2016 Convective mass transfer from a submerged drop in a thin falling film. *J. Fluid Mech.* **789**, 630–668.

LEE, K.Y. & MOONEY, D.J. 2001 Hydrogels for tissue engineering. *Chem. Rev.* **101** (7), 1869-1880.

LI, J., HU, Y., VLASSAK, J.J. & SUO, Z. 2012 Experimental determination of equations of state for ideal elastomeric gels. *Soft Matt.* **8** (31), 8121–8128.

LI, R., SHI, Y., ALSAEDI, M., WU, M., SHI, L. & WANG, P. 2018 Hybrid hydrogel with high water vapor harvesting capacity for deployable solar-driven atmospheric water generator. *Environ. Sci. Technol.* **52** (19), 11367-11377.

LOGG, A., MARDAL, K.-A. & WELLS, G.N. 2012 *Automated Solution of Differential Equations by the Finite Element Method*. Springer.

MACMINN, C.W., DUFRESNE, E.R. & WETTLAUFER, J.S. 2015 Fluid-driven deformation of a soft granular material. *Phys. Rev. X* **5**, 011020.

MADHUSUDANAN, P., RAJU, G. & SHANKARAPPA, S. 2020 Hydrogel systems and their role in neural tissue engineering. *J. R. Soc. Interface* **17**, 20190505.

NORDBOTTEN, J.M. & CELIA, M.A. 2006 Similarity solutions for fluid injection into confined aquifers. *J. Fluid Mech.* **561**, 307–327.

ONUKI, A. 1988 Pattern formation in gels. *J. Phys. Soc. Japan* **57** (3), 703–706.

OZCELIK, B. 2016 Degradable hydrogel systems for biomedical applications. In *Biosynthetic Polymers for Medical Applications* (ed. L. Poole-Warren, P. Martens & R. Green), pp. 172–188. Elsevier.

PHADNIS, A., MANNING, K.C., SANDERS, I., BURGIN, T.P. & RYKACZEWSKI, K. 2018 Droplet-train induced spatiotemporal swelling regimes in elastomers. *Soft Matt.* **14**, 5869–5877.

PHILLIPS, O.M. 2009 *Geological Fluid Dynamics: Sub-Surface Flow and Reactions*. Cambridge University Press.

PORTER, T., STEWART, R., REED, J. & MORTON, K. 2007 Models of hydrogel swelling with applications to hydration sensing. *Sensors* **7** (9), 1980–1991.
RICHTER, A., KLATT, S., PASCHEW, G. & KLENKE, C. 2009 Micropumps operated by swelling and shrinking of temperature-sensitive hydrogels. *Lab on a Chip* **9** (4), 613–618.

RUTQVIST, J. 2012 The geomechanics of CO2 storage in deep sedimentary formations. *Geotechnol. Geol. Engng* **30**, 525–551.

SEMINARA, A., ANGELINI, T.E., WILKING, J.N., VLAMAKIS, H., EBRAHIM, S., KOLTER, R., WEITZ, D.A. & BRENNER, M.P. 2012 Osmotic spreading of *Bacillus subtilis* biofilms driven by an extracellular matrix. *Proc. Natl Acad. Sci.* **109**, 1116–1121.

SPENCER, A. 2004 *Continuum Mechanics*. Dover.

STUBBE, B., MIGNON, A., VAN DAMME, L., CLAES, K., HOEKSEMA, H., MONSTRY, S., VAN VЛИERBERGHE, S. & DUBRUEL, P. 2021 Photo-crosslinked gelatin-based hydrogel films to support wound healing. *Macromol. Biosci.* **21**, 2100246.

TANAKA, T. & FILLMORE, D.J. 1979 Kinetics of swelling of gels. *J. Chem. Phys.* **70** (3), 1214–1218.

TANG, G., ZHOU, B., LI, F., WANG, W., LIU, Y., WANG, X., LIU, C. & YE, X. 2020 Advances of naturally derived and synthetic hydrogels for intervertebral disk regeneration. *Front. Bioeng Biotechnol.* **8**, 745.

TOKITA, M. & TANAKA, T. 1991 Friction coefficient of polymer networks of gels. *J. Chem. Phys.* **95** (6), 4613–4619.

TORTORELLA, S., INZALACO, G., DAPPORTO, F., MATURI, M., SAMBRI, L., BURATTI, V.V., CHIARIELLO, M., FRANCHINI, M.C. & LOCATELLI, E. 2021 Biocompatible pectin-based hybrid hydrogels for tissue engineering applications. *New J. Chem.* **45**, 22386.

TRELOAR, L.R.G. 1975 *The Physics of Rubber Elasticity*. Oxford University Press.

VAN DE VELDE, P., DERVAUX, J., PROTIÈRE, S. & DUPRAT, C. 2022 Spontaneous localized fluid release on swelling fibres. *Soft Matt.* **18** (24), 4565–4571.

VAN DE VELDE, P., PROTIÈRE, S. & DUPRAT, C. 2021 Dynamics of drop absorption by a swelling fiber. *Soft Matt.* **17**, 6168–6175.

VARADY, M.J., PEARL, T.P., STEVENSON, S.M. & MANTOOTH, B.A. 2016 Decontamination of VX from silicone: characterization of multicomponent diffusion effects. *Ind. Engng Chem. Res.* **55** (11), 3139–3149.

VERVOORT, S. 2006 Behaviour of hydrogels swollen in polymer solutions under mechanical action. PhD thesis, École Nationale Supeérieure des Mines de Paris, France.

VIRTANEN, P., et al. 2020 SciPy 1.0: fundamental algorithms for scientific computing in Python. *Nat. Meth.* **17**, 261–272.

WALSH, C., STRIDE, E., CHEEMA, U. & OVENDEN, N. 2017 A combined three-dimensional in vitro – in silico approach to modelling bubble dynamics in decompression sickness. *J. R. Soc. Interface* **14**, 20170653.

VAN DER WALT, S., SCHÖNBERGER, J.L., NUNEZ-IGLESIAS, J., BOULOGNE, F., WARNER, J.D., YAGER, N., GOUILLART, E. & YU, T. 2014 scikit-image: image processing in Python. *PeerJ* **2**, e453.

WANG, H.F. 2001 *Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology*. Princeton University Press.

WICHTERLE, O. & LIM, D. 1960 Hydrophilic gels for biological use. *Nature* **185**, 117–118.

YOON, J., CAI, S., SUO, Z. & HAYWARD, R.C. 2010 Poroelastic swelling kinetics of thin hydrogel layers: comparison of theory and experiment. *Soft Matt.* **6** (23), 6004–6012.

YU, C., MALAKPOOR, K. & HUYGHE, J.M. 2020 Comparing mixed hybrid finite element method with standard FEM in swelling simulations involving extremely large deformations. *Comput. Mech.* **66**, 287–309.

ZHANG, Y.Y., ETZOLD, M.A. & LEFAUVE, A. 2021 Growth of gas-filled penny-shaped cracks in decompressed hydrogels. *Soft Matt.* **17** (4), 815–825.