2.2. Suppression of vapour transport

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

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

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

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

2.3. Imaging system and data analysis

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

DigiFlow was used for camera control (Dalziel Reference Dalziel2017). 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. Reference van der Walt, Schönberger, Nunez-Iglesias, Boulogne, Warner, Yager, Gouillart and Yu2014). Suitable parameters for the edge detector (radius of Gaussian filter, upper and lower gradient thresholds for edge detection) were manually set for each dataset using an interactive interface. To mitigate edge effects arising from the finite size of the hydrogel sheet, we truncated each time series when the swelling from the blister approached the edges. For each dataset, we extracted the position of the upper edge of the unswollen hydrogel before the droplet was placed (see figure 2a) and the upper edge of hydrogel and droplet once the droplet was placed (see figure 2c–i). For analysis, the blister height (vertical displacement) was defined as

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

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

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

2.4. Derived quantities

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

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

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

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

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

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

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

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

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

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

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

2.5. Observations

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

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

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

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

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

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

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

3.1. Problem description

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

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

3.2. Governing equations

Conserving mass in both the solid and fluid phases gives

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

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

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

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

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

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

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

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

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

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

where $K$ and $G$ are the osmotic and shear moduli, respectively (Doi Reference Doi2009; Etzold et al. Reference Etzold, Linden and Worster2021).

As shown by both Doi (Reference Doi2009) and Bouklas & Huang (Reference Bouklas and Huang2012), the linear elastic equation (3.5) can be obtained by linearisation of nonlinear constitutive equations similar to those used by Engelsberg & Barros (Reference Engelsberg and Barros2013), Bertrand et al. (Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016) 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 (Reference Doi2009) and Etzold et al. (Reference Etzold, Linden and Worster2021), 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 Reference Biot1941). In a hydrogel, $K$ represents the combination of osmotic (solution) effects and a contribution from the elasticity of the polymer matrix (Doi Reference Doi2009). It is therefore appropriate to see (3.5) as a simplified version of physically more encompassing nonlinear elastic models. Note that the form of (3.5) implies that, when the system is in its reference state, $\boldsymbol {\sigma }=0$. While this is not always physically appropriate, it is immaterial for the current analysis in this paper; we will return to this point in § 6.

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

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

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

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

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

where $\beta$ is an empirical parameter chosen by matching to experimental data, recalling that $\phi _{0p}$ is defined as the polymer volume fraction at the start of an experiment before the water drop has been added, i.e. at ambient laboratory conditions. For a purely cuboidal network swelled isotropically, previous researchers have proposed $\beta = 2/3$ (Etzold et al. Reference Etzold, Linden and Worster2021). Al-Amodi & Hill (Reference Al-Amodi and Hill2022) 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. (Reference Etzold, Linden and Worster2021). In contrast, experiments have measured values ranging between $\beta =1.5$ (Tokita & Tanaka Reference Tokita and Tanaka1991) and 1.85 (Grattoni et al. Reference Grattoni, Al-Sharji, Yang, Muggeridge and Zimmerman2001).

3.3. Dimensionless shallow-layer scalings

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

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

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

3.4. Non-dimensional governing equations

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

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

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

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

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

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

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

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

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

where

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

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

3.5. Boundary conditions

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.1. Self-similar solutions: small deformation limit

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

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

with volume conservation condition

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

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

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

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

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

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

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

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

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

where

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

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

4.2. Self-similar solutions: large deformation limit

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

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

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

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

with volume conservation condition

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

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

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

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

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

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

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

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

where $C$ satisfies

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

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

Thus, the volume conservation condition (4.21b) becomes

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

and we find

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

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

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

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

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

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

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

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

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

4.3. Numerical exploration

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

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

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

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

where $Q$ is a test function.

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

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

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

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

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

5.1. Fitting procedure

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

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

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

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

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

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

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

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

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

Whilst $\beta =2.25$ is considerably larger than the suggestion of $2/3$ made by Etzold et al. (Reference Etzold, Linden and Worster2021), it is more comparable to many previous experimental studies which found $1.5\leq \beta \leq 1.85$ (Tokita & Tanaka Reference Tokita and Tanaka1991; Grattoni et al. Reference Grattoni, Al-Sharji, Yang, Muggeridge and Zimmerman2001; Bertrand et al. Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016). 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. Reference Grattoni, Al-Sharji, Yang, Muggeridge and Zimmerman2001), $\varOmega$ would have decreased slightly to $5.82\times 10^{9}\,{\rm s}\,{\rm m}^{-2}$ (a change of around 15 %). This small change does not influence the characteristic trends that are observed as the system evolves over time. This insensitivity of the fit might be due to only modest changes of the observables over the observation period (factor 2). While we have not tried expressions for $\kappa ^{*}$ in terms of $\phi _p$ other than (3.7), at small $\phi _p$ the leading-order term in the Taylor expansion will be of the polynomial form we have chosen in our model (3.7).

5.2. Profile comparison

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

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

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

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

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

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

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

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

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