††thanks: A.P. and C.A. contributed equally to this work.††thanks: A.P. and
C.A. contributed equally to this work.††thanks: To whom correspondence may be
addressed.
Email<EMAIL_ADDRESS>or<EMAIL_ADDRESS>To whom
correspondence may be addressed.
Email<EMAIL_ADDRESS>or<EMAIL_ADDRESS>
# Obstructed swelling and fracture of hydrogels

Abigail Plummer  Princeton Center for Complex Materials, Princeton University,
Princeton, NJ 08544  Caroline Adkins  Department of Civil and Environmental
Engineering, Princeton University, Princeton, NJ 08544  Jean-François Louf
Department of Chemical and Biological Engineering, Princeton University,
Princeton, NJ 08544  Department of Chemical Engineering, Auburn University,
Auburn, AL 36849  Andrej Košmrlj  Department of Mechanical and Aerospace
Engineering, Princeton University, Princeton, NJ 08544 Princeton Materials
Institute, Princeton University, Princeton, NJ 08544 Sujit S. Datta Department
of Chemical and Biological Engineering, Princeton University, Princeton, NJ
08544

###### Abstract

Obstructions influence the growth and expansion of bodies in a wide range of
settings---but isolating and understanding their impact can be difficult in
complex environments. Here, we study obstructed growth/expansion in a model
system accessible to experiments, simulations, and theory: hydrogels swelling
around fixed cylindrical obstacles with varying geometries. When the obstacles
are large and widely-spaced, hydrogels swell around them and remain intact. In
contrast, our experiments reveal that when the obstacles are narrow and
closely-spaced, hydrogels unexpectedly fracture as they swell. We use finite
element simulations to map the magnitude and spatial distribution of stresses
that build up during swelling at equilibrium, providing a route toward
predicting when this phenomenon of self-fracturing is likely to arise.
Applying lessons from indentation theory, poroelasticity, and nonlinear
continuum mechanics, we also develop a theoretical framework for understanding
how the maximum principal tensile and compressive stresses that develop during
swelling are controlled by obstacle geometry and material parameters. These
results thus help to shed light on the mechanical principles underlying
growth/expansion in environments with obstructions.

Many growth and expansion processes are sculpted through confinement by rigid
obstructions. Familiar examples include muffins rising into their
characteristic shape during baking Ding  _et al._ (2019), trees growing around
boulders Taylor  _et al._ (2021), and even cities expanding around
inhospitable geographic features Borsdorf and Haller (2020). Obstructed growth
and expansion also play pivotal roles---both harmful and beneficial---in many
practical applications. For example, excessive tissue growth around metal mesh
tubes inserted into blood vessels is a common, but life-threatening,
complication of stenting Byrne  _et al._ (2015); Kuhl  _et al._ (2007); Cheng
and Zhang (2021); conversely, the expansion of spray foam into cracks and in
between walls underlies the thermal insulation of many energy-efficient
buildings Cabeza  _et al._ (2010). More broadly, obstructed growth and
expansion critically influence the emergence of form and function across
diverse non-living and living systems, ranging from hydrogels added to soil
for water retention to biofilms and biological tissues in complex environments
Louf  _et al._ (2021); Beebe  _et al._ (2000); Alben (2022); Chu  _et al._
(2018); Streichan  _et al._ (2014); Bengough  _et al._ (1997); Alric  _et al._
(2022). Therefore, we ask: are there general principles that dictate how
obstructions influence growth and expansion? And if so, how do we discover
them?

Figure 1: Hydrogels swelling around obstacles remain intact at equilibrium
when obstacles are further apart, but fracture when obstacles are closer
together. (A) Schematic of the experimental platform, showing a hydrogel
(yellow) surrounded by cylindrical obstacles with radius $r_{\rm{obs}}$
separated according to $r_{\rm{ctr}}$ and confined vertically by parallel
plates separated by $\Delta z$. (B-C) Top view images taken over the course of
swelling, with $r_{\rm{obs}}=5~{}\rm{mm}$. The approximate borders of the
obstacles and hydrogels are outlined for clarity. As each hydrogel swells, the
yellow dye fixed in its polymer network becomes more dilute; thus, the color
intensity serves as a proxy for the local polymer concentration. (B) When the
obstacles are further apart ($r_{\rm{ctr}}=20~{}\rm{mm}$), a hydrogel reaches
a four-lobed clover-like shape at equilibrium. (C) When the obstacles are
closer together ($r_{\rm{ctr}}=10~{}\rm{mm}$), cracks appear at the surface of
the hydrogel as it swells, driving the repeated production of discrete
fragments.

When the growth/expansion of a body is resisted by surrounding obstructions,
large and spatially-nonuniform stresses can develop, influencing subsequent
growth/expansion in turn Ambrosi  _et al._ (2019); Goriely (2017). Being able
to understand the distribution and magnitude of these stresses is thus a
necessary step in the development of widely-applicable, predictive models of
growth and expansion. Unfortunately, due to the complexity of real-world
materials and environments, model systems in which the coupling between growth
and stress can be systematically studied are scarce. Here, we address this
issue using studies of spherical hydrogel beads swelling in 3D-printed
obstacle arrays with defined geometries. Hydrogels are cross-linked networks
of hydrophilic polymers that can absorb large amounts of water and swell while
still retaining integrity. As a result, they are extensively used in
biomedical, environmental, and manufacturing applications, and have well-
characterized and highly-tunable properties such as their degree of swelling
and elasticity Zhang and Khademhosseini (2017); Oyen ; Deligkaris  _et al._
(2010); Guo  _et al._ (2020). Indeed, the comprehensive theoretical literature
on hydrogel swelling makes computations of shape and internal stresses
accessible for a variety of boundary conditions. While much work has focused
on the case of a hydrogel swelling while adhered to another material Tanaka
_et al._ (1987); Hong  _et al._ (2008, 2009); Kang and Huang (2010); Marcombe
_et al._ (2010); Amar and Ciarletta (2010); Dervaux and Ben Amar (2012); Kim
_et al._ (2010); Cho  _et al._ (2019), non-adhered swelling around
obstructions has received limited attention, primarily being considered in the
context of indentation testing Hui  _et al._ (2006); Hong  _et al._ (2008); Hu
_et al._ (2010); Bouklas  _et al._ (2015).

Our study extends this body of work. First, we use experiments to directly
visualize how hydrogel swelling is altered by obstacles of systematically-
varying geometries. When obstacles are positioned further apart, the hydrogel
swells through the spaces between them and maintains its integrity. By
contrast, when the obstacles are closer together, we observe a surprising
phenomenon: the hydrogel fractures, repeatedly tearing itself apart as it
swells! We then use theory and simulations to rationalize these observations,
and importantly, quantify and understand the distribution of stresses. Taken
together, our work provides a prototypical example of obstructed
growth/expansion and uncovers unexpected swelling behaviors and mechanical
instabilities that can result during this process---highlighting the rich
physics waiting to be explored in this area of soft mechanics.

##  I Experiments

Our experimental platform is schematized in Fig. 1A and detailed in Materials
and Methods. To define the obstacles, we 3D-print four rigid cylindrical
columns of radius $r_{\rm{obs}}$ to be placed an equal distance $r_{\rm{ctr}}$
from a central point; in the experiments, we vary $r_{\rm{obs}}$ and
$r_{\rm{ctr}}$ between $2.5-15$ mm and $3.5-25$ mm, respectively. The
cylinders are securely attached to horizontal, parallel laser-engraved acrylic
plates spaced vertically by a fixed amount, $\Delta z=40$ mm. Importantly,
these plates are transparent, permitting direct visualization of a hydrogel as
it swells between the cylindrical obstacles and parallel plates. Hence, at the
beginning of each experiment, we place a spherical polyacrylamide hydrogel
bead of initial radius $\sim 6$ mm (characterized in SI Appendices A, B) in
the center and submerge the entire apparatus in a bath of ultrapure milli-Q
water---thereby initiating swelling, which we track using a camera focused on
the top plate.

We first examine the case of obstacles that are spaced further apart. As the
hydrogel swells, it contacts the top and bottom plates, as well as the
cylinders, and continues to swell through the space between them (Fig. 1B,
Movie S1). It eventually reaches an unchanging four-lobed equilibrium shape
that reflects the balance between the osmotic pressure driving swelling and
the elastic stresses arising from the polymer network and hydrogel-obstacle
contact Louf and Datta (2021); Louf  _et al._ (2021).

The case of closer-spaced obstacles of the same size is dramatically
different. We observe similar behavior to the previous case of less
confinement at early times: the hydrogel contacts the surrounding surfaces and
swells through the space between them. As these lobes continue to swell,
however, cracks abruptly form at the hydrogel surface (48 h in Fig. 1C),
reflecting the development of large stresses during obstructed swelling.
Remarkably, this process then continues, resulting in elaborate, multi-step
fracturing of the hydrogel as it swells, repeatedly ejecting fragments of the
hydrogel outward (Fig. 1C, Movie S2). This process eventually stops as the
hydrogel reaches its final equilibrium degree of swelling. Several other
representative examples of fracturing hydrogels with different obstacle
geometries are shown in Movies S3 and S4. The fracturing process varies
significantly between samples, reflecting the acute sensitivity of crack
formation and propagation on the random imperfections in the hydrogel and the
complex topography arising from previous fracturing Gdoutos (2020); Wang  _et
al._ (2022); Li  _et al._ (2023). To our knowledge, this is the first report
of fracture induced by obstructed swelling.

Figure 2: Experiments reveal a hydrogel fracture threshold that depends on
obstacle radius and spacing. Each $\square$ indicates an experiment in which
the hydrogel reached an intact equilibrium shape as in Fig. 1B, while each
$\times$ indicates an experiment that resulted in fracture as in Fig. 1C.
Schematics show a top view of the obstacle geometries for the indicated
points. The grey excluded region in the top left shows parameters for which
the obstacles would overlap. The green background shading provides a guide to
the eye.

These observations suggest that, when sufficiently obstructed, a
growing/expanding body can tear itself apart. To characterize the dependence
of this phenomenon on confinement, we repeat these experiments with obstacles
of varying sizes $r_{\rm{obs}}$ and spacings $r_{\rm{ctr}}$. We observe a
similar phenomenon in all cases, as summarized in the state diagram in Fig. 2.
When the obstacles are far apart, the hydrogel swells and retains its
integrity ($\square$), while when the obstacles are closer i.e.,
$r_{\rm{ctr}}$ is below a threshold value, the hydrogel repeatedly self-
fractures as it swells ($\times$). This threshold value also depends on the
size of the obstacles. One might expect that fracturing is promoted when the
hydrogel has less free space to swell into i.e., larger $r_{\rm{obs}}$; in
stark contrast to this expectation, we find that for a given $r_{\rm{ctr}}$,
fracturing occurs below a threshold value of $r_{\rm{obs}}$. We rationalize
this intriguing observation using theory and simulations in the following
sections. The two geometric parameters $r_{\rm{ctr}}$ and $r_{\rm{obs}}$ thus
delineate a boundary between swelling without fracture and obstructed swelling
causing fracture, as shown by the light and dark green regions, respectively,
in Fig. 2.

To demonstrate the generality of this phenomenon, we repeat our experiments in
fully 3D granular packings, which more closely mimic the obstructions
experienced by hydrogels in many applications Misiewicz  _et al._ (2022);
Krafcik  _et al._ (2017); Kapur  _et al._ (1996); Dolega  _et al._ (2017); Lee
_et al._ (2019); Woodhouse and Johnson (1991); Demitri  _et al._ (2013);
Guilherme  _et al._ (2015); Lejcuś  _et al._ (2018); Wei and Durian (2013).
Our previous experiments Louf  _et al._ (2021) using this platform used a
solvent that promoted only slight hydrogel swelling, and therefore accessed
the case of weak confinement; however, repeating these experiments with a
solvent that promotes more hydrogel swelling indeed gives rise to self-
fracturing, as shown in Fig. S4 (SI Appendix C). Thus, this phenomenon of
obstructed growth/expansion causing fracture arises not just in idealized
geometries, but also in more realistic complex spaces.

##  II Theory and simulations

What are the essential physical principles that govern this phenomenon? To
address this question, we develop finite element simulations that incorporate
the energetic penalty of contacting obstructions into a model of hydrogel
swelling based on classic Flory-Rehner theory. Our goal is to better
understand the complex distribution of stresses that arises during obstructed
swelling, not to quantitatively capture all the features of the experiments;
as such, we use a simplified two-dimensional (2D) model that permits
straightforward visualization of stresses and thereby enables us to develop an
intuitive and analytic description of the underlying physics.

We describe stretching of the hydrogel polymer network and solvent-polymer
interactions with the commonly-used free energy density Flory and Rehner
(1943); Flory (1953); Hong  _et al._ (2008, 2009); Kang and Huang (2010);
Bouklas  _et al._ (2015):

$\displaystyle\frac{F_{\rm{en}}}{A_{0}k_{B}T}$
$\displaystyle=\frac{n_{p}}{2}\left(F_{iK}F_{iK}-2-2\ln(\det(\mathbf{F}))\right)$
$\displaystyle+\frac{1}{\Omega}\left(\Omega C\ln\left(\frac{\Omega C}{1+\Omega
C}\right)+\chi\frac{\Omega C}{1+\Omega C}\right).$ (1)

The first term describes the elastic free energy: $n_{p}$ is the number of
polymer chains per unit dry reference area $A_{0}$ and $F_{iK}=\frac{\partial
x_{i}}{\partial X_{K}}$ is the deformation gradient tensor, with
$F_{iK}F_{iK}=\mathrm{tr}(\mathbf{F}^{T}\mathbf{F})=\sum_{i}\lambda_{i}^{2}$
in terms of principal stretches $\lambda_{i}$. Following standard conventions,
deformed/current configurations are denoted by lowercase letters and
reference/initial configurations are denoted by capital letters unless
otherwise noted Goriely (2017). The second term describes the mixing free
energy: $\Omega$ is the area of a solvent molecule, $C$ is the nominal
concentration of solvent (number of solvent molecules per unit reference
area), and $\chi$ is the Flory-Huggins interaction parameter. Both terms are
scaled by the product of the Boltzmann constant and temperature, $k_{B}T$.

Given that the hydrogel network is held together by permanent cross-links
between its polymer chains, we additionally assume that it only changes volume
by uptake of solvent, which allows us to express concentrations in terms of
the deformation: $\det(\mathbf{F})=1+\Omega C$. We require the chemical
potential $\mu$ to be zero on the boundary of the hydrogel to mimic submerging
it in pure water, as in the experiments. To impose this boundary condition, we
perform the standard Legendre transform of Eq. (1) and derive a new free
energy $\hat{F}_{\rm{en}}$ in terms of $\mu$ as
$\hat{F}_{\rm{en}}(x_{i},\mu)=F_{\rm{en}}(x_{i},C)-\mu C$ Hong  _et al._
(2009). Finally, we model contact by imposing an energy penalty for overlap
between the hydrogel and obstacles; as detailed in Materials and Methods, our
results are insensitive to the choice of this parameter. Note that the
chemical potential boundary condition $\mu=0$ is enforced in the regions of
contact between the hydrogel and the obstacles (see SI Appendix E for further
discussion).

Figure 3: Finite element simulations quantify how the equilibrium strains and
stresses that develop in the hydrogel depend on obstacle geometry. (A-D) Maps
of the solvent concentration per unit reference area with the finite element
grid superimposed (left column) and the principal stresses (right column) in
the hydrogel at equilibrium, normalized by the area of the solvent molecule
$\Omega$ and the fully-swollen hydrogel shear modulus $G_{0}$, respectively.
Confinement increases from top to bottom, in this case by changing the
obstacle spacing with fixed $r_{\rm{robs}}/r^{*}=0.3$; thus,
$\Delta/r^{*}=1-r_{\rm{ctr}}/r^{*}=0.02,0.1,0.5,0.6$ from top to bottom, where
$\Delta$ is the difference between the unconfined swollen hydrogel radius
$r^{*}$ and the distance to an obstacle from the center $r_{\rm{ctr}}$. Each
mesh point in the right column bears two perpendicular lines, oriented and
scaled according to the direction and magnitude of the principal stresses
$\sigma_{1}$, $\sigma_{2}$ at that point, and colored such that compressive
stresses are blue and tensile stresses are red. The stresses exceed the range
of the color bar for (C-D). (E-F) Points show the maximum principal tensile
and compressive stresses obtained from the simulations, taken along the line
connecting the hydrogel center to an obstacle center, at this same fixed
$r_{\rm{obs}}/r^{*}=0.3$ as a function of $\Delta/r^{*}$. Stresses are again
normalized by $G_{0}$. The values of $\Delta/r^{*}$ corresponding to (A--D)
are marked by the grey diamonds. The predictions of linear indentation theory
(solid lines) agree well with the simulation data for small $\Delta/r^{*}$
(purple shaded region). With increasing $\Delta/r^{*}$, the maximum tension
exceeds linear predictions, but can be reproduced by a _geometrically_
nonlinear elastic theory with a linear constitutive law (blue shaded region,
SI Appendix H). For compression, the linear theory is accurate over a larger
range of $\Delta/r^{*}$, but the geometrically nonlinear theory cannot explain
the deviations; instead, an elastic model incorporating _material_
nonlinearities (Eq. (10), dashed line) better captures the data. As
$\Delta/r^{*}$ increases even further, the hydrogel exhibits a symmetry-
breaking instability (D, red shaded region).

To visualize and track hydrogel deformation during obstructed swelling, along
with the concomitant development of internal stresses, we implement this model
in FEniCS Alnæs  _et al._ (2015), an open-source finite element platform. Our
simulations consider a circular hydrogel swelling around four fixed obstacles.
Although our primary focus is the stress distribution in the hydrogel at
equilibrium, we simulate the full dynamics of obstructed swelling to provide
numerical stability and ensure that there is a realistic path from the initial
to the final equilibrium swollen state (SI Appendix E).

We begin by examining hydrogel swelling around obstacles that have the same
size, but varying spacing---just as in the experiments shown in Fig. 1(B,C).
Four representative examples, varying from the case of weak to strong
confinement, are shown in Fig. 3(A-D). The corresponding maximum tensile and
compressive principal stresses on the line connecting the hydrogel center to
an obstacle center are plotted in panels (E,F). We quantify these variations
using $\Delta/r^{*}\equiv\left(r^{*}-r_{\rm{ctr}}\right)/r^{*}$, where $r^{*}$
is the equilibrium, fully-swollen, unconfined radius of the hydrogel and
$\Delta$ is the difference between $r^{*}$ and the distance from the center to
the closest obstacle $r_{\rm{ctr}}$. The case of weak confinement in (A) has
$\Delta/r^{*}=0.02$, increasing to $\Delta/r^{*}=0.6$ for the case of strong
confinement in (D). As seen by the color maps of solvent concentration and
principal stresses in panel (A), in this limit of weak confinement, the
hydrogel is barely deformed and the resultant stresses are not visible when
plotted on the same scale as (B--D). Thus, as we describe in Sec. II.1, the
hydrogel stresses can be captured analytically using linear theories in this
regime. Increasing hydrogel confinement causes the magnitudes of the stresses
to increase (panels E,F), but they still remain primarily localized near the
obstacles; as discussed in Sec. II.2 and Sec. II.3, the largest compressive
stresses can still be described by linear theory, but the largest tensile
stresses require consideration of a geometric nonlinearity. As the separation
between obstacles decreases further, large compressive stresses span the
entirety of the hydrogel (panel C), eventually triggering a symmetry-breaking
instability (panel D) as discussed in Sec. II.4. Finally, while our model is
not suitable to treat fracture directly, we discuss the connection between our
calculations of stresses and the experimental observations of swelling-induced
self-fracture in Sec. III.

###  II.1 Weak confinement

Consider a hydrogel disk that has swollen around obstacles to reach
equilibrium. The chemical potential is spatially uniform and therefore all
solvent transport has stopped. Nonetheless, due to contact with the obstacles,
the distribution of solvent is inhomogeneous through the hydrogel: solvent
preferentially enters the uncompressed lobes of the hydrogel between the
obstacles, as shown by the yellow color in Fig. 1B Louf and Datta (2021). Now,
consider another hydrogel that was first swollen, unobstructed, to its
equilibrium size, and then slowly and incrementally squeezed by an identical
set of obstacles moving towards its center, acting as four indenters. Solvent
must exit the hydrogel where it is indented by the obstacles; recall our
condition $1+\Omega C=\det(\mathbf{F})$. For the same final obstacle geometry,
these two scenarios must have identical solvent distributions and stress
profiles at equilibrium. Thus, the long time limit of obstructed swelling can
be treated as an indentation problem, which is well-studied in the limit of
small deformations. Making this analogy between obstructed swelling and
indentation allows us to apply lessons from a large body of literature on
linear contact mechanics and poroelasticity to derive expressions for the
stress tensor in the hydrogel Hui  _et al._ (2006); Hu  _et al._ (2010);
Goriely  _et al._ (2016); Weickenmeier  _et al._ (2016).

Assuming that deformations relative to the fully swollen state are small, we
linearize Eq. (1) to find effective linear elastic parameters (SI Appendix F,
Hu  _et al._ (2011); Bouklas and Huang (2012)). In particular, for a 2D
hydrogel, the effective Poisson’s ratio $\nu$ and Young’s modulus $E$ are
given by

$\displaystyle\nu$
$\displaystyle=1-2n_{p}\Omega\left(n_{p}\Omega\left(1+\frac{1}{\lambda_{0}^{2}}\right)+\frac{1}{\lambda_{0}^{2}(\lambda_{0}^{2}-1)}-\frac{2\chi}{\lambda_{0}^{4}}\right)^{-1},$
(2) $\displaystyle E$ $\displaystyle=2(1+\nu)n_{p}k_{B}T,$ (3)

where $\lambda_{0}$ is the principal stretch corresponding to the fully
swollen state. The expression for the Poisson’s ratio can be understood as its
value for an incompressible 2D solid, $\nu=1$, minus a correction---reflecting
the fact that the compressibility of the swollen hydrogel is generated via
solvent transport. This linearization also yields the shear modulus,
$G_{0}=n_{p}k_{B}T$, which we use to normalize stresses throughout this paper.

Given these effective elastic parameters, we solve for the stresses in the
hydrogel as a 2D linear contact mechanics problem (SI Appendix G, Barber
(2018); Johnson (1987)). This approach provides expressions for the stress
tensor $\sigma_{ij}$ along the line directly beneath the top obstacle as a
function of $y$ as shown in Fig. 4B:

$\displaystyle\sigma_{xx}$
$\displaystyle=\frac{2\zeta}{\pi}\Bigg{(}\frac{2r^{*3}}{(r^{*2}+y^{2})^{2}}-\frac{1}{r^{*}}-\frac{2(r^{*}-y)}{a^{2}}$
$\displaystyle+\frac{2(r^{*}-y)^{2}+a^{2}}{a^{2}\sqrt{(r^{*}-y)^{2}+a^{2}}}\Bigg{)},$
(4) $\displaystyle\sigma_{yy}$
$\displaystyle=\frac{2\zeta}{\pi}\left(\frac{1}{r^{*}+y}-\frac{1}{r^{*}}+\frac{2r^{*}y^{2}}{(r^{*2}+y^{2})^{2}}+\frac{1}{\sqrt{(r^{*}-y)^{2}+a^{2}}}\right),$
(5) $\displaystyle\sigma_{xy}$ $\displaystyle=0,$ (6)

where $\zeta<0$ is the force applied to the indenters and the half contact
width $a$ is

$a=\sqrt{-\frac{4\zeta}{E\pi\left(\frac{1}{r_{\rm{obs}}}+\frac{1}{r^{*}}\right)}},$
(7)

as defined in Fig. 4B. Note that $y$ is defined with respect to the hydrogel’s
unobstructed, fully-swollen state and ranges from $-r^{*}$ to $r^{*}$.

Figure 4: When a hydrogel is weakly confined, linear indentation theory can be
used to predict the stresses, but misses key features in stronger confinement.
(A, B) Stress components $\sigma_{xx}$ and $\sigma_{yy}$, respectively,
determined from simulations (points) and linear theory (lines) along a line
running from the center of the hydrogel ($x=y=0$) to the point of contact with
the top obstacle ($x=0,~{}y=r^{*}$) with fixed $r_{\rm{obs}}/r^{*}=0.3$. The
inset to (B) defines the variables used: the hydrogel-obstacle contact width
is $2a$ and the indenter displacement is $\Delta$, defined relative to the
undeformed hydrogel radius $r^{*}$ (dashed orange line). As the indentation
depth $\Delta/r^{*}\equiv 1-r_{\rm{ctr}}/r^{*}$ increases, the discrepancy
between the linear theory and simulations increases, as expected. With
increasing confinement, tension ($\sigma_{xx}>0$) builds up at
$r_{\rm{ctr}}/r^{*}\approx 0.75$, as shown by the magnified inset in A.

Given these analytical expressions, we next ask: Under what confinement
regimes (as parameterized by $\Delta/r^{*}$) does this linearized theory work,
and when does it break down? And since we would like to gain an intuitive
understanding of obstructed swelling-induced fracture, how does the maximal
principal tensile stress --- which can be used to approximate material
strength Gdoutos (2020)\--- vary with confinement? To address these questions,
we first re-cast Eqs. (4)--(5) in terms of $\Delta/r^{*}$ to facilitate direct
comparison with the results of the nonlinear simulations. To do so, we
integrate the strain $u_{yy}=\frac{\sigma_{yy}}{Y}-\frac{\nu}{Y}\sigma_{xx}$
using Eqs. (4)--(5) to find the displacement at the surface of the indenter
relative to the center of the hydrogel, and expand the result to linear order
in $a/r^{*}$. This procedure gives

$\Delta=-\frac{\zeta}{E\pi}\left(\ln\left(\frac{16r^{*2}}{a^{2}}\right)+\frac{1}{2}(\pi-6-\pi\nu)\right),$
(8)

which we then invert, apply Eqs. (2)--(3), and substitute the resulting
expression for $\zeta(\Delta)$ into Eqs. (4)--(5). The resulting expressions
for $\sigma_{xx}(\Delta)$ and $\sigma_{yy}(\Delta)$ cannot be expressed in
terms of elementary functions, so we omit them, but are shown by the solid
lines in Fig. 4 for the illustrative case of $r_{\rm{obs}}/r^{*}=0.3$; for
comparison, the symbols show the results of the full nonlinear simulations.

As expected, when $\Delta/r^{*}\ll 1$, the linearized indentation solution
agrees well with the nonlinear simulation results, while as $\Delta/r^{*}$
increases, the discrepancy between the two becomes more apparent. In
particular, as exemplified by the data for $\Delta/r^{*}=0.16$ (dark blue
circles) and $\Delta/r^{*}=0.13$ (yellow squares) in Fig. 4, the linear
solution underestimates the compression at both the center $y=0$ and boundary
$y=r^{*}$ of the hydrogel, as well as the tension that builds up at
$y/r^{*}\approx 0.75$ (Fig. 4A, inset). Indeed, though tension
($\sigma_{xx}>0$) does appear in the linear solutions for small
$\Delta/r^{*}$, it disappears with increasing $\Delta/r^{*}$, in contrast to
the simulation results (see arrow in Fig. 3C, for example).

Since we are interested in fracture behavior, we thus focus our attention on
the maximal value of this tensile stress for the same illustrative case of
$r_{\rm{obs}}/r^{*}=0.3$. By symmetry, we expect the largest stresses to lie
beneath each obstacle, since stresses must go to zero at the hydrogel boundary
away from obstacles in the weak confinement regime. Moreover, because the
straight edges of the finite element mesh can introduce spurious tensile
forces at the edge of the hydrogel-obstacle contact (described further in SI
Appendix I), we plot the maximum and minimum values of the principal stresses
along the $x=0$ line shown in the schematic inset of Fig. 4B. The results are
displayed in Fig. 3(E,F). As noted in Fig. 4A, the linear theory only predicts
the presence of tension for small confinement before deviating from the
nonlinear simulation results at $\Delta/r^{*}\gtrsim 0.02$, as shown by the
solid line and points in Fig. 3E, respectively. Interestingly, however, linear
theory captures the maximum compressive stress over a broader range of
confinement, shown by the solid line in Fig. 3F, which agrees well with the
simulation data up to $\Delta/r^{*}\approx 0.2$.

Thus, while linear indentation theory can predict both tensile and compressive
stresses during obstructed swelling in weak confinement, it underpredicts both
for larger deformations---suggesting that the assumptions made in the linear
theory are no longer valid. We revisit these assumptions for both tension and
compression in the next two sections, respectively.

###  II.2 Tension beyond the linear regime

The linear theory in the previous section relies on a number of assumptions
that can fail as deformations increase:

  1. 1.

The effective elastic parameters [Eqs. (2)--(3)] are independent of strain,

  2. 2.

The stress is linearly related to the strain,

  3. 3.

The strain tensor is linear in the displacements.

To assess the validity of these assumptions, we compare the results of our
full nonlinear simulations to those of more complex elastic models that
incorporate _material_ /_geometric_ nonlinearities.

First, we explore the limits of assumption 1 by relaxing assumptions 2 & 3\.
Specifically, we compare the hydrogel simulation results to those of a
compressible neo-Hookean elastic material with elastic parameters given by
Eqs. (2)--(3). As detailed in SI Appendix F, the neo-Hookean model closely
reproduces the stress profiles of the hydrogel simulations over a broad range
of $\Delta/r^{*}$ up to $\approx 0.4$, well beyond the limits of the linear
theory at $\sim 0.02$. Therefore, nonlinearities due to the effective elastic
parameter mapping can be neglected up to this point.

Next, we explore the limits of assumption 2 by relaxing assumption 3.
Specifically, we use a St. Venant-Kirchoff elastic model with strain tensor
$u_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial
u_{j}}{\partial x_{i}}+\frac{\partial u_{k}}{\partial x_{i}}\frac{\partial
u_{k}}{\partial x_{j}}\right)$; thus, the strain tensor is nonlinear in
displacements, but we still require that the hydrogel material follows a
linear constitutive law. Note that derivatives are taken with respect to
coordinates in the unobstructed, fully-swollen state, denoted here with
lowercase letters for simplicity of presentation. Intriguingly, as detailed in
SI Appendix H, the St. Venant-Kirchhoff model quantitatively reproduces the
maximum principal tensile stress in the hydrogel simulations up to
$\Delta/r^{*}\approx 0.1$, well beyond the limit of the linear theory at
$\approx 0.02$. Hence, the excess tension $\sigma_{xx}$ that develops beneath
the obstacles just beyond the linear regime is driven by _geometric_
nonlinearity, and does not require a nonlinear constitutive relationship. At
even larger displacements $\Delta/r^{*}>0.1$, the St. Venant-Kirchoff
simulations are unstable and we expect that both geometric and material
nonlinearities contribute to the tensile stress.

How exactly does geometric nonlinearity generate tension during obstructed
swelling? We answer this question using an illustrative argument reminiscent
of the derivation of the Föppl-von Kármán equation Landau  _et al._ (1986). As
detailed in _Materials and Methods_ , we first find the variation of the
integrated St. Venant-Kirchhoff strain-energy function with respect to
displacements, which can be written in terms of the second Piola-Kirchhoff
stress tensor, $\sigma_{ij}^{PK}$. In 2D, this quantity gives the stress
component in a material direction $i$ perpendicular to a line that has unit
length and normal $j$ in the reference configuration Audoly and Pomeau (2010).
We then make approximations specific to our obstacle geometry. Ultimately, we
find

$\sigma_{xx}^{PK}\approx-\frac{\partial\sigma_{yy}^{PK}}{\partial
y}\left(1+\frac{\Delta}{r^{*}}\right)r_{\rm{obs}}.$ (9)

Since $\sigma_{yy}^{PK}$ becomes more negative as $y$ increases,
$\frac{\partial\sigma_{yy}^{PK}}{\partial y}<0$ (Fig. 4B), and therefore
$\sigma_{xx}^{PK}>0$ indicating tension. Thus, geometric nonlinearity
generates tension beneath an indenter perpendicular to the indentation
direction, qualitatively matching our simulations.

###  II.3 Compression beyond the linear regime

A notable result shown in Fig. 3 is that while linear indentation theory
predicts tensile stresses for $\Delta/r^{*}<0.02$, it captures the compressive
stresses over a broader range, up to $\Delta/r^{*}\approx 0.2$. Which
nonlinearities drive the deviations that arise at even larger displacements?
We answer this question by following the same procedure as in the previous
section, detailed further in SI Appendix F & H. In contrast to the case of
tension, we do not find any parameters for which the St. Venant-Kirchoff model
is more accurate than the linear model. Furthermore, as shown in Fig. 3F, past
the linear regime, the compressive (Cauchy) stress underneath the top obstacle
scales like that of a compressible neo-Hookean elastic material experiencing
uniaxial compression, with the principal stretch parallel to indentation set
to $\lambda_{1}=1-\Delta/r^{*}$ and the principal stretch perpendicular to the
indentation set to 1:

$\sigma_{yy}\sim
G_{0}\left(1-\frac{\Delta}{r^{*}}\right)+\frac{2G_{0}\nu}{1-\nu}\frac{\ln\left(1-\frac{\Delta}{r^{*}}\right)}{\left(1-\frac{\Delta}{r^{*}}\right)}-\frac{G_{0}}{\left(1-\frac{\Delta}{r^{*}}\right)}.$
(10)

Thus, unlike the case of tension, deviations from the linear theory do not
arise from geometric nonlinearities and can instead be attributed to
_material_ nonlinearities.

###  II.4 Symmetry-breaking instability

A striking phenomenon arises in our simulations as the separation between
obstacles decreases further: as shown in Fig. 3D, the hydrogel displays a
symmetry-breaking instability and swells preferentially along a diagonal. Why
does this instability arise? Inspecting the spatial distribution of
compressive stresses provides a clue. As the hydrogel swells in increasing
amounts of confinement, its central core becomes increasingly compressed (see,
e.g., Fig. 3C--D). Compressing this circular core along a _single_ axis,
forming an ellipse, requires less energy than does compressing this core
isotropically. Thus, as confinement increases, one expects the case of
asymmetric swelling to be energetically preferred, leading to this instability
---as described further in SI Appendix J.

Figure 5: Simulations predict that hydrogel stresses grow as obstacles are
made smaller and brought closer together in a manner similar to the
experimental findings. (A, B) Most positive (tensile) principal stress
$\sigma_{\rm{max}}$ and magnitude of the most negative (compressive) principal
stress $|\sigma_{\rm{min}}|$, respectively, again normalized by the fully-
swollen shear modulus $G_{0}$. Each (light or dark) green square corresponds
to a simulation. As in Fig. 2, the grey region corresponds to overlapping
obstacles. Note that we take the maxima/minima over the entire mesh, rather
than just beneath an obstacle, to avoid non-monotonic behavior in the
instability regime.

##  III Comparison between simulations and experiments

Our theoretical analyses and simulations capture the essential features of the
hydrogel deformation observed in experiments (e.g., Figs. 1B and 3C). They
also enabled us to explore how stresses develop during obstructed swelling
more generally (Figs. 3\--4). While our model is not suitable to directly
treat the swelling-induced self-fracture observed experimentally, the
simulated stresses help rationalize this phenomenon. To this end, we compare
the simulations to the experimental state diagram shown in Fig. 2 by plotting
the maximum principal tensile and compressive stresses as a function of
$r_{\rm{obs}}$ and $r_{\rm{ctr}}$. The results, shown in Fig. 5, bear a
compelling resemblance to the experimental results. In particular, the
convexity and shape of the experimental fracture boundary are similar to the
simulated contours of maximum principal stress. Indeed, appealing to a
commonly-used fracture criterion for brittle materials Gdoutos (2020), we
expect that hydrogel fracture occurs when the maximum tensile stress exceeds a
threshold --- whose exact value would establish the position of the
experimental fracture boundary in Fig. 5A. We conjecture that this threshold
is reached prior to the symmetry-breaking instability revealed by the
simulations, as we do not observe it in our experiments; experiments using
tougher hydrogels than those studied here may be a useful way to probe this
deformation mode in future work.

##  IV Discussion

Despite the ubiquity of obstructed growth and expansion in our everyday lives,
how exactly this process generates large, spatially non-uniform stresses in a
body --- and how these stresses influence its subsequent growth/expansion in
turn --- has remained challenging to systematically study. One reason is the
lack of model systems in which this intricate coupling between growth and
stress can be probed both experimentally and theoretically. We addressed this
need by studying the swelling of spherical hydrogel beads in obstacle arrays
of tunable geometries. Our experiments revealed a striking phenomenon: under
weak confinement, a hydrogel retains its integrity and assumes a symmetric,
four-lobed clover-like shape, while in stronger confinement, it repeatedly
tears itself apart as it swells. We elucidated the underlying physics by
adopting established models of hydrogel swelling to map the tensile and
compressive stresses arising during swelling. In particular, we found that
when a hydrogel is weakly deformed, stresses are well-described by linear
indentation theory, while as a hydrogel is increasingly deformed, geometric
and material nonlinearities engender large tensile and compressive stresses
tangential and normal to the obstacles, respectively, driving the hydrogel
towards fracture.

Because our study represents a first step toward fully unravelling the
mechanics of obstructed growth and expansion, it necessarily has limitations.
For example, while the experimental results shown in Fig. 2 revealed a
fracture threshold that varies with obstacle geometry, quantitative comparison
to theory and simulation will require more precise control over the system
geometry and dimensionality Joshi  _et al._ (2022, 2023), both in experiments
and in the model, along with a more detailed treatment of the microscopic
processes underlying fracturing Lai and Hu (2018); Kuna (2013); Mao and Anand
(2018); Wang and Hong (2012); Mao and Anand (2018); Baumberger and Ronsin
(2020); Yang  _et al._ (2022). Moreover, many more experimental trials near
this threshold and with hydrogels of varying mechanical properties, along with
high-resolution imaging of crack propagation Leslie  _et al._ (2021), will be
useful in characterizing the details of the fracturing process, which likely
depend on the presence of microscopic imperfections in the hydrogel. Finally,
although here we restricted our attention to the case of rigid obstacles, many
scenarios involve growth/expansion around _compliant_ constraints---e.g., the
development of biofilms, tissues, and organs in the body Helmlinger  _et al._
(1997); Zhang  _et al._ (2021); Fortune  _et al._ (2022); Goodwin  _et al._
(2019), with potential implications for biological function Möbius  _et al._
(2015); Atis  _et al._ (2019). Extending our study to the case of deformable
obstacles would therefore be a useful direction for future work.

Our results may be especially relevant to diverse applications of hydrogels,
and other swellable materials, that frequently involve their confinement in
tight and tortuous spaces. For example, driven by growing demands for food and
water, hydrogels are increasingly being explored as additives to soils to
absorb and release water to plants and therefore reduce the burden for
irrigation Misiewicz  _et al._ (2022); Woodhouse and Johnson (1991); Demitri
_et al._ (2013); Guilherme  _et al._ (2015); Lejcuś  _et al._ (2018). They are
also widely adopted in other applications, such as oil recovery, construction,
mechanobiology, and filtration, all of which involve hydrogel swelling in
confinement Krafcik  _et al._ (2017); Kapur  _et al._ (1996); Dolega  _et al._
(2017); Lee  _et al._ (2019). A common assumption made in all these cases is
that the hydrogel retains integrity as it swells; however, our study indicates
that these applications should be evaluated for the possibility of swelling-
induced self-fracture. Indeed, fracture could lead to the production and
dispersal of many small hydrogel fragments, potentially reducing their utility
and leading to environmental contamination. This process should therefore be
carefully considered in a wide range of real-world contexts.

## Materials and Methods

#### Experimental details

To create the obstacle array, we 3D print cylindrical columns in Clear v4
resin using a Form3+ industrial 3D printer (Formlabs), and cut the acrylic
plates using an Epilog Laser Mini 24 laser cutter and engraving system. We
secure the columns to the acrylic plates using a twist-and-lock mechanism. The
hydrogels are polyacrylamide beads (‘‘water gel beads" obtained from
Jangostor) and are stored in a screw cap container prior to experiments; as
such, they experience some slight swelling due to ambient humidity. The
hydrogel beads have varying sizes and colors, but all appear to have the same
swelling behavior and beads of similar sizes are used for experiments (SI
Appendix A). These hydrogel beads were extensively characterized in our
previous experiments Louf  _et al._ (2021); Louf and Datta (2021), which
provide additional detail.

Early in the swelling process (e.g. 2 h in Fig. 1B,C), each hydrogel appears
out of focus since it has not yet made contact with the top plate (focal plane
for imaging), and cusps are visible on its surface due to differential
swelling as water enters the hydrogel from its outer surface Tanaka (1986);
Tanaka  _et al._ (1987); Bertrand  _et al._ (2016); Curatolo  _et al._ (2017).
We verify that the hydrogel beads swell to an equilibrium shape without
rupturing when no obstacles are present (either with or without the 40 mm
$z$-confinement), indicating the transient stresses that engender these cusps
are insufficient to drive fracture as in other less tough gels de Silva and
Lapitsky (2016); Leslie  _et al._ (2021). Because the plates and obstacles are
made of acrylic or a polymeric resin, respectively, we do not expect or see
evidence of adhesion between the polyacrylamide hydrogel surface and the
confining surfaces. For experiments in which we image the entire process of
obstructed swelling, we verify that the hydrogel reaches an equilibrium shape
(in the cases of no fracturing) when its size/shape does not noticeably change
for at least three hours. In other experiments where we do not image the
entire process of obstructed swelling, we verify that equilibrium is reached
by waiting an additional 12 hours. Each symbol in Fig. 2 represents a single
experiment.

#### Simulation details

We create meshes using FEniCS’s built-in mesh generation function with cell
size set such that there are 30 vertices along the radius. We confirmed that
this mesh resolution was sufficient for numerical convergence. We set
$n_{p}\Omega=0.001$ and $\chi=0.3$. The penalty function is integrated over
the hydrogel boundary, and is given by

$\frac{F_{\rm{pen}}}{2\pi r_{0}k_{B}T}=\frac{p}{4\pi r_{0}}\sum_{i}\langle
r_{\rm{obs}}^{2}-(\mathbf{x}-\mathbf{x}_{i})^{2}\rangle_{+}^{2},$ (11)

where $p$ is the penalty strength, $\mathbf{x}_{i}$ is position of the center
of the $i^{\rm{th}}$ obstacle, $r_{0}$ is the dry reference radius of the
hydrogel, and the sum is over all obstacles. The brackets
$\langle\cdots\rangle_{+}$ take the positive part of the argument, defined as
$\langle x\rangle_{+}=\frac{x+|x|}{2}$. Thus, when evaluated at positions away
from any obstacles, the penalty function is zero, but takes a large positive
value inside the obstacles. We set $p=6.25\pi/r_{0}^{4}$ to generate the data
shown in this text, and have verified that using $p=62.5\pi/r_{0}^{4}$
produces the same results.

Following the suggestions of Refs. Zhang  _et al._ (2009); Bouklas  _et al._
(2015); Ang  _et al._ (2020), we use the backwards Euler scheme for time
integration (see SI Appendix E for further discussion of dynamics), Taylor-
Hood mixed elements (quadratic elements for the displacement field and linear
elements for the chemical potential field), and early time ramping of the
chemical potential boundary condition to ensure numerical stability. The
Newton-Raphson method is used at each time step, and equilibrium is defined by
when zero iterations are required for a step to complete.

To find the maximum/minimum principal Cauchy stresses, we first calculate the
eigenvalues of the stress tensor for a given displacement field. We project
these eigenvalue fields onto a function space of discontinuous Lagrange
elements of order 1. We then compare the eigenvalues defined on this mesh. The
largest positive eigenvalue is the maximum principal (tensile) stress
$\sigma_{\rm{max}}$, and the most negative eigenvalue is the minimum principal
(compressive) stress, $\sigma_{\rm{min}}$. To find the minimum and maximum
stresses beneath an obstacle (data in Fig. 3(E,F)), we instead project the
stress tensor onto the vertical line directly beneath the top obstacle. The
minimum value of $\sigma_{yy}$ and the maximum value of $\sigma_{xx}$ along
this line are plotted as the below top obstacle minimum and maximum stresses
respectively.

#### Tension generated by geometric nonlinearity

We describe the argument leading up to Eq. (9) in more detail, which
demonstrates how tension can appear beneath an obstacle when a nonlinear
strain tensor is used. The variation of the integrated St. Venant-Kirchhoff
strain-energy function with respect to displacements in terms of the second
Piola-Kirchhoff stress tensor is

$\delta W=\int\sigma_{ij}^{PK}\delta
u_{ij}dA=\int\sigma_{ij}^{PK}\left(\frac{\partial\delta u_{i}}{\partial
x_{j}}+\frac{\partial u_{k}}{\partial x_{i}}\frac{\partial\delta
u_{k}}{\partial x_{j}}\right)dA.$ (12)

Upon integrating by parts, assuming we cannot vary the displacements at the
boundaries due to the presence of obstacles, we find

$\displaystyle\delta W=\int\left(\frac{\partial\sigma_{ij}^{PK}}{\partial
x_{j}}\delta u_{i}+\frac{\partial}{\partial
x_{j}}\left(\sigma_{ij}^{PK}\frac{\partial u_{k}}{\partial x_{i}}\right)\delta
u_{k}\right)dA,$
$\displaystyle=\int\left(\frac{\partial\sigma_{ij}^{PK}}{\partial x_{j}}\delta
u_{i}+\frac{\partial\sigma_{ij}^{PK}}{\partial x_{j}}\frac{\partial
u_{k}}{\partial x_{i}}\delta
u_{k}+\sigma_{ij}^{PK}\frac{\partial^{2}u_{k}}{\partial x_{i}x_{j}}\delta
u_{k}\right)dA.$ (13)

Next, in order to make approximations specific to our geometry, we explicitly
list all the terms that appear for a two-dimensional solid. Since we will
examine stresses directly beneath an obstacle, we set $\sigma_{xy}^{PK}=0$ by
symmetry, yielding:

$\displaystyle\delta W$ $\displaystyle\approx\int dA\delta
u_{x}\Bigg{(}\frac{\partial\sigma_{xx}^{PK}}{\partial x}\left(1+\frac{\partial
u_{x}}{\partial x}\right)+\frac{\partial\sigma_{yy}^{PK}}{\partial
y}\frac{\partial u_{x}}{\partial y}$
$\displaystyle+\sigma_{xx}^{PK}\frac{\partial^{2}u_{x}}{\partial
x^{2}}+\sigma_{yy}^{PK}\frac{\partial^{2}u_{x}}{\partial y^{2}}\Bigg{)}$
$\displaystyle+\int dA\delta
u_{y}\Bigg{(}\frac{\partial\sigma_{yy}^{PK}}{\partial y}\left(1+\frac{\partial
u_{y}}{\partial y}\right)+\frac{\partial\sigma_{xx}^{PK}}{\partial
x}\frac{\partial u_{y}}{\partial x}$
$\displaystyle+\sigma_{xx}^{PK}\frac{\partial^{2}u_{y}}{\partial
x^{2}}+\sigma_{yy}^{PK}\frac{\partial^{2}u_{y}}{\partial y^{2}}\Bigg{)}.$ (14)

In equilibrium, the coefficients of $\delta u_{x}$ and $\delta u_{y}$ must be
zero (in the absence of body forces). We focus our attention on the
coefficient of $\delta u_{y}$. If we consider the internal stresses directly
beneath the obstacle, we can set $\partial u_{y}/\partial x$ to zero by
symmetry. We approximate the deformation as an affine contraction in the $y$
direction, which sets $\partial u_{y}/\partial y\approx\Delta/r^{*}$ and
$\partial^{2}u_{y}/\partial y^{2}\approx 0$. We assume that the curvature of
the $y$ displacement, $\frac{\partial^{2}u_{y}}{\partial x^{2}}$, is
determined by the curvature of the obstacle, $1/r_{\rm{obs}}$. With these
substitutions, the equilibrium condition becomes

$\frac{\partial\sigma_{yy}^{PK}}{\partial
y}\left(1+\frac{\Delta}{r^{*}}\right)+\frac{\sigma_{xx}^{PK}}{r_{\rm{obs}}}\approx
0.$ (15)

Solving for $\sigma_{xx}^{PK}$, we find

$\sigma_{xx}^{PK}\approx-\frac{\partial\sigma_{yy}^{PK}}{\partial
y}\left(1+\frac{\Delta}{r^{*}}\right)r_{\rm{obs}}.$ (16)

Just beyond the linear regime, $\sigma_{yy}^{PK}$ should follow the same
trends as $\sigma_{yy}$ predicted by the linear theory. In Fig. 4B, we observe
that $\sigma_{yy}$ decreases as $y$ increases ($\partial\sigma_{yy}/\partial
y<0$) until it plateaus at $y/r^{*}\approx 1$. Thus,
$\partial\sigma_{yy}/\partial y$ reaches its minimum a small distance away
from the obstacle boundary, and our argument predicts that the largest
geometric nonlinearity-generated tensions will appear there. This location is
indeed where the greatest tensile stresses appear in simulations (arrow in
Fig. 3C). We can also compare the magnitude of the tension predicted by Eq.
(16) using $\sigma_{yy}$ from the linear theory with the tension measured in
hydrogel simulations--- however, we find an estimate of the tension that is
approximately five times larger than the maximum tension found via simulations
and scales more gradually with increasing $\Delta/r^{*}$. This discrepancy is
not surprising given the approximations made in this calculation.

###### Acknowledgements.

It is a pleasure to acknowledge John Kolinski, Amaresh Sahu, and Carolyn Bull
for helpful discussions, as well as funding from the NSF through the Princeton
MRSEC DMR-2011750, the Project X Innovation Fund, a Princeton Materials
Science Postdoctoral Fellowship (A.P.), and the Camille Dreyfus Teacher-
Scholar Program of the Camille and Henry Dreyfus Foundation (S.S.D).

## References

  * Ding  _et al._ (2019) M. Ding, X. Han, S. Wang, T. F. Gast, and J. M. Teran, ACM Transactions on Graphics (TOG) 38, 1 (2019).
  * Taylor  _et al._ (2021) I. Taylor, K. Lehner, E. McCaskey, N. Nirmal, Y. Ozkan-Aydin, M. Murray-Cooper, R. Jain, E. W. Hawkes, P. C. Ronald, D. I. Goldman, _et al._ , Proceedings of the National Academy of Sciences 118, e2018940118 (2021).
  * Borsdorf and Haller (2020) A. Borsdorf and A. Haller, in  _The Elgar Companion to Geography, Transdisciplinarity and Sustainability_ (Edward Elgar Publishing, 2020) pp. 140--154.
  * Byrne  _et al._ (2015) R. A. Byrne, M. Joner, and A. Kastrati, European heart journal 36, 3320 (2015).
  * Kuhl  _et al._ (2007) E. Kuhl, R. Maas, G. Himpel, and A. Menzel, Biomechanics and modeling in mechanobiology 6, 321 (2007).
  * Cheng and Zhang (2021) J. Cheng and L. T. Zhang, Journal of Biomedical Science and Engineering 14, 33 (2021).
  * Cabeza  _et al._ (2010) L. F. Cabeza, A. Castell, M. Medrano, I. Martorell, G. Pérez, and I. Fernández, Energy and Buildings 42, 630 (2010).
  * Louf  _et al._ (2021) J.-F. Louf, N. B. Lu, M. G. O’Connell, H. J. Cho, and S. S. Datta, Science Advances 7, eabd2711 (2021).
  * Beebe  _et al._ (2000) D. J. Beebe, J. S. Moore, J. M. Bauer, Q. Yu, R. H. Liu, C. Devadoss, and B.-H. Jo, Nature 404, 588 (2000).
  * Alben (2022) S. Alben, Proceedings of the Royal Society A 478, 20210742 (2022).
  * Chu  _et al._ (2018) E. K. Chu, O. Kilic, H. Cho, A. Groisman, and A. Levchenko, Nature communications 9, 1 (2018).
  * Streichan  _et al._ (2014) S. J. Streichan, C. R. Hoerner, T. Schneidt, D. Holzer, and L. Hufnagel, Proceedings of the National Academy of Sciences 111, 5586 (2014).
  * Bengough  _et al._ (1997) A. Bengough, C. Croser, and J. Pritchard, in  _Plant Roots-From Cells to Systems_ (Springer, 1997) pp. 107--116.
  * Alric  _et al._ (2022) B. Alric, C. Formosa-Dague, E. Dague, L. J. Holt, and M. Delarue, Nature Physics 18, 411 (2022).
  * Ambrosi  _et al._ (2019) D. Ambrosi, M. Ben Amar, C. J. Cyron, A. DeSimone, A. Goriely, J. D. Humphrey, and E. Kuhl, Journal of the Royal Society Interface 16, 20190233 (2019).
  * Goriely (2017) A. Goriely,  _The mathematics and mechanics of biological growth_ , Vol. 45 (Springer, 2017).
  * Zhang and Khademhosseini (2017) Y. S. Zhang and A. Khademhosseini, Science 356, eaaf3627 (2017).
  * (18) M. L. Oyen,  59, 44, publisher: Taylor & Francis _eprint: https://doi.org/10.1179/1743280413Y.0000000022.
  * Deligkaris  _et al._ (2010) K. Deligkaris, T. S. Tadele, W. Olthuis, and A. van den Berg, Sensors and Actuators B: Chemical 147, 765 (2010).
  * Guo  _et al._ (2020) Y. Guo, J. Bae, Z. Fang, P. Li, F. Zhao, and G. Yu, Chemical Reviews 120, 7642 (2020).
  * Tanaka  _et al._ (1987) T. Tanaka, S.-T. Sun, Y. Hirokawa, S. Katayama, J. Kucera, Y. Hirose, and T. Amiya, Nature 325, 796 (1987).
  * Hong  _et al._ (2008) W. Hong, X. Zhao, J. Zhou, and Z. Suo, Journal of the Mechanics and Physics of Solids 56, 1779 (2008).
  * Hong  _et al._ (2009) W. Hong, Z. Liu, and Z. Suo, International Journal of Solids and Structures 46, 3282 (2009).
  * Kang and Huang (2010) M. K. Kang and R. Huang, Journal of Applied Mechanics 77 (2010).
  * Marcombe  _et al._ (2010) R. Marcombe, S. Cai, W. Hong, X. Zhao, Y. Lapusta, and Z. Suo, Soft Matter 6, 784 (2010).
  * Amar and Ciarletta (2010) M. B. Amar and P. Ciarletta, Journal of the Mechanics and Physics of Solids 58, 935 (2010).
  * Dervaux and Ben Amar (2012) J. Dervaux and M. Ben Amar, Annu. Rev. Condens. Matter Phys. 3, 311 (2012).
  * Kim  _et al._ (2010) J. Kim, J. Yoon, and R. C. Hayward, Nature materials 9, 159 (2010).
  * Cho  _et al._ (2019) H. J. Cho, N. B. Lu, M. P. Howard, R. A. Adams, and S. S. Datta, Soft matter 15, 4689 (2019).
  * Hui  _et al._ (2006) C.-Y. Hui, Y. Y. Lin, F.-C. Chuang, K. R. Shull, and W.-C. Lin, Journal of Polymer Science Part B: Polymer Physics 44, 359 (2006).
  * Hu  _et al._ (2010) Y. Hu, X. Zhao, J. J. Vlassak, and Z. Suo, Applied Physics Letters 96, 121904 (2010).
  * Bouklas  _et al._ (2015) N. Bouklas, C. M. Landis, and R. Huang, Journal of the Mechanics and Physics of Solids 79, 21 (2015).
  * Louf and Datta (2021) J.-F. Louf and S. S. Datta, Soft matter 17, 3840 (2021).
  * Gdoutos (2020) E. E. Gdoutos,  _Fracture mechanics: an introduction_ , Vol. 263 (Springer Nature, 2020).
  * Wang  _et al._ (2022) M. Wang, M. Adda-Bedia, J. M. Kolinski, and J. Fineberg, Journal of the Mechanics and Physics of Solids 161, 104795 (2022).
  * Li  _et al._ (2023) C. Li, X. Wei, M. Wang, M. Adda-Bedia, and J. M. Kolinski, Journal of the Mechanics and Physics of Solids , 105330 (2023).
  * Misiewicz  _et al._ (2022) J. Misiewicz, S. S. Datta, K. Lejcuś, and D. Marczak, Materials 15, 4465 (2022).
  * Krafcik  _et al._ (2017) M. J. Krafcik, N. D. Macke, and K. A. Erk, Gels 3, 46 (2017).
  * Kapur  _et al._ (1996) V. Kapur, J. C. Charkoudian, S. B. Kessler, and J. L. Anderson, Industrial & engineering chemistry research 35, 3179 (1996).
  * Dolega  _et al._ (2017) M. E. Dolega, M. Delarue, F. Ingremeau, J. Prost, A. Delon, and G. Cappello, Nature communications 8, 14056 (2017).
  * Lee  _et al._ (2019) W. Lee, N. Kalashnikov, S. Mok, R. Halaoui, E. Kuzmin, A. J. Putnam, S. Takayama, M. Park, L. McCaffrey, R. Zhao, _et al._ , Nature communications 10, 144 (2019).
  * Woodhouse and Johnson (1991) J. Woodhouse and M. Johnson, Agricultural water management 20, 63 (1991).
  * Demitri  _et al._ (2013) C. Demitri, F. Scalera, M. Madaghiele, A. Sannino, and A. Maffezzoli, International Journal of Polymer Science 2013 (2013).
  * Guilherme  _et al._ (2015) M. R. Guilherme, F. A. Aouada, A. R. Fajardo, A. F. Martins, A. T. Paulino, M. F. Davi, A. F. Rubira, and E. C. Muniz, European Polymer Journal 72, 365 (2015).
  * Lejcuś  _et al._ (2018) K. Lejcuś, M. Śpitalniak, and J. Dabrowska, Polymers 10, 271 (2018).
  * Wei and Durian (2013) Y. Wei and D. J. Durian, Physical Review E 87, 053013 (2013).
  * Flory and Rehner (1943) P. J. Flory and J. Rehner, The Journal of Chemical Physics 11, 521 (1943).
  * Flory (1953) P. J. Flory,  _Principles of polymer chemistry_ (Cornell University Press, 1953).
  * Alnæs  _et al._ (2015) M. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells, Archive of Numerical Software 3 (2015).
  * Goriely  _et al._ (2016) A. Goriely, J. Weickenmeier, and E. Kuhl, Physical review letters 117, 138001 (2016).
  * Weickenmeier  _et al._ (2016) J. Weickenmeier, E. Kuhl, and A. Goriely, Journal of the Mechanics and Physics of Solids 96, 572 (2016).
  * Hu  _et al._ (2011) Y. Hu, X. Chen, G. M. Whitesides, J. J. Vlassak, and Z. Suo, Journal of Materials Research 26, 785 (2011).
  * Bouklas and Huang (2012) N. Bouklas and R. Huang, Soft Matter 8, 8194 (2012).
  * Barber (2018) J. R. Barber,  _Contact mechanics_ (Springer, 2018).
  * Johnson (1987) K. L. Johnson,  _Contact mechanics_ (Cambridge university press, 1987).
  * Landau  _et al._ (1986) L. D. Landau, A. Kosevich, L. P. Pitaevskii, and E. M. Lifshitz,  _Theory of Elasticity_ (Butterworth-Heinemann, Oxford, 1986).
  * Audoly and Pomeau (2010) B. Audoly and Y. Pomeau,  _Elasticity and geometry: from hair curls to the non-linear response of shells_ (Oxford university press, 2010).
  * Joshi  _et al._ (2022) C. Joshi, D. G. Goldstein, C. Wennerholm, E. Downey, E. Hamilton, S. Hocking, A. Andrei, J. H. Adler, and T. J. Atherton, arXiv preprint arXiv:2208.07859 (2022).
  * Joshi  _et al._ (2023) C. Joshi, M. Q. Giso, J.-F. Louf, S. S. Datta, and T. J. Atherton, arXiv preprint arXiv:2304.04252 (2023).
  * Lai and Hu (2018) Y. Lai and Y. Hu, Soft Matter 14, 2619 (2018).
  * Kuna (2013) M. Kuna, Solid Mechanics and Its Applications 201, 153 (2013).
  * Mao and Anand (2018) Y. Mao and L. Anand, Journal of the Mechanics and Physics of Solids 115, 30 (2018).
  * Wang and Hong (2012) X. Wang and W. Hong, Soft Matter 8, 8171 (2012).
  * Baumberger and Ronsin (2020) T. Baumberger and O. Ronsin, Mechanics of Soft Materials 2, 1 (2020).
  * Yang  _et al._ (2022) Y. Yang, H. Guo, Z. Du, W. Hong, T. Lu, and T. Wang, Journal of the Mechanics and Physics of Solids , 105007 (2022).
  * Leslie  _et al._ (2021) K.-A. Leslie, R. Doane-Solomon, S. Arora, S. J. Curley, C. Szczepanski, and M. M. Driscoll, Soft Matter 17, 1513 (2021).
  * Helmlinger  _et al._ (1997) G. Helmlinger, P. A. Netti, H. C. Lichtenbeld, R. J. Melder, and R. K. Jain, Nature biotechnology 15, 778 (1997).
  * Zhang  _et al._ (2021) Q. Zhang, J. Li, J. Nijjer, H. Lu, M. Kothari, R. Alert, T. Cohen, and J. Yan, Proceedings of the National Academy of Sciences 118, e2107107118 (2021).
  * Fortune  _et al._ (2022) G. T. Fortune, N. M. Oliveira, and R. E. Goldstein, Physical Review Letters 128, 178102 (2022).
  * Goodwin  _et al._ (2019) K. Goodwin, S. Mao, T. Guyomar, E. Miller, D. C. Radisky, A. Košmrlj, and C. M. Nelson, Development 146, dev181172 (2019).
  * Möbius  _et al._ (2015) W. Möbius, A. W. Murray, and D. R. Nelson, PLoS computational biology 11, e1004615 (2015).
  * Atis  _et al._ (2019) S. Atis, B. T. Weinstein, A. W. Murray, and D. R. Nelson, Physical review X 9, 021058 (2019).
  * Tanaka (1986) T. Tanaka, Physica A: Statistical Mechanics and its Applications 140, 261 (1986).
  * Bertrand  _et al._ (2016) T. Bertrand, J. Peixinho, S. Mukhopadhyay, and C. W. MacMinn, Physical Review Applied 6, 064010 (2016).
  * Curatolo  _et al._ (2017) M. Curatolo, P. Nardinocchi, E. Puntel, and L. Teresi, Journal of Applied Physics 122, 145109 (2017).
  * de Silva and Lapitsky (2016) U. K. de Silva and Y. Lapitsky, ACS applied materials & interfaces 8, 29015 (2016).
  * Zhang  _et al._ (2009) J. Zhang, X. Zhao, Z. Suo, and H. Jiang, Journal of Applied Physics 105, 093522 (2009).
  * Ang  _et al._ (2020) I. Ang, Z. Liu, J. Kim, C.-Y. Hui, and N. Bouklas, Journal of the Mechanics and Physics of Solids 145, 104132 (2020).