Seismic fragility assessment of geotechnical seismic isolation (GSI) for bridge configuration

The seismic vulnerability of bridges may be reduced by the application of Geotechnical Seismic Isolation (GSI) below the foundations of the columns and the abutments. However, the role of GSI on the seismic response of bridges has been limitedly examined in literature. Therefore, this research has been conducted to study the effect of applying GSI on the seismic response of bridges to address the aforementioned gap in knowledge. Advanced nonlinear dynamic three-dimensional finite element analyses have been conducted using OpenSees to study the influence of the GSI. The cases of traditional and isolated bridges subjected to earthquakes have been considered to assess the GSI effects. The results showed that the GSI reduces the seismic effect on the column while its effect seems to be less significant for the abutments. In addition, fragility curves for the traditional and isolated cases have been developed and compared to provide insights with a probabilistic-based approach. The results of this paper provide a useful benchmark for design considerations regarding the use of GSI for bridges.


Background
Geotechnical Seismic Isolation (GSI) may be considered as an innovative method to reduce the seismic vulnerability of civil structures and infrastructures. The use of the GSI has been examined by several researchers (e.g., Tsang 2009;Tsang et al. 2012). Tsang (2009) introduced the definition of the GSI itself. In addition, Tsang et al. (2012) examined the use of rubber-soil mixtures on the acceleration and inter-story drift of low to medium rise 1 3 buildings. Equivalent liner model was used to model the dynamic response of the soil and the GSI in Tsang et al. (2012) study, where the dynamic response of the soil was characterized using the secant shear modulus and the damping ratio. Forcellini (2017) investigated the effect of the thickness of the GSI on the response of a bridge using different earthquake records with a PGA range of 0.3-0.9 g. Forcellini (2017) used granulated rubber-soil mixtures as the GSI material.
More recently, Forcellini (2020) examined the efficiency of the GSI for several configurations of buildings using three-dimensional finite element analysis. The paper focused specifically on the effect of the flexibility of the building, which was simulated by changing the number of stories of the building. Thus, buildings with a number of stories of 3, 5, and 7 was modelled in this study. The study results were presented in terms of the acceleration change due to the use of the GSI. Tsang et al. (2021) performed centrifuge studies to examine the efficiency of the GSI system for the case of a building subjected to seismic effect. The case of a building subjected to earthquake effect and resting on GSI with a thickness of 2.0 m was modelled in the centrifuge experiments. Rubber-soil mixture was used as the GSI system with two different rubber percentages (30% and 40%). The study analyzed the effect of the GSI on the displacement of the building, deformation of the columns of the building, and the acceleration of the building. Banović et al. (2019) conducted series of small-scale laboratory tests to examine the response of a building resting on a stone pebble as a GSI system. They studied the effect of the thickness of the stone pebble layer, fraction of pebbles, pebble compaction and moisture, pressure applied by the foundation, and repeated shakes on the response of the building. Banović et al. (2020) examined the efficiency of the stone pebble as a GSI system for small scale buildings with foundation width of 0.7 and 1.2 m, respectively, using small scale laboratory models. In addition, the influence of the stiffness of the building combined with the GSI system is also examined. They noticed that increasing the size of the foundation for this GSI system reduced the effect of rocking and the bearing capacity of the foundation. Tsiavos et al. (2019) examined the effect of the size of the rubber on the engineering properties of the sand-rubber mixture. The authors then used the sand-rubber mixture in a 1-g shaking table test to the study the influence of the size of the rubber and the thickness of the sand-rubber mixture layer (GSI) on the performance of a rigid block made from steel and subjected to seismic shake. They also investigated the effect of the foundation material (i.e., bottom of the foundation) which is in contact with the GSI layer. Furthermore, Tsiavos et al. (2020) proposed a PVC 'sand-wich' (PVC-s) seismic isolation system, which consisted of a thin film of sand grains between two PVC sheets (sandwiching). The authors performed large-scale shaking table tests to demonstrate the effectiveness of the low frictional resistance of sand particles against polymers. In other study, Tsiavos et al. (2021) considered two other seismic protection mechanisms; these mechanisms involved the use a steel wire mesh and steel ties attached to the surface of the outer walls of the structure to improve seismic protection against in-plane and out-of-plane failure. These mechanisms allowed the development of seismic mitigation strategies in developing countries using robust and low-cost design solutions. Dhanya et al. (2020) proposed geogrid reinforcements to improve the low bearing capacity of the GSI layer for the case of a low rise building subjected to seismic effect. Two-dimensional finite element analysis was conducted to address the aim of the study. The study demonstrated the efficiency of using a GSI layer with a thickness of 0.1-0.2 of the width of the building with two layers of geogrid reinforcement to overcome the low bearing capacity and high compressibility of the GSI layer. Forcellini (2021a) developed fragility curves for the case of a building resting on layers of soil and GSI and subjected to seismic shake. Pitilakis et al. (2021) examined the effect of the percentage of rubber of a gravel-rubber mixture on the dynamic response of a vibrating steel frame using field-based tests. The layer of the rubber-gravel mixture was placed directly underneath the foundation of the steel frame. The thickness of gravel-rubber layer was considered equal to 0.5 m with rubber percentages of 10 and 30%, respectively. The authors found that a percentage of 30% of rubber was enough to cut off all of the waves transmitted from the steel frame to the soil due to vibration effect. Gatto et al. (2021) studied the effectiveness of a proposed GSI system for the case of a building subjected to seismic effect. The proposed GSI system involved the use of a Polyurethane injected in the soil. A building resting on sand and subjected to seismic effect was considered as the case study to examine the effectiveness of the proposed GSI system. The research concluded that the polyurethane reduced the acceleration induced in the building due to seismic effect. Table 1 shows the existing literature in order to make a comparison on the bases of the type of the study, the materials that have been used, the constitutive models used for both the soil and GSI, the investigated parameters, and considered intensities of the performed earthquakes.
More recently, Kuvat and Sadoglu (2020) investigated the dynamic properties of asphalt-sand mixtures as a GSI material using cyclic triaxial tests. In addition, Sarajpoor et al. (2021) performed a series of large size dynamic hollow cylinder tests to assess the possibility of increasing damping of sandy soils by adding bitumen. The effects of different factors such as bitumen content, loading frequency, and bitumen viscosity on the dynamic properties of the mixture were examined in this research. Furthermore, several contributions investigated the effect of soil-structure-interaction on the seismic vulnerability of buildings, such as Petridis and Pitilakis (2021) and Forcellini (2021b, c).
As illustrated in the aforementioned detailed literature review, there has been limited studies on the application of the GSI in bridges. This is due to the relatively low bearing capacity and high compressibility of the GSI system. However, this problem could be solved by placing the GSI layer in depth wherever the loads imposed by the bridge distributed so that it becomes very low and do not impose stability or settlement problems; this approach has been considered in this paper. Hence, this paper aims to assess the beneficial effect of GSI system for the case of a bridge subjected to seismic effect using fragility curves. Two models have been developed to address the aim of the paper; the first model is a benchmark structure representing a typical Californian highway bridge founded on a rigid soil. The second model represents a realistic case of the same configuration where GSI has been applied at the foundation soil. The comparisons are conducted employing the probability of exceedance for different intensities.

FEM methodology
This paper aims to propose analytical fragility curves for a benchmark bridge configuration that represents a typical Californian highway bridge. Three-dimensional finite element analysis (FEM) has been employed to accurately model the complicated interaction among the soil, the foundation, and the structure of the bridge. The use of the three-dimensional FEM has been deemed essential as it was aimed to take into account the effect of the rotations along the vertical and transversal axes which is not possible to be simulated in the two-dimensional FEM. Nonetheless, the three-dimensional modelling of this problem is quite challenging as it involves modelling the complex response of the soil when subjected to seismic effect which requires accurate simulation of the soil non-linearity, soil plasticity, and hysteretic damping. In addition, the response of N/A means that the information is not applicable N/P means that the information has not been provided in the paper Table 1 (continued) the GSI should also be modelled accurately to ensure robust results. Thus, the wellknown open access software OpenSees (Mazzoni et al. 2009) has been employed in the numerical modelling of the soil because this software contains soil constitutive model (pressure-independent multi-yield (PIMY) model) which has the capability to accurately model the dynamic response of the soil (Yang et al. 2003;Alzabeebee and Forcellini 2021).
The developed finite element model is shown in Fig. 1a, b. Furthermore, the details of the bridge modelling are presented in Figs. 2 and 3. The following subsections discuss the construction of the finite element models.  based on the wavelength of the seismic shake and the ultimate frequency (i.e., the frequency beyond which there was an insignificant spectral content). As suggested in Forcellini (2017Forcellini ( , 2018Forcellini ( and 2019, the size of the elements has been gradually raised as the soil domain distances from the bridge (center of the model) to the model boundaries to ensure the consideration of the pure shear which develops in locations away from the bridge. In addition, the accuracy of the developed finite element mesh has been further examined by ensuring that the acceleration of the surface in the locations near the model boundaries is identical to those of free field condition. The boundary conditions have been modeled as transmitting boundaries in order to dissipate the radiating waves and to accurately consider the damping by preventing the reflection of the seismic waves back into the soil medium. Shear deformations have been allowed by leaving the longitudinal and transverse directions unconstrained at lateral boundaries which have been modelled by adopting the Penalty method. A tolerance of 10 −4 has been selected in order to avoid problems associated with the equations system conditions (see . Additionally, the base of the finite element model has been restrained in all of the directions to simulate the rock layer below the soil deposits as suggested in previous studies (Alzabeebee 2019a(Alzabeebee , b, 2020Kampas et al. 2020). The developed finite element model is shown in Fig. 1a, b. Furthermore, the seismic effect has been modelled as a prescribed acceleration applied at the bottom of the finite element model in a methodology similar to that employed in previous studies (Alzabeebee 2019a(Alzabeebee , b, 2020Kampas et al. 2020).

Soil modelling
The soil has been assumed to be homogenous over the whole depth to simplify the analysis. Furthermore, a stiff soil has been modelled in the analysis to simulate the case of bridge sit on a stiff soil deposit or a soil layer which has been improved by one the soil improvement methods (e.g., compaction). The Pressure Independ Multi Yield (PIMY) model has been utilized to model the response of the soil in the analysis. This constitutive model has been considered a robust option because it is able to simulate the key features of the soil response when subjected to seismic load (non-linear response to load, hysteretic damping, and permanent deformation) and has been used in many previous studies to simulate the seismic response of the soil (Su et al. 2017;Mina and Forcellini 2020). The PIMY model also utilizes the Von Mises failure criteria (Yang et al. 2003). The parameters of the soil employed in the analysis are taken from Forcellini (2017) and are shown in Table 2.

Foundations
The three-dimensional modelling has been conducted to robustly model the soil-structure interaction. This model allows to account for the seismic settlement, title, and displacement of the foundation. Some assumptions were introduced to reduce the computational demand of the soil-structure interaction of this problem; these assumptions are as follows: (1) The foundations of the abutments have been simulated as rigid slabs using an equivalent linear elastic model. (2) The foundation of the column was designed as a Type 1 Caltrans pile shaft, with a diameter similar to that of the column on top of it (Forcellini 2017). Thus, the longitudinal reinforcements extend continuously from the pile shaft into the column. In addition, the thickness of the foundation is considered equal to 0.5 m for the abutments; this foundation has been designed using the most detrimental loading condition (i.e., highest applied pressure and the highest possible bending moment). In addition, the interaction of foundations of the abutments and the soil has been modelled using rigid link elements similar to the approach adopted in Forcellini (2017Forcellini ( , 2021a. The connection between the points of the soil and the rigid links have been modelled using two points with equal degrees of freedom (Mazzoni 2009) to ensure similar displacement due to static and dynamic loads as shown in Fig. 3. The pile shaft was assumed to have a diameter of 1.22 m and a length of 10 m and was linked with the supported column using equal degree of freedom (Mazzoni et al. 2009).

Traditional bridge model (fixed based model)
The reference case with no GSI has been named as fixed based (FB) model (Fig. 2). The followings summarize the assumptions made in the FB model.
(1) A linear elastic model has been employed to model the concrete of the column and the deck slab. The modulus of elasticity, Poisson's ratio and the unit weight of the concrete are considered equal to 2.8 × 10 4 MPa, 0.20, and 24 kN/m 3 , respectively. This simplification is based on the assumption that the designed GSI and the isolators on the abutments perform correctly and thus the bridge elements (column and deck) may be assumed to be capacity designed so that it is able to respond in the elastic range.
The cross-sectional area, and moment of inertial in the transverse and longitudinal directions are considered equal to 0.950 m 2 , 0.108 m 4 and 0.108 m 4 , respectively. The connection between the deck slab and the column has been modeled with a fixed link.
(3) The deck (90 m long) has been assumed to behave following the linear elastic model since the connections at the ends consist of isolation devices that guarantee capacitydesigned conditions. Therefore, the deck has been modelled using beam elements with linear elastic properties. The cross-sectional area, and moment of inertial in the transverse and longitudinal directions are considered equal to 5.72 m 2 , 2.81 m 4 and 53.90 m 4 , respectively. The dead and live loads above the deck slab have been modelled using a distributed load with a magnitude of 130.3 kN/m. (4) The deck slab has been connected with the abutments using isolation devices in the longitudinal direction as shown in Figs. 2 and 3 In particular, the deck has been linked with a rigid link that connects two nodes to which the equivalent linear springs (modeling the isolators) have been connected to (detail in Fig. 2c). The isolators consist of seismic isolators that simulate ISOSISM HDRB-S (http:// www. fpcit alia. it/ freys sinet) which were designed following the design procedure described in the Eurocode 8, parts 1 and 2 (UNI EN-1998EN- -1:2004UNI EN-1998-2:2005) and consisted of: (1) Determination of the equivalent elastic stiffness of the device by assuming a target value for the period of the isolated bridge (T BI = 0.8 s): where, M is the mass of the system and TBI is the fundamental period of the isolated system. (2) The vertical loads under static and seismic conditions were calculated in order to select the diameter, shear modulus, damping, rubber width, and total height of the isolators (equal to 650 mm, 0.4 MPa, 10%, 105 mm, and 206 mm, respectively). They were modeled with the model named M1 in Forcellini (2021b) to simulate the isolation system with an equivalent linear formulation that considered equivalent linear springs of 1260 kN/m and 10% damping (Fig. 2c) along longitudinal direction.
(3) The design displacement was then considered to verify the safety of the response.
The number of designed devices (2 for each abutment) were determined in order to minimize the eccentricity between the center of stiffness of the isolation system and the center of mass of the superstructure.
The other directions have been modeled with high stiffness zero-length elements to eliminate the deformations. The height of the embankment of the abutments have been considered equal to 25 m and have been modelled as a dead load applied on the abutment foundation. The applied load is estimated to be equal to 30,000 kN. The fundamental period of the bridge in the longitudinal direction is 0.811 s.

The GSI model
The GSI model is structurally similar to the FB model with the only difference is that the GSI has been simulated at a depth of 10 m below the base of the foundation slabs. This depth has been chosen as a compromise between the economic costs of the application and the efficiency of the GSI (more details in Forcellini 2021b). The location of the GSI is shown in Fig. 4. In particular, the first 10 m of soil was verified to distribute the vertical loads and thus reducing the contact pressures on the foundation level due to the vertical loads. In this regard, the first 10 m of the soil below the foundation slabs needs to be sufficiently stiff with high shear strength parameters to guarantee that bearing capacity and compressibility will not be a potential limitation to the use of the GSI in bridge applications. In addition, the limited thickness of the GSI (i.e., 0.5 m) also helped to reduce the potential effect of the high compressibility and low shear strength parameters of the GSI layer. The GSI has been simulated using the PIMY model and the parameters (Table 2) have been calibrated by considering the values proposed by Anastasiadis et al. (2012) for a GSI layer consisting of a Rubber-soil mixture material (RSM) with a rubber percentage of 30%. In particular, the value of G (5.0 MPa) has been demonstrated to be appropriate by Tsang et al. (2021).

Seismic scenario
The seismic scenario is described by 100 input motions, selected from the PEER NGA database (http:// peer. berke ley. edu/ nga/), to significantly the effect of the dynamic characteristics of the system. These input motions have been applied at the base of the models, along the x-axis (longitudinal direction). The motions have been divided into 5 series of 20 motions each with characteristics of moment magnitude (Mw) and closest distance (R), as already applied in Forcellini (2017): (1) Mw = 6.5-7.2, R = 15-30 km, (2) Mw = 6.5-7.2, R = 30-60 km, (3) Mw = 5.8-6.5, R = 15-30 km, (4) Mw = 5.8-6.5, R 30-60 km, (5) Mw = 5.8-7.2, R = 0-15 km. of the selected records. These measures (PGA and SA) will enable having multiple levels of intensities, which is essential to obtain distributed outputs and develop fragility curves. It should be mentioned that the fragility curves are used to consider a probabilistic approach that implicitely considers the effect of the repeated seismic motions without considering the uncertainties connected with materials (such as deterioration, age, fatigue,... etc).

Methodology of fragility curves
Fragility curves consist of a probabilistic-based approach to express the relationship between the demand and the capacity. In this paper, analytical fragility curves have been built to represent the relationship between the likelihood of exceeding threshold limits of damage by FEMA (2003) and applied in Xie and DesRoches (2019). The paper refers to two EDPs (Engineering Demand Parameters): 1) the drift of the column and 2) the deck displacement. Both of these EDPs were herein considered linear because of the assumption that both of the isolators at the abutments and the GSI allow to consider the bridge elements (column and deck) to be capacity designed so that they respond in the elastic range. Several limit states have been considered for the drift of the column: 0.80, 1.51, 3.00 and 4.00 for slight (SL1), moderate (SL2), extensive (SL3), and complete (SL4) damage state, respectively. Two limit states have been considered for the longitudinal displacement between the deck and the abutments: 102 mm and 305 mm as slight (SL1) and moderate (SL2) damage states, respectively (more details in Xie and DesRoches 2019). The results from the analyses have been obtained for the two models (FB and GSI models) and for the 100 selected input motions by considering a representative intensity measure (Im). The assumptions herein are that the uncertainties of the results may be represented by lognormal distributions and two parameters were calculated: the logarithmic mean (µ) and standard deviation (β) of the lognormal seismic intensity measure. Therefore, the results were used to develop linear regressions and also to determine the values of the mean and the log-standard deviation. The probability of exceedance (PE) was then calculated as: where PE is the probability of the structural damage (D) to exceed the i-th damage state (C), while ϕ is the standard normal cumulative distribution function (more details in Forcellini (2021a, c). Figure 6a, b present the effect of increasing the seismic intensity on the response of the bridge for both of the considered configuration (with and without GSI). PGA (peak ground acceleration) and SA (spectral acceleration) have been considered as intensity measures for the column drifts and the displacements of the deck, respectively. The aforementioned limit states are also included in the figures to provide an insight into the demand and the capacity. It is obvious based on the figures that increasing the seismic

Seismic effects on the response of the bridge
intensity rises the column drift and the deck displacement due to the development of strain in the column and also in the deck-abutment links. However, it is evident from the figures that the GSI reduces the effect of the seismic shake compared with the FB. The comparison among the two intensity measures is fundamental, in order to consider the best fitting variable in terms of quality of interpolation. In this regard, Table 3 shows the comparison between the values of the coefficient of determination (R 2 ) that represents the proportion of variance in the drift/deck displacement that can be explained by the various intensity measures, or the quality of the regression to fit the obtained results for various models (FB and GSI). The table shows that SA best fits the relationship between the seismic intensity and the drift/deck displacement and thus, SA will be applied in the next sections to develop fragility curves.

Fragility curves for the FB model
This section presents the results for the FB model in terms of the fragility curves for both the column and the deck displacement, as defined in the previous section. In particular, Tables 4 and 5 show the values of the logarithmic mean (µ) and standard deviation (β) calculated from the results shown in Fig. 6a, b. In particular, the slope of the curves depends on the lognormal standard deviation (β) that defines the dispersion about the mean value. Figure 7a, b show the fragility curves for the column and the abutment, respectively. For both of the structural elements, it is worth noticing that the probabilities of exceedance (PE) connected with the various limit states reach high values even for small intensities and this occurs for all the limit states. If the column ductility and the deck displacements are compared for the correspondent limit states, it can be said that the vulnerability of the columns is higher than that of the abutments. In this regard, for the same level of spectral acceleration (i.e., 0.10 g), PE reaches 0.726 and 0.663 for LS1 for the column drift and the deck displacement. Therefore, the contribution of the abutments to the bridge system fragility is less significant if compared with the column fragility, confirming what demonstrated in Xie and DesRoches (2019. In other words, the column is shown to drive the fragility of the entire system since isolation devices are concentrated on the connections between the deck and the abutments while the column-deck connection does not allow dissipations.

Fragility curves for the GSI model
This section shows the results for the GSI model in terms of the column ductility and deck displacement. In particular, Tables 6 and 7 show the values of the logarithmic mean (µ) and standard deviation (β) derived using the results presented in Fig. 6a, b. Figure 8a, b present the obtained fragility curves for the column and abutment, respectively.
For the column (Fig. 8a), it is worth noticing that the probabilities connected with various limit states are smaller than those resulted for the FB model. In particular, for PGA of 0.40 g, PE are equal to 0.821, 0.773, 0.712 and 0.637 for LS1, LS2, LS3 and LS4, respectively, demonstrating the role of the GSI in reducing the vulnerability of the system. For the abutment (Fig. 8b), the PE are more similar to those resulted for the FB model. For example, at PGA of 0.10 g, the probabilities of exceedance are equal to 0.687 and 0.410 for LS1 and LS2, respectively. This demonstrates that in the case of the GSI, the abutments are  significantly more vulnerable than the column. In particular, for the same level of spectral acceleration (i.e., 0.10 g), PE reaches 0.136 and 0.587 for LS1 for the column drift and the deck displacement, respectively. These results demonstrate the importance of the GSI in protecting the column and hence, in this case, the abutments drive the fragility of the entire system (contrarily to model FB, Sect. 4.2).

Comparison of the GSI and FB fragility curves
The fragility curves correspondent with the FB and GSI are compared for the column for different LSs in Fig. 9. The differences among the PE resulted from the two models depend on the considered limit state. For example, for PGA of 0.30 g, the PE for the GSI model is 23.03%, 28.11%, 31.07% and 36.07% smaller than that calculated for the FB model (more details in Table 8), meaning that 1) GSI has beneficial effect in reducing the vulnerability of the system for all the level of damage, 2) GSI effectiveness increases for high level of damage. When the fragility curves for the abutments are compared (Fig. 10), the beneficial effect of the GSI is less pronounced. For example, for a PGA of 0.20 g, for GSI model there are reductions of 6.89% and 7.87% than the values calculated for the FB model (more details in Table 9). It is worth noting that at the abutments GSI is coupled with base isolation. In particular, following CALTRANS specifications (criterion 6.3.1.3-1), the seismic response of the bridge may be dominated by the abutments when the effective abutment longitudinal displacement is large (for high intensities). Therefore, for our case study, the role of the isolation on the abutments may become predominant on the GSI and thus the presence of the isolators improves the seismic behaviour of the bridge mainly at high intensities. In addition, it is worth noting that GSI is significantly effective when applied to the columns, confirming that the deck displacement adds a minor increase to the fragility of the bridge system if compared with the column fragility; similar observations are also noted by Xie andDesRoches (2019 andForcellini (2021a). The column is shown to drive the fragility of the entire system and thus, the application of the GSI below the column contributes to reduce the seismic vulnerability of the system.
Overall, the presented results are limited to the application that is based on the interposition of the GSI layer between two strata (Fig. 4). It is also worth noting the importance of the first 10 m layer in order to reduce the pressures applied by the bridge and thus    guarantee that the bridge configuration will not encounter stability or settlement problems due to low bearing capacity and high compressibility of the GSI layer.

Conclusions
This paper investigated the influence of the GSI on the seismic vulnerability of a benchmark bridge configuration representing a typical Californian highway two-span concrete bridge. Three-dimensional soil-foundation-structure numerical models were performed with OpenSees to derive analytical fragility curves that describe the vulnerability of the columns and the abutments. Several limit-states were considered according to the multihazard loss estimation methodology proposed by HAZUS (FEMA 2003). The results showed the beneficial effect of applying the GSI in reducing the seismic vulnerability of the bridge, especially at lower intensities particularly below the column. The application of the GSI at the abutments seems not to be efficient due to the presence of the isolators in the connection between the abutments and the deck which helps in dissipating the seismic energy. Therefore, applying the GSI only below the column is an important outcome in order to reduce the costs of massive design. Even if the present study is conducted with one specific case study, the general conclusions are applicable to other similar bridge configurations. In this regard, the findings may be applied by several stakeholders, such as designers, experts, and bridge owners to refine the isolation modelling and develop reliable fragility curves for other classes of bridges.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions. The authors have not disclosed any funding.

Conflict of interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.

Availability of data and material
The data used to support the findings of this study are available from the corresponding author upon request.

Code availability
The data used to support the findings of this study are available from the corresponding author upon request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.