A Bayesian estimation method for variational phase-field fracture problems

In this work, we propose a parameter estimation framework for fracture propagation problems. The fracture problem is described by a phase-field method. Parameter estimation is realized with a Bayesian approach. Here, the focus is on uncertainties arising in the solid material parameters and the critical energy release rate. A reference value (obtained on a sufficiently refined mesh) as the replacement of measurement data will be chosen, and their posterior distribution is obtained. Due to time- and mesh dependencies of the problem, the computational costs can be high. Using Bayesian inversion, we solve the problem on a relatively coarse mesh and fit the parameters. In several numerical examples our proposed framework is substantiated and the obtained load-displacement curves, that are usually the target functions, are matched with the reference values.

phase-field formulation for quasi-brittle fracture is used. The variational phase-field formulation is a thermodynamically consistent framework to compute the fracture failure process. This formulation emanates from the regularized version of the sharp crack surface function, which was first modeled in a variational framework in [1]. Regularized fracture phenomena are described with an additional auxiliary smooth indicator function [2], which is denoted as crack phase-field (here indicated by d). Along with a mechanical field (denoted by u), a minimization problem for the multi-field problem (u, d) can be formulated. The main feature of such a variational formulation is to approximate the discontinuities in u across the lower-dimensional crack topology with the phasefield function d.
The resulting, regularized formulation leads to a diffusive transition zone between two phases in the solid, which corresponds to the fractured phase (i.e., d = 0) and intact phase (i.e., d = 1), respectively. The transition zone is determined by the phase-field regularization parameter , also well-known as the length-scale parameter. The parameter is related to the element size h and specifically h ≤ (e.g., = 2h). Therefore, sufficiently small length-scales are computationally demanding. To date, the focus in such cases was on local mesh adaptivity and parallel computing in order to reduce the computational cost significantly; see for instance [3][4][5][6][7][8][9][10][11]. Another recent approach is a global-local technique in which parts of the domain are solved with a simplified approach [12,13] that also aims to reduce the computational cost.
Generally, material parameters fluctuate randomly in space. In fact, the mechanical material parameters are spatially variable and, therefore, the uncertainty related to spatially varying properties can be represented by random fields. For instance, the material stiffness property has spatial variability. In fact, there are several sources of uncertainty including the class of extensometer or strain gauge resolution, uncertainty in the dimensional measurements, the classification and resolution of the load cell, misalignment of the specimen or strain measurement device, temperature effects, operator-dependent factors, data fitting routines and analysis methods, etc [14]. Therefore, in order to provide a reliable model, the uncertainty effect must be taken into account.
The main goal in this work is to identify such uncertain parameters for phase-field fracture problems. The underlying framework of parameter estimation using Bayesian inference is described in the following. Bayesian inference is a probabilistic method used to estimate the unknown parameters according to the prior knowledge. The observations (experimental or synthetic measurements) can be used to update the prior data and provide the posterior estimation. The distribution provides useful information about the possible range of parameters and their variations and mean. Markov chain Monte Carlo (MCMC) [15] is a common computational approach for extracting information of the inverse problem (posterior distribution). Metropolis-Hastings (MH) algorithm [16] is the most popular MCMC method to generate a Markov chain employing a proposal distribution for new steps. In practice, a reliable estimation of influential parameters is not possible or needs significant efforts. In [17,18], the authors used the Metropolis-Hastings algorithm to estimate the unknown parameters in field-effect sensors. It enables authors to estimate probe-target density of the target molecules which can not be experimentally estimated. We refer interested readers to [19,20] for more applications of Bayesian estimation. In the same line, other optimization approaches can be used to determine intrinsic material properties of the specimen from experimental load-displacement curves, see e.g., [21].
As previously mentioned, we consider fractures in elastic solids in this work. The principal material parameters are the shear modulus μ and the effective bulk modulus, K = λ + 2μ 3 (here λ denotes Lamé's first parameter) and Griffith's critical energy release rate G c . Using Bayesian inversion, the objective is to determine the unknown elasticity parameters.
For a homogeneous material, the stability requires positivedefiniteness of the elasticity tensor. For an externally unconstrained homogeneous solid, the conditions of structural stability needs that the fourth-order stiffness tensor is positive-definite. The condition for an isotropic, linear elastic medium gives rise to the shear modulus μ and the effective bulk modulus, K be strictly positive [22]. Regarding λ, the bound λ > − 2μ 3 may relate it to the shear modulus. Also, for the isotropic materials (as used in this paper) Poisson's ratio ν satisfies the condition −1 < ν < 1 2 [23]. These two elasticity parameters (λ and ν) are not well-suited for the estimation due to their bounds and dependency. Therefore, for the elasticity parameterization, we chose the eigenvalues, i.e., K and μ, and strive to estimate the joint probability, being updated jointly using MCMC.
Griffith's theory describes that crack propagation occurs if a certain reduction of the potential energy due to the change of surface energy associated with incremental crack extension reaches to its critical value [24]. Here, Griffith's critical energy rate G c measures the amount of energy dissipated in a localized fracture state [25], thus has units of energyper-unit-area. In case G c is unknown, one possibility is to employ the Bayesian setting for its identification. Physically speaking, there is a direct relation between G c and material stiffness, which means that in stiffer materials more energy is needed for the crack initiation. Computationally speaking, this value is independent of the elasticity parameters. Finally, since we should deal with three positive values (μ, K ) and G c , in order to remove the positivity constraints, we transfer these parameters and estimate the transfered values μ * = log(μ), K * = log(K ), and G * c = log(G c ). In our Bayesian framework, a reference value (obtained on a sufficiently refined mesh which termed here to the virtual observation) as the replacement of measurement will be chosen. Then, the posterior density of the elasticity parameters (joint probability) and the critical energy release rate is obtained. The computational costs can be high, specifically when an appropriate estimation is required inside multiphysics frameworks, see e.g. [3,[26][27][28]. Using Bayesian inversion, we strive to solve such problems with a coarser mesh and fit the parameters. The obtained load-displacement curve (as an important characteristic output) is matched with the reference value. In spite of using coarser meshes and therefore significantly lower computational costs (in terms of CPU timings), the accuracy of the solution is reliable (crack initiation and material fracture time estimated precisely).
The paper is organized as follows: In Sect. 2, we describe the variational isotropic phase-field formulation for the brittle fracture that is a thermodynamically consistent framework to compute the fracture failure process. In Sect. 3, the Bayesian inference is explained. We describe how the MH algorithm will be used to estimate the unknown parameters in phasefield fracture. Also we point out the critical points in the loaddisplacement curve, which must be estimated precisely with the Bayesian approach. In Sect. 4, the Bayesian framework is adopted to estimate unknown parameters in the phasefield fracture approach. In Sect. 5, three specific numerical Fig. 1 a Geometric setup: the intact region indicated by R and C is the crack phase-field surface. The entire domain is denoted by . The crack phase-field is approximated in the domain F . The fracture boundary is ∂ F and the outer boundary of the domain is ∂ . F is represented by means of d such that the transition area is 0 < d < 1 with thickness 2 . b Regularized crack phase-field profile for a different length scale. A smaller value for the length scale lets the crack phase-field profile converge to a delta distribution examples with different parameters and geometry will be given. We will use two proposal distributions (uniform and normal distribution) to sample the candidates and estimate the unknown parameters with different mesh sizes. Finally, in 6 we will draw paper conclusions and explain our future planes for employing Bayesian inversion in heterogeneous materials.

The primary fields for the variational phase-field formulation
We consider a smooth, open and bounded domain ⊂ R δ (here δ = 2). In this computational domain, a lower dimensional fracture can be indicated by C ⊂ R δ−1 . In the following, Dirichlet boundaries conditions indicated as ∂ D := ∂ , and Neumann boundaries conditions are given on ∂ N := N ∪ ∂C, where N is the outer boundary of and ∂C is the crack boundary. The geometric setup including notations is illustrated in Fig. 1a. The surface fracture C is estimated in F ⊂ ⊂ R δ . A region without any fracture (i.e., an intact region) is indicated by The variational phase-field formulation is a thermodynamically consistent framework to compute the fracture process. Due to the presence of the crack surface, we formulate the fracture problem as a two-field problem including the displacement field u(x) : → R δ and the crack phase-field d(x) : → [0, 1]. The crack phase-field function d(x) interpolates between d = 1, which indicates undamaged material, and d = 0, which indicates a fully broken material phase.
For stating the variational formulations, the spaces are used. Herein, W in denotes a closed, non-empty and convex set which is a subset of the linear function space W = H 1 ( ) (see e.g., [29]).

Variational formulation for the isotropic mechanical contribution
In the following, a variational setting for quasi-brittle fracture in bulk materials with small deformations is formulated. To formulate the bulk free energy stored in the material, we define the first and second invariants as with the second-order infinitesimal small strain tensor defined as The isotropic scalar valued free-energy function reads where K = λ + 2 including the stored internal energy and the imposed external energy is where τ is the imposed traction traction vector on ∂ N C := N ∪ C and the body-force is neglected. Following [1], we define the total energetic functional which includes the stored bulk-energy functional and fracture dissipation as where G c is the so called the Griffith's critical elastic-energy release rate. Also, H δ−1 refers to the (δ − 1)-dimensional Hausdorff measure (see e.g. [2]). Following [2], H δ−1 is regularized (i.e. approximated) by the crack phase-field d(x) (see e.g. [2]). Doing so, a second-order variational phase-field formulation is employed; see Sect. 2.3. Additional to that, a second-order stress degradation state function (intactedfractured transition formulation) is used as a monotonically decreasing function which is lower semi-continuous order; see Sect. 2.5.

Crack phase-field formulation in a regularized setting
Let us represent a regularized (i.e., approximated) crack surface for the sharp-crack topology (which is a Kronecker delta function) thorough the exponential function d(x) = 1 − exp −|x|/l , which satisfies d(x) = 0 at x = 0 as a Dirichlet boundary condition and d(x) = 1 as x → ±∞. This is explicitly shown in Fig. 1b for different length scales. Here, x is a position variable in the Cartesian coordinate system, meaning u and d have a certain value at each position within the geometry. The first observation through the explicit formulation is that, the crack phase-field d constituting a smooth transition zone dependent on the regularization parameter . In engineering or physics, is often a so-called characteristic length-scale parameter. This may be justified since this zone weakens the material and is a physical transition zone from the unbroken material to a fully damaged state. In practice, choices such as = 2h or = 4h are often employed. Following [30,31], a regularized crack surface energy functional for the second term in Eq. 8 reads based on the crack surface density function γ (d, ∇d) per unit volume of the solid. This equation is the so-called AT-2 model because of the quadratic term in PDE.
We set sharp crack surfaces as Dirichlet boundary conditions in C ⊂ . Hence, the crack phase-field d(x, t) is obtained from the minimization of the regularized crack density function as Figure 2 gives the numerical solution that arises from the minimization Eq. 10 and demonstrates the effect of different regularized length scales on the numerical solution. Clearly, a smaller length scale leads to a narrower transition zone (see Fig. 2c). That is also in agreement with the crack phase-field profile shown in Fig. 1b.

Strain-energy decomposition for the bulk free-energy
Fracture mechanics is the process which results in the compression free state. As a result, a fracture process behaves differently in the positive phase and in negative phase, see e.g. [32]. In the following, an additive split for the strain energy density function to distinguish the positive and negative phases is used. Instead of dealing with a full linearized strain tensor ε(u), the additive decomposition of the strain tensor based on its eigenvalues is used [5,31].
Herein, x ± := x±|x| 2 refers to the a Macaulay brackets for x ∈ R ± . Furthermore, ε + and ε − refer to the positive and negative parts of the strain, respectively. The {ε i } are the principal strains (i.e., the eigenvalues of the ε(u)) and the {N i } are the principal strain directions (i.e., the eigenvectors of the ε(u)). To determine the positive and negative parts of total strain ε, a positive-negative fourth-order projection tensor is such that the fourth-order projection tensor P ± ε projects the total linearized strain ε onto its positive-negative counterparts, i.e., ε ± = P ± ε : ε. Hence an additive formulation of the strain-energy density function consisting of the positive and the negative parts reads tension term Here, the scalar valued principal invariants in the positive and negative modes are Here, the first positive/negative invariant I 1 (ε) is strictly related to the tension/compression mode, respectively, meaning that if tr(ε) > 0 requires that we are in tension mode otherwise we are in compression state. The second invariant I 2 (ε) is mainly due to the positive and negative eigenvalue of the strain tensor, where its positive value requires that we are either in shear or in tension mode otherwise it is in compression. Thus, we distinguished between tension/compression and also a isochoric mode of our constitutive model, and only the positive part of the energy is degraded.

Energy functional for the isotropic crack topology
Due to the physical response of the fracture process, it is assumed that the degradation of the bulk material due to the crack propagation depends only on the tensile and isochoric counterpart of the stored bulk energy density function. Thus, there is no degradation of the bulk material in negative mode, see [31]. Hence, the degradation function denoted as g(d + ) acts only on the positive part of bulk energy given in Eq. 12, i.e., This function results in degradation of the solid during the evolving crack phase-field parameter d. Due to the transition between the intact region and the fractured phase, the degradation function has flowing properties, i.e., Following [31], the small residual scalar 0 < κ 1 is introduced to prevent numerical instabilities. It is imposed on the degradation function, which now reads The stored bulk density function is denoted as w bulk . Together with the fracture density function w f rac , it gives the the total density function with Formulation 2.1 (Energy functional for isotropic crack topology) We assume that K and μ are given as well as initial conditions u 0 = u(x, 0) and d 0 = d(x, 0). For the loading increments n ∈ {1, 2, . . . , N }, find u := u n ∈ V and d := d n ∈ W in such that the functional Herein, to make sure that phase-field quantity d lies in the interval [0, 1], we define d + to map negative values of d to positive values. In Formulation 2.1, the stationary points of the energy functional are determined by the first-order necessary conditions, namely the Euler-Lagrange equations, which can be found by differentiation with respect to u and d.

Formulation 2.2 (Euler-Lagrange equations)
Let K > 0, μ > 0 be given as well as the initial conditions u 0 = u(x, 0) and d 0 = d(x, 0). For the loading increments n ∈ {1, 2, . . . , N }, find u := u n ∈ V and d := d n ∈ W in such that Herein, E u and E d are the first directional derivatives of the energy functional E given in Formulation 2.1 with respect to the two fields, i.e., u and d, respectively. Also, D is a crack driving state function which depends on a state array of strain-or stress like quantities and δu ∈ {H 1 ( ) 2 : δu = 0 on ∂ D } is the deformation test function and δd ∈ H 1 ( ) is the phase-field test function.
Furthermore, the second-order constitutive stress tensor with respect to Eq. 18 reads with

Crack driving forces for brittle failure
Following [33,34], we determine the crack driving state function to couple between two PDEs. Hence, crack driving state function acts as a right hand side for the phase-field equation.
To formulate the crack driving state function, we consider the crack irreversibility condition, which is the inequality constraintḋ ≤ 0 imposed on our variational formulation. The first variation of the total pseudo-energy density with respect to the crack phase-field given in (17) reads Herein, the functional derivative of γ l (d, ∇d) with respect to d is Maximization the inequality given in Eq. 22 with respect to the time history s ∈ [0, t n ] reads We multiply Eq. 24 by l G c . Then Eq. 24 can be restated as Here, H := H(ε, t) denotes a positive crack driving force that is used as a history field from initial time up to the current time. Note that the crack driving state function D is affected by the length-scale parameter and hence depends on the regularization parameter.
The multi-field problem given in Formulation (2.3) depending on u and d implies alternately fixing u and d, which is a so called alternate minimization scheme, and then solving the corresponding equations until convergence. The alternate minimization scheme applied to the Formulation (2.3) is summarized in Algorithm 1.

The influence of the Ä on the stress-strain curve
In this part, the influence of the κ on the stress-strain curve is taken into account. Following [26], the homogeneous solu- Fig. 3 The influence of the κ on the stress-strain curve; left plot represent κ = κ( ) and right plot presents κ as a numerical parameter which is sufficiently small Initialization of alternate minimization scheme (k = 1): Algorithm 1: Alternate minimization scheme for Formulation (2.3) at a fixed loading step n.
tion at the quasi-static stationery state of the phase-field partial differential equation in the loading case takes the following form which results from the free Laplacian operator (•) = 0 assumption in Eq. 24 without any source terms (zero lefthand sides). Here, the crack driving state function D is given in Eq. 25. Because, we are in the elastic limit, prior to the onset of fracture, then no split is considered. We now aim to relate a stress state σ with the isotopic phase-field formulation. To do so, a non-monotonous function in the one-dimensional setting for the degrading stresses takes the following form by To see the influence of the κ in Eq. 29, the concrete material is considered. Following, [35] for a concrete material which has a brittle response, a typical values for material parameters reads, E = 29 GPa, σ c = 4.5 MPa and G c = 70 N/m. (30) We set = 0.0105 m. Thus, we can do a plot for the stressstrain curve through Eq. 29 by considering the material set given above. Figure 3 shows the effect of stress state for different strain loading. The black curve represents the stressstrain curve while κ = 0. Evidently, it can be grasped through Fig. 3 with κ = 0 the σ c is exactly σ c = 4.5 MPa as it is required for the concrete material, see [35]. If we consider κ = 0 as a function of characteristic length-scale, see Fig. 3 left, we can observe a good agreement with κ = 0 up to the peak point while after some strain value it becomes different as κ changes. Unfortunately, we can not observe any converged response if we consider κ as a function of . In contrast, if we chose κ sufficiently small, see Fig. 3 right, as much as κ reduced, in here less than κ ≤ 10 −4 , we observed a very identical response with κ = 0, thus it behaves as numerical parameters rather than material parameters.

Stochastic model for Bayesian inversion
In this section, we explain how we use Bayesian inversion to identify parameters. Then, we introduce a computationally effective numerical technique to estimate the unknown parameters.
In the phase-field model, the uncertainties arise from the elasticity parameters including the shear modulus μ and the bulk modulus K as well as Griffith's critical elastic energy release rate (material stiffness parameter) G c , which are assumed to be random fields. Specifically, we represent the parameters uncertainty (spatial variability) by a spatiallyvarying log-normal random field.
The Karhunen-Loéve expansion (KLE) expansion method is used to reduce the dimensionality of the random field. The field representing the elasticity parameters and the energy release rate can be characterized by its expectation and covariance using the expansion. Considering the probability density function P, the covariance function is which leads to the KL-expansion Here the first term is the mean value, k n are the orthogonal eigenfunctions, ψ n are the corresponding eigenvalues of the eigenvalue problem [36] D Cov (x, y)k n ( y) d y = ψ n k n (x), and the {ξ n (ω)} are mutually uncorrelated random variables satisfying where E indicates the expectation of the random variables. The infinite series can be truncated to a finite series expansion (i.e., an N KL -term truncation) by [36] For the Gaussian random field, we employ an exponential covariance kernel as where ζ is the correlation length as well as σ is the standard deviation.
For a random field, we describe the parameters using a KL-expansion. Considering the Gaussian field ξ(x), a lognormal random field can be generated by the transformationξ (x) = exp(ξ(x)). For instance, for the parameter K , the truncated KL-expansion can be written as

Bayesian inference
We consider Formulation 2.3 as the forward model y = G( (x)), where G : L 2 ( ) → L 2 ( ). The forward model explains the response of the model to different influential parameters (here μ, K , and G c ). We can write the statistical model in the form [37] where M indicates a vector of observations (e.g., measurements). The error term ε arises from uncertainties such as measurement error due experimental situations. More precisely, it is due to the modeling and the measurements and is assumed to have a Gaussian distribution of the form N (0, H ) with known covariance matrix H . The error is independent and identically distributed and is independent from the realizations. Here, for sake of simplicity, we assume H = σ 2 I (for a positive constant σ 2 ). For a realization θ of the random field corresponding to a realization m of the observations M, the posterior distribution is given by where π 0 (θ ) is the prior density (prior knowledge) and W m is the space of parameters m (the denominator is a normalization constant) [38]. The likelihood function can be defined as [37] π(m|θ) := 1 As an essential characteristic of the phase-field model, the load-displacement curve (i.e., the global measurement) in addition to the crack pattern (i.e., the local measurement) are appropriate quantities to show the crack propagation as a function of time. Figure 4 indicates the load-displacement curve during the failure process. Three major points are the following.
1 First stable position. This point corresponds to the stationary limit such that we are completely in elastic region (d(x, 0) = 1 ∀x ∈ \C).
2 First peak point. Prior to this point crack nucleation has occurred and now we have crack initiation. Hence, this Fig. 4 The schematic of load-deflection response for the failure process including primary path (prior to the crack initiation, i.e., between point 1 and 2) and secondary path (during crack propagation, i.e., between point 2 and 3) peak point corresponds to the critical load quantity such that the new crack surface appears (i.e., there exist some elements which have some support with d = 0). 3 Failure point. At this point, failure of the structure has occurred and so increasing the load applied to the material will not change the crack surface anymore.
The interval between point 1 and point 2 in Fig. 4 typically refers to the primary path where we are almost in the elastic region. The secondary path (sometimes referred to as the softening damage path) starts with crack initiation occurring at point 2. The whole process recapitulates the load-deflection curve in the failure process.
The main aim of solving the inverse problem followed here is to determine the random field to satisfy (38). We strive to find a posterior distribution of suitable values of the parameters μ, K , and G c in order to match the simulated values (arising from (27)) with the observations. The distribution provides all useful statistical information about the parameter.

Remark 3.1
Note that the principal parameters h, κ, and are mathematically linked in Formulation 2.3. Here, we use = 2h and κ is sufficiently small which is compatible with Sect. 2.7. In Sect. 5, the values of κ in the computations will be specified. Further, a sufficiently small h is chosen to obtain the reference solution.
The crack pattern is a time-dependent process (more precisely in a quasi-static regime, the cracking process is load-dependent), i.e., after initiation it is propagated through time. In order to approximate the parameters precisely, we estimate the likelihood during all time steps. Therefore, the posterior distribution maximizes the likelihood function for all time steps, and therefore we have an exacter curve for all crack nucleation and propagation times.
MCMC is a suitable technique to calculate the posterior distribution. When the parameters are not strongly corre-lated, the MH algorithm [39] is an efficient computational technique among MCMC methods. We propose a new candidate (so-called θ , i.e., a value of (μ, K , G c )) according to a proposal distribution (for instance uniform or normal distributions) and calculate its acceptance/rejection probability. The ratio indicates how likely the new proposal is with regard to the current sample. In other words, by using the likelihood function (40), the ratio determines whether the proposed value is accepted or rejected with respect to the observation (here the solution of Formulation 2.3 with a very fine mesh). As mentioned, fast convergence means that the parameters are fully correlated. A summary of the MH algorithm is given below.
Initialization: set prior data θ 0 and number of samples N . for i = 1 : N do 1. Propose a new candidate based on the proposal distribution θ * ∼ K(θ * | θ i−1 ).

Compute the acceptance/rejection probability
) .

Bayesian inversion for phase-field fracture
In this key section, we combine the phase-field algorithm from Sect. 2 with the Bayesian framework presented in Sect. 3.
First, we define two sampling strategies as follows: • One-dimensional Bayesian inversion. We first use N samples (according to the proposal distribution) and extract the posterior distribution of the first set e.g., (μ * , K * ) where other parameter is according to the mean value. Then obtained information is used to estimate the posterior distribution of next unknown (i.e., G * c ). In order to employ the estimated values, the exponential of the estimated parameters is used in the AT-2 model (see Algorithm 2).
• Multi-dimensional Bayesian inversion. A three-dimensional candidate (μ * , K * , G * c ) is proposed and the algorithm computes its acceptance/rejection probability.
To make the procedures more clarified we explain the multi-dimensional approach in Algorithm 3. Clearly, for the one-dimensional setting; for each parameter (e.g., θ * = (μ * , K * )), it can be reproduced separately. We will study both techniques in the first example and the more efficient method will be used for other simulations.
Here, n max is the sufficiently large value that is set by the user. Also, tol load is a sufficiently small value to guarantee that the crack phase-field model reached to the material failure time. Note, in part (iv) for the while-loop step, the criteria F n < tol load in the secondary path (i.e., during crack propagation state) guarantees that reaction force under imposed Dirichlet boundary surface is almost zero. Hence, no more force exists to produce a fractured state. We now term this as a complete failure point. But, in some cases, e.g., shear test as reported in [31], by increasing the monotonic displacement load,F n is not reached to zero. For this type of problem, if n < n max holds, then phase-field step (i.e., while-loop step) in Algorithm 3 will terminate.
The physical aim of using Bayesian inversion in phasefield fracture is adjusting the effective parameters to fit the solution with the reference values (see Remark 3.1). With (future) experiments (experimental load-displacement until the failure point), these can be used as reliable reference values.  The units are in kN/mm 2

Numerical examples
In this section, we consider three numerical test problems to determine the unknown parameters using given Bayesian inference. Specifically, we propose: • Example 1 the single edge notch tension (SENT) test; • Example 2 double edge notch tension (DENT) test; • Example 3 tension test with two voids.
The observations can be computed by very fine meshes (here the reference values) as an appropriate replacement of the measurements (see Remark 3.1). Regarding the observational noise, σ 2 = 1 × 10 −3 is assumed. The main aim here is to estimate the effective parameters (μ, K , and G c ) in order to match the load-displacement curve with the reference value. To characterize the random fields, we can use the KL-expansion with N KL = 100 and the correlation length ζ = 2 as well.
In all examples, the phase-field parameters set by κ = 10 −8 , and regularized length scale = 2h (respecting the condition h < l). The stopping criterion for the iterative Newton method scheme, i.e. the relative residual norm that is 3. Calculate the likelihood function (40) forF (during all n-steps, untiln) with respect to θ * where m n indicates the reference value at the n-th loading step. 4. Compute the acceptance/rejection probability ν(θ * | θ i−1 ). 5. Use Algorithm 2 to determine θ i (i.e., θ * is accepted/rejected). end Algorithm 3: The multi-dimensional Bayesian inversion for phase-field fracture. In the examples, the random fields modeled as a log-normal random field. For the numerical simulations, all variables are discretized by first-order quadrilateral finite elements.

Example 1: The single edge notch tension (SENT) test
This example considers the single edge notch tension. The specimen is fixed at the bottom. We have traction-free conditions on both sides. A non-homogeneous Dirichlet condition is applied at the top. The domain includes a predefined single notch (as an initial crack state imposed on the domain) from the left edge to the body center, as shown in Fig. 5a. We set A = 0.5 mm hence = (0, 1) 2 mm 2 , hence the predefined notch is in the y = A plane and is restricted to 0 ≤ |C| ≤ A. This numerical example is computed by imposing a monotonic displacementū = 1 × 10 −4 at the top surface of the specimen in a vertical direction. The finite ele- Fig. 8 The load-displacement curve for the one-dimensional (black) and three-dimensional (red) posterior distributions in addition to the ones for the prior distribution (green) and the reference value (blue) for the SENT example (Example 1) with h = 1/160 Fig. 9 The autocorrelation function for one-and multidimensional Bayesian inference in the SENT example ment discretization corresponding to h = 1/80 is indicated in Fig. 5b.
For the shear modulus, we assume the variation range (60 kN/mm 2 , 100 kN/mm 2 ). Regarding the the bulk modulus K , the parameter varies between 140 kN/mm 2 and 200 kN/mm 2 . Finally, we consider the interval between 2.1 × 10 −3 kN/mm 2 and 3.3 × 10 −3 kN/mm 2 for G c . Furthermore, we assume that in this example, the variables are spatially constant random variables (they are not random fields).
We solved the PDE model (Formulation 2.3) with μ = 80 kN/mm 2 , K = 170 kN/mm 2 , and G c = 2.7 × 10 −3 kN/mm 2 [31] and the displacement during the time (as the reference solution) with h = 1/320 was obtained. The main goal is to obtain the suitable values of μ, K , and G c such that the simulations match the reference value.
For this example, we use a uniformly distributed prior distribution and the uniform proposal distribution where χ indicates the characteristic function of the interval [θ 1 , θ 2 ] (where θ denotes a set of parameters). First, we describe the effect of each parameter on the displacement. As the elasticity constants (i.e., μ and K ) become larger, the material response becomes stiffer; crack initiation takes longer to occur. Additionally, a larger crack release energy rate (as an indicator for the material resistance against the crack driving force) delays crack nucleation and hence crack dislocation. All these facts are illustrated in Fig. 6.
The joint probability density of the elasticity parameters and the marginal probability of the posterior are shown in Figure 7 including one-and three-dimensional Bayesian inversions. The mean values of the distributions are μ = 84.2 kN/mm 2 and μ = 85.1 kN/mm 2 are obtained for the shear modulus. Here an acceptance rate of 27% is obtained. Regarding the material stiffness parameter G c , the acceptance rates are near 29%. The values are summarized in Table 1.
To verify the parameters obtained by the Bayesian approach, we solved the forward model using the mean values of the posterior distributions. Figure 8 shows the load-displacement diagram according to prior and posterior distributions. As expected, during the nucleation and propagation process, using Bayesian inversion results in better agreement (compared to the prior). Furthermore, a better estimation is achieved by simultaneous multi-dimensional Bayesian inversion. From now onward, this approach will be used for Bayesian inference.

The convergence of MCMC
A customary method to assess the convergence of the MCMC is the calculation of its autocorrelation. The lag-τ autocorrelation function (ACF) R : N → [−1, 1] is defined as where θ n is the n-th element of the Markov chain andθ indicates the mean value. For the Markov chains, R(τ ) is positive and strictly decreasing. Also, a rapid decay in the ACF indicates the samples are not fully correlated and mixing well. Figure 9 shows the convergence of the MCMC where the elasticity parameters and the crirical elastic energy. Also, we estimated the convergence observed in the multi-dimensional approach. As expected, the multi-dimensional approach converges slower than the one-dimensional one. As noted above, the phase-field solution depends on h and . A detailed computational analysis was for instance performed in [5,6]. In general, for smaller h (and also smaller ) the crack path is better resolved, but leads to a much higher computational cost. Figure 10 illustrates the crack pattern using different mesh sizes varying between h = 1/20 and h = 1/320. For these mesh sizes, we show the load-displacement diagram in Fig.  11 and the corresponding CPU time in Table 3.
Here we strive to solve the problem using a coarse mesh and employ MCMC to find parameters that make the solution more precise compared with the reference value. Figure 12 shows the obtained displacement with both prior and posterior distribution for h = 1/20. The efficiency of the Bayesian estimation is pointed out here since the peak point and the failure point are estimated precisely. The posterior distribu-tions are shown in the right panel as well. The estimation can also be performed for finer meshes: Figs. 13 and 14 illustrate the load-displacement curves for h = 1/40 and h = 1/80, respectively. In both cases, in addition to the precise estimation of the crack-initiation point and the material-failure point, the curve is closer to the reference value. Again, the posterior distributions are shown on the right panels. Finally, the mean values of the posterior distributions in addition to their acceptance rates are indicated in Table 2.

Example 2: Double edge notch tension (DENT) test
This numerical example is a fracture process that occurs through the coalescence and merging of two cracks in the domain. We consider the tension test with a double notch located on the left and right edge. The specimen is fixed on the bottom. We have traction-free conditions on both sides. A non-homogeneous Dirichlet condition is applied to the topedge. The domain has a predefined two-notch located in the left and right edge in the body as shown in Fig. 15a. We set A := 20 mm and B := 10 mm hence = (20, 10) 2 mm 2 .
For the double-edge-notches, let H 1 := 5.5 mm and H 2 := 3.5 mm with the predefined crack length of l 0 := 5 mm (Fig.  15a). This numerical example is computed by imposing a monotonic displacementū = 1 × 10 −4 at the top surface of the specimen in a vertical direction. The finite element discretization that uses h = 1/80 is indicated in Fig. 15b. According to the truncated KL-expansion, for the bulk modulus K , Eq. 37 gives the mean value ofK = 23.58 and the standard deviation of σ K = 0.28. Therefore, the parameter varies between 10 kN/mm 2 and 14 kN/mm 2 . For the shear modulus, the expectation ofμ = 22.8 and the standard deviation of σ μ = 0.23 leads to the variation Similarly, by using a KLexpansion for G c , we obtained the variation range between 8 × 10 −5 kN/mm 2 and 12 × 10 −5 kN/mm 2 . Figure 16 illustrates the effect of their different values including shear and bulk modulus and G c on the curve.
We assumed the uniform proposal distribution, namely the normal distribution As we aforementioned, the random field can be represented using the KL-expansion. In this example (32) is employed to parameterize the elasticity and energy rate parameters. The random perturbations are imposed on the ξ coefficients in the KL expansion. According to the proposal, the mean of the KL-expansion is updated.
Here we plan to study the effect of the number of samples N on the posterior distribution. Figure 17 shows the joint distributions of (K , μ), and the marginal distribution of G c using N = 3 000, N = 12 000, and N = 50 000. The calculations are done with h = 1/80, and h = 1/320 is used as the reference. As shown, with a larger number of samples, the distribution is close to a normal distribution. Table 4 points out the mean values in addition to the acceptance rate of all influential parameters.
As the next step, we use different mesh sizes for the Bayesian inversion using 15 000 samples. Figure 18 shows the crack pattern using different mesh sizes changing from h = 1/10 to h = 1/320. Finer meshes lead to a smoother  and more reliable pattern. Figure 19 depicts the loaddisplacement diagram using the prior values. With coarse meshes, the curve is significantly different from the reference including crack initiation. Using Bayesian inversion (see Fig. 17) enables us to predict the crack propagation and initiation more precisely. As the figure shows, even for the coarsest mesh (compare h = 1/10 to h = 1/320) the peak and fracture points are estimated precisely. For finer meshes (e.g., h = 1/80) the diagram is adjusted tangibly compared to the reference value. Finally, a summary of the mean values (of posterior distributions) and their respective acceptance rate is given in Table 5. The significant advantage of the developed Bayesian inversion is a significant computational cost reduction. As shown, for SENT and DENT, by using Bayesian inference for coarser meshes, the estimated load-displacement curve is very close to the reference values. We should note that the needed CPU time for h = 1/320 is approximately 4.5 hours; however the solution with h = 1/80 is obtained in less than 10 minutes. This fact pronounces the computational efficiency provided by Bayesian inversion, i.e., obtaining a relatively precise solution in spite of using much coarser meshes.

Example 3: Tension test with two voids
Here we consider the tension test where two voids are located in the domain as a more complicated example. The voids are used to weaken the material and to lead to crack nucleation/initiation without an initial singularity (i.e., a preexisting crack). The specimen is fixed on the bottom. We have traction-free conditions on both sides. A non-homogeneous  Dirichlet condition is applied to the top. Domain includes a predefined two voids in the body, as depicted in Fig. 20a. We set A = 0.5 mm hence = (0, 1) 2 mm 2 . The radius of left void is r 1 := 0.247 with the center c 1 := (0. 21, 0.197). The radius of the right void is r 2 := 0.0806 with the center c 2 := (0.7, 0.197). This numerical example is computed by imposing a monotonic displacementū := 1 × 10 −4 at the top surface of the specimen in vertical direction. The finiteelement discretization corresponding h = 1/40 is shown in Fig. 20b. Due to the resemblance to the first example (SENT) we use the same range of parameters and the variables are again spatially constant random variables. The load-displacement curves obtained from different values of μ, K , and G c are illustrated in Fig. 21.
This numerical example includes two voids results in multi-stage crack propagation. Hence, in the load-displacement curve, two peak points exist to demonstrate multi-stage crack propagation, see Fig. 21. Figure 22 shows the proposal distribution where a uniform prior distribution is used for Bayesian inversion with 10 000 samples. Here we use h = 1/160 as the reference solution and h = 1/80 is employed to estimate the parameters. In summary, the mean values are μ = 63 KN/mm 2 , K = 162 KN/mm 2 , and G c = 0.003 KN/mm 2 , and the acceptance rates are 28% (the elasticity parameters) and 21% (the critical energy rate).
We solve the forward model with the mean values of the estimated parameters. As Fig. 23 shows, the difference between the prior distribution and the reference solution is significantly large. By using Bayesian inversion, we could compensate this difference; crack initiation and material failure points are estimated precisely. Although multidimensional Bayesian inversion increases the computational costs (CPU time), the estimated solution is closer to the reference value.
Finally, we show the crack patterns obtained by different meshes varying from h = 1/10 to h = 1/160. We use Bayesian inference to estimate the unknown parameters with the three-dimensional approach for h = 1/20 and h = 1/40. As Fig. 25 illustrates, although the solution based on the posterior distribution is more precise (i.e., a better estimation of crack initiation and the fracture point) compared to the one based on the prior distribution, there is still a difference compared to the reference value. These results conform to Fig. 24, since the estimated crack pattern is considerably larger than the reality.

Conclusions and future works
In this work, we proposed a Bayesian approach to estimate material parameters for propagating fractures in elastic solids. For the fracture model, we adopted a phase-field approach. For the parameter estimation, we employed a Bayesian framework. We studied three phase-field fracture settings, and in each one, bulk and shear modulus as well as the critical elastic energy release rate were estimated with respect to a reference solution.
The developed Bayesian framework enabled us to provide useful knowledge about unknown parameters. By using Bayesian inversion, we could estimate the load-displacement curve precisely even with coarse meshes. For instance, in the first example (SENT), the diagram for h = 1/320 and h = 1/80 are essentially same, although a noticeable CPU time reduction is achieved. Interestingly, using even coarser meshes, the crack initiations and material fracture times can be estimated very well in all examples.
As one future application, the Bayesian approach will be used in multiscale problems to study crack propagation in heterogeneous materials, e.g., in composites. Due to their complexities, Bayesian inference will be employed to estimate material properties when the fiber-reinforced structures have a random distribution.