Experimental evaluation and uncertainty quantification for a fractional viscoelastic model of salt concrete

This study evaluates and analyzes creep testing results on salt concrete of type M2. The concrete is a candidate material for long-lasting structures for sealing underground radioactive waste repository sites. Predicting operational lifetime and security aspects for these structures requires specific constitutive equations to describe the material behavior. Thus, we analyze whether a fractional viscoelastic constitutive law is capable of representing the long-term creep and relaxation processes for M2 concrete. We conduct a creep test to identify the parameters of the fractional model. Moreover, we use the Bayesian inversion method to evaluate the identifiability of the model parameters and the suitability of the experimental setup to yield a reliable prediction of the concrete behavior. Particularly, this Bayesian analysis allows to incorporate expert knowledge as prior information, to account for limited experimental precision and finally to rigorously quantify the post-calibration uncertainty.

tainment is guaranteed and a contamination of the environment becomes extremely unlikely. The sealing of repositories in salt rock is performed by salt concrete and the concrete of type M2 (DBE 2004) is a candidate material for this objective in German disposal sites. The safety aspect and the extremely long operational lifetime of corresponding concrete structures require a profound knowledge of the material behavior including the mechanical properties. Particular phenomena of concrete, which have to be considered especially on long time scales, are creep and relaxation processes, which are generalized as viscoelastic material behavior (Christensen 2013;Lakes 1999;Creus 1986). The quantification of these phenomena is problematic and subject to large uncertainties.
Our study aims to model and predict the viscoelastic behavior of concrete M2. Therefore, accurate experiments, an appropriate model, experimental evaluation and model calibration with uncertainty analysis are required. Accordingly, we analyze creep tests with M2 specimens that run for about two years. To model the measured creep behavior, a fractional viscoelastic constitutive law (Bagley and Torvik 1985;Schmidt and Gaul 2002b;Mainardi 2010) was assessed, i.e. the material is modelled using fractional differential and integral operators (Samko et al. 1993;Oldham and Spanier 1974;Podlubny 1999;Diethelm 2010) in order to account for long-term creep over decades. Fractional calculus was used by several authors in the thirties and forties of the last century as an empirical method for describing viscoelastic material behavior, see the references in Mainardi (2010, Chap. 3.5). Particularly, Scott Blair (1947) used a fractional derivative of strain to model an intermediate material behavior interpolating between Hooke's and Newton's laws. Later, Koeller (1984) introduced the springpot 1 as a new rheological element that can be used in mechanical network models besides spring and dashpot, leading to a valid description of viscoelastic materials. A well-known model of this type is the fractional Zener model, see Mainardi (2010, Chap. 3.2), which is used in this article. The associated model parameters for concrete M2 are determined from the experimental data by a least-squares calibration and the identifiability of these parameters from the chosen experiments is evaluated by an informed Bayesian inversion (Congdon 2003;Besag et al. 1995).
The main question in this paper is whether a parameter identification for the considered model is possible from the given experimental data, particularly in view of the experimental duration which is short compared to the characteristic time scales of creep. More specifically, we examine the identifiability of the model parameters in terms of probability distributions. Thereby, we formulate a sufficient amount of constraints and reasonable soft information on the parameters in order to see how this influences or improves the parameter inference and whether it leads to plausible combinations of parameters with acceptable uncertainty intervals. Moreover, the differences between the results of least-squares calibration and Bayesian inversion are studied.
Well-known results on fractional calculus and the considered viscoelastic model are introduced in Sect. 2. The experimental setup, the data and their correspondence to the model parameters are illustrated in Sect. 3. After an introduction to Bayesian analysis and an explanation of the choice of prior parameter distributions in Sect. 4, the results of calibration and inference are discussed in Sect. 5 and the paper is concluded in Sect. 6.

Introduction to linear viscoelasticity
The mechanical behavior of a deformable body, aside from the general principles of mechanics (i.e. equilibrium and compatibility equations), is determined by a characteristic material behavior specified by constitutive equations, which relate forces and deformations. Assuming only small displacements of the body, forces and deformations can be described by the Cauchy stress tensor σ and the infinitesimal strain tensor ε as defined in any textbook on continuum mechanics (Gurtin 1981). The characteristic property of a viscoelastic constitutive equation is that the stress at a certain time instant depends on the entire history of the strain and vice versa. Besides this memory hypothesis, three other assumptions lead to a quite general viscoelastic constitutive equation. These assumptions are • causality (no future stress (resp. strain) state can affect the current strain (resp. stress) state), • linearity (the principle of superposition holds), • non-ageing (the material behavior is independent of shifts on the time scale).
The resulting viscoelastic constitutive equation is given by see Christensen (2013), Creus (1986), Lakes (1999). Here, the functions J ij kl are known as creep functions and describe the strain increase of a material under a unit step stress (given by the Heaviside function (t)). Particularly, it follows from (1) by integrating the Dirac delta function An alternative representation of the constitutive law in terms of stress relaxation can be obtained by reversing the roles of stress and strain in (1). This so-called relaxation representation, see e.g. Christensen (2013, Chap. 1.2), is not considered here, as we are interested in modelling creep. The general tensorial constitutive equations (1) can be simplified by using intrinsic symmetries of the stress and strain tensors as well as material symmetries. For the applications in this article, we consider the most simple case of isotropic material behavior. This leads to a description by only two remaining creep functions J h , J d , such that where are the hydrostatic (subscript h) and deviatoric (subscript d) components of the strain tensor ε and the stress tensor σ , see Gurtin and Sternberg (1962) for a detailed derivation. Hence, (3) and (4) are two one-dimensional independent functional equations that fully represent the constitutive law. Subsequently, we will omit the subscripts h and d of the creep functions as long as their general properties are discussed. From energy considerations and using the principle of fading memory (Christensen 2013, Chap. 3), one can obtain monotonicity conditions for creep functions and their first and second derivatives, i.e.
Moreover, as we are interested in describing viscoelastic solids, we assume an initial elastic material response and finite asymptotic creep, such that The most simple creep function that satisfies (6) and (7) is given by Christensen (2013, Chap. 1.5) where τ ε is a characteristic time quantity called retardation time. It determines the time scale in which creep takes place, see Fig. 1. As a mechanistic idea, the creep function (8) is associated with a mechanical network model consisting of two springs (modulus E 1 , E 2 ) and a dashpot (viscosity η), named Zener model (Zener 1948), see Fig. 2. Connecting Hooke's law for a spring σ = Eε and Newton's law for a dashpot σ = ηε, the Zener model is represented by an ordinary differential equation (Mainardi 2010) which leads, under the assumption of a unit stress step at t = 0 and the initial condition J (0) = 1 E 1 , to the creep function representation which is equivalent to (8) when defining the respective parameters by comparing the coefficients.
The applicability of the Zener model to fit real material data is limited as the creep process usually proceeds over a large time range. A possible improvement is given by the Prony series approach, considering several retardation times τ i , i = 1, . . . , n, leading to a creep function However, a creep function of the form (11) fits the experimental creep data only in the measured time ranges and predicts a fast decay of creep outside the measured time span, which results in a bad extrapolation (Bagley and Torvik 1985;Schmidt and Gaul 2002a). A different approach, related to fractional calculus, leads to a solution of this problem, see Sect. 2.3.

Introduction to fractional calculus
The following terms and properties are well-known results from fractional calculus, see e.g. Diethelm (2010), Oldham and Spanier (1974), Podlubny (1999), Samko et al. (1993). Let T > 0 and f ∈ L 1 (−∞, T ] be an integrable function. The Liouville-Weyl fractional integral of f of order α ∈ (0, 1) is introduced as (Mainardi 2010) where represents the Euler Gamma function It can easily be checked that (12) simplifies to a plain primitive of f at the limit α → 1.
To examine the limit for α → 0, consider the following. If, additionally, the function f is absolutely continuous (in the sense of Samko et al. (1993, Chap. 2, Def. 6.1)) and f and f vanish in an appropriate manner in the negative time limit, we obtain by partial integration of (12) The representation (14) shows, that (12) is a plain identity for α → 0. Together, this provides the intuition that the fractional integral for α ∈ (0, 1) is something between the identity and the first-order primitive. Furthermore, the Liouville-Weyl fractional derivative of Caputo type of order α ∈ (0, 1) of an absolutely continuous function f is given as (Mainardi 2010) which is (one possible definition for) a left inverse of the operator I α , i.e.
see Samko et al. (1993, Chap. 2, § 6.2). Similar to the fractional integral, the fractional derivative interpolates between the identity and the first-order derivative. To see this, consider the fractional derivatives of the function x(t) = t χ (0,∞) (t) in Fig. 3, where χ (0,∞) is the unit function vanishing on (−∞, 0]. Further details regarding fractional integrals and derivatives are given, e.g., in Oldham and Spanier (1974), Podlubny (1999), Diethelm (2010) and especially for functions with unbounded domains in Samko et al. (1993). Before we explain the role of fractional calculus in viscoelasticity, we want to emphasize one more theoretical aspect, which is important for the further reasoning, i.e. we want to investigate the problem which represents a linear functional equation with a given initial function. (17) is equivalent to the fractional initial value problem where is the fractional Caputo derivative with zero initialization time. The difference from (15) is merely the lower integration bound. The solution of (18) contains the so-called one-

Fractional Zener model
Comparing (1) and (14), one can see, as in Mainardi (2010, Chap. 3), that a slowly increasing creep function of the form leads to a (one-dimensional) constitutive law which is in view of (16) equivalent to The constitutive law (24) defines another rheological element besides spring and dashpot, called springpot (Koeller 1984), which is depicted as a rhombus, see Fig. 4a. However, a springpot is not sufficient to describe a viscoelastic solid, because (22) satisfies (6) but not (7). In view of the classical Zener model, one can consider a spring (modulus E 1 ) in series with a springpot (parameters p, α) with another spring (modulus E 2 ) in parallel, which is known as fractional Zener model (Mainardi 2010), see Fig. 4b. For the first spring one can obtain the creep function whereas the second spring parallel to the springpot leads to the problem  (17). Finally, using (21) leads to the overall creep function The creep behavior according to (27) is shown in Fig. 5. The limit case α → 1 results in the classical Zener model (10), which shows creep in a small time range, whereas smaller values of α ∈ (0, 1) lead to an extended time window of creep.
Comparing (8) with (27), one obtains a generalized retardation time which again determines the characteristic time scale in which creep occurs. Finally, to obtain a full model for the three-dimensional material behavior, one can consider a long-term creep function of the form (27) for hydrostatic (3) and deviatoric (4) strain and stress components (Schmidt and Gaul 2002a), i.e. one can choose in (3) and (4). The model in (29) and (30) . All these parameters have to be chosen or calibrated to obtain a predictive model, and classical literature values from simpler models cannot easily be transferred. Therefore, experimental data and a calibration procedure (possibly taking into account uncertainties and soft information) must be used.

Experimental design and parameter fitting
For calibrating the parameters in (29) and (30), we use data from an experiment described in Sect. 3.1, followed by a discussion on parameter identification in Sect. 3.2 and a reparametrization for Bayesian inversion in Sect. 3.3.

Test record and data
We ran a series of tests at the Materials Testing Institute (MPA) of the University of Stuttgart from January 2019 to January 2021. These experiments consist of creep and shrinkage tests for three specimens made of concrete M2. The cylinder-shaped concrete bodies were produced at the Materials Testing Institute in Brunswick and the tests at the MPA Stuttgart started when the specimens were 56 days old.
Just before starting the tests, we determined the compressive strength of M2 as 25 MPa on average. The creep tests are driven at a constant compressive axial load of 8.3 MPa, which is about 1 3 of the compressive strength and hence small enough that no cracks occur during the tests.
The shortening of the specimens in the axial direction is measured with the help of three dial gauges, while the lateral strain is determined using two strain gauges, see Fig. 6. The deformation is measured in the central part of the specimens, where the boundary effects have decayed and homogeneous stress and strain states can be assumed.
In parallel with the creep tests, we measure the autogenic shrinkage of the unloaded specimens with dial gauges, as shrinkage is not covered by the constitutive model and has to be subtracted from the measured data to obtain the elastic and creep strains. The deformation data are recorded in logarithmic time steps, which is important for parameter identification. All details regarding the test series are summarized in Table 1. The results of the measurements are depicted in Fig. 7.
For each of the three specimens, the strains measured by the multiple dial and strain gauges were averaged and shown together with error bars that represent an uncertainty of one empirical standard deviation. The bottom right graph shows the averages of the axial and lateral strain of the three specimens as well as the mean shrinkage.

Parametrization and parameter identification
In the given creep tests, we consider an approximately constant load that is applied (almost) instantaneously at the beginning of the test. This results in the stress history where σ 0 is the initial stress with hydrostatic part σ 0,h and deviatoric part σ 0,d , and (t) is still the Heaviside function. For this case, in view of (2) and using the constitutive laws (29) and (30), we obtain the strain components for t ≥ 0. For the specific compression test performed at the MPA, we consider stress and strain tensors that can be expressed in the principal coordinate system I of the test cylinders where σ zz , ε zz < 0 and ε xx = ε yy > 0. The assumptions in (33) are justified by the uniaxial loading in e I z -direction, the axisymmetry and the purely axial stretch in the principal directions of the specimens observed in the experiments, at least in the central area of the concrete bodies, where the measurements are performed. All quantities σ zz , ε xx = ε yy , ε zz that occur in (33) are measured within the experiments. A connection between (32) and (33) may be derived using (5) as Particularly, using (34), the constitutive law for the deviatoric components reduces to a scalar equation with Finally, using (32) and (34) together with the measured data, we are able to identify the parameters in (29) and (30).
To illustrate the influence of its parameters, we show the creep function (27) of the fractional Zener model in a log-log plot (Fig. 5, Fig. 9) and study the resulting S-shaped graphs. Clearly, the initial value of the creep function corresponding to the purely elastic component of the material response is given by E 1 as The slope of the central part of the graph, i.e. the velocity of the creep, is determined by the parameter α (Fig. 5). A change in the springpot coefficient p results in a shift of the graph along the time axis (Fig. 9a), whereas the parameter E 2 determines the asymptotic behavior as see Fig. 9b. The parameter identification is performed in two qualitatively different steps. First, a best-fit calibration is obtained from a nonlinear least-squares optimization. Second, the parameters are considered as random variables and their (joint) probability distributions are determined depending on prior assumptions and the experimental results using Bayesian inversion. The results of the first method can be found in Sect. 5.1. The second method is introduced in Sect. 4 and evaluated and compared with the first method in the rest of Sect. 5.

Reparametrization for Bayesian inversion
For the Bayesian inversion, we consider a slightly different, reorganized parameter set than in (29) and (30). The reorganization does not change the model equations, but makes the problem more easily accessible to expert knowledge to obtain more informed prior probability distributions. Particularly, the elastic material response can be described by Young's modulus E and the Poisson's ratio ν, which are related to the original elastic parameters E 1,h and E 1,d as (Gurtin 1973, Sect. 22) Moreover, instead of using the parameters p h , p d to identify the time range in which creep occurs, we consider characteristic retardation times τ h , τ d as in (28), viz.
such that (32) can be reformulated as

Introduction to Bayesian inversion
Suppose we have the following forward model: where M is the deterministic model (40) together with (34), θ denotes the vector of uncertain parameters E, ν, E 2,h , α h , τ h , E 2,d , α d , τ d for the model M, d represents the measurement data (ε xx (t i ), ε zz (t i )) i=1,...,n , and e represents unavoidable measurement errors. The purpose of Bayesian inversion is to infer the posterior distribution of the uncertain model parameters θ given prior knowledge about θ and the measurement data d.
The posterior distribution of the uncertain parameters θ can be obtained through Bayes' rule as (Congdon 2003) where p(θ ) is the prior distribution of the unknown parameters, which describes the prior knowledge (expert knowledge) about θ independent of the data d; p(d | θ ) is the likelihood function, which quantifies how probable are the measured data d for a given realization of uncertain parameters θ ; and p(d) is the so-called marginal likelihood. The marginal likelihood p(d) is a normalization factor for the posterior distribution. Since p(d) is not a function of the parameters θ , it is often ignored when making inference of parameters. This allows the posterior distribution to be written as proportional to the product of prior and likelihood, as reflected in the last part of the equation (42) (Cary and Chapman 1988).
When representing experimental errors and imprecisions via a multivariate Gaussian distribution for e in (41), the likelihood function can be expressed as where e is the covariance matrix of the error term e and n d is the number of measurements.
In this work, the number of measurements is n d = 2n, where n is the number of time steps in which we obtain strain data from the experiments. The variance 2 of the error term e is set as s 2 e,x = 0.1 mm m 2 for (ε xx (t i )) i=1,...,n and s 2 e,z = 0.5 mm m 2 for (ε zz (t i )) i=1,...,n according to the maximal empirical variances in the experiments (see Fig. 7). The covariance matrix e is a diagonal matrix with the diagonal elements s 2 e,x (first element to nth element) and s 2 e,z ((n + 1)th element to the last element), which assumes the independence of the measurement errors. A powerful tool for solving Bayesian inversion problems is the Markov chain Monte Carlo (MCMC) method (Besag et al. 1995;Worden and Hensman 2012). In this work, we will use the popular Metropolis-Hastings algorithm (Hastings 1970) to solve our problem.

Choice of prior distributions of model parameters
Choosing proper prior distributions for the model parameters is important for Bayesian inversion. In this work, we choose the prior distributions as follows. Associated samplings and plots are given in Sect. 5 in Figs. 11 and 13.

Young's modulus E: log-normal distribution
Young's modulus E is purely positive and can vary by orders of magnitude (i.e., it makes sense to look at ln E). We have expectations (mean value) of E from external data (see below), and these data show roughly symmetric scattering (standard deviation). To accommodate these pieces of information, we need a purely positive 2-parameter distribution on the log-scale that can be roughly symmetric. On the log-scale with a given mean and variance, the log-normal distribution is the max-entropic distribution (least information), and it is relatively symmetric after the back-transformation if the standard deviation is small compared to the mean.
From the literature (DBE 2004, Sect. 6.2), we get the mean and variance of E as μ E = 13.5 GPa and s 2 E = (1.5 GPa) 2 . 3 Then, we can calculate the corresponding parameters for the log-normal distribution as (Krishnamoorthy 2015) Poisson's ratio ν: beta distribution Poisson's ratio ν is bounded between 0 and 0.5, and we expect it to be mostly below 0.4, and typically around 0.3 (DBE 2004, Sect. 6.3). Therefore, we need a two-side-bounded distribution with two parameters (to tune the mean and the standard deviation). We see no specific reason for a symmetric distribution of ν. As the beta distribution (Krishnamoorthy 2015) is a bounded 2-parameter distribution with flexible shape, it is a good choice. Since the original beta distribution is defined on the interval [0, 1], we will use a scaled beta distribution for ν as 2ν ∼ Beta(α, β) with α = 8 and β = 4, which leads to a desired maximum probability density at about 0.3.

E 2,h and E 2,d : improper distribution
All we know about E 2,h and E 2,d is that they are positive. We have made the observation that over geological time scales and geological pressures rocks can act like fluids. And as design engineers, we fear that concrete may do the same. Hence, we choose the prior PDF p(θ) ∼ 1/θ for both E 2,h and E 2,d , because it puts a lot of the probability mass close to zero (which is the limiting case of asymptotically viscous behavior). This distribution is often suggested for variables where no other information than positivity exists (Robert et al. 2009), and it is the limiting case of the log-normal distribution for s E 2 → ∞. Although this is a so-called improper distribution (without normalization), it is enough for MCMC simulation.
α h and α d : uniform distribution In the given model, the admissible intervals for α h and α d are [0, 1]. We know nothing more than the interval over α h and α d . Because the uniform distribution is the max-entropic (least information) distribution under these conditions (Robert et al. 2009), we choose the prior distribution U(0, 1) for both α h and α d .
For the retardation time τ , we know that it is positive, and we believe that all orders of magnitude are equally probable. Hence, we start from a uniform prior on the real axis for log τ (which is similar to the assumptions for E 2,h and E 2,d ), but then softly exclude the cases that would be implausibly fast. This leads to a soft lower bound expressed by the CDF of the standard normal distribution. Thus, our PDF of the chosen prior distribution for both log τ h and log τ d can be expressed as where is the CDF of the standard normal distribution, τ 0 = 1000 d, and s log(τ ) = 1.
Besides the univariate choice of parameters described above, there are some further constraints relating several of the model parameters. These joint constraints result from the assumption that the lateral strain of the specimens is positive and monotonically increasing. An opposite behavior is only known from materials with a special micro-structure and cannot be supposed to occur for concrete without phase transition (which is implausible in our experiments). Necessary conditions for increasing lateral strain can be derived by considering the relations according to (34) using (29) and (30). For ε xx to be positive, the limits t → 0 and t → ∞ in (29) and (30), using (46), result in the conditions and respectively. The first condition (47) is, in view of (38), equivalent to ν > 0 and hence automatically fulfilled by the chosen univariate prior distribution for ν. The second condition (48) provides an additional constraint. To obtain conditions forε xx to be positive, we consider the approximation of the Mittag-Leffler function for t → 0 and t → ∞, viz. and see Erdélyi (1953, p. 210, (21)). Using these properties in (29) and (30) (omitting subscripts) leads to the expansions and where (52) together with (46) finally reveals the necessary condition forε xx > 0. The conditions (48) and (53) are used as additional constraints on Bayesian inversion, leading to joint prior distributions. A summary of the chosen prior distributions is given in Table 2.

Quality of fit
and the resulting predicted curves for hydrostatic and deviatoric components together with experimental data are shown in Fig. 10. The fitted curves provide a good representation of the experimental results and it is shown that the slow creep process of concrete M2 is (at least over the experimental duration) well described by the fractional Zener model. Small deviations, plausibly due to measurement errors just after loading the specimens, are visible. The fit of the deviatoric strain component is slightly better than for the hydrostatic part. For the visibly deviant data points we included ±2s e error bars (∼ 95% confidence level) into the figure to demonstrate the likely explanation by measurement inaccuracies. We omitted the other intervals to avoid cluttering the figures.

Identified parameter distributions
A Bayesian inversion was performed in view of the experimental data and assumed prior information of Sect. 4.2 using an MCMC with N = 5 · 10 7 samples. The resulting marginal histograms of the posterior distributions for the eight model parameters are shown in Fig. 11 together with the sampled prior distributions. Furthermore, the parameters (54) from the least-squares calibration (best fit) as well as the samples with the maximal likelihood value (ML) and with the maximal posterior probability (MAP) are visualized in the histograms. The histograms show comparably smooth unimodal distributions. The best fit and the ML sample are close, which indicates a good mixing in the MCMC. To understand this argument, note that the (weighted) least squares calibration and the maximum likelihood fit are, under the condition of a Gaussian likelihood, analytically identical. Marginal posterior distributions are in some cases widespread (Fig. 11g, Fig. 11h) and partly similar to the prior distributions (Fig. 11a, Fig. 11b). This indicates weakly informative data, i.e. the prior information influences the posterior distributions more than the knowledge from the experimental results. This assertion is supported by the location of the ML and MAP samples. The distance between ML and MAP is quite large and the MAP is typically located in regions of high prior probability. Possible reasons for the weak inference  (Fig. 7) and a remaining compensation among the parameters. The latter statement becomes visible in the scatter plots in Fig. 12, showing two-dimensional projections of parameter correlations. Particularly, the triples (E 2,h , α h , log(τ h )) and (E 2,d , α d , log(τ d )) show strong correlations. This indicates that a plain least squares optimization is inappropriate because it fails to quantify the remaining uncertainties that are substantial and nonlinear in this case.

Comparing constrained and unconstrained priors
The univariate prior distributions given in Sect. 4.2 are additionally constrained by the conditions (48) and (53), leading to joint priors for E 2,h , E 2,d , α h and α d . The constrained and unconstrained prior distributions as well as the associated posterior distributions are shown in Fig. 13. As (48) suggests, the univariate prior distribution of E 2,h is shifted to larger values, whereas lower values of E 2,d show a higher probability. Similarly, (53) renders larger values of α h and smaller values of α d more likely. Particularly, the uniform marginal priors of α h and α d become linear distributions due to the constraint (53). The posterior distributions of E 2,h and E 2,d substantially change under the influence of the additional constraints, whereas α h and α d are only slightly affected.

Comparison to calibration optimization
The graphs in Fig. 11 show a relatively wide gap between the parameter set (54) obtained from the least squares calibration and the MAP E = 13.48 GPa, ν = 0.31, E 2,h = 0.20 GPa, E 2,d = 0.16 GPa, α h = 0.16, α d = 0.12, τ h = 1.346 · 10 10 d, τ d = 4.461 · 10 11 d (55) resulting from the MCMC. As already mentioned in Sect. 5.2, this indicates weakly informative experimental data. This reconfirms the importance of a Bayesian analysis, already stated in Sect. 5.2. However, the argument this time is that the expert knowledge, literature values and plausibility constraints that can be used in Bayesian analysis significantly influence the result, which can once again be decision relevant.

Overall assessment
Another remarkable result of the MCMC is the gap between the global maximum of the multivariate posterior distribution (i.e. of MAP) and the maxima in the marginal histograms of E 2,h (Fig. 11c), log(τ h ) (Fig. 11g) and log(τ d ) (Fig. 11h). This indicates a highly nonlinear interrelation of the parameters in the joint posterior distribution. In order to check the plausibility of the location of the MAP, we consider bar plots in Fig. 14 that show the maximal posterior density within the classes of the posterior histogram. The bar plots maximize, where the MAP lies (as expected). Particularly, the plots for log(τ h ) (Fig. 14g) and log(τ d ) (Fig. 14h) are relatively flat, which shows that the associated parameter values with maximal probability are weakly identifiable.
For a plausibility check of the location of the ML, additional similar bar plots showing the maximal likelihood found within the histogram classes are shown in Fig. 14. Analogously, these bar plots should maximize near the ML, which is fulfilled in all graphs of Fig. 14. The bar plots for E (Fig. 14a) and ν (Fig. 14b) show a relatively large variation.
The plots in Fig. 14 emphasize that the nonlinear parameter interactions (and the discussed offsets of global to univariate maxima) are not merely an artefact of our chosen implementation, but facts of poor identifiability.
Finally, we demonstrate the consequences of parameter interactions and measurement uncertainties for the prediction of the future material behavior in Fig. 15. The graphs show simulations of the lateral and axial strains of our concrete specimens for 10 9 days under the same constant loading conditions as in the experiment, but for different parameter sets. Particularly, we plotted the simulated strains for 250 realizations of the prior and posterior distributions (satisfying ε xx (t) > 0 andε xx (t) > 0 for all t > 0). The results for all posterior realizations, least squares and MAP parameter sets remain quite close to the experimental data and to each other during the experimental time range. However, they become more widespread afterwards due to parameter compensation effects and uncertainties. Specifically, it can be seen that the least squares and MAP parameter sets do not satisfy the monotonicity condition for ε xx . This confirms that (53) is necessary but not sufficient.
Note that the forward uncertainty quantification shown here does not account for all possible sources of uncertainty. While we capture parametric uncertainties and measurement uncertainties within the fractional Zener model, our study does not include uncertainties pertaining to the choice of model equations (such as alternatives to the Zener model), possible long-term ageing effects of concrete, or uncertainties in boundary conditions. This  Long-term prediction of strains within the creep tests for different parameter sets. The prediction for the ML parameter set is omitted as the graphs of ML and the best fit realization are really close means that the posterior realizations of the parameters visualize a lower bound on the evolution of the uncertainties. Still, the identified spreading of the predicted strains on large time scales shows that the uncertainty quantification is really necessary. However, a comparison with the prior realizations (which fill the graphs almost completely) illustrates that the considered experiments yield a lot of information and that a longer experimental time could improve the prediction.

Conclusion
This article provides and analyzes experimental data from long-term creep tests on a salt concrete of type M2 and extensively discusses whether the material behavior can be modelled by a fractional viscoelastic constitutive law. It is shown that a fractional Zener model for hydrostatic and deviatoric stress and strain components yields a beautiful description of the experimental results over the full trial period of two years (Fig. 10). Furthermore, the applied Bayesian inversion yields a probability distribution with plausible combinations of model parameters (Fig. 11). However, the uncertainty quantification reveals that the parameters of the fractional model are only weakly identifiable from the current experimental setup, i.e. the (posterior) probability distributions of the parameters are comparably widespread and their shapes are mainly influenced by prior information (instead of experimental data). As a consequence, the parameter set obtained from the least squares calibration is quite different to the set with maximal posterior probability. Particularly, this effect occurs more severely for the parameters related to the hydrostatic stress and volumetric strain components and for both generalized retardation times. Possible reasons for the weak identifiability are given by large variances of the experimental results (which are not unusual in concrete creep tests, though), the strongly nonlinear interrelation of the model parameters within the model behavior and the experimental setup, which does not allow separate measurement of volumetric and deviatoric effects, which should be considered in prospective measurements.
The computational cost and scientific effort for a Bayesian analysis are admittedly much higher than for a least-squares calibration. However, we are convinced that the benefits of such a Bayesian analysis outweigh these costs, as it measures the parameter identifiability and provides a quantification of the post-calibration uncertainty, which represents the foundation for the computation of predictive uncertainties and the optimal design of experiments. Moreover, the Bayesian analysis reveals nonlinear parameter interactions and allows information from expert elicitation.
The results of this paper can be used in a next step to obtain a forward uncertainty quantification, i.e. an estimation of the transfer of parameter uncertainty to the uncertainty of deformation of concrete structures for long time spans. For an even more complete uncertainty quantification, uncertainties from the choice between competing models, ageing effects and uncertain boundary conditions can be included. Additionally, future experiments can be designed, that are optimal with respect to a certain criterion, e.g. a failure limit of a concrete structure. These tasks will be considered in our future work.
Funding Open Access funding enabled and organized by Projekt DEAL. This work is supported by the Federal Ministry of Education and Research of Germany within the project "ProVerB" under Grant Number 01IS17096B.
We thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for supporting this work by funding -EXC2075 -390740016 under Germany's Excellence Strategy.

Data Availability
The experimental data sets analyzed during this study are included in this published article and its supplementary information files (Online Resource 1).

Code Availability
The MATLAB code to perform the MCMC sampling is included in this published article and its supplementary information files (Online Resource 2).

Conflict of Interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.