<figure><img src="image_1.png" /></figure>

HAL
open science

# Inverse analysis of the impact response of a rockfall protection structure: Application towards warning and survey

Ritesh Gupta, Franck Bourrier, Stéphane Lambert

► To cite this version:

Ritesh Gupta, Franck Bourrier, Stéphane Lambert. Inverse analysis of the impact response of a rockfall protection structure: Application towards warning and survey. 2024. hal-04582717

HAL Id: hal-04582717

https://hal.science/hal-04582717

Preprint submitted on 22 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Bayesian inference based inverse analysis of the impact response of a rockfall protection structure: Application towards warning and survey

Ritesh Gupta, Franck Bourrier, Stéphane Lambert

► To cite this version:

Ritesh Gupta, Franck Bourrier, Stéphane Lambert. Bayesian inference based inverse analysis of the impact response of a rockfall protection structure: Application towards warning and survey. Engineering Structures, 2024, 321, pp.118800. 10.1016/j.engstruct.2024.118800. hal-04582717v2

HAL Id: hal-04582717

https://hal.science/hal-04582717v2

Submitted on 6 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Inverse analysis of the impact response of a rockfall protection structure: Application towards warning and survey

Ritesh Guptaª, Franck Bourrierª,ᵇ, Stéphane Lambertª,\*

ªUniv. Grenoble Alpes, INRAE, CNRS, IRD, Grenoble INP, IGE, 38000 Grenoble, France

ᵇUniv. Grenoble Alpes, Inria, CNRS, Grenoble INP, Institute of Engineering, LJK, 38000 Grenoble, France

## Abstract

The viability of a rockfall protection structure is vital for the hazard mitigation of habitations and infrastructures. This article investigates the feasibility and potential of inverse analysis applied to data collected on on-site rockfall protection structures exposed to real events. As an application case, a rockfall protection wall made from interconnected concrete blocks which are piled up in a zig-zag pattern is considered. The numerical model of this structure is created via the use of a Non-Smooth Contact Dynamics (NSCD) method-based modelling technique. The NSCD model is developed using the SICONOS software and spatiotemporally calibrated from two real-scale impact experiments with 520kJ and 1020kJ projectile energy. This model is used to investigate the variability in wall mechanical response against different impact conditions. The simulation results served as input data for developing the inverse analysis method. As a first application, it is proposed to use the inverse analysis to aid in remote decision-making shortly after an event, based on real-time measurements. Then, the use of inverse analysis to retrieve the impact condition characteristics (energy, location) from data collected after the event is addressed. The proposed approach appeared efficient for back-analysing (i.e., output to input) data related to the wall response for being used as a warning based on its displacement with respect to the protected element at risk and damage to the wall with root mean square errors (RMSE) of 16 cm and 82 kJ respectively and for a rockfall site survey with RMSE of 71 kJ energy transferred to the wall and position of impact with 37 cm.

**Keywords:** rockfall, impact, meta-model, statistical learning, NSCD modelling, inverse analysis

<footnote>\*Corresponding author
Email address: stephane.lambert@inrae.fr (Stéphane Lambert)</footnote>
# Bayesian inference based inverse analysis of the impact response of a
rockfall protection structure: Application towards warning and survey

Ritesh Guptaª, Franck Bourrierª,ᵇ, Stéphane Lambertª,\*

## Abstract

The viability of a rockfall protection structure is vital for the hazard mitigation of habitations and infrastructures. This article investigates the feasibility and potential of inverse analysis applied to data collected on on-site rockfall protection structures exposed to real events. As an application case, a rockfall protection wall made from interconnected concrete blocks which are piled up in a zig-zag pattern is considered. The corresponding numerical model is developed in a python-based open source software SICONOS which implements the Non-Smooth Contact Dynamics (NSCD) method. The numerical model was previously spatiotemporally calibrated from two real-scale impact experiments with 520kJ and 1020kJ projectile energy and used in this work to investigate the variability in wall mechanical response against different impact conditions. The simulation results served as input data to develop the Bayesian inference based inverse analysis method aided by the polynomial chaos expansion based metamodelling technique for two purposes. First, to aid in remote decision-making shortly after an event, based on real-time measurements and second, to retrieve the impact condition characteristics (energy, location) from data collected after the event is addressed. The proposed approach appeared efficient for back-analysing (i.e., output to input) data related to the wall response for being used as a warning based on its displacement with respect to the protected element at risk and damage to the wall with root mean square error (RMSE) of 16 cm and 82 kJ respectively and for a rockfall site survey with RMSE of 71 kJ energy transferred to the wall and position of impact with 37 cm.

**Keywords:** rockfall, impact, meta-model, statistical learning, NSCD modelling, inverse analysis

<footnote>ªUniv. Grenoble Alpes, INRAE, CNRS, IRD, Grenoble INP, IGE, 38000 Grenoble, France</footnote>

<footnote>ᵇUniv. Grenoble Alpes, Inria, CNRS, Grenoble INP, Institute of Engineering, LJK, 38000 Grenoble, France</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Email address: stephane.lambert@inrae.fr (Stéphane Lambert)</footnote>
1. Introduction

Passive rockfall protection mitigation measures such as embankments [1], galleries [2] or flexible barriers [3] are exposed to severe loading when intercepting rock blocks threatening elements at risk. During their normal operation, these structures thus experience significant deformation, displacement, and damage and possibly fail to withstand the impact or to satisfactorily control the rock blocks trajectory [4].

In this context, and similarly, as for other gravitational natural hazards, increasing use is being made of equipment for monitoring on-site rockfall protection structures, to collect the data upon the impact of a rockfall [5]. On-site monitoring of structures with accelerometers or force sensors, for example, in particular aim at serving as a warning system for supporting decision-making (e.g., road closure) or for providing information regarding the loading amplitude or the structure response in real situations.

On-site monitoring is becoming increasingly common in the field of landslides due to the improvement, miniaturisation and cost reduction of sensors, data acquisition and transmission systems [6, 7].

Numerical models are widely used for addressing rockfall protection structures' response as a forward problem. Published research demonstrates the accuracy of the model predictions as compared to experimental data (e.g. [8, 9, 10]) or addressing the structure response under some specific impact conditions (e.g. [11, 12]) or considering sets of realistic impact conditions (e.g. [13, 14]). By contrast, inverse analysis conducted (i.e., finding input from the output) based on numerical models have rarely been used in this field. To the best of the authors' knowledge, the only exception concerns the work presented by Escallón et al. [15] who proposed an inverse optimisation process for determining the parameters of macroscopic FE models of steel wire-rope cables and steel wire-rings.

This article investigates the feasibility of conducting an inverse analysis of the structure response for deriving information with added value for stakeholders and engineers. It is for example proposed to analyse the measured data pertinent to the structure's performance to retrieve the impact conditions. The proposed inverse analysis consists of comparing data collected on-site during or after an event with the structure impact response as obtained from numerical simulations. This comparison relies on the use of statistical learning methods enhanced by meta-modelling tools. A statistically modelled series of different impact conditions and a learning-enabled link between input and output responses provide us with the necessary components to do the inverse analysis. Further, the establishment of a mathematical link between the structure's response and the corresponding impact condition is enabled through the use
## 1. Introduction

Passive rockfall protection mitigation measures such as embankments [1], galleries [2] or flexible barriers [3] are exposed to severe loading when intercepting rock blocks threatening elements at risk. During their normal operation, these structures thus experience significant deformation, displacement, and damage and possibly fail to withstand the impact or to satisfactorily control the rock block's trajectory [4].

In this context, and similarly, as for other gravitational natural hazards, increasing use is being made of equipment for monitoring on-site rockfall protection structures, to collect the data upon the impact of a rockfall [5]. On-site monitoring of structures with accelerometers or force sensors, for example, in particular aims at serving as a warning system for supporting decision-making (e.g., road closure) or for providing information regarding the loading amplitude or the structure response in real situations.
On-site monitoring is becoming increasingly common in the field of landslides due to the improvement, miniaturisation and cost reduction of sensors, data acquisition and transmission systems [6, 7].

Numerical models are widely used for addressing rockfall protection structures' response as a forward problem. Published research demonstrates the accuracy of the model predictions as compared to experimental data (e.g. [8, 9, 10]) or addressing the structure response under some specific impact conditions (e.g. [11, 12]) or considering sets of realistic impact conditions (e.g. [13, 14]). By contrast, inverse analysis conducted (i.e., finding input from the output) based on numerical models has rarely been used in this field. To the best of the authors' knowledge, the only exception concerns the work presented by Escallón et al. [15] who proposed an inverse optimisation process for determining the parameters of macroscopic FE models of steel wire-rope cables and steel wire-rings.

In this work, we investigate the feasibility of conducting an inverse analysis of the structure response to derive the information with added value for stakeholders and engineers. It is for instance proposed to analyse the measured data pertinent to the structure's performance to retrieve the impact conditions. The proposed inverse analysis presents means to correlate the data retrievable from an on-site structure during and after an rockfall event, with the corresponding impact response obtained from the numerical model simulation of that structure. This correlation relies on the use of statistical learning methods enhanced by meta-modelling tools. A statistically modelled series of different impact conditions and a learning-enabled link between input and output responses provide us with the necessary components to do
of meta-models having negligible computation costs.

Here, we focus on an innovative passive rockfall protection structure made of pile-up concrete blocks to form a wall with a zig-zag pattern, for which a Non-Smooth Contact Dynamics (NSCD) method-based numerical model was previously developed. The calibration of this model was conducted against experimental data following a complex method involving statistical learning, resulting in improved confidence in the model's predictive capacities, reported in [16]. In parallel, specific instrumentation including sensors (i.e. accelerometers and inclinometers), remote data acquisition and processing has been developed to equip these structures. In principle, the proposed approach for conducting an inverse analysis could be applied to any type of structure, and in particular structures exposed to dynamic loading such as flexible barriers.

The paper is organised as follows: First, the NSCD model definition and its calibration are briefly presented. Then the mechanical response of the wall to impact is addressed thanks to NSCD model simulations. First, varying the impact location and then via six parameters collectively describing the rock block trajectory, kinematics and impact point on the wall. The impact simulation results highlight the intricacy of the structure response and suggest the simulation outputs which best allow differentiating the impact cases. Then, a method for conducting the inverse analysis of the collected data is developed to support the real-time decision and rockfall activity monitoring purpose of a given site.

## 2. Structure and model presentation

The considered structure is a composition of concrete blocks, reinforced with internal rebars and placed together in a staggered pattern to form a wall whose geometry can be adapted to specific requirements (e.g. linear or zig-zag conformation along its longitudinal axis). The blocks are connected thanks to metallic tubes and cables passing through vertical cylindrical holes in the blocks, providing mechanical continuity to the structure [17]. There exist mechanical plays between concrete blocks of the same row as well as between the metallic tubes and the blocks.

This type of articulated structure constitutes an alternative to other massive passive rockfall protection structures with reduced footprint [18, 19] with the advantages of high versatility in the design and deformability under impact. For such an application where the structure is exposed to a localised dynamic loading, tubes and cables increase the number of blocks involved in the structure's response proportional
the inverse analysis. Further, the establishment of a mathematical link between the structure's response and the corresponding impact condition is enabled through the use of meta-models having negligible computation costs.

In this work, we focus on an innovative passive rockfall protection structure made of pile-up concrete blocks to form a wall with a zig-zag pattern, for which a Non-Smooth Contact Dynamics (NSCD) method-based numerical model was previously developed. The calibration of this model was conducted against experimental data following a complex method involving statistical learning, resulting in improved confidence in the model's predictive capacities, reported in Gupta et al. [16]. In parallel, specific instrumentation including sensors (i.e. accelerometers and inclinometers), remote data acquisition and processing has been developed to equip these structures. In principle, the proposed approach for conducting an inverse analysis could be applied to any type of structure, and in particular structures exposed to dynamic loading such as flexible barriers.

The article is organised as follows: First, the type of the protection structure considered in this study is described followed by the presentation of its numerical equivalent developed in the non-smooth contact dynamics (NSCD) framework - named the NSCD model. A brief description of the previously conducted NSCD model calibration using the displacement data from the impact response for two full-scale structures is presented alongside. Then the mechanical response of the calibrated NSCD model wall to the model projectile impact is addressed via an extensive campaign of simulations. An illustrative structural response is first presented at three distinct sections of wall geometry followed by the generalised impact assessment via varying the six parameters collectively describing the rock block trajectory, kinematics and impact point on the wall. The impact simulation results highlight the intricacy of the structure response and suggest the simulation outputs which best allow differentiating the impact cases. Then, a method for conducting the inverse analysis from the acquired data is developed using a combination of meta-modelling and Bayesian inference techniques to support the real-time decision-making and rockfall activity monitoring purposes for a given site. A discussion on of predictive capability of the inverse analysis method, its limitations and the perspective on improvements towards future research concludes the work.
to the impact loading.

The dynamic response of this structure is evaluated thanks to two full-scale impact experiments with energies of 520kJ and 1020kJ on a wall 3.2 m in height and 14.1 m in length (Figure 1a). The projectile used for this purpose is 2600kg in mass and 1.1 m in dimension, corresponding to approximately 1/3rd the wall height, conforming with the requirements of flexible barrier testing [20]. The wall displacement evolution data is recorded at a couple of points on the impact axis (at points Top and Base) and on the distant axis (at points C and D) as further detailed in Furet et al. [17].

<figure><img src="image_3.png" /><figcaption>Figure 1: (a) The experimental real-scale wall structure and (b) its numerical equivalent developed using SICONOS software based on the framework of non-smooth contact dynamics (NSCD) method, hereby referred to as NSCD model.</figcaption></figure>

Based on these impact experiments, a numerical model of the structure was developed under the Non-Smooth Contact Dynamics (NSCD) framework implemented in SICONOS [21] software. The NSCD model of the structure is a collection of blocks and connectors (Figure 1b). The detailed description of the model development and calibration is presented in Gupta et al. [16].

The same block geometry is considered in the model as that in the real structure. The cables and slings are modelled with connectors. The interaction between the structure model components, as well as the interaction with the projectile, abides with the framework of the NSCD method.
## 2. Experimental Structure and NSCD model

### 2.1. Presentation

The considered structure is a composition of concrete blocks, reinforced with internal rebars and placed together in a staggered pattern to form a wall whose geometry can be adapted to specific requirements (e.g., linear or zig-zag conformation along its longitudinal axis). Each block weighs about 1800 kg, is cuboid in shape with dimensions 1.56m×0.76m×0.8m having rounded extremities in the horizontal plane and is traversed by two cylindrical holes of diameter 154 mm along the vertical axis. The blocks are connected thanks to metallic tubes and cables passing through these holes, providing mechanical continuity to the structure [17]. There exist mechanical plays between concrete blocks of the same row as well as between the metallic tubes and the blocks.

This type of articulated structure constitutes an alternative to other massive passive rockfall protection structures with reduced footprint [18, 19] with the advantages of high versatility in the design and deformability under impact. For such an application where the structure is exposed to a localised dynamic loading, tubes and cables increase the number of blocks involved in the structure's response proportional to the impact loading.

The dynamic response of this structure is evaluated thanks to two full-scale impact experiments carried out at the pendulum testing facility of the Université Gustave Eiffel test site (Montagnole, France). This testing facility enables the control of the impact position and constraining the projectile incidence angle and rotational velocity to zero. The blocks are arranged following a zig-zag pattern (at 45° angle) to improve the structure stability against tilting. The wall is made up of 38 blocks and four half-blocks stacked in four layers resulting in 3.2 m height and 14.1 m length, presented in Figure 1a. Notably, the layers from the ground up are named ‘base’, ‘second’, ‘third’ and ‘top’ for all future mentions. The projectile used for this purpose is 2600kg in mass and 1.1 m in dimension, corresponding to approximately $1/3^{rd}$ the wall height, conforming with the requirements of flexible barrier testing [20]. The released projectile impacted the wall at its mid-length and at a height of about 1.7 m from the ground. The velocity at impact was 20m/s and 28m/s during the first and second impacts, corresponding to kinetic energies of 520 kJ and 1020 kJ respectively.

The wall displacement evolution data is recorded at a couple of points on the impact axis (at points Top and Base) and on the distant axis (at points C and D). Besides, the data from two points A and B
The model parameters governing the interactions between the various bodies of the system were calibrated against the experimental database describing the spatio-temporal response of the wall during the two impact tests. Notably, the calibration was based on the displacement measured with time (at three-time instants) at the four locations previously mentioned. The model calibration was supported by the Bayesian statistical learning method, improving the confidence in the derived set of model parameters, and thereby enabling the replication of the structure response with time and space, over a range of impact energies up to 1 MJ. In addition, the model being simple and developed under the NSCD method framework provides on average about 30 times faster computation times in comparison to conventional FE-based models [17] and allows envisaging running a large number of computations.

In the following, all simulations are conducted on a wall 3.2m high and about 28.2m long (Fig. 2). The increase in length as compared to that in Gupta et al. [16] aims at avoiding any potential boundary effects when varying the impact location along the wall longitudinal axis. For these simulations, the same projectile is considered. Moreover, the comparison of spatio-temporal displacement response for both small and large wall lengths (i.e., 14.1 and 28.2m respectively) for the two calibrated impact cases are reported qualitatively and quantitatively similar. Therefore, the calibrated model constitution parameters as reported in Gupta et al. [16] are directly applied to the long wall model.

### 3. Structure response investigation

Due to its zig-zag conformation, discrete nature and complex design where concrete blocks are interconnected thanks to metallic components and including mechanical plays, the structure response is expected to be highly sensitive to impact conditions, for example varying the impact location or the projectile pre-impact trajectory.

First, the structural response variability is illustrated by considering three impact cases in similar conditions to that during the experiments but with different impact locations. Then the investigation is generalised varying 6 parameters describing the impact conditions.

#### 3.1. Response for the illustrative impact conditions

The wall response is addressed by varying the projectile impact location to illustrate its complexity and variability. Three impacts with a 520 kJ energy and located in the centre of each of the distinct
in the impact axis is also available which is used for the cross-reference purpose only. Measures along the impact axis were derived from video records while measures from cable extensometers were considered along the distant axis. Notably, due to the absence of this later data for the 520 kJ impact test at point C, the evolution obtained from simulations presented in Furet et al. [17] are considered. This exception is thought to have a minor influence due to the very small amplitude in displacement observed in this specific case. The experimentally recorded evolution of displacement with time at these four points for two impact tests can be referred to in Figure 2. Further information on the two real-scale impact experiments setup, instrumentation and data acquisition shall be referred to in Furet et al. [17].

<figure><img src="image_3.png" /><figcaption>Figure 1: (a) The experimental real-scale wall structure and (b) its numerical equivalent developed using SICONOS software based on the framework of non-smooth contact dynamics (NSCD) method, hereby referred to as NSCD model.</figcaption></figure>

Based on these impact experiments, a numerical model of the structure was developed under the Non-Smooth Contact Dynamics (NSCD) framework implemented in a python-based software package SICONOS [21]. The NSCD model of the structure is a collection of blocks and connectors, presented in Figure 1b. The detailed description of the model development and calibration is presented in Gupta et al. [16] and briefed hereafter.

The NSCD method was developed to solve multi-body multi-contact problems with rigid and/or
sections of the zig-zag pattern were simulated. These sections are referred to as convex-wedge, angled-
wedge and concave-wedge, when viewed from the impacted side, as presented in Figure 2. The wall
response is addressed in terms of block displacement, velocity, inclination and plastic damage.

<figure><img src="image_4.png" /><figcaption>Figure 2: Top view of the wall model highlighting the repetitive pattern and an identifier for different sections.</figcaption></figure>

3.1.1. *Displacement*

The displacement at the top of the wall is extremely sensitive to the impact location, presented for three illustrative impact conditions (IC), via the top view of the wall in Figure 3. Here, the three simulated impacts one each at convex, angled and concave wedges are named as IC-1, IC-2 and IC-3 respectively. This is particularly noticeable in terms of structure conformation. The impact on the concave-wedge results in aligned blocks (Fig. 3c), while an impact on the convex-wedge amplifies the amplitude of the zig-zag (Fig. 3a). The length of the movement zone along the wall longitudinal axis varies from approximately 6 m in the case of an impact on the angled-wedge to 7 m for the two other impact cases, corresponding to four and six number of top-row blocks respectively. Figure 3 also suggests a difference in the displacement of the individual blocks. Larger block displacement is observed for impacts on the concave-wedge.

As this zig-zag patterned wall aims at protecting a given element at risk it appears much more
relevant to consider the residual distance to the so-called safety line shown in Figure 3 in addition to of
the displacement of the blocks. The safety line may be defined based on the position of the protected
element at risk. It is here arbitrarily located at 1.5 m from the convex-wedge block's rear face. From
Figure 3, it can be deduced that the minimum residual distance to the safety line is much smaller in the
case of the impact in the convex-wedge (0.5 m. approx.), as compared to the two other impact cases (1.3
m approx.). With respect to this criterion, an impact on the convex-wedge thus reveals comparatively
more critical than others.
deformable bodies with rigid contact laws such as Signorini's model of unilateral contact and Coulomb's dry friction without any kind of regularization (viscous friction) or compliance. It finds many modelling applications in the field of civil engineering and geomechanics such as masonry and stone structures exposed to static, cyclic, and dynamic loading, cohesive and non-cohesive granular materials and rockfall propagation on slopes to name a few.

In the present work, the system is modelled as a collection of rigid blocks connected by unilateral constraints with Coulomb friction. The finite-freedom dynamics of rigid bodies with unilateral constraints are known to be non-smooth, in the sense that the velocities of the system possess jumps when a contact is closing with a positive relative velocity. These velocity jumps are described through the introduction of an impact law. Thereby, the NSCD method can perform the numerical time integration of the multi-body system under impacts in a dynamic condition.

The NSCD model of the experimental structure is a collection of the model objects - i.e., blocks, connectors, projectile and ground, as shown in Figure 1b. The same block and projectile geometry is considered in the model as that in the real structure and the tubes and slings are collectively modelled as connectors. In the NSCD framework, these model objects are defined as a combination of different contactors, collectively making the model equivalent to the real structure. These contactors are defined with a uniquely assigned collision group to identify them and differentiate one from another. The interaction between model bodies is assigned through Newton impact friction non-smooth law [22]. This law governs the interaction between a pair of contactors (identified through their respective collision groups) via a user-defined coefficient of friction ($\mu$) and coefficient of restitution ($e$).

## 2.2. Calibration

The model parameters governing the interactions between the various bodies of the system were calibrated against the experimental database describing the spatio-temporal response of the wall during the two impact tests. The calibration was based on the displacement measured with time (at three-time instants) at the four locations previously mentioned. These three-time instances correspond to the evolving displacement from the start of the impact at (1) the initial nearly linear slope with a representative instance of 0.1 seconds for the impact axis and 0.25 seconds for the distant axis, (2) maximum block displacement and, (3) displacement at rest i.e., at 1.0 second.

The model calibration was carried out for five parameters. The first two take into account the imposed
<figure><img src="image_5.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c)</figcaption></figure>

Figure 3: Top view of the structure before and after a 520-kJ impact located in the centre of (a) the convex-wedge (IC-1), (b) the angled-wedge (IC-2) and (c) the concave-wedge (IC-3).

### 3.1.2. Velocity

The variability in response is here illustrated by the time evolution of the velocity of the concrete blocks. For this purpose, the velocity computed at the mid-bottom of the block rear vertical face ‘$v_{face}$’ (i.e., on the side opposite to impact) is considered which represents the potential position for the sensor deployment in the in-situ structure.

Instrumentation of on-site structures is likely to consist of 1D accelerometers installed on the concrete block faces for measuring the acceleration normal to the block with time. For a simplified analysis, the desired velocity profile from the NSCD model shall also be computed normal to the block face at all time instances. For this, the ‘$v_{face}$’ vector is projected onto the unit normal vector ($\hat{n}_b$) to the block face. The computation process is presented in Equation 1.
constraints for the block-connector ($d_z$) and connector-connector ($v_p$) interaction which collectively controls the relative flexibility of the wall and the relative looseness in the sling. The other three parameters concern the interaction laws between contacting bodies in the NSCD framework. They are defined as the friction coefficient between blocks i.e., concrete-concrete ($\mu_{cc}$) and between blocks and ground i.e., concrete-soil ($\mu_{cs}$) and the restitution coefficient (e). The calibration was supported by the Bayesian inference statistical learning method, improving the confidence in the derived set of model parameters as 6.8 cm, 7.1 cm, 0.316, 0.307 and 0.222 chronologically as mentioned above. Consequently, it enables the replication of the structure response with time and space, over a range of impact energies up to 1 MJ.

In addition, the model being simple and developed under the NSCD method framework provides on average about 30 times faster computation times in comparison to conventional FE-based models [17] and allows envisaging running a large number of computations.

<figure><img src="image_4.png" /><figcaption>Figure 2: Comparison of the wall displacement response for (a) 520 kJ and (b) 1020 kJ impact tests at ‘Top’ and ‘Base’ points (for impact axis) and ‘C’ and ‘D’ points (for distant axis) recorded experimentally (exp) with its counterpart evaluated from the calibrated NSCD model of the short-wall (SW) and long-wall (LW).</figcaption></figure>

In the present work, all simulations are conducted on a wall 3.2m high and about 28.2m long (Fig. 3). The increase in length as compared to that in Gupta et al. [16] aims at avoiding any potential boundary effects when varying the impact location along the wall longitudinal axis. For these simulations, the same projectile is considered. Moreover, the comparison of spatio-temporal displacement response for both small wall (SW) and large wall (LW) lengths (i.e., 14.1 and 28.2m respectively) for the two calibrated
$$v_{\text{face}} = v_g + (p\Omega p^{-1}) \times q \qquad (1)$$

$$v_{\text{nb}} = (v_{\text{face}} \cdot \hat{n}_{\text{b}}) \hat{n}_{\text{b}}$$

where $v_g$ is the velocity at the block gravity centre, $\Omega$ is the rotational velocity, $p$ is the quaternion (presenting block's 3D orientation) and $q$ is the relative position vector of the block face (envisaged sensor location) to its centre of gravity.

The time evolution of velocity normal to block face computed for all three cases is presented in Figure 4 (all left). In addition, the peak velocity of the blocks and the time to reach this peak are plotted in Figure 4 (all right).

A general trend where the peak velocity localises at the impacted block is observed. Besides, the number of blocks experiencing rapid displacement from the impact beginning is comparatively higher when the projectile impacts the concave-wedge (Fig. 4c). Indeed, three blocks close to the impact location experience very similar curves. The time for the wall to come to rest is smaller in case of an impact on the convex-wedge with a duration of about 0.5 as compared to 0.8 seconds after an impact on the angled-wedge in particular (Fig. 4a vs Fig. 4b). Overall, the velocity of all the blocks appears to provide rich information, with high amplitude as well as significant differences from one case to the other.

The blocks with significant maximum velocities are located within a distance along the wall length ranging from 8 to 12 m approximately depending on the impact location. This observation is in line with the observation made in Figure 3. The shape of the maximum velocity pattern significantly differs from one impact case to the other. These observations suggest that the distribution of peak velocity along the wall could be specific to the impact case in a rather univocal manner.

### 3.1.3. Inclination

In addition to the displacement and velocity, the inclination of the wall blocks, i.e., their rotation around a horizontal axis, seems relevant to address as it relates to wall post-impact stability. In addition, inclination may be measured on-site with rather low-cost sensors. The post-impact inclination of all the blocks in the wall is presented in Figure 5. An impact on the concave-wedge results in block inclination over a wall length as large as 8 m (Fig. 5). The maximum value of about 17° approx. is observed at the upper row in the impact vicinity. Results concerning the impact on the convex-wedge and angled-wedge are not presented as they induced nearly zero wall inclination.
impact cases are reported qualitatively and quantitatively similar, as presented in Figure 2. Therefore, the calibrated model constitution parameters as reported in Gupta et al. [16] are directly applied to the long wall model. Notably, the SW displacements are quantitatively slightly different (up to 4cm) from the one reported in Gupta et al. [16] due to the current computations in an updated version of SICONOS.

## 3. Impact response investigation

Due to its zig-zag conformation, discrete nature and complex design where concrete blocks are interconnected thanks to metallic components and including mechanical plays, the structure response is expected to be highly sensitive to impact conditions, for example varying the impact location or the projectile pre-impact trajectory. Therefore, an extensive campaign of numerical simulations is carried out on the NSCD model of the structure to understand, evaluate and exploit the response under different impact conditions.

First, the structural response variability is illustrated by considering three impact cases in similar conditions to that during the experiments but with different impact locations. Then the investigation is generalised varying the six parameters describing the impact conditions.

### 3.1. Response for the illustrative impact conditions

The wall response is addressed by varying the projectile impact location to illustrate its complexity and variability. The projectile with an impact energy of 520kJ is simulated where the point of impact is located in the centre of each of the distinct sections of the zig-zag pattern. The 520kJ impact energy is significantly below the wall nominal capacity of 1 MJ where severe block damage, loss of mechanical continuity and very large wall displacements are observed. Therefore, the 520kJ is considered to represent the structure's serviceable impact energy levels, for which the complexity of the wall response is hereby illustrated. These sections are referred to as convex-wedge, angled-wedge and concave-wedge, when viewed from the impacted side, as presented in Figure 3. The origin of the wall structure for the 'along length position' (i.e., Y=0) is assigned at the wall centre and for the 'wall movement due to impact' (i.e., X=0) at the rear face of the wall. The wall response is addressed in terms of block displacement, velocity, deflection and plastic damage hereafter.
<figure><img src="image_8.png" /><figcaption>Figure 4: Evolution of the velocity of the concrete block in the movement zone (all left) and the corresponding mapping of the maximum velocity for blocks in all four rows i.e., top, third, second and base (all right) for the illustrative impact cases at (a) convex wedge (IC-1), (b) angled wedge (IC-2) and concave wedge (IC-3) respectively. Here, 'Y' refers to the global position of the mid-bottom of the rear face of each block along the wall length.</figcaption></figure>

Conclusively, the inclination doesn't provide sufficiently rich information to envisage relating it to the impact conditions, in particular, because some impact conditions result in null inclination angle values and hence are not deemed appropriate to explore further.
<figure><img src="image_5.png" /><figcaption>Figure 3: Top view of the wall model highlighting the repetitive pattern and an identifier for different sections. Here, the abscissa refers to the block position along the wall length where Y=0 is located at the wall centre. And, the ordinate refers to the direction of wall movement when a projectile hits the wall from the impacted side where X=0 is located at the wall's rear face.</figcaption></figure>

### 3.1.1. Displacement

The displacement of the wall is extremely sensitive to the impact location, presented for three illustrative impact conditions (IC), via the top view of the wall in Figure 4. Here, the three simulated impacts one each at convex, angled and concave wedges are named IC-1, IC-2 and IC-3 respectively. The sensitivity to the impact location is particularly noticeable in terms of structure conformation. The impact on the concave-wedge results in aligned blocks (Fig. 4c), while an impact on the convex-wedge amplifies the amplitude of the zig-zag pattern (Fig. 4a). The length of the movement zone along the wall longitudinal axis varies from approximately 6 m in the case of an impact on the angled-wedge to 7 m for the two other impact cases, corresponding to four and six top-row blocks respectively. A significant difference in the displacement of the individual blocks is also observed in Figure 4, where the displacement is evaluated at the block's centre of gravity. Larger block displacement is observed for impacts on the concave-wedge in comparison to the other two ICs.

As this zig-zag patterned wall aims at protecting a given element at risk, it appears much more relevant to consider the residual distance to the so-called 'Safety line' shown in Figure 4 in addition to the displacement of the blocks. The safety line may be defined based on the position of the protected element at risk. It is here arbitrarily located at 1.5 m from the convex-wedge block's rear face. From Figure 4, it can be deduced that the minimum residual distance to the safety line is much smaller in the case of the impact in the convex-wedge (0.5 m. approx.), as compared to the two other impact cases (1.3 m approx.). With respect to this criterion, an impact on the convex-wedge thus reveals comparatively more critical than others.
<figure><img src="image_9.png" /><figcaption>Figure 5: Inclination of all blocks in four layers of the wall (i.e., top, third, second and base) with respect to the vertical at rest for the impacts at concave wedge (IC-3) under 520kJ energy impact. The results for the impact at the convex-wedge (IC-1) and angled-wedge (IC-2) are not presented because of their negligible magnitude.</figcaption></figure>

### 3.1.4. Plastic damage

The rockfall impact on the structure causes plastic damage to the blocks, dissipating a portion of the incident projectile kinetic energy. Notably, the NSCD model is created as an accumulation of rigid bodies where the contact law controls the interaction. Despite that, the quantitative estimation of the energy dissipation mechanisms is possible following the work reported by Acary [22] in the NSCD framework. Subsequently, the relation presented in Equation 2 is implemented in the present work to compute the energy dissipation due to plasticization ($D_p$).

$$D_p = \sum_{k=0}^{N} \sum_{\alpha \in I} \frac{1}{2} (v_{N,k+1}^{\alpha} + v_{N,k}^{\alpha}) P_{N,k+1}^{\alpha} \quad (2)$$

Here, $v_N$ is the normal component of the velocity vector and $P_N$ is the normal impulse at a given time step 'k' (out of total N) at a contact point $\alpha \in I$. The details of velocity and impulse computations can be referred to at Gupta et al. [16].

The time evolution of $D_p$ for the three illustrative cases with different impact locations is presented in Figure 6. A significant contribution of plastic damage is observed at the early time instance when the projectile impacts the wall. With reference to the zoomed-in window (see Figure 6b), the damage evolves relatively fast for the impact at convex-wedge (IC-1) and reaches a constant magnitude of about
<figure><img src="image_6.png" /><figcaption>Figure 4: Top view of the displaced component of the whole length of the wall structure before and after a 520-kJ impact located in the centre of (a) the convex-wedge (IC-1), (b) the angled-wedge (IC-2) and (c) the concave-wedge (IC-3).</figcaption></figure>

### 3.1.2. Velocity

The variability in the wall's response is here illustrated by the time evolution of the velocity of the concrete blocks during the impact. For this purpose, the velocity computed at the mid-bottom of the block rear vertical face ‘$v_{face}$’ (i.e., on the side opposite to impact) is considered which represents the potential position for the sensor deployment in the in-situ structure. In this position, the instrumentation shall likely be with 1D accelerometers to capture the time evolution of the acceleration normal to the block. For a simplified analysis, the desired velocity profile from the NSCD model shall also be computed normal to the block face at all time instances. For this, the ‘$v_{face}$’ vector is projected onto the unit normal vector ($\hat{n}_{b}$) to the block face. The computation process is presented in Equation 1.
<figure><img src="image_10.png" /><figcaption>Figure 6: (a) Time evolution of the energy dissipation due to plastic damage for the three illustrative cases, with impact at the convex-wedge (IC-1), the angled-wedge (IC-2) and concave-wedge (IC-3), (b) zoom-in window between 300 and 400 kJ</figcaption></figure>

380kJ at about 0.35 seconds post-impact. The other two impact cases evolve for a relatively longer time (about 0.6 seconds) before reaching a constant value. This observation is in line with the duration over which the wall moves, as illustrated in Figure 4. Conclusively, the plastic damage is taken as a quantity of interest for structural response assessment.

Notably, the remaining part of the energy transferred by the projectile to the wall is dissipated by friction at the contacts, as detailed in [23] where it is also demonstrated that the computation scheme complies with the fundamental principle of energy conservation.

## 3.2. Response to a large set of impact conditions

Following the illustrative structural response description in the previous section, the numerical simulation framework is extended to address the response considering a large set of impact conditions. In a similar way as previously done for flexible barriers [13, 24, 25], the parameters describing the incident projectile trajectory, kinematics and location, hereafter referred to as ICPs (Impact Conditions Parameters) are varied over realistic ranges to account for nearly all possible distinct impact conditions which may occur in the field. The ICPs and their magnitude ranges are given in Table 1.

The first two parameters (i.e., translational and rotational velocities) account for the incident kinetic energy of the projectile. Given the projectile geometric and mechanical characteristics considered in this
$$v_{\text{face}} = v_g + (p\Omega p^{-1}) \times q \qquad (1)$$

$$v_{\text{nb}} = (v_{\text{face}} \cdot \hat{n}_{\text{b}}) \hat{n}_{\text{b}}$$

where $v_g$ is the velocity at the block gravity centre, $\Omega$ is the rotational velocity, $p$ is the quaternion (presenting block's 3D orientation) and $q$ is the relative position vector of the block face (envisaged sensor location) to its centre of gravity. The time evolution of velocity normal to block face computed for all three cases is presented in Figure 5 (all left). In addition, the peak velocity of the blocks and the time to reach this peak are plotted in Figure 5 (all right).

A general trend where the peak velocity localises at the impacted block is observed. Besides, the number of blocks experiencing rapid displacement from the impact beginning is comparatively higher when the projectile impacts the concave-wedge (Fig. 5c). Indeed, three blocks close to the impact location experience very similar curves. The time for the wall to come to rest is smaller in case of an impact on the convex-wedge with a duration of about 0.5 as compared to 0.8 seconds after an impact on the angled-wedge in particular (Fig. 5a vs Fig. 5b). Overall, the velocity of all the blocks appears to provide rich information, with high amplitude as well as significant differences from one case to the other.

The blocks with significant maximum velocities are located within a distance along the wall length ranging from 8 to 12 m approximately depending on the impact location. This observation is in line with the observation made in Figure 4. The shape of the maximum velocity pattern significantly differs from one impact case to the other. These observations suggest that the distribution of peak velocity along the wall could be specific to an impact case in a rather univocal manner.

### 3.1.3. Deflection

In addition to the displacement and velocity, the deflection of the wall blocks, i.e., their rotation around a horizontal axis, seems relevant to address as it relates to wall post-impact stability. In addition, deflection may be measured on-site with rather low-cost sensors. The post-impact deflection of all the blocks in the wall is presented in Figure 6. An impact on the concave-wedge results in block deflection over a wall length as large as 8 m. The maximum value of about 17° is observed at the upper row in the impact vicinity. Notably, the deflection results concerning the impact on the convex-wedge and angled-wedge are not presented as they induced nearly zero wall deflection.

Conclusively, the deflection doesn't provide sufficiently rich information to envisage relating it to the
<table><caption>Table 1: Parameters considered for generating close-to-reality impact conditions</caption><thead><tr><td>Parameter</td><td>Possible range</td><td>Distribution</td><td>Unit</td></tr></thead><tbody><tr><td>Translational velocity (v)</td><td>10 - 25</td><td>Gaussian</td><td>m/s</td></tr><tr><td>Rotational velocity (Ω)</td><td>0.0 - 5.6</td><td>Uniform</td><td>rot/sec</td></tr><tr><td>Impact position - along length (offset) (y)</td><td>0.0 - -3.53</td><td>Uniform</td><td>m</td></tr><tr><td>Impact position - along wall height (z)</td><td>0.55 - 2.10</td><td>Uniform</td><td>m</td></tr><tr><td>Impact inclination (α)</td><td>-60 - +60</td><td>Uniform</td><td>°</td></tr><tr><td>Impact deviation (β)</td><td>-45 - +45</td><td>Gaussian</td><td>°</td></tr></tbody></table>

190 study, the translation velocity results in a kinetic energy ranging from 130kJ to 800 kJ. The rotational velocity range is defined as per the work reported by Bourrier et al. [26] providing the rotational kinetic energy up to 240kJ for the projectile geometry used in this work.

193 The impact location refers to the projectile's centre of gravity at impact. The range for the impact locations along the wall longitudinal axis (Y axis) is defined considering that the structure consists of the repetition of a pattern. Thereby, the impact locations were restricted to a representative segment of the wall conformation. By definition, the response observed for this segment is extendable to other parts of the wall, except for the wall extremities. Indeed, impacts at a close distance from the wall extremities shall result in a different impact response, in particular depending on some design choice with influence on the mobility of the blocks at the wall extremity (such as retaining cables or abutments). The case of impacts close to the wall extremities is thus not considered in the present work. The impact location along the vertical axis (z) is defined considering the wall height and the block size, which is 1.1m. The minimum value for the impact location along the vertical axis is half the projectile dimension. The maximum value corresponds to the distance between the projectile centre of gravity and the wall crest equalling the projectile diameter.

205 The inclination and deviation angles are varied over the ranges reported by Toe et al. [24]. The inclination angle (α) represents the relation of the z-component of projectile impact path with the normal to the wall face [13]. Negative α values thus correspond to upward trajectories (i.e., away from the ground). The β angle accounts for the deviation in the horizontal plane of the incident rock block
<figure><img src="image_7.png" /><figcaption>Figure 5: Evolution of the velocity of the concrete block in the movement zone (all left) and the corresponding mapping of the maximum velocity for blocks in all four rows i.e., top, third, second and base (all right) for the illustrative impact cases at (a) convex wedge (IC-1), (b) angled wedge (IC-2) and concave wedge (IC-3) respectively. Here, 'Y' refers to the global position of the mid-bottom of the rear face of each block along the wall length.</figcaption></figure>

226 impact conditions, in particular, because some impact conditions result in null deflection angle values and hence are not deemed appropriate to explore further.
trajectory with respect to the normal to the wall longitudinal axis (and not with respect to the impacted block uphill face). The range considered for $\beta$ implies that the wall longitudinal axis is considered perpendicular to the line of the maximum gradient of the uphill slope.

In this study, the kinetic energy is considered sufficient for describing the projectile kinematics in view of addressing the structure response. This obviously constitutes a simplification as it is evidenced that, for a given projectile kinetic energy, the response of a rockfall protection structure may be significantly influenced by the projectile mass-to-velocity ratio [27, 28]. Besides, the upper limit of the translational kinetic energy is kept less than the wall nominal capacity as observed during the experiments, which is above 1 MJ. This is motivated by two observations: (1) The representation of the model to correctly model the structure displacement and damage is lower at impact kinetic close to the nominal value (e.g., due to loss of mechanical continuity in connectors [16]) and (2) The damage and displacement are severe for translational kinetic energy at impact higher than 800 kJ, imposing visit and repair/maintenance works without any doubts.

### 3.2.1. Statistical sample

A statistical set of 300 model computations, based on the Latin hypercube sampling (LHS) method [29] is retrieved to represent the model response under multiple impact conditions. The created statistical sample is presented in Figure 7. Here, the filled input space assures that possible combinations of the input variables are well considered in the finite set of ICPs.

It is worth highlighting that the ICPs were defined ignoring some correlations that may exist between the six parameters. For example, the velocity in the case of a downward projectile trajectory is higher on average than that for an upward trajectory. This results from the fact that a downward trajectory is associated with a rock block-free fall while an upward projectile trajectory immediately succeeds in a rebound in the ditch which induces energy dissipation. This means that some combinations of parameter values are not realistic (for example the lowest $\alpha$ value with the highest translational velocity, $v$).

Moreover, a probability distribution is assigned to each parameter to account for the influence of rock-fall trajectory, as suggested in Lambert et al. [14]. A Gaussian distribution is chosen for the translational velocity ($v$) and deviation angle ($\beta$) indicating that the extreme limits of both these parameters are less probable. The parameter $\beta$ distribution abides by the idea that the structure installation on site is such that the projectile is most likely to impact normally to its face. Similarly, the parameter $v$ distribution
<figure><img src="image_8.png" /><figcaption>Figure 6: Deflection of all blocks in four layers of the wall (i.e., top, third, second and base) with respect to the vertical at rest for the impacts at concave wedge (IC-3) under 520kJ energy impact. The results for the impact at the convex-wedge (IC-1) and angled-wedge (IC-2) are not presented because of their negligible magnitude.</figcaption></figure>

### 3.1.4. Plastic damage

The rockfall impact on the structure causes plastic damage to the blocks, dissipating a portion of the incident projectile kinetic energy. Notably, the NSCD model is created as an accumulation of rigid bodies where the contact law controls the interaction. Despite that, the quantitative estimation of the energy dissipation mechanisms is possible following the work reported by Acary [23] in the NSCD framework. Subsequently, the relation presented in Equation 2 is implemented in the present work to compute the energy dissipation due to plastic damage ($D_p$).

$$D_p = \sum_{k=0}^{N} \sum_{\alpha \in I} \frac{1}{2} (v_{N,k+1}^{\alpha} + v_{N,k}^{\alpha}) P_{N,k+1}^{\alpha} \quad (2)$$

Here, $v_N$ is the normal component of the velocity vector and $P_N$ is the normal impulse at a given time step 'k' (out of total N) at a contact point $\alpha \in I$. The details of velocity and impulse computations can be referred to at Gupta et al. [16].

The time evolution of $D_p$ for the three illustrative cases with different impact locations is presented in Figure 7. A significant contribution of plastic damage is observed at the early time instance when the projectile impacts the wall. With reference to the zoomed-in window, the damage evolves relatively fast for the impact at convex-wedge (IC-1) and reaches a constant magnitude of about 380kJ in about 0.35 seconds post-impact. The other two impact cases evolve for a relatively longer time (about 0.6 seconds)
<figure><img src="image_11.png" /><figcaption>Figure 7: Input sample of size 300 (based on Latin hypercube sampling (LHS) sampling method) comprising six impact condition parameters (ICPs) and their probability distribution in the respective ranges of variability.</figcaption></figure>

takes into account the most likely impacting projectile energy and keeps the low and high energy cases as rare events. The uniform distribution is assigned for all the remaining parameters as no prior information is available on their probability of occurrence.

### 3.2.2. Displacements under impact

The position of the top row blocks at the end of impact (i.e., at rest) for all 300 simulations is presented in Figure 8. The hollow circles correspond to the location of the extremity of the vertical connectors position. The blocks represented in blue correspond to the impact case resulting in the least distance from the wall to the safety line. Impacts in the convex-wedge resulted in a lower distance to the safety line which is consistent with observations made based on Figure 3. The movement zone for all simulated cases stays within approximately 7m on either side from the centre of the wall (i.e., $Y = 0$). The cloud of points, together with the mean positions and the variability (estimated as twice the standard deviation), reveals that the amplitude of the displacement in the X- and Y-axis directions are much larger at a distance from the wall centre (i.e. at the concave wedge) and that it is significantly dependent on the impact location along the Y-axis.
<figure><img src="image_9.png" /><figcaption>Figure 7: Time evolution of the energy dissipation due to plastic damage for the three illustrative cases, with impact at the convex-wedge (IC-1), the angled-wedge (IC-2) and concave-wedge (IC-3) and the zoom-in window between 300 and 400 kJ</figcaption></figure>

243 before reaching a constant value towards the end of the one-second simulation time. This observation is in line with the duration over which the wall moves, as illustrated in Figure 5. Conclusively, the plastic damage is taken as a quantity of interest for structural response assessment as it presents the variability in response with different impact conditions.

247 Notably, the remaining part of the energy transferred by the projectile to the wall is dissipated by friction at the contacts, as detailed in Lambert et al. [24] where it is also demonstrated that the computation scheme complies with the fundamental principle of energy conservation.

## 3.2. Response to a large set of impact conditions

Following the illustrative structural response description in the previous section, the numerical simulation framework is extended to address the response considering a large set of impact conditions. In a similar way as previously done for flexible barriers [13, 25, 26], the parameters describing the incident projectile trajectory, kinematics and location, hereafter referred to as ICPs (Impact Conditions Parameters) are varied over realistic ranges to account for nearly all possible distinct impact conditions which may occur during natural rockfall event. The ICPs and their magnitude ranges are given in Table 1.

The first two parameters (i.e., translational and rotational velocities) account for the incident kinetic
<figure><img src="image_12.png" /><figcaption>Figure 8: Displacement of the wall's top row blocks under 300 sampled impact conditions. The position of the blocks holes position after impact is presented for each simulation and the impact case resulting in the least distance to the safety line is highlighted for better visibility.</figcaption></figure>

### 3.2.3. Sensitivity analysis

The variability in structure response is further addressed by investigating the influence of each pa- rameter on the different quantities of interest (QoI). In accordance with Section 3.1, the considered QoIs are the minimum distance to the safety line ($u^{safety}$), the maximum displacement of the wall ($u_{max}$), the maximum concrete block velocity recorded ($v_{max}$) during impact and plastic damage ($D_p$).

The influence of each parameter on the QoIs is investigated through the Sobol sensitivity method also known as analysis of variance [30]. This method decomposes the variance of the output parameters as the sum of the contributions of the different input parameters including the possible interaction between input parameters. Each contribution is characterised by the ratios of the partial variance to the total variance, called Sobol sensitivity indices.

The accurate computation of Sobol indices demands a large number of model computations (of order $10^6$). This is highly impractical in the present study if the NSCD model is used directly. This limitation is circumvented with the help of the meta-modelling technique which enables to create a surrogate of the NSCD model, allowing direct computation of Sobol indices at a highly reduced cost [31, 32]. In this case, a meta-model (or surrogate) can be defined as a mathematical operator describing the response envelope of the wall in the 6D space corresponding to the six variables defining the impact conditions related to the rock projectile.

Here, a database of 300 model simulations is processed for each QoI to acquire the model output set and
<table><caption>Table 1: Projectile model impact condition parameters (ICPs), their magnitude range and the probabilistic distributions for close-to-reality wall impact computations. Here, the Uniform (U) distribution characterised by the minimum (a) and maximum (b) parameters, and Gaussian/Normal (N) distribution characterised by mean (μ) and standard deviation (σ) parameters are used.</caption><thead><tr><td>Parameter</td><td>Possible range</td><td>Unit</td><td>Distribution</td></tr></thead><tbody><tr><td>Translational velocity (v)</td><td>10 - 25</td><td>m/s</td><td>N(μ = 17.5, σ = 2.5)</td></tr><tr><td>Rotational velocity (Ω)</td><td>0.0 - 5.6</td><td>rot/sec</td><td>U(a = 0.0, b = 5.6)</td></tr><tr><td>Impact position - along length 'offset' (y)</td><td>0.0 - -3.53</td><td>m</td><td>U(a = 0.0, b = 3.53)</td></tr><tr><td>Impact position - along wall height (z)</td><td>0.55 - 2.10</td><td>m</td><td>U(a = 0.55, b = 2.10)</td></tr><tr><td>Impact inclination (α)</td><td>-60 - +60</td><td>°</td><td>U(a = -60, b = 60)</td></tr><tr><td>Impact deviation (β)</td><td>-45 - +45</td><td>°</td><td>N(μ = 0, σ = 15)</td></tr></tbody></table>

258

259

260

261

Given the projectile geometric and mechanical characteristics considered in this study, the translation velocity results in a kinetic energy ranging from 130kJ to 800 kJ. The rotational velocity range is defined as per the work reported by Bourrier et al. [27] providing the rotational kinetic energy up to 240kJ for the projectile geometry used in this work.

262

263

264

265

266

267

268

269

270

271

272

273

274

The impact location refers to the projectile's centre of gravity at impact. The range for the impact locations along the wall length (Y axis) is defined considering that the structure consists of the repetition of a pattern. Thereby, the impact locations were restricted to a representative segment of the wall conformation inclusive of all three wedges (i.e., convex, angled and concave), defined as 'y' offset or simply 'y' for future mentions. By definition, the response observed for this segment is extendable to other parts of the wall, except for the wall extremities. Indeed, impacts at a close distance from the wall extremities shall result in a different impact response, in particular depending on some design choice with influence on the mobility of the blocks at the wall extremity (such as retaining cables or abutments). The case of impacts close to the wall extremities is thus not considered in the present work. The impact location along the vertical axis (z) is defined considering the wall height and the block size, which is 1.1m. The minimum value for the impact location along the vertical axis is half the projectile dimension. The maximum value corresponds to the distance between the projectile centre of gravity and the wall crest equalling the projectile diameter.
270

271

272

273

274

275

276

277

its corresponding ICPs. Then, the polynomial chaos expansion (PCE) based meta-modelling technique (as per UQlab PCE module [33], see Appendix A for details) is used to formulate a generalised link between the input ICPs and the output QoI. For each meta-model creation, a mathematical relation (as presented in Equation A.2), is established between the 300 distinct sets of input parameters and the corresponding displacement output set. The mathematical relation here can be analogically referred to as obtaining a regression for a 2D database. The accuracy of the PCE-based meta-model is presented by the leave-one-out (LOO) error (Equation A.5), reported of order $10^{-2}$ for $u^{safety}$, $u_{max}$ and $D_p$ and of order $10^{-1}$ for $v_{max}$.

278

279

280

The UQlab sensitivity analysis module [34] is used for such computations. The first and total order Sobol indices for the QoI are presented in Figure 9. The Sobol indices are computed to evaluate the influence of each of the six ICPs on the four aforementioned QoI.

<figure><img src="image_13.png" /><figcaption>Figure 9: Sobol sensitivity analysis of the ICPs on different QoI towards the representation of the generalised impact assessment.</figcaption></figure>

281

282

283

284

285

The first-order Sobol indices reflect the main effect of each ICP and the total order indices reflect the main effect inclusive of the contribution from the interaction between different ICPs. When there is no interaction between variables, the first and total index magnitudes shall be the same. However, the contrary is observed here which reflects the significant interaction between ICPs in terms of their influence on the QoIs.

286

287

288

289

Moreover, the main effect of the rotational velocity (except for $D_p$), the vertical position of impact ($z$) and deviation ($\beta$) are observed to have the lowest sensitivity to all QoIs. On the contrary, all QoIs are reported to be highly sensitive to the other three input ICPs. The $y$ position of projectile impact is observed to be highly influential confirming the structural response variability from convex to angled to
The inclination and deviation angles are varied over the ranges reported by Toe et al. [25]. The inclination angle (α) represents the relation of the z-component of projectile impact path with the normal to the wall face [13]. Negative α values thus correspond to upward trajectories (i.e., away from the ground). The β angle accounts for the deviation in the horizontal plane of the incident rock block trajectory with respect to the normal to the wall longitudinal axis (and not with respect to the impacted block uphill face). The range considered for β implies that the wall longitudinal axis is considered perpendicular to the line of the maximum gradient of the uphill slope. Also, β determines the planer (x,y) components of Ω while the z-component is considered to have no magnitude such that the incoming projectile spin is always in line with its trajectory. Notably, the sketch in the Figure 11 visually illustrates the considered inclination and deviation angles.

In this study, the kinetic energy is considered sufficient for describing the projectile kinematics in view of addressing the structure response. This obviously constitutes a simplification as it is evidenced that, for a given projectile kinetic energy, the response of a rockfall protection structure may be significantly influenced by the projectile mass-to-velocity ratio [28, 29]. Besides, the upper limit of the translational kinetic energy is kept less than the wall nominal capacity as observed during the experiments, which is above 1 MJ. This is motivated by two observations: (1) The representation of the model to correctly model the structure displacement and damage is lower at impact kinetic close to the nominal value (e.g., due to loss of mechanical continuity in connectors [16]) and (2) The damage and displacement are severe for translational kinetic energy at impact higher than 800 kJ, imposing visit and repair/maintenance works without any doubts.

### 3.2.1. Statistical sample

A statistical set of 300 model computations, based on the Latin hypercube sampling (LHS) method [30] is retrieved to represent the model response under different impact conditions. The created statistical sample is presented in Figure 8, where the magnitudes of all parameters are in the range reported in Table 1. Here, the filled input space assures that possible combinations of the input variables are well considered in the finite set of ICPs. In general, 50-500 model runs are considered a 'reasonable budget' to create the metamodels of the NSCD model presented later on. The NSCD model computation time is about 30-40 minutes for one projectile impact simulation and hence the computation cost of 300 runs is deemed sufficient. Such sample size also enables the creation of all requisite metamodels to exploit
concave wedges, as reported in Figures 3 and 8.

The $v$ and $\alpha$ parameters collectively account for the projectile velocity component parallel to the ground (xy plane) and the component parallel to the wall face (z-axis) as presented in Equation 3.

$$v_z = v \sin \alpha$$

$$v_{xy} = v \cos \alpha \qquad (3)$$

$$v_x = v_{xy} \cos \beta, \quad v_y = v_{xy} \sin \beta$$

It is thus evident that these three input ICPs control all resulting QoI directly. Notably, the parameter $\beta$ determines the further division of projectile velocity ($v$) in x-y plane (i.e., $v_x$ and $v_y$). Since $\beta$ is observed to have a low influence on the QoI, the collective effect of projectile velocity in xy-plane is sufficient for further consideration.

### 3.3. Towards inverse analysis

The investigation of the structure response under close-to-reality conditions is of paramount importance, first, for evaluating its on-site efficiency and, second, for defining the best strategy for developing the inverse analysis method.

The wall response to impact is revealed extremely complex, with a significant dependence on parameters describing the impact conditions (ICPs). Differences in response were revealed by the concrete blocks' kinematics, in terms of displacement and velocity. The sensitivity to the impact point location along the wall's longitudinal axis results from the wall's zig-zag conformation and its discrete and articulated nature. In the case of an impact with a zero deviation and at a 520 kJ impact energy, the length of the wall experiencing significant displacement is up to 12m. This is confirmed by the distinct patterns of the peak concrete block velocity along the wall (e.g., see Figure 4). In addition, the wall response depends on the projectile orientation with respect to the wall face ($\alpha$ and $\beta$) which is attributed to the amount of energy that is transferred to the wall during impact (which will be further addressed in section 5.2).

The complexity of this response justified considering a specific approach for developing a relevant inverse analysis method. In addition, it is concluded that the inclination of the wall is not relevant for conducting an inverse analysis. On the contrary, the distribution of the concrete block's peak velocity seemed to be impact-conditions-specific and is thus considered a good input data candidate for conducting
304

various response mechanisms without optimising the sample size for each mechanism separately.

<figure><img src="image_10.png" /><figcaption>Figure 8: Input sample of size 300 (based on Latin hypercube sampling (LHS) sampling method) comprising six impact condition parameters (ICPs) and their probability distribution in the respective ranges of variability.</figcaption></figure>

305

306

307

308

309

310

311

312

Moreover, a probability distribution is assigned to each parameter to account for the influence of rockfall trajectory, inspired by the work reported in Lambert et al. [14]. A Gaussian distribution is chosen for the translational velocity (v) and deviation angle (β) indicating that the extreme limits of both these parameters are less probable. The parameter β distribution abides by the idea that the structure installation on-site is such that the projectile is most likely to impact normally to its face. Similarly, the parameter v distribution takes into account the most likely impacting projectile energy and keeps the low and high-energy cases as rare events. The uniform distribution is assigned for all the remaining parameters as no prior information is available on their probability of occurrence.

313

314

315

316

317

318

It is worth highlighting that the ICPs were defined ignoring some correlations that may exist between the six parameters. For example, the velocity in the case of a downward projectile trajectory is higher on average than that for an upward trajectory. This results from the fact that a downward trajectory is associated with a rock block-free fall while an upward projectile trajectory immediately succeeds the rebound from the ditch which induces energy dissipation. This means that some combinations of parameter values are not realistic (for instance, the lowest α value with the highest translational velocity, v).
inverse analysis based on real-time measurements made on-site. These conclusions suggest that the priority in terms of on-site structure instrumentation should be placed on accelerometers, from which the concrete block velocity can be derived. By contrast, it is less relevant to install inclinometers for the purpose of conducting inverse analysis.

In this work, it is proposed to conduct inverse analysis based on real-time measurements, during impact, for the purpose of warning and decision-making upon impact. This requires installing sensors in different locations on the structure, connected to data acquisition and transfer equipment. Then, it is also proposed to conduct inverse analysis based on data collected on the structure after impact, such as a cloud of points revealing the structure envelope obtained from scanning tools. These data sets may be used for conducting post-impact inverse analysis in view, for example, of quantifying the consequence of the impact on the structure or retrieving information related to the event.

In both cases, the inverse analysis basically consists of exploiting the data collected on-site to statis-
tically establish a link with the various impact conditions. The inverse analysis relies on meta-models
created based on a large set of NSCD model simulations of the structure response under close-to-reality
conditions. The Bayesian interface statistical learning method is used where the data collected on-site
shall be fed as evidence into the meta-models based forward model to find out the likelihood of obtaining
the same result (a single value or a pattern) from a particular set of input parameters. A brief descrip-
tion of the Bayesian interface is presented in Appendix B and its implementation in the present work is
detailed in section 4.2.

4. Inverse analysis based on real-time measurements

4.1. *General considerations*

The instrumented protective structures e.g., installed sensors, most often aim at warning the structure
owner or the protected infrastructure manager of any event. Real-time measurements could also be used
to improve decision-making over the next few moments after impact. Such a process could be based on
images from cameras, but sensors have the advantage of providing reliable information independently of
the light and scene conditions.

Two decisions could be taken remotely immediately after an impact on a structure protecting a traffic line: traffic interruption and the need for a structure inspection. Both are related to the consequences of
319

320

321

322

323

A specific site-investigation-based ICP correlation can be established for a better representation of the input sample. However, this is out of the scope of the present work and hence we proceed further with the uncorrelated ICPs. Moreover, the Gaussian distribution for *v* and *β* parameters results in the representation of each ICP set as not equiprobable. Subsequently, it enables having less number of simulated cases corresponding to rare events.

324

### 3.2.2. Displacements under impact

325

326

327

328

329

330

331

332

333

334

The position of the top row blocks at the end of impact (i.e., at rest) for all 300 simulations is presented in Figure 9. The hollow circles correspond to the location of the extremity of the vertical connectors position. The blocks represented in blue correspond to the impact case resulting in the least distance from the wall to the safety line. Impacts in the convex-wedge resulted in a lower distance to the safety line which is consistent with observations made based on Figure 4. The movement zone for all simulated cases stays within approximately 7m on either side from the centre of the wall (i.e., $Y = 0$). The cloud of points, together with the mean positions and the variability (estimated as twice the standard deviation), reveals that the amplitude of the displacement in the X- and Y-axis directions are much larger at a distance from the wall centre (i.e. at the concave wedge) and that it is significantly dependent on the impact location along the Y-axis.

<figure><img src="image_11.png" /><figcaption>Figure 9: Top view of the displaced component of the whole length of the wall as per sampled 300 sets of projectile impact condition parameters. The position of the block holes position after impact is presented for each simulation and the impact case resulting in the least distance to the safety line is highlighted for better visibility.</figcaption></figure>
343

344

345

346

347

the impact on the structure, in terms of displacement and damage. An excess in one of these latter may be critical to the safety of the traffic line or imply rapid repair or maintenance work. In this objective, the inverse analysis is intended to rapidly and remotely provide an indication on the magnitude of the consequences of the impact, for example using consequence classes, rather than to provide a precise value of the displacement or damage.

348

349

350

351

352

353

354

355

356

The inverse analysis is based on measurements retrieved from the on-site structure. The devices used for these measurements should consider cost constraints to be of interest in an operational context. In particular, the number and type of sensors should be optimised. The previous section suggested that a reliable inverse analysis could be conducted from the block velocity, which can be derived from the on-site acceleration measurements recorded during the impact. Based on this, the proposed inverse analysis is developed considering acceleration measurements from uniaxial sensors, because of their much lower cost compared to triaxial ones. Also, it is developed to constrain the number of required data, meaning the number of measuring points on the structure (i.e., sensors). The number of measuring points is kept at a minimum while allowing a sufficiently good prediction accuracy of the inverse analysis.

357

358

359

360

361

362

363

364

365

366

367

In this perspective, the sensor deployment scheme is defined for a representative wall length, in such a manner it can be replicated along the full wall length. Considering the minimum length of the movement zone observed in Figure 4, the choice is made to have at least 3 sensors per wall segment. More precisely, sensors shall be located in the centre of the vertical face of concrete blocks experiencing significant displacements whatever the ICPs. This number of sensors per representative length is presumed to be the minimum number for allowing a reliable inverse analysis. It corresponds to 6 concrete blocks equipped with sensors over the maximum movement zone length, as shown in Figure 10. Here, the illustrative process of acceleration data acquisition followed by its integration to obtain velocity evolution with time for each sensor is also presented. The maximum velocity magnitude is identified and consequently, the pattern of maximum velocity for each sensor processed data is retrieved which is deemed unique for each impact condition.

368

369

370

371

The inverse analyses presented in the following make use of the maximum concrete block velocity along the direction normal to the block face at these 6 locations. The set of maximum block velocities recorded along the wall at these locations is in the following referred to as the maximum velocity pattern (MVP).
### 3.2.3. Sensitivity analysis

The variability in structure response is further addressed by investigating the influence of each pa- rameter on the different quantities of interest (QoI). In accordance with Section 3.1, the considered QoIs are the minimum distance to the safety line ($u^{safety}$), the maximum displacement of the wall ($u_{max}$), the maximum concrete block velocity recorded ($v_{max}$) during impact and plastic damage ($D_p$).

The influence of each parameter on the QoIs is investigated through the Sobol sensitivity method also known as analysis of variance [31]. This method decomposes the variance of the output parameters as the sum of the contributions of the different input parameters including the possible interaction between input parameters. Each contribution is characterised by the ratios of the partial variance to the total variance, called Sobol sensitivity indices.

The accurate computation of Sobol indices demands a large number of model computations (of order $10^6$). This is highly impractical in the present study if the NSCD model is used directly. This limitation is circumvented with the help of the meta-modelling technique which enables to create a surrogate of the NSCD model, allowing direct computation of Sobol indices at zero cost [32, 33]. In this case, a meta-model (or surrogate) can be defined as a mathematical operator describing the response envelope of the wall in the 6D space corresponding to the six variables defining the impact conditions related to the rock projectile.

Here, a database of 300 model simulations is processed for each QoI to acquire the model output set and its corresponding ICPs. Then, the polynomial chaos expansion (PCE) based meta-modelling technique (as per UQlab PCE module [34], see Appendix A for details) is used to formulate a generalised link between the input ICPs and each output QoI. A total of four metamodels are created corresponding to four QoI, named Meta-$u^{safety}$, Meta-$u_{max}$, Meta-$v_{max}$ and Meta-$D_p$. For each meta-model creation, a mathematical relation (as presented in Equation A.2), is established between the 300 distinct sets of input parameters and the corresponding displacement output set. The mathematical relation here can be analogically referred to as obtaining a regression for a 2D database. The accuracy of the PCE-based meta-model for each QoI is presented by the leave-one-out (LOO) error (Equation A.5), reported of order $10^{-2}$ for Meta-$u^{safety}$, Meta-$u_{max}$ and Meta-$D_p$ and of order $10^{-1}$ for Meta-$v_{max}$. Notably, these PCE-based meta-models are also used for the inverse analysis work, detailed in the section 4.

The UQlab sensitivity analysis module [35] is used for the computation of the Sobol indices. The first
<figure><img src="image_14.png" /><figcaption>Figure 10: The schematics of the wall model where some blocks are equipped with sensors for data acquisition</figcaption></figure>

## 4.2. Inverse analysis method

The inverse analysis relied on the confrontation between the MVP obtained from on-site measurements to that obtained from a large number of simulation-based results. The developed method for conducting the inverse analysis comprised a combination of meta-models and Bayesian interface statistical analysis numerical methods. The flow diagram is presented in Figure 11 and detailed as follows.

<figure><img src="image_15.png" /><figcaption>Figure 11: Computational workflow of the inverse analysis using a combination of meta-models, event data and Bayesian inversion methods</figcaption></figure>

The rich database of NSCD model simulations for 300 different ICPs (Figure 7) is processed (as per Equation 1) to obtain the maximum velocity normal to the block face for all blocks where a sensor is planned to be installed. This provides a database of MVP for 300 model computations. This database is
and total order Sobol indices for the QoI are presented in Figure 10, computed to evaluate the influence of each of the six ICPs on the four aforementioned QoI.

<figure><img src="image_12.png" /><figcaption>Figure 10: Sobol sensitivity analysis of the ICPs on different QoI towards the representation of the generalised impact assessment.</figcaption></figure>

The first-order Sobol indices reflect the main effect of each ICP and the total order indices reflect the main effect plus the contribution from the interaction between different ICPs. Thereby, when there is no interaction between variables, the first and total index magnitudes shall be the same. However, the contrary is observed here which reflects the significant interaction between ICPs in terms of their influence on the QoIs. The main and the total effects of each ICP are briefly discussed as follows.

The main effect of the rotational velocity (except for $D_p$), the vertical position of impact (z) and deviation (β) are observed to have the lowest sensitivity to all QoIs. On the contrary, all QoIs are reported to be highly sensitive to the other three input ICPs. The y position of projectile impact is observed to be highly influential confirming the structural response variability from convex to angled to concave wedges, as reported in Figures 4 and 9. The v and α parameters collectively account for the projectile velocity component parallel to the ground (xy plane) and the component parallel to the wall face (z-axis) as presented in Equation 3.

$$v_z = v \sin \alpha$$

$$v_{xy} = v \cos \alpha \qquad (3)$$

$$v_x = v_{xy} \cos \beta, \quad v_y = v_{xy} \sin \beta$$

It is thus evident that these three input ICPs influence all resulting QoI directly. Notably, the parameter β determines the further division of projectile velocity (v) in x-y plane (i.e., $v_x$ and $v_y$). Since
processed to create a PCE-based meta-model associating each ICP value set to a distinct MVP.

The accuracy of these meta-models is estimated via leave-one-out-error which is estimated of order $10^{-2}$ to $10^{-1}$. The accuracy of the meta-model representing the MVP, is validated by using the same simulation set complemented with an independent set of 300 simulations, for a total of 600 simulations, and comparing the predicted vs true value from the NSCD model at all six sensor locations, collectively presented in Figure 12.

<figure><img src="image_16.png" /><figcaption>Figure 12: Validation of the PCE-based meta-model of the maximum velocity by comparison of the reported magnitude from 600 NSCD model simulations and the corresponding estimation by the meta-model.</figcaption></figure>

The prediction accuracy of the created meta-model is quantified via the predictability coefficient ($Q^2$), estimated considering the results from the 600 simulations and 6 points, for a total of 3600 observations, as:

$$Q^2 = 1 - \frac{\sum_{i=1}^{N} (y_i^{\text{pred}} - y_i^{\text{true}})^2}{\sum_{i=1}^{N} (y_i^{\text{true}} - y_{\text{mean}})^2} \quad (4)$$

where $y_i^{\text{pred}}$ is the meta-model predicted value and $y_i^{\text{true}}$ is the value processed from the NSCD model for block (i) and $y_{\text{mean}}$ is the arithmetic mean of all 'true' predictions. The interest of using the meta-models is to generalise the maximum concrete block velocity data obtained from 300 (i.e., limited number of) simulations, to any possible ICP combination (i.e., a large number), which is a prerequisite for the next step.

The Bayesian inversion process is deployed where the processed PCE-based meta-model is used as
β is observed to have a low influence on the QoI, the collective effect of projectile velocity in xy-plane is sufficient for further consideration. The negligible influence of ‘z’ infers the stability of the wall structure throughout its height. Lastly, the negligible influence of ‘Ω’ on block displacement and velocity indicates that the projectile spin induces a reaction from the block predominately in the vertical direction. A relatively significant influence on $D_p$ further emphasises this interaction inferring that Ω contributes to the block damage with minimal induced displacement.

The total effect of each ICP is observed to reflect the interaction between them as their total Sobol index is observed higher than the corresponding first-order index. It is an indication of the non-linearity in the model response. The $v$ and $\alpha$ parameters collectively account for the energy transferred to the wall, which in turn influence all four evaluated QoIs. Thereby, an interaction between them is justified. Similarly, the parameter $y$, despite being independent of the velocity parameter, accounts for the effect of the wall pattern (i.e., convex, angled and concave wedges) on the wall response. Therefore, abiding with the illustrative observations in Figure 4, the same energy impact at different locations results in a different response, i.e., sensitive to all four QoIs. The effect on $D_p$ is observed least in comparison with the other three QoI which indicates that the damage to the block is predominantly influenced by the impact energy with a minor influence of the impact location. The reported $D_p$ from the three illustrative impact cases (see Figure 7) supports this observation. Lastly, similar to the main effect, the contribution of $\Omega$ and $z$ and $\beta$ parameters in the total effect is observed relatively less.

### 3.3. Towards inverse analysis

The investigation of the structure response under close-to-reality conditions is of paramount importance, first, for evaluating its on-site efficiency and, second, for defining the best strategy for developing the inverse analysis method.

The wall response to impact is revealed extremely complex, with a significant dependence on parameters describing the impact conditions (ICPs). Differences in response were revealed by the concrete blocks’ kinematics, in terms of displacement and velocity. The sensitivity to the impact point location along the wall’s longitudinal axis results from the wall’s zig-zag conformation and its discrete and articulated nature. For instance, in the case of IC-1 (i.e., an impact with a zero deviation and at a 520 kJ impact energy), the length of the wall experiencing significant displacement is about 12m (from -6 to +6), and similarly, it is about 8m for IC-2 and 9m for IC3 (see Figure 4). This is confirmed by the
the forward model (or 'prior', see Equation B.2). This approach is extensively used by many researchers (e.g., [35, 36, 16]) to accelerate Bayesian computations. The Bayesian interface is inspired by Bayes' theorem [37] - a representation of the changing beliefs - simply demonstrating that the probability of a 'hypothesis' being correct becomes more reliable with supporting 'evidence'. In our work, the 'hypothesis' states that the 300 distinct sets of input parameters present a set of 300 distinct MVPs. The idea of reliability increase in our hypothesis is that, if we have the on-site recorded MVP as 'evidence', then there exists at least one set of ICPs (i.e., 'posterior', see Equation B.3) such that the recorded evidence is reproduced. The data obtained from on-site measurements shall provide the MVP which serves as evidence for inverse analysis. The Bayesian process estimates the likelihood of this MVP amongst all possible MVPs recorded from the meta-models and predicts the corresponding set of ICPs.

In the UQlab Bayesian interface framework [38], the uncertainty in the model prediction is assigned via added Gaussian discrepancy (see Equation B.5) when correlated with the recorded on-site data. In the present work, as the inverse analysis shall rely on a single MVP recorded on-site, deemed less precise, the discrepancy with a known residual variance of order $10^{-3}$ is manually assigned.

The output from the Bayesian interface-based inverse analysis is a set of ICP values retrieved as point estimates 'maximum a posteriori' and 'mean' (see Equation B.13). This predicted set of ICP values is then used to estimate the QoI (in this case damage to the wall or distance to the safety line). As a conventional approach, running a new NSCD model simulation with the predicted ICPs shall provide detailed information on the wall response from which the QoI can be extracted as a subset. However, this operation shall require about 30 minutes of computation time and added time for post-processing. This added delay time between on-site observation and warning system activation is unfavourable.

A procedure by which the QoI is instantly computed is developed to reduce the reaction time. The computational workflow components for this process is highlighted in the dashed box of Figure 11. The database of NSCD model simulations results is processed to evaluate the QoI as outputs. The resulting database of QoI is then used to create the PCE-based meta-model. It serves as a surrogate of the NSCD model where the input ICPs are processed to compute the corresponding QoI output. The leave-one-out (LOO) error of order $10^{-2}$ is reported for the two PCE-based meta-models created for the two QoI (i.e., distance to safety and wall damage) mentioned above. The predicted ICPs are then fed into this surrogate of NSCD model as input which is then evaluated to eventually predict the QoI and the whole
distinct patterns of the peak concrete block velocity along the wall (see Figure 5). In addition, the wall response depends on the projectile orientation with respect to the wall face (α and β) which is attributed to the amount of energy that is transferred to the wall during impact (which will be further addressed in see section 5.2).

The complexity of this response justified considering a specific approach for developing a relevant inverse analysis method. In addition, it is concluded that the deflection of the wall is not relevant for conducting an inverse analysis. On the contrary, the distribution of the concrete block's peak velocity seemed to be impact-conditions-specific and is thus considered a good input data candidate for conducting inverse analysis based on real-time measurements made on-site. These conclusions suggest that the priority in terms of on-site structure instrumentation should be placed on accelerometers, from which the concrete block velocity can be derived. By contrast, it is less relevant to install inclinometers for the purpose of conducting inverse analysis.

In this work, it is proposed to conduct inverse analysis based on real-time measurements, during impact, for the purpose of warning and decision-making upon impact. This requires installing sensors in different locations on the structure, connected to data acquisition and transfer equipment. Then, it is also proposed to conduct inverse analysis based on data collected on the structure after impact, such as a cloud of points revealing the structure envelope obtained from scanning tools. These data sets may be used for conducting post-impact inverse analysis in view, for example, of quantifying the consequence of the impact on the structure or retrieving information related to the event.

In both cases, the inverse analysis consists of exploiting the data collected on-site to statistically establish a link with the various impact conditions. The inverse analysis relies on meta-models created based on a large set of NSCD model simulations of the structure response under close-to-reality conditions. The Bayesian inference statistical learning method is used where the data collected on-site shall be fed as evidence into the meta-models based forward model to find out the likelihood of obtaining the same result (a single value or a pattern) from a particular set of input parameters. A brief description of the Bayesian inference is presented in Appendix B and its implementation in the present work is detailed in section 4.2.
424 computational process terminates.

425 Notably, in the absence of an exhaustive database concerning real structures, the essential ‘evidence’
426 database (i.e., recorded MVP) to test the inverse analysis process is not available. Therefore, the concrete
427 block velocity patterns from the same set of 300 model simulations and a new independent set of 300
428 simulations with different impact conditions are used as pseudo-evidences to test the process.

429 Each of these 600 known block velocity patterns is individually fed as input into the inverse analysis
430 using the multiple model output feature of Bayesian analysis in UQlab (demonstrated in Equation B.8),
431 providing the predicted ICP values. The evaluation of the reliability of the inverse analysis is based
432 on the comparison between the QoI reported via the ICP used for running the simulation, and the one
433 estimated from evaluating the ICP predicted by the inverse analysis onto the meta-model created for the
434 same QoI.

435 **4.3. Inverse analysis for warning purpose**

436 First, the inverse analysis addresses the traffic interruption issue for which the post-impact minimum
437 distance of the wall to the safety line constitutes the QoI. The higher the distance, the higher the safety.
438 The predictive capability of the inverse analysis workflow presented in Figure 11 is tested for the maximum
439 concrete block velocity database obtained from the 600 NSCD model simulations. The MVPs from the
440 three impact cases mentioned in Section 3.1 are also used to predict the minimum distance to the safety
441 line.

442 The comparison between the predictions from the inverse analysis + meta-model and the simulation
443 results is presented in Figure 13. The relative proximity of the predictions to the ideal diagonal reflects
444 the potential of the presented methodology. Points above the diagonal concern cases for which the inverse
445 analysis overestimates the distance to the safety line. On the contrary, the inverse analysis underestimates
446 the distance to the safety lines for points below the diagonal, which is on the safe side.

447 The root mean square (RMS) error of the deviation from the expected response for the 600 simulations
448 is 16 cm. About 58% of this RMS error is due to the surrogate of the NSCD model (i.e., dashed box
449 component of the workflow). In fact, in some sectors of the 6D space of the six ICPs, this surrogate
450 locally fails in precisely capturing the wall response, in particular due to the presence of non-linearities.
451 As an alternative, the reliability of the inverse analysis could be improved by running an NSCD model
452 simulation with the predicted ICPs in place of using the surrogate of the NSCD model, at the expense
4. Inverse analysis based on real-time measurements

4.1. General considerations

The instrumented protective structures e.g., installed sensors, most often aim at warning the structure
owner or the protected infrastructure manager of any event. Real-time measurements could also be used
to improve decision-making over the next few moments after impact. Such a process could be based on
images from cameras, but sensors have the advantage of providing reliable information independently of
the light and scene conditions.

Two decisions could be taken remotely immediately after an impact on a structure protecting a traffic line: traffic interruption and the need for structure inspection. Both are related to the consequences of the impact on the structure, in terms of displacement and damage. An excess magnitude may be critical to the safety of the traffic line or imply rapid repair or maintenance work. In this objective, the inverse analysis is intended to rapidly and remotely indicate the magnitude of the consequences of the impact, for example using consequence classes, rather than to provide a precise value of the displacement or damage.

The inverse analysis is based on measurements that shall be retrieved from the on-site structure
exposed to the natural rockfall impact event. The devices used for these measurements should con-
sider cost constraints to be of interest in an operational context. In particular, the number and type
of sensors should be optimised. The previous section suggested that a reliable inverse analysis could
be conducted from the block velocity, which can be derived from the on-site acceleration measurements
recorded during the impact. Based on this, the proposed inverse analysis is developed considering accel-
eration measurements from uniaxial sensors, because of their much lower cost compared to triaxial ones.
Also, it is developed to constrain the number of required data, meaning the number of measuring points
on the structure (i.e., sensors). The number of measuring points is kept at a minimum while allowing a
sufficiently good prediction accuracy of the inverse analysis.

In this perspective, the sensor deployment scheme is defined for a representative wall length, in such a
manner it can be replicated along the full wall length. Considering the minimum length of the movement
zone observed in Figure 4, the choice is made to have at least 3 sensors per wall segment (i.e., component
of the wall length comprising of one convex, angled and concave wedge each). More precisely, sensors
shall be located in the centre of the vertical face of concrete blocks experiencing significant displacements
whatever the ICPs. This number of sensors per representative length is presumed to be the minimum
of a longer reaction time upon event.

<figure><img src="image_17.png" /><figcaption>Figure 13: Comparison of the minimum distance from the safety line, estimated from the NSCD model and the corresponding prediction from the inverse analysis tool. Here, the predictions are divided into three zones reflecting high (H), medium (M) and Low (L) levels of severity inclusive of the corresponding zonal false positive (+) and false negative (x) predictions. Results are based on a set of 600 impact conditions complemented with the impact cases presented in section 3.1 (IC-1 to IC-3).</figcaption></figure>

For practical purposes, it is proposed to classify the impact response into three classes depending on the distance to the safety line: [0-0.75 m], [0.75-1.25 m] and [1.25-1.5m]. These classes correspond to decreasing criticality concerning transport corridor safety, respectively high (H), medium (M) and low (L) criticality. NSCD model simulation results indicate that 15%, 59% and 26% of the impact cases are classified in class H, M and L respectively. Predictions based on the inverse analysis + metamodel give values of 12%, 63% and 25% for classes H, M and L respectively. Besides, the false negative predictions are estimated to be about 24% and 10% among all in the H and M classes respectively. The remaining predictions i.e., good and false positives are in favour of the decision-making by the user, and hence are deemed acceptable.

In the end, the predictions are considered sufficiently good to be used in an operational context for estimating the consequences of the rockfall event in terms of structure displacement. The undeniable interest in this approach is that it automatically and remotely provides stakeholders with information of great value for decision-making upon impact. No quantitative evaluation would be derived rapidly from data measured on-site in the absence of such an approach.
number for allowing a reliable inverse analysis. It corresponds to 6 concrete blocks equipped with sensors over the maximum movement zone length, as shown in Figure 11. Here, the illustrative process of acceleration data acquisition followed by its integration to obtain velocity evolution with time for each sensor is also presented. The maximum velocity magnitude is identified and consequently, the pattern of maximum velocity for each sensor processed data is retrieved which is deemed unique for each impact condition.

<figure><img src="image_13.png" /><figcaption>Figure 11: The schematics of the wall model indicating the blocks equipped with sensors for data acquisition. The sketch on the top right presents the notation of the impacting projectile inclination (α) and deviation (β)</figcaption></figure>

The inverse analyses presented in the following make use of the maximum concrete block velocity along the direction normal to the block face at these 6 locations to represent the impact response corresponding to all possible sets of ICPs. The set of maximum block velocities recorded along the wall at these locations is in the following referred to as the maximum velocity pattern (MVP).

## 4.2. Inverse analysis method

The inverse analysis relied on the confrontation between the MVP that shall be obtained from an on-site measurements to that obtained from a large number of simulation-based results. The developed method for conducting the inverse analysis comprised a combination of meta-models and Bayesian inference statistical learning methods. The flow diagram is presented in Figure 12 and detailed as follows.

The rich database of NSCD model simulations for 300 different ICPs (Figure 8) is processed (as per Equation 1) to obtain the maximum velocity normal to the block face for all blocks where a sensor is planned to be installed. This provides a database of MVP for 300 model computations. This database
4.4. Inverse analysis for structure damage quantification purpose

Second, the inverse analysis is used to remotely estimate damage to the wall. It is defined as the cumulative energy dissipated by the plasticization of all the wall concrete blocks, $D_p$, after the wall is at rest. Damage to all concrete blocks of the wall is computed from NSCD model simulations results as per section 3.1.4. The approach considered for this inverse analysis is the same as that for the previous warning criterion. Figure 14 provides a comparison between $D_p$ values obtained from NSCD model simulations and the predicted ones. The RMS error of the deviation from the expected response is 82 kJ, which may be considered acceptable when considering all uncertainties associated with the quantification of energy dissipation by plasticization due to real rockfall impacts. Similar to the previous prediction of distance, about 35% of the contribution to the RMS error is attributed to the surrogate of the NSCD model. Besides, a similar prediction for the three illustrative cases is also presented which reflects a good prediction except for the impact at the angled-wedge (i.e., IC-2).

<figure><img src="image_18.png" /><figcaption>Figure 14: Comparison of the predicted damage to the wall, $D_p$, estimated from the NSCD model and the corresponding prediction from the inverse analysis tool. Here, the predictions are divided into three zones reflecting Low (L), medium (M) and high (H) levels of damage inclusive of the corresponding zonal false positive (+) and false negative (×) predictions. Results are based on a set of 600 impact conditions and impact cases presented in section 3.1 (IC-1 to IC-3).</figcaption></figure>

The cases in the right upper corner constitute the most critical impact cases. Points above the diagonal concern cases for which the inverse analysis overestimates the damage, which is on the safe side. For practical purposes, it is proposed to classify the damage into three classes depending on the value of $D_p$: [0-100 kJ], [100-300 kJ] and [300-800 kJ] for low (L), medium (M) or high (H) damage. These limits were
<figure><img src="image_14.png" /><figcaption>Figure 12: Computational workflow of the inverse analysis using a combination of meta-models, event data and Bayesian inversion methods</figcaption></figure>

is processed to create a PCE-based meta-model, named Meta-MVP, associating each set of ICPs to a distinct MVP.

The accuracy of this meta-model is estimated via leave-one-out-error which is estimated of order $10^{-2}$ to $10^{-1}$ for 6 locations collectively making a MVP. Besides, the predictive capability of the Meta-MVP is presented by comparing the predicted vs true value from the NSCD model at all six sensor locations. For this, the same 300 simulations set complemented with an independent set of 300 simulations, where ICP sets generated using Halton sampling method [30] within the bounds reported in Table 1 are used. The collective outcome from a total of 600 simulations is presented in Figure 13.

The prediction accuracy of the Meta-MVP is quantified via the predictability coefficient ($Q^2$), estimated considering the results from the 600 simulations and 6 points, for a total of 3600 observations, as:

$$Q^2 = 1 - \frac{\sum_{i=1}^{N} (y_i^{\text{pred}} - y_i^{\text{true}})^2}{\sum_{i=1}^{N} (y_i^{\text{true}} - y_{\text{mean}})^2} \quad (4)$$

where $y_i^{\text{pred}}$ is the meta-model predicted value and $y_i^{\text{true}}$ is the value processed from the NSCD model for block (i) and $y_{\text{mean}}$ is the arithmetic mean of all 'true' predictions. The interest of using the meta-models is to generalise the maximum concrete block velocity data obtained from 300 (i.e., limited number of) simulations, to any possible ICP combination (i.e., a large number), which is a prerequisite for the
arbitrarily defined from the distribution of $D_p$ over all cases.

The class limits refer to the level of damage to the wall corresponding to the fracture energy of one concrete block when a fracture propagates along a vertical plane in the middle of the concrete block. A classification based on the energy dissipation associated with the fracture of the considered concrete blocks would be more appropriate. Nevertheless, this would require specific and complex tests which are beyond the scope of the present study. In the absence of specific dynamic test results for these large and highly reinforced concrete blocks, results obtained by Guo et al. [39] for a reinforced beam exposed to impact were extrapolated accounting for the difference in block section and resulted in a value of 100 kJ approx. This value is a rough estimate of the first limit for the considered concrete blocks.

According to NSCD model simulations results presented in Figure 14, 8%, 72% and 20% of the impact cases result in damage classified as H, M and L respectively. Values of 6%, 76% and 18% for classes H, M and L are obtained from the predictions respectively. Besides, the false negative predictions are estimated to be about 38% and 5% among all in the H and M classes respectively. The remaining predictions i.e., good and false positives are in favour of the decision-making by the user, and hence are acceptable.

The differences between simulation results and predictions are associated with points at a large distance from the diagonal (Fig. 14). Further method developments could improve the accuracy of these predictions. Nevertheless, the predictions following the described inverse analysis are considered sufficiently good to be used in an operational context to aid in remote decision-making, the RMS error is 82 kJ.

### 4.5. Comments concerning the use of inverse analysis based on real-time measurements

The results presented in section 4.3 and 4.4 demonstrate the feasibility of conducting an inverse analysis of data collected on-site during impact, focusing on two specific purposes.

The proposed method, making use of numerical modelling, meta-modelling and statistical approaches, revealed rather efficient in providing useful information which could not be derived from the direct observation of the data received from the site. The duration of the whole process for treating the input data typically lasts about five minutes with the developed script making it compatible with rapid reaction upon event. The analysis scripts are written in Matlab as per UQLab framework [40] allowing running this complex process in two operations (one for inverse analysis and the other for output of interest prediction). The accuracy of the predictions is considered sufficiently good for remotely classifying the
<figure><img src="image_15.png" /><figcaption>Figure 13: Validation of the PCE-based meta-model of the maximum velocity pattern (Meta-MVP) by comparison of the reported magnitude from 600 NSCD model simulations and the corresponding estimation by the meta-model.</figcaption></figure>

next step.

The Bayesian inversion process is deployed where the processed PCE-based meta-model is used as the forward model (or ‘*prior*’, see Equation B.2). This approach is extensively used by many researchers (e.g., [36, 37, 16]) to accelerate Bayesian computations. The Bayesian inference is inspired by Bayes’ theorem [38] - a representation of the changing beliefs - simply demonstrating that the probability of a ‘hypothesis’ being correct becomes more reliable with supporting ‘evidence’. In our work, the ‘hypothesis’ states that the 300 distinct sets of input parameters present a set of 300 distinct MVPs. The idea of reliability increase in our hypothesis is that, if we have the on-site recorded MVP as ‘evidence’, then there exists at least one set of ICPs (i.e., ‘*posterior*’, see Equation B.3) such that the recorded evidence is reproduced. The data that shall be obtained from on-site measurements for a real projectile impact on the wall event shall provide the MVP which serves as evidence for inverse analysis. The Bayesian process estimates the likelihood of this MVP amongst all possible MVPs recorded from the meta-models and predicts the corresponding set of ICPs.

In the UQlab Bayesian inference framework [39], the uncertainty in the model prediction is assigned via added Gaussian discrepancy (see Equation B.5) when correlated with the data that shall be recorded from the real on-site event of the projectile impact on the wall. In the present work, as the inverse analysis relies on a single MVP that shall be recorded on-site, deemed less precise, the discrepancy with a known residual variance of order $10^{-3}$ is manually assigned.
events and aiding in decision-making. If required for another purpose, this accuracy could be improved further, for example by also accounting for other measure types, such as inclination, or by improving the used methodologies for creating the meta-models.

In real situations, the probability of occurrence of an impact in given conditions is extremely variable. By contrast, and as mentioned in section 3.2, the ICP sets considered in this work were defined by not considering any possible correlation between the various parameters and, more generally, considering that the ICP sets were equiprobable. The statistical assessment of the inverse analysis predictive capacities is related to the considered ranges and distributions of the six ICPs. It should thus be considered indicative at this stage. It should be improved by accounting for the probability of each scenario considered in this evaluation, for example considering site-specific data describing the rock projectile trajectory.

It is worth noting that, as the second step of the inverse analysis provides the set of ICP values resulting in an MVP similar to the observed one, then a preliminary estimate of the projectile velocity at impact and consequently the impact energy can also be retrieved complementarily.

# 5. Inverse analysis based on post-impact measurements

## 5.1. General considerations

In addition to real-time measurements, the wall conformation after impact may be easily obtained by employing techniques such as photogrammetry or Lidar. Impact-induced displacements with very good predictions may be obtained if the same type of data is collected before the impact. In such a situation, the number of data from the site would be much larger than that for real-time response. Typically, the position of the centre of gravity of all blocks can be derived from these data sets and will be considered in the following.

In this context, the aim of an inverse analysis could be to retrieve information related to the rock projectile kinematics at impact, such as its kinetic energy. This information could thus contribute to better quantifying the rockfall hazard at the wall location.

The development of the inverse analysis method is based on the following considerations. In case there is damage to concrete blocks, notably in the impact location vicinity, the determination of the blocks' gravity centre will not be possible. This eventuality should be accounted for in the development of the inverse analysis method. In addition, it is considered that the rock projectile volume and the impact
516

517

518

519

520

521

522

The output from the Bayesian inference-based inverse analysis is a set of ICP values retrieved as point estimates ‘maximum a posteriori’ and ‘mean’ (see Equation B.13). This predicted set of ICP values is then used to estimate the QoI (in this case damage to the wall or distance to the safety line). As a conventional approach, running a new NSCD model simulation with the predicted ICPs shall provide detailed information on the wall response from which the QoI can be extracted as a subset. However, this operation shall require about 30-40 minutes of computation time and added time for post-processing. This added delay time between on-site observation and warning system activation is unfavourable.

523

524

525

526

527

528

529

530

531

A procedure by which the QoI is instantly computed is developed to reduce the reaction time. The computational workflow components for this process are highlighted in the dashed box of Figure 12. The database of NSCD model simulation results is processed to evaluate the QoI as outputs. The resulting database of QoI is then used to create the PCE-based meta-model for each QoI. It serves as a surrogate of the NSCD model where the input ICPs are processed to compute the corresponding QoI output. The leave-one-out (LOO) error of order $10^{-2}$ is reported for the two PCE-based meta-models created for the two QoI (i.e., distance to safety and wall damage), named as Meta-$u^{safety}$ and Meta-$D_p$ (see section 3.2.3). The predicted ICPs are then fed into this surrogate of the NSCD model as input which is then evaluated to eventually predict the QoI and the whole computational process terminates.

532

533

534

535

536

537

Notably, in the absence of an exhaustive database of the MVP that shall be retrieved from the structure exposed to the natural rockfall events, the essential ‘evidence’ database (i.e., recorded MVP) to test the inverse analysis process is not available. Therefore, the concrete block velocity patterns from the same set of 300 model simulations and the new independent set of 300 simulations with different impact conditions are used as pseudo-evidences to test the process as the realisation of the natural rockfall events on the structure.

538

539

540

541

542

543

Each of these 600 known block velocity patterns is individually fed as input into the inverse analysis using the multiple model output feature of Bayesian analysis in UQlab (demonstrated in Equation B.8), providing the predicted ICP values set. The evaluation of the reliability of the inverse analysis is based on the comparison between the QoI reported via the ICP set used for running the simulation and the one estimated from evaluating the ICP set predicted by the inverse analysis onto the meta-model created for the same QoI.
location can be determined during the on-site visit. In this process, we assume that the NSCD model which is calibrated for a unique projectile mass is also able to predict the structure response for another rock projectile mass. In other words, this implies that for a given projectile kinetic energy, the wall response is the same whatever the mass and velocity. Under this assumption, the kinetic energy provided by the inverse analysis combined with the observed rock projectile mass gives an estimate of its velocity at impact.

The ICP comprise six physical parameters which describe the projectile kinematics and the point of impact on the wall. An accurate prediction of all individual parameters is of limited interest as the user is more focused on the estimation of the incident rock projectile kinetic energy. Therefore, it is proposed to reduce the dimension of the space of parameters to improve the predictions’ reliability. The set of six ICPs is reduced to three by combining $v$, $\Omega$, $\alpha$ and $\beta$ to calculate the total kinetic energy of the projectile ($KE_{proj}$) as per the Equation below and keeping position parameters $y$ and $z$ intact.

$$KE_{proj} = \frac{1}{2}m[v_{xy}]^2 + \frac{1}{2}m[v_z]^2 + \frac{1}{2}I[2\pi\Omega]^2 \quad (5)$$

Further, the total kinetic energy of the projectile can be divided into two components: one accounting for the energy on the plane parallel to the ground ($KE_p$) and the other accounting for the rest of the energy ($KE_{np}$). This is done by using only the planer component of the translational velocity ($v_{xy}$) component to compute $KE_p$ and the remaining two components (i.e., $v_z$ and $\Omega$) to compute $KE_{np}$.

This division is motivated by the observation that $KE_p$ mainly governs the wall displacement, as it is integrally transferred to the wall (except if the projectile rebound velocity is significant, which is not the case here). On the contrary, the rest of the energy (which is associated with projectile rotation and velocity along the vertical direction) is marginally transferred to the wall and thus has very little influence on the wall response.

This proposition is qualitatively explored by the Sobol indices evaluated on reduced ICPs for the maximum 3D displacement, presented in Figure 15. Here, $KE_p$ is observed highly influential for the resulting maximum displacement of the wall and as suspected, $KE_{np}$ provides negligible influence. The y-position of the projectile impact is observed as influential and supports the conclusion presented in the previous section. The position-z is observed to have a negligible influence which supports the conclusions in the previous section as well. Thereby, observing from the Sobol indices, the inverse analysis is conducted
544

4.3. Inverse analysis for warning purpose

545

First, the inverse analysis addresses the traffic interruption issue for which the post-impact minimum
distance of the wall to the safety line constitutes the QoI. The higher the distance, the higher the safety.

The predictive capability of the inverse analysis workflow presented in Figure 12 is tested for the maximum
concrete block velocity database obtained from the 600 NSCD model simulations. The MVPs from the
three impact cases mentioned in Section 3.1 are also used to predict the minimum distance to the safety
line.

546

547

548

549

550

551

The comparison between the predictions from the inverse analysis + meta-model and the simulation
results is presented in Figure 14. The relative proximity of the predictions to the ideal diagonal reflects
the potential of the presented methodology. Points above the diagonal concern cases for which the inverse
analysis overestimates the distance to the safety line. On the contrary, the inverse analysis underestimates
the distance to the safety lines for points below the diagonal, which is on the safe side.

552

553

554

555

556

The root mean square (RMS) error of the deviation from the expected response for the 600 simulations
is 16 cm. About 58% of this RMS error is due to the surrogate of the NSCD model (i.e., dashed box
component of the workflow). Besides, the presented process is also tested for the three illustrative impact
cases, reflecting a good prediction. In fact, in some sectors of the 6D space of the six ICPs, this surrogate
locally fails in precisely capturing the wall response, in particular due to the presence of non-linearities.
As an alternative, the reliability of the inverse analysis could be improved by running an NSCD model
simulation with the predicted ICPs in place of using the surrogate of the NSCD model, at the expense
of a longer reaction time upon event.

557

558

559

560

561

562

563

564

For practical purposes, it is proposed to classify the impact response into three classes depending on
the distance to the safety line: [0-0.75 m], [0.75-1.25 m] and [1.25-1.5m]. These classes correspond to
decreasing criticality concerning transport corridor safety, respectively high (H), medium (M) and low
(L) criticality. NSCD model simulation results indicate that 15%, 59% and 26% of the impact cases are
classified in class H, M and L respectively. Notably, the relatively high number of simulated cases are in
the middle category as a consequence of the Gaussian distribution of the magnitude of the velocity input
parameter. Predictions based on the inverse analysis + metamodel give values of 12%, 63% and 25% for
classes H, M and L respectively. Besides, the false negative predictions are estimated to be about 24%
and 10% among all in the H and M classes respectively. The remaining predictions i.e., good and false

565

566

567

568

569

570

571

572
<figure><img src="image_19.png" /><figcaption>Figure 15: Sobol sensitivity analysis of the impact condition parameters (ICPs) of derived and reduced dimensions (i.e., from 6 to 4) for the maximum 3D displacement</figcaption></figure>

taking into account the variability in the $KE_p$ and $y$ parameters only hereby referred to as reduced dimension ICPs.

## 5.2. Inverse analysis method

The process of inverse analysis for displacement-based prediction of the impact conditions follows nearly the same workflow as presented in the previous section and illustrated in Figure 11. Here, the set of 300 simulations is processed to generate the meta-model of the pattern of the displacement at rest of the gravity centre of all blocks (similar to the MVP).

The accuracy of the created meta-model was estimated of order $10^{-1}$ via leave-one-out-error. Besides, the created meta-model for 3D displacement is validated by using the same simulation set and an independent set of 300 simulations (thus a total of 600 simulations) and comparing the predicted vs true value (from NSCD model) at all 26 blocks in the movement zone, collectively presented in Figure 16. Here, the predictability coefficient ($Q^2$) presents the prediction accuracy, estimated for $N = 600 \times 26 = 15600$ observations as per Equation 4.

This created meta-model is used as a forward model for the Bayesian interface-based inverse analysis. The displacement data retrieved from real events shall serve as evidence for the inverse analysis, similar to the ‘on-site MVP’ mentioned in the previous section. Notably, here, two threshold levels of displacement magnitude are set at 0.2m and 0.05m for the recorded data from all 26 blocks. The first threshold of 0.2m is assigned to limit the number of blocks in the displacement pattern to avoid overfitting. The
<figure><img src="image_16.png" /><figcaption>Figure 14: Comparison of the minimum distance from the safety line, estimated from the NSCD model and the corresponding prediction from the inverse analysis tool accompanied with Meta-u$safety$. Here, the predictions are divided into three zones reflecting high (H), medium (M) and Low (L) levels of severity inclusive of the corresponding zonal false positive (+) and false negative (×) predictions. Results are based on a set of 600 impact conditions complemented with the impact cases presented in section 3.1 (IC-1 to IC-3).</figcaption></figure>

positives are in favour of the decision-making by the user, and hence are deemed acceptable.

In the end, the predictions are considered sufficiently good to be used in an operational context for estimating the consequences of the rockfall event in terms of structure displacement. The undeniable interest in this approach is that it automatically and remotely provides stakeholders with information of great value for decision-making upon impact. No quantitative evaluation would be derived rapidly from data measured on-site in the absence of such an approach.

### 4.4. Inverse analysis for structure damage quantification purpose

Second, the inverse analysis is used to remotely estimate damage to the wall. It is defined as the cumulative energy dissipated by the plastic damage of all the wall concrete blocks, $D_p$, after the wall is at rest. Damage to all concrete blocks of the wall is computed from NSCD model simulation results as per section 3.1.4. The approach considered for this inverse analysis is the same as that for the previous warning criterion. Figure 15 provides a comparison between $D_p$ values obtained from NSCD model simulations and the predicted ones. The RMS error of the deviation from the expected response is 82 kJ, which may be considered acceptable when considering all uncertainties associated with the quantification of plastic energy dissipation due to real rockfall impacts. Similar to the previous prediction of distance,
<figure><img src="image_20.png" /><figcaption>Figure 16: Validation of PCE-based meta-model of the 3D displacement of 26 blocks by comparison of the reported magnitude from 600 NSCD model simulations and the corresponding estimation by the meta-model</figcaption></figure>

second threshold of 0.05m is assigned to take into account a sufficient number of blocks for a small wall displacement situation where the first threshold limits the total number of moving blocks to maximum three.

The corresponding reduced dimension ICPs are predicted and the inverse analysis process terminates. Notably, the workflow corresponding to the evaluation of the QoI and the associated process to create a meta-model serving as a surrogate of the NSCD model is therefore not required.

### 5.3. Evaluation of the method

Following the guidelines of the previous section, the predictions of the energy transferred to the wall and the position of impact are compared with the corresponding NSCD model simulation results. These comparisons are presented in Figure 17 for the simulated 600 distinct impact conditions and the three illustrative impact conditions.

The energy transferred to the wall is predicted with an overall RMS error of 71 kJ (which is below the fracture energy of the concrete block) based on 600 simulations. The overall prediction is satisfactory for a wide range of energy magnitudes supporting the potential of the presented methodology. Nevertheless, two out of three illustrative cases are badly predicted with up to 120 kJ difference from the expected 520kJ for all cases (IC-1 and IC-3).

Further, the position of impact along the wall length is rather well predicted for impacts at convex and concave-wedges, including the illustrative impact tests in these two locations (IC-1 and IC-3). By
about 35% of the contribution to the RMS error is attributed to the surrogate of the NSCD model.

Besides, the process is tested for the three illustrative cases, reflecting a good prediction except for the impact at the angled-wedge (i.e., IC-2).

<figure><img src="image_17.png" /><figcaption>Figure 15: Comparison of the predicted damage to the wall, $D_p$, estimated from the NSCD model and the corresponding prediction from the inverse analysis tool accompanied with Meta-$D_p$. Here, the predictions are divided into three zones reflecting Low (L), medium (M) and high (H) levels of damage inclusive of the corresponding zonal false positive (+) and false negative (×) predictions. Results are based on a set of 600 impact conditions and impact cases presented in section 3.1 (IC-1 to IC-3).</figcaption></figure>

The cases in the right upper corner constitute the most critical impact cases. Points above the diagonal concern cases for which the inverse analysis overestimates the damage, which is on the safe side. For practical purposes, it is proposed to classify the damage into three classes depending on the value of $D_p$: [0-100 kJ], [100-300 kJ] and [300-800 kJ] for low (L), medium (M) or high (H) damage. These limits are arbitrarily defined from the distribution of $D_p$ over all cases.

The class limits refer to the level of damage to the wall corresponding to the fracture energy of one concrete block when a fracture propagates along a vertical plane in the middle of the concrete block. A classification based on the energy dissipation associated with the fracture of the considered concrete blocks following different patterns and mechanisms (compression or tension) and in different locations in each block (e.g., middle of the blocks and corners) would be more appropriate for the practitioners. However, such specific estimations are not possible in the framework of the NSCD method used for the definition of the NSCD model in the present work. Also, the experimental estimation of these specific energy dissipation magnitudes through dynamic tests is not in the scope of the present work. Thereby, in
<figure><img src="image_21.png" /><figcaption>Figure 17: Comparison of displacement-scan inverse analysis-based predictions of (a) energy transferred to the wall and (b) position along the wall length with the corresponding magnitude deduced from 600 different sets of the impact condition parameters.</figcaption></figure>

contrast, the predictions for the angled wedge zone (between 1.2 to 2.2 m approx.) are significantly less accurate thereby contributing to the overall RMS error of 37 cm. This significant difference in predictions accuracy is attributed to the zig-zag conformation of the wall, which amplifies the influence of the deviation angle ($\beta$) on the structure response depending on the impact location, as illustrated by the Sobol indices. In fact, considering a reduced ICPs input space dimension suppresses the distinct contribution of $\beta$ which in turn suppresses the variability of the impact at the angled wedge.

## 6. Conclusions

This work explored the potential of inverse analysis to obtain the information of interest in the operational context of a rockfall protection structure. The inverse analysis relies on simulation-based predictions of the structure response and makes use of the Bayesian interface and meta-modelling techniques. This work in particular focuses on rockfall protection walls made of piled-up and articulated concrete blocks for which an NSCD model was previously developed and calibrated against real-scale impact experiments.

The main conclusions drawn from this research are:

* Provided a well-calibrated model is available, it is possible to carry out reliable inverse analysis of
604

605

606

607

the absence of specific dynamic test results for these large and reinforced concrete blocks, results obtained
by Guo et al. [40] for a reinforced beam exposed to impact were extrapolated accounting for the difference
in block section and resulted in a value of 100 kJ approx. This value is a rough estimate of the first limit
for the considered concrete blocks.

608

609

610

611

612

613

614

According to NSCD model simulations results presented in Figure 15, 8%, 72% and 20% of the impact cases result in damage classified as H, M and L respectively. Notably, the relatively higher number of cases in the middle category are a consequence of the Gaussian distribution of the magnitude of the velocity input parameter. Values of 6%, 76% and 18% for classes H, M and L are obtained from the predictions respectively. Besides, the false negative predictions are estimated to be about 38% and 5% among all in the H and M classes respectively. The remaining predictions i.e., good and false positives are in favour of the decision-making by the user, and hence are acceptable.

615

616

617

618

619

The differences between simulation results and predictions are associated with points at a large dis-
tance from the diagonal (Fig. 15). Further method developments could improve the accuracy of these
predictions. Nevertheless, the predictions following the described inverse analysis are considered suffi-
ciently good to be used in an operational context to aid in remote decision-making, the RMS error is 82
kJ.

620

4.5. *Comments concerning the use of inverse analysis based on real-time measurements*

621

622

The results presented in section 4.3 and 4.4 demonstrate the feasibility of conducting an inverse
analysis of data that shall be collected on-site during impact, focusing on two specific purposes.

623

624

625

626

627

628

629

630

631

632

The proposed method, making use of numerical modelling, meta-modelling and statistical approaches,
revealed rather efficient in providing useful information which could not be derived from the direct
observation of the data that shall be received from the site. The duration of the whole process for
treating the input data typically lasts about five minutes with the developed script making it compatible
with rapid reaction upon event. The analysis scripts are written in Matlab as per UQLab framework
[41] allowing running this complex process in two operations (one for inverse analysis and the other for
output of interest prediction). The accuracy of the predictions is considered sufficiently good for remotely
classifying the events and aiding in decision-making. If required for another purpose, this accuracy could
be improved further, for example by also accounting for other measure types, such as deflection, or by
improving the used methodologies for creating the meta-models.
the response of complex structures subjected to localised dynamic loading.

* Investigating the structure response considering a wide range of impact conditions, representative of those encountered in situ is a prerequisite in view of defining the strategy for developing the inverse analysis method.

* The inverse analysis may be based on real-time or post-impact measurements and may serve various purposes, in particular for remote warning upon impact or rockfall activity survey.

* The presented inverse analysis framework is transferable to other types of such structures (such as flexible barriers) via the usage of the relevant recorded database.

As a perspective, the accuracy of the predictions by the inverse analysis will be improved, focusing on the meta-models which are key components in the proposed method. As a supplementary outcome, the recorded concrete block velocity database for different impact conditions can be further exploited towards the optimisation of the sensor location.

## References

[1] S. Lambert, F. Bourrier, Design of rockfall protection embankments: A review, Engineering Geology 154 (2013) 77–88.

[2] F. Calvetti, Rockfall shelters covered by granular layers, European Journal of Environmental and Civil Engineering (2011).

[3] C. Gentilini, G. Gottardi, L. Govoni, A. Mentani, F. Ubertini, Design of falling rock protection barriers using numerical models, Engineering Structures 50 (2013) 96–106. Engineering Structures: Modelling and Computations (special issue IASS-IACM 2012).

[4] F. Bourrier, J. Baroth, S. Lambert, Accounting for the variability of rock detachment conditions in designing rockfall protection structures, Natural Hazards (2016).

[5] A. Segalini, A. Valletta, A. Carri, R. Savi, Impact identification on flexible rockfall barriers: on site test of a wireless monitoring system, IOP Conference Series: Earth and Environmental Science 1124 (2023) 012125.
633

In real situations, the probability of occurrence of an impact in given conditions is extremely variable.

634

By contrast, and as mentioned in section 3.2.1, the ICP sets considered in this work were defined by not considering any possible correlation between the various parameters. The statistical assessment of the inverse analysis predictive capacities is related to the considered ranges and distributions of the six ICPs. It should thus be considered indicative at this stage. It should be improved by accounting for the probability of each scenario considered in this evaluation, for example considering site-specific data describing the rock projectile trajectory.

635

636

637

638

639

640

It is worth noting that, as the second step of the inverse analysis provides the set of ICP values resulting in an MVP similar to the observed one, a preliminary estimate of the projectile velocity at impact and consequently the impact energy can also be retrieved complementarily.

641

642

643

# 5. Inverse analysis based on post-impact measurements

644

## 5.1. General considerations

645

In addition to real-time measurements, the wall conformation after impact can be easily obtained by employing techniques such as photogrammetry or Lidar. An accurate impact-induced displacements may be obtained if the same type of data is collected before the impact. In such a situation, the number of data that shall be from the site would be much larger than that for real-time response. Typically, the position of the centre of gravity of all blocks can be derived from these data sets and will be considered in the following.

646

647

648

649

650

651

In this context, the aim of an inverse analysis could be to retrieve information related to the rock projectile kinematics at impact, such as its kinetic energy. This information could thus contribute to better quantifying the rockfall hazard at the wall location.

652

653

654

The development of the inverse analysis method is based on the following considerations. In case there is damage to concrete blocks, notably in the impact location vicinity, the determination of the blocks' gravity centre will not be possible. This eventuality is accounted for in the development of the inverse analysis method. In addition, it is considered that the rock projectile volume and the impact location can be determined during the on-site visit. For the simplicity of executing the process, we assume that the NSCD model which is calibrated for a unique projectile mass is also able to predict the structure response for another rock projectile mass. In other words, this implies that for a given projectile kinetic

655

656

657

658

659

660
[6] M. Ragnoli, A. Leoni, G. Barile, G. Ferri, V. Stornelli, Lora-based wireless sensors network for rockfall and landslide monitoring: A case study in pantelleria island with portable lorawan access, J. Low Power Electron. Appl. 12 (2022).

[7] P. Mayer, R. Rogge, A. Caviezel, J. Munch, A. Ringenbach, M. Magno, L. Benini, Design and evaluation of a lora controlled rugged multisensor unit for induced rockfall experiments, in: 2023 9th International Workshop on Advances in Sensors and Interfaces (IWASI), 2023, pp. 52–57. doi:10.1109/IWASI58316.2023.10164375.

[8] L. Castanon-Jano, E. Blanco-Fernandez, D. Castro-Fresno, D. Ferreño, Use of explicit fem models for the structural and parametrical analysis of rockfall protection barriers, Engineering Structures (2018).

[9] L. Dugelas, J. B. Coulibaly, F. Bourrier, S. Lambert, M.-A.nut, I. Olmedo, F. Nicot, Assessment of the predictive capabilities of discrete element models for flexible rockfall barriers, International Journal of Impact Engineering 133 (2019) 103365.

[10] Z.-H. Zhu, J.-H. Yin, J.-Q. Qin, D.-Y. Tan, A new discrete element model for simulating a flexible ring net barrier under rockfall impact comparing with large-scale physical model test data, Computers and Geotechnics (2019).

[11] J. P. Hambleton, O. Buzzi, A. Giacomini, M. Spadari, S. W. Sloan, Perforation of flexible rockfall barriers by normal block impact, Rock Mechanics and Rock Engineering (2013).

[12] Z. Yu, L. Luo, C. Liu, L. Guo, X. Qi, L. Zhao, Dynamic response of flexible rockfall barriers with different block shapes, Landslides (2021).

[13] A. Mentani, L. Govoni, G. Gottardi, S. Lambert, F. Bourrier, D. Toe, A new approach to evaluate the effectiveness of rockfall barriers, Procedia Engineering 158 (2016) 398–403. VI Italian Conference of Researchers in Geotechnical Engineering, CNRIG2016 - Geotechnical Engineering in Multidisciplinary Research: from Microscale to Regional Scale, 22-23 September 2016, Bologna (Italy).

[14] S. Lambert, D. Toe, A. Mentani, F. Bourrier, A Meta-Model-Based Procedure for Quantifying the On-Site Efficiency of Rockfall Barriers, Rock Mechanics and Rock Engineering 54 (2021) 487–500.
energy, the wall response is the same whatever the mass and velocity. Notably, this assumption is not tested via the NSCD simulations of varying mass of the projectile and hence we acknowledge that it can be considered a strong assumption. Under this assumption, the kinetic energy provided by the inverse analysis combined with the observed rock projectile mass gives an estimate of its velocity at impact. Thereby, the reliability of the prediction shall be relatively higher when the impacting rock mass in a natural event is close to the simulated projectile mass.

The ICP comprise six physical parameters which describe the projectile kinematics and the point of impact on the wall. The accurate prediction of all parameters is of limited interest for post-impact measurements as the user is more focused on the estimation of the incident rock projectile kinetic energy. Therefore, it is proposed to use a reduced dimension of the input parameters space to create the metamodel from the 3D displacement database. The set of six ICPs is reduced to three by combining $v$, $\Omega$, $\alpha$ and $\beta$ to calculate the total kinetic energy of the projectile ($KE_{proj}$) as per the Equation 5 below and keeping position parameters $y$ and $z$ intact.

$$KE_{proj} = \frac{1}{2}m[v_{xy}]^2 + \frac{1}{2}m[v_z]^2 + \frac{1}{2}I[2\pi\Omega]^2 \quad (5)$$

Further, the total kinetic energy of the projectile can be divided into two components: one accounting for the energy on the plane parallel to the ground ($KE_p$) and the other accounting for the rest of the energy ($KE_{np}$). This is done by using only the planer component of the translational velocity ($v_{xy}$) component to compute $KE_p$ and the remaining two components (i.e., $v_z$ and $\Omega$) to compute $KE_{np}$.

This division is motivated by the observation that $KE_p$ mainly governs the wall displacement, as it is integrally transferred to the wall (except if the projectile rebound velocity is significant, which is not the case here). On the contrary, the rest of the energy (which is associated with projectile rotation and velocity along the vertical direction) is marginally transferred to the wall and thus has very little influence on the wall response.

This proposition is qualitatively explored by the Sobol indices evaluated on reduced ICPs for the maximum 3D displacement, presented in Figure 16. Here, $KE_p$ is observed highly influential for the resulting maximum displacement of the wall and as suspected, $KE_{np}$ provides negligible influence. The y-position of the projectile impact is observed as influential and supports the conclusion presented in the previous section. The position-z is observed to have a negligible influence which supports the conclusions
670

671

[15] J. Escallón, C. Wendeler, E. Chatzi, P. Bartelt, Parameter identification of rockfall protection barrier components through an inverse formulation, Engineering Structures 77 (2014) 1–16.

672

673

674

[16] R. Gupta, F. Bourrier, V. Acary, S. Lambert, Bayesian interface based calibration of a novel rockfall protection structure modelled in the non-smooth contact dynamics framework, Engineering structures 297 (2023).

675

676

677

[17] A. Furet, S. Lambert, P. Villard, J. P. Jarrin, Experimental and numerical impact responses of an innovative rockfall protection structure made of articulated concrete blocks, Rock Mechanics and Rock Engineering (2022).

678

679

[18] S. Lambert, F. Bourrier, P. Gotteland, F. Nicot, An experimental investigation of the response of slender protective structures to rockfall impacts, Canadian geotechnical journal (2019).

680

681

[19] O. Korini, M. Bost, J.-P. Rajot, Y. Bennani, N. Freitag, The influence of geosynthetics design on the behavior of reinforced soil embankments subjected to rockfall impacts, Engineering Geology (2021).

682

[20] EOTA, Ead 340059-00-0106: Falling rock protections kits, 2018.

683

684

685

[21] V. Acary, F. Perignon, Siconos: A Software Platform for Modeling, Simulation, Analysis and Control of Nonsmooth Dynamical Systems, SIMULATION NEWS EUROPE, ArgeSIM/ASIM 17 (2007) 19–26.

686

687

688

[22] V. Acary, Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact, ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik 96 (2016) 585–603.

689

690

691

[23] S. Lambert, R. Gupta, F. Bourrier, V. Acary, What can we learn from simulation-based quantification of energy dissipation in rockfall protection structures? case of an articulated wall modelled with the nscd method, Rock Mechanics and Rock Engineering (submitted).

692

693

694

[24] D. Toe, A. Mentani, L. Govoni, F. Bourrier, G. Gottardi, S. Lambert, Introducing meta-models for a more efficient hazard mitigation strategy with rockfall protection barriers, Rock Mechanics and Rock Engineering 51 (2018) 1097–1109.

695

696

[25] M. Previtali, M. O. Ciantia, S. Spadea, R. Castellanza, G. Crosta, Assessing rockfall barrier performance through block propagation codes and meta-models, in: M. Barla, A. Di Donna, D. Sterpi,
<figure><img src="image_18.png" /><figcaption>Figure 16: Sobol sensitivity analysis of the impact condition parameters (ICPs) of derived and reduced dimensions (i.e., from 6 to 4) for the maximum 3D displacement</figcaption></figure>

in the previous section as well. Thereby, observing from the Sobol indices, the inverse analysis is conducted
taking into account the variability in the input $KE_p$ and $y$ parameters only hereby referred to as reduced
dimension ICPs.

## 5.2. Inverse analysis method

The process of inverse analysis for displacement-based prediction of the impact conditions follows the same workflow as illustrated in Figure 12, however by replacing some components detailed hereafter. The recorded database from 300 NSCD model simulations is associated with the reduced dimension ICPs, instead of 6 ICPs. The 3D displacement post-impact is evaluated at all blocks in the movement zone and the database of the displacement pattern (as a replacement for the MVP mentioned in the previous section) corresponding to each set of reduced dimension ICPs is obtained. A PCE-based metamodel is created to associate each reduced dimension ICP set to a distinct 3D-displacement pattern, named Meta-3DP.

The accuracy of the Meta-3DP is estimated via the leave-one-out error, reported of order $10^{-1}$ for all 26 blocks in the movement zone collectively making a 3D-displacement pattern. Besides, the created meta-model for 3D displacement is validated by using the same simulation set and the independent set of 300 simulations (presented in Section 4.2), thus a total of 600 simulations, and comparing the predicted vs true value (from NSCD model) at all 26 blocks in the movement zone, collectively presented in Figure 17. Here, the predictability coefficient ($Q^2$) presents the prediction accuracy, estimated for $N = 600 \times 26 = 15600$ observations as per Equation 4.
697

A. Insana (Eds.), Challenges and Innovations in Geomechanics, Springer International Publishing, Cham, 2023, pp. 291–298.

698

699

[26] F. Bourrier, F. Berger, P. Tardif, L. Dorren, O. Hungr, Rockfall rebound: comparison of detailed field experiments and alternative modelling approaches, Earth Surface Processes and Landforms 37 (2012) 656–665.

700

701

702

[27] M. Spadari, A. Giacomini, O. Buzzi, J. Hambleton, Prediction of the bullet effect for rockfall barriers: a scaling approach, Rock Mechanics and Rock Engineering (2012).

703

704

[28] A. Breugnot, S. Lambert, P. Villard, P. Gotteland, A discrete/continuous coupled approach for modeling impacts on cellular geostructures, Rock Mechanics and Rock Engineering (2016).

705

706

[29] C. Lataniotis, S. Marelli, B. Sudret, UQLab user manual – The Input module, Technical Report, Chair of Risk, Safety & Uncertainty Quantification, ETH Zurich, 2015. Report UQLab-V0.9-102.

707

708

[30] I. M. Sobol, A screening design for factorial experiments with interactions, Mathematical and Computer Modelling 1 (1993) 407–414.

709

710

[31] B. Sudret, Global sensitivity analysis using polynomial chaos expansions, Reliability Engineering & System Safety 93 (2008) 964–979. Bayesian Networks in Dependability.

711

712

[32] E. Rouzies, C. Lauvernet, B. Sudret, A. Vidard, How is a global sensitivity analysis of a catchment-scale, distributed pesticide transfer model performed? application to the peshmelba model, Geosci-entific Model Development 16 (2023) 3137–3163.

713

714

715

[33] S. Marelli, B. Sudret, Uqlab user manual - polynomial chaos expansions, 2015. doi:10.13140/RG.2.1.3778.7366.

716

717

[34] S. Marelli, C. Lamas, K. Konakli, C. Mylonas, P. Wiederkehr, B. Sudret, UQLab user manual – Sensitivity analysis, Technical Report, Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland, 2022. Report UQLab-V2.0-106.

718

719

720

[35] M. A. Hariri-Ardebili, B. Sudret, Polynomial chaos expansion for uncertainty quantification of dam engineering problems, Engineering Structures 203 (2020) 109631.

721

722

[36] X. Guo, D. Dias, C. Carvajal, L. Peyras, P. Breul, Reliability analysis of embankment dam sliding stability using the sparse polynomial chaos expansion, Engineering Structures 174 (2018) 295–307.

723
<figure><img src="image_19.png" /><figcaption>Figure 17: Validation of PCE-based meta-model of the 3D displacement pattern (Meta-3DP) of 26 blocks by comparison of the reported magnitude from 600 NSCD model simulations and the corresponding estimation by the meta-model</figcaption></figure>

This created meta-model is used as a forward model for the Bayesian inference-based inverse analysis.

The displacement data retrieved from real events shall serve as evidence for the inverse analysis, by replacing the ‘on-site MVP’ mentioned in the workflow of Figure 12. Notably, here, two threshold levels of displacement magnitude are considered at 0.2m and 0.05m based on the observation from processed displacement data from 300 simulations for all 26 blocks. The first threshold of 0.2m is assigned to limit the number of blocks in the displacement pattern to avoid overfitting. The second threshold of 0.05m is assigned to take into account a sufficient number of blocks for a small wall displacement situation where the first threshold limits the total number of moving blocks to a maximum of three.

The corresponding reduced dimension ICPs are predicted and the inverse analysis process terminates.

Notably, the workflow (in the dashed box of Figure 12) corresponding to the evaluation of the QoI and the associated process to create a meta-model serving as a surrogate of the NSCD model is therefore not required.

### 5.3. Evaluation of the method

Following the guidelines of the previous section, the predictions of the energy transferred to the wall and the position of impact are compared with the corresponding NSCD model simulation results. These comparisons are presented in Figure 18 for the simulated 600 distinct impact conditions and the three illustrative impact conditions.
[37] T. Bayes, N. Price, LII. An essay towards solving a problem in the doctrine of chances. by the late Rev. Mr. Bayes, F. R. S. communicated by Mr. Price, in a letter to John Canton, A. M. F. R. S, Philosophical Transactions of the Royal Society of London 53 (1763) 370–418.

[38] P.-R. Wagner, J. Nagel, S. Marelli, B. Sudret, UQLab user manual – Bayesian inversion for model calibration and validation, Technical Report, Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland, 2022. Report UQLab-V2.0-113.

[39] J. Guo, J. Cai, Q. Chen, X. Liu, Y. Wang, Z. Zuo, Dynamic behaviour and energy dissipation of reinforced recycled aggregate concrete beams under impact, Construction and Building Materials (2019).

[40] S. Marelli, B. Sudret, UQLab: A Framework for Uncertainty Quantification in Matlab, 2014, pp. 2554–2563. doi:10.1061/9780784413609.257.

## Appendix A. PCE based meta-model

Consider the articulated structure model represented by a $M(\mathbf{X})$ as an equivalent mathematical model. Here, $\mathbf{X} \in \mathbb{R}^M$ is a random vector with independent components described by the joint probability density function (PDF) $f_{\mathbf{X}}$. Consider also a finite variance computational model as a map $Y = M(\mathbf{X})$, with $Y \in \mathbb{R}$ such that:

$$ \mathbb{E}[Y^2] = \int_{D_{\mathbf{X}}} M(\mathbf{x})^2 f_{\mathbf{X}}(\mathbf{x}) d\mathbf{x} < \infty \quad (A.1) $$

Then, under the assumption of Equation A.1, the PCE of $M(\mathbf{X})$ is defined as:

$$ Y = M(\mathbf{X}) = \sum_{\alpha \in \mathbb{N}^M} y_{\alpha} \Psi_{\alpha}(\mathbf{X}) \quad (A.2) $$

where, the $\Psi_{\alpha}(\mathbf{X})$ are multivariate polynomials orthonormal with respect to $f_{\mathbf{X}}$, $\alpha \in \mathbb{N}^M$ is a multi-index that identifies the components of the multivariate polynomials $\Psi_{\alpha}$ and the $y_{\alpha} \in \mathbb{R}$ are the corresponding coefficients. In practical applications, the sum in Equation A.2 needs to be truncated to a finite sum by introducing the truncated polynomial chaos expansion:

$$ M(\mathbf{X}) \approx M^{PC}(\mathbf{X}) = \sum_{\alpha \in \mathcal{A}} y_{\alpha} \Psi_{\alpha}(\mathbf{X}) \quad (A.3) $$
<figure><img src="image_20.png" /><figcaption>Figure 18: Comparison of displacement-scan inverse analysis-based predictions of (a) energy transferred to the wall and (b) position along the wall length with the corresponding magnitude deduced from 600 different sets of the impact condition parameters.</figcaption></figure>

The energy transferred to the wall is predicted with an overall RMS error of 71 kJ (which is below the fracture energy of the concrete block) based on 600 simulations. The overall prediction is satisfactory for a wide range of energy magnitudes supporting the potential of the presented methodology. Nevertheless, two out of three illustrative cases are badly predicted with up to 120 kJ difference from the expected 520kJ for all cases (IC-1 and IC-3).

Further, the position of impact along the wall length is rather well predicted for impacts at convex and concave-wedges, including the illustrative impact tests in these two locations (IC-1 and IC-3). By contrast, the predictions for the angled wedge zone (between 1.2 to 2.2 m approx.) are significantly less accurate thereby contributing to the overall RMS error of 37 cm (which is about 24% of the length of a block). This significant difference in prediction accuracy is attributed to the zig-zag conformation of the wall, which amplifies the influence of the deviation angle ($\beta$) on the structure response depending on the impact location, as illustrated by the Sobol indices. Considering a reduced ICPs input space dimension suppresses the distinct contribution of $\beta$ which in turn suppresses the variability of the impact at the angled wedge.

Overall, the inverse analysis method based on the post-impact measurements is deemed acceptable to retrieve the incident projectile impact position and the energy transferred to the wall. The presented
where, $\mathcal{A} \subset \mathbb{N}^M$ is the set of selected multi-indices of multivariate polynomials.

In this work, the least-angle regression (LARS) method is used to create the PCE meta-model truncated to the maximum polynomial degree ($p$) ranging from 1 to 20, and using hyperbolic truncation scheme ($q$) ranging from 0.5 to 1.

$$ \mathcal{A}^{M,p,q} = \{\alpha \in \mathcal{A}^{M,p} : \|\alpha\|_q \le p\}, \quad \text{where, } \|\alpha\|_q = \left(\sum_{i=1}^{M} \alpha_i^q\right)^{1/q} \qquad (A.4) $$

The accuracy of the constructed PCE is estimated by computing the leave-one-out (LOO) cross-validation error ($\epsilon_{LOO}$). It consists in building $N$ meta-models $\mathcal{M}^{PC\setminus i}$, each one created on a reduced experimental design $\mathcal{X}\setminus\mathbf{x}^{(i)} = \{\mathbf{x}^{(j)}, j = 1, \dots, N, j \ne i\}$ and comparing its prediction on the excluded point $\mathbf{x}^{(i)}$ with the real value $y^{(i)} = \mathcal{M}(\mathbf{x}^{(i)})$ [33]. The leave-one-cross-validation error can be written as:

$$ \epsilon_{LOO} = \frac{\sum_{i=1}^{N} \left( \mathcal{M}(\mathbf{x}^{(i)} - \mathcal{M}^{PC\setminus i}(\mathbf{x}^{(i)})) \right)^2}{\sum_{i=1}^{N} \left( \mathcal{M}(\mathbf{x}^{(i)} - \hat{\mu}_Y) \right)^2} \qquad (A.5) $$

where, $\hat{\mu}_Y$ is the mean of the experimental design sample.

## Appendix B. Bayesian interface for inverse analysis

Consider the computational model $\mathcal{M}$ that allows the analyst to predict certain quantities of interest gathered in a vector $\mathbf{y} \in \mathbb{R}^{N_{out}}$ as a function of input parameters $\mathbf{x}$:

$$ \mathcal{M}: \mathbf{x} \in \mathcal{D}_{\mathbf{X}} \subset \mathbb{R}^M \mapsto \mathbf{y} = \mathcal{M}(\mathbf{x}) \in \mathbb{R}^{N_{out}} \qquad (B.1) $$

The Bayesian interface for inverse analysis focuses on identifying the input parameters of a computational model to recover the observations in the collected output data set. It comprises of a computational forward model $\mathcal{M}$, a set of input parameters $\mathbf{x} \in \mathcal{D}_{\mathbf{X}}$ that need to be inferred, and a set of experimental data $\mathcal{Y}$. Here, $\mathcal{Y} \stackrel{\text{def}}{=} \{\mathbf{y}_1, \dots, \mathbf{y}_N\}$ is a global data set of $N$ independent measured quantities of interest ($\mathbf{y}_i$).

The forward model $\mathbf{x} \mapsto \mathcal{M}(\mathbf{x})$ is a mathematical representation of the system under consideration. The lack of knowledge on the input parameters is modelled by considering them as a random vector, denoted by $\mathbf{X}$ which is assumed to follow a so-called prior distribution (with support $\mathcal{D}_{\mathbf{X}}$), as presented in Figure 7 in the present work.
740

741

method well captures the variability in the retrieved information. Further efforts can be made to improve
the predictive capability of the inverse analysis method.

742

# 6. Conclusions

743

744

745

746

747

This work explored the potential of inverse analysis to obtain the information of interest in the opera-
tional context of a rockfall protection structure. The inverse analysis relies on simulation-based evaluation
of the structure response and makes use of the Bayesian inference and meta-modelling techniques. This
work in particular focuses on rockfall protection walls made of piled-up and articulated concrete blocks
for which an NSCD model was previously developed and calibrated against real-scale impact experiments.

748

The main conclusions drawn from this research are :

749

750

* Provided a well-calibrated model is available, it is possible to carry out reliable inverse analysis of the response of complex structures subjected to localised dynamic loading.

751

752

753

* Investigating the structure response considering a wide range of impact conditions, representative of those encountered in situ is a prerequisite in view of defining the strategy for developing the inverse analysis method.

754

755

* The inverse analysis may be based on real-time or post-impact measurements and may serve various purposes, in particular for remote warning upon impact or rockfall activity survey.

756

757

758

* The presented inverse analysis framework is transferable to other types of such structures (such as flexible barriers) via the usage of the representative numerical model and the relevant recorded database on-site.

759

760

761

* The inverse analysis process can be implemented to other wall configurations of the same constituent blocks and connectors considering the same calibrated model parameters and following the workflow presented in Figure 12.

762

763

764

765

766

The purpose of the two presented inverse analysis processes (i.e., for real-time impact and post-impact conditions) is different and hence a direct comparative assessment is not representative. Nonetheless, as both analyses yield predictions of one distance-based parameter and one energy-based parameter, an indirect comparative assessment is made that can be useful in future research. The displacement pattern metamodel is about 8.3% less accurate than the velocity pattern metamodel accessed via the reported $Q^2$
$$X \sim \pi(x) \qquad (B.2)$$

The Bayesian statistics combine this prior knowledge of the parameters with the few observed data points to obtain a statistical model called posterior distribution ($\pi(x | y)$) of the input parameters, using Bayes' theorem [37], expressed as:

$$\pi(x | y) = \frac{\pi(y | x)\pi(x)}{\pi(y)} \qquad (B.3)$$

Now, considering the available data set ($\mathcal{Y}$) as independent realizations of $Y | x \sim \pi(y | x)$, the collected measurements result in the definition of the likelihood function $\mathcal{L}(x; \mathcal{Y})$, which is a function of input parameters $x$:

$$\mathcal{L}: x \mapsto \mathcal{L}(x; \mathcal{Y}) \stackrel{\text{def}}{=} \prod_{i=1}^{N} \pi(y_i | x) \qquad (B.4)$$

This implicitly assumes in/dependence between individual measurements in $\mathcal{Y}$. Intuitively the likelihood function for a given $x$ returns the relative likelihood of observing the data at hand, under the assumption that it follows the prescribed parametric distribution $\pi(y | x)$.

As all models are simplifications of the real world, a discrepancy term ($\epsilon$) is introduced to connect real-world observations ($\mathcal{Y}$) to the predictions by the model. In practice, the discrepancy term represents the effects of the measurement error and model inaccuracy. The discrepancy term introduced here reads:

$$y = \mathcal{M}(x) + \epsilon \qquad (B.5)$$

Here, the $\epsilon$ is assumed as an additive Gaussian discrepancy [38] with a zero mean and given covariance matrix ($\Sigma$):

$$\epsilon \sim \mathcal{N}(\epsilon | 0, \Sigma) \qquad (B.6)$$

Taking insights from the discrepancy term definition, a particular measurement point ($y_i \in \mathcal{Y}$, is a realization of the Gaussian distribution with mean value $\mathcal{M}(x)$ and covariance matrix $\Sigma$. This distribution is named as discrepancy model and is expressed as:

$$\pi(y | x) = \mathcal{N}(y | \mathcal{M}(x), \Sigma) \qquad (B.7)$$
values. The resulting energy-based prediction is 13% more accurate for the former whereas the distance-based prediction relative to the corresponding block dimension (i.e., length for position of impact and width for minimum distance to the safety line) is reported 11% less accurate for the former when accessed via the reported RMS errors.

As a perspective, the predictive capability of the proposed methodology shall be improved by focusing on the accuracy of the individual components of the whole process. The improvement in the surrogate model via an increase in the $Q^2$ coefficient and improvement in the inverse analysis method via (1) reducing the RMS error of the predictions and (2) reducing the number of false indicators cases ‘especially negative’ shall be adopted. As a supplementary outcome, the recorded concrete block velocity database for different impact conditions can be further exploited towards the optimisation of the sensor location.

Moreover, the present method is limited to considering a single rock impact and therefore the assessment of multiple rocks impacting the wall at the same time is subjected to future research.

## References

[1] S. Lambert, F. Bourrier, Design of rockfall protection embankments: A review, Engineering Geology 154 (2013) 77–88.

[2] F. Calvetti, Rockfall shelters covered by granular layers, European Journal of Environmental and Civil Engineering (2011).

[3] C. Gentilini, G. Gottardi, L. Govoni, A. Mentani, F. Ubertini, Design of falling rock protection barriers using numerical models, Engineering Structures 50 (2013) 96–106. Engineering Structures: Modelling and Computations (special issue IASS-IACM 2012).

[4] F. Bourrier, J. Baroth, S. Lambert, Accounting for the variability of rock detachment conditions in designing rockfall protection structures, Natural Hazards (2016).

[5] A. Segalini, A. Valletta, A. Carri, R. Savi, Impact identification on flexible rockfall barriers: on site test of a wireless monitoring system, IOP Conference Series: Earth and Environmental Science 1124 (2023) 012125.

[6] M. Ragnoli, A. Leoni, G. Barile, G. Ferri, V. Stornelli, Lora-based wireless sensors network for
In application, the discrepancy model defines the connection between the supplied data (Y) and the forward model. In the present work, as the inverse analysis is conducted on the pseudo on-site measurement (as evidence), the discrepancy model with known residual variance is assigned.

Afterwards, the N independent available measurements (in our case, maximum velocity pattern and 3D displacement databases) gathered in the data-set (i.e., $\mathcal{Y} = \{\mathbf{y}_1, ..., \mathbf{y}_N\}$) are used to define the likelihood function as:

$$ 
\begin{aligned}
\mathcal{L}(\mathbf{x}; \mathcal{Y}) &= \prod_{i=1}^{N} \mathcal{N}(\mathbf{y}_i | \mathcal{M}(\mathbf{x}), \mathbf{\Sigma}) \\
&= \prod_{i=1}^{N} \frac{1}{\sqrt{(2\pi)^{N_{\text{out}}} \det(\mathbf{\Sigma})}} \exp \left( -\frac{1}{2} (\mathbf{y}_i - \mathcal{M}(\mathbf{x}))^{\top} \mathbf{\Sigma}^{-1} (\mathbf{y}_i - \mathcal{M}(\mathbf{x})) \right)
\end{aligned}
\quad (B.8) $$

Thereby, following Bayes' theorem, the posterior distribution $\pi(\mathbf{x} | \mathcal{Y})$ of the parameters ($\mathbf{x}$) given the observations in $\mathcal{Y}$ can be written as:

$$ \pi(\mathbf{x} | \mathcal{Y}) = \frac{\mathcal{L}(\mathbf{x}; \mathcal{Y})\pi(\mathbf{x})}{Z} \quad (B.9) $$

Here, $Z$ is a normalizing factor, known as the marginal likelihood or evidence, is added to ensure that this distribution integrates to 1:

$$ Z \stackrel{\text{def}}{=} \int_{D_{\mathbf{x}}} \mathcal{L}(\mathbf{x}; \mathcal{Y}) \pi(\mathbf{x}) d\mathbf{x} \quad (B.10) $$

The closed-form solutions do not exist in practice, and therefore the posterior distribution is obtained through Markov chain Monte Carlo (MCMC) simulations. In the present work, among many proposed algorithms in [38], the Adaptive Metropolis (AM) algorithm is used and 100 parallel chains with 1000 steps are assigned to the MCMC solver. The start of the covariance adaptation in AM algorithm is assigned at the 100th step (see Wagner et al. [38] for more details).

The output ($\mathbf{y}$) predictive capabilities of the Bayesian interface is assessed through the comparison of prior ($\pi(\mathbf{y})$) and posterior ($\pi(\mathbf{y} | \mathcal{Y})$) output distributions as:

$$ \pi(\mathbf{y}) = \int_{D_{\mathbf{x}}} \pi(\mathbf{y} | \mathbf{x})\pi(\mathbf{x})d\mathbf{x} \quad (B.11) $$

$$ \pi(\mathbf{y} | \mathcal{Y}) = \int_{D_{\mathbf{x}}} \pi(\mathbf{y} | \mathbf{x})\pi(\mathbf{x} | \mathcal{Y})d\mathbf{x} \quad (B.12) $$
rockfall and landslide monitoring: A case study in pantelleria island with portable lorawan access,
J. Low Power Electron. Appl. 12 (2022).

[7] P. Mayer, R. Rogge, A. Caviezel, J. Munch, A. Ringenbach, M. Magno, L. Benini, Design and evaluation of a lora controlled rugged multisensor unit for induced rockfall experiments, in: 2023 9th International Workshop on Advances in Sensors and Interfaces (IWASI), 2023, pp. 52–57. doi:10.1109/IWASI58316.2023.10164375.

[8] L. Castanon-Jano, E. Blanco-Fernandez, D. Castro-Fresno, D. Ferreño, Use of explicit fem models for the structural and parametrical analysis of rockfall protection barriers, Engineering Structures (2018).

[9] L. Dugelas, J. B. Coulibaly, F. Bourrier, S. Lambert, M.-A.nut, I. Olmedo, F. Nicot, Assessment
of the predictive capabilities of discrete element models for flexible rockfall barriers, International
Journal of Impact Engineering 133 (2019) 103365.

[10] Z.-H. Zhu, J.-H. Yin, J.-Q. Qin, D.-Y. Tan, A new discrete element model for simulating a flexible ring net barrier under rockfall impact comparing with large-scale physical model test data, Computers and Geotechnics (2019).

[11] J. P. Hambleton, O. Buzzi, A. Giacomini, M. Spadari, S. W. Sloan, Perforation of flexible rockfall barriers by normal block impact, Rock Mechanics and Rock Engineering (2013).

[12] Z. Yu, L. Luo, C. Liu, L. Guo, X. Qi, L. Zhao, Dynamic response of flexible rockfall barriers with different block shapes, Landslides (2021).

[13] A. Mentani, L. Govoni, G. Gottardi, S. Lambert, F. Bourrier, D. Toe, A new approach to evaluate the effectiveness of rockfall barriers, Procedia Engineering 158 (2016) 398–403. VI Italian Conference of Researchers in Geotechnical Engineering, CNRIG2016 - Geotechnical Engineering in Multidisciplinary Research: from Microscale to Regional Scale, 22-23 September 2016, Bologna (Italy).

[14] S. Lambert, D. Toe, A. Mentani, F. Bourrier, A Meta-Model-Based Procedure for Quantifying the On-Site Efficiency of Rockfall Barriers, Rock Mechanics and Rock Engineering 54 (2021) 487–500.

[15] J. Escallón, C. Wendeler, E. Chatzi, P. Bartelt, Parameter identification of rockfall protection barrier components through an inverse formulation, Engineering Structures 77 (2014) 1–16.
Lastly, in the present studies, the purpose of the Bayesian interface is to obtain the ‘best set of impact condition parameters (ICPs)’. Given the posterior distribution ($\pi(x | \mathcal{Y})$), we are interested in finding a suitable set among the posterior computed set i.e., $X | \mathcal{Y}$. This is done through a point estimator ($\hat{x}$) computed from:

$$ \pi(y | \mathcal{Y}) \stackrel{\text{def}}{=} \pi(y | \hat{x}) \qquad (\text{B.13}) $$

This point estimator can be a mean or mode (maximum a posteriori ‘MAP’) [38] of the posterior distribution as per the user’s choice.
[16] R. Gupta, F. Bourrier, V. Acary, S. Lambert, Bayesian interface based calibration of a novel
rockfall protection structure modelled in the non-smooth contact dynamics framework, Engineering
structures 297 (2023).

[17] A. Furet, S. Lambert, P. Villard, J. P. Jarrin, Experimental and numerical impact responses of an
innovative rockfall protection structure made of articulated concrete blocks, Rock Mechanics and
Rock Engineering (2022).

[18] S. Lambert, F. Bourrier, P. Gotteland, F. Nicot, An experimental investigation of the response of
slender protective structures to rockfall impacts, Canadian geotechnical journal (2019).

[19] O. Korini, M. Bost, J.-P. Rajot, Y. Bennani, N. Freitag, The influence of geosynthetics design on the
behavior of reinforced soil embankments subjected to rockfall impacts, Engineering Geology (2021).

[20] EOTA, Ead 340059-00-0106: Falling rock protections kits, 2018.

[21] V. Acary, F. Perignon, Siconos: A Software Platform for Modeling, Simulation, Analysis and Control
of Nonsmooth Dynamical Systems, SIMULATION NEWS EUROPE, ArgeSIM/ASIM 17 (2007) 19–
26.

[22] F. Dubois, V. Acary, M. Jean, The Contact Dynamics method: A nonsmooth story, Comptes
Rendus Mécanique 346 (2018) 247–262.

[23] V. Acary, Energy conservation and dissipation properties of time-integration methods for nonsmooth
elastodynamics with contact, ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift
für Angewandte Mathematik und Mechanik 96 (2016) 585–603.

[24] S. Lambert, R. Gupta, F. Bourrier, V. Acary, What can we learn from simulation-based quantifi-
cation of energy dissipation in rockfall protection structures? case of an articulated wall modelled
with the nscd method, Rock Mechanics and Rock Engineering (submitted).

[25] D. Toe, A. Mentani, L. Govoni, F. Bourrier, G. Gottardi, S. Lambert, Introducing meta-models for
a more efficient hazard mitigation strategy with rockfall protection barriers, Rock Mechanics and
Rock Engineering 51 (2018) 1097–1109.

[26] M. Previtali, M. O. Ciantia, S. Spadea, R. Castellanza, G. Crosta, Assessing rockfall barrier perfor-
mance through block propagation codes and meta-models, in: M. Barla, A. Di Donna, D. Sterpi,
847

A. Insana (Eds.), Challenges and Innovations in Geomechanics, Springer International Publishing, Cham, 2023, pp. 291–298.

848

849

[27] F. Bourrier, F. Berger, P. Tardif, L. Dorren, O. Hungr, Rockfall rebound: comparison of detailed field experiments and alternative modelling approaches, Earth Surface Processes and Landforms 37 (2012) 656–665.

850

851

852

[28] M. Spadari, A. Giacomini, O. Buzzi, J. Hambleton, Prediction of the bullet effect for rockfall barriers: a scaling approach, Rock Mechanics and Rock Engineering (2012).

853

854

[29] A. Breugnot, S. Lambert, P. Villard, P. Gotteland, A discrete/continuous coupled approach for modeling impacts on cellular geostructures, Rock Mechanics and Rock Engineering (2016).

855

856

[30] C. Lataniotis, S. Marelli, B. Sudret, UQLab user manual – The Input module, Technical Report, Chair of Risk, Safety & Uncertainty Quantification, ETH Zurich, 2015. Report UQLab-V0.9-102.

857

858

[31] I. M. Sobol, A screening design for factorial experiments with interactions, Mathematical and Computer Modelling 1 (1993) 407–414.

859

860

[32] B. Sudret, Global sensitivity analysis using polynomial chaos expansions, Reliability Engineering & System Safety 93 (2008) 964–979. Bayesian Networks in Dependability.

861

862

[33] E. Rouzies, C. Lauvernet, B. Sudret, A. Vidard, How is a global sensitivity analysis of a catchment-scale, distributed pesticide transfer model performed? application to the peshmelba model, Geosci-entific Model Development 16 (2023) 3137–3163.

863

864

865

[34] S. Marelli, B. Sudret, Uqlab user manual - polynomial chaos expansions, 2015. doi:10.13140/RG.2.1.3778.7366.

866

867

[35] S. Marelli, C. Lamas, K. Konakli, C. Mylonas, P. Wiederkehr, B. Sudret, UQLab user manual – Sensitivity analysis, Technical Report, Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland, 2022. Report UQLab-V2.0-106.

868

869

870

[36] M. A. Hariri-Ardebili, B. Sudret, Polynomial chaos expansion for uncertainty quantification of dam engineering problems, Engineering Structures 203 (2020) 109631.

871

872

[37] X. Guo, D. Dias, C. Carvajal, L. Peyras, P. Breul, Reliability analysis of embankment dam sliding stability using the sparse polynomial chaos expansion, Engineering Structures 174 (2018) 295–307.

873
[38] T. Bayes, N. Price, LII. An essay towards solving a problem in the doctrine of chances. by the late Rev. Mr. Bayes, F. R. S. communicated by Mr. Price, in a letter to John Canton, A. M. F. R. S, Philosophical Transactions of the Royal Society of London 53 (1763) 370–418.

[39] P.-R. Wagner, J. Nagel, S. Marelli, B. Sudret, UQLab user manual – Bayesian inversion for model calibration and validation, Technical Report, Chair of Risk, Safety and Uncertainty Quantification, ETH Zurich, Switzerland, 2022. Report UQLab-V2.0-113.

[40] J. Guo, J. Cai, Q. Chen, X. Liu, Y. Wang, Z. Zuo, Dynamic behaviour and energy dissipation of reinforced recycled aggregate concrete beams under impact, Construction and Building Materials (2019).

[41] S. Marelli, B. Sudret, UQLab: A Framework for Uncertainty Quantification in Matlab, 2014, pp. 2554–2563. doi:10.1061/9780784413609.257.

## Appendix A. PCE based meta-model

Consider the articulated structure model represented by a $\mathcal{M}(\mathbf{X})$ as an equivalent mathematical model. Here, $\mathbf{X} \in \mathbb{R}^M$ is a random vector with independent components described by the joint probability density function (PDF) $f_{\mathbf{X}}$. Consider also a finite variance computational model as a map $Y = \mathcal{M}(\mathbf{X})$, with $Y \in \mathbb{R}$ such that:

$$ \mathbb{E}[Y^2] = \int_{\mathcal{D}_{\mathbf{X}}} \mathcal{M}(\mathbf{x})^2 f_{\mathbf{X}}(\mathbf{x}) d\mathbf{x} < \infty \quad (A.1) $$

Then, under the assumption of Equation A.1, the PCE of $\mathcal{M}(\mathbf{X})$ is defined as:

$$ Y = \mathcal{M}(\mathbf{X}) = \sum_{\alpha \in \mathbb{N}^M} y_{\alpha} \Psi_{\alpha}(\mathbf{X}) \quad (A.2) $$

where, the $\Psi_{\alpha}(\mathbf{X})$ are multivariate polynomials orthonormal with respect to $f_{\mathbf{X}}$, $\alpha \in \mathbb{N}^M$ is a multi-index that identifies the components of the multivariate polynomials $\Psi_{\alpha}$ and the $y_{\alpha} \in \mathbb{R}$ are the corresponding coefficients. In practical applications, the sum in Equation A.2 needs to be truncated to a finite sum by introducing the truncated polynomial chaos expansion:

$$ \mathcal{M}(\mathbf{X}) \approx \mathcal{M}^{\text{PC}}(\mathbf{X}) = \sum_{\alpha \in \mathcal{A}} y_{\alpha} \Psi_{\alpha}(\mathbf{X}) \quad (A.3) $$
where, $\mathcal{A} \subset \mathbb{N}^M$ is the set of selected multi-indices of multivariate polynomials.

In this work, the least-angle regression (LARS) method is used to create the PCE meta-model truncated to the maximum polynomial degree (p) ranging from 1 to 20, and using hyperbolic truncation scheme (q) ranging from 0.5 to 1.

$$ \mathcal{A}^{M,p,q} = \{\alpha \in \mathcal{A}^{M,p} : \|\alpha\|_q \le p\}, \quad \text{where, } \|\alpha\|_q = \left(\sum_{i=1}^{M} \alpha_i^q\right)^{1/q} \qquad (A.4) $$

The accuracy of the constructed PCE is estimated by computing the leave-one-out (LOO) cross-validation error ($\epsilon_{LOO}$). It consists in building N meta-models $\mathcal{M}^{PC\setminus i}$, each one created on a reduced experimental design $\mathcal{X}\setminus\mathbf{x}^{(i)} = \{\mathbf{x}^{(j)}, j = 1, ..., N, j \neq i\}$ and comparing its prediction on the excluded point $\mathbf{x}^{(i)}$ with the real value $y^{(i)} = \mathcal{M}(\mathbf{x}^{(i)})$ [34]. The leave-one-cross-validation error can be written as:

$$ \epsilon_{LOO} = \frac{\sum_{i=1}^{N} \left( \mathcal{M}(\mathbf{x}^{(i)} - \mathcal{M}^{PC\setminus i}(\mathbf{x}^{(i)})) \right)^2}{\sum_{i=1}^{N} \left( \mathcal{M}(\mathbf{x}^{(i)} - \hat{\mu}_Y) \right)^2} \qquad (A.5) $$

where, $\hat{\mu}_Y$ is the mean of the experimental design sample.

## Appendix B. Bayesian inference for inverse analysis

Consider the computational model $\mathcal{M}$ that allows the analyst to predict certain quantities of interest gathered in a vector $\mathbf{y} \in \mathbb{R}^{N_{out}}$ as a function of input parameters $\mathbf{x}$:

$$ \mathcal{M}: \mathbf{x} \in \mathcal{D}_{\mathbf{X}} \subset \mathbb{R}^M \mapsto \mathbf{y} = \mathcal{M}(\mathbf{x}) \in \mathbb{R}^{N_{out}} \qquad (B.1) $$

The Bayesian inference for inverse analysis focuses on identifying the input parameters of a computational model to recover the observations in the collected output data set. It comprises of a computational forward model $\mathcal{M}$, a set of input parameters $\mathbf{x} \in \mathcal{D}_{\mathbf{X}}$ that need to be inferred, and a set of experimental data $\mathcal{Y}$. Here, $\mathcal{Y} \stackrel{\text{def}}{=} \{\mathbf{y}_1, ..., \mathbf{y}_N\}$ is a global data set of N independent measured quantities of interest ($\mathbf{y}_i$).

The forward model $\mathbf{x} \mapsto \mathcal{M}(\mathbf{x})$ is a mathematical representation of the system under consideration. The lack of knowledge on the input parameters is modelled by considering them as a random vector, denoted by $\mathbf{X}$ which is assumed to follow a so-called prior distribution (with support $\mathcal{D}_{\mathbf{X}}$), as presented in Figure 8 in the present work.
$$X \sim \pi(x) \qquad (B.2)$$

The Bayesian statistics combine this prior knowledge of the parameters with the few observed data points to obtain a statistical model called posterior distribution ($\pi(x | y)$) of the input parameters, using Bayes' theorem [38], expressed as:

$$\pi(x | y) = \frac{\pi(y | x)\pi(x)}{\pi(y)} \qquad (B.3)$$

Now, considering the available data set ($\mathcal{Y}$) as independent realizations of $\mathbf{Y} | x \sim \pi(y | x)$, the collected measurements result in the definition of the likelihood function $\mathcal{L}(x; \mathcal{Y})$, which is a function of input parameters $x$:

$$\mathcal{L}: x \mapsto \mathcal{L}(x; \mathcal{Y}) \stackrel{\text{def}}{=} \prod_{i=1}^{N} \pi(y_i | x) \qquad (B.4)$$

This implicitly assumes in/dependence between individual measurements in $\mathcal{Y}$. Intuitively the likelihood function for a given $x$ returns the relative likelihood of observing the data at hand, under the assumption that it follows the prescribed parametric distribution $\pi(y | x)$.

As all models are simplifications of the real world, a discrepancy term ($\epsilon$) is introduced to connect real-world observations ($\mathcal{Y}$) to the predictions by the model. In practice, the discrepancy term represents the effects of the measurement error and model inaccuracy. The discrepancy term introduced here reads:

$$y = \mathcal{M}(x) + \epsilon \qquad (B.5)$$

Here, the $\epsilon$ is assumed as an additive Gaussian discrepancy [39] with a zero mean and given covariance matrix ($\mathbf{\Sigma}$):

$$\epsilon \sim \mathcal{N}(\epsilon | \mathbf{0}, \mathbf{\Sigma}) \qquad (B.6)$$

Taking insights from the discrepancy term definition, a particular measurement point ($y_i \in \mathcal{Y}$, is a realization of the Gaussian distribution with mean value $\mathcal{M}(x)$ and covariance matrix $\mathbf{\Sigma}$. This distribution is named as discrepancy model and is expressed as:

$$\pi(y | x) = \mathcal{N}(y | \mathcal{M}(x), \mathbf{\Sigma}) \qquad (B.7)$$
In application, the discrepancy model defines the connection between the supplied data ($\mathcal{Y}$) and the forward model. In the present work, as the inverse analysis is conducted on the pseudo on-site measurement (as evidence), the discrepancy model with known residual variance is assigned.

Afterwards, the $N$ independent available measurements (in our case, maximum velocity pattern and 3D displacement databases) gathered in the data-set (i.e., $\mathcal{Y} = \{\mathbf{y}_1, ..., \mathbf{y}_N\}$) are used to define the likelihood function as:

$$ 
\begin{aligned}
\mathcal{L}(\mathbf{x}; \mathcal{Y}) &= \prod_{i=1}^{N} \mathcal{N}(\mathbf{y}_i | \mathcal{M}(\mathbf{x}), \mathbf{\Sigma}) \\
&= \prod_{i=1}^{N} \frac{1}{\sqrt{(2\pi)^{N_{\text{out}}} \det(\mathbf{\Sigma})}} \exp \left( -\frac{1}{2} (\mathbf{y}_i - \mathcal{M}(\mathbf{x}))^{\top} \mathbf{\Sigma}^{-1} (\mathbf{y}_i - \mathcal{M}(\mathbf{x})) \right)
\end{aligned}
\quad (B.8) $$

Thereby, following Bayes' theorem, the posterior distribution $\pi(\mathbf{x} | \mathcal{Y})$ of the parameters ($\mathbf{x}$) given the observations in $\mathcal{Y}$ can be written as:

$$ \pi(\mathbf{x} | \mathcal{Y}) = \frac{\mathcal{L}(\mathbf{x}; \mathcal{Y})\pi(\mathbf{x})}{Z} \quad (B.9) $$

Here, $Z$ is a normalizing factor, known as the marginal likelihood or evidence, is added to ensure that this distribution integrates to 1:

$$ Z \stackrel{\text{def}}{=} \int_{\mathcal{D}_{\mathbf{x}}} \mathcal{L}(\mathbf{x}; \mathcal{Y}) \pi(\mathbf{x}) d\mathbf{x} \quad (B.10) $$

The closed-form solutions do not exist in practice, and therefore the posterior distribution is obtained through Markov chain Monte Carlo (MCMC) simulations. In the present work, among many proposed algorithms in [39], the Adaptive Metropolis (AM) algorithm is used and 100 parallel chains with 1000 steps are assigned to the MCMC solver. The start of the covariance adaptation in AM algorithm is assigned at the 100th step (see Wagner et al. [39] for more details).

The output ($\mathbf{y}$) predictive capabilities of the Bayesian inference is assessed through the comparison of prior ($\pi(\mathbf{y})$) and posterior ($\pi(\mathbf{y} | \mathcal{Y})$) output distributions as:

$$ \pi(\mathbf{y}) = \int_{\mathcal{D}_{\mathbf{x}}} \pi(\mathbf{y} | \mathbf{x})\pi(\mathbf{x})d\mathbf{x} \quad (B.11) $$

$$ \pi(\mathbf{y} | \mathcal{Y}) = \int_{\mathcal{D}_{\mathbf{x}}} \pi(\mathbf{y} | \mathbf{x})\pi(\mathbf{x} | \mathcal{Y})d\mathbf{x} \quad (B.12) $$
Lastly, in the present studies, the purpose of the Bayesian inference is to obtain the ‘best set of impact condition parameters (ICPs)’. Given the posterior distribution ($\pi(x | \mathcal{Y})$), we are interested in finding a suitable set among the posterior computed set i.e., $X | \mathcal{Y}$. This is done through a point estimator ($\hat{x}$) computed from:

$$ \pi(\mathbf{y} | \mathcal{Y}) \stackrel{\text{def}}{=} \pi(\mathbf{y} | \hat{\mathbf{x}}) \qquad (\text{B.13}) $$

This point estimator can be a mean or mode (maximum a posteriori ‘MAP’) [39] of the posterior distribution as per the user’s choice.