A Machine Learning Surrogate Framework for Bayesian Calibration of Nonlinear Concrete Damage Models

1  Introduction

Concrete is a material of paramount importance within the engineering field [1]. Due to its microstructure containing inherent defects such as micropores and microcracks, concrete structures exhibit pronounced nonlinear and stochastic mechanical behavior under mechanical loading and environmental actions. Reliable characterization of such nonlinear responses is therefore essential for structural safety and rational design.

Systematic research through both experimental investigations and numerical modeling has been conducted to characterize the nonlinear structural behavior of damaged concrete structures [2,3]. Among these approaches, nonlinear damage models are widely incorporated into finite element frameworks to capture the complex mechanical responses. However, practical finite element analyses of real structures typically involve nonlinear constitutive relationships and iterative solution procedures, which result in substantial computational costs.

Damage model parameters govern the nonlinear material behavior and significantly influence the structural responses computed by finite element simulations. Because the mapping from model parameters to structural responses is implicit and highly nonlinear, reliable parameter estimation must rely on calibration against experimental observations. Bayesian calibration provides a rigorous framework for updating uncertain parameters using experimental data and quantifying predictive uncertainty [4,5]. Formulated as an inverse problem, Bayesian calibration requires repeated evaluations of the forward model to explore posterior distributions. For nonlinear finite element damage models, this results in significant computational burden.

To reduce the computational cost of repeated forward simulations, surrogate modeling techniques have been widely developed for approximating expensive nonlinear models. Representative data-driven approaches include Gaussian process regression (GPR), neural networks, polynomial chaos expansions (PCE), and support vector regression (SVR) [6–10]. These methods provide efficient computational tools for uncertainty analysis and, when integrated with Bayesian inference, can substantially accelerate posterior estimation. Related studies have further demonstrated the value of surrogate-based approaches in improving computational efficiency and supporting uncertainty aware approximation in computational mechanics [11,12]. Meanwhile, Bayesian calibration has been applied to parameter identification and uncertainty quantification in engineering systems [13,14]. Further studies have considered model-bias correction, reliability assessment, and durability-related calibration in concrete structures [15–17]. However, existing surrogate-based calibration frameworks primarily emphasize approximation accuracy and sampling efficiency, while the treatment of structural model bias is often simplified. In practice, calibration uncertainty arises from two conceptually distinct sources: observational noise and structural model bias. Observational noise originates from measurement inaccuracies and experimental variability and is typically modeled as a stochastic error term in the likelihood function. In contrast, structural model bias corresponds to the systematic deviation between model predictions and experimental observations that persists even when the model parameters are optimally calibrated. Such bias may result from simplified modeling assumptions, incomplete physical representation, or idealized constitutive descriptions in computational models [18]. This distinction is illustrated in Fig. 1.

Figure 1: The distinction between test error and model bias.

While observational noise is routinely incorporated into Bayesian calibration, structural model bias is often neglected or implicitly absorbed into parameter uncertainty [19,20]. Moreover, surrogate approximation errors are sometimes conflated with structural model bias, leading to ambiguity in uncertainty interpretation and potentially biased posterior estimates [21]. Surrogate-assisted Bayesian calibration can substantially reduce the computational burden associated with repeated model evaluations [22]. These considerations highlight the need for an integrated probabilistic framework that simultaneously addresses computational efficiency and systematic model bias within Bayesian inference. In the proposed framework, the surrogate model is introduced to approximate the finite element response for computational acceleration, whereas the response-level model bias term represents the systematic discrepancy between the finite element model and the experimental response. Therefore, the model bias considered in this study refers to the finite element model bias rather than the surrogate approximation error.

To address these limitations, this study introduces a machine learning surrogate framework for Bayesian calibration of nonlinear concrete damage models. The main methodological contribution lies in the explicit separation between data-driven surrogate approximation for accelerating Bayesian inference and response-level model bias representation for quantifying systematic finite element model inadequacy. A data-driven surrogate is constructed to approximate the nonlinear relationship between model parameters and structural responses obtained from nonlinear finite element simulations, and is incorporated into the Bayesian sampling loop to accelerate posterior inference. Furthermore, an explicit probabilistic model bias term is introduced at the response level to represent structural model inadequacy while distinguishing it from observational noise. The proposed framework is demonstrated through the calibration of a nonlinear finite element damage model for a reinforced concrete column.

2  Machine Learning Surrogate Framework

Bayesian calibration of nonlinear concrete damage models requires repeated evaluations of the forward operator, which is defined in the present study through a nonlinear finite element procedure involving constitutive integration and iterative equilibrium solution. When such a high-fidelity model is embedded into a sampling-based inference framework, the computational cost of posterior estimation becomes prohibitive. To overcome this difficulty, a machine learning surrogate framework is developed by introducing a data-driven approximation of the forward mapping and integrating it with Bayesian inverse analysis. The proposed framework combines nonlinear damage modeling, surrogate-based approximation, model bias representation in the likelihood formulation, and adaptive posterior sampling within a unified computational procedure. The overall calibration workflow is illustrated in Fig. 2.

Figure 2: Bayesian model calibration process.

2.1 Nonlinear Damage Forward Model

The nonlinear mechanical behavior of concrete is characterized by stiffness degradation and strength softening, which cannot be adequately captured by traditional elastic or plastic models [23]. In damage mechanics, internal variables are introduced to represent the progressive deterioration of micro defects such as cracks and pores, thereby enabling the simultaneous modeling of elastic degradation and plastic deformation [24–28].

In this study, a bi-scalar elastoplastic damage model for concrete [29] is adopted as the forward constitutive model. Two independent damage variables are incorporated to distinguish tensile softening and compressive weakening mechanisms. For the uniaxial case, the stress–strain relationship is expressed as:

σ=(1−dt)σ¯t+(1−dc)σ¯c=(1−dt/c)σ¯t/c(1)

where dt and dc denote the tensile and compressive damage variables, respectively; σ¯t and σ¯c are the tensile and compressive components of the effective stress; and σ represents uniaxial stress. The subscript t/c, as appropriate, denotes tensile and compressive quantities. The effective stress σ¯ is defined by:

σ¯=E(ε−εp)=Eεe(2)

where E is the elastic modulus of the undamaged material, while ε, εp and εe represent total strain, plastic strain, and elastic strain, respectively. A one-dimensional empirical damage evolution equation can be obtained by simultaneously considering the following equation:

dt/c={1−ρt/cnt/cnt/c−1+(xt/c)nt/cxt/c≤11−ρt/cαt/c(xt/c−1)2+xt/cxt/c>1(3)

where ρt/c=ft/cEεt/c,nt/c=Eεt/cEεt/c−ft/c,xt/c=εeεt/c. Here, ft, εt and αt denote the tensile strength, peak strain, and descending branch parameter, respectively, while fc, εc and αc denote the corresponding compressive quantities. These parameters are selected according to relevant design provisions [30].

To account for the coupling between plasticity and damage, the plastic strain is decomposed into positive and negative components:

εp=εtp+εcp(4)

where εtp represents the tensile plastic strain, and εcp represents the compressive plastic strain. Since tensile plastic strain values are typically negligible, yielding εp=εcp. The corresponding plastic evolution function is given as:

εp=∅pmax(εce)=ξp[exp⁡(npdc)−1]max(εce)(5)

where ∅p is the plastic strain coefficient, ξp and np are parameters governing the evolution of plastic deformation. The uniaxial constitutive law is illustrated in Fig. 3, comprising four stages: linear, nonlinear, descending, and convergence [23].

Figure 3: Uniaxial constitutive law of concrete.

The confinement effect of transverse reinforcement is incorporated through a confinement coefficient AC, which is calculated based on the equivalent lateral confinement stress Fl as follows [31]:

AC=2.2541+7.94Flfc−1.254−2Flfc(6)

This coefficient is used to account for the effect of transverse reinforcement by enhancing the effective compressive strength and adjusting the corresponding strain and damage-related parameters in compression.

Based on the above constitutive formulation, the nonlinear structural response is governed by a set of material parameters. In this study, eight constitutive parameters are treated as uncertain and subject to Bayesian calibration. The uncertain parameter vector is defined as:

θ={E,ft,εt,αt,fc,εc,αc,Fl}(7)

Given θ, the finite element solver produces the corresponding structural response:

y=M(θ)(8)

where M(⋅) denotes the nonlinear mapping defined by the finite element solver. Each evaluation of M(θ) requires a full nonlinear finite element analysis. The physical meanings of these parameters are summarized in Table 1.

2.2 Data-Driven Surrogate Modeling of the Forward Mapping

The forward operator M(θ) is evaluated through a full nonlinear finite element analysis involving iterative equilibrium solution and constitutive integration. When embedded into sampling-based Bayesian calibration, the large number of required operator evaluations leads to substantial computational cost. To alleviate this difficulty, a data-driven surrogate model is constructed to approximate the implicit operator:

M⌢(θ)=M(θ)(9)

The surrogate model replaces the computationally intensive finite element simulation during posterior sampling, while retaining the nonlinear relationship between model parameters and structural response. The construction procedure comprises training data generation, regression-based operator learning, hyperparameter optimization, and predictive error assessment.

Beyond its role as a computational accelerator, the surrogate model is treated in this study as an integral component of the inverse problem formulation. In particular, the surrogate directly influences the structure of the likelihood function and thus affects the geometry of the posterior distribution. Consequently, surrogate model selection must be aligned not only with approximation accuracy but also with the requirements of stable and reliable Bayesian inference.

2.2.1 Support Vector Machines Formulation

The surrogate is trained using a dataset generated from the FE model. The admissible parameter domain is defined as θ∈Ωθ⊂Rd, where d denotes the number of material parameters. A Latin Hypercube Sampling (LHS) strategy is employed to generate parameter samples within the predefined parameter ranges [32,33]. This sampling design provides sufficient coverage of the admissible parameter domain while reducing clustering effects in the high dimensional parameter space.

For each sampled parameter vector the FE solver is evaluated to obtain the corresponding structural response. The resulting dataset, consisting of parameter samples and corresponding structural responses, denoted by

D={(θi,yi)}i=1N(10)

is used as the training set for surrogate construction. The sample size is determined by balancing surrogate accuracy and computational cost. Its sufficiency is further examined through learning curves constructed from separate training and test sets.

To select an appropriate surrogate model, several candidate methods, including Kriging, PCE, and SVR, were preliminarily evaluated under identical training and testing settings. This comparison was based on the finite element simulation dataset generated for the reinforced concrete column case. Details of the loading protocol, the force-displacement response at the column top, and the train-test split are provided in Section 3.3. Their prediction accuracies were compared using normalized root mean square error (NRMSE), the mean absolute percentage error (MAPE) and the coefficient of determination (R2), as summarized in Table 2. The NRMSE was calculated by normalizing the RMSE with the range of the true response values. This preliminary comparison indicates that SVR provides the most accurate approximation for the present dataset.

SVR is employed to approximate the nonlinear mapping from model parameters to structural response [34]. Given the training dataset, SVR seeks a regression function that achieves a balance between prediction accuracy and model flatness. Based on the ε-insensitive loss function, deviations lying within a prescribed tolerance band are ignored, while larger deviations are penalized. In the kernel-based representation, the resulting regression function can be written as [35]:

f(θ)=wTφ(θ)+b=∑i=1N(αi−αi∗)φ(θi)Tφ(θ)+b(11)

where w is the weight coefficient vector, b represents the bias parameter, αi and αi∗ are the Lagrange multipliers. By using a kernel function to define the inner product, the expression becomes:

K(θi,θ)=φ(θi)φ(θ)(12)

and the regression function becomes:

f(θ)=∑i=1N(αi−αi∗)K(θi,θ)+b(13)

In this study, a Gaussian kernel with anisotropic scaling is adopted to account for the different sensitivities of the structural response to individual parameters. The associated hyperparameters are optimized using a span leave-one-out error estimator in combination with the Covariance Matrix Adaptation Evolution Strategy (CMA-ES).

2.2.2 Justification of Surrogate Model Selection

The choice of surrogate model plays a critical role in surrogate-assisted Bayesian calibration, as it directly affects predictive accuracy, numerical stability, and the efficiency of posterior sampling.

First, the nonlinear mapping from constitutive parameters to structural response is characterized by a high-dimensional input space and a limited number of training samples due to the computational cost of FE simulations. In such small sample regimes, SVR provides strong generalization capability through the principle of structural risk minimization, which enables effective control of model complexity and reduces the risk of overfitting.

Second, the ε-insensitive loss formulation in SVR enhances robustness to localized irregularities and numerical noise in the training data. This property is particularly beneficial for nonlinear FE responses, where small numerical inconsistencies may arise during nonlinear finite element simulations. By filtering out minor fluctuations, SVR yields a stable approximation of the underlying physical mapping.

Third, compared with GPR, which requires inversion of covariance matrices and may become computationally demanding or numerically sensitive for large datasets or high-dimensional outputs, SVR provides a more scalable alternative while maintaining sufficient approximation accuracy for the present application. Although GPR provides predictive uncertainty, SVR is adopted here as a deterministic surrogate for computational acceleration, while model bias is represented separately by the response-level bias term.

More importantly, from the perspective of Bayesian inverse problems, a smooth and numerically stable surrogate approximation is helpful for obtaining stable likelihood evaluations and efficient posterior sampling. Kernel-based SVR provides relatively smooth mappings in the parameter space, which is beneficial for stable posterior sampling in the present application.

Therefore, the selection of the surrogate model in the present study considers not only approximation accuracy, but also numerical stability and suitability for Bayesian calibration.

2.2.3 Surrogate Validation and Error Assessment

For Bayesian calibration, it is necessary to assess whether the approximation error introduced by the surrogate model significantly affects likelihood evaluation and posterior inference. To this end, the surrogate predictive accuracy is evaluated using separate training and test sets through statistical error metrics, including NRMSE, MAPE and R2. In addition, the maximum relative deviation between surrogate predictions and FE responses is examined to assess the worst-case approximation behavior.

The magnitude of surrogate error is then compared with the experimental measurement uncertainty adopted in the likelihood formulation. Since the surrogate approximation error is found to be substantially smaller than the observational noise level, its contribution to the total uncertainty is considered negligible in the present study. Accordingly, surrogate error is not explicitly included in the Bayesian calibration model, and the learned operator M⌢(θ) is used directly in place of the FE model during posterior inference.

2.3 Bayesian Model Calibration as an Inverse Problem

Bayesian calibration formulates parameter identification as an inverse problem in which uncertain model parameters are updated by combining prior knowledge with experimental observations. The observed structural response is denoted by d, and the corresponding model prediction is given by the forward operator M(θ), where θ is the vector of uncertain constitutive parameters. In the subsequent surrogate-assisted implementation, M(θ) is replaced by the trained surrogate operator M⌢(θ). In structural engineering applications, discrepancies between model predictions and experimental data may arise from both random observational noise and systematic model inadequacy. To obtain reliable parameter estimates, these two sources of uncertainty should be clearly distinguished and incorporated into the likelihood formulation.

2.3.1 Bayesian Model Calibration without Model Bias

The typical Bayesian forward model is formulated as:

d=M⌢(θ)+εe(14)

where εe denotes the observational noise vector. It is commonly modeled as a realization of a Gaussian random vector with zero mean and covariance matrix ∑e, i.e., εe∼N(0,∑e). In practice, when multiple experimental observations from the same test series are available, the covariance matrix can be estimated as:

∑e=1n−1∑i=1n(di−d¯)(di−d¯)T(15)

where di denotes the i-th experimental observation vector and d¯ is the corresponding sample mean vector.

The goal of Bayesian inversion is to infer the posterior distribution of θ given prior knowledge and measurement data. According to Bayes’ theorem, the posterior distribution of the uncertain parameters is described as follows:

p(θ|d)=p(d|θ)p(θ)p(d)(16)

Here, p(θ) represents the prior distribution, while p(d|θ) is the likelihood function, i.e., the probability of observing d given θ. p(d) is the marginal distribution function of d, which is usually omitted during posterior inference. Therefore, the posterior is proportional to the product of the prior and likelihood. The likelihood function can be derived as:

p(d|θ)=exp⁡(−12(d−M⌢(θ))T∑e−1(d−M⌢(θ)))(2π)nddet(∑e)(17)

This formulation implicitly assumes that the computational model is structurally adequate and that discrepancies between model predictions and experimental data are solely attributed to random observational noise.

2.3.2 Bayesian Model Calibration with Model Bias

In practice, systematic discrepancies may remain between model predictions and experimental observations even when the parameters are optimally calibrated. Such discrepancy, commonly referred to as model bias, reflects the effects of structural simplifications, numerical approximations, and constitutive idealizations that cannot be eliminated through parameter updating alone.

To account for this effect, the forward model in Eq. (14) is extended as:

d=M⌢(θ)+εM+εe(18)

where the enhanced term εM denotes model bias introduced at the response level.

In this study, the model bias is treated as an effective statistical discrepancy rather than a fully latent model form uncertainty jointly inferred with the model parameters. For a dataset consisting of N experimental samples, the residual vector between experiments and finite element predictions is defined as:

r(k)=d(k)−yFEM(k),k=1,…,N(19)

where yFEM(k) denotes the corresponding finite element prediction. In this study, the residual samples are constructed using finite element predictions from the simulation dataset generated before Bayesian calibration, rather than from posterior calibrated responses. The same set of experimental responses is used as the reference response for residual construction and parameter calibration. The estimated bias statistics are introduced as fixed hyperparameters in the likelihood function.

The collection of residuals is regarded as empirical realizations of the random vector εM, which provides an empirical probabilistic representation of response-level model bias. Similarly, εM is assumed to follow a normal distribution with a mean vector μM and covariance matrix ∑M, i.e., εM∼N(μM,ΣM). The mean vector is estimated from the sample average of the residuals, while the marginal variances are estimated from the residual samples. If output correlations are neglected, the covariance matrix can be taken as diagonal. In the present study, however, correlation between response quantities is represented using a squared exponential kernel:

(∑M)ij=σM(xi)σM(xj)exp⁡(−(xi−xj)22l2)(20)

where (∑M)ij denotes the covariance between the model bias components at response coordinates xi and xj, σM(xi) and σM(xj) are the marginal standard deviations of the model bias estimated from the residual samples, and l is the kernel length scale controlling the correlation range. The length scale is optimized using the residual samples.

Assuming independence between the model bias and observational noise, the total covariance is given by:

ΣM+e=ΣM+Σe(21)

Accordingly, the likelihood function is written as:

p(d|θ)=N(M(θ)+μM,ΣM+e)(22)

2.4 Computational Efficiency and Acceleration Analysis

The computational cost of Bayesian calibration is primarily governed by the number of evaluations of the forward model. In the present problem, each evaluation of the nonlinear forward operator requires a full finite element simulation involving incremental loading, constitutive updates, and iterative equilibrium solution. As a consequence, a single model evaluation may be computationally expensive.

Bayesian inference is performed by sampling from the posterior distribution of the uncertain parameters. Since the posterior distribution generally does not admit a closed form solution, sampling-based methods are required, and a large number of model evaluations is typically needed to ensure stable posterior estimation. When the forward model is directly represented by a nonlinear finite element solver, repeated evaluations can render the calibration procedure computationally prohibitive.

To overcome this difficulty, the surrogate model introduced in Section 2.1 is employed to approximate the nonlinear mapping between model parameters and structural responses. Once trained using a limited number of finite element simulations, the surrogate can provide rapid predictions for newly proposed parameter samples at negligible computational cost. In the proposed framework, the expensive finite element solver is therefore used only for training data generation, while posterior sampling is carried out using the surrogate model.

Posterior inference is performed using the Metropolis Hastings algorithm [36,37]. Candidate samples are generated from a proposal distribution and accepted or rejected according to the posterior density. To improve numerical stability, the likelihood function is evaluated in logarithmic form. In addition, an adaptive step size strategy is adopted to regulate the proposal scale according to the current acceptance rate. The step size is automatically adjusted to maintain the acceptance rate within a target interval of 0.2 to 0.4, thereby improving sampling efficiency, reducing manual tuning, and enhancing convergence stability.

Through this surrogate-assisted inference strategy, the number of expensive finite element evaluations is restricted to the offline training stage, whereas the online Bayesian updating procedure is conducted using the computationally efficient surrogate. This substantially reduces the overall computational burden and makes Bayesian calibration feasible for nonlinear damage models with repeated likelihood evaluations.

To provide a clear overview of the complete computational procedure, a detailed flowchart of the proposed surrogate-assisted Bayesian calibration methodology is presented in Fig. 4.

Figure 4: Flowchart of the proposed surrogate-assisted Bayesian calibration methodology.

3  Case Study of Reinforced Concrete Columns

This section presents a case study of reinforced concrete columns to demonstrate the applicability of the proposed machine learning surrogate-based Bayesian calibration framework. The experimental program is used as the benchmark, and the corresponding nonlinear finite element model is established to reproduce the structural response. Experimental observations are then incorporated to calibrate the constitutive parameters of the concrete damage model.

The implementation procedure consists of four main parts. First, a finite element benchmark model is established for the selected column specimens. Second, a surrogate model is constructed from a dataset generated by finite element simulations. Third, sensitivity analysis is performed to examine the relative influence of uncertain parameters on the structural response. Finally, Bayesian calibration is carried out to infer the posterior distributions of the constitutive parameters, with and without explicit consideration of model bias.

3.1 Finite Element Benchmark Model

The applicability of the proposed method was examined through simulation of reinforced concrete column experiments [38]. Eight groups of rectangular columns were tested under low cycle reversed loading. Among them, specimens T1, T2 and T4 from the same test series were selected for Bayesian calibration, while T3 was excluded due to premature failure. These specimens share the same main structural and loading conditions, and mainly differ in transverse reinforcement anchorage details. The concrete compressive strength was 25.6 MPa, and the yield strengths of longitudinal and transverse reinforcement were 474 and 333 MPa, respectively. The specimen dimensions, reinforcement details, and corresponding experimental force-displacement curves are shown in Fig. 5. Each column was fixed at the base and free at the top, subjected first to axial load and then to displacement controlled loading up to 100 mm.

Figure 5: Column case study (a) loading procedure and reinforcement layout; (b) experimental force-displacement curves.

The experimental results exhibit five distinct stages: linear elastic response, cracking accompanied by stiffness reduction, yielding with peak load attainment, softening with rapid crack development and structural degradation, and a residual stage, where load capacity stabilizes at a relatively low level despite continued deformation.

A fiber beam based finite element model was established using the same concrete constitutive model described in Section 2.1. The concrete section was modeled as a rectangular beam section, while the reinforcement was represented by an equivalent box section based on the principles of equal area and equal position. The model employed B31 beam elements [39] and was discretized into five segments along the height. Consistent with the experimental loading protocol described above, the column base was fixed, the top was left free, and the specimen was first subjected to axial load followed by displacement-controlled lateral loading up to 100 mm.

Concrete behavior was described using the damage plasticity model with confinement effects incorporated through equivalent transverse confining stress [31]. The initial material parameters of the damage model were taken as the prior distribution means listed in Section 3.2. The reinforcement was modeled using a bilinear constitutive relation with a hardening coefficient of 0.01. This model was used to generate the simulation dataset required for surrogate construction and subsequent Bayesian calibration.

3.2 Random Variables and Prior Distributions

Eight parameters in the concrete constitutive model were treated as uncertain random variables in the Bayesian calibration framework. Their prior distributions were defined based on available experimental data and reported values in literature [40,41]. Normal prior distributions were assigned to all parameters as summarized in Table 3. Since these parameters must remain positive on physical grounds, samples taking negative values were excluded during sampling.

The uncertain parameter vector was selected according to the physical meaning of the bi-scalar damage constitutive model, rather than as purely numerical tuning variables. Although some parameters correspond to intrinsic material properties, their effective values in structural-scale simulations may still be affected by material heterogeneity, specimen variability, scale effects, and modeling idealizations. Therefore, they are treated as uncertain constitutive parameters within physically constrained prior ranges.

For the elastic modulus, a coefficient of variation of 0.077 was adopted. For the tensile and compressive strengths, a coefficient of variation of 0.18 was used. The prior means of εt and αt were determined from the tabulated parameters of the uniaxial tensile stress-strain curve in GB 50010 [35] according to ft. Similarly, the prior means of εc and αc were determined from the tabulated parameters of the uniaxial compressive stress-strain curve according to fc. The prior mean of Fl was calculated using the Mander confinement model. For these parameters, moderate coefficients of variation were adopted based on engineering judgment because these parameters are not directly measured in standard material tests. The selected prior ranges were also used to generate the finite element simulation dataset for surrogate training. In this way, the subsequent surrogate construction and Bayesian calibration were performed consistently within the same admissible parameter domain.

3.3 Training Dataset and Surrogate Validation

To construct the surrogate model described in Section 2.2, LHS was performed over the parameter space. A total of 2000 parameter sets were generated and evaluated through finite element simulations, and the resulting dataset is illustrated in Fig. 6. For each sampled parameter vector, the force-displacement response at the column top was obtained from the finite element benchmark model. Therefore, in this study, the structural response used for surrogate construction specifically refers to the column-top force-displacement response, where the output components are the reaction forces corresponding to prescribed displacement points. The multi-output force-displacement response was reconstructed by training independent SVR models at different displacement points and assembling the predicted reaction forces along the displacement axis.

Figure 6: Distribution of parameters sampled for finite element simulation of column specimens.

Based on the generated data, SVR surrogate models were established to approximate the nonlinear mapping from the constitutive parameters to the structural response. Consistent with the SVR setting described in Section 2.2, the models were built using the l2 - ε loss function. The optimized hyperparameters include the penalty parameter, the insensitive loss parameter, and the anisotropic kernel scale parameters associated with the eight input variables. Both input and output scaling were enabled, and the input parameters and output responses were standardized using the mean and standard deviation of the training data. Among the 2000 samples, 1800 were used for training and the remaining 200 were used for testing, and the same split was used for all displacement points.

The adequacy of the dataset size was examined through learning curves constructed from separate training and test sets, as shown in Fig. 7. The results indicate that the learning curve becomes stable after approximately 1500 samples, confirming that the adopted split of 1800 training samples and 200 testing samples is sufficient and that no evident overfitting is observed.

Figure 7: Surrogate dataset and validation in the column case (a) finite element simulation dataset; (b) learning curve for train-test split.

The predictive performance of the surrogate models was further evaluated by NRMSE, MAPE and R2. The force-displacement response was further divided into different critical stages for stage-wise validation, as summarized in Table 4. The trained surrogate models are therefore considered suitable for the subsequent sensitivity analysis and Bayesian calibration. In addition, the same scale comparison with the observational noise and model bias indicates that the SVR approximation error is relatively small. Therefore, the SVR surrogate models were used as deterministic computational accelerators in the subsequent Bayesian calibration.

3.4 Sensitivity Analysis

Before performing Bayesian calibration, sensitivity analysis was conducted based on the LHS and the established SVR surrogate model. Sobol sensitivity indices were calculated to quantify the relative influence of the eight uncertain parameters on the structural response [42]. In the sensitivity analysis, 2000 parameter samples were used, and the structural reaction force at each displacement point was taken as the output response for evaluating the Sobol sensitivity indices along the force-displacement curve. The results are shown in Fig. 8.

Figure 8: Sensitivity analysis results in the column case (a) first-order sensitivity index; (b) total sensitivity index.

Among the eight parameters, the most influential ones are the concrete compressive strength, the peak compressive strain, and the stirrup confinement-related parameter, whereas the remaining variables also contribute to the response evolution to varying degrees. In the elastic stage, the elastic modulus and compressive strength mainly govern the structural stiffness and load response. As the displacement increases, the compressive strength becomes increasingly dominant and largely determines the ultimate load capacity. In the plastic stage, the nonlinear behavior of concrete becomes more pronounced, and the influence of stirrup confinement is significantly enhanced, leading to increased toughness and improved crack resistance.

Overall, the sensitivity evolution reflects the transition from elastic response to nonlinear deterioration, as well as the changing interaction between concrete and reinforcement. These results confirm that the selected parameters are identifiable from the structural response and provide support for the subsequent Bayesian calibration.

3.5 Bayesian Calibration Results

Bayesian calibration was performed to infer the posterior distributions of the eight constitutive parameters and quantify the associated uncertainty in concrete behavior. Four parallel MCMC chains with 500,000 samples per chain were generated. A Gaussian random walk proposal with an adaptive step size was adopted. A burn-in length of 40% and a thinning interval of 250 were used, resulting in 4800 retained posterior samples in total. The average acceptance rate of the four chains was 0.35. Two calibration processes were conducted, namely calibration without model bias and calibration with model bias.

The trace plots shown in Fig. 9 exhibit stable mixing behavior without evident drift. This was further supported by the Gelman Rubin statistics listed in Table 5, while the effective sample sizes (ESS) were used to quantify the sampling efficiency. Although some parameters exhibited relatively small ESS values due to posterior autocorrelation, the Gelman Rubin statistics remained close to 1, indicating satisfactory convergence. In addition to posterior stability, substantial computational acceleration was achieved through the proposed surrogate-assisted framework. The computational cost was reduced from about 10 min per finite element simulation to less than 1 s per surrogate prediction, corresponding to a speedup of at least 600 times.

Figure 9: Trace plots for the four chains of the eight parameters in the column case.

The pairwise correlation plots in Fig. 10 show three main categories of relationships among the eight parameters, namely strong correlation, weak correlation, and nonlinear dependence. In particular, the strength related parameters and the corresponding peak strain parameters exhibit strong positive correlations, indicating their coupled influence on the load carrying capacity of concrete. By contrast, several parameter pairs show bent and asymmetric contours, indicating nonlinear coupling effects related to tensile softening behavior.

Figure 10: Pairwise correlation plot of eight parameters in the column case. The diagonal panels show marginal posterior distributions, while the off-diagonal panels present bivariate posterior density contours. Warmer colors (red) denote higher probability density, while cooler colors (blue) denote lower density.

The posterior distributions compared with priors are shown in Fig. 11. In all cases, the posterior distributions become more concentrated after calibration, indicating that the experimental data provide informative updating of the constitutive parameters. A comparison between the two calibration schemes shows that the posterior distributions obtained without model bias are narrower for some parameters, whereas those obtained with model bias are generally broader and more conservative. In particular, explicit treatment of model bias increases the conservativeness of some parameter estimates and improves the robustness of the inferred uncertainty range, while its influence on other parameters remains limited. The corresponding posterior statistics are summarized in Table 6.

Figure 11: Comparison of prior and two posterior distributions of eight parameters in the column case.

It should be noted that the eight calibrated parameters are identifiable to different degrees from the global force-displacement response. The Sobol sensitivity analysis indicates that the global response provides stronger information for response dominant parameters, such as the compressive strength, peak compressive strain, and stirrup confinement related parameter. Other parameters can still be updated through Bayesian calibration, but their posterior uncertainty reduction may be less pronounced due to weaker sensitivity or stronger posterior correlation. Additional local physical measurements, such as strain, crack development, or local deformation responses, would be helpful for further improving parameter identifiability.

The posterior response in Fig. 12 was obtained using the experimental curves included in the Bayesian calibration. The calibration results are shown through the 95% credible interval of the posterior response. Good agreement is observed between the posterior mean predictions and experimental data. The improvement is mainly driven by the updating of the most sensitive parameters, while the introduction of model bias further reduces the discrepancy between simulation and experiment. Overall, the proposed Bayesian calibration framework effectively refines parameter distributions, improves the agreement between model predictions and calibration data, and provides a more reliable uncertainty representation for nonlinear structural responses.

Figure 12: Credible interval with posterior distribution in the column case.

4  Conclusion

This study presents a machine learning surrogate-based Bayesian calibration framework for nonlinear concrete damage models. The framework combines SVR-based surrogate modeling, response-level model bias representation and adaptive MCMC sampling within a unified probabilistic inference procedure. By replacing repeated finite element evaluations with a trained surrogate model during posterior sampling, the proposed method substantially improves computational efficiency while retaining predictive accuracy.

The reinforced concrete column case study demonstrates that proposed framework effectively calibrates constitutive parameters from experimental force-displacement data. The surrogate model accurately reproduces the nonlinear responses generated by the finite element benchmark model, and Bayesian updating reduces prior uncertainty in the calibrated parameters. Sensitivity analysis further identifies the compressive strength, peak compressive strain, and stirrup confinement related parameter as the dominant variables governing the structural response.

The comparison between calibration results with and without model bias shows that explicit bias treatment yields broader but more realistic posterior uncertainty, whereas neglecting model bias leads to comparatively overconfident parameter estimates. The close agreement between the posterior responses and the experimental observations used for calibration further demonstrates the effectiveness of the proposed framework within the present calibration setting.

Overall, the proposed method provides an efficient and robust strategy for parameter calibration and uncertainty quantification in nonlinear concrete damage models. Nevertheless, the present case study is mainly applicable to reinforced concrete column and beam specimens modeled using fiber beam elements and calibrated against global force-displacement responses. When the framework is extended to shear walls, three-dimensional solid models, or full cyclic hysteretic responses, additional experimental observations, more refined finite element descriptions, and higher dimensional parameter representations may be required. In particular, local physical information, such as strain, stress and crack development, can help improve parameter identifiability and better constrain the calibration process. Future work should therefore focus on incorporating multi source response information and extending the framework to more complex structural systems, loading conditions, and surrogate construction strategies.

Acknowledgement: Not applicable.

Funding Statement: This research was funded by Shanghai Municipal Commission of Science and Technology (Grant No. 24TS1416400) and National Natural Science Foundation of China (Grant No. 52478198).

Author Contributions: The authors confirm contribution to the paper as follows: Conceptualization, methodology and writing—original draft preparation: Yi Chen; supervision, funding acquisition and writing—review and editing: Xiaodan Ren. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data available on request from the authors.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Teigen JG, Frangopol DM, Sture S, Felippa CA. Probabilistic FEM for nonlinear concrete structures. I: theory. J Struct Eng. 1991;117(9):2674–89. doi:10.1061/(asce)0733-9445(1991)117:9(2674). [Google Scholar] [CrossRef]

2. Genikomsou AS, Polak MA. Finite element analysis of punching shear of concrete slabs using damaged plasticity model in ABAQUS. Eng Struct. 2015;98(7):38–48. doi:10.1016/j.engstruct.2015.04.016. [Google Scholar] [CrossRef]

3. Husek M, Kala J. Material structure generation of concrete and its further usage in numerical simulations. Struct Eng Mech. 2018;68(3):335–44. [Google Scholar]

4. Beck JL, Katafygiotis LS. Updating models and their uncertainties. I: bayesian statistical framework. J Eng Mech. 1998;124(4):455–61. doi:10.1061/(asce)0733-9399(1998)124:4(455). [Google Scholar] [CrossRef]

5. Walters DJ, Biswas A, Lawrence EC, Francom DC, Luscher DJ, Fredenburg DA, et al. Bayesian calibration of strength parameters using hydrocode simulations of symmetric impact shock experiments of Al-5083. J Appl Phys. 2018;124(20):205105. doi:10.1063/1.5051442. [Google Scholar] [CrossRef]

6. Gao J, Du K, Qi J. Quantifying the uncertainty of structural parameters using machine learning–based surrogate models. ASCE-ASME J Risk Uncertain Eng Syst Part A Civ Eng. 2025;11(2):04025023. doi:10.1061/ajrua6.rueng-1550. [Google Scholar] [CrossRef]

7. Song Z, Xu W, Valdebenito MA, Faes MGR. Efficient forward and inverse uncertainty quantification for dynamical systems based on dimension reduction and Kriging surrogate modeling in functional space. Mech Syst Signal Process. 2025;235:112898. doi:10.1016/j.ymssp.2025.112898. [Google Scholar] [CrossRef]

8. Chen B, Shen L, Zhang H. Gaussian process regression-based material model for stochastic structural analysis. ASCE-ASME J Risk Uncertain Eng Syst Part A Civ Eng. 2021;7(3):04021025. [Google Scholar]

9. Funk S, Airoud Basmaji A, Nackenhorst U. Globally supported surrogate model based on support vector regression for nonlinear structural engineering applications. Arch Appl Mech. 2023;93(2):825–39. doi:10.1007/s00419-022-02301-3. [Google Scholar] [CrossRef]

10. Ni P, Xia Y, Li J, Hao H. Using polynomial chaos expansion for uncertainty and sensitivity analysis of bridge structures. Mech Syst Signal Process. 2019;119:293–311. doi:10.1016/j.ymssp.2018.09.029. [Google Scholar] [CrossRef]

11. Li W, Liu Z, Ma Y, Liu W, Meng Z, Ma J, et al. Sensitivity analysis of structural dynamic behavior based on the sparse polynomial chaos expansion and material point method. Comput Model Eng Sci. 2025;142(2):1515–43. doi:10.32604/cmes.2025.059235. [Google Scholar] [CrossRef]

12. Guan Y, Zhi P, Wang Z. A novel variable-fidelity Kriging surrogate model based on global optimization for black-box problems. Comput Model Eng Sci. 2025;144(3):3343–68. doi:10.32604/cmes.2025.069515. [Google Scholar] [CrossRef]

13. Faroz SA, Ghosh S, Pushkaran T. Optimum calibration of a corrosion rate instrument using information gain criterion within a Bayesian framework. Struct Saf. 2023;104(1):102354. doi:10.1016/j.strusafe.2023.102354. [Google Scholar] [CrossRef]

14. Kučerová A, Sýkora J, Havlásek P, Jarušková D, Jirásek M. Efficient probabilistic multi-fidelity calibration of a damage-plastic model for confined concrete. Comput Meth Appl Mech Eng. 2023;412(2):116099. doi:10.1016/j.cma.2023.116099. [Google Scholar] [CrossRef]

15. Rossat D, Baroth J, Briffaut M, Dufour F, Monteil A, Masson B, et al. Bayesian inference with correction of model bias for Thermo-Hydro-Mechanical models of large concrete structures. Eng Struct. 2023;278(5):115433. doi:10.1016/j.engstruct.2022.115433. [Google Scholar] [CrossRef]

16. Simwanda L, Babafemi AJ, De Koker N, Viljoen C. Bayesian calibration and reliability analysis of ultra high-performance fibre reinforced concrete beams exposed to fire. Struct Saf. 2023;103:102352. doi:10.1016/j.strusafe.2023.102352. [Google Scholar] [CrossRef]

17. Yuan Z, Li Q, Li K. Measurement plan targeting the accuracy of calibrated chloride ingress model for concrete structures in marine environment. Struct Saf. 2024;106(6):102405. doi:10.1016/j.strusafe.2023.102405. [Google Scholar] [CrossRef]

18. Arendt PD, Apley DW, Chen W. Quantification of model uncertainty: calibration, model discrepancy, and identifiability. J Mech Des. 2012;134(10):100908. doi:10.1115/1.4007390. [Google Scholar] [CrossRef]

19. Hegde A, Weiss E, Windl W, Najm HN, Safta C. A Bayesian calibration framework with embedded model error for model diagnostics. Int J Uncertain Quantif. 2024;14(6):37–70. doi:10.1615/int.j.uncertaintyquantification.2024051602. [Google Scholar] [CrossRef]

20. Ghahari F, Sargsyan K, Taciroglu E. Quantification of modeling uncertainty in the Rayleigh damping model. Earthq Engng Struct Dyn. 2024;53(9):2950–6. doi:10.1002/eqe.4143. [Google Scholar] [CrossRef]

21. Reiser P, Aguilar JE, Guthke A, Bürkner PC. Uncertainty quantification and propagation in surrogate-based Bayesian inference. Stat Comput. 2025;35(3):66. doi:10.1007/s11222-025-10597-8. [Google Scholar] [CrossRef]

22. Taflanidis AA, Aakash BS, Yi SR, Conte JP. Surrogate-aided Bayesian calibration with adaptive learning strategies. Mech Syst Signal Process. 2025;237:113014. doi:10.1016/j.ymssp.2025.113014. [Google Scholar] [CrossRef]

23. Chen J, Ren X, Feng DC, Kohler J, Sørensen JD, Wu JY, et al. Recent developments in mechanical and uncertainty modelling of concrete. Struct Saf. 2025;113(6):102526. doi:10.1016/j.strusafe.2024.102526. [Google Scholar] [CrossRef]

24. Dougill JW. On stable progressively fracturing solids. Z Für Angew Math Und Phys ZAMP. 1976;27(4):423–37. doi:10.1007/BF01594899. [Google Scholar] [CrossRef]

25. Jie L. A review on the constitutive model for static and dynamic damage of concrete. Adv Mech. 2010;40(3):284–98. [Google Scholar]

26. Kachanov LM. Rupture time under creep conditions. Int J Fract. 1999;97(1):11–8. doi:10.1023/A:1018671022008. [Google Scholar] [CrossRef]

27. Li J, Ren X. Stochastic damage model for concrete based on energy equivalent strain. Int J Solids Struct. 2009;46(11–12):2407–19. doi:10.1016/j.ijsolstr.2009.01.024. [Google Scholar] [CrossRef]

28. Mazars J, Pijaudier-Cabot G. Continuum damage theory—application to concrete. J Eng Mech. 1989;115(2):345–65. doi:10.1061/(asce)0733-9399(1989)115:2(345). [Google Scholar] [CrossRef]

29. Li J, Zeng S, Ren X. A stochastic rate dependent damage model for concrete. J Tongji Univ. 2014;42(12):1783–9. [Google Scholar]

30. GB 50010-2010. Code for design of concrete structures. Beijing, China: China Architecture & Building Press; 2010. [Google Scholar]

31. Mander JB, Priestley MJN, Park R. Theoretical stress-strain model for confined concrete. J Struct Eng. 1988;114(8):1804–26. doi:10.1061/(asce)0733-9445(1988)114:8(1804). [Google Scholar] [CrossRef]

32. Loh WL. On Latin hypercube sampling. Ann Statist. 1996;24(5):2058–80. doi:10.1214/aos/1069362310. [Google Scholar] [CrossRef]

33. McKay MD, Beckman RJ, Conover WJ. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics. 2000;42(1):55–61. doi:10.1080/00401706.2000.10485979. [Google Scholar] [CrossRef]

34. Vapnik V. The nature of statistical learning theory. New York, NY, USA: Springer Science & Business Media; 2013. [Google Scholar]

35. Smola AJ, Schölkopf B. A tutorial on support vector regression. Stat Comput. 2004;14(3):199–222. doi:10.1023/B:STCO.0000035301.49549.88. [Google Scholar] [CrossRef]

36. Hastings WK. Monte Carlo sampling methods using Markov chains and their applications. Biometrika. 1970;57(1):97. doi:10.2307/2334940. [Google Scholar] [CrossRef]

37. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. Equation of state calculations by fast computing machines. J Chem Phys. 1953;21(6):1087–92. doi:10.1063/1.1699114. [Google Scholar] [CrossRef]

38. Tanaka H. Effect of lateral confining reinforcement on the ductile behaviour of reinforced concrete columns [dissertation]. Canterbury, New Zealand: University of Canterbury; 1990. [Google Scholar]

39. Timoshenko SP. On the correction for shear of the differential equation for transverse vibrations of prismatic bars. Phil Mag Ser. 1921;41(245):744–64. doi:10.1080/14786442108636264. [Google Scholar] [CrossRef]

40. Du X-L, Jin L. Research on the heterogeneous statistical properties of elastic modulus of a concrete meso-scale unit. Eng Mech. 2012;29(10):106–15. [Google Scholar]

41. Ellingwood BR, Galambos TV, MacGregor JG. Development of a probability based load criterion for american national standard a58: building code requirements for minimum design loads in buildings and other structures. Washington, DC, USA: US Department of Commerce, National Bureau of Standards; 1980. [Google Scholar]

42. Sobol’ IYM. On sensitivity estimation for nonlinear mathematical models. Mat Model. 1990;2(1):112–8. [Google Scholar]