Obstructed swelling and fracture of hydrogels

,


I. Introduction
Many growth and expansion processes are sculpted through confinement by rigid obstructions.Familiar examples include muffins rising into their characteristic shape during baking [1], trees growing around boulders [2], and even cities expanding around inhospitable geographic features [3].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 [4][5][6]; conversely, the expansion of spray foam into cracks and in between walls underlies the thermal insulation of many energy-efficient buildings [7].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 [8][9][10][11][12][13][14].Therefore, we ask: are there general principles that dictate how obstructions influence growth and expansion?And if so, how do we discover them?
When the growth/expansion of a body is resisted by surrounding obstructions, large and spatiallynonuniform stresses can develop, influencing subsequent growth/expansion in turn [15,16].Being able to understand the distribution and magnitude of these stresses is thus a necessary step in the development of widelyapplicable, predictive models of growth and expansion.However, model systems in which the coupling between growth and stress can be systematically studied in structured environments 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 [17][18][19][20].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 [21][22][23][24][25][26][27][28][29][30], non-adhered swelling around obstructions has received more limited attention [22,[31][32][33][34][35], often in the context of indentation testing.
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 dramatic 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 complex 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.

II. 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 obs to be placed an equal distance r ctr from a central point; in the experiments, we vary r obs and r ctr between 2.5 − 15 mm and 7.5 − 25 mm, respectively.The cylinders are securely attached to horizontal, parallel laser-engraved acrylic plates spaced vertically by a fixed amount, ∆z = 40 mm (see SI Sec.B for a discussion on the selection of this value).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 ∼ 6 mm (initial state characterized in SI Sec.A) 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 fo-cused 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 fourlobed equilibrium shape.As the 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.However, the deformed 3D shape of the hydrogel makes it challenging to extract quantitative information about relative expansion from a purely top-down view.We leave this effort for future work.
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 sequence 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).The fracturing process is dy- namic: As cracks propagate through compressed areas, releasing stresses, those regions are then able to increase their solvent content and swell.The process eventually stops as the central hydrogel body reaches its final equilibrium degree of swelling.Another representative example of a fracturing hydrogel with a different obstacle geometry is shown in Movie S3.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 [36][37][38].
These observations suggest that, when appropriately 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 obs and spacings r 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 (□), while when the obstacles are closer i.e., r ctr is below a threshold value, the hydrogel repeatedly self-fractures as it swells (×).The fracturing threshold also depends on the size of the obstacles: For a given r ctr , fracturing occurs below a threshold value of r obs .The finding that smaller obstacles promote fracture is reminiscent of wire cutting tests in fracture mechanics, in which thin wires pushed into soft materials induce cutting [39,40].The two geometric parameters r ctr and r 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. Intuitively, this non-trivial dependence of the onset of fracturing on both obstacle spacing and size reflects the balance between the osmotic pressure driving swelling-which is set by the difference between the solvent chemical potential in the unswollen hydrogel compared to at its boundary-and the tensile stresses that develop in the hydrogel as it swells against narrow obstacles [8,41].We theoretically model and computationally quantify this balance in the following sections.
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 [42][43][44][45][46][47][48][49][50][51].Our previous experiments [8] 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 Sec.C).Thus, this phenomenon of obstructed growth/expansion causing fracture arises not just in idealized geometries, but also in more realistic complex spaces.

III. 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.The dimensional reduction is discussed at length in SI Sec.D.
We describe stretching of the hydrogel polymer network and solvent-polymer interactions with the commonly-used free energy density [22-24, 33, 52, 53]: 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 = ∂xi ∂X K is the deformation gradient tensor, with F iK F iK = tr(F T F) = i λ 2 i in terms of principal stretches λ i .Following standard conventions, deformed/current configurations are denoted by lowercase letters and dry reference configurations are denoted by capital letters unless otherwise noted [16].The second term describes the mixing free energy: Ω is the area of a solvent molecule, C is the nominal concentration of solvent (number of solvent molecules per unit dry refer-ence area), and χ 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(F) = 1 + ΩC.We require the chemical potential µ 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 Fen in terms of µ as Fen (x i , µ) = F en (x i , C) − µC [23].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 µ = 0 is enforced in the regions of contact between the hydrogel and the obstacles.It would be interesting to also consider the impact of dynamic boundary conditions in this problem: for example, a boundary condition requiring no normal solvent flux in the regions where the hydrogel contacts obstacles could be imposed [33].Since we focus on equilibrium quantities in this work, we use a constant chemical potential boundary for simplicity.
To visualize and track hydrogel deformation during obstructed swelling, along with the concomitant development of internal stresses, we implement this model in FEniCS [54], 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 Sec.E).
We begin by examining hydrogel swelling around obstacles that have the same size, but varying spacingjust 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 ∆/r * ≡ (r * − r ctr ) /r * , where r * is the equilibrium, fully-swollen, unconfined radius of the hydrogel and ∆ is the difference between r * and the distance from the center to the closest obstacle r ctr .The case of weak confinement in (A) has ∆/r * = 0.02, increasing to ∆/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.III A, the hydrogel stresses can be captured ana-lytically using linear theories in this regime.Increasing hydrogel confinement to the point shown in panel (B) causes the magnitudes of the stresses to increase (panels E,F), but they still remain primarily localized near the obstacles; as discussed in Sec.III B and Sec.III C, the largest compressive stresses can still be described by linear theory, but the largest tensile stresses require consideration of a geometric nonlinearity due to their connection to the rotation of material as it conforms to the obstacle boundary.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.III D. 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.IV.

A. 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 [41].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 + ΩC = det(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 [31,32,55,56].
Assuming that deformations relative to the fully swollen state are small, we linearize Eq. ( 1) to find effective linear elastic parameters in equilibrium (SI Sec.F, [57,58]).In particular, for a 2D hydrogel, the effective equilibrium Poisson's ratio ν and Young's modulus E are given by where λ 0 is the principal stretch corresponding to the fully swollen state.The expression for the Poisson's ra- With increasing ∆/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 Sec.H).For compression, the linear theory is accurate over a larger range of ∆/r * , but the geometrically nonlinear theory cannot explain the deviations; instead, an elastic model incorporating material nonlinearities better captures the data (Eq.( 10), SI Fig. S9) .As ∆/r * increases even further, the hydrogel exhibits a symmetry-breaking instability (D, red shaded region).
tio can be understood as its value for an incompressible 2D solid, ν = 1, minus a correction-reflecting the fact that the compressibility of the swollen hydrogel in equilibrium is generated via solvent transport (i.e., the hydrogel responds to an instantaneous deformation like an incompressible solid before solvent is able to equilibrate).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 equilibrium elastic parameters, we solve for the stresses in the hydrogel as a 2D linear contact mechanics problem (SI Sec.G, [59,60]).This approach provides expressions for the stress tensor σ ij along the line directly beneath the top obstacle as a function of y as shown in Fig. 4B: where ζ < 0 is the force applied to the indenters and the half contact width a is 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 * .Given these analytical expressions, we next ask: Under what confinement regimes (as parameterized by ∆/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 [36]-vary with confinement?To address these questions, we first re-cast Eqs. ( 4)-( 5) in terms of ∆/r * to facilitate direct comparison with the results of the nonlinear simulations.To do so, we integrate the strain u yy = σyy Y − ν Y σ 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 which we then invert, apply Eqs. ( 2)-( 3), and substitute the resulting expression for ζ(∆) into Eqs.( 4)-( 5).The resulting expressions for σ xx (∆) and σ yy (∆) 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 obs /r * = 0.3; for comparison, the symbols show the results of the full nonlinear simulations.
As expected, when ∆/r * ≪ 1, the linearized indentation solution agrees well with the nonlinear simulation results, while as ∆/r * increases, the discrepancy between the two becomes more apparent.In particular, as exemplified by the data for ∆/r * = 0.16 (dark blue circles) and ∆/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 * ≈ 0.75 (Fig. 4A, inset).Indeed, though tension (σ xx > 0) does appear in the linear solutions for small ∆/r * , it disappears with increasing ∆/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 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 Sec.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 ∆/r * ≳ 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 ∆/r * ≈ 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.

B. 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.The effective equilibrium elastic parameters [Eqs.
(2)-( 3)] are independent of strain, 2. The stress is linearly related to the strain, 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 Sec.F, the neo-Hookean model closely reproduces the stress profiles of the hydrogel simulations over a broad range of ∆/r * up to ≈ 0.4, well beyond the limits of the linear theory at ∼ 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-Kirchhoff elastic model with strain tensor u ij = 1 2 ∂ui ∂xj + ∂uj ∂xi + ∂u k ∂xi ∂u k ∂xj ; 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 Sec.H, the St. Venant-Kirchhoff model quantitatively reproduces the maximum principal tensile stress in the hydrogel simulations up to ∆/r * ≈ 0.1, well beyond the limit of the linear theory at ≈ 0.02.Hence, the excess tension σ xx that develops beneath the obstacles just beyond the linear regime is driven by geometric nonlinearity, related to the rotations of the hydrogel material as it accommodates the obstacles during swelling, and does not require a nonlinear constitutive relationship.At even larger displacements ∆/r * > 0.1, the St. Venant-Kirchhoff 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 [61].As detailed in Materials and Methods, we first find the variation of the integrated St. Venant-Kirchhoff strainenergy function with respect to displacements, which can be written in terms of the second Piola-Kirchhoff stress tensor, σ P K ij .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 [62].We then make approximations specific to our obstacle geometry.Ultimately, we find Since σ P K yy becomes more negative as y increases, ∂σ P K yy ∂y < 0 (Fig. 4B), and therefore σ P K xx > 0 indicating tension.Thus, geometric nonlinearity generates tension beneath an indenter perpendicular to the indentation direction, qualitatively matching our simulations.

C. Compression beyond the linear regime
A notable result shown in Fig. 3 is that while linear indentation theory predicts tensile stresses for ∆/r * < 0.02, it captures the compressive stresses over a broader range, up to ∆/r * ≈ 0.2.Intuitively, this robustness is due to the fact that the large compressive stresses appear perpendicular to the obstacles and are relatively unaffected by rotations of the hydrogel material around the obstacles.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 Sec.F & H.In contrast to the case of tension, we do not find any parameters for which the St. Venant-Kirchhoff model is more accurate than the linear model.Furthermore, as shown in Fig. S9, 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 λ 1 = 1 − ∆/r * and the principal stretch perpendicular to the indentation set to 1: (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.

D. 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 symmetrybreaking 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 Sec.J.

IV. 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 swellinginduced 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 obs and r 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 [36], 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 ob-  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.
serve 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.

V. 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 [63,64], both in experiments and in the model, along with a more detailed treatment of the microscopic processes underlying fracturing [65][66][67][68][69][70].Moreover, many more experimental trials near this threshold and with hydrogels of varying mechanical properties, along with high-resolution imaging of crack propagation [71], 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 [72][73][74][75], with potential implications for biological function [76,77].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 [42,[47][48][49][50].They are also widely adopted in other applications, such as oil recovery, construction, mechanobiology, and filtration, all of which involve hydrogel swelling in confinement [43][44][45][46].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.

Experimental design
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 twistand-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 Sec.A).These hydrogel beads were extensively characterized in our previous experiments [8,41], 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 [21,[78][79][80].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 are insufficient to drive fracture as in other less tough gels [71,81].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 the trial shown in Fig. 1, for example, equilibration took 88 h.In other trials where we do not image the entire process of obstructed swelling, we let the hydrogel swell either until fracture has occurred or for at least 40 h past the time to fracture for a trial with the same r obs and a smaller r ctr .If we already recorded an outcome of 'no fracture' for a trial with the same r obs and a smaller r ctr , we waited a minimum of 8 h beyond the equilibration time for an imaged reference trial with similar geometry.Each symbol in Fig. 2 represents a single experiment.

Simulation design
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 Ω = 0.001 and χ = 0.3.The penalty function is integrated over the hydrogel boundary, and is given by where p is the penalty strength, x i is position of the center of the i th obstacle, r 0 is the dry reference radius of the hydrogel, and the sum is over all obstacles.The brackets ⟨• • • ⟩ + take the positive part of the argument, defined as ⟨x⟩ + = 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π/r 4 0 to generate the data shown in this text, and have verified that using p = 62.5π/r 4 0 produces the same results.
Following the suggestions of Refs.[33,82,83], we use the backwards Euler scheme for time integration (see SI Sec.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 σ max , and the most negative eigenvalue is the minimum principal (compressive) stress, σ 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 σ yy and the maximum value of σ 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 (12) Upon integrating by parts, assuming we cannot vary the displacements at the boundaries due to the presence of obstacles, we find 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 σ P K xy = 0 by symmetry, yielding: In equilibrium, the coefficients of δu x and δu y must be zero (in the absence of body forces).We focus our attention on the coefficient of δu y .If we consider the internal stresses directly beneath the obstacle, we can set ∂u y /∂x to zero by symmetry.We approximate the deformation as an affine contraction in the y direction, which sets ∂u y /∂y ≈ −∆/r * and ∂ 2 u y /∂y 2 ≈ 0. We assume that the curvature of the y displacement, ∂ 2 uy ∂x 2 , is determined by the curvature of the obstacle, 1/r obs .With these substitutions, the equilibrium condition becomes Solving for σ P K xx , we find Just beyond the linear regime, σ P K yy should follow the same trends as σ yy predicted by the linear theory.In Fig. 4B, we observe that σ yy decreases as y increases (∂σ yy /∂y < 0) until it plateaus at y/r * ≈ 1.Thus, ∂σ yy /∂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 σ yy from the linear theory with the tension measured in hydrogel simulations-however, we find an estimate of the tension that is approximately twenty times larger than the maximum tension found via simulations and scales incorrectly with increasing ∆/r * .These discrepancies are not surprising given the many approximations made in this calculation.

Material properties
To assess the consistency of material properties between different hydrogels, we measured the density of five white beads in three different hydration states: ambient, dry, and swollen.The "ambient" state is the condition of the beads prior to any experimentation.This state is characterized by a small amount of swelling due to ambient humidity, evidenced by slight pliability but high elastic resistance.The "dry" state was created by placing hydrogel beads in an oven at 100 • C for approximately 24 hours, until the material became completely rigid.These dry beads were then placed in deionized water for 96 hours to reach their fully "swollen" state.Volumes of the spherical beads for each state were extrapolated from the diameters measured with ImageJ.The results of these experiments and calculations are displayed in Figure S1.Dry Ambient Swollen FIG.S1: Density is reasonably consistent between different hydrogels in the same state, despite scatter in their mass and volume.The three sets of points correspond to the dry, ambient, and swollen states of the hydrogel.The slope of the line for each set of points gives the average density for each set of five trials.The average densities are 1.12 g/cm 3 , 0.993 g/cm 3 , and 0.510 g/cm 3 for the dry, ambient, and swollen states respectively.arXiv:2307.11827v2[cond-mat.soft]22 Jan 2024

Size distribution
The hydrogel spheres used in the experiments have slightly different initial masses prior to swelling.This variation in size certainly impacts the location of the reported fracture threshold, as differently sized hydrogels experience different stresses in the same obstacle geometry, and limits our ability to make quantitative statements.In Fig. S2, we display the fracture data from Fig. 2 of the main text, now scaling each point by the initial mass of the hydrogel.The variation in initial mass between all trials is relatively modest, and the masses of the hydrogels in the trials close to the fracture boundary in particular are similar.

B. Fracture due to vertical confinement
Hydrogel swelling in the chambers shown in Fig. 1A of the main text is limited both by the presence of obstacles as well as the top and bottom plates.In order to explore the role played by vertical confinement, we allowed hydrogels to swell in chambers with varying ∆z and no obstacles.
We observed that hydrogels swelling with top and bottom plates separated by ∆z < 20 mm reliably fracture even in the absence of obstacles.In contrast, for ∆z = 40 mm, the separation used for all of the experiments presented in this work, we did not observe fracture in the absence of obstacles.One such trial with ∆z = 40 mm is shown in Fig. S3-note that this particular hydrogel had an initial mass of 1.19 g, larger than any of the hydrogels used to create Fig. 2 of the main text.
Understanding fracture as a function of vertical confinement is an interesting direction for future research.Such a study could build on previous theoretical work on flat plate compression of a sphere [1].

C. 3D packing experiments
Following the work of Louf et al. [2], we allow hydrogel beads to swell in a 3D disordered granular medium composed of borosilicate glass beads of radius 1.5 mm.As described in the original work, beads of this size mimic coarse, unconsolidated soil.
We saturate the packing with deionized water and subject the entire column to an applied load of 12 pounds to prevent restructuring of the granular matrix in response to hydrogel swelling.Since the refractive index of water does not match that of the glass beads, we are unable to collect clear images of the hydrogel during swelling.However, we can disassemble the packing for visual analysis post-swelling.
After a short period of 6 hours, we see that pieces of the hydrogel have already broken off and are distributed throughout the packing (Fig. S4).When this experiment was replicated with extended swelling periods of 12, 24, and 72 hours, increasing degrees of fragmentation were observed.
Hydrogel swelling can clearly generate stresses large enough to cause rupture in 3D granular media, and hydrogel fragmentation therefore may be relevant in agriculture applications.Further studies of hydrogels exposed to water in 3D granular media are an interesting and important extension of our current work.

D. Dimensional reduction
In this section, we provide further discussion on our choice to use a 2D model and discuss the corrections we expect to have due to the 3D nature of our experiments.First, we argue that corrections due to the spherical shape of the hydrogel can be neglected near the midplane and thus locally approximate the spherical hydrogel as a cylinder.Then, we describe how an axially compressed hydrogel cylinder can be modeled in 2D via the generalized plane strain approximation.

Neglect of spherical features
To frame our discussion, we refer to Fig. S5 which shows the side view of a fully-swollen hydrogel (top view shown in Fig. 1B of the main text).From this viewpoint, we clearly see that the spherical shape preferred by the hydrogel leads to z-dependent obstacle-induced stresses-the cylindrical obstacles penetrate further into the hydrogel at the xy plane halfway between the plates compared to the xy plane at the top or bottom plate.In other words, if we consider horizontal slices of the hydrogel, each slice has a preferred radius that changes as a function of z due to the spherical geometry.We can estimate this preferred radius by neglecting the shape distortion caused by the top and bottom plates (i.e., assuming the hydrogel swells to a perfect sphere in the absence of obstacles).Given a preferred sphere radius R * and defining z = 0 to be halfway between the plates, the preferred disk radius r * of a slice of the hydrogel in the xy plane as a function of z is given by Indenter displacement is defined as the difference between the preferred hydrogel radius r * (z) and the leading edge of the obstacle.Since the obstacles are cylinders and their shape/location is independent of z, the z-dependent preferred hydrogel radius generates a z-dependent effective indenter displacement.For a fixed value of r ctr , the dimensionless indenter displacement as a function of z is When corrections of order z 2 /R * 2 can be neglected, we can ignore the z-dependence and locally approximate the hydrogel as a cylinder.This approximation is reasonable close to the midplane of the hydrogel, where the confinement is greatest.Thus, we assume our system meets the criteria of a generalized plane strain approximation near the midplane of the hydrogel at equilibrium.We next describe this approximation, first in the context of linear elasticity and then in the context of the nonlinear hydrogel model.We emphasize that this dimensional reduction only holds in equilibrium, as the dynamic features observed in experiments, including initial swelling as well as fracture itself, have 3D features that cannot be straightforwardly described in 2D.

Generalized plane strain: Hookean elasticity
Generalized plane strain approximations can describe a cylinder in frictionless contact with walls that impose a uniform axial contraction, experiencing tractions independent of z on its sides [3].Accordingly, we assume that the strain tensor elements u xx , u xy and u yy are independent of z, u xz = u yz = 0, and u zz = α, a constant.Assuming Hooke's law holds, we can write down the standard equations for the stress tensor σ ij in terms of the strain tensor u ij , bulk Young's modulus E 3D , and bulk Poisson's ratio ν 3D [4].
Since ∂σzz ∂z = σ xz = σ yz = 0, there are only two equations of equilibrium.Note that α drops out of the equilibrium equations entirely.
To make comparisons with 2D elasticity, it will also be useful to express the strains in terms of stresses.In order to do this, we write σ zz in terms of σ xx and σ yy using the relation for u zz .
Using this relation, the elements of the strain tensor are When α = 0, we recover the standard plane strain result.

Relation to 2D elastic solid
We compare the stress and strain tensors derived using the generalized plane strain approximation to those of a 2D elastic solid.For a 2D material, there is no z direction.The two equilibrium equations are therefore identical to those for generalized plane strain.However, the relationship between stress and strain differs.Strain tensor elements are given in terms of stress tensor elements as where E 2D = 4µ(µ+λ) 2µ+λ and ν 2D = λ 2µ+λ are the 2D Young's modulus and Poisson's ratio defined in terms of 2D Lamé parameters µ and λ.We observe that if we make the substitutions 1−ν 3D and subtract ν 3D αδ ij , we can transform the 2D strain tensor into the generalized plane strain tensor.
Therefore, if we solve for the stresses using a 2D model, we can create solutions to a corresponding generalized plane strain model using the procedure described above.

Generalized plane strain: nonlinear hydrogel theory
Thus far, we have only discussed the generalized plane strain approximation in the context of Hookean elasticity.We show here how the same kinematic assumptions affect the nonlinear hydrogel model when deformations are not necessarily small.The free energy density of a 3D hydrogel is [5,6]: where F 3D is the 3D deformation gradient tensor, v is the volume of a solvent molecule, n 3D p is the number of polymer chains per unit reference volume V 0 , and C 3D is the 3D nominal concentration (number of solvent molecules per unit reference volume).
The Cauchy stress is given by where we define the function We now apply our generalized plane strain assumptions.The deformation gradient tensor and its inverse can be written where J is the determinant of the 2D deformation gradient tensor F iK and λ z = 1 + α is the axial stretch.Thus, J 3D = λ z J. Constrained hydrogel swelling is well-studied-Ref.[7] provides a particularly relevant analysis of indentation on swollen constrained hydrogels.With these assumptions, σ xz = σ yz = 0, and σ zz is given by As in Hookean elasticity, this stress is independent of z and fully determined by the strains in the x and y directions and the prescribed stretch in the z direction.Therefore, ∂σzz ∂z = 0 and there are again only two equilibrium equations to solve.In the absence of body forces, these are ∂σxj ∂xj = 0 and ∂σyj ∂xj = 0.The other nonzero elements of the Cauchy stress tensor can be written In the absence of obstacles, the hydrogel will swell isotropically in x and y, reaching an equilibrium swelling stretch ratio λ 0 that will be a function of λ z .We find λ 0 by setting σ ij = 0 with F iJ = λ 0 δ iJ , J = λ 2 0 .This gives the relation This equation must be solved numerically.The relationship between λ z and λ 0 for µ 3D = 0, n 3D p v = 0.001, χ = 0.3 is shown in Fig. S6.When λ z = λ 0 , we regain the standard 3D result [8].

Relation to 2D hydrogel
We compare the stress tensor for a hydrogel in generalized plane strain to that of the purely 2D hydrogel described in the main text.For the energy functional The corresponding Cauchy stresses are where We can use 2D solutions to construct generalized plane strain solutions with arbitrary λ z .If we make the following substitutions in Eq. (S24), we regain Eq. (S21).
Note that, as in the Hookean elasticity example, these substitutions change the units of various terms-this occurs because 2D and 3D stress have different units.

Summary
Our experimental system can be reasonably modeled in 2D using generalized plane strain assumptions near the hydrogel midsection in equilibrium.In the main text, we work with 2D elastic models to simplify our calculations and simulations.As described above, it is straightforward to transform 2D elasticity solutions to generalized plane strain solutions and vice versa.The imposed stretch in the z direction enters as an effective chemical potential in the hydrogel model, so chemical potential boundary conditions in the 2D simulations would need to be modified to make more quantitative comparisons with experiments.

E. Dynamics
We simulate hydrogel swelling dynamically using the kinetic law proposed in Hong et al. [5].
where j i and c are the solvent flux and concentration in the current configuration respectively, µ is the chemical potential, and D is solvent diffusivity.To evolve the concentration field, we discretize the continuity equation using the concentration and solvent flux in the reference configuration: ∂C ∂t + ∂J K ∂X K = 0, with J K = det(F) ∂X K ∂xi j i [9].We do not expect the dynamical behavior of our simulations to closely reproduce what we observe in our experiments for several reasons.First, the kinetic law we use is simple and does not account for effects such as changes to the diffusion constant as the polymer network expands.See, e.g., [10][11][12], for further discussions about appropriate expressions for solvent flux.Second, before the hydrogel makes contact with the top and bottom plates, it swells equally in all three dimensions.A 2D simulation therefore only provides a reasonable approximation of equilibrium experimental behavior, as discussed in ESI Sec.D.
Nevertheless, the simulated 2D dynamics are interesting and complex, and are relevant to hydrogel geometries that remain quasi-2D throughout the swelling process, such as a thin hydrogel disk between two flat plates.We therefore show some characteristic results and provide a brief discussion on simulated dynamical behavior here.In all cases, we initialize the model at a radius r 0 and impose an isotropic initial stretch λ x = λ y = 1.5 to avoid the singularity that appears for the dry state with λ x = λ y = 1 [8].As described in Materials and methods, simulations are run until they reach equilibrium which occurs at different times for different trials.In Fig. S7, we show the maximum principal tensile (top set of curves) and compressive (bottom set of curves) stresses as a function of time for different values of ∆/r * with r obs /r * = 0.3.These data are a subset of the data displayed in Fig. 5 of the main text, with maxima/minima taken over the entire mesh.
To understand these data, we use the ∆/r * = 0 curve as a reference.This curve corresponds to the unconstrained case.We observe a peak in both the maximum compressive and tensile principal stresses at dimensionless time t ≡ tD r 2 0 ≈ 4, followed by a decay to zero.This peak in stress occurs because the outer regions of the hydrogel swell faster than the inner regions, and can lead to the formation of transient cusp structures, as described in the main text and elsewhere [13].For ∆/r * ≤ 0.5, simulations follow the free swelling curve until they make contact with the obstacles.Following contact, the hydrogel develops stresses that do not disappear at equilibrium.The data for ∆/r * = 0.6, the strongest confinement shown here, has a number of interesting differences from the other trials.This hydrogel encounters obstacles prior to the time points shown in Fig. S7, and thus sustains larger stresses earlier compared to the other trials.At t ≈ 700, we see a sharp increase in both the compressive and tensile stresses in this trial, not observed in other trials-this is the symmetry-breaking instability discussed in Sec.IIID of the main text.We compare these dynamics to those in Fig. S8.Here, we show principal stresses as a function of time with ∆/r * held fixed and r obs /r * varied, again with maxima/minima taken over the entire mesh.These data correspond to vertical slices of Figs.5(A,B).Since ∆/r * is the same for all trials, the hydrogels all encounter obstacles at the same time, t ≈ 6.After that time, trials behave differently due to the different curvature of the obstacles.

F. Elastic moduli
In this section, we provide the derivation for the effective Poisson's ratio and Young's modulus of a 2D hydrogel in equilibrium (Eqs.( 2) and (3) of the main text) by adapting the argument presented in Bouklas and Huang [14].We begin with the Cauchy stress for a 2D hydrogel model, with We apply a uniaxial stress in the x direction and require σ yy = 0.This sets We assume that the resulting stretches are small deformations relative to the stress-free equilibrium state with stretch λ 0 , defined for the 2D hydrogel via the equation This allows us to define the strain tensor elements in terms of stretches as  Even for these large deformations, we observe reasonably good agreement between the two models.Note that we plot both the maximum/minimum principal stress over the entire hydrogel domain (labeled global max./min.)as well as the maximum/minimum principal stress along the line connecting the hydrogel center to an obstacle center (labeled below obstacle max./min.).For further discussion of the differences between these two measures, see ESI Sec.I.In panel B, the dashed line shows Eq. 10 of the main text for comparison.

G. Indentation solutions
Here we derive the indentation theory equations in Sec.IIIA of the main text.We start with the Flamant solution applied to a point force ζ acting normal to an elastic half space in 2D [15].We assume the half space fills the region y < 0 and use the convention that ζ < 0 corresponds to a force that compresses the plane.In radial coordinates with the origin at the location at which the point force acts, this is a simple radial stress distribution Directly beneath the point source, θ = −π/2 and σ rr = ζ πr < 0. The displacements due to a point force grow logarithmically in 2D.Therefore, we must be mindful of boundaries, as the finite size of our deformable material will influence our calculations.We can gain an appreciation for this subtlety by considering two diametrically opposed point forces acting on the top and bottom of a 2D elastic disk, as in Goodier and Timoshenko [16] p 107.In this case, all points on the boundary of the disk experience isotropic compression due to the pair of point forces, and a uniform tension σ ij = − ζ πr * δ ij must be added to maintain a stress-free boundary.By the same argument, σ ij = − 2ζ πr * δ ij must be added to the stress tensor to free the boundaries when four perpendicular point forces are acting on the elastic disk.
We now consider the specific geometry we are interested in: four circular indenters of radius r obs acting on an elastic disk of radius r * .In the weak confinement regime, we expect the greatest stresses to develop on the line connecting the center of the elastic disk and the center of an indenter, as stresses go to zero at the points at which the disk and indenter lose contact.Therefore, by symmetry we only need to solve for stresses along this line (x = 0 line in inset to Fig. 4B in the main text).Appealing to Saint-Venant's principle, we simplify our calculation by modeling one obstacle as a 2D Hertzian indenter and the other three obstacles as point forces.We proceed by generalizing the calculation in Johnson [17] p129.
We place the origin at the center of the contact between the elastic body and the Hertzian indenter, and let b > 0 be the distance from this origin along the centerline with x = 0.The stress tensor contributions from the point indenters along the line of interest are ) The Hertzian indenter contributes ) where the contact length 2a is given according to The stress tensor is the sum of these contributions, plus the isotropic tension from the corrective solution, σ ij = − 2ζ πr * δ ij .With the substitution y = r * − b, we find Eqs.( 4)- (7) of the main text.The strain u yy can be found as We find the displacement at the surface directly beneath the indenter by integrating We expand this expression in the limit a/r * ≪ 1, neglecting terms quadratic in a/r * to find This expression can be inverted using the Lambert W-function or product logarithm ( [18], §4.13).Note that this function is multivalued for small indenter force and the k = −1 branch must be chosen for consistent results.

H. St. Venant-Kirchhoff model
As described in Sec.IIIB of the main text, we simulate St. Venant-Kirchhoff materials surrounded by four circular indenters in FEniCS.We use the same mesh as in the hydrogel trials (approximately 30 vertices across the radius; see Fig. 3, left column).The strain-energy density function is where 1−ν are Lamé coefficients, set to match the effective elastic properties of the hydrogel as described in ESI Sec.F, u i is a displacement in the ith direction, and indices run over x and y.Displacements are defined relative to the zero stress reference configuration with radius r * .
As the indenter displacement increases, St. Venant-Kirchhoff simulations can become unstable.For example, due to the linear constitutive law, there is a finite energetic cost for compressing material to a point [19].To maintain stability when possible, we incrementally increase both indenter displacement and penalty strength.Nonetheless, we cannot simulate values of ∆/r * > 0.115 for r obs /r * = 0.3.Results from these simulations prior to this threshold are shown in Fig. S11  As discussed in the main text, the St. Venant-Kirchhoff model is able to capture the maximum tensile stress reasonably well.However, the compressive stresses in the hydrogel model (and neo-Hookean model, see Fig. S10B) are significantly larger than those in the St. Venant-Kirchhoff model for large deformations.

I. Maximum tensile stresses in weak confinement
In Figs.S9 and S11, the global maximum of principal tensile stress is larger than the theoretical expectation at small indentation depths.These deviations are the result of isolated cells experiencing large tensile stresses on the boundary of the hydrogel where the material loses contact with the obstacle, and we consider it to be an unphysical effect originating with our finite mesh resolution.In Fig. S13, we show how the tension changes as we increase the resolution.At the relatively coarse resolution used for the simulations in this work with 30 vertices along the radius, the anomalous tension appears close to the boundary, while at higher resolutions it is more reliably located on the boundary itself.We compare the maximum principal stress excluding the boundary points (blue crosses in Fig. S13) to the maximum principal stress including the boundary and the maximum tensile stress beneath the obstacle.We find good agreement between the maximum stress excluding the boundary and the theoretical expectation at high resolution.However, since the maximum compression occurs close to the boundary as well, excluding these points systematically creates disagreement between theory and simulation for compressive stresses.Thus, in Fig. 3 we display the maximum and minimum stresses taken along a line connecting the obstacle center to the hydrogel center beneath the top obstacle, as decribed in Materials and methods.

J. Symmetry-breaking instability
When obstacles are very close together (∆/r * > 0.56 for r obs /r * = 0.3), simulations show that the hydrogel disk swells primarily along a diagonal, rather than equally in all four pore spaces.This effect is quantified in Fig. S14.
As discussed in the main text, this instability can be understood as the hydrogel prioritizing an elliptical shape over a circular shape.In the weak confinement regime, deformations are relatively localized.Therefore, by symmetry, the hydrogel will swell evenly into all four pores.However, as the confinement increases and the deformation becomes global, maintaining symmetric swelling requires isotropic compression of the center of the hydrogel.Assuming uniform deformations and reasonable elastic parameters, it is less costly for the hydrogel to compress an amount ∆ along a single axis, forming an ellipse, than to compress an amount ∆ along two axes, forming a smaller circle.For concreteness, we demonstrate this with a compressible neo-Hookean model (Eq.S45).Compressibility results from the migration of solvent molecules, as shown in ESI Sec.F. To transform a circle into a smaller circle, we impose  .

FIG. 1 :
FIG. 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 obs separated according to r ctr and confined vertically by parallel plates separated by ∆z.(B-C) Top view images taken over the course of swelling, with r obs = 5 mm.The approximate borders of the obstacles and hydrogels are outlined for clarity.(B) When the obstacles are further apart (r ctr = 20 mm), a hydrogel reaches a four-lobed clover-like shape at equilibrium.(C) When the obstacles are closer together (r ctr = 10 mm), cracks appear at the surface of the hydrogel as it swells, driving the repeated production of discrete fragments.

FIG. 2 :
FIG.2: Experiments reveal a hydrogel fracture threshold that depends on obstacle radius and spacing.Each □ indicates an experiment in which the hydrogel reached an intact equilibrium shape as in Fig.1B, while each × 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.

FIG. 3 :
FIG. 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 dry 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 Ω 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 obs /r * = 0.3; thus, ∆/r * = 1 − r ctr /r * = 0.02, 0.1, 0.5, 0.6 from top to bottom, where ∆ is the difference between the unconfined swollen hydrogel radius r * and the distance to an obstacle from the center r 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 σ 1 , σ 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 obs /r * = 0.3 as a function of ∆/r * .Stresses are again normalized by G 0 .The values of ∆/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 ∆/r * (purple shaded region).With increasing ∆/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 Sec.H).For compression, the linear theory is accurate over a larger range of ∆/r * , but the geometrically nonlinear theory cannot explain the deviations; instead, an elastic model incorporating material nonlinearities better captures the data (Eq.(10), SI Fig.S9) .As ∆/r * increases even further, the hydrogel exhibits a symmetry-breaking instability (D, red shaded region).

9 FIG. 4 :
FIG.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 σ xx and σ 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 obs /r * = 0.3.The inset to (B) defines the variables used: the hydrogel-obstacle contact width is 2a and the indenter displacement is ∆, defined relative to the undeformed hydrogel radius r * (dashed orange line).As the indentation depth ∆/r * ≡ 1 − r ctr /r * increases, the discrepancy between the linear theory and simulations increases, as expected.With increasing confinement, tension (σ xx > 0) builds up at r ctr /r * ≈ 0.75, as shown by the magnified inset in A.

FIG. 5 :
FIG.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 σ max and magnitude of the most negative (compressive) principal stress |σ 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.
FIG.S2: Fracture threshold as a function of obstacle geometry as in Fig.2of the main text.The size of the points gives the initial mass of the hydrogel according to the key in the top right, with radius scaled according to mass.The dark green/light green points correspond to trials that fracture/do not fracture respectively, with the dotted line showing the approximate boundary between the two behaviors.The triangular region on the top left shows parameters for which obstacles would overlap.
FIG.S3: A fully swollen hydrogel remains intact when compressed between two parallel plates spread 40 mm apart.
FIG. S5: Side view of the hydrogel and obstacles pictured in Fig. 1B of the main text at ∼ 90 hours.
FIG.S7: Maximum tensile and compressive principal stresses as a function of time (rescaled by r 0 , the dry hydrogel radius, and D, the diffusivity) for r obs /r * = 0.3, varying ∆/r * .
FIG. S8: Maximum tensile and compressive principal stresses as a function of time (rescaled by r 0 , the dry hydrogel radius, and D, the diffusivity) for ∆/r * = 0.5, varying r obs /r * .The light blue curves in this figure and Fig. S7 display the same data.
FIG. S9:Comparison between maximum principal Cauchy tensile and compressive stresses in simulations with the nonlinear hydrogel model and neo-Hookean elastic model.Results for the neo-Hookean model are shown through ∆/r * = 0.48.Even for these large deformations, we observe reasonably good agreement between the two models.Note that we plot both the maximum/minimum principal stress over the entire hydrogel domain (labeled global max./min.)as well as the maximum/minimum principal stress along the line connecting the hydrogel center to an obstacle center (labeled below obstacle max./min.).For further discussion of the differences between these two measures, see ESI Sec.I.In panel B, the dashed line shows Eq. 10 of the main text for comparison.

47 FIG
FIG.S10: Comparison between Cauchy stress profiles in simulations with the nonlinear hydrogel model and neo-Hookean elastic model, plotted against undeformed fully-swollen coordinates (corresponding to an unobstructed equilibrium state for the hydrogel).Neo-Hookean results are shown with large colored markers, while hydrogel results are shown in grey.Data corresponding to the same value of ∆/r * have the same marker shape.Although we clearly see differences between the two models at large ∆/r * , the two models produce similar stress profiles.

11 FIG
FIG. S11:Comparison between maximum principal Cauchy tensile and compressive stresses in simulations with the nonlinear hydrogel model and St. Venant-Kirchhoff model.In (A), we observe that the St. Venant-Kirchhoff model produces approximately the same tensile stresses as the hydrogel model, implying that a geometric nonlinearity can explain tensile stresses up to ∼ ∆/r * = 0.115.In (B), we see that the St. Venant-Kirchhoff model displays compressive stresses that are lower than the hydrogel model, implying that a material nonlinearity is responsible for the additional compressive stress.Note that we plot both the maximum/minimum principal stress over the entire hydrogel domain (labeled global max./min.)as well as the maximum/minimum principal stress along the line connecting the hydrogel center to an obstacle center (labeled below obstacle max./min.).
FIG. S13: Large tensile stresses appear at isolated points near where the hydrogel loses contact with an obstacle.(A) Positive values of the principal stresses for ∆/r * = 0.98, r obs /r * = 0.3, plotted on top of the deformed hydrogel mesh and rescaled by the fully swollen shear modulus G 0 .(B) We compare the maximum principal stresses shown in Fig. 5 of the main text (purple circles) to the maximum principal stress excluding the boundary cells, taken by excluding cells at positions ≥ 95% of the hydrogel radius.Panels (C, D) show the same calculations at twice the resolution.

r 1 r 2 δr = r 1 * 2 − 2 −
FIG. S14:As ∆/r * increases (the obstacles are brought closer together), the hydrogel begins to swell preferentially along a diagonal.The difference in the maximum distance from the origin of the top right and top left hydrogel lobe, δr, provides an order parameter for this transition.Insets show the boundary of the hydrogel at the indicated points (obstacles not pictured). ) )