Buckling initiation in layered hydrogels during transient swelling

Subjected to compressive stresses, soft polymers with stiffness gradients can display various buckling patterns. These compressive stresses can have different origins, like mechanical forces, temperature changes, or, for hydrogel materials, osmotic swelling. Here, we focus on the influence of the transient nature of osmotic swelling on the initiation of buckling in confined layered hydrogel structures. A constitutive model for transient hydrogel swelling is outlined and implemented as a user-subroutine for the commercial finite element software Abaqus. The finite element procedure is benchmarked against linear perturbation analysis results for equilibrium swelling showing excellent correspondence. Based on the finite element results we conclude that the initiation of buckling in a two-layered hydrogel structure is highly affected by transient swelling effects, with instability emerging at lower swelling ratios and later in time with a lower diffusion coefficient. In addition, for hard-on-soft systems the wavelength of the buckling pattern is found to decrease as the diffusivity of the material is reduced for gels with a relatively low stiffness gradient between the substrate and the upper film. This study highlights the difference between equilibrium and transient swelling when it comes to the onset of instability in hydrogels, which is believed to be of importance as a fundamental aspect of swelling as well as providing input to guiding principles in the design of specific hydrogel systems.


Introduction
Hydrogels are polymer networks swollen in an aqueous solution. Natural entities, e.g. extracellular matrix or the vitreous body of the eye, can be understood within this concept. Similar, various synthetic hydrogels are exploited in commercial products due to their particular properties (e.g. disposable diapers and contact lenses), and have a significant potential for applications like smart valves [1], tissue engineering [2,3,4], drug delivery systems [5,6,7], and biological sensors [8,9]. For many of these applications, the capability of hydrogels to swell or shrink as a response to stimuli (e.g. changes in temperature, mechanical forces, pH, salinity level, electric field, specific molecules recognized by included capture moieties) is exploited. During this swelling process, instabilities can occur at the surface of the gels causing wrinkling (i.e. global harmonic waves) and/or creasing (i.e. localized sharp folds) [10,11,12,13,14,15,16,17,18,19,20], as illustrated in Figure 1. While creasing mainly occurs as the first mode of instability when homogeneous gels are exposed to large swelling [12,18], wrinkling patterns are first and foremost found in gels having a gradient stiffness through the thickness, caused by a variation in the crosslinking density within the gel [14,21,22], or by the deposition of a thin and stiff film at the outer surface of the gel, possibly to alter properties like permeability, stability or biocompatibility [23,24,25]. In addition, the impact of the mechanically layered or anisotropic character of natural entities on morphogenesis attains growing interest [26,18] with structure formation in aging skin as an example [27].
A phenomenon closely related to swelling induced instability in layered hydrogels is the evolution of buckling patterns in layered polymer plates under mechanical loading [28,29,17,30]. Assuming linear elastic material behavior, analytical expressions can be found for the critical strain at the onset of buckling and the wavelength of the resulting wrinkling pattern based on the stiffness ratio between the film and the substrate and the thickness of the film [31]. Along the same lines, linear perturbation analysis (LPA) frameworks have been developed for studying the onset of instability in layered hydrogel systems exposed to equilibrium swelling [32,33], i.e. assuming a homogeneous field for the chemical potential in the gel at all times. While these analytical or semi-analytical approaches are computationally efficient, they neglect the transient nature of the osmotic swelling process in hydrogels, giving a time-dependent and inhomogeneous field of the chemical potential through the gel, which could possibly alter buckling initiation. Accounting for the effects of transient swelling in layered gels calls upon numerical procedures like the finite element method (FEM).
For finite element simulations of the equilibrium swelling of gels, numerous constitutive models are available in the literature, both general formulations defining the gel by its basic properties [34,35] and formulations that relate the parameters in the simulations more explicitly to the network properties of the specific gels [36,24]. To also describe the transient nature of hydrogel swelling there are generally two methods that are used in the literature. The most available approach is to utilize the analogy between gel diffusion and heat transfer and implement the swelling process using the thermal modeling capabilities already available in commercial finite element codes [37,38]. The other and more general approach is to formulate diffusion-deformation specific elements [39,40,41], however, the implementation of such elements can be considered challenging and time-consuming.
The studies in the literature using finite element simulations for predicting the onset of instability mainly focus on buckling caused by equilibrium swelling (i.e. omitting transient effects) using 2D [42,32,43,44] or 3D models [45,46,47]. Transient effects related to swelling induced buckling are less described. Bouklas et al. [41] briefly discussed simulations of transient swelling, however, they did not demonstrate the difference between equilibrium and transient swelling for the initiation of buckling. Toh et al. [48] only considered plates with homogeneous material properties, in addition, they used a multi-point constraint routine to force a sinusoidal buckling pattern and hence neglected creasing as a possible mode of instability. In their recent study, Yu et al. [49] investigated buckling initiation in a 2D model of a thin constrained plate exposed to transient swelling, however, they did not consider gels with a stiffness gradient. Finally, Dortdivanlioglu and Linder [50] studied the initiation of buckling in layered confined hydrogels during transient swelling, though, they did not discuss the dependence of the diffusion coefficient for the initiation of instability nor the effective swelling ratio at the onset of buckling. In addition, none of the mentioned studies discuss the onset of instability during transient swelling of soft-on-hard layered gels.
Although buckling caused by mechanical compression or equilibrium swelling of soft layered materials is well studied using experimental, analytical, or numerical procedures, the effect of the transient nature of hydrogel swelling for the onset of buckling remains insufficiently understood. Hence, we address the effects of transient swelling on the critical swelling ratio, the time to buckling, and the wavelength of the resulting buckling pattern for confined hydrogels. Furthermore, we study both hard-on-soft systems, with stiffness ratios between the film and the substrate in the moderate (20, 100) and low range (2,5), and a soft-on-hard system with a stiffness ratio of 0.5. Due to its importance for the presented simulation results, we give a detailed discussion on the introduction of surface imperfections, the quantification of onset of instability, and the use of a plane strain assumption. For the case of equilibrium swelling, we benchmark our model against an LPA framework available in the literature [32], obtaining excellent correspondence. The results of this study are believed to be of importance for the development of models that can predict hydrogel swelling behavior in a general manner, as well as providing valuable insight into the nature of swelling induced buckling processes in layered hydrogels, which can affect guiding principles for the design of hydrogel systems.
The article is organized as follows: In the next section, the geometrical properties of the problem studied herein are described. In Section 3, the constitutive formulation used for transient hydrogel swelling is presented. Section 4 briefly describes the LPA approach used to predict the onset of buckling during equilibrium swelling as a benchmark for the finite element approach. Thereafter, in Section 5, the finite element modeling is described in detail, including the implementation of the constitutive model for the commercial finite element software Abaqus and a validation case. Section 6 presents and discusses the results obtained from the finite element procedure and compares them to the linear perturbation results for the case of homogeneous swelling. Finally, conclusions from the study are presented in Section 7.

Problem definition
To study the onset of swelling induced instability in constrained layered hydrogel plates we define a square plate as illustrated in Figure 2, having edges of length L and a total thickness of H. The upper part of the plate consists of a film with a contrast stiffness and a thickness T . In the present study, the lengths of the plate are set to L = 8 mm, the total plate thickness is set to H = 0.5 mm, while the thickness of the film is set to T =50 µm. We introduce the film to plate thickness ratio η = T /H, and note that a value of η = 0.1 is used herein. The influence of the chosen dimensions on the obtained results is discussed in Section 6.6.
All lateral edges of the plate are constrained from in-plane deformation, while there are no restrictions to swelling in the out-of-plane direction. At its upper surface, the plate is exposed to a solvent. Hence, a change in the chemical potential in the surrounding of the gel will cause an inhomogeneous and timedependent profile for the chemical potential through the thickness of the gel. This chemical potential can cause swelling and hence give the plate a new and time-dependent total height of h. The chemical potential profile through the thickness of the gel will depend on the diffusion coefficient of the solvent molecules. As the goal of this study is to investigate the effect of transient swelling on the onset of buckling, we study the described problem with various values for the diffusion coefficient of the solvent molecules. Calculations are performed for diffusion coefficients between 10 −11 m 2 /s and 1 m 2 /s. While 10 −11 m 2 /s represents the lower range of physical values for a gel swelling in water [51], 1 m 2 /s is an unrealistic value, however, it is used to obtain buckling results when the chemical potential is approximately homogeneous through the gel.  Figure 2: Illustration of a layered hydrogel plate in the coordinate system x where only the upper free surface is exposed to a solvent.

Equilibrium swelling
The constitutive modeling of hydrogel behavior applied herein is based on the work by Hong et al. [52,34], Kang and Huang [35], and Toh et al. [37]. The free-energy function for the hydrogel is assumed to originate from the additive contributions of stretching of the polymer network and mixing of the polymer and the solvent molecules [53,54,55] where N is the number of polymeric chains per reference volume, kT is the temperature in the unit of energy, F is the deformation gradient tensor, I 1 = trb is the first invariant of the left Cauchy-Green tensor b = FF T , J = det F is the volume ratio of the material, v is the volume per solvent molecule, C is the nominal concentration of solvent molecules, and χ is the Flory-Huggins parameter.
By assuming that both the polymer network and the solvent molecules retain their volumes through the swelling process, we find that the volume increase of the gel only can come from an increase in the number of solvent molecules inside the gel, hence we can write Further, to enable implementation for the finite element method, we can introduce a new free-energy function W by the use of a Legendre transformation aŝ and hence ensure the deformation gradient F and the chemical potential µ to be the two independent variables of the model. The Cauchy stress tensor σ can then be obtained from the free-energy potential function as where I is the second-order identity tensor. Normalized quantities are applied in the following and are defined according toW =Ŵ /N kT ,σ = σv/kT , andμ = µ/kT .

Kinetics
In order to capture transient swelling in the constitutive model, the migration of solvent molecules into the hydrogel network must be accounted for. Here, we adopt the modeling procedure by Hong et al. [52] and assume that the small solvent molecules diffuse in the hydrogel network. A kinetic law for the gel can be written as where Φ (X, t) is the nominal flux vector for the solvent molecules, X denotes the position of the material points in the reference configuration, while M is the mobility tensor. Further, the diffusion coefficient of the solvent molecules D is taken as isotropic and independent of the deformation gradient F and the nominal concentration of solvent molecules C [52]. The true flux vector φ(x, t) is assumed to follow Fick's first law [56] where x = x (X, t) is the position of the material points in the current configuration. The nominal and true flux vectors can be related through Combining Equations (2), (5), (6), and (7) we find that the mobility tensor can be written as It can be noted that the outlined material model neglects possible viscoelastic effects in the polymer network. The benefit of this approach is that the transient swelling is the only time-dependent feature of the model, and it is thereby easy to isolate the transient effects. The impact of this simplification on the obtained results is discussed in Section 6.4.

Linear perturbation analysis
Based on the work by Wu et al. [32] for a bilayer hydrogel we implemented an LPA framework for studying the onset of mechanical instability during equilibrium swelling (i.e. assuming a homogeneous chemical potential). The procedure is based on a 2D perturbation analysis where the perturbed deformation gradientF is assumed to take the form letting λ be the out-of-plane deformation, while the assumed perturbations from the equilibrium state in the x 1 -and x 3 -direction are given by u 1 = U 1 (x 3 ) sin (ωx 1 ) and u 3 = U 3 (x 3 ) cos (ωx 1 ). The further steps of the method are thoroughly described in Wu et al. [32] and its implementation culminates in solving the determinant of an 8×8 matrix equal to zero, where all nonzero elements are given between Equations (43) and (44) in Wu et al. [32]. A stability plot, showing the estimated critical swelling ratio at the onset of instability for different wavelengths of the initial harmonic perturbation, is given in Figure 3. Here, n denotes the stiffness ratio between the film and the substrate, n = N v f /N v s , with N v f and N v s being the stiffness of the film and the substrate, respectively. For readability, only results with the stiffness of the substrate set to N v s = 0.001 is shown, however, similar trends would be found for N v s = 0.01. The analysis is based on the plate height and film thickness as given in Section 2, while the plate width is assumed infinite in the LPA framework.
For the soft-on-hard system (n = 0.5) we can see that the minimum point for the critical swelling ratio occurs at the zero-wavelength limit. This indicates creasing as the dominating mode of instability and a similar response was found for soft-on-hard systems with stiffness ratio values of 0.1 and 0.8 (results not shown here due to readability). For the hard-on-soft systems (n > 1) on the other hand, the minimum point of the critical swelling ratio is found at a defined wavelength, indicating wrinkling to be the dominating mode of instability. This difference between soft-on-hard and hard-on-soft systems is in accordance with previously reported results [32].
From Figure 3 it can also be noted that for hard-on-soft systems the critical swelling ratio increases as the stiffness ratio is reduced. Further, the instability curves for the hard-on-soft gels display a wider plateau in buckling wavelength around the critical point as the stiffness ratio is increased. Hence, it can be expected that the critical wavelength obtained in experiments would be more sensitive to experimental uncertainties, like the initial surface geometry, as the stiffness ratio is increased. Along the same lines, we can expect that the critical wavelength obtained using numerical methods as FEM could be more sensitive to numerical approximations as the stiffness ratio is increased.
To study the effect of transient swelling for the initiation of buckling, we will make use of the finite element method as described in the following section. The swelling ratio and the wavelength obtained by the LPA will be used as a benchmark for the finite element simulations with a large diffusion coefficient (i.e. 1 m 2 /s), causing a nearly homogeneous chemical potential through the thickness of the gel, to mimic the assumptions of the LPA calculations.

. Implementation
To accommodate for finite element simulations of the defined swelling problem, the constitutive model described in Section 3.1 was implemented as a Fortran routine to be used with the commercial finite element software Abaqus/Standard. The equilibrium behavior of the model is implemented as a user-defined material, UMAT [57], defining the equations for the normalized stressσ and the Jacobian tangent stiffness. As the chemical potential is singular in the dry state of the hydrogel Ω 0 , an intermediate configuration Ω 1 is introduced. Consequentially, the full deformation tensor F can be split in a multiplicative manner as F = F 1 F 0 . F 0 maps the reference configuration to the intermediate configuration, and F 1 maps the intermediate configuration to the current configuration, as illustrated in Figure 4. In the intermediate configuration, the gel is assumed stress-free and in a state of isotropic strain such that F 0 = λ 0 I. In the present study, all finite element simulations start in the intermediate configuration Ω 1 with a homogeneous distribution of the normalized chemical potentialμ 1 in the gel and dimensions as given in Section 2. The value of λ 0 is found by numerically solvingσ = 0 at the start of the finite element simulation for the initial chemical potential and the material parameters of the hydrogel.
To account for the transient behavior of the gel swelling process, a coupled solvent diffusion and large deformation procedure is utilized. This procedure follows the main ideas presented by Toh et al. [37] and a detailed description is included in Appendix A.
The implemented Fortran code and an example input file for Abaqus is made available as a Mendeley dataset linked to this work [58].

Validation
To validate the implemented model, we study the case of 1D swelling of a homogeneous perfect plate. We start from an intermediate state Ω 1 where the gel is in a state of isotropic swelling (F 0 = λ 0 I) so that it equilibrates at a homogeneous chemical potential ofμ 1 = −2. The plate is then confined from in-plane expansion but can still freely swell out-of-plane. The chemical potential at the surface of the plate is changed toμ (X 3 = H, t) = 0 and the gel gradually swells in the out-of-plane direction giving a deformation gradient with the non-zero elements F 11 = F 22 = λ 0 and F 33 = λλ 0 where λ = λ(X 3 , t). The swelling ratio of the gel can then be expressed as J = λ 3 0 λ. As the swelling in the x 3 -direction is free, the normalized Cauchy stress componentσ 33 must be zero throughout the gel. From Equation (4) we then get Intermediate configuration Current configuration X x=φ(X) Equation 10 yields an expression for the normalized chemical potential that varies through the thickness of the gel asμ From the condition of molecular incompressibility, we can express the nominal concentration of solvent molecules as By combining Equations (5) and (8) the nominal flux in the x 3 -direction can be written as where we can express ∂μ ∂X3 through use of the chain rule ∂μ The conservation of solvent molecules requires that where ∂C ∂t can be obtained by the use of Equation (12 Finally, by inserting Equation (14) into Equation (13), we can with the use of Equation (16) write Equation (15) as The total out-of-plane swelling ratio (h/H) obtained through time solving Equation (17) with the method of lines is compared to the FEM solution using the implemented Fortran code in Figure 5, yielding a good validation of the Fortran implementation of the constitutive model.

Boundary conditions
To induce swelling in the hydrogel plate described in Section 2, the chemical potential is changed at the upper surface fromμ 1 toμ = 0. To avoid an abrupt change in boundary condition in the finite element simulation, the chemical potential is changed during a smooth step with a duration of 0.2 seconds (assumed to be short compared to the time to buckling initiation for transient swelling). The potential at the exposed surface is then kept constant at zero and the hydrogel swelling is driven by the diffusion process. Through this process, the plate is constrained from in-plane deformation at all lateral faces, while the lower surface is constrained from out-of-plane deformation. To comply with the LPA results and also accommodate for large swelling ratios, we want to mimic a nearly dry state at the start of the simulation, and hence set the initial homogeneous chemical potential toμ 1 = −2, giving an initial free swelling of J 0 = λ 3 0 ≈ 1.03.

Material parameters
The material parameters used in this study are given in Table 1. N v s defines the normalized stiffness of the substrate, while the parameter n defines the ratio between the stiffness in the film N v f and the stiffness of the substrate N v s according to n = N v f /N v s . For the diffusion coefficient D, a value in the range between 10 −9 and 10 −11 m 2 /s is reported as representative for the diffusion of water molecules in a gel [52,51]. However, we perform simulations for every decade of D in the range from 10 −11 to 1 m 2 /s to study the change in the buckling behavior of hydrogels as we gradually go from a transient to an equilibrium swelling process. The set of parameters as given in Table 1 leads to a total of 120 simulations to cover the parameter space.
In addition, simulations were performed for n = 0.1 and n = 0.8. These simulations confirm the general trends for soft-on-hard systems as presented in the following. Hence, the results for these two stiffness ratios are not included in the current presentation to enhance readability.

Buckling initiation
To obtain a measure for the global swelling ratio at the point of buckling and beyond, a plane is optimized with respect to the position of the uppermost layer of nodes in the film as illustrated in Figure 6. The height of the plane, h, is found by minimizing the residual r given as the sum of the squared distances between the nodes in their current position and the plane, i.e. r = n i=1 (∆ i ) 2 . The minimization procedure is performed using the SciPy package of Python. The global swelling ratio is defined as the plane height h divided by the initial plate thickness H = 0.5 mm.
There are multiple approaches that could be used to quantify the point where buckling initiates in a finite element simulation, e.g. abrupt changes in stress, strain energy or geometry. For practical applications, we consider the topology of the surface to be the most relevant measure. Hence, we define the point of onset of buckling by the use of a threshold value for the difference between the highest and lowest point on the surface, denoted ∆x 3 , i.e. buckling occurs if ∆x 3 ≥ ∆x crit 3 . Obviously, the obtained results for the onset of buckling will be influenced by the choice of the buckling limit, and the value to be used for ∆x crit 3 should depend on the application of the simulated gel. For the purpose of this study, we choose to use ∆x crit 3 =1 µm. Although the absolute values for the swelling ratio and time to initiation of buckling would be altered with a different threshold value, all trends presented in Section 6 would be preserved if the buckling limit would be set to a lower value like 0.5 µm.

Nodes
Optimized plane

Initial imperfection
To accommodate buckling, an initial imperfection is introduced in the upper surface of the stiff film. A Python script is used to change the x 3 position of all the surface nodes in the Abaqus input file. To ensure that the method of introducing this imperfection is not dominating the obtained results, we use two different modes of the initial imperfection. First, we introduce a single sinusoidal wave over the width of the gel where x i 3 is the new x 3 -position of the node i, H is the x 3 node coordinate in the perfect model, is the maximum perturbation given as a fraction of the initial height, and x i 1 is the x 1 -position of node i. Second, we use a random perturbation given by where p ∈ [−1, 1] is a pseudorandom number. While the expression in Equation (18) gives a smooth surface, it assumes a specific shape of the surface of the gel. The random imperfection in Equation (19), on the other hand, avoids a priori assumptions of the shape of the initial surface but can produce highly uneven surfaces.
The effect of the initial imperfection size, , on the maximum height difference in the plate surface, ∆x 3 , for slow and fast diffusion processes, D = 10 −11 and D = 1 m 2 /s, respectively, is shown in Figure 7. Clearly, a larger imperfection size would lead to faster growth of the height difference in the plane, however, the same asymptote is reached for all simulations where buckling is initiated. It can be seen that for the random initial imperfections, there is a significant increase in ∆x 3 at the beginning of the simulations with a slow diffusion process. This effect arises as the random imperfection introduces peak and valley nodes in the mesh, where the peak nodes would be exposed to more solvent and hence swell faster than the valley points (similar to free swelling of a cube where the corners would swell faster than the center of the faces [39]). This initial effect is not seen for D = 1 m 2 /s, as this would resemble an equilibrium process with equal swelling ratios in both peak and valley points, nor for the sinusoidal imperfection as this surface will be initially smooth. To demonstrate the effect of the randomness introduced in the mesh by Equation (19), results from three subsequent random imperfections applied to an initially perfect mesh using = 0.01% are shown (i.e. overlapping blue curves in Figure 7). Finally, the plot also shows that the imperfection must be larger than a critical size to trigger buckling, as seen for the curve for a random imperfection with = 0.001% where no buckling is initiated.
In all further calculations, we set to be 0.01%, meaning that the plate height at each node will be in the range between 0.49995 mm and 0.50005 mm and that we initially have ∆x 3 ≤ 0.1 µm. It is important to note that the initial imperfection must be smaller than the threshold value used to quantify the onset of buckling. In the following, we perform simulations using either the harmonic imperfection (Equation (18)) or the random imperfection (Equation (19)). Buckling limit D = 1 m 2 /s D = 10 11 m 2 /s = 0.05% rand = 0.025% rand = 0.01% rand = 0.01% sin = 0.001% rand Figure 7: Effect of the size of the initial imperfection for the evolution of ∆x 3 vs swelling ratio, using N vs = 0.01 and n = 5.

Plane strain assumption
To capture both the gradient of the chemical potential through the thickness of the plate in the transient swelling analyses and the evolving in-plane buckling pattern requires a relatively fine mesh, putting significant demands on computational resources.
However, it can be noted that for the problem at hand, a plane strain assumption would be correct up to the point of buckling. After buckling has initiated though, a plane strain assumption would restrain the model to a 1D buckling pattern, in contradiction with experimental and theoretical studies [59] having shown that other buckling patterns, like a hexagonal or herringbone mode, would be energetically favorable.
To investigate if a 2D plane strain assumption would predict buckling at the same swelling ratio as a full 3D model, despite the difference in bucking mode, a simple case of equilibrium swelling with N v s = 0.001 and n = 100 was used. Due to the homogeneous chemical potential and the large wavelength of the resulting buckling pattern, a relatively coarse mesh could be used. The 2D and 3D problems were both discretized in the same manner through the thickness using fully integrated higher order elements. The obtained buckling patterns and the increase in ∆x 3 during swelling are shown in Figures 8 and 9 respectively. The modes shown in Figure 8 are obtained from the first increment after initiation of buckling according to the threshold level ∆x 3 ≥ 1 µm. The initial buckling pattern in the 3D simulation can be seen to have a checkerboard pattern in the central region of the plate, this gradually expands to cover the plate before it would transition to a herringbone mode at larger swelling ratios [59]. For the plane strain simulations, the 2D nature of the model forces the plate to buckle in the less advantageous 1D buckling pattern more commonly observed when one of the in-plane stresses dominates [60]. From the fringe plots, it can be seen that the value of ∆x 3 is larger in the 3D simulation. However, from Figure 9 it is seen that the swelling ratio (i.e. the height of the optimized plane) at the onset of buckling is similar between the 2D and 3D simulations. Hence, to reduce the computational demands, a plane strain model is used in the further to study the point of onset of instability. This limits the study from considering detailed post-buckling analysis as the correct buckling pattern cannot be obtained.

3D simulation
Plane strain simulation Figure 8: Illustration of buckling mode for an equilibrium simulation with N vs = 0.001 and n = 100 obtained in a 3D and a 2D simulation. The fringe plot gives the out-of-plane deformation. The 2D result is extruded in the x 2 -direction to visualize its 3D resemblance.

Element type and size
For the further 2D simulations, the plate is discretized by fully integrated linear plane strain temperaturedisplacement elements (denoted CPE4T in Abaqus) where the temperature field is used to mimic the distribution of the normalized chemical potential through the plate (see Appendix Appendix A for the analogy between gel diffusion and heat transfer). Linear elements were preferred over second-order elements as the first use a lumped heat capacity matrix and hence avoid spurious oscillations in the temperature field during changes in the boundary conditions [61]. Thereby, with our use of the heat transfer and gel swelling analogy, oscillations in the chemical potential during transient swelling is avoided. See [41,62] and the references therein for a further discussion on oscillations in the chemical potential during finite element simulations of transient gel swelling.
The effect of the element size on the obtained results was studied by running simulations with N v s = 0.01, n = 2 and D = 10 −11 m 2 /s using between 1 and 30 square elements over the film thickness. To avoid a change in the surface topology as the mesh is refined, the harmonic imperfection (Equation (18)) was used. As the gradient in the chemical potential is largest close to the surface of the gel, the total size of the model Onset of buckling 2D model 3D model Figure 9: Comparison of the predicted onset of buckling for a 2D and a 3D model during equilibrium swelling.
was reduced by gradually increasing the mesh size going from the film and into the substrate as illustrated in the enlarged view in Figure 10. The resulting swelling ratio at the onset of instability as a function of the number of elements over the thickness of the stiff film is shown in Figure 11. The red squares refer to the swelling ratio, while the blue circles refer to the relative change in the swelling ratio compared with the previous (larger) element size. The same trend is observed with respect to the time to onset of instability. As the swelling ratio at the onset of buckling seems to be converged for a mesh with 30 square elements over the film thickness, this mesh size, leading to a total of 194 400 2D elements, is used in all further calculations. It is worth noting that for this element size, the random imperfection with set to 0.01% produces a maximum perturbation of the upper nodes of 3% of the initial element height. A sketch of the final mesh with a random imperfection can be seen in Figure 10.  can be observed between the two analysis methods. The predicted wavelength, on the other hand, shows a perfect correspondence between the two methods for n = 2 and n = 100, for the two other stiffness ratios a reasonably good correspondence is seen. A reason for the discrepancy in the wavelength results can stem from the finite width of the FEM model, while the LPA method assumes an infinitely wide plate. However, the generally good correspondence between the FEM and LPA results indicates that the length of the plate L in the FEM model is sufficient to yield results representative for infinitely wide plates. It can be noted that a comparison between FEM and LPA results using N v s = 0.01 shows a similar correspondence, and were omitted from Figure 12 merely to enhance readability.
For the soft-on-hard system (n = 0.5) both the LPA and the FEM approaches predict creasing as the initial instability mode. For the swelling ratio at the onset of buckling, on the other hand, the FEM and LPA results deviate with about 8%. However, it is observed that for n = 0.5 and N v s = 0.001 the FEM results have not fully reached the equilibrium plateau level at the diffusion coefficient of 1 m 2 /s (seen in Section 6.3, Figure 14), explaining a part of the discrepancy.

Stress and chemical potential profile at the initiation of buckling
The distribution of the normalized chemical potential,μ, and the normalized in-plane stress,σ 11 , through the thickness of the plate at the onset of buckling are shown in Figure 13 for all values of the diffusion coefficient using stiffness parameters of N v s = 0.01 and n = 2 (note that the presented trends are representative also for other parameter combinations). Here, X 3 is used to denote the x 3 -position in the plate in the intermediate state (i.e. the starting point for the simulations) such that the curves are easily comparable. The data are extracted from the right edge of the plane strain plate, however, the profiles would be representative for the whole plate at the initiation of buckling.
For the normalized chemical potential shown in Figure 13a, three distinct regions of the diffusion coefficient can be found. First, for all simulations with D ≥ 10 −2 m 2 /s, we see that the chemical potential at the onset of buckling is nearly homogeneous through the thickness of the plate and we have overlapping curves, i.e. equilibrium swelling is represented. We denote this region I. Second, we have a transition region in the range 10 −6 m 2 /s ≤ D ≤ 10 −3 m 2 /s where the chemical potential through the plate is increasingly inhomogeneous as the diffusivity is reduced. We denote this region II. Finally, the results for a diffusion coefficient D ≤ 10 −7 m 2 /s produces overlapping curves as the profile of the chemical potential through the thickness at the onset of buckling is nearly independent of the diffusion coefficient. We denote this region III. For the profiles of the normalized in-plane stresses at the onset of buckling seen in Figure 13b, we see the same three regions as for the normalized chemical profile, with overlapping curves for fast (I) and slow (III) diffusion and a transition zone in-between (II). For the equilibrium swelling simulations (i.e. D ≥ 10 −2 m 2 /s), the chemical potential is homogeneous through the thickness of the plate, and consequentially the compressive stress changes through the plate only due to the stiffness difference between the film and the substrate.
For the slow diffusion processes (i.e. D ≤ 10 −7 m 2 /s), on the other hand, the compressive stress is gradually reduced (i.e. becomes less negative) going downwards through the film thickness. Then there is a discontinuous reduction as one goes from the stiff film to the soft substrate. In the substrate, there is a further gradual decrease towards zero as one approaches the bottom of the plate. It can also be seen that the maximum compressive stress, located at the top of the film, increase as the diffusion coefficient is reduced. Figure 14 shows the swelling ratio of the plate (i.e. the global amount of swelling h/H) at the onset of buckling for values of the diffusion coefficient D in the range from 10 −11 to 1 m 2 /s, n between 100 and 0.5, and a normalized stiffness of the substrate of 0.001 or 0.01. The plotted data were obtained using a random initial imperfection in the upper surface of the plate (i.e. Equation 19), however, a similar plot would be produced if a harmonic initial imperfection was used (i.e. Equation 18). In accordance with the results for the normalized chemical potential and compressive stress profiles at the onset of buckling (discussed in Section 6.2) the swelling ratio at the onset of buckling can be divided into the three distinct regions I, II, and III. In region I, there is a relatively large and stable swelling ratio at buckling initiation that equals the equilibrium solution. Then, in region II there is a transition zone where the swelling ratio at buckling initiation is gradually reduced as the diffusivity is lowered. Finally, in region III a relatively low and stable swelling ratio at buckling initiation is obtained. In the plot, it can also be seen that the range of the diffusion coefficient in which the three regions occur depends on the value of n, with the transitions between the regions moving towards higher values of D as n is reduced.

Swelling ratio at buckling initiation
Further, the swelling ratio at the onset of instability can also be seen to depend on the absolute stiffness of the substrate and the film, and not only the ratio between the two. For the soft-on-hard system (n = 0.5), a larger critical swelling ratio is obtained at high values of D for the softer gels (i.e. N v s = 0.001) compared to the stiffer gels (i.e. N v s = 0.01). For hard-on-soft gels (n > 1) on the other hand, a larger swelling ratio at the onset of buckling is obtained for the stiffer gels (i.e. N v s = 0.01) compared to the softer gels (i.e. N v s = 0.001). This is especially seen in the results for n = 2 combined with low values of D, where a swelling ratio close to unity is present at the onset of buckling if N v s = 0.001, considerably below a swelling ratio of 1.2 as found for N v s = 0.01. A further discussion on the large difference between these parameter combinations is given in Section 6.5 discussing the obtained buckling profiles.  Figure 15 shows the time to the onset of instability on a log-log plot for the simulations where wrinkling was obtained as the mode of instability (see Section 6.5 for a discussion on buckling profiles). Simulations predicting creasing were excluded from the plot as the time to onset of creasing is known to show a strong mesh dependence [41].

Time to buckling initiation
In Figure 15, it is seen that there is a negligible effect on the time to the onset of buckling in region I and II (i.e. D ≥ 10 −6 m 2 /s), however, reducing the diffusion coefficient even further (i.e. into the physically reasonable region for water molecules diffusing in a gel) gives a clear rise in the time needed for buckling to occur, with an increasing effect as n is reduced for hard-on-soft gels. For all parameter combinations, a linear log-log relation between the time to instability and the diffusion coefficient can be seen in region III, corresponding to the stable region for swelling at the onset of instability for low values of D as shown in Figure 14. It is also seen that there is a clear dependence on the absolute value of the stiffness in the plate and the substrate and not only the ratio between the two, with longer times to buckling initiation for stiffer gels.
From the time to instability results it can be noted that the time used on the smooth increase of the chemical potential at the surface of the gel in the FEM model, i.e. 0.2 seconds, is negligible compared to the time to the onset of buckling for low values of the diffusivity. In addition, with the timescale needed to obtain buckling in region III, it is assumed that the contribution from viscoelastic effects in the polymer network would be negligible. For region I and II, on the other hand, further studies are needed to investigate how viscous effects might alter the initiation of buckling.

Buckling profiles
To study the dependence on the diffusion coefficient for the obtained buckling pattern we extract the wavelength found in the FEM simulations for all values of D. Figure 16 shows the obtained wavelengths for n = 0.5, n = 2 and n = 5. For n = 20 and n = 100, no clear trend in the critical wavelength could be seen.
For the soft-on-hard gel (n = 0.5), a critical wavelength of zero was obtained for all values of the diffusion coefficient as creasing was found to be the first mode of instability. This result was also found in the simulations with n = 0.1 and n = 0.8. For the hard-on-soft gels on the other hand (n = 2 and n = 5), wrinkling was found as the first mode of instability in most of the simulations, with a tendency of a reduction in the critical wavelength Λ c as the diffusion coefficient D was reduced. For n = 2 and N v s = 0.001, the critical wavelength drops to zero when the diffusivity goes below 10 −5 m 2 /s, as creasing (rather than wrinkling) was found as the critical mode of instability in these simulations. This is illustrated in Figure 17 showing the instability mode for the fast and the slow diffusion processes for N v s = 0.001 and n = 2. It can be seen that for fast swelling, a structured wrinkling pattern is obtained where the size of the initial imperfection is negligible compared to the length scale of the wrinkles. For the slow swelling simulation, on the other hand, it can be noted that the out-of-plane deformation is limited to the upper part of the film and that a localized folding process has started at the surface of the gel. The location of the developing crease depends on the initial surface imperfection in the model. The results shown in Figures  16 and 17 were obtained with a random imperfection in the upper surface of the plate (i.e. Equation (19)), however, similar results were obtained using the harmonic imperfection (i.e. Equation (18)). In addition, the shift from wrinkling to creasing as the diffusion coefficient decrease from 10 −5 m 2 /s to 10 −6 m 2 /s was found to be insensitive of the mesh size (tested in the range between 20 and 50 elements over the film thickness).
The wrinkling to creasing transition for a hard-on-soft gel clearly demonstrates how the time-dependent nature of swelling can cause a discrepancy between the buckling mode obtained in experiments compared to that predicted by an equilibrium analysis. Further, the fact that N v s = 0.001 and n = 2 reach instability by creasing before wrinkling at low values of D explains the large difference in the critical swelling ratio between N v s = 0.001 and N v s = 0.01 during slow swelling for n = 2 seen in Figure 14.
The authors hypothesize that creasing can occur during transient swelling of hard-on-soft gels with low film to substrate stiffness ratios due to a small difference between the strain energy in the creased and wrinkled configurations. For a sufficiently slow swelling process (i.e. a high ratio between the film thickness and the diffusion coefficient) in a sufficiently soft gel, instability will be triggered while the compressive stresses are focused in the upper part of the film, causing creasing to occur rather than wrinkling. However, further experimental and theoretical research is needed to fully understand the wrinkling to creasing transition during transient swelling for a hard-on-soft system and to quantify the critical conditions for this transition to occur. 6.6. Dependence on initial dimensions 6.6.1. Preliminaries For the length of the plate L, the generally good agreement between the LPA results, assuming an infinitely wide plate, and the FEM results, giving the plate a finite length, indicates that the ratio H/L is small enough for the results herein to be representative for infinitely wide plates. To discuss the influence on the obtained results by the chosen values for the film vs plate ratio η and total plate thickness H, we will consider equilibrium and transient swelling separately.

Equilibrium swelling
For a discussion on the geometrical dependencies under equilibrium swelling, the LPA framework is an effective method. First, looking at the film vs plate thickness ratio η, the dependence on the critical swelling ratio and the normalized critical wavelengthΛ c = Λ c /H are shown in Figures 18a and 18b, respectively. The data in the plots were obtained using N v s = 0.001, although similar trends could be obtained with N v s = 0.01. Starting with the soft-on-hard case (n = 0.5), both the critical swelling ratio and the normalized critical wavelength can be seen to be nearly independent of η. For the hard-on-soft cases (n > 1) on the other hand, we see that near the film to plate height ratio of η = 0.1 used herein the critical swelling ratio (Figure 18a) is relatively stable with respect to the film fraction. However, there is a clear tendency of an increasing critical swelling ratio as η gets larger than 0.25, with a stronger dependence observed for lower stiffness ratios. The predicted normalized critical wavelength (Figure 18b) for the hard-on-soft cases is changing significantly near η = 0.1, with increasing wavelengths as the film to plate thickness ratio is increased.
Second, for a change in the initial height of the plate (keeping the ratio η constant), both the critical swelling ratio and the normalized critical wavelength (as presented in Figure 18) were found to be independent of H under the assumption of equilibrium swelling of an infinitely wide plate.

Transient swelling
For a discussion on how the results for transient swelling depend on the initial thickness of the plate H (keeping η constant), it is worth noting that the defined problem introduces a relation between the length scale of the swelling plate and the diffusion coefficient of the gel material. This relation means that scaling D with a factor x 2 would be equivalent to scaling the dimensions of the plate with a factor 1/x. Hence, the results presented in the previous sections for varying diffusivity and constant dimensions could have been obtained with a constant diffusivity and varying the plate dimensions. This scaling is demonstrated in Figure 19, plotting the critical swelling ratio as a function of the initial plate thickness for a constant diffusivity of D = 10 −9 m 2 /s. It is important to note the difference in the dependence on the initial plate height H between equilibrium and transient swelling. This difference highlights the importance of the length scale of a swelling hydrogel system for whether simple equilibrium analyses would yield precise predictions for the initiation of buckling, with improved precision of equilibrium estimates for thinner gels. In addition, the length scale dependence in transient swelling means that the regions I, II, and III and their respective diffusivity values as given in Section 6.2 are specific for the plate dimensions used.
A further discussion on how a change in the ratios H/L and T /H would influence the initiation of buckling during transient swelling is considered out-of-scope for the present work.    Figure 19: Effect of the varying the initial plate thickness (retaining the shape of the plate) on the critical swelling ratio for a constant diffusivity of D = 10 −9 m 2 /s.

Conclusion
In this paper, we show that the onset of instability in hydrogels with gradient stiffness is highly influenced by the kinetic nature of the swelling process. Both the critical swelling ratio and the time to buckling initiation were found to increase as the diffusion coefficient of the material was reduced. In addition, these measures were found to depend on the absolute stiffness of the film and the substrate and not only the ratio between the two.
For the geometrical configuration studied herein, changing the diffusivity within the physical region for hydrogels in water 10 −9 − 10 −11 m 2 /s would have a negligible effect on the swelling ratio at the onset of buckling, while the time to the onset of buckling could change significantly, with an increasing effect as the stiffness ratio between the stiff film and the soft substrate is reduced.
For soft-on-hard gels, creasing was found as the first mode of instability independent of the diffusion coefficient. For hard-on-soft gels on the other hand, the diffusivity of the gel material could alter the buckling pattern by reducing the wavelength as the diffusivity was reduced. This is most evident for the combination of material parameters N v s = 0.001 and n = 2, where a stable wrinkling pattern was predicted by an equilibrium analysis, while creasing was triggered as the first mode of instability when using a physical value for the diffusion coefficient of the gel.
The results presented herein for a change in the diffusivity with constant plate dimensions could also be read as a change in plate dimensions (retaining the shape) with a constant value for the diffusivity. Whether simplified analyses based on an assumption of equilibrium swelling would be predictive for the onset of buckling in a layered hydrogel system would hence depend on the length scale of the problem.
For further work, we suggest that the transition from wrinkling to creasing in hard-on-soft gels is studied both theoretically and experimentally. Further, viscoelastic effects of the polymer network should be included in the constitutive modeling to study the influence this will have on the initiation of buckling. In addition, analyses of the post-buckling behavior of confined swelling gels where the three-dimensional nature of the buckling pattern is accounted for could increase the understanding of the fundamental mechanisms governing swelling induced buckling and pattern evolution.

Conflicts of interest
There are no conflicts to declare.
To find an analytical expression for ∂C ∂μ to implement as the specific heat capacity is non-trivial. A strategy suggested by Toh et al. [37] is to derive an expression forμ from Equation (4), using Equation (2) to relate J and Cμ = ln vC 1 + vC where trσ denotes the trace of the normalized Cauchy stress tensor andĪ 1 = (1 + vC) −2/3 I 1 . Assuming C = C(μ), the right-hand side of Equation (A.4) is differentiated with respect toμ by use of symbolic differentiation. Using ∂μ/∂μ = 1, an expression for ∂C ∂μ is obtained as Note that Equation (A.5) is slightly different from the expression given in Toh et al. [37]. The good benchmark results shown in Figure 5 indicates that Equation (A.5) yields good accuracy for the problem studied herein.