Seismic Vulnerability Assessment of Masonry Structures in Archaeological Sites
The main aim of this dissertation, by the author Jinwoo Kim, is to develop the modeling strategies for masonry structures with the Discrete Element Method (DEM), which has been verified for its high suitability for the dynamic analysis of ancient masonry structures.
Seismic risk management for monuments
Cultural heritage is a valuable asset in our society, but many of them are located in earthquake-prone areas throughout the world. This often imposes significant seismic risks to the monuments and historic constructions. Especially, unreinforced masonry (URM) is one of the most common yet vulnerable types of heritage structures.
To safeguard the priceless cultural heritage in our society, it is essential to better understand the seismic behavior of historic masonry structures and monuments.
Dynamic analysis of masonry structures
To date, such an important task still remains challenging to engineers despite the advanced modern engineering technology. This usually stems from the unconventional characteristics of monumental structures, often found incomplete as free-standing columns or partially remaining walls. It is predictable that such incomplete structures will present different mechanical properties than structural elements of buildings in function do.
The computational tool of choice for structural analysis has been probably the Finite Element Method. It is based on continuum idealizations and capable for the analysis of 'strength and elasticity', solving a variety of engineering problems in practice. However, the fact that masonry structures are inherently discontinuous makes the popular method less suitable for masonry analysis.
The DEM, on the other hand, is based on discontinuum idealizations and adequately addresses the discontinuous nature of masonry structures (José V. Lemos, 2007). Moreover, the equation of motion allows relatively large displacements with separation and rotation compared to the FEM. This is especially beneficial for predicting complex failure modes of block structures due to 'stability'.
Problem Statement
Despite the clear advantages, there seem to be major challenges in the DE analysis of monumental masonry structures, such as:
- The necessary modeling parameters are difficult to estimate;
- The monumental structures in incomplete forms have different mechanical behavior than structural elements of buildings in function.
Objectives and scope of the study
This dissertation aims to advance the DE modeling strategies of masonry structures in archaeological sites. The new strategies must be able to determine appropriate modeling parameters without having to harm the existing structures. Also, the DE analysis will focus on the instability of masonry structures due to ground motions. Finally, the outcomes are expected to help the seismic vulnerability assessment procedures of our cultural heritage.
Two cases of masonry structures in archaeological sites are studied:
1) the Pompeii Water Towers in Pompei, Italy
2) the Apollo Nymphaeum in Hierapolis, Turkey.
These case studies are the continuation of the ongoing research projects at the Department of Civil, Environmental, and Architectural Engineering at the University of Padua. The scope of this study encompasses the extension of the previous work, in addition to a general review of the DE modeling parameters derivation. For the DE analysis, the commercial three-dimensional DEM software 3DEC is used.
Discrete Element Method (DEM) for masonry structures
The Discrete Element Method (DEM), also called as Distinct Element Method, is a numerical method to solve the motion of discontinuous medium. The fundamental algorithms of the DEM were based on the work by Cundall an Strack (1979). Initially, it was mainly developed for the analysis of geomaterials, such as jointed rock, soil, or granular materials. In the following decades, there has been a great advancement in computational power, and the DEM has become far more versatile in various disciplines.
Rocking motion of rigid body
The first analytical solution to the rocking response of a single rigid block subjected to horizontal ground motions was established by Housner (1963). The author derived an equation of motion for pure rocking of a single rigid block and studied the scale effect of the slenderness ratio.
Out-of-plane behavior of masonry walls
In addition to free-standing columns, masonry walls with out-of-plane behavior are one of the typical structures of concern. Peña & Lourenço (2006) modeled the rocking walls as rigid blocks in the DEM software UDEC and compared results with the experimental data. Although it was a relatively simplified two-dimensional model, the model was valid to predict the out-of-plane behavior and yielded good agreement.
DE MODELING STRATEGIES
Input parameters
The main modeling parameters required in the DE analysis are mainly broken down into two types:
1) Block properties:
- Density;
- Elastic modulus;
- Shear modulus.
2) Interface properties:
- Joint normal stiffness;
- Joint shear stiffness;
- Friction angle;
- Dilatancy angle;
- Cohesion;
- Tensile strength.
Joint stiffness parameters derivation
In the DEM, elastic materials are modeled as equivalent springs in the interface, which is already a common practice in the field. The first systematic approach to such a modeling method for masonry structures was established by Lourenço (1996). The author describes a homogenization technique to calculate the joint stiffness of masonry walls from the elastic modulus of the unit and mortar. Figure 8 provides a schematic drawing of the homogenization technique for composite material.

Damping constants derivation
The DEM allows implementing equivalent viscous damping in the equation of motion. The numerical result could be very sensitive to the damping parameters, which makes determining appropriate damping parameters a critical factor to the accuracy of the model. In 3DEC, damping is introduced as equivalent viscous damping, in the form of Rayleigh damping. Rayleigh damping consists of two components; mass-proportional and stiffness-proportional.
Critical time step
The solution scheme for the DEM is based on the central difference algorithm and is conditionally stable. A critical time step that satisfies the numerical stability criterion is determined as the minimum between the required time step for block deformation computation and inter-block displacement calculation. When the rigid block model is used, and only the critical time step for inter-block displacement is considered:
Evaluation of the modeling parameters derivation
To determine the joint stiffness parameters, specifically, the elastic and shear modulus are required for an estimation. However, the mechanical properties of historic structures in archaeological sites are oftentimes not available and also can be highly variable. This leaves high uncertainties in the required properties.
This dissertation, therefore, includes the novel approach by which the parameters in the computational model are calibrated by the experimental results from the ambient vibration tests. The calibrated values will be compared with the parameters by the formulation, to discuss its effectiveness. In addition, zero damping will be used for the cases introduced herein, because of its capability to predict the seismic response during strong motion. Nevertheless, the effect of different damping constants suggested in the previously mentioned researches will be compared for a sensitivity study.
THE FREE-STANDING MASONRY COLUMNS
The Pompeii Water Towers are part of the urban water supply system in Roman Pompeii. In the Roman age, water was carried by public aqueducts to Porta Vesuvio and stored in the water distribution building Castellum Aquae (Figure 11).

The stored water would then run through the fourteen Water Towers located in the city. Olsson carried out an extensive archaeological study on the water distribution system and reconstructed a plausible original operational scenario. The Water Towers, connected at different elevations, are supposed to be connected by lead pipes and distributed water to street fountains and private houses in the city.
In this dissertation, only the four Water Towers are studied and will be referred to as WT1, WT2, WT3, and WT4.
In addition to the archaeological literature, a survey campaign was carried out with the aim to better understand the current conservation state and inner composition of the WTs. The survey included:
- Non-destructive testing (NDT), such as Sonic Tomographic Tests and Ground Penetrating Radar (GPR) investigations, to characterize the mechanical properties of the materials;
- Ambient vibration tests
- Output-only identification techniques
The results provided the modal parameters of the Water Towers, such as natural frequencies, damping ratio, and mode shapes.
Modeling parameters
Table 2 summarizes the modeling parameters determined for the WTs.

Calibration of joint stiffness parameters
The joint stiffness parameters are calibrated using the experimental data from ambient vibration tests. Figure 16 illustrates the mode shapes of the first five natural frequencies, obtained by the eigenmodes analysis in 3DEC and the experimental data.

Evaluation of calibrated parameters
To evaluate the effectiveness of the calibration, the elastic modulus and shear modulus are estimated by back-calculation and compared with typical ranges of values recommended by the Italian code NTC 2008.
Would you like to know more? Leave your email below to receive the full academic dissertation behind this article:
Artificial earthquakes
A series of artificial earthquakes are generated from the elastic spectrum for time history analysis.
Table 6 summarizes the seismic parameters necessary for characterizing the site effect to generate code complying accelerograms. First, the soil type and topography were determined based on the results from the geological survey campaign and the classifications in the Italian code NTC 2008.
Subsequently, Ultimate Limit States (SLV; Life Limit State and SLC; Collapse Limit State) with various life spans and classes are considered, as summarized in Table 7.

For historic structures, it is subjectable to determine the design requirements due to different interpretations on the usage and importance of the structure. Generating artificial earthquakes from various spectra is also desired since the Italian code recommends using different shapes with varying corner periods.
Time history analysis
A series of time history analyses were carried out in 3DEC using the generated artificial earthquakes. The seismic load is applied as a velocity history prescribed at the base. For each earthquake scenario, sixteen different combinations of accelerograms are made in regards to the directionality.
First, the current state of the Water Towers was analyzed. Subsequently, additional time history analyses were carried out to study the original state of the WTs and the sensitivity of results to the damping parameters.
Fragility Curves
The overturning stability of a rigid body subjected to horizontal ground motions must be understood from a probabilistic point of view. The results of the displacement and rotational history confirm that the rocking response of the water towers is in no systematic way and not consistent. The fragility curves are constructed on the basis of the multiple stripe analysis results from the Time history analysis.
Seismic vulnerability
The fragility curves of all WTs are compared to identify the most vulnerable WT (Figure 25). WT3 shows the fragility curve surpassing the others, while WT2 is the most stable and found in a lower range. WT1 and WT4 show a similar curve, while WT4 appears slightly more vulnerable to moderate earthquake motion (up to PGA = 0.5 g).

Finally, the collapse probability of Water Towers is assessed according to the local seismic condition. The Pompeii Water Towers are located in a moderate seismic zone, where the seismic hazard map in Figure 26 defines a probable earthquake with a PGA of 0.125 to 0.200. At a PGA of 0.2 g, for example, the collapse probability is approximately 10 % for WT3 and nearly zero for the other WTs.

The original state
According to the archaeological studies on the site, water containers once existed on top of the Water Towers to store water. The containers are modeled accordingly in 3DEC to study the ancient state of the water towers as shown in Figure 27.

Damping parameters
In this section, three different damping constants are compared with zero damping:
- D1) using the empirical equation;
- D2) specifying 50 % damping at the edge frequency;
- D3) specifying 100 % damping at the edge frequency
THE DRY-JOINTS MASONRY WALLS
The Apollo Nymphaeum studied in this dissertation is a subset of the Temple of Apollo in Hierapolis, Turkey (Figure 33). To the current archaeological knowledge, the Apollo Nymphaeum was built at the beginning of the 3rd century, while its function has long been misunderstood since its recognition. The maps in Figure 33 display the Nymphaeum located at the western edge of the temple and surrounding buildings. As can be seen in the images, the temple is now in ruin with the only remaining part of three continuous walls. Stone remnants found on the ground imply the original structure may have collapsed.

Figure 34 displays a front view from one of the LiDAR images of the Nymphaeum. Currently seen, however, is only the remaining part of the original Nymphaeum. The reconstructive drawing by archaeological study in Figure 35 reveals that the Nymphaeum used to have two floors and narrow terraces, supported by a series of columns. The small hole found in the niche is believed to accommodate a water foundation in the inner space.

Structural features
The Apollo Nymphaeum consists of three continuous walls, North Wall, East Wall, and West Wall resting on multiple layers of base and foundation, making a C-shaped section in plan. The base consists of three masonry layers made of the same stone. Found at the bottom of the Nymphaeum is a foundation made of pebble and cobble aggregates mixed with lime mortar.
The Apollo Nymphaeum is made of dry-joints stone masonry walls as shown in Figure 38. Such construction technique is typical of the Greek age. The masonry layout reveals that the horizontal bed joints are continuous around the walls, while the vertical bed joints are partially staggered. They are perfectly cut, and there is no mortar or very thin bed joint.

Geometrical modeling
Modeling a complex geometry is oftentimes painstaking and time-consuming when applied to large and irregular structures. 3DEC provides a useful command that generates random joints using a statistical model, without having to draw the exact geometry of units and interfaces one by one. With this command, it is possible to assign joints with normally distributed spacing by inputting a mean and a standard deviation of the joint spacing.
Modeling parameters
Table 10 summarizes the modeling parameters for the Hierapolis nymphaeum in 3DEC.

Artificial earthquakes
Table 11 summarizes the site-specific seismic parameters to define the target spectrum. Artificial earthquakes are generated from the elastic spectrum in accordance with Eurocode 8 (Table 12).

Time history analysis
Time history analysis is carried out in 3DEC to assess the seismic behavior of the Hierapolis Nymphaeum. Figure 45 displays the DE model at the end of ground motions (t = 25 s) at each PGA level.

Seismic vulnerability
Figure 50 (a) displays the maximum normalized horizontal displacement with a clear trend of monotonic increase. Overall, the out-of-plane behavior at all sections is below the threshold value 1.

CONCLUSIONS
The free-standing columns: the Pompeii Water Towers
To an earthquake with a PGA of 0.2 g, WT3 has approximately 10 % overturning probability and is identified to be the most vulnerable. The other WTs have nearly zero probability of an earthquake of the same PGA.
At the PGA same to the seismic capacity calculated by the kinematic analysis, WTs have very low overturning probabilities below 10 %. Also, WT4 was identified as the most vulnerable structure.
It can be said the DE analysis provided more detailed information about the dynamic stability and better addressed the actual 3D behavior.
The dry-joints masonry walls: the Apollo Nymphaeum
The seismic vulnerability of the Apollo Nymphaeum is assessed based on the qualitative observations of the failure modes and the displacement history. At a PGA of 0.78 g, the Nymphaeum showed residual displacements that leaves gaps between blocks. The first visible failure occurred at an extremely high PGA of 1.28 g. At the highest PGA of 1.28 g, the maximum normalized horizontal displacement approached 0.8. This is equivalent to the center of mass displaced by 80 % of the thickness of one leaf (≈ 1/3 of the wall thickness).
It is concluded that the Apollo Nymphaeum has sufficient seismic capacity to an earthquake with a PGA of 0.3 g, which is defined in the national hazard map. However, although the monument will not collapse, there is a concern with residual deformation.
DE Modeling Strategies
Overall, it is concluded that the DEM code can predict the dynamic response efficiently and provide valuable insight for seismic vulnerability assessment.
Future works
The main challenges of the DE analysis and the seismic vulnerability assessment were the unconventional characteristics of monumental structures. This made it difficult to determine the engineering parameters. If the modeling strategy were to be applied to other general cases, there is a great need for further research to fill this gap.
The currently available performance criteria are mostly based on the strength, not the stability thus may not be relevant to the overturning stability. In this dissertation, the threshold value to separate the collapse and the non-collapses was set based on the scatter plots of PGA against normalized rotation. This can be repeated for other cases, but once robust performance criteria are established, the seismic vulnerability can be assessed for various levels of damage.
A modeling technique to replicate a complex irregular masonry layout is introduced. In 3DEC, it is possible to generate a hundred cases of random joints from one normal distribution model. However, further verification may be necessary to study the difference resulting from the different random generations.
One of the greatest challenges of assessing the seismic behavior of the Hierapolis Nymphaeum was the lack of information regarding its engineering properties and construction details. During the modeling procedure, it was difficult to find sufficient references for material properties of local travertine as an ancient building material. Moreover, the unknown construction details inside the wall forced to make the three-leaf wall assumption, which left extra uncertainties in reproducing the actual mechanical behavior. To tackle these challenges, experimental data or an additional survey campaign may be required.