A 3D Numerical Approach to Assess the Temporal Evolution of Settlement Damage to Buildings on Cavities Subject to Weathering

The goal of this paper is to show how a recently developed advanced chemo-hydro-mechanical (CHM)-coupled constitutive and numerical model for soft rocks can be applied to predict the temporal evolution of settlement damage to buildings on cavities subject to weathering. In particular, a Building Damage Index (BDI) and its evolution with time is proposed. The definition of the BDI is inspired by the work of Boscardin and Cording (J Geotech Eng 115:1–21, 1989) and uses the surface differential settlements obtained by finite element (FE) analyses to assess how far a building is from a non-acceptable service condition. By modelling the reactive transport of chemical species in 3D and using a coupled CHM constitutive and numerical model, it is possible to simulate weathering scenarios and monitor the temporal evolution of surface settlements making the BDI time dependent. This approach is applied to evaluate the damage evolution of two buildings lying on two anthropic caves in a calcarenite deposit belonging to the Calcarenite di Gravina Formation. Standard and advanced experimental tests are performed on the in situ material, and the results are used to calibrate the constitutive model. The soundness of both constitutive relationship and reactive transport solver is subsequently tested by simulating two laboratory-scale boundary value experiments. The first is a model footing test on dry and wet calcarenite, while the second is a small-scale pillar that, after the saturation-induced short-term water weakening, fails due to a long-term dissolution weathering process. Finally, both 2D and 3D coupled FE analyses simulating different weathering scenarios and corresponding settlements affecting the buildings above the considered cavities are presented. Particular attention is placed on assessing the BDI and its temporal evolution.


Introduction
Sinkholes, of both natural and anthropogenic origin, represent an important geological hazard in Europe and, worldwide and every year, cause significant property losses and causalities. Sinkhole formation is often the end result of complex physical hydraulic chemical phenomena (weathering) inducing continuous changes of the soils and rocks forming the earth crust causing a progressive reduction in mechanical properties including strength and deformability and variations of permeability. From an engineering perspective, the weathering-induced degradation of mechanical properties makes such phenomenon quite relevant as it generates deformation in the considered geotechnical systems and can compromise its mechanical stability in the worst case. For instance, leaching chemical-aggressive species can degrade the bond strength of cemented soils, causing settlements of structures founded on these materials (Lollino et al. 2013). In evaporites, the interaction of water with the gypsum pillars of an abandoned mine can transform this rock into a much weaker and compressible material, which can eventually collapse under the weight of the overburden (Castellanza et al. 2010). Also the simple saturation of a dry porous carbonate rock, limestone or other karstic formations leads to considerable changes in strength and stiffness that may cause large surface settlements unacceptable for the service life of existing structures, or even collapse in the form of sinkholes (Parise and Lollino 2011;Parise and Vennari 2017;Sowers 1996). The substantial increase in sinkhole formation in the Apulian region in the last decades (Parise and Lollino 2011;Parise and Vennari 2017;Sowers 1996) has pushed the regional authorities to introduce new regulations regarding the procedures and requirements needed for risk assessment (L. R. 2002/19). In particular, in the annex introduced in 2006 specifically regarding the risk associated with underground mines (A.d.b. Puglia 2006), it is specified that for all existing structures in high-risk regions (PG3) two options are possible: the first is to fill the cavity with materials with similar mechanical properties of the host Calcarenite, the second is to perform numerical analyses to assess and predict the mechanical behaviour of the geosystem and propose eventual consolidation measures to fulfil required levels of serviceability and safety. Also a simple planning permission to change the use class of a construction now requires the application of one of these two options. When a building is resting on the top of an abandoned cavity (Fig. 1), the engineering question to be addressed is: how long will the building remain in serviceability conditions?
Despite considerable progress made in recent years in the study of the interactions between the soil (or rock) mass and the environment (including hydraulic, thermal and chemical coupling phenomena) (Alonso et al. 1990;Cekerevac and Laloui 2004;Della Vecchia and Romero 2013;Gajo et al. 2015;Gens 2010;Hueckel and Borsetto 1990;Jommi 2000;Uchaipichat and Khalili 2009), the application of advanced models to perform 3D numerical simulations to address such question is rare. In particular, chemo-hydro-mechanical (CHM)-coupled advanced numerical simulations are usually performed for 2D idealized problems to show the robustness of integration schemes or to simulate laboratory-scale boundary value problems (Fernandez-Merodo et al. 2007;Tamagnini et al. 2002;Tamagnini and Ciantia 2016).
In engineering practice, according to Eurocode 7 (EC7-BSI 2004), it is necessary to demonstrate that building serviceability requirements can be met at the serviceability limit state (SLS). Typically, building settlements (either gross or distortional) have to be smaller than limiting values which depend on the building type, ground type, and design approach. The choice of the analytical or numerical method used to assess the expected values of settlements is left to Fig. 1 a, b Map of Canosa di Puglia with the highlighted high-, medium-and low-risk zones and case study location (http://93.51.158.165/gis/ map_defau lt.phtml ), c section view of cave interacting with building; d picture of Cave A the engineer, and to this end, elastic methods, Winklerbased analyses and standard finite element (FE) analyses may typically be used (Knappett and Craig 2012). However, the effect of weathering on serviceability of buildings is not contemplated by the design code, and in particular, no specific requirements for considering this are given in EC7. As shown by Conte et al. (2010), a simple way to model weathering in FE analyses is to assume that a portion of the domain consists of a formation with lower mechanical characteristics. As described therein, the authors show how weathering changes stability of a weathered slope modelled in 2D. However, limit state analyses do not provide reliable deformation results as the deformations associated with c-ϕ analysis are not representative of reality.
This work proposes to exploit a CHM-coupled numerical model to study the serviceability conditions of buildings lying on cavities affected by environmental weathering by introducing a practical scalar index: the building damage index (BDI). To answer to the question 'how long will the building remain in safety conditions', instead of assuming a strength profile and then performing a strength reduction analysis, hydro-chemical weathering scenarios are postulated and, by means of the reactive transport of contaminants, they are simulated. In this way, the spatial-temporal evolution of rock strength and stiffness will develop, and the surface settlements will induce damage to the structure until it is no longer safe. The method hence requires the use of both a coupled CHM constitutive model and a FE code which solves the hydro-mechanical problem and the transport of chemical species in porous medium. Small-scale coupled weathering experiments (needed to calibrate and validate the model) are then a prerequisite of fundamental importance for the proposed approach. In fact, in order to make reliable future predictions of the BDI, the model must first capture the temporal evolution of weathering at the laboratory scale.
In the following, after describing the method used to calculate the temporal evolution of the BDI, a brief overview of 3D contaminant transport modelling concepts is presented. Subsequently, the main elements of a recently developed advanced multiscale coupled CHM constitutive model (Ciantia and di Prisco 2015) used to describe the degradation of the soft carbonate rock induced by changes in degree of saturation (short-term debonding) and chemical dissolution (long-term debonding) (Ciantia et al. 2015a, b;Ciantia and Hueckel 2013) are described. Standard mechanical and advanced CHM experimental tests on the calcarenite from site are also performed and used to calibrate the model. Experimental results of two boundary value experiments reproducing weathering scenarios at the laboratory scale are also presented and used for validation. Finally, 2D and 3D coupled FE analyses simulating long-term effect of the weathering-induced settlements and the temporal evolution of BDI of a building in Canosa di Puglia built on top of an abandoned cavity (Fig. 1) are presented. Attention is placed on comparing how the BDI evolves with time for the 2 and 3D cases.

A Time-Evolving Building Damage Index, BDI(t)
Despite the fact that final stages of sinkhole formation in soft rocks is usually characterized by a rapid brittle type of global failure (Lollino and Andriani 2017), recent interferometric synthetic aperture radar (inSAR) measurements show that ground movements start to accumulate with time several months before the collapse (Chang and Hanssen 2014;Intrieri et al. 2015;Jones and Blom 2014;Nof et al. 2013). In the Marsala Panatletico 2011 sinkhole, signs of instability, such as open fractures and detachment blocks, were observed already 1 year prior to collapse (Fazio et al. 2017). Leucci et al. (2002) show that due to leakage of water into the ground can cause differential surface settlements causing the opening of fractures on the overlying buildings. The formation of cracks on buildings, implying ground level subsidence, months before the Gallipoli sinkhole of 2007 ( Fig. 2) were also observed (Delle Rose 2007). These surface settlements, if identified within time, can avoid disasters and allow for actions for risk mitigation (De Giorgi and Leucci 2014). All these observations lead to the need of introducing an approach able to estimate the level of safety of the buildings lying on rocks susceptible to weathering prior to the brittle  (Polimeno 2007) collapse of the geosystem, namely the building damage index, BDI. The formulation of the BDI is inspired by the Eurocode 7 serviceability limit states (which are defined as states that correspond to conditions beyond which specified service requirements for a structure or structural member are no longer met) and the work of Boscardin and Cording (1989) that related the level of damage of buildings to surface settlements. Accordingly, the surface differential settlements obtained by FE analyses are here used to assess how far a building is from a non-acceptable service condition. By modelling the reactive transport of chemical species and using a coupled CHM constitutive model where physical time is an intrinsic variable and rock damage is accounted for, it is possible to simulate weathering scenarios and monitor the temporal evolution of surface settlements and hence a BDI function of time. The BDI(t) is clearly strongly dependent on the weathering scenario assumed and on the constitutive model used to reproduce the mechanical behaviour. It gains in reliability as both the weathering mechanism reproduced, and the coupled constitutive models become more realistic. While advanced constitute relationships for calcarenites are now well established (Ciantia and di Prisco 2015), realistic hydrochemical weathering scenarios for underground cavities such as the ones reported by Ghabezloo and Pouya (2006) are harder to validate due to complexity of the processes involved. Nevertheless, in the absence of in situ experimental data on the hydro-chemical initial and boundary conditions, one may assume worst-case scenarios, that would in any case result more realistic than a classical stability analysis trough simultaneous homogenous strength reduction in all the whole domains (Castellanza et al. 2008). Temporal evolution of surface settlements are then used to determine the corresponding evolution of horizontal strain and angular distortion at the base of the building used to employ the well-known Boscardin and Cording abacus and to calculate BDI(t) that, referring to Fig. 3, is defined as: where ||OA(t)|| is the distance of the current state from the origin, while ||OB|| is the distance from level of damage 3 to the origin obtained as an extension of OA. When the BDI < 1, the settlement-induced damage on the building is slight or negligible, while if the BDI > 1, the damage is moderate or severe.
In order to employ Boscardin and Cording's abacus to define the subsidence induced damage to the building and use Eq. (1), the following assumptions are required: (i) the weathering effect causes surface settlements (subsidence) like the ones induced by tunnelling excavation, (ii) the building stiffness is not considered as it is assumed not to affect the subsidence profiles, and the building displacement and deformations are assumed to coincide with the ones of the rock and (iii) the building is assumed to comply with the Boscardin and Cording (1989) mechanical and geometrical assumptions.
The first assumption is deemed to be reasonable as the bell-shaped subsidence basin that characterizes sinkholes (Gutiérrez et al. 2008) has 2D sections which are very similar to the typical Gaussian-shaped transverse settlement observed in tunnelling excavation (Peck 1969). The latter two also limit the use of this approach to a certain typology of buildings, the one used to construct the abacus.
Despite the fact that this work uses a simplified method to define the level of structural damage, the methodological approach would remain unchanged also when running full soil-structure interaction simulations and can also be applied using other abaci or building damage criteria (e.g. Burland 1995;Burland and Wroth 1974;Franzius et al. 2006;Potts and Addenbrooke 1997;Son and Cording 2005).
The key aspect needed to introduce the BDI(t) is to link the evolution of surface settlements to the hydro-chemical time-dependent weathering processes, and this is only possible if the mechanical behaviour is described by means of advanced coupled CHM constitutive and numerical models.

Constitutive and Numerical Modelling of Weathering in Soft Porous Carbonate Rocks
The FEM code used in this work is the GeHoMadrid code (GHM) developed by prof. M. Pastor and co-workers (Fernandez-Merodo et al. 2007). GHM is a non-commercial research code which solves the hydro-mechanical problem, following the u − p w Biot coupled formulation (Zienkiewicz  Boscardin andCording (1989) plot andShiomi 1984), and hence the hydraulic problem (new seepage velocity, pore-pressure and saturation variations) and the mechanical response are solved simultaneously.

Modelling the Transport of Chemical Species in Porous Medium
From a mathematical point of view, the transport of a reactive compound is derived by fulfilling a mass balance equation for an elementary volume of the porous medium (see Zheng and Bennett 2002). The effects of chemical reactions on solute transport are generally included in the advection-diffusion equation through a chemical sink/source term, which can be formulated for each chemical compound of interest and added to the general advective-dispersive partial differential equation as follows: where C k is the concentration of the kth compound, n is the effective porosity, q i is the component of the real seepage velocity along the i axis, D ij is a diffusion tensor, incorporating the effects of dispersion, q S is the volumetric flow rate of sink/source per unit volume and C k S is the concentration associated with the source or sink and ∑ N n=1 R n is the chemical sink/source term representing the rate of change in solute mass of a particular species due to N chemical reactions. In this work, as in Fernandez-Merodo et al. (2007), it is assumed that the chemical reaction governing the dissolution of calcium carbonate forming the porous carbonate rock is given by the global equation: (2) where the subindexes (g), (aq) and (s) refer to gas, aqueous and solid phases, respectively. The calcium carbonate species in the solid phase (CaCO 3(s) ) are assumed to be fixed, while the acid ions H 3 O + (aq) are mobile. The transport equation for both species can be written as: where K 1 and K 2 are two constants characterizing the reaction rate and " [.]'' indicates the amount of chemical concentration (i.e. solute mass per unit volume of water). The algorithm concerning the convection-diffusion reactive problem, developed for 2D problems by Fernandez-Merodo et al. (2007), was extended to the 3D case and tested here for the first time. Details of the numerical integration procedure are given in Appendix 1.

Modelling the Mechanical Behaviour
As shown in Fig. 4a where the first sketches of Prof. Nova are reported (Nova 1992(Nova , 1997, rock behaviour can be mathematically described by extending the strain-hardening plasticity framework used to model soils by increasing the size of the yield locus giving the material tensile strength. In particular, the intersections of the yield locus (f = 0) with the positive and negative p' axis are called p c = p s + p m , and p t; p s acts as preconsolidation pressure, as it is for soils, and it is assumed to depend on both plastic volumetric and deviatoric strain rates. In contrast, p t and p m account for the effects of interparticle bonding. The existence of bonds therefore implies a nonzero tensile strength (p t > 0) (Di Prisco et al. 1992) and an increase in the yield stress along radial loading paths. p m (4) Fig. 4 a Concept of weathering for bonded geomaterials induced by shrinkage of the yield locus proposed by Nova (1992Nova ( , 1997 and b, c current version of the weatheinrg model driven by the STD and LTD mechanisms (Ciantia and di Prisco 2015)

Chemo-Hydro-Mechanical FEM Coupling
Despite that the mechanical part uses a Lagrangian approximation and the transport part an Eulerian one, the key point of the coupling is that both solvers use the same finite element mesh. Coupled formulation in GeHoMadrid uses for instance elements of kind T6P3 (triangle), Q8P4 (quadrilateral) in two dimensions or H8P4 (hexahedra), B20P4 (bricks) in three dimensions. The unknowns of the ''transport solver'' are the concentrations of CaCO 3 and H 3 O + , Eq. (4). The coupling uses a ''staggered'' algorithm which works as follows: knowing the seepage velocity and the concentration value of CaCO 3 and H 3 O + for each node at time t, the transport solver computes the new values at time t + dt. The normalized dissolved mass, ξ dis , is calculated using Eq. (26) in Appendix 2, and it is used as an input in the hydro-mechanical solver affecting the mechanical response (see Eq. 9). For a partially saturated sample, a water retention curve needs to be defined, which, knowing the pore pressure, allows the saturation degree S r , that governs the STD to be determined. As described in Fernandez-Merodo et al. (2007), since the two phenomena analysed have different time steps, the computations are arranged in such a way that there is an inner loop with small time steps for transport within the hydro-mechanical main loop, and several time steps are computed with the transport solver for each hydro-mechanical step. Another simplifying assumption can be made by not considering changes in seepage velocity; if the seepage velocity is constant, there is no change in pore pressure. The mechanical behaviour can be computed with a simple displacement formulation. For sake of simplicity, in this work the weathering path considered is as follows: during saturation, because the inundation phase is assumed to be rapid enough to prevent any dissolution, only the STD process is assumed to take place; consequently, ξ dis is constant, while S r evolves from 0 to 1. As detailed in Ciantia et al. (2015b), STD evolves for a limited range of degree of saturation. Once S r > S r,cr dissolution is assumed to start, and hence LTD starts evolving with ξ dis . Although features such as pressure solution (Rutter and Elliott 1976), opening of cracks and discontinuities (Lollino and Andriani 2017), corresponds to this increase in the case of isotropic compression path. Although p t and p m can generally be considered as two independent hardening variables, they are usually assumed to be linearly related, p t = kp m (Lagioia and Nova 1995). To incorporate weathering, Nova and co-workers (Nova et al. 2003) introduced a scalar measure (going from 1 to 0) of weathering, namely the weathering function, Y. The latter described the shrinkage of the yield locus induced by nonmechanical processes (Lumb 1962) in a phenomenological way. Recently, micro-and macro-experimental evidence has shown the existence of two types of bonding, depositional and diagenetic, Ciantia et al. (2015a, b) and Ciantia and di Prisco (2015) redefined the weathering function by means of a multiscale approach (Fig. 4b, c). This allowed the weathering to be described using physical quantities namely the degree of saturation (S r ) and the normalized dissolved mass ( dis ) of calcium carbonate (Fig. 5). A detailed description of the Nova model is given in Ciantia and di Prisco (2015), and hereafter the hardening laws: and the constitutive relationship: are reported to ease the discussion on the numerical results. Y and Y E represent the strength and stiffness weathering functions, which are the key functions that link the strength and stiffness evolution to the weathering driving variables, S r and dis . In Appendix 2 the full expression of the weathering functions is reported. Fig. 5 Microscopic weakening processes: a saturation-induced short-term debonding of depositional bonds; b, c dissolutiondriven long-term debonding and grain dissolution process (Ciantia et al. 2015a) chemical enhanced crack propagation (Hueckel et al. 2017) and strain localization cannot be captured by GHM and the adopted constitutive model, the definition and use of the BDI(t) would remain unchanged even if more advanced numerical methods such as FDEM (Geomechanica Inc 2016) or P-FEM (Rodriguez et al. 2016) were to be used. Note that in the latter case a coupled constitutive relationship would still be required.

Experimental Results and Model Calibration
The material tested in the laboratory was retrieved from the man-made cave of the selected case study (represented in Fig. 1) that was excavated in a calcarenite deposit in the urban area of Canosa di Puglia. The calcarenite is characterized by a porosity of 59% and a unit weight in a dry condition of 11.1 kN/m 3 . The experimental tests performed to adequately characterize the soft rock and calibrate the constitutive model are: (i) Thin sections, scanning electron microscope (SEM) imaging with punctual energy-dispersive X-ray spectrometer (EDS) measurements to identify the chemi-cal composition and the presence of depositional and diagenetic bonds (Fig. 6). (ii) Uniaxial (UC), Brazilian (BT) and drained triaxial tests (TXD), on saturated samples to define the shape of the yield locus (Fig. 7). (iii) UC and BT tests on dry (S r = 0), partially saturated (0 < S r < 1) and saturated (S r ~ 1) samples to calibrate the strength (Y) and stiffness (Y E ) weathering functions for the STD mechanism (Fig. 8). (iv) The post-peak responses of the UC and TXD tests are used to calibrate the mechanical hardening parameters, s , t and t . (v) UC and BT tests on saturated intact and artificially (with acid) pre-weathered samples to calibrate the strength and stiffness weathering function for the LTD mechanism, are performed (Fig. 9).
Following Ciantia et al. (2015a), by analysing the thin sections with both a polarized light microscope and SEM, it was possible to identify the presence of both depositional and diagenetic bonds (Fig. 6). Punctual EDS analyses (Fig. 6) performed during SEM imaging of the thin sections confirm that both type of bonds are formed by CaCO 3 microcrystals. While UC and BT tests are very useful to identify the yield locus at low confinement pressures, the TXD test are needed to identify the yield locus shape at high pressures (Fig. 7). In Fig. 8a typical axial stress-strain curves at different degrees of saturation are shown; the strength and stiffness reduction associated with higher saturation degrees are clearly visible. Figure 8b, c shows the normalized trend of σ c and E versus S r . These results have been obtained by performing more than 40 UC tests on specimen with 38 mm and Young modulus with S r diameter using a displacement-controlled loading frame at 0.05 mm/min. The LTD process has been analysed by means of a two-step experimental test: following the same procedure of Castellanza and Nova (2004), the samples were first weathered by means of an acid solution (pH = 2.8) in order to accelerate the weathering process. In this initial phase, the dissolved mass ΔM is the measured variable and ξ dis is evaluated by means of Eq. (26). Subsequently, the samples are tested, and in Fig. 9b, c, the dependence of both strength and stiffness on ξ dis is illustrated, while in Fig. 9a

Model Validation at the Laboratory Scale
As it was already mentioned, the reliability of the approach is strongly affected by various aspects, but the most critical is the correct numerical reproduction of weathering mechanism taking place. The capability of the model to capture rock weathering from a quantitative point of view at the laboratory scale is then a prerequisite. As the CHM initial and boundary conditions in the laboratory can be controlled with reasonable accuracy, before using the approach as a predicting tool, the calibrated model should first capture the laboratory-scale response. For this reason, two BVP laboratory tests were performed and simulated with the already calibrated model. The first (BVP #1) is a model footing test on dry and wet calcarenite, while the second (BVP #2) is a loaded small-scale pillar subject to both the STD and LTD weathering mechanisms. For this latter the experimental results are taken from Ciantia et al. (2015b). The calcarenite samples of the two BVP's were prepared from a Calcarenite di Gravina Formation block extracted few miles from the analysed cavity, hence have different initial strength conditions (p m0 ), but the rest of the model parameters are the same.

BVP#1 Model Footing Tests on Dry and Wet Calcarenite
BVP#1 is aimed to prove the capability of both the constitutive model and numerical code to properly simulate the short-term water-weakening effects (STD) on the bearing capacity of a small-scale footing. Two experimental tests were run, and the results are shown Figs. 10 and 11a, while the FE simulation results are reported in Fig. 11b. Both experiment and numerical simulation are performed controlling the vertical displacement of the rigid footing. It is worth noting the formation of a punching mechanism transforming the calcarenite into an unbonded soil just located in a bulb below the foundation (see Fig. 10d-f). The experimental load-displacement curve is characterized by an initial linear path until a specific point described as the generalized yielding point (see Fig. 11a) followed by a nonlinear hardening characterized by an increase in bearing capacity as the vertical displacement, s v, increases (Nova and Parma 2011). Concerning the numerical simulation, it can be observed that the dry and wet behaviour of the calcarenite is very well captured by the model in terms of load displacement curves for the different environmental conditions (i.e. dry and wet) which are considered by assuming a value of S r = 0 and 1 for the dry and wet foundation, respectively. The FE results also show that, as the prescribed vertical settlement of the rigid footing increases, the stresses below the foundation increase and that the yielding of the rock mainly develops in tabular zone that propagates in depth forming a plastic bulb. This can be observed in Fig. 11b, c, in terms of evolution of the bonding internal variable p m and volumetric plastic strains, respectively. Interestingly, there is a very good match in terms of size and shape of the plastic bulb determined numerically (Fig. 11b, c) and experimentally at the end of the test (Fig. 10e, f).

BVP#2: STD and LTD Processes Inducing Model Pillar Failure
The second BVP simulation is devoted to reproduce the failure of a pillar induced by a LTD process occurring after a STD process (see Fig. 12a). The experimental test simulates the degradation processes that might occur in a rock pillar within an underground cavity, when the rock is first soaked by water and then subjected to chemical weathering (Castellanza et al. 2008). As the applied vertical stress is smaller than the saturated strength of the calcarenite, the STD will not induce failure, whereas similar experiments where the pillar failure is induced by STD are reported in Ciantia et al. (2015b). The testing layout and an illustration of the different stages of the test are shown in Fig. 12a, b, respectively. Figure 12c shows the problem geometry, the spatial discretization adopted in the FE simulation and the assumed boundary conditions. In the first stage of the test (OA), the specimen is loaded in uniaxial conditions, starting from a dry state. Then the specimen is saturated at constant axial stress by imposing 1e−6 0.9 1.657 1.657 10 − 0.1 5 0 Table 3 Weathering parameters ξ dis,cr (-) S r,cr (-) cd,0 (kPa) cw,0 (kPa) E d,0 (MPa) E w,0 (MPa) 0.4 0.25 1600 ± 100 740 ± 100 315.0 ± 100 150.0 ± 50 Table 4 Elastic parameters 45 ± 5 0.09 ± 0.02 1000 5 ± 1 Fig. 10 a Experimental apparatus; b-d typical punching behaviour; e, f pictures of typical debonded bulbs just below the foundation an upward water flow (stage AB), up to full saturation (stage BC) during a period of 2 min. It can be observed how the sample deforms due to the stiffness reduction experienced during the STD (Fig. 13a). In the final stage of the test (CD), the specimen is brought to failure by a long-term exposure to an acid solution (pH = 3.7) on the lateral surface. The isochrones of the computed S r , ξ dis and p m are shown in Fig. 14.
The first three columns (0 to 10 min) show the effect of the STD, driven by an increase in S r from the bottom of the specimen, on the boding-related internal variable p m . The last three columns (20 to 70 min) show the effect of the LTD, driven by an increase in ξ dis from the lateral boundary of the specimen, on p m . The response of the material during the first loading stage-performed in steps-is almost linearly elastic. The subsequent STD produces additional deformations during the wetting stage. The progressive chemical dissolution of intergranular bonds gives rise to progressively increasing deformations with time at constant stress, eventually leading to the failure of the specimen when the reduction in uniaxial rock strength due to the LTD is such that the axial load cannot be sustained any longer. The FE simulations, as demonstrated by the comparison between experimental observations and numerical predictions, are used to calibrate the diffusion tensor D ij to correctly capture this behaviour. In the simulation, the dissolution constants K 1 and K 2 are taken from the work of Ciantia (2013) that obtained them performing dissolution tests on the same calcarenite. The computed isochrones of the axial stress (σ y ) and of ξ dis along a radius at the specimen mid-height are shown in Fig. 13b, c, respectively. As the acid flows into the specimen from the lateral boundary, the values of ξ dis are not uniform along the section, but minimum at the specimen axis and maximum at the boundary. The inhomogeneous degradation process resulting from the horizontal acid flow causes a progressive weakening of the external portion of the specimen that induces a reduction in axial stress close to the lateral boundary and a corresponding stress increase at the specimen axis, due to stress redistribution within the specimen. This phenomenon is qualitatively similar to the one observed in the pillars of abandoned mines (Fig. 12d) or in ancient temple columns (Fig. 12e) when the intact rock is exposed to long-term chemical degradation (Parise and Lollino 2011;Tamagnini et al. 2002).

Damage Analysis of Buildings Built on Top of an Abandoned Cavity in Canosa di Puglia
The BDI(t) of the complex 3D geometry concerning the interaction of a cave (excavated about two centuries ago) and two buildings in a PG3 zone (Fig. 1) is here addressed. It should be noted that such situations are extremely common  . 11 a Experimental results and numerical simulation in (S r = 0) and wet (S r = 1) conditions; b evolution of vertical plastic strains and c bonding variable p m under the model footing in the Apulian region (Parise and Vennari 2017), but this particular case study is considered as the authors were directly involved with the risk assessment of this cavity (Ciantia et al. 2015c) and hence disposed of in situ material needed for the model calibration (see Sect. 4). The current global stability of the same geosystem determined by means of the standard c-ϕ reduction method (Potts et al. 1990) resulted to be 3.4 and 1.5 for the dry and saturated states, respectively (Ciantia et al. 2015c). Despite this current stable condition in both dry and saturated conditions, it is known that the long-term weathering of the rock could lead to subsidence and eventually to a failure mechanism such as the one in Fig. 2 (Ghabezloo and Pouya 2006). The BDI(t) as introduced in Sect. 2 is hence used to tentatively predict the evolution of the damage on the buildings under study. The weathering scenario assumed consists of two phases: 1. a saturation phase (STD process) that could be induced by heavy rainfalls or leakages from water supplies. 2. a long-term dissolution phase (LTD process) assuming the cavity remains saturated for long periods of time in an open-system scenario (dissolved calcite ions are washed away).
The laboratory tests on the calcarenite collected in situ were used to calibrate the initial mechanical properties of the material which, assuming them to be the same for all the domains, were assigned as the model's initial state. During phase 1, simulated by setting S r = 1 and hence using the 'wet' mechanical properties for the calcarenite, the geosystem remains in equilibrium. For this reason, in what follows, only the long-term chemical damage is considered (phase 2). In the absence of environmental hydraulic field data and boundary conditions, inspired by Ghabezloo and Pouya (2006), the cave boundaries were exposed to a different acidic pH environment. Such boundary conditions trigger the development of a coupled chemo-mechanical diffusion and deformation process. To consider the real geometry of the problem and realistic building loads, a 3D numerical analysis is performed, and the mesh and problem definition are shown in Fig. 15a, b. To highlight the importance of the three-dimensionality of the problem, a 2D plain strain analysis in correspondence of section A-A′ is also performed. The mesh of section A-A′ and the problem definition are shown in Fig. 15c. For the 3D FEM model discretization, 46,961 quadratic tetrahedral H10P4C4 elements (67,581 nodes and 169,270 degrees of freedom) are used, while the 2D model is meshed with 1879 triangular T6P3C3 elements (3920 nodes and 8642 degrees of freedom). Standard displacement boundary conditions have been applied to the external boundaries, with fixed horizontal and vertical displacements at the bottom and fixed horizontal displacements on the lateral boundaries. The external loads (i.e. the foundation loads) acting at the ground surface are kept constant, while the LTD process evolves from the inner boundary of the cave. To represent building A (see Fig. 1), a pressure of 150 kPa is imposed under each of the 9 square footings, while a pressure to represent building B a uniform pressure of 40 kPa is imposed under the area covered by it. For the 2D analysis, the same pressure of building A is applied despite the 2D geometrical limitations (plain strain hypothesis). On the cavity surfaces, a fixed constant concentration of [H 3 O + ] ions, corresponding to a pH of 6.7, is imposed. They will propagate into the rock by means of a diffusive phenomenon causing the long-term dissolution of CaCO 3 . The reactive transport model parameters used for the simulations are those of BVP#2 and are summarized in Table 5. To highlight the importance of the cavity environmental conditions, two further 3D simulations were run: one imposing a pH of 7.5, while in the other, the pH imposed at the cavity boundaries was set to 6. Figures 16 and 17 report the main results as a function of physical time (days) from the beginning of the weathering process assumed (i.e. the LTD process) for the 2D and 3D Fig. 14 FE simulation of the small-scale pillar weathering test: isochrones of S r (first row), ξ d (second row)and p m (third row) with time. Each snapshot represents half a sample as represented in Fig. 12c simulations, respectively, for the 6.7 pH simulations. The dissolution process that evolves from the inner boundary is shown in the first column of Fig. 16 by the contour in time and space of the contaminant ion [H 3 O + ]. As the contaminant concentration increases, the dissolution process evolves causing a mass reduction and a consequent increase in the scalar ξ dis leading to a reduction in the bonding-related variables p m and p t (2nd column of Figs. 16 and 17). This strength reduction causes yielding of the rock to concentrate in correspondence of the partition wall for both the 2D and 3D simulations. As the exposure time with the [H 3 O + ] ions increases, the plastic strains localize into a shear band (3rd column) and the displacements increase (4th column). Nevertheless, before plastic strains localize into a shear band, as clearly shown in Fig. 18a, material degradation developing from the wall boundaries induces a stress redistribution in the partition wall Fig. 18b. Although non-local approaches (Conte et al. 2010) should be adopted to avoid mesh dependency when looking at the failure mechanism, it was observed that further mesh refinement did not change the results of the simulation in terms of surface displacement with time. Interestingly, the failure of the partition wall does not cause a global failure mechanism, and this is due to the residual frictional strength within the shear band combined with the rooftop geometry. In fact, as represented in Fig. 18c, d, the stress corresponding to a point within the shear band starts from an initial elastic state, yields after about 120 days of acid exposure and evolves until critical state is reached. Figure 19 shows comparison of the displacement field for the 2D and 3D numerical simulations used to calculate the BDI for the 6.7 pH environment. The chemo-mechanical coupled induced displacements are well documented in terms Fig. 15 a, b 3D mesh and sections considered for the BDI calculation, and c section, mesh, boundary and initial conditions and environmental loads for 2D simulation     Fig. 19d, e. As expected, the 2D simulation deforms more than the 3D one thanks to the lateral support provided by the rock not intersected by the tunnel. The change in curvature of the displacement time curve (more visible for the 2D case) can be identified as the sinkhole inception time as from this point onward the displacement will accelerate until the simulations stops converging. Interestingly, the surface subsidence observed in the 3D simulation acquires a circular shape, and it is caused by the failure of the partition wall. Finally, in Fig. 20a the angular distortion horizontal strain paths for the 2 and 3D analyses with different pH conditions are reported. By means of Eq. (1), in Fig. 20b, the BDI evolving with time is also reported. The time needed to attain a BDI of one results to be 1080, 200 and 60 days for the 7.5, 6.7 and 6 pH case, respectively. The plot also shows how the three-dimensional effect delays the BDI attaining a value of unity from 165 to 200 days for the 6.7 pH simulation. In this contribution, only two sections were considered for the BDI assessment of building A; however, any section can be analysed and the final BDI(t) would be the first one (section A-A′ in this simulation) reaching a value of 1.

Conclusions
Many abandoned man-made caves induce sinkholes because of the degradation of the constituent soft rock. In this paper, the problem superficial subsidence prior to sinkhole formation in urban area is tackled and in particular: (i) A procedure to determine the temporal evolution of settlement damage to buildings by means of a building damage indicator, namely the building damage index, BDI is proposed. The weathering-induced temporal evolution of horizontal strain and angular distortion are used to introduce physical time to the BDI which is a practical scalar index of the evolution of serviceability of buildings. (ii) Despite all the limitations intrinsic to standard FE and local constitutive relationships, for the first time, a 3D numerical approach including the reactive transport of chemical species and an advanced strain-hardening plasticity constitutive model was used to tentatively predict the BDI(t). (iii) To introduce a time-dependent building damage index, a weathering scenario must be assumed and consequently modelled. For this reason, the reactive transport of chemical species is fundamental as it enables to describe the spatial-temporal evolution of chemical species that reproduce numerically the assumed weathering scenario. The reliability of the  approach increases as the weathering mechanism reproduced is more realistic and with the soundness of the constitutive model which has to be hydrochemo-mechanically coupled. It is shown how BDI is strongly dependent on the weathering scenario assumed. (iv) In order to determine the BDI(t) of a real building on an abandoned cavity in Southern Italy, standard mechanical and advanced chemo-hydro-mechanical (CHM) experimental tests on the in situ rock are performed and used to calibrate the constitutive model which accounts for the hydro-chemical damage effects induced by weathering.
(v) Two model footing tests on a dry and wet calcarenite were performed, and a further small-scale boundary value experimental test from the literature is simulated to test the soundness of the model and the FE code to deal with CHM-coupled boundary value problems. (vi) For all simulations of the LTD-induced damage of the abandoned anthropic cavity, the BDI attained a value of 1 well before the progressive failure mechanisms starts to develop. For the same weathering scenario (6.7 pH case), smaller displacements are observed in the 3D model with respect to the 2D analysis and the BDI attains a value of one after a longer time. (vii) A plane strain 2D analysis followed by three more realistic 3D FEM stability analyses is performed. In all cases, the numerical results show that weathering of the porous rock causes displacements at the ground level that, before accelerating, slowly increases with time. Plastic strains begin to develop at the corners of the cavity by the partition wall well before they start to localize into a shear band. This plastic strain evolution induces a lot of stress redistribution within the rock. Despite the mesh dependency limitations, this result may suggest that such local failures can be very useful as collapse premonitory signs as in the model they start developing few months before failure. This is in accordance with what observed in the 2007 Gallipoli sinkhole where evident signs of damage on pillars were observed months before the collapse (Delle Rose 2007). (viii) Despite the non-symmetric geometry of the two cavities, the subsidence basin in the 3D analysis assumes a circular shape, which is typical for the sinkholes observed in these rocks (Parise and Lollino 2011). As recently observed by Nof et al. (2013) and Jones and Blom (2014) with interferometric synthetic aperture radar (inSAR), subsidence basins can be used a precursory signs of sinkhole formation. Interestingly, also the 3D FEM simulation presented here captures this feature.
(ix) The numerical predictions show that, for open-system scenarios, it results that serviceability conditions are lost after 65 days of exposure in acidic environment (pH 6), whereas the same level of damage is achieved after 3 years if the pH is set to 7.5. From an engineering perspective, the precautions to be taken to maintain the buildings in good serviceability conditions would be to monitor flood occurrence and to avoid the closing of ventilation raises. In fact, for this particular case, rock saturation an the corresponding STD would still not cause collapse and proper air ventilation would avoid the formation of aggressive acidic environments that, in the presence of water, would initiate the LTD.
Limitations of the BDI include aspects inherent to most deterministic numerical models of geotechnical problems. In particular the reliability of the approach will depend on: (i) the amount (and associated uncertainty) of information required concerning the geometry of the problem, (ii) the CHM material properties, the chemical reactions and processes, (iii) local heterogeneities and (iv) the assumption of the weathering scenarios. Many of these inconveniences may be addressed by employing a probabilistic approach (Griffiths et al. 2002;Griffiths and Fenton 2004;Zhang and Ng 2005). We also believe that in a near future remote sensing techniques such as the interferometric synthetic aperture radar (inSAR) from satellite images will provide temporal evolution of subsidence that will allow to test and improve the reliability of the proposed methodology. described in (Donea 1984;Löhner et al. 1984) already applied to both fluid dynamics and solid mechanics (Mabssout and Pastor 2003a, b;Peraire et al. 1986;Quecedo and Pastor 2002) and a system of ordinary differential equations solved with a fourth-order Runge-Kutta algorithm (Hirsch 1988;Lambert 1973). In particular, denoting with ∇ the vector differential operator, Eq. (4) is written in a more compact manner as: where we have introduced the vector of unknowns The mass balance equation includes convection, diffusion and retardation effect due to the chemical reaction between calcite and acid ions. Equation (11) combines both parabolic and hyperbolic PDEs with a source term S, and it can be solved using a splitting operator technique in which each of the two operators, convective-diffusive transport and sources, is treated separately (Toro 2001): We want to evolve n from time t = t n to the new value n+1 at t = t n+1 by a time step Δt = t n+1 − t n . Therefore, two problems will be considered: (i) A diffusion-advection problem of the type PDE ∶ t + ∇ − ∇ 2 = 0 initial condition, IC ∶ ( , t n ) = n ( ) Δt ⇒ adv−dif (ii) and the source problem, which is an ordinary differential equation (ODE): The initial condition IC of the advection-diffusion problem (13) is the initial condition n of the full problem, and the solution after a time step Δt is adv−dif , which may be regarded as a predicted solution. ODE (14) is then solved by a time Δt and with initial condition adv−dif , and the solution of Eq. (14) is then regarded as the approximate solution to the full problem. One may choose the best scheme for each type of sub-problem.

Adopted Taylor-Galerkin Scheme for Advective-Diffusive Transport
where = ∫ T ϕ ϕ d is the classical mass matrix and Δ̂ = n+1 − n . It is important to remark that the system of equations can be written in the form and can be solved using an iterative Jacobi or "accelerated viscous relaxation scheme". Both are described in (Zienkiewicz and Taylor 2000). The iterations are: where L is the lumped mass matrix, which is diagonal. The method can be considered as explicit, and provides an accurate solution in less than 5 iterations.

Adopted Runge-Kutta, Fourth-Order Accurate Scheme for the Source Term
The second part of the splitting concerns the source term: Runge-Kutta (RK) algorithms provide a high accuracy in the evaluation of ODEs (Hirsch 1988;Lambert 1973). The general form of a RK algorithm is with where the consistency condition ∑ k k=1 k = 1 has to be fulfilled.
Here we have chosen a fourth-order RK algorithm, as it provides an excellent combination of accuracy and computational effort. The coefficients are given by: from where we obtain As ( ) = ⋅ with constant:

Appendix 2
As described in Ciantia and di Prisco (2015) the two driving variables of weathering of calcarenite are the normalized dissolved mass ξ dis and the degree of saturation S r . ξ dis , is the dissolution reaction progress variable calculated as where ΔM and M 0 are calcium carbonate's dissolved and initial mass, respectively. The weathering function Y governs the size of the yield locus when the material is subject to environmental loads and is a scalar that goes from 1 (intact) to 0 (completely weathered) and can be written as: where While the stiffness weathering function, Y E , governs the stiffness of the material when subject to environmental loads and results: where (24) 1 = n 2 = n + 1 2 Δt ( 1 ) 3 = n + 1 2 Δt ( 2 ) 4 = n + Δt ( 3 ) n+1 = n + Δt 6 ( 1 ) + 2 ( 2 ) + 2 ( 3 ) + ( 4 ) In Eqs. (28)-(32), w c0 and d c0 represent the uniaxial compression strengths of calcarenite rock under wet and dry conditions, respectively; E d 0 and E w 0 are the Young's moduli of the intact material in dry and wet states, respectively; S r,cr is the minimum degree of saturation necessary for all the depositional bonds to suspend, while ξ dis,cr is the normalized dissolved mass that corresponds to the ideal threshold between the bonded and the granular material state (all bonds have been dissolved).