Parameter identification and uncertainty propagation of hydrogel coupled diffusion-deformation using POD-based reduced-order modeling

This study explores reduced-order modeling for analyzing the time-dependent diffusion-deformation of hydrogels. The full-order model describing hydrogel transient behavior consists of a coupled system of partial differential equations in which chemical potential and displacements are coupled. This system is formulated in a monolithic fashion and solved using the finite element method. We employ proper orthogonal decomposition as a model order reduction approach. The reduced-order model performance is tested through a benchmark problem on hydrogel swelling and a case study simulating co-axial printing. Then, we embed the reduced-order model into an optimization loop to efficiently identify the coupled problem's material parameters using full-field data. Finally, a study is conducted on the uncertainty propagation of the material parameter.


Introduction
Hydrogels, known for their biocompatibility, high water content, and adaptable mechanical properties, are at the forefront of innovation in the fabrication of functional materials for biomedical applications (Zhang and Khademhosseini, 2017).The fabrication of hydrogels has been supported by the integration of quantitative experimental techniques with computational models describing the material and mechanical properties of the fabricated hydrogel.A comprehensive review of available experimental methods to investigate hydrogel properties, along with an overview of some common mathematical modeling approaches describing hydrogels' mechanics and transport properties, can be found in Caccavo et al (2018).A recent review of the theories capturing hydrogel's diffusion-deformation coupled mechanisms can be found in Urrea-Quintero et al (2024).
Diffusion-deformation models for hydrogels combine continuum mechanics and thermodynamics with solvent diffusion physics through the balance of mass and Fick's law of diffusion (Chester and Anand, 2010;Bouklas et al, 2015;Liu et al, 2016), but lack an analytical solution in general (Chester et al, 2015;Urrea-Quintero et al, 2024).Hence, the model's solution requires numerical methods, such as the Finite Element Method (FEM) (Anand and Govindjee, 2020;Urrea-Quintero et al, 2024).Running such a fullorder model (FOM) finite element simulation can become costly (Bonatti et al, 2021;Chirianni et al, 2024).This has hindered the widespread adoption of multiphysics modeling in hydrogel manufacturing, where many-query tasks are often performed, and several hundreds or thousands of similar runs are often necessary.Typical many-query tasks are material parameter identification (Zou et al, 2018), uncertainty quantification (Pagani and Manzoni, 2021), process optimization, model-based 3D printing monitoring, or feedback printing error compensation technologies.(Lassila et al, 2014a;Strazzullo et al, 2020).
Model-order reduction (MOR) offers a promising solution to the challenges posed by the computational burden of multiphysics models.By using a ROM, we can save on computation resources while still capturing the dominant behaviors of the system.ROMs can be tailored to specific regions of interest, ensuring accuracy where it is most crucial.ROMs can bridge the gap between the detailed insights provided by multiphysics models and the practical necessities of the many-query tasks mentioned above; see, for example, Chen et al (2015), Strazzullo et al (2020), and Pagani and Manzoni (2021).
In the MOR literature, Proper Orthogonal Decomposition (POD) and its variants are particularly popular across multiple disciplines; see, for example, Lu et al (2019) for a review of POD-based methods for MOR in a variety of research areas.Additionally, we refer to Benner et al (2015aBenner et al ( , 2020) ) and Gräßle and Hinze (2018) for a comprehensive introduction to POD applied to MOR and Benner et al (2015b) for a survey in the context of parametrized PDEs.In particular, POD-based reduced-order modeling has been successfully applied in engineering, in the context of fluid mechanics (Girfoglio et al, 2022;Ballarin et al, 2015;Nonino et al, 2021), image processing (Vaccaro, 1991), rotor dynamics (Lu et al, 2015a(Lu et al, , 2016)), for structural damage detection (Eftekhar Azam et al, 2017), and in the area of structural dynamics (Kerschen et al, 2005).POD-based methods can be also used to solve optimal control problems (Negri et al, 2013;Brunton and Noack, 2015;Ballarin et al, 2022;Volkwein, 2001).Moreover, POD-based ROMs can be physically interpreted based on the concept of proper orthogonal modes (POMs) (Kerschen et al, 2005).However, while the theory is solid for linear and some non-linear systems (Benner et al, 2020), applying ROMs to tasks such as hydrogels fabrication by printing is largely unexplored.In this work, we exploit the capabilities of POD to predict the transient diffusion-deformation of hydrogels in co-axial printing and its efficiency in tasks many-query such as material parameter identification and uncertainty propagation.
Different POD-based algorithms have been proposed in the literature.Global POD methods focus on the estimation of global POD basis functions (Benner et al, 2015a(Benner et al, , 2020)).In contrast, local POD methods are based on snapshot clustering (Sahyoun and Djouadi, 2013).Adaptive POD methods have been developed to improve POD basis functions construction by choosing additional snapshot locations in an optimal way (Lass and Volkwein, 2014).Incremental POD has also been proposed, allowing an adaptive enrichment of the POD basis in case of unforeseen changes in the solution behavior (Fareed et al, 2018;Fischer et al, 2024).Transient POD methods focus on snapshots of time-dependent problems (Lu et al, 2015b).The so-called nested-POD method applies to time-dependent partial differential equations (PDEs) and enables a two-step dimensionality reduction: first in the temporal domain and then in the spatial domain, which also includes parameter information (Kadeethum et al, 2021(Kadeethum et al, , 2022)).This approach has been shown to be particularly useful when the PDE requires many steps to be approximated.In this work, we adopt the nested-POD approach and compare it to the traditional POD.
In this paper, we consider a specialized theory for small deformations superposed on a previously homogeneously swollen gel as presented by Anand and Govindjee (2020).This model is derived following rigorous thermodynamic concepts and can be considered the linearized version of the more general diffusion-deformation theory presented by Chester and Anand (2010) and Chester et al (2015).This diffusion-deformation model was recently validated by Chen et al (2020) and showed good agreement with experimental data to study the diffusion-deformation process within the first three hours.If long-time deformations are present, a nonlinear material model accounting for the hyperelastic deformation of the polymer network and describing the mixing of the solvent with the polymer network should be considered (Chester and Anand, 2010;Bouklas et al, 2015;Chester et al, 2015;Urrea-Quintero et al, 2024).We implement two POD-based ROMs, study their accuracy, and evaluate the computation speed-up to assess its suitability for many-query scenarios.
The main contributions of the paper can be summarized as follows: 1. We present the first work that considers POD-based parametric ROMs for the diffusion-deformation of hydrogels.2. We demonstrate the effectiveness of the obtained parametric ROMs for material parameter identification and uncertainty propagation for hydrogel printing processes.
The paper is organized as follows: Section 2 summarizes the coupled diffusion-deformation theory of gels in the small deformation regime.Section 3 introduces the POD framework.Section 4 states the material parameter identification problem enabled by the reduced-order model.In Section 5, several simulations are conducted and computationally analyzed.Concluding remarks are given in Section 6.

Small deformation theory for the diffusion-deformation of hydrogels
Diffusion-deformation theories of gels integrate fluid transport and solid mechanics principles to explain how gels swell or shrink due to fluid movement and how these processes influence their shape and mechanical properties (Chester and Anand, 2010;Bouklas et al, 2015;Chester et al, 2015;Urrea-Quintero et al, 2024).In a diffusion-deformation model of gels, the driving forces are fluid diffusion into or out of the gel and mechanical deformation due to external loads or internal stresses.The model should allow for predicting swelling or shrinking behaviors, stress distributions, and changes in mechanical properties over time.The primary outcomes of this model include the gel's deformation patterns and mechanical responses under varying conditions.Essential parameters for model calibration encompass the diffusion coefficient, the gel matrix's mechanical properties such as shear and bulk modulus, and fluid-gel interaction parameters such as Flory's interaction parameter.
This section summarizes the theory describing the diffusion-deformation mechanisms in elastomeric gels that undergo small deformations under isothermal conditions.The model is derived as a special case of a more general theory of chemoelasticity (Anand, 2015;Anand and Govindjee, 2020).This linear theory is suitable for capturing swelling and drying for the first few hours of the diffusion-deformation process (Anand and Govindjee, 2020;Chen et al, 2020).For example, the theory can describe the diffusiondeformation of hydrogels in bioprinting applications where the hydrogel is mixed in the printer's nozzle with living cells and nutrients necessary for cell growth and proliferation, or in co-axial printing where the hydrogel is directly crosslinked in the printer's nozzle (Kjar et al, 2021).
As notation rules, we denote gradient grad(•) and divergence div(•).The time derivative of any field is denoted by ∂ t (•).The operator tr(A) refers to the trace of the second-order tensor A. We denote the spatial dimension with d, and in this section, we exclusively work with d = 3.Finally, let I := {0, T f } be the time interval with end time value T f > 0. (1)

Kinematics of the deformation
For small deformations, the strain tensor can be approximated by the linearized strain tensor, ε : which represents the symmetric part of the displacement gradient.

Governing partial differential equations
The two governing PDEs for the quasi-static vector-valued displacements u : B → R d and scalar-valued, time-dependent, fluid content in terms of the concentration c : B × I → R consist of: 1.The local form of the balance of linear momentum: (3) Here, σ : B → R d×d denotes the Cauchy stress tensor, specified below, and σ 0 : B → R d×d its initial value at t = 0. External actions consist of body forces per unit deformed volume b : B → R d .Moreover, boundary conditions prescribe displacements ū : B → R d and traction t : ∂Bt → R d on separate portions of the boundary.Notice that inertial effects have been neglected due to the considerably slow dynamics of the fluid diffusion evolution w.r.t. the time scale of the wave propagation.2. The local form of the mass balance of fluid content inside the hydrogel: Boundary conditions are defined by a Robin condition on the boundary ∂B R , which prescribes a relationship between the fluid flux j and the chemical potential µ.Specifically, α R is a proportionality constant, µ is the chemical potential at the boundary, and µ ∞ : ∂B R × I → R represents a reference chemical potential.This condition models the proportional flux response to the difference in chemical potential at the boundary and is considered to be more consistent with experimental observations of fluid absorption in hydrogels (Chen et al, 2020).Furthermore, µ 0 : B × {t = 0} → R refers to the initial value of the chemical potential inside the hydrogel.Notice that the mass balance of fluid content is written in terms of c and j, but its corresponding boundary and initial conditions involve the chemical potential µ.The connection of equation ( 4) with µ becomes clear by introducing a constitutive relation for j, for example, through Fick's laws of diffusion.

Thermodynamically consistent constitutive theory for linear isotropic gels
In this subsection, we present a thermodynamically consistent constitutive theory describing the coupled diffusion-deformation mechanisms for linear isotropic gels.Only the main ingredients are provided, but the step-by-step derivation of the constitutive model is detailed in A.
For linear isotropic hydrogels, the following thermodynamically-consistent constitutive relation can be established for Cauchy's stress with G representing the elastic shear modulus, λ d the drain Lamé parameter, 1 the identity tensor, χ the stress-chemical modulus, and Λ the chemical modulus.
Additionally, the following constitutive relation can be established for the species concentration as a function of the chemical potential Moreover, we append a thermodynamically motivated choice for the fluid flux j based on Fick's law for species diffusion.That is, j is proportional to the gradient of the chemical potential, namely, with M as the species mobility tensor defined as where D is a constant for fluid diffusivity, k B is the Boltzmann constant, and T is the absolute temperature.
The final model is characterized by the four independent material parameters, two from isotropic linear elasticity [G, λ d ] and the other two from the fluid equations [Ω, D], where Ω is the molar volume of the fluid as defined in A.
It is worth mentioning that the theory introduced here can be linked with the classical theories of chemoelasticity and linear poroelasticity; see Anand (2015) and Anand and Govindjee (2020) [Chapters 15 and 16] for more details.The related transformations are explained in Appendix A.

Coupled diffusion-deformation model in normalized form
To facilitate the analysis of the model derived in A, Problem 4, the physical quantities are normalized as follows In equation ( 9), l is the reference length of the hydrogel, λ * is the dimensionless first Lamé constant, which indicates how much volumetric versus shear deformation contributes to the material response, and A is the dimensionless scaling factor for the chemical potential's influence on the stress.This normalization closely follows the methodology presented by Chen et al (2020).
The reader should note that the normalized model only contains two non-dimensional parameters, compared to six parameters in the original formulation.Moreover, as the time-dependent chemical potential enters into the displacements, the governing PDEs become quasi-static, indirectly depending also on time.Thus, we write u : B × I → R d in the remainder of the paper.
The resulting coupled PDE model with normalized quantities reads: Problem 1 (Strong form linear diffusion-deformation of gels).Given ū, t * , µ * ∞ as boundary data and σ * 0 , µ * 0 as initial data, find u : with the constitutive law for stresses in normalized form given by and the strain tensor defined as in equation (2), which now contains only two dimensionless material parameters, namely [λ * , A], defined in equation (9).
The reader should notice that, for simplicity in the notation, the asterisks in the normalized quantities will be removed in the subsequent sections.

Weak formulation
The solution of the coupled PDE system given by equations ( 10) -( 17) consists of a vector-valued field of displacements u and a scalar-valued field of the chemical potential µ.
In the following, we derive the weak formulation.We adopt standard notation for the usual Lebesgue and Sobolev spaces, see, for example, Wloka (1987).The functional space H 1 (B) d is a Sobolev space that consists of functions defined on a bounded domain B ⊂ R d , with square integrable partial derivatives up to the first order.To this end, we define the trial and test spaces as follows To formulate both problem statements in an abstract fashion, we introduce for the displacement system the bilinear form a((µ, u))(v).Furthermore, let b(v) be the given right-hand side data.Next, for the balance of fluid concentration, we use c((µ, u))(q) and d(q).Then, the weak formulation can be written as Problem 2 (Weak form).Find (u, µ) ∈ V ū × Q, with µ(0) = µ 0 and σ(0) = σ 0 , such that for t ∈ I it holds where with " : " denoting the double contraction of the second-order tensors σ and grad(v), where σ(µ, u) is defined by equation (17), which couples the two balance equations.

Discrete weak formulation and numerical solution
In order to deal with efficient POD-based ROMs, we are required to fulfill the assumption of affine parametric dependence on the operators appearing in the weak form (18) at the FOM level.This allows us to decouple the construction stage of the reduced-order space (offline) from the parametric evaluation stage (online).
To meet this requirement, let us assume that the computational domain B is decomposed into triangular elements.We employ Taylor-Hood elements such that the inf-sup condition is fulfilled, that is, V h,0 × Q h ⊂ V 0 × Q, with quadratic functions in V h and linear functions in Q h .Then, the parametric semi-discrete variational monolithic formulation for the diffusion-deformation model reads: Problem 3 (Parametric spatial semi-discrete weak form).Given the material parameter θ where equation (23) 1 can be reformulated as Additionally, equation (23) 2 can be written analogously as Here, M((µ h , u h ), t; θ) is the semi-discrete operator that captures the transient behavior of the coupled diffusion-deformation problem and is affine parameterized by the material parameters θ = [λ, A].
The time-dependent operator defined in equation ( 28) 1 is further discretized in time with the firstorder A-stable implicit Euler scheme.The arising linear system of equations in Problem 3 is solved with a sparse direct solver.The implementation is done in FEniCS.All details of this FOM realization can be found in Urrea-Quintero et al (2024).

Parametric POD-based reduced order modeling
In this section, we describe parametric POD-based and parametric nested-POD approaches as reducedorder modeling strategies applied to the transient diffusion-deformation of hydrogels.Such a ROM construction allows to speed up the solution of the coupled model presented in Section 2, which enables its use in many-query scenarios.
The computation in the reduced-order modeling strategy is divided into two phases: offline and online.
(1) The offline phase refers to the parametric ROM construction.It mainly consists of two steps: (1a) Collect a set of snapshots that are FOM solutions for different parameter values and time instances.The simulation time interval for each parameter instance remains the same.(1b) Find an optimal reduced-order approximation of the primary variables (µ h , u h ) by projecting the snapshot matrix into a reduced linear subspace using Galerkin approach.Figure 1 illustrate the main steps of the offline phase for the parametric ROM construction.
(2) The online phase refers to the use of the ROM as a surrogate of the FOM to predict the diffusiondeformation of hydrogels in many-query tasks.Here, the parametric ROM is used for parameter identification; see Section 5.
The offline steps are explained in detail below.The starting point is the assembly of the snapshot matrix in Section 3.1, followed by the singular value decomposition (SVD) in Section 3.2.We introduce the nested-POD approach as a computationally efficient variant in Section 3.3.The Galerkin projection is outlined in Section 3.4.The implementation of the parametric POD and nested-POD-based ROMs is inspired by Lassila et al (2014b) and Kadeethum et al (2021).

Assembly of snapshots matrix
The snapshots are obtained by solving the model This results in a global snapshot matrix h the number of chemical potential degrees of freedom, iii.N t the number of time steps, and iv.N θ the number of parameter instances used for training.
Different assembly strategies arise from splitting the subset global snapshots into subsets in order to enhance computational efficiency.Here, we consider partitioning in the primary variables (µ h , u h ) as well as in space and time.The reader should note that the latter approach is also referred to as the nested-POD approach.
Partitioning of primary variables: To this end, we consider a partitioned approach to derive reduced bases for µ h and u h separately.We focus on the displacement field in what follows, although a similar procedure applies to the chemical potential.Given the parameter instances θ i , ∀i ∈ N θ of the training dataset, the corresponding temporal snapshot matrix S u (θ i ) ∈ R N u h ×Nt is assembled as: for the displacement field u h (., t n ; θ i ) at a specific time t n , for n = 0, . . ., N t , and parameter instance θ i .
The complete snapshots matrix S u ∈ R N u h ×(Nt•N θ ) , is then constructed by aggregating the individual temporal snapshot matrices S u (θ i ) as: Consequently, each row in S u corresponds to a degree of freedom in the spatial discretization, and each column corresponds to the state of the system at a particular time for a particular parameter sample.

POD for the displacement field
In this subsection, we describe the key parts of the POD for the displacement field.The chemical potential POD is computed in an analogous fashion.In synthesis, POD reduces the dimensionality of S u (32) by computing a matrix V ∈ R N u h ×r , which represents the best approximation in the least-squares sense to the snapshot matrix S u .Here, r ≪ N t • N θ is the number of retained modes based on analyzing the singular values that retain the most significant information in the snapshots matrix.We refer the reader to equation ( 34).
The POD solution is then derived by solving the optimization problem: where ∥•∥ F denotes the Frobenius norm, which measures the distance between S u and its approximation VV T S u in terms of all entries of the snapshots matrix.Moreover, σ i refers to the i-th singular value of S u and VV T yield the identity matrix with dimension of N u h × N u h .The singular values σ i provide quantitative guidance for choosing the size of the POD basis.That is, the number of modes r required to capture a desired fraction of energy of the FOM, measured by the fraction of retained energy where E total is the total energy of the system and E r is the retained energy by r modes.Consequently, the number of modes r is chosen such that the retained energy ratio meets a predefined target fraction

Nested-POD approach for the displacement field
For a large number of parameters N θ and/or time steps N t , applying SVD on the snapshot matrix (32) may become inefficient, requiring large memory storage.An alternative is to apply the nested-POD approach; see, for instance, Kadeethum et al (2021), and also our own study in Figure 5 (Section 5).
The snapshot matrix S u (32) can be further decomposed into temporal and space parts to apply the nested-POD approach.That is, the time dimension is reduced first, followed by the spatial dimension for the parameter samples.The dimensionality reduction steps can be summarized as follows: 1. Temporal dimensionality reduction: POD is applied individually for every parameter sample θ i on the snapshots matrix S u (θ i ) (31).In this way, we obtain the basis functions for the first n t ≪ N t temporal modes and for each θ i .The snapshot matrix after the temporal reduction is represented as Su Then, in analogy to the optimization problem (33), POD is applied again on Su to retain the r first modes and construct spatial bases functions.
For both the temporal and spatial dimensionality reduction, the criterion ( 34) is employed to select the number of required basis functions.
It is worth remarking here that the SVD algorithm is at the core of dimensionality reduction.In this work, we adopt the SVD algorithm provided by the RBniCS library.Other methods can also be used to improve the SVD efficiency further.For example, the method of snapshots (Sirovich, 1987) and the streaming method of snapshots (Hemati et al, 2014), or modern variants such as randomized (Martinsson et al, 2011) and incremental (Brand, 2006;Kühl et al, 2024) SVD.

Galerkin projections
In this subsection, we finally construct the reduced-order solutions using Galerkin projections.In both the POD and nested-POD cases, only the first r modes are retained during the dimensionality reduction and are employed to construct the basis functions for the reduced order space U r for the displacement field u.Thus, we map the temporal and spatial dimensions to their representation in a linear manifold, producing the best approximation in the reduced space U r for u(•; t, θ).
The optimal representation is obtained on a reduced basis space by means of Galerkin projection.Let (w u 1 , . . ., w u r ) denote the basis functions spanning U r .Given a time t n and a parameter instance θ i in the training set, the best approximation Ũh ( where the coefficients ξ U k are solutions to the Galerkin projection problem.Specifically, for given The operator (w u k , w u j ) Ur forms a matrix that can be precomputed and stored, facilitating the efficient computation of the coefficients ξ U k for each new instance of the parameters or time.The reduced basis approximation Ûh (•; t n , θ i ) is computed as: In the same fashion, we obtain the reduced basis approximation for the chemical potential: Here, z µ 1 , . . ., z µ r are the basis functions spanning the reduced order space Γ r .Furthermore, ξ Γ k are the coefficients resulting from the Galerkin projection.
We can numerically approximate the parametric discrete FOM M((µ h , u h ), t; θ) in Problem 3 by the POD-based ROM M(( Γh , Ûh ), t; θ) as to approximate µ h and u h in the online phase.
The reader should notice that, in the online phase, the ROM should be only used for predictions of µ h and u h within the parameter range for which it has been trained.Particularly, the ROM model in equation ( 41) needs to be evaluated at the new desired θ values.If, in the online phase, it becomes apparent that predictions outside this parameter range are necessary, additional training is recommended, as extrapolation is prone to errors.Then, one can look into various methods proposed, for instance, by Lange et al (2024), Irick and Brown (2019), or Fischer et al (2024).Furthermore, some research has also been focusing on using machine learning combined with ROM to tackle this issue, as done by Haasdonk et al (2023).However, in the present study, we did not consider such cases.
Finally, we utilize the L 1 , L 2 and L ∞ norms to evaluate the approximation errors, for example, e u := u F OM − u ROM between the FOM u h F OM and ROM u h ROM solutions.

Material parameter identification from full-field data
We exemplary consider material parameter identification from full-field observations to evaluate the ROM's performance in a many-query scenario.A recent review on parameter identification from full-field data in computational mechanics can be found in Römer et al (2024).Here, our objective is to identify the optimal set of dimensionless material parameters θ opt = [λ opt , A opt ] ∈ R 2 that minimize a metric L(θ) for the discrepancy between model predictions (parametric POD-based ROM or FOM) and some observational data capturing the diffusion-deformation process of the hydrogel.To this end, material parameter identification can be cast as the below optimization problem: The discrepancy metric L(θ) is often referred to as loss or objective function.In order to ensure proper convergence of the optimization problem stated above, the loss contributions from the different physical fields are computed and scaled individually.This leads to the following composite loss function whose contributions are computed from the relative L 2 norms of the difference between observed and computed displacement fields u obs and u(θ) as well as of the difference between observed and computed concentration fields µ obs and µ(θ): The solutions µ(θ) and u(θ) can be obtained either from the FOM predictions or its ROM approximation given by equation ( 41) at those discrete time steps t i at which observational data is available.The solution of the optimization problem stated in equation ( 42) requires iterative solution techniques, which makes the use of a ROM appealing.In this work, the L-BFGS-B optimizer is used, a variant of the Limited Memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm, which belongs to the family of quasi-Newton methods (Zhu et al, 1997).This algorithm is particularly suited for solving large-scale optimization problems with bound constraints.Its hyperparameters are provided in Section 5.1.4.

Numerical simulations
In this section, we conduct two extensive numerical studies to evaluate the effectiveness of parametric POD-based ROMs in capturing the main dynamics of the FOM describing the diffusion-deformation process of hydrogels.The motivation behind the numerical computations is to capture the step-by-step process of manufacturing hydrogels using bioprinting as functional materials in biomedical applications.Typically, hydrogel fabrication is done in two stages.The first stage involves characterizing the material to be manufactured, whereas, in the second stage, a bioprinter can be used to fabricate more intricate structures, possibly capturing complexities closer to those in living tissues.
Figure 2 illustrates the simulation setup of the two problems.First, a benchmark problem is considered, namely, a free-swelling of a 2D square hydrogel block (Figure 2a).We implement parametric POD and nested-POD-based ROMs, test their accuracy, compare the simulation speed-up with respect to the FOM, and perform a material parameter identification to evaluate the ROMs' effectiveness in this manyquery scenario.The second test is a case study mimicking co-axial hydrogel bioprinting (Figure 2b).For the case study, we present the results for the parametric nested-POD-based ROM and test its accuracy and computation speed-up with respect to the FOM.The code for the POD-based ROM is provided for interested readers.Finally, we conduct a parameter uncertainty propagation and exploit the simulation speed-up achieved by the nested-POD-based ROM.The FOM is solved using FEniCS1 (Alnaes et al, 2015) and the ROM is implemented in RBniCS2 (Rozza et al, 2024).All codes to reproduce and verify the results can be accessed in the following online repository: Agarwal et al (2024).

Benchmark: free-swelling of a 2D square hydrogel block
This subsection refers to the first stage of hydrogels manufacturing.The first stage involves characterizing the material to be manufactured.This is done through a series of experiments in the laboratory where molding techniques are used to fabricate simple hydrogel blocks under controlled conditions.These fabricated hydrogels are then subjected to, for example, free-swelling experiments to i. characterize material properties such as bulk and shear moduli and absorption capabilities of diffusing species and ii.validate mathematical models that can be used later as predictive tools for the material behavior.

Simulation setup
As a benchmark problem, we consider the transient free-swelling of a gel in a two-dimensional setting as illustrated in Figure 2a.The 2D hydrogel block initially has a square cross-section.The hydrogel block is immersed into a non-reactive solvent with a reference chemical potential µ 0 = 0.Only a quarter of the whole model is considered because of the symmetry of the block.In terms of mechanical boundary conditions, the bottom and left edges are subject to Dirichlet boundary conditions u y = 0 and u x = 0, respectively.The top and right edges are considered to be traction-free.In terms of chemical boundary conditions, the bottom and left edges are subject to zero fluid flux.The chemical potential is subject to a Robin boundary condition as defined in equation ( 15) on the top and right edges.A similar setup to this example can be found in, e,g., Chester and Anand (2011), Liu et al (2016), or Chen et al (2020).
The following values are selected as the nominal dimensionless material parameter 3 : λ = 1558 and A = 4000, and the Robin boundary coefficient is defined as α = 0.66.The chemical potential's initial condition is µ o = −0.3124,and the domain is considered to be initially stress-free.These material parameters taken from Chen et al (2020) represent a gelatin-based hydrogel.The reader is referred to Chen et al (2020) for more details on the gel preparation and its mechanical characterization.
The domain is discretized into triangular elements with quadratic interpolation for displacement and linear interpolation for the chemical potential to obtain the Taylor-Hood elements as introduced in Section 2. The dimensionless coupled system of equations ( 23) is solved as a variational monolithic problem using the FEniCS package as a PDE solver.We remark here that we solve the FOM in a monolithic fashion, but later for the POD, we partition the snapshot matrix into two parts in order to apply the partitioned approach as described in Section 3.
Following the convergence analysis reported in B, a mesh density of N h = 50, that is, 5000 triangular elements, is specified to obtain what we denote as the FOM nominal simulation setup.The total dimensionless simulation time is chosen as T = 0.25, and the time interval {0, T } is discretized into 100 time steps.The reader should notice that T = 0.25 corresponds to about t = 3 hours of the hydrogel diffusion-deformation before normalization based on equation ( 9).In their experimental work, Chen et al (2020) found that the linear model is valid up to 3 h of fluid absorption.For larger simulation times, a nonlinear theory should be adopted (Urrea-Quintero et al, 2024).

Numerical results for the full-order model
Figure 3 shows the primary variables for the FOM at three different simulation times.Figure 3a displays the evolution of the chemical potential within the 2D hydrogel block.It can be observed that the chemical potential transitions over time from its initial condition to its boundary value.This increase in chemical potential indicates that as the gel absorbs the solvent, its internal fluid concentration rises, driving the absorption process towards equilibrium.
Figure 3b illustrates the evolution of the displacement magnitude.It is evident that displacement increases in tandem with solvent absorption, with the highest displacement occurring in the corner fully exposed to the solvent.The symmetry in the displacement field reflects the isotropic swelling of the gel.Figure 4 shows the time evolution of normalized chemical potential (Figure 4a), displacement (Figure 4b), and normalized stress (Figure 4c) at different points in the 2D hydrogel block.The figure reveals that the four selected points undergo different diffusion-deformation pathways, highlighting the spatial dependency of the problem.Figure 4 complements observations from Figure 3. Figure 4c shows that while a region within the domain experiences compressive stresses, another region experiences the opposite.This stress gradient is the main driver of the fluid's diffusion.The figure also reaffirms that the points on the right side remain traction-free throughout the simulation.

Numerical results for the reduced-order model
In this section, we provide numerical experiments to validate and demonstrate the ROMs' accuracy and computational efficiency.All the following computations are conducted in RBniCS.The same mesh density and time steps as for the nominal simulation are taken.Material parameters are sampled within the ranges λ = [1000, 2000] and A = [2000,6000].
In order to minimize the overall computational effort, the immediate task is now to identify the minimum number of basis functions and samples.For this purpose, we follow the procedure presented in Section 3 and summarized in Figure 1.This procedure is followed for both the traditional parametric POD-based ROM and the nested variant, both of which have been introduced in Section 3.
Figure 5a depicts the generated training and testing samples from the material parameters used in the construction of the ROM.The samples were generated using Monte Carlo sampling from a uniform distribution.We noticed that 10 samples are the minimum to obtain an acceptable approximation of the FOM with both POD and nested-POD-based ROM.We trained the POD-based ROMs using 30 samples. Figure 5b reports the wall time required for training the ROMs adopting the POD and nested-POD algorithms for an increasing number of training samples.The wall time is composed of the time required to create the snapshots and to perform the SVD on the resultant snapshots matrix.
POD-based approach: Figure 6 shows the singular values and energy retained by the SVD on the snapshots matrix.From Figures 6a and c, it is observed that the singular values magnitude decays very fast for both chemical potential and displacement.Such a behavior is advantageous, as it facilitates the development of an efficient low-dimensional ROM.In fact, from Figures 6b and d, it is possible to conclude that only six basis functions are required to retain up to η = 99.9999% of the energy contained in the FOM.
An illustration of the error for the selected number of basis functions can be found in Figure 7. Here, the POD is assembled, and the predictions are compared to the nominal FOM model (A = 4000, λ = 1558, and α = 0.66).The reader should notice that these parameter values were not included in the training dataset; therefore, they can be used to validate the POD-based ROM.
Figure 7 depicts a comparison of the temporal evolution of chemical potential (Figure 7a) and displacement (Figure 7b) at different points inside the 2D hydrogel block for the POD-based ROM.The figures show that the root mean square error (RMSE) between the POD-based ROM and FOM has its maximum at the beginning of the simulation.However, for the chemical potential, the RMSE reaches 3 × 10 −3 in about T = 0.05 and keeps decreasing as time evolves.For the displacement, the RMSE is smaller than 10 −3 for almost the whole simulation time.
Nested-POD-based approach: Figures 8 and 9 report the singular values and fraction of energy retained for the chemical potential and displacement for the nested-POD approach, respectively.On the one hand, Figures 8a, b  temporal dimension, for which SVD is locally applied for each material parameter sample.The behavior of the singular values and retained energy are presented for three out of the 30 samples taken for the parameters.The figures show that the singular values and retained energy behavior are very similar for the temporal dimension in each case.From the analysis, it is concluded that only four temporal modes are necessary to retain the desired fraction of energy η.On the other hand, Figures 8c, d and 9c, d display the behavior of the singular values and retained energy for the spatial dimension, for which SVD is applied over the snapshots matrix containing the temporal basis functions and spatial information for the 30 training parameter samples.From this analysis, it is concluded that only six spatial basis functions are required to retain the desired fraction of energy η contained in the FOM for the chemical potential and displacement.
As we found that the results obtained with the nested-POD-based ROM were identical to those reported in Figure 7, we instead plot a full-field visualization of the simulation outcome for [λ, A] = [1558,4000].To this end, Figure 10 displays the chemical potential (Figure 10a) and displacement (Figure 10c) obtained with the six basis functions at three different time steps for the nested-POD-based ROM, namely, T = {0.02,0.15, 0.25}.Errors with respect to the FOM simulation results reported in Figure 3 are displayed for the chemical potential in Figure 10b and for the displacement in Figure 10d.As seen from Figure 10, the nested-POD-based ROM can predict the FOM to a satisfactory extent.From Figures 10b and d, it can be observed that the maximum discrepancy between the FOM and ROM is less than 0.55% for the chemical potential and 0.3% for the displacement field, respectively.It is also clear that the maximum discrepancy occurs at the beginning of the simulation but fades as time evolves.This result is expected because, at the beginning of the simulation, the chemical potential field has a very sharp gradient close to the exposed boundaries, but it is almost flat towards the center.It is well known that global surrogate models suffer when approximating flat functions.
POD and nested POD-based ROMs require the same number of basis functions to retain the desired fraction of energy η, and both approaches lead to the same level of accuracy; see also C for the convergence error analysis.However, the demand for computational resources significantly increases with the number of parameter samples.The advantage of the nested-POD algorithm becomes evident as reported in  Figure 11 demonstrates the significant computational speed-up achieved with the ROM.For instance, the FOM computation time for N h = 10 is 0.33s, but for N h = 160, it is already 213.4s.This corresponds to an exponential increase of the computational cost w.r.t. the mesh density.In contrast, the ROM computation time for N h = 10 is 0.018 s, but for N h = 160, it is still only 0.025 s.That is, the ROM simulation is about a factor of 8536 faster than the FOM simulation for N h = 160.
In the following, we restrict our analysis to the nested-POD approach only.This is because the same number of basis functions is required to retain the desired fraction of energy η, the same level of accuracy has been obtained with POD and nested POD-based ROMs, and given the computational advantage of nested-POD algorithm for a large number of training samples.However, the code is provided to replicate results with POD and nested-POD-based ROM for all upcoming examples; see Agarwal et al (2024).

Material parameter identification from full-field data
This section shows the potential of using the developed nested-POD-based ROM to identify the model's material parameters from full-field data.In the lack of experimental data, we take the solution of the nominal FOM (λ = 1558, A = 4000) at T = {0.15,0.25} as synthetic data.It is worth recalling that the ROM did not encounter these values in the training dataset.
The optimization problem (42) to identify the optimal value of the material parameters was implemented in Python and solved using the minimize function from the Scipy library (Virtanen et al, 2020).The hyperparameters of the L-BFGS-B algorithm are the maximum number of function evaluations and iterations ('maxfun' and 'maxiter'), both set to 15, 000.The function tolerance ('ftol') is set to 2.22 × 10 −9 , and the gradient tolerance ('gtol') to 1 × 10 −5 .The function tolerance dictates the precision of the objective function and the convergence criterion based on the gradient's magnitude, respectively.The step size for numerical approximation ('eps') is 1 × 10 −8 , and the maximum number of line search steps ('maxls') per iteration is 20.This setup allowed for efficient parameter space exploration while adhering to the constraints imposed by the problem's physical context.
The optimization process starts with an initial guess for the material parameters, set to θ 0 = [λ 0 , A 0 ] = [1800, 3800], but we verified that the same result is obtained if different initial conditions are chosen.The bound constraints [λ, A] = [(1000, 2000), (2000,6000)] ensure that the solution remains within the parameter range the ROM was trained for.
Table 1 summarizes the results obtained from the parameter identification with ROM and FOM, as well as the related computational cost.From the results, it is possible to observe the nested-POD-based ROM's effectiveness in performing many-query tasks.The material parameters identified with the ROM have a relative error of 0.95 % and 1.95 % for µ and u, respectively.This discrepancy is attributed to the low number of temporal snapshots (T = {0.15,0.25}) used for the identification problem.If more snapshots were considered for the identification problem, the error between the identified and correct material parameters is expected to decrease, but especially with the FOM, the computational cost would also increase significantly.Considering both the ROM training and the actual parameter identification, only a third of the time needed with the FOM was consumed.It is important to remark that the material parameter values identified using the FOM also deviate from the nominal values used to create the synthetic data.
It is worth remarking that the optimization in equation ( 42) converges within 6 steps using the FOM, but 84 model calls were necessary to approximate the Jacobian.Thus, the computational time could be reduced significantly when using adjoints.However, our focus is on repeated calibration in laboratory contexts, where the ROM would still be beneficial.

Case study: co-axial hydrogel bioprinting
Co-axial bioprinting is one technique that can be used to automate the fabrication process (Kjar et al, 2021).It uses a dual-nozzle in a core/sheath configuration system to extrude a hydrogel precursor and a crosslinking agent simultaneously; see Figure 2b.This core/sheath configuration, which enables immediate gelation and structural stability, can be adapted for specific hydrogel composition designs, allowing for multi-material deposition and improved resolution.It is crucial in tissue engineering for creating complex 3D structures with controlled properties like stiffness and porosity (Kjar et al, 2021).
Model validation in the first stage is then employed to make predictions of the co-axial printing process and understand the effect of material parameter uncertainties on the printed hydrogel.In this case, the model might need to be evaluated many times to evaluate different printing conditions.More precisely, the material parameters λ and A could present a significant source of uncertainty.Consequently, in Subsection 5.2.3, we exploit the ROM's efficiency and propagate the uncertainty in the material parameters onto two quantities of interest (QoIs) in the co-axial process, namely, the chemical potential at the nozzle's tip and the induced stresses within the domain due to swelling.
On the one hand, the chemical potential is directly related to the species concentration and can inform about the degree of crosslinking (DoC) of the extruded material (Hajikhani et al, 2021) that influences the gel's properties and determines its possible application.On the other hand, induced stresses due to fluid absorption at the co-axial printer's nozzle are often overlooked in the literature on hydrogel bioprinting.The focus has primarily been on induced shear stresses from extrusion (Chirianni et al, 2024).However, this does not capture the full picture of co-axial printing, where additional stresses are induced in the nozzle due to fluid absorption and hydrogel crosslinking, resulting in swelling and shrinking stresses.The mean and maximum induced stresses could be monitored to evaluate if the printed hydrogel has the right conditions for cells viability after printing (Chirianni et al, 2024).

Simulation setup
We consider a vertical hydrogel bar where the height is four times the length as marked by the dashed white box in Figure 2b.We assume that the diffusion time of the fluid into the hydrogel is relatively small compared to the velocity at which the hydrogel is extruded out of the printer's nozzle.This allows us to neglect the physics of the hydrogel transport through the co-axial nozzle.We study the transient free-swelling of the hydrogel bar in a two-dimensional setting, assuming a homogeneous mix of gel and solvent at the nozzle's tip.The material parameters are the same as in the benchmark problem mentioned earlier.Therefore, if the characteristic length in the previous example is l = 0.0015 m, the dimensions of the hydrogel bar are 0.0015 × 0.0060 m 2 or, equivalently, 1.5 × 6.0 mm 2 .This size is representative of a filament being crosslinked in the nozzle of a co-axial bioprinter.
We assume that only 75% of the gel bar is exposed to the solvent on the left and right edges.The remaining 25% of the left and right edges and the top surface have zero fluid flux for the solvent concentration boundary conditions.Due to the symmetry of the problem, only half of the domain is considered.The axisymmetric axis is marked in red in Figure 2b.A Robin boundary condition imposes the chemical potential on the exposed portion on the left edge as defined in equation (15).Symmetric (Neumann) boundary conditions are defined at the center of the domain.In terms of mechanical boundary conditions, the axisymmetric edge is subject to a Dirichlet boundary condition u x = 0, while all other boundaries are considered to be traction-free.
The domain is discretized into triangular elements with quadratic interpolation for displacement and linear interpolation for the chemical potential, resulting in Taylor-Hood elements as introduced in Section 2. Equation ( 23) is solved as a variational monolithic problem using the FEniCS package as a PDE solver.

Numerical results for the reduced-order model
In this section, we provide numerical experiments to validate and demonstrate the parametric nested-POD-based ROM's accuracy and computational efficiency, as we did for the benchmark case.The procedure in Figure 1 is also followed in this section to assemble the nested-POD-based ROM.The material parameters are varied within the same range as in the benchmark problem to create the training and testing sets; see Figure 5.
Figures 12 and 13 show the singular values and fraction of energy retained for the nested-POD algorithm for the chemical potential and displacement, respectively.On the one hand, Figures 12a, b and 13a, b report the behavior of the singular values and retained energy for the temporal dimension.They show the decay of the singular values and the retained energy for three out of the 30 parameter samples, as we did for the benchmark problem.It is seen from the figures that the singular values and retained energy behavior are very similar for the temporal dimension.We conclude from the analysis that six temporal modes are enough to retain η = 99.9999% of the energy.On the other hand, Figures 12c, d  and 13c, d display the decay of the singular values and retained energy for the spatial dimension.This analysis concludes that we can retain the desired fraction of energy η contained in the FOM with eight spatial basis functions.The nested-POD-based ROM is assembled with six and eight basis functions for the temporal and spatial domains.Four locations at the hydrogel bar domain are selected to verify the ROM's accuracy as marked in Figure 14.The root mean square error (RMSE) is measured over time at each location.
The ROM for the hydrogel bar developed using the nested-POD method is compared with the nominal FOM model (A = 4000, λ = 1558 and α = 0.66).Figure 14  figures show that the maximum RMSE between both models remains lower than 10 −2 for the chemical potential and displacement.
Figure 15 displays the chemical potential (Figure 15a) and displacement (Figure 15c) obtained with eight basis functions for the temporal and spatial domains at three different time steps, namely, T = {0.01,0.12, 0.25}.The corresponding errors with respect to the FOM are illustrated for the chemical potential in Figure 15b and for the displacement in Figure 15d, respectively.
As seen from Figures 14 and 15, the nested-POD-based ROM can predict the FOM to a satisfactory extent.From Figure 15, it can be observed that the maximum discrepancy between the FOM and ROM is less than 0.40% for the chemical potential and 0.24% for the displacement field, respectively.
Notice that the displacement field is not symmetric.The asymmetry in the displacement field suggests uneven swelling of the hydrogel bar during the printing process.This is relevant because this uneven swelling could impact the mechanical integrity of the final printed structure, particularly in complex multi-layered scaffold designs, where stability is crucial.

Uncertainty propagation of the material parameter
In this subsection, we exploit the ROM's efficiency to understand the effect of material parameters on two QoIs.We consider λ and A as random variables and propagate their uncertainty onto i. the chemical potential at the nozzle's tip and ii. the induced stresses within the domain due to swelling.
In principle, the probability density function (PDF) of the two material parameters can be determined as discussed below: 1. Bayesian updating of the material parameters from measurement data; the reader can, for instance, refer to Noii et al (2022), Grashorn et al (2023), andAnton et al (2024).However, acquiring species concentration-related and displacement measurement data during printing is a technological challenge.2. Bayesian inference from benchmark: Alternatively, one could quantify the uncertainty for the benchmark problem and then propagate this uncertainty in the case study.However, although both problems involve the same hydrogel material, the distinct boundary conditions and domain dimensions result in different chemical potential profiles and displacement magnitudes, thereby influencing the material response differently, as evidenced when comparing Figures 14, 15  2020)).The reader should notice that assuming a Gaussian distribution implies that λ and A could take unphysical negative values, but the probability of getting them at the given mean and standard deviation is approximately 7 × 10 −24 , i.e., essentially zero.Figure 16 depicts the probability density function (PDF) of along with the normalized histogram of 1000 samples generated from the PDF for λ (Figure 16a) and A (Figure 16b).It is observed from the figure that, indeed, the two random variables follow a normal distribution, and the generated samples remain within the range the ROM was trained (λ = [1000, 2000] and A = [2000,6000]).We solve the nested-POD-based ROM iteratively and observe how the uncertainty in λ and A propagates to the QoIs. Figure 17 shows the effect of the uncertain material parameters on the temporal evolution of the selected QoIs for 1000 ROM realizations.It is worth remarking that one iteration of the ROM took an average of 15 seconds, which led to a total computation time of around 4 hours on a standard laptop.The same 1000 iterations adopting the FOM instead will take up to 24 hours on the same laptop (an average of 90 seconds for each FOM iteration), which shows the advantage of adopting the nested-POD-based ROM in this many-query task.
Figures 17a and b report the evolution of the normalized chemical potential at the nozzle's tip. Figure 17a refers to the corner at (X = 0, Y = 0) and Figure 17b to the nozzle's tip center at the point (X = 0.5, Y = 0).From Figure 17 can be observed that variations up to 10% in the material parameters can lead to considerable differences in the diffusion-deformation process evolution.The effect of varying λ and A becomes more pronounced at the center of the nozzle's tip.
The normalized mean induced stress in the XX-direction at each time instance is computed as the average of the stress values over the domain as with A 2D−block the simulated hydrogel bar area and dA = dxdy representing the differential area element.
The maximum normalized induced stress in the XX-direction at each time instance is computed as the maximum value of the stress over the domain as σ max XX (t) = max  and σ max Y Y .It is observed from the figures that the level of induced stress has pronounced peaks at the beginning of the simulation, which can lead to some damage to the cells contained in bio-inks.
These induced stresses are normalized by the shear modulus G, which Chen et al (2020) found to be equal to 33 kPa from the experiments.Hence, the maximum induced stress can be amplified 80 times and reaches a maximum of 2600 kPa based on the results reported in Figure 17d.Just as a reference example, Chirianni et al (2024) studied the level of damage caused by the shear stresses during hydrogels extrusion in bioprinting applications and reported that cells could suffer damage up to 50% for induced stresses between 100 and 500 kPa.This threshold is clearly exceeded at the beginning of the performed simulation and can have significant consequences on the cell viability at the end of the process.One alternative to preventing such high-stress fluctuations could be implementing a feedback control mechanism that can monitor the mean and maximum induced stresses informed by the ROM and dynamically adjust the solvent concentration at the mixing zone (see Figure 2b) to keep the induced stresses under a desired limit.

Conclusions
In this paper, we demonstrated that POD-based algorithms are feasible for developing a ROM of a coupled diffusion-deformation model.The accuracy of the ROM was tested in two different scenarios: i. free-swelling of a 2D gel block, a well-established benchmark in the literature, and ii. a case study resembling a bioprinting process using a co-axial bioprinter.In both cases, we highlighted the computational benefits of adopting the ROM, where the computational cost of the ROM online evaluation was only 0.01% of running the full-order model.Furthermore, when using a finer mesh, we found that the ROM is even 5172 times faster than the FOM in simulating the time-dependent diffusion-deformation evolution in the co-axial printing application.
We demonstrated the benefits of adopting a ROM in many-query tasks, that is, parameter identification from full-field data and uncertainty propagation.In parameter identification, we showed that the computational resources spent in the offline phase for training data generation and ROM construction are repaid by a relatively large number of evaluations (90 on average) of the trained model in an online phase.Additionally, we exploited the simulation speed-up provided by the ROM and propagated the uncertainty in the material parameters to two quantities of interest in co-axial bioprinting: the chemical potential at the nozzle's tip and the stresses due to the fluid's absorption that lead to the hydrogel's swelling.The significant simulation speed-up provided by the ROM paves the way for more advanced applications, such as online feedback error compensation during the bioprinting process.This would ensure that the printed products meet design and functional requirements, ultimately aiming to produce consistent, high-quality hydrogel products to advance the reliability of bioprinting technologies.
Consider a continuum homogeneous elastomeric body B within the Euclidean space R d and its boundary ∂B = ∂Bu ∪ ∂B t = ∂B R .Here, ∂Bu denotes the displacement (Dirichlet) boundary, ∂B t the traction (Neumann) boundary, and ∂B R the fluid flux (Robin) boundary.The outward normal vector to the domain boundaries is denoted by n ∈ R d .The displacement field is defined as u : B → R d and the deformation is described by the displacement gradient tensor, H : B → R d×d , H = grad(u).

Fig. 1
Fig. 1 Illustration of the offline phase workflow for the POD-based ROM's construction.

Fig. 2
Fig. 2 Illustration of the simulated problems: a. benchmark: 2D square hydrogel block and b. case study: coaxial printing.Notes: i.The white dashed boxes mark the domains extracted for the simulation setup.ii.The shadowed area in b refers to the neglected part due to the axisymmetric nature of the case study.
Fig. 4 Benchmark: transient FOM simulation.a. normalized chemical potential, b displacement in x-direction, and b. normalized stresses at different points in the 2D hydrogel block.

Fig. 5
Fig. 5 Benchmark: material parameter sampling for the POD and nested-POD-based ROM construction.

FOMFig. 7
Fig. 7 Benchmark: POD-based ROM simulation results with six basis functions.a. normalized chemical potential and square root error with respect to the FOM; b. displacement magnitude at different simulation times, and square root error with respect to the FOM.

Fig. 8
Fig. 8 Benchmark: nested-POD-based ROM.SVD of the chemical potential snapshots matrix: a. singular values and b. retained energy for the temporal domain; and c. singular values and d. retained energy for the spatial domain.

Figure 5b .
Figure 5b.While both algorithms perform similarly with a small number of training samples, the nested-POD algorithm becomes more efficient as the number of training samples increases.

Fig. 9
Fig. 9 Benchmark: nested-POD-based ROM.SVD of the displacement snapshots matrix: a. singular values and b. retained energy for the temporal domain; and c. singular values and d. retained energy for the spatial domain.
Fig. 10 Benchmark: nested-POD-based ROM simulation results for primary variables.a. normalized chemical potential field and b. absolute ROM error with respect to the FOM, as well as c. displacement magnitude field and d. absolute ROM error with respect to the FOM at time steps T = {0.01,0.15, 0.25}.
Fig. 12 Case study: hydrogel bar -nested-POD-based ROM.SVD of the chemical potential snapshots matrix: a. singular values and b. retained energy for the temporal domain; c. singular values and d. retained energy for the spatial domain.
Fig. 13 Case study: hydrogel bar -nested-POD-based ROM.SVD of the displacement snapshots matrix: a. singular values and b. retained energy for the temporal domain; c. singular values and d. retained energy for the spatial domain.
Fig. 14 Case study: hydrogel bar.Nested-POD-based ROM simulation results with eight basis functions.a. normalized chemical potential and square root error with respect to the FOM; b. displacement magnitude and square root error with respect to the FOM.
Fig. 15 Case study: hydrogel bar.Nested-POD-based ROM simulation results.a. normalized chemical potential and b. absolute error of the ROM with respect to the FOM; c. displacement magnitude and d. absolute error of the ROM with respect to the FOM at time steps T = {0.01,0.15, 0.25} Fig. 16 Case study: hydrogel bar.Uncertainty propagation analysis.Probability distribution function and normalized histogram of the generated samples a. of λ and of A.

Figures
Figures 17c and d show the normalized mean σ meanXX and maximum σ max XX induced stresses in the XX-direction due to the diffusion of the solvent into the hydrogel, but identical results are obtained for σ mean Y Y . The DoC is a QoI Benchmark: nested-POD-based ROM v.s. the FOM.Computational analysis of the model parameters identification.