
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)11745-8
10.1016/j.heliyon.2024.e35714
e35714
Research Article
Integrated approach for seismic wave prediction and structural evaluation of Oya tuff quarry underground spaces: Field observations and numerical modeling
Cheng Chuantao Chengchuantao@fuji.waseda.jp
a⁎
Shizuo Noguchi b
Dintwe Tumelo K.M. c
Seiki Takafumi d
Iwanami Motoi a
a Department of Social and Environmental Engineering, Waseda University, Tokyo, 169-8555, Japan
b Public Interest Incorporated Association Oya Regional Maintenance Corporation, Utsunomiya, 321-0043, Tochigi, Japan
c Department of Civil Engineering, Faculty of Engineering, University Botswana, Gaborone, 999-106, Botswana
d Department of Construction Engineering, Utsunomiya University, Utsunomiya, 321-0904, Tochigi, Japan
⁎ Corresponding author. Chengchuantao@fuji.waseda.jp
06 8 2024
30 8 2024
06 8 2024
10 16 e3571419 3 2024
6 7 2024
1 8 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
This study evaluates the structural stability of large underground spaces in seismic conditions, represented by the Oya underground stone mining plant. By directly monitoring the seismic response of the underground mining site, significant earthquake activities at the plate boundaries of the Tokyo region and Ibaraki Prefecture offshore area were observed. Additionally, through an in-depth analysis of seismic records from different locations within the underground structure, the dynamic characteristics and motion patterns of the Oya underground stone mining plant were revealed, revealing its movement trajectory during earthquakes. Additionally, this study innovatively applied seismic waves measured at the original site as input parameters and artificially generated seismic waves based on their response spectra. A numerical analysis was performed after ensuring the model's high correlation with the original site was met. The findings demonstrate that the results of both parameter input methods are confirmable and valuable. Under severe seismic conditions, instability was observed in some regions of the underground mining site. The study also discusses the location and damage mechanisms of the mining site's structure under seismic effects, providing valuable insights for the safety assessment of similar large underground spaces and proposing new approaches for selecting input parameters in seismic analysis.

Highlights

• Tokyo area and the plate boundary (the offshore area of Ibaraki Prefecture) are seismically active.

• The vibrations followed a logarithmic decline rate with the constant distance divergence.

• Instability may happen in some parts of the underground quarry in a seismic solid state.

• During the seismic state, some instability may be evident in some parts of the underground quarry.

Keywords

Incremental dynamic analysis
Seismic wave prediction
Dynamic characteristics of underground structures
And safety evaluation standards for underground space
==== Body
pmc1 Introduction

Most earthquakes occur in zones near the boundary of tectonic plates (massive plate-shaped bedrock covering the earth's surface). The offshore area of Ibaraki Prefecture is a typical area where plate movement occurs (the ground of seismic source). The seismic probability of an earthquake striking this area is more than 90% in the next 30 years, and the possibility of a strong earthquake cannot be ruled out, still [1]. If a strong earthquake occurs, Utsunomiya City, Tochigi Prefecture, may be significantly affected, and this is where our study subjects are located. There are more than 200 underground quarries in Oya town, Utsunomiya City, which together form a large-scale underground cavity group. To reuse the abandoned resources, the local government has renovated some underground quarries, and some have become sightseeing places in Utsunomiya. In 2018, they were chosen as Japan's material and cultural heritage. However, many large-scale collapse accidents of underground quarry spaces have occurred in Oya town since 1989 [2]. Therefore, it is necessary to explore and understand the influence of ground motion on the stability of underground structures from the viewpoint of effective utilization.

Large-scale underground spaces like the Oya Underground Quarry are widely distributed globally, with Japan being particularly abundant in such sites, including but not limited to Ikeshima Coal Mine Remains (Nagasaki City), Iwami Ginzan Silver Mine (Oda City, Shimane Prefecture), Sado Kinzan Remains (Sado City, Niigata Prefecture), and Super-Kamiokande (Hida City, Gifu Prefecture). Since the March 11, 2011 seismic event, commonly known as the Great East Japan Earthquake, there has been an increasing emphasis on the safety of large underground spaces. With the continuous publication of studies on solid earthquake prediction, local governments, and residents are increasingly concerned about the seismic capacity of underground spaces. They aim to ensure underground structures can safely withstand future strong earthquakes.

When assessing the seismic bearing capacity of underground structures, reliance primarily rests on three scientific research methods: in-situ observation, model testing, and theoretical analysis. Specifically, in-situ observation encompasses seismic Monitoring and disaster investigations, with seismic observation being particularly crucial for continuous Monitoring of the underground structure's motion over multiple seismic cycles. This study focuses extensively on seismic observation to accurately capture the dynamic motion state of underground structures and collect raw seismic data. As for model testing, it involves shake table experiments and artificial source tests. However, this method is not considered in this study due to uncertainties in boundary conditions, doubts about model reliability, and challenges in data acquisition. On the theoretical analysis front, numerical methods such as the finite element, finite difference, and discrete element methods are predominantly utilized. Despite significant advancements in these methods in recent years, research on the seismic input mechanisms in underground space seismic response analysis remains insufficient. Previous studies often employed “typical” seismic waves, randomly selected model boundaries, or directly applied seismic motion time histories at the bedrock surface as inputs, placing them at the bottom boundary of numerical models. However, due to the presence of site amplification effects [3] and site depth attenuation effects [4], the application of these methods is not rigorous [5]. Therefore, this study, considering seismic input mechanisms comprehensively, conducts a more realistic dynamic analysis to ensure the reliability and accuracy of research results.

When delving into the seismic response characteristics of underground structures, the primary and fundamental task is to establish appropriate seismic input parameters. These parameters' time history inputs primarily originate from two primary sources: the naturally recorded seismic data accurately captured and recorded by actual seismographs and the seismic data synthesized artificially based on scientific principles and methods. Together, these two sources form the basis of time history inputs, providing a solid foundation for the theoretical analysis of seismic response characteristics in underground structures.

(1) The naturally recorded seismic data was accurately captured and recorded by actual seismographs.

In existing seismic response analyses, a common practice is directly applying seismic acceleration records from databases such as Ki-net and PEER-NGA to the bottom of computational models, often based on ground acceleration amplitudes corresponding to seismic intensities specified in seismic design codes. However, this approach has significant limitations and may lead to overestimating seismic response results [6]. Similarly, direct use of seismic acceleration records from the bedrock surface is also not ideal, as it confuses the concept of bedrock seismic motion with boundary input seismic motion in numerical computations. Arbitrary boundary selection may lead to a series of issues [5], as seismic wave propagation to the ground surface results in radiation-damping effects [7], with some energy dissipating into the far-field foundation. Furthermore, in addition to the amplification effects inherent to the strata, the presence of shallow subsurface structures also contributes to the amplification of seismic signals [8]. These factors significantly influence the computed results of seismic response in large underground spaces, resulting in poor comparability between different seismic analysis results for the same site.

To address these issues, the Incremental Dynamic Analysis (IDA) method has emerged [9][10]. This method is a Performance-Based Earthquake Engineering (PBEE) analysis approach [11], capable of effectively analyzing the entire process from stability to instability of engineering structures under seismic actions while considering the randomness of seismic motion. Through the IDA method, we can more accurately assess the performance of underground structures under seismic actions [12][13], providing robust support for engineering design and safety assessment.

(2) seismic data synthesized artificially based on scientific principles and methods.

The seismic motion process exhibits typical non-stationary characteristics, dynamically transitioning from weak to strong and gradually weakening. Early researchers employed stabilization techniques, such as uniform modulation or equivalent stabilization, to simulate this non-stationary behavior. Pioneering work by Bollotnill et al. in 1960 [14] introduced the concept of capturing the non-stationary intensity characteristics of ground motion through the product of a stable Gaussian process and a deterministic envelope function. Subsequently, researchers have developed various intensity envelope models [15][16] and conducted systematic and in-depth studies on model parameters' determination and decay patterns.

Among various methods for synthesizing non-stationary acceleration, the approach centered on fitting target response spectra is widely employed. The uniform modulation triangle series and model proposed by Scanlan and Sachs in 1974 have emerged as the primary means of seismic motion simulation due to their extensive use and significant value in engineering practice. With the continuous accumulation of experiences from major earthquakes, non-stationary frequency seismic wave input has become a crucial factor in ensuring the accuracy of computational results. While researchers such as Nigam et al. [17] and Yokoyama [18] have conducted in-depth discussions on the phase characteristics of seismic motion time series, existing methods still struggle to accurately capture the dynamic changes in frequency content across different periods. Boore et al. [19] proposed using phase derivatives to synthesize seismic motion, aiming to enhance simulation accuracy. Additionally, Shrikhande et al. [20] introduced arbitrarily varying phase angles to generate piecewise linear curves, further enriching the means of simulating phase unwrapping or phase curve changes in seismic motion.

As research progresses, artificial seismic wave synthesis has been achieved by precisely matching seismic records to target site response spectra using parameters and criteria set by multiple databases such as PEER-NGA [21][22], offering novel perspectives in earthquake engineering. In recent years, with the widespread adoption of machine learning techniques, more advanced matching technologies have made the seismic wave input process increasingly efficient [23][24][25].

However, despite these methods' flexibility in selecting seismic wave search parameters, differences exist between the selected seismic and original site records. No study has utilized completely accurate original site seismic records for artificial seismic wave synthesis, which may be attributed to a shortage of seismic records, lack of records near the engineering site, and mismatched geological conditions between the recorded site and the engineering site.

Therefore, this study comprehensively evaluates the seismic performance of underground quarry spaces by integrating both actual seismic recordings and artificially synthesized seismic waves. A comprehensive assessment of the seismic performance of underground quarry spaces was conducted by directly Monitoring seismic records at the original site and using actual seismic instrument data combined with Incremental Dynamic Analysis (IDA) and related theoretical analysis methods. Additionally, spectral information from the original site seismic records was further utilized to synthesize artificial seismic waves as input parameters for numerical analysis, aiming to reconcile the contradiction between accurate structural analysis methods and rough design seismic input. The study found that the results obtained from both methods corroborated each other, achieving satisfactory evaluation outcomes and providing valuable references for stability assessments of similar structures.

2 Prototype site

The quarry underground space is in Oya Town, Utsunomiya City, Tochigi Prefecture, north of Tokyo. It is about 8 km from east to west, and its area is 4.13 km2, has a depth of 200 m, and has a reserve of about 1 billion tons. There are more than 200 underground quarries in the area. The mined geological feature is the soft rock called the Oya tuff, a sedimentary rock. It is characterized by lightweight, high porosity, and low hardness. The rhyolitic ash was deposited in the Miocene epoch around the Aquitanian stage, about 20 Ma. The stratigraphic sequence of Oya (Fig. 1) is composed of Paleozoic basement sandstone and chert strata covered with tuff, which can be divided into three horizons according to the rock strength: upper (Oya layer II-III), middle (Oya layer IV-V), and lower (S2,VI). The underground space of the tuff quarry aims to extract quality stone, which is widely used in house construction and interior decoration. The subterranean exploitation of tuff rock formations primarily employs the vertical shaft and residual pillar mining technique, prominently observed within the quarry's underground space. This subterranean void presents itself as a vertical pit measuring 60 meters in depth (Fig. 2(a)), 30 meters in length, and 28 meters in width. Access is established from the base of the vertical pit (Fig. 2(b)), giving rise to an underground expanse spanning 160 meters in length and 28 meters in width (Fig. 2(c)). The upper portion is surrounded by soil and upheld by 14 residual pillars (Fig. 2(d)).Figure 1 The stratigraphic sequence of the Oya formation near the study site.

Figure 1

Figure 2 The field map and design drawing of the Oya tuff quarry underground space: (a) The vertical shaft; (b) The entrance; (c) The entry to the entrance to the underground quarry; and (d) The design drawing.

Figure 2

3 Methods

The seismic performance of underground structures is studied through field observations, experiments, and numerical analysis methods. The primary aim of the research is to assess whether the subject of the study can withstand a major earthquake that may occur in the future. Seismic observation techniques are used to analyze the movement patterns of underground structures during earthquakes and to record the spectral characteristics of natural seismic waves. These spectral data are further utilized to generate artificial seismic waves, providing a scientific basis for earthquake prediction. Field experiments validate the reliability of the model through stress wave attenuation experiments and explore the stress wave attenuation patterns in the underground space of the Oya Tuff quarry. In terms of numerical analysis, the finite difference method effectively handles large deformations and complex nonlinear problems by substituting derivatives in partial differential equations with difference approximations and integrating over time steps. Therefore, this study employs the finite difference method to establish a numerical model and uses artificial seismic waves to perform dynamic analysis of the model. Additionally, a dynamic stability assessment benchmark is established using the IDA method to evaluate the impact of artificial seismic waves of varying intensities on the model analysis results (as shown in Fig. 3), thereby predicting whether the underground structures' seismic performance is sufficient to cope with potential future major earthquakes.Figure 3 The methodology.

Figure 3

3.1 Observations on the prototype

To facilitate the observations of the seismic records on the wall of the underground space, the seismographs, monitoring the three perpendicular directions investigated, were set from each other at the upper part, or -11.5 m from the ground level, GL, the central part, or -34.5 m below the GL, and the lower part, or -59.5 m below the GL of the shaft (Fig. 4).Figure 4 Ground motion observation points in the underground space.

Figure 4

3.2 Analysis of the model

In the field of geology, Oya stone is classified as a pumiceous lapilli tuff. It is widely believed that during the first half of the Miocene, most of the Japanese archipelago remained below sea level, with frequent volcanic activity. The volcanic ash and sand emitted during these eruptions settled in the marine environment and gradually solidified, eventually forming the Oya stone seen today [26]. To systematically investigate the main mechanical properties of the pumiceous Oya stone, previous research [27] has conducted numerous indoor rock mechanics tests in the Oya area. Samples were collected from the surrounding rock and residual pillars of the Oya area, from the surface to the bedrock, using drilling methods. The tests conducted included physical tests, uniaxial compression tests, fracture tensile tests, elasticity tests, and triaxial compression tests. Due to space limitations, detailed results of some tests are presented in Table 1, while average results are shown in Table 2. To ensure the uniformity of geological parameters, this study selected the average values of each test parameter. These parameters were then converted to the International System of Units and calculated as input parameters for the numerical model, with the specific calculation process detailed below. The determination of the bulk modulus K and shear modulus G is based on the secant modulus of elasticity E and Poisson's ratio ν, obtained from the indoor tests shown in Table 2. These are calculated using the following mechanical formulas, as shown in Eq. (1) and Eq. (2):(1) K=E3(1−2ν)

(2) G=E2(1+ν)

Additionally, the selection of cohesion C and the internal friction angle ϕ was also based on specific data from the indoor tests shown in Table 2. Specifically, this study used the results of the consolidated undrained shear test Ccu, setting the value of cohesion C at 2.10 MPa. while the internal friction angle ϕ was chosen as an integer value between ϕcu and ϕcd, which is 30∘. As for the tensile strength σt and density ρt, they were directly taken from the data obtained in the indoor tests of Table 2. The specific parameters are shown in Table 3.Table 1 Detailed data table.

Table 1	Apparent	
specific	
gravity	

	Effective	
porosity	

	Uniaxial	
compressive	
strength	

	Fracture	
tensile	
test	

	Tangential modulus	
of elasticity	
(×104)	

	Secant modulus	
of elasticity	
(×104)	

	Poisson's	
ratio	

	
%	kgf/cm2	kgf/cm2	kgf/cm2	kgf/cm2	
Inside the quarry	
⋅ Residual pillar	
(40∼60 m)	

	1.77	37.3	78.0	11.5	3.16	2.31	0.20	
1.78	39.8	69.8	8.3	3.22	2.28	0.24	
1.79	40.8	84.4	11.8	3.38	2.47	0.24	
		79.0	12.5				
	
Inside the quarry	
⋅ wall	
Block sampling	
(40∼60 m)	

	1.82	39.3	81.4	10.6	3.22	2.44	0.23	
1.85	37.7	78.4	13.4	3.32	2.39	0.21	
1.85	37.0	82.6	8.8	3.43	2.38	0.21	
		83.2		3.01	2.30	0.22	
		82.7		2.74	2.23	0.24	
		81.8		2.92	1.93	0.24	
	
Boring from	
ground surface	
⋅ shallow	
(15∼23 m)	

	1.65	45.1	46.4					
1.78	38.8	84.5					
1.80	35.6	79.8					
1.73	49.3	57.0					
1.83	41.9	80.8					
1.84	41.3	75.6					
	
Boring from	
ground surface	
⋅ Residual pillar	
(30∼50 m)	

	1.66	44.7	87.2					
1.72	38.6	90.6					
1.73	37.9	75.4					
1.78	44.4	86.4					
1.78	43.4	94.7					
1.79	42.7	93.0					
	
average	1.77	40.8	82.5	11.0	3.11	2.28	0.23	
	
Standard deviation	

	0.06	3.5	5.9	1.7	0.24	0.17	0.01	
	
Number of samples	

	18	18	20	7	6	6	6	

Table 2 The engineering parameters of the underground area.

Table 2Physical and	
mechanical constants	

	Test results	Test item	
Apparent specific gravity	1.77	Physical testing	
Effective porosity (%)	

	40.08	
Uniaxial compressive	
strength (MPa)	

	8.09	Uniaxial compression	
test	

	
Poisson's ratio	0.23	
Tangential modulus (MPa)	3.05Image 1 103	
Secant modulus (MPa)	

	2.24 Image 1 103	
P wave velocity (km/s)	2.40	Elastic test	
S wave velocity (km/s)	1.23	
Fracture tensile strength (MPa)	1.08	Tensile strength	
test	

	
Adhesion C (MPa)	Ccu	2.10	Triaxial compression	
test	

	
Ccd	2.31	
Internal friction	
angle Φ (deg.)	

	ϕcu	31.3	
ϕcd	29.6	

Table 3 The parameters of the analysis of the model.

Table 3Bulk	
modulus	

	Shear	
modulus	

	Cohesion	Internal	
friction	
angle	

	Tensile	
strength	

	Density	Poisson's	
ratio	

	
(GPa)	(GPa)	(MPa)	(deg.)	(MPa)	(kg/m3)		
1.38	0.91	2.1	30	1.08	1770	0.23	

This study utilized FLAC 3D software, based on the finite difference method, to construct an underground space analysis model measuring 500 m x 500 m x 100 m (Fig. 5) for subsequent analysis. According to related reports [28], when the confining pressure on Oya tuff exceeds approximately 5.88 MPa, the peak strength does not significantly decrease, exhibiting elastic-perfectly plastic behavior. Based on the data from Table 3, this study first conducted a numerical simulation of a uniaxial compression test on a specimen measuring 50 mm x 50 mm x 100 mm, using elastic-perfectly plastic behavior as the failure criterion. The simulation results (Fig. 6) indicate that the material entered a fully plastic state upon reaching stress of 7.2 MPa, consistent with previous reports. Consequently, in subsequent analyses, a dynamic analysis was performed using an elastic-perfectly plastic model that integrates both the Mohr-Coulomb criterion and tensile failure standards.Figure 5 Ground motion observation points in the underground space.

Figure 5

Figure 6 The elastic-perfectly plastic model.

Figure 6

To optimize computational efficiency, hexahedral elements measuring 2 m x 2 m x 2 m were strategically placed near the underground space, while finer hexahedral elements of 0.1 m x 0.1 m x 0.1 m were utilized in the vibration attenuation test area. The rest of the grid was set to 10 m x 10 m x 10 m, totaling 373,930 elements. Considering the vertical shaft and residual pillar mining technique employed in the Oya tuff underground quarry, the voids created by mining are supported by residual pillars formed from unmined areas, as illustrated in Fig. 2 and Fig. 5. These residual pillars along with the surrounding rock constitute a continuous whole, thus eliminating the need for additional contact models in this simulation. Furthermore, an infinite boundary condition was implemented to negate the effects of reflected waves at the model's sides. This setup and parameter selection not only ensure the accuracy and computational efficiency of the analysis but also provide a robust tool for subsequent assessments of the seismic performance of underground structures.

3.3 The stress wave decay experiments

The replicability for investigating the propagation characteristics of elastic waves vibrating on the inhomogeneous Oya tuff by numerical analysis was evaluated in this study. Experiments were conducted to measure the propagation characteristics of the stress waves in the Oya tuff quarry underground space. The device for measuring elastic wave acceleration and the distance between them are shown in Fig. 7. To study and understand the distance attenuation characteristics of acceleration propagation, the distance between accelerometers was 0.3 m, 0.6 m, 1.2 m, 2.4 m, and 4.8 m, respectively. The direction of the connecting line between the accelerometer and the vibration-adding point (the direction of vibration propagation) was set as X-axis +, and the vertically upward direction was set as Z-axis +. The target frequency range of measurement and recording was less than 500 Hz, and the sampling frequency was 10 kHz. During the measurements, about 1.5 seconds was recorded with a sampling frequency of 10 kHz. Since the vibration of the impact hammer was artificially vibrated, the measurements were done 10 times under the same conditions, and the best data set was selected according to the vibration load waveform and vibration propagation. The numerical simulation was done by finding the position corresponding to the experiment's site in the model established earlier. The stress waves were observed by the elastic wave measurement devices one after another after the acceleration of the Impact hammer. In this study, the stress waves observed by the first stress wave measurement device were used as the input condition of the corresponding node in the numerical simulation. As the stress wave propagated, the energy decayed continuously, and the stress wave peaks of the five observation points correspond to the experiment (Fig. 7 Model of the experiment).Figure 7 The experiment plot (reappear measuring line).

Figure 7

3.4 The seismic resistance evaluation based on Incremental Dynamic Analysis (IDA)

3.4.1 Overview of IDA

This method evaluates the properties of the structure continuously mined from the nonlinear time history response of the gradually increasing input ground motion intensity. The evaluation criteria for the stability of large-scale underground space can be obtained according to the determined ground motion level and property level (Table 4) [29].Table 4 The seismic resistance evaluation criteria.

Table 4Seismic level	Property level	
Operating basis earthquake (OBE)	Serviceability	
Safety evaluation earthquake (SEE)	Safety	

3.4.2 Performance parameter selection

IDA curve is a nonlinear curve generally depicted in the IM-DM two-dimensional coordinate system. The seismic intensity coefficient IM represents the seismic intensity in the two-dimensional IM-DM coordinate system in IDA. In this study, PGA is the multiple of IM, and DM is the index of seismic performance parameters in the IM-DM coordinate system. The choice of DM depends on the purpose and characteristics of the structure and is often expressed in terms of displacement and deformation. The performance of the underground space needs to be considered to maintain sufficient stability for continuously mining Oya stone in the underground quarry space. In addition, relative deformation occurs in the surrounding bedrock of the underground space of the Oya quarry during earthquakes. Therefore, considering the safety of continuous mining and the sense of security for surrounding residents, the maximum relative displacement angle, θ max (Fig. 8), Eq. (3) of the upper and lower layers of underground space is used as DM:(3) θmax=δh=δ40−δ60h

Figure 8 Diagram of the interlayer displacement angle.

Figure 8

θmax: maximum interlayer displacement Angle

δ40: maximum layer displacement at 40 m underground

δ60: maximum layer displacement at 60 m underground

h: height

DM: [(maximum relative displacement angle) = (maximum displacement of the upper layer-maximum displacement of the lower layer)/layer difference in underground space (Fig. 8)]. In addition, the numerical analysis model determined the position where the displacement difference was the largest.

3.4.3 Selection of the input seismic wave based on field measurements

The seismic records with offshore plate boundaries were selected as the focus (Table 5) based on the ground motion records measured at the bottom of the quarry space. The records were multiplied by the multiple changes with PGA (peak ground acceleration) at all intensities (Table 6) according to the relationship between magnitude and PGA. The upper limit of PGA (Grade 8) was set to amplitude modulation to the underground space failure state.Table 5 The set of epicenters as the target of the analysis.

Table 5	Latitude (deg.)	Longitude (deg.)	Dep (km)	Mag	
1	37.251	141.501	33	1.9	
2	37.701	141.619	45	6.4	
3	36.869	141.134	49	3.6	
4	34.836	139.851	51	3.6	

Table 6 Target values of the magnitude and peak acceleration.

Table 6Intensity	Acceleration (m/s2)	Acceleration/gravity ((m/s2)/g)	
1	0∼0.008	0∼0.000815	
3	0.025∼0.08	0.002548∼0.008155	
4	0.08∼0.25	0.008155∼0.025484	
5-	0.25∼0.8	0.025484∼0.081549	
5+	0.8∼2.5	0.081549∼0.254842	
6-	2.5∼4	0.254842∼0.407747	
6+	4∼5	0.407747∼0.51	
7+	5∼	0.51∼	

3.5 Analysis of the ground motion made by the trigonometric series superposition method

Various studies have used open-source data, such as PEER NGA-West2 in the United States and Ki-NET in Japan, measured on the surface as the input time history of ground motion. However, this will lead to the input of seismic time history, which the structural characteristics of the study object will not determine. It was found that the seismic performance evaluation results obtained with different input seismic time histories were very different, even if the seismic intensity was the same. Due to the amplification characteristics of the site, when the (near) surface observation data is input into the ground motion time history as the base, it does not reflect the accurate seismic performance evaluation results to determine whether the seismic capacity of the Oya quarrying underground space can cope with the strong ground motion that may occur in the future, In this study, the average spectrum of the seismic waves observed at the Oya quarry underground space was analyzed as an element. The synthetic trigonometric series superposition method was used to synthesize artificial seismic waves as the input ground motion time history, which effectively avoided the field amplification effect and fully matched the structural characteristics of the research object. This was combined with the near-field seismic environment and site conditions to meet the requirements of the most unfavorable conditions.

3.5.1 Overview

The seismic records with offshore plate boundaries were selected as the focus (Table 5) based on the ground motion records measured at the bottom of the quarry space. The records were multiplied by the multiple changes with PGA (peak ground acceleration) at all intensities (Table 6) according to the relationship between magnitude and PGA. The upper limit of PGA (Grade 8) was set to amplitude modulation to the underground space failure state.

The basic idea of this method is to calculate the data following a stationary Gaussian approximation by using trigonometric series and then multiply it by the seismic intensity envelope to obtain the non-stationary acceleration time history waveform [14]. The specific calculation process is shown in Eq. (4):(4) X=f(t)⋅a(t)

f(t): seismic intensity envelope.

a(t): constant acceleration moments waveforms.(5) a(t)=∑k=0nCkcos⁡(ωk+ϕk)

In the equation for the steady-state acceleration time-course waveform a(t), presented in Eq. (5), ϕk represents the uniformly distributed random phase within the interval (0,2π), while Ck and ωk denote the amplitude and frequency of the k-th frequency component, respectively. Let the amplitude of the component closest to the two highest frequencies be equal to the amplitude of the two largest frequency components.

where Ck is obtained from the given power spectral density function, as shown in Eq. (6):(6) Ck=[4S(ωk)⋅Δω]1/2

where ω=2πT and ωk=2πkT. where T is the total holding time, and the conversion relationship between the reaction spectrum and the power spectrum is given by Eq. (7):(7) S(ω)=ξπω⋅[SaT(ω)]2ln⁡[−πωTln⁡(1−P)]

where SaT(ω) is the response spectrum for a given target acceleration, and ξ and P represent the damping ratio and the probability of exceeding the response, respectively. Generally, it is assumed that P≤0.1512, and in this study, P=0.15 is used.

f(t) is the envelope function proposed to account for the non-stationarity of the acceleration, denoted as Eq. (8):(8) f(t)={(tt1)20≤t≤t11t1≤t≤t2e−c(t−t2)t2≤t≤T

where t1 is the rising interval of the peak, t1 ∼ t2 is the smooth segment of the peak, and T represents the duration. c is the decay constant.

The initial seismic wave obtained by the above method is bound to have a large error between its response spectrum and the target response spectrum. This discrepancy necessitates a method to correct the initial wave so that its response spectrum gradually approaches the target response spectrum.

The relative error between the calculated response spectrum S(ω) and the observed seismic response spectrum SaT(ω) is corrected according to the methodology proposed in the related study. The correction equation is as follows, as shown in Eq. (9):(9) E(ωk)=|sa(ωk)−saT(ωk)saT(ωk)|

The amplitude adjustment factor is given by Eq. (10):(10) Rω(ωk)=|SaT(ωk)Sa(ωk)|

When E(ωk) is larger than the allowable error ϵ, the power spectrum and amplitude spectrum of the k-th fitting point must be corrected. The correction formula for the power spectrum is given by Eq. (11):(11) Si+1(ωk)=[SaT(ωk)Sa(ωk)]2⋅Si(ωk)

The amplitude spectrum is corrected as shown in Eq. (12):(12) Ci+1=Ci[R(ωi)]q

In the equation, Ci represents the amplitude corresponding to the points around the k-th control point, including the k-th point itself. The variable q is determined based on the contribution of the individual components of the Fourier spectrum to the maximum response; q takes the value of 1 or -1, depending on whether the contribution is positive or negative, respectively.

The target peak acceleration was calculated from the prediction map of seismic probability and seismic intensity [1]. The mean value of the response spectrum recorded in seismic observations with the offshore plate boundary as the earthquake epicenter was taken as the target spectrum (Fig. 9). The prediction graph of earthquake probability and seismic intensity established the equivalent magnitude. Concerning several focal areas occurring in the offshore area of Ibaraki Prefecture, the coordinates of the earthquake epicenter were assumed to be about 36°6.5' north latitude and 141°15.9' east longitude, and the horizontal distance from the Oya underground space was 133.4 km. After establishing the equivalent epicenter distance, a regression equation was derived using earthquake observation data statistics to determine the time parameters required by the earthquake intensity envelope function [24]. The specific formula is shown in Eq. (13):(13) log10⁡t1=−1.074+1.005log10⁡(D+10)log10⁡(t2−t1)=−2.268+0.326M+0.5815log10⁡(D+10)log10⁡C=1.941−0.28M−0.5670log10⁡(D+10)

where M is the equivalent earthquake magnitude, and D is the equivalent epicenter distance.Figure 9 The bottom three-component spectrum of the underground space.

Figure 9

3.5.2 The unsteady ground motion in Oya

To construct a seismic time history waveform, it is necessary, in general, to evaluate the three characteristics of ground motion: the amplitude (peak acceleration), the spectrum, and time. The target peak acceleration was calculated from the prediction map of seismic probability and seismic intensity [1]. Here, the response spectrum as the target spectrum was obtained from the analysis, and the equivalent magnitude and epicenter distance were calculated through the prediction map of seismic probability and seismic intensity. The regression equation was established according to the firm ground motion observation data statistics, and the parameters required for the ground motion intensity envelope function were determined. Considering the most unfavorable state of the structure, the dynamic stability evaluation benchmark was established, and the magnitude of 7.2 (Table 7) and the probability of exceeding 2% in 50 years (Table 8) were selected. The ground motion made according to this step is shown in Fig. 10.Table 7 Parameters of the ground motion intensity envelope.

Table 7D (km)	133.4	
M	6.7	6.8	6.9	7	7.1	7.2	
C	0.069	0.065	0.06	0.06	0.054	0.05	
t1 (s)	12.34	
t2 (s)	27.17	28.12	29.35	30.98	32.32	34.24	

Table 8 Horizontal peak acceleration of the foundation of the site.

Table 8	30 years	50 years	
Value exceeding	
probability (%)	

	3%	6%	2%	5%	10%	39%	
Seismic intensity	
Less than Grade(<)	
More than Grade(>)	

	<6	<6	<6	<6	<6	>5	
Value of maximum	
acceleration of	
foundation	
(cm/s2)	

	272.81	256.17	296.56	272.63	255.8	223.84	

Figure 10 Synthetic acceleration time history waveform of the input ground motion (along the X (a), Y (b), Z (c) components).

Figure 10

4 Results and discussion

4.1 Analysis of monitoring results

4.1.1 Monitoring of epicenter of the earthquake

The offshore area of Ibaraki Prefecture, where the plates meet, has the highest probability of earthquakes in Japan, with a probability of over 90% of a major earthquake occurring in the next 30 years. The seismic intensity of the plates is predicted to be 6.7 to 7.2. There is also the possibility of a giant earthquake. In the event of such an earthquake, there is a 14% chance that the Oya region will be affected by a magnitude six or higher. During more than a month of field monitoring (Fig. 11), it was observed that the seismic activity along the boundary of the Pacific plate and the North American plate in the Japan Trench area of Japan was highly active and was mainly concentrated in the area of latitude 36° N to 37° N and longitude 140.5° E to 142° E. Thus, it verifies the speculation that there will be a massive earthquake in the offshore area of Ibaraki Prefecture, where the plates meet, in the next 30 years. This study also attempted to reproduce the artificial seismic waveforms centered in this region. According to the monitoring data, seismic activity under Tokyo is also widespread, and the possibility of Tokyo being affected by a strong earthquake in the future is very high.Figure 11 Seismic epicenter location and magnitude observed in the Oya tuff quarry underground space.

Figure 11

4.1.2 Monitoring of epicenter of the earthquake

Fig. 12 illustrates the variation in the ratio of velocity amplitudes at each observation point along the vertical direction of the shaft wall. This ratio is calculated based on the average maximum amplitude recorded at each observation point compared to the monitoring values of the seismic instruments in the upper part. The research findings indicate that, compared to the vertical dynamic component, the shaking of the horizontal dynamic component significantly intensifies when approaching the middle height of the vertical shaft wall. However, no observed phenomenon of site amplification extending to the surface was observed, mainly due to concrete reinforcement measures taken by local safety authorities at a depth of approximately 11 meters underground to prevent collapse or landslide at the shaft exit. Nevertheless, due to the age of these measures, relevant data cannot be verified.Figure 12 The ratio of the velocity amplitude of three directions at observation points.

Figure 12

Assuming complete restraint of movement in the upper part of the excavation wall as a reference, we compared the relative movements of the middle and lower parts. Specifically, the horizontal shaking in the middle part is approximately 1.95 times that of the upper part, while in the lower part, it is about 1.5 times that of the upper part. Additionally, we found that the shaking magnitude of the vertical dynamic component is roughly half that of the horizontal dynamic element.

It is worth noting that the seismic monitoring platforms Ki-net and Kik-net in Japan primarily monitor seismic activity at around 30 meters below the surface and the underground bedrock. However, according to the data obtained from the seismic monitors installed at the middle to upper positions in this study, using the surface and bedrock data as input parameters for numerical simulation would result in significant errors. Therefore, careful selection of seismic input parameters is essential when assessing underground structures' seismic performance to ensure the evaluation results' accuracy.

4.1.3 Monitoring of epicenter of the earthquake

To consider the influence of seismic wave propagation on residual pillars and walls, the 1-second waveform trajectory of the central vibration part of the S wave was drawn. From Fig. 13, it can be seen that the shaking of the residual pillar (Fig. 13(d)∼(f), amplitude 17) is more significant than that of the wall (Fig. 13(a)∼(c), amplitude 9), with the pillar's shaking nearly twice that of the wall, and there is also a phase difference. Particularly, from Fig. 13 (f), it is found that in the same earthquake, the trajectory of the horizontal vibration component of the residual pillar is close to an arc. This preliminary confirmation confirms that although the underground space of the Oya Tuff quarry shakes as a whole, there is also a deviation in the phase and degree of shaking of each part.Figure 13 Comparison of the three-component trajectories between residual pillar and wall.

Figure 13

4.2 The reliability test of the model

4.2.1 Monitoring of epicenter of the earthquake

The study selected seismic wave data obtained from monitoring at the bottom of the Oya Quarry, which underwent meticulous processing to remove extremely long-period components, ensuring the accuracy of its use as input parameters for the analysis model. We applied these processed data to the dynamic analysis of the underground quarry to explore its seismic response characteristics in depth. The response parameters extracted from the model results, including the positions of residual pillars and the middle and bottom of the pit wall (as shown in Fig. 4), underwent a detailed comparison with actual observations. Through comparative analysis (Fig. 14), we found that the trend of the velocity response spectrum output by the model closely matched the actual observed values, and the prominent periods of the two were essentially consistent, demonstrating the high reliability of the model. Despite the absence of historical data (This is described in section 4.1.2), which hinders the complete numerical reproduction of the concrete retaining wall at the upper part of the pit, causing discrepancies between observed and analyzed values at the upper wall, these differences are primarily due to the proximity (only 0.5 m) of the seismometer to the retaining wall. However, considering the entire pit's depth of over 60 m, length of 30 m, and width of 28 m, the presence or absence of the retaining wall has a negligible impact on the overall analysis results. The data shown in Fig. 14(d)∼(i) indicate a high consistency between the analysis values and observed values at the middle and bottom of the pit wall, confirming that the impact of the retaining wall on pit performance is minimal. The focus of this study is the interior space of the limestone underground quarry, not the auxiliary structures such as the pit, entering the vast underground space described in Fig. 2(b) and Fig. 2(c), which measures 28 m wide and 169 m long. Data from the same vertical face of the pit at -34.5 m (Fig. 14(d)∼(f)) were unaffected by the retaining wall. Furthermore, the analysis of critical control parameters such as the interlayer displacement angles (Fig. 8), and the displacement measurements of the deeper layers within the underground quarry at -40 m and -60 m, also showed no impact from the retaining wall. The comparison of internal pillars within the underground quarry (Fig. 14(a)∼(c)) further verified that the local effects of the retaining wall do not extend to deeper areas of the quarry. Given the relative stability of the geology of limestone underground quarries, not all pits are equipped with retaining walls. Therefore, the results of this study are universally significant for understanding the seismic resistance of similar limestone quarry structures. From a safety perspective, even though the impact of the retaining wall is minimal, opting not to include the retaining wall at the pit entrance in the analysis is a better practice considering the most unfavorable factors. Thus, despite the limitations of historical records, the numerical model analysis results in this study remain reliable and meaningful, taking into account broad applicability and the most adverse factors.Figure 14 Comparison of the observed and analytical values of the pillars and the middle and lower parts of the shaft.

Figure 14

4.2.2 Monitoring of epicenter of the earthquake

As shown in Fig. 15(a)∼(c), the elastic waves of the three components measured at the first point of the experimental observation points were input to the corresponding position of the first point of the analysis model. Then, the accelerations of three components were obtained at the following four observation points. Their maximum acceleration was taken, and the resulting data was analyzed. The results of field tests and numerical simulations show that the maximum acceleration tends to decrease logarithmically during vibration propagation. Due to the continuous energy divergence in the propagation process, the vibration velocity at each measuring point decreased with the increased distance from the impact point. In addition, with the increase in distance, the attenuation speed of the waveform gradually slowed down, showing fast attenuation at the proximal end and slow attenuation at the distal end, and then tended to be flat. The pair of X-axis distances was taken to show the downward trend of the field test and numerical simulation, which showed the straight-line results more clearly. The accuracy and downward trend of the simulation were evaluated according to the slope of the fitting function. The determination coefficient (R2) of the vibration attenuation-related analysis value relative to the observed value was more than 0.9 (Fig. 15(d)). The numerical simulation results agree with the experimental results. To verify the reliability of the analysis model, the analysis values of the seismic wave spectrum and stress wave attenuation were compared with the experimental values, which showed that the model used in this study had high consistency. It may be concluded that it can be applied to the dynamic stability evaluation of the underground space with this model as the target object.Figure 15 Comparison of field measurements and numerical analysis of distance attenuation of maximum acceleration: (a) X direction; (b) Y direction; (c) Z direction and (d) comparison of field measurements and numerical analysis of distance attenuation of maximum accelerations.

Figure 15

4.3 Results of IDA quartile analysis

The IDA quantile curve is obtained in this study by analyzing the model results. The curve shows in detail the whole process of the analysis model from a safe state to complete collapse (Fig. 16). Following the study by Vamvatsikos and Cornell [9], the point of 0.8 times the elastic slope, Ke, as the function guarantee point, and the point of 0.2 times Ke as the safety guarantee point was taken [27]. The results showed that if gravitational acceleration, g, was 9.81 m/s2 and the PGA of the engineering foundation surface reached 0.29 g, the acceleration was 2.88 m/s2, while the corresponding plastic zone was from 6.3% to 6.9%. This was treated as the safety guarantee point. Although the underground space would become unstable over these values, it would not collapse immediately. When PGA reaches the safety guarantee point, reinforcement supports must be implemented. The IDA curve of 75% quartile (Fig. 16) concluded that the underground space had a PGA value of magnitude greater than Grade 5. This classification is based on the seismic intensity (maximum acceleration of earthquake) of Japan, which is respectively 1, 2, 3, 4, 5-, 5+, 6-, 6+, and 7+, and the corresponding values are shown in Table 9, has a plastic zone rate of 2.22% ∼ 6.3% corresponding to an acceleration of 1.98 m/s2 at 0.2 g. Before this, it was considered that it had no significant influence on the function of the underground space. Besides, it was concluded that safety inspection and repair work were needed. Based on the results, the target of seismic grade (Table 9) was set. Therefore, when the probability of collapse is 75%, the analysis model of underground space of the tuff shows that the function (Such as mining operations by workers, tourism for visitors, and storage of critical items, among others) of the underground space cannot be guaranteed under the assumed seismic magnitude of the weak standard magnitude five and horizontal maximum acceleration of 1.98 m/s2. In addition, when the seismic wave with a standard magnitude of 6 and weak maximum acceleration of 2.88 m/s2 was used as the engineering foundation, the underground space of the tuff quarry began to lose structural stability; also, to manage the deformation of the underground space, the maximum interlayer displacement angle was adopted in this study. The maximum displacement of the quarry angle can be defined from the underground space, IDA 4.5 Image 1 10−4 rad repaired damage. With up to 0.18 Image 1 10−2 rad, the whole structure of the underground space of the tuff quarry may begin to lose stability.Figure 16 The IDA quartile curve.

Figure 16

Table 9 Target of seismic grade determination based on the analytical results of 75% quartile.

Table 9Seismic	
resistance	
sorting	

	Intensity	Peak	
Acceleration	

	Maximum	
interlayer	
displacement	
angle	

	Plastic	
state	

	
(m/s2)	(rad)	(%)	
AA (OBE)	5+	<1.98	
<0.2 g	

	4.5×10−4	<2.22	


	
Countermeasures	The damage is small, and no additional	
countermeasures are required.	

	


	
A (SEE)	5+∼6-	1.98∼2.88	
0.2 g∼0.29 g	

	4.5 ×10−4 ∼0.0018	2.22∼6.3	


	
Countermeasures	If the structure of the underground space is slightly or	
repairable damaged without affecting the function,	
its function can be improved by support.	
It needs to be inspected and repaired.	

	


	
B (Unstable)	>Grade 6-	>2.88	>0.0018	>6.3	


	
Countermeasures	Adverse effects on function. At this time,	
the underground space begins to become unstable,	
but it will not collapse immediately.	
It needs to be reinforced by backfilling or support	

	

In the numerical analysis, under the influence of a strong earthquake, cracks and collapses occurred in the target area's underground space, both in the plastic region of the deformation regime. The size and distribution of the plastic region represent the unstable state of the target underground space. Fig. 17(a)∼(d) shows the plastic zone after applying the seismic wave at the bottom of the model. Based on the field measurements of the ground motion in the underground space, the strong motion was inferred by the superposition of the trigonometric series. The results show that there is a plastic region, and the central yield morphology is caused by tension. The possibility of tensile stress failure of the roof near the residual pillar was found to be relatively high. During the earthquake, the underground space's failure position was mainly the roof's intersection, the residual pillar, and the surrounding bedrock. Fig. 17(a)∼(d) shows the distribution of the plastic zone after the gradually increasing amplitude of peak acceleration of the input artificial seismic wave just below the roof around the underground space. The distribution of the plastic zone state and the maximum displacement angle at this point are shown in Table 10. The evaluation criteria in Table 8 can be used to determine which grade the management value belongs to. According to the predicted dynamic response of Magnitude 7.2 earthquake ground motion in the offshore area of Ibaraki Prefecture to the model, all indices (Table 9) of seismic intensity of more than Grade 5 of the underground space are within the safe range. After reaching the magnitude of less than Grade 6, it can be seen from the 75% quartile evaluation table that the underground space has lost stability. When the peak acceleration was 2.5 m/s2, the plastic zone did not reach 2.22% in the 75% quartile table. For the sake of safety, it is recommended that the worst-case index be used to determine the seismic grade when evaluating the magnitude, the maximum acceleration, the maximum interlayer displacement angle, and the plastic zone.Figure 17 The plastic zone with the artificial seismic wave peak acceleration with different accelerations: (a) 0.8 m/s2, (b)1.8 m/s2, (c)2.5 m/s2 and (d)2.96 m/s2.

Figure 17

Table 10 The seismic resistance evaluation of the artificial seismic wave in the Oya quarry.

Table 10Intensity	Peak	
Acceleration	

	Maximum	
interlayer	
displacement	
angle	

	Plastic	
state	

	Seismic	
resistance	
Sorting	

	
(m/s2)	(rad)	(%)	
5-	0.8	1.04×10−3	4.93×10−3	AA(OBE)	
5+	1.8	2.34×10−3	2.65x×10−1	A(SEE)	
6-	2.5	3.19×10−3	1.72	(B)Unstable	
>6-	2.96	3.30×10−3	45.07	(B)Unstable	

4.4 Discussion

In this study, the horizontal dynamic component/vertical dynamic component of the velocity of the Oya stone wall is 1.4 ∼ 1.95. This suggests that the shaking degree of the underground space structure of Oya tuff quarry in the horizontal direction is significantly greater than that in the vertical direction. This study showed that near-focus earthquakes often start with ups and downs, followed by a horizontal swing. Far-focus earthquakes have fewer up-and-down jolts but are mainly characterized by horizontal swing. Therefore, the rocking degree of the surrounding rock wall in the horizontal direction is significantly greater than that in the vertical direction. The horizontal swing aggravates the shear action of the residual pillar, the surrounding rock wall, and the roof and floor. The shear failure of underground structures should be considered in seismic Monitoring, damage inspection, and support reinforcement.

In simultaneous earthquakes, the shaking degree of the residual pillar was twice that of the wall. Here, the residual pillar was only constrained by two ends; the horizontal direction showed an elliptical movement trend, and the surrounding rock wall was constrained in three directions. During the earthquake, the residual pillar was relatively free and had more significant deformation than otherwise. Therefore, it may also be concluded that the intersection position of the roof and the residual pillar was more likely to be damaged than the intersection position of the surrounding rock wall and ceiling. As it is a rigid connection, it is vulnerable to shear damage. Thus, when supporting reinforcement, priority should be given to the intersection of the roof slab and the residual pillar, followed by the junction of the roof slab and the surrounding bedrock.

This study compares the spectrum of seismic waves and the attenuation of stress waves with the experimental values to illustrate to what extent the results obtained by this research method are limited. The analytical model used in this study is highly consistent, which makes the dynamic evaluation criteria derived from the model dynamic analysis more convincing. It is important to note that improving the model accuracy is still a long-term goal.

Since there is no suitable seismic solid record in the area close to the site, the impact of the largest earthquake expected to occur at the plate junction in the next 30 years (with the source located in offshore Ibaraki Prefecture) on the Oya tuff quarry underground space was simulated by creating artificial seismic waves using the trigonometric superposition method. At the same time, it was found that the Oya tuff quarry underground space would be destabilized. This work uses the Fourier spectrum of seismic records of the target site itself directly in the synthesis. This paper uses the Fourier spectrum of the seismic records of the target site itself directly in the synthesis. As a result, the input ground-motion time interval is more in line with the structure's dynamic characteristics. It can better meet the requirements of the most unfavorable conditions. During more than one month of field monitoring, analysis of the seismic records revealed extremely active seismic activity at the plate boundary (especially offshore Ibaraki Prefecture) and in the Tokyo area, verifying the speculation that strong earthquakes at the Okinawa plate boundary in Ibaraki Prefecture and strong direct-down type earthquakes in Tokyo will occur in the next 30 years. In the future, the impact on the Oya Tuff quarry underground space can be predicted according to the research method in this paper.

For IDA, the research method is mainly applied to reinforced concrete and other above-ground structures. According to the Chinese code, the interlayer displacement Angle is unstable at the limit value of about 0.1Image 110−2 rad. According to the Japanese code, the interlayer displacement angle is limited to less than 0.5Image 110−2 rad [30], but it can be relaxed to less than 0.83Image 110−2 rad according to the conditions [31]. Therefore, it was proved that the displacement angle of the underground pure rock structure was almost concordant with that of the aboveground pure rock structure. The most important thing when using the IDA method is to select suitable seismic records. In the discussion of the solid seismic record selection method based on the PEER-NGA database, Kun et al. [21] considered the search for appropriate seismic records according to the seismic, geological characteristics of specific projects to provide the optimal input ground shaking for carrying out incremental dynamic analysis. However, the data from the PEER-NGA and Ki-NET databases are seismic records observed at (near) the surface. Due to the site's amplification effect, the database's seismic intensity is much larger than the bedrock seismic intensity. In contrast, the seismometer in this study was installed at the most profound depth of 60 m below ground near the bedrock. The measured ground vibrations used for input not only conform to the structural properties but also avoid the bias of the results caused by the amplification effect of the site. It makes the management values used to establish the dynamic stability evaluation criteria more credible. Eventually, the damage caused by artificial seismic waves can be well evaluated. It can be considered that the results obtained by this research method are more advantageous. It should be noted that the limitations of this study include the fact that not all engineering projects have the conditions to set up seismometers to obtain seismic waves at the original site. However, this approach is still quite applicable for representative or influential structural facilities that require model reliability.

Based on the magnitude of the earthquake, the maximum acceleration, and the inter-story displacement angle in the dynamic stability evaluation benchmark we have established for assessing dynamic stability, it is possible to determine and develop a disaster warning system based on the specific values observed in real-time by the existing seismic monitoring system at the Oya Observatory.

5 Conclusions

This study comprehensively analyzed the seismic response characteristics of large underground spaces such as underground quarries. Direct simulation using seismic waves from the ground surface and bedrock was inaccurate and may lead to biased conclusions. By integrating actual seismic records and artificially synthesized seismic waves, a more accurate simulation of the seismic impact on underground structures was achieved, thereby enhancing the accuracy and reliability of the evaluation.

The research results indicate that seismic solid events significantly damage underground quarry spaces. The residual pillars exhibit elliptical motion in the horizontal direction, with much greater displacement than the surrounding rock walls, exacerbating shear forces and posing a severe threat to the stability of underground structures. Due to the high correlation between the model and the site, this method can help identify damaged locations quickly after seismic events, providing valuable reference for significant engineering projects.

It is recommended that a seismic monitoring system be established to capture real-time seismic data and combine dynamic analyses of artificial seismic waves and Incremental Dynamic Analysis (IDA) to promote seismic resistance evaluations of all large underground structures. Moreover, personalized seismic hazard warning systems and reinforcement schemes should be established based on site characteristics and engineering importance, especially in seismically active areas. Additionally, real-time assessment of the safety of critical structures using numerical twin technology is a future research direction.

However, this study has certain limitations due to restrictions on installing seismometers. Therefore, other means may be necessary for selecting seismic input parameters in practical applications. Nevertheless, this research provides new insights and methods for evaluating the seismic performance of underground structures, offering valuable reference for research and practice in related fields.

In conclusion, this study comprehensively evaluated the seismic performance of underground quarry spaces by integrating actual seismic records and artificially synthesized seismic waves. It provides an essential basis for stability assessments of similar structures and offers new perspectives and methods for research and application in earthquake engineering.

Funding

This research was supported by 10.13039/501100002241 JST SPRING under Grant No. JPMJSP2128 . A donation from Chuden Engineering Consultants Co. Ltd, Japan, and a subsidy from Utsunomiya City Hall for Oya area development.

CRediT authorship contribution statement

Chuantao Cheng: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Investigation, Conceptualization. Noguchi Shizuo: Writing – original draft, Visualization, Investigation, Data curation. Tumelo K.M. Dintwe: Writing – original draft, Software. Takafumi Seiki: Writing – review & editing, Validation, Supervision, Methodology, Investigation. Motoi Iwanami: Visualization, Validation, Supervision, Investigation.

Declaration of Competing Interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: CHENG CHUANTAO reports financial support was provided by 10.13039/501100002241 Japan Science and Technology Agency . Takafumi SEIKI reports financial support was provided by Chuden Engineering Consultants Co. Ltd, Japan. Takafumi SEIKI reports financial support was provided by Utsunomiya City Hall for Oya area development. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

All data to support the conclusions have been either provided or are otherwise publicly available.
==== Refs
References

1 Disaster Prevention Movement Disaster prevention movement http://www.bousai.go.jp/kohou/kouhoubousai/r03/101/news_12.html 2020
2 Remembrance of natural disasters: the fall of the valley stone mining site http://www.ooyakousya.o0o0.jp/sisutem/pamphlet-2021.pdf
3 Celebi M. Topographical and geological amplifications determined from strong-motion and aftershock records of the 3 March 1985 Chile earthquake Bull. Seismol. Soc. Am. 77 4 1987 1147 1167 10.1785/BSSA0770041147
4 Jin-Jun H. Li-Li X. Variation of earthquake ground motion with depth Acta Seismol. Sin. 18 2005 72 81 10.1016/0148-9062(87)92282-0
5 Chen H. Overview of the seismic input at dam sites in China Earthq. Sci. 35 5 2022 410 425 10.1016/j.eqs.2022.05.006
6 Mendoza M. Auvinet G. The Mexico earthquake of September 19, 1985—behavior of building foundations in Mexico city Earthq. Spectra 4 4 1988 835 853 10.1193/1.1585505
7 Zhang C. Jin F. Seismic safety evaluation of high concrete dams: part 1: state-of-the-art design and research - sciencedirect Seism. Saf. Eval. Concr. Dams 2013 67 78 10.1016/B978-0-12-408083-6.00003-9
8 Zucca M. Crespi P. Tropeano G. Simoncelli M. On the influence of shallow underground structures in the evaluation of the seismic signals Ingegn. Sism. 38 1 2021 23 35
9 Vamvatsikos D. Cornell C.A. Incremental dynamic analysis Earthq. Eng. Struct. Dyn. 31 3 2002 491 514 10.1002/eqe.141
10 Vamvatsikos D. Seismic Performance, Capacity and Reliability of Structures as Seen Through Incremental Dynamic Analysis 2002 Stanford University 10.1002/eqe.141
11 Court A. Simonen K. Webster M. Trusty W. Morris P. Linking next-generation performance-based seismic design criteria to environmental performance (atc-86 and atc-58) Structures Congress 2012 2012 922 928 10.1061/9780784412367.082
12 Pathak S. Watt S. Khennane A. Al-Deen S. Dynamical systems approach for the evaluation of seismic structural collapse and its integration into pbee framework Soil Dyn. Earthq. Eng. 135 2020 106184 10.1016/j.soildyn.2020.106184
13 Li J. Zhong Z. Wang S. Bi K. Hao H. Seismic fragility analysis of water supply pipelines retrofitted with corrosion-protection liner buried in non-uniform site Soil Dyn. Earthq. Eng. 176 2024 108333 10.1016/j.soildyn.2023.108333
14 Bolotin V. Statistical theory of the aseismic design of structures Proceedings of the 2nd World Conference on Earthquake Engineering, vol. 2 Tokyo 1960 1365 1374
15 Shinozuka M. Sato Y. Simulation of nonstationary random process J. Eng. Mech. Div. 93 1 1967 11 40 10.1061/JMCEA3.0000822
16 Scanlan R.H. Sachs K. Earthquake time histories and response spectra J. Eng. Mech. Div. 100 4 1974 635 655 10.1061/JMCEA3.0001911
17 Nigam N. Phase properties of a class of random processes Earthq. Eng. Struct. Dyn. 10 5 1982 711 717 10.1002/eqe.4290100508
18 Yokoyama T. Theofanopoulos N. Watabe M. Distribution of phase differences in relation to the earthquake magnitude, distance to the fault and local soil conditions Proc. of 9th WCEE, Abstract 1 (A06-04) 1988 55
19 Boore D.M. Phase derivatives and simulation of strong ground motions Bull. Seismol. Soc. Am. 93 3 2003 1132 1143 10.1785/0120020196
20 Shrikhande M. Gupta V.K. On the characterisation of the phase spectrum for strong motion synthesis J. Earthq. Eng. 5 04 2001 465 482 10.1142/S1363246901000571
21 Kun J. Ruizhi W. Yefei R. Ground motion recordings selection for seismic design code J. Build. Struct. 38 12 2017 57 67 10.14006/j.jzjgxb.2017.12.007
22 Montejo L.A. Response spectral matching of horizontal ground motion components to an orientation-independent spectrum (rotdnn) Earthq. Spectra 37 2 2021 1127 1144 10.1177/8755293020970981
23 Fayaz J. Galasso C. A deep neural network framework for real-time on-site estimation of acceleration response spectra of seismic ground motions Comput.-Aided Civ. Infrastruct. Eng. 38 1 2023 87 103 10.1111/mice.12830
24 Lin X. Liu X. Hui J. Shan W. Assessment on detailed regional seismic damage risk of buildings based on time-history dynamic analyses Bull. Earthq. Eng. 2024 1 21 10.1007/s10518-024-01883-3
25 Kim T. Kwon O.-S. Song J. Deep learning-based response spectrum analysis method for building structures Earthq. Eng. Struct. Dyn. 2024 10.1002/eqe.4086
26 T. Dintwe, T. Seiki, O. Aydan, N. Tokashiki, H. Shimada, Experimental Study on the Effects of Large Single-Joints on Oya Tuff Pillar Strength All Days, 2019, ISRM–14CONGRESS–2019–240.
27 Shimizu C. Examination Report on Structural Stability of Haruhiko Confectionery Underground Factory 1989 Shimizu Construction Publisher Tokyo
28 Society J.G. Methods of Determining Design Soil Constants for Rock Foundations 2007 Japanese Geotechnical Society
29 Iai S. International standard (iso) on seismic actions for designing geotechnical works–an overview Soil Dyn. Earthq. Eng. 25 7–10 2005 605 615
30 China GB50011-2010, Pdf, URL http://www.jianbiaoku.com/webarbs/book/276/4455736.shtml
31 The building standard act of Japan Articles, URL http://hourei.net/law/325AC0000000201
