Visualizing water seepage dynamics in grotto relics via atom-based representative model

Water seepage in grotto relics, i.e., Yungang Grottoes, Dazu Rock Carvings, is a key issue to accurately describe the deterioration and weathering process of grotto rock mass. Considering rainfall infiltration, Finite element simulation was performed for studying the water flow through macro-channel of fractured rock in the 4th cave of Yungang Grottoes, where a group of joints with directions of S62°E and N5°W are widely developed. A 3D atom-based representative model was derived from X-ray diffraction (XRD) patterns and the related semi-quantitative calculation of grotto rock powders, for visualizing the associated seepage characteristics through micro-channel by means of molecular dynamics simulation, for the first time. By analyzing various properties, ranging from the configuration and energetic behaviors to the dynamic characteristics, the calculated water flux and mass flow rate were equal to 270 ns−1 and 8.10 × 10–12 g s−1, respectively. A dynamic process of water transport from the entrance region to the exit region was examined and it is consistent with the relative concentration profiles at the corresponding stage. The tagged O atoms experienced a zigzag movement instead of linear motion as expected, roughly exhibited the same target direction. The seepage characteristics in grotto relics experienced a complex evolution process and three types can be summarized: water infiltrates through micro-channels with a low flow rate; it flows through fracture with a relatively high flow rate; it turned into a kind of analogous pipe flow in inter-connected fracture network, resulting in water seepage hazard. Current simulation studies provide helpful insights for understanding the water flow-infiltration behavior of fractured rock in grotto relics.


Introduction
Grotto relics in China are widely distributed with the unique artistic features and historical value. As a global tangible cultural heritage, Yungang Grottoes and Dazu Rock Carvings have been inscribed into the World Heritage List by the United Nations Educational, Scientific and Cultural Organization (UNESCO) [1][2][3]. However, in the past few decades, the Grotto cultural relics have experienced different degrees of damage due to the longterm adverse effects of natural forces, human influences and environmental erosion, among which water seepage and weathering were responsible for the terrible deterioration and greatly threaten to the conservation of grotto relics in China [4][5][6]. Affected by complex geological structures and hydrological conditions, they have been a serious challenge. Previous studies indicated that water is one of the main causes of weathering to the grottoes and the phenomena such as swelling, disintegration and softening occurred after the water penetrated and interacted with the rock, resulting in the reduction of stability. The various forms of occurrence and special action *Correspondence: Junxia Wang jxwang@whrsm.ac.cn State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China of slight seepage produced a cyclical effect to the Grotto rock mass, and turned to be a dominant factor, accelerating the deterioration and weathering destruction of rock mass [7][8][9][10][11]. The existence of fracture network usually functions as an underground water flow path, thus posing a serious threat to safety of grotto relics. Therefore, the formation and mechanism of water seepage of fractured rock in grotto relics is a key issue to accurately describe the erosion and weathering process.
Fluid flow through macro-scale channels, i.e., faults, joints, fracture and so on, is crucial for evaluating the stability of geotechnical engineering and the mechanical behaviors in rock mass because they are recognized as the main channel for water migration. The finite element (FE) method and finite difference (FD) method were commonly used methods in solving complex problems in seepage mechanics and engineering, evolving from the seepage path, seepage process to hydraulic fracturing and coupling implementation [12][13][14][15][16]. The FE and FD method proved successful in evaluating water seepage from macro-scale, but there is little relevant research associated with the water seepage of fractured rock in grotto relics. Furthermore, the flow process among micro-channels associated with H 2 O molecules migration was ignored.
Fluid flow through micro-scale channels (small pores), driven by external fields, is critical to many fields of application, i.e., membrane separation, drug delivery, sensors, fuel cells, gas separation, desalination and so on [17][18][19][20][21][22]. Molecular dynamics (MD) simulation has been proved to be an efficient technique to simulate microscopic motions of molecules based on the known physics and the interactions between the molecules [23][24][25]. Using molecular dynamics approach, it was able to tackle the question of the interstitial fluid flow at micro-scale and describe properly the confined transport of H 2 O molecules through the channels. Pressure difference from either an atomic pressure or hydrostatic pressure gradient, is generated in MD simulation by applying an external constant force on each atom of the liquid flow and has been used to derive liquid transport in many materials, such as carbon nanotubes, graphenes, zeolites, metaorganic frames, polymers and others [26][27][28]. These reports reveal that the pressure difference, pore shape, pore size, temperature, hydrogen-bond (H-bond) and orientation play a vital role in diffusions and ordering of H 2 O molecules. However, the structure and dynamics of confined water from bulk phase going through the seepage channel of fractured rock in grotto relics obtained from MD simulations has never been reported. The distribution, orientation and motion trajectories of H 2 O molecules passing through the seepage channel are still unclear.
Because of the particularity of cultural relics, many non-destructive and in-situ methods are restricted, except for field observations [3]. Numerical modeling is the first choice for quantifying water seepage in grotto relics. In this study, we will be confining ourselves to perform FE and MD simulation targeting at visualizing the water seepage through macro-channel and micro-channel related to the direct infiltration of natural rainfall, in which an atomistic configuration of parallel feldsparquartz platelets mimicking skeleton micro-channel was derived and then the Lennard-Jones(LJ) potentials and Einstein relations was employ, thus to provide a better understanding of the driving processes of internal water transport in grotto relics and a theoretical basis for the conservation of the World Cultural Heritage.

Geological background
Yungang Grottoes (113°20‵E, 40°04‵N) is located at Datong city of Shanxi Province, China, as displayed in Fig. 1. The grottoes are known as one of the largest grotto groups in China and the world-famous stone carving art, including 45 major caves, 252 shrines and approximately 51,000 sculptures. There is an average of 423 mm rain per year with rainfall distribution mostly in July, August and September. The annual average evaporation is 1748 mm, with the maximum of 801 mm in June.
The whole grottoes are divided into three parts: the east (1st-4th Cave), the middle (5th-20th Cave) and the west (21th-53th Cave). The 4th Cave with typical features of seepage hazards ( Fig. 2a-f ) was chosen as example of seepage channel survey where the blue box shows the actual observed water seepage and the bottom of the grottoes is 10 m higher than the groundwater level, completely in the aeration zone. Precipitation is Fig. 1 The geological location of Yungang Grottoes in China the only source of replenishment for groundwater, which continues seeping via fractures until arriving in groundwater. Tensile joints were the most developed joints in the east of Yungang Grottoes, appearing in groups, with dip direction nearly erect and extension direction nearly east-west, roughly parallel to the faults. A total of 97 joints inside the 4th Cave were measured and the stereographic diagram for the joint patterns analysis was displayed in Fig. 2g, using upper hemisphere projection, as evidence from the rose diagram for the joint patterns analysis in Fig. 2h. The pole of each facet was shown by a triangle and the average plane was drawn with a circle line. Two groups of the joints in the 4th Cave were obvious: one group was 208°∠68° with direction of S62°E and the other is 265°∠75°with direction of N5°W。

Measurement of chemical composition
The grotto rock powders were prepared and measured using a X-ray diffraction detector (XRD, D8 Advance) with a scanning range from 3 to 70°. The XRD patterns of the grotto rock powders were analyzed in Fig. 3a for identifying the chemical composition existed in the grotto rocks. It is characteristic of obvious peaks at 20°, 26°, 36°, 39°, 50° and 60°, which could indexed to the standard lattice parameters of quartz. The typical diffraction peaks of feldspar were displayed at 28°, 35° and 42°. The peaks at 8°, 12°, 19° and 24° indicated the existence of biotite and kaolinite. It is mainly composed of 54.39% quartz and 31.20% feldspar according to semi-quantitative calculation as illustrated in Fig. 3b and as a result, feldspar-quartz slab was generated for the current computational studies.

Computational methods and models Finite element computational details
Various studies have shown that the critical Reynolds number (Re) for a flow between a Darcy and a non-Darcy flow is 10 and the flow satisfies Darcy`s law when Re ≤ 10 [29]. Re is defined as (Eq. 1) [30] (1) Re = ρQ µb where ρ is the density of water, Q is the flow rate, μ is the viscosity of water, b is the aperture of the fracture. In this paper, ρ is 1.0 × 10 3 kg m −3 , μ is 1.005 mPa s [30], Q is 1.0 × 10 -5 -2.0 × 10 −6 m 3 /s [31] and b is 0.001 m. So, the Re equals to 1.99-9.95 and it was governed by Darcy`s law.
The macro-scale FE model was conducted using commercial available ANSYS software as depicted in Fig. 4. The mesh was generated for all the volume entities with 10-node tetrahedral element to better mesh quality and the converging rate. As a consequence, the FE model was discretized using approximately 8078 tetrahedral  where v is the seepage rate, k is the permeability h is the total head,l is the length of seepage path.
The flow rate Q can be calculated according to Eq. 3 [32][33][34]: where w is the cross-section area.
The general governing differential equation for threedimensional seepage can be expressed mathematically as following (Eq. 4) [35]: where k x , k y and k z are the permeability in the X-, Y-and Z-direction, γ w is the specific weight of water and m w is storage curve slope.

MD computational details
The external pressure was applied to drive the fluid flow. All MD simulations are performed by Accerlery Materials Studio software. Figure 5 depicts the initial configuration of parallel feldspar-quartz slabs mimicking skeleton nanochannel with the inter-space (d) of 2.0 nm. A simple point charge-extended (SPC/E) model was applied for H 2 O molecules because of its excellent description for bulky water. 737 H 2 O molecules with density equal to = 1.0 g/cm −3 was randomly filled into the left reservoir with dimension of 2.00 × 2.51 × 4.40 nm 3 . The model consists of a 3 × 2 supercell of quartz slab and a 8 × 5 supercell of feldspar slab, for which position constraints were employed. One-layer Graphene sheet with size of 2.46 × 5.04 nm 2 is assigned to motion group and acts as moveable wall (piston) for creating the driven force toward H 2 O molecules in the left reservoir and the right reservoir (2.00 × 2.51 × 4.40 nm 3 ) was kept empty. Geometry optimization of 5,000 iterations was achieved for energy minimization. Following this, a MD simulation via the isothermal-isometric(NVT) ensemble with the Nose thermostat is performed using Drieding force field at 298 K with a timestep of 1 fs. The intermolecular potential energy includes the LJ 6-12 type potential (E vdW ) and the Coulomb potential (E ele ) to describe van der Waals and electrostatic interactions, respectively. Their mathematical expressions are displayed in Eq. 5 [36].

Results and discussions
The main source of the seepage water in the grottoes was the rainfall infiltration and seepage recharge through joint fissures or joint fractures, mainly depending on the rainfall intensity and permeability [37]. The seepage field calculated from FE method was displayed in the form of color-coded contours of the water head and flow rate distribution in response to the rainfall intensity of 100 mm d −1 (Fig. 6). Joints with a direction of S62°E were recognized to be an effective channel for water seepage. In the fracture zone, water infiltrates quickly from the top surface into the joint fracture with a high value of water heads, in accordance with the actual seepage  Fig. 2a-f, where the flow rate was ranging from 1.67 × 10 -5 to 4.67 × 10 -8 m 3 s −1 . The direct infiltration of natural rainfall acted as the driving force of water migration, providing a framework to construct an atomic model of the pore flow in the subsequent section.
The flow direction of water seepage is generally from the low potential area to the high potential area. H 2 O molecules suffered from a large energy barrier and they should overcome the large occupancy fluctuations to occupy the empty vacancies, which comes from larger energy barrier with strong interactions between water and channel. Figure 7a displayed the time course of energy variation evolving potential energy and its components during NVT simulation. It is found that they showed stable energy levels, a relatively strong electrostatic interaction (attractive) and a relatively weak van der Waal (vdW) interaction (repulsive), leading to a high potential energy barrier. Figure 7b illustrated the variation in the number of H 2 O molecules entering into seepage channel versus time. The water occupancy increased   [26,38], equal to 270 ns −1 and 8.10 × 10 -12 g s −1 , respectively. The magnitude of mass flow rate was close to the mass flow rate of graphenebased nanochannel with the pore width of 0.7-1 nm [28].
The flow behavior of H 2 O molecules in seepage channel involved three steps: (1) permeate into the seepage channel; (2) moving along the channel; (3) flow out of the channel. In order to probe the transport phenomena in detail, three stages (stage I, stage II and stage III) were defined as clarification of H 2 O molecules in the entrance region, the center region and the exit region, respectively. The typical feature at various stages was demonstrated  Fig. 8a, giving a dynamic process of water transport from the entrance region to the exit region. We divided the simulation box into many slabs along X-axis direction for statistical analysis. Hence, the relative concentration is given by the ratio of concentration in the slab to its average concentration across the entire system, as demonstrated in Fig. 8b. Significant discrepancies emerge at the number and magnitude of the peak values existing in the relative concentration distribution curve. Specifically, at stage I, H 2 O molecules are uniformly dispersed in the left reservoir with uniform distribution as verified by the multi-peaks. Part of H 2 O molecules entered into the channel because the relative concentration in the range of 2-4 nm changed from 0 to 2.0, as can be seen at stage II. When the whole channel was filled by H 2 O molecules, the relative concentration among the channel fluctuated around 1.0, indicating a uniform distribution at stage III.
These findings are consistent with the morphologies shown in Fig. 8a. The orientation distribution of O-H bonds in H 2 O molecules in a specified stage was also obtained to gain insight into ordering of the confined water, as plotted in Fig. 9. The orientation was defined as the angle between the flow direction (X-axis) and O-H bonds in H 2 O molecules. We have also included the corresponding orientation distribution at each stage for comparison. As can be seen, for each stage, there was a broader angular distribution and no ordered configuration, suggesting that the orientation behavior of confined H 2 O molecules responded to pressure-driven flow field without any disruption. The various distributions of may be derived from the fluctuation of H 2 O molecules and random thermal motion [39]. It was unfavorable for forming H-bonds among confined H 2 O molecules. A detailed discussion on the H-bond inspection is presented in the following section.
The aforementioned orientation behavior may be relevant to the H-bonds formation. When the H 2 O molecules orient themselves parallel to the flow direction, it is in favor of the formation of H-bonds among the confined H 2 O molecule. For a more detailed discussion, the number and structure of H-bonds for confined H 2 O molecules were characterized, as shown in Fig. 10. When the distance between the O atom of one H 2 O molecule, and the H atom of another H 2 O molecule was less than 3.0 Å, one can assume that a H-bond was formed [18,40]. Having a precise look at Fig. 10, it has experienced a nearly decreasing trend in the number of H-bonds with ignorance of some fluctuations. The H 2 O molecules should break some of H-bonds to enter the narrow channel due to steric crowding and large energy barrier, leaving a reduction in the number of H-bonds. This finding agrees well with the results of orientation distribution O-H bonds in H 2 O molecules. In addition, it is further interpreted as the strong repulse interaction between H 2 O molecules and feldspar because some reduction has been occurred may be attributed to the disappearance of the H-bonds between H 2 O molecules and feldspar surface (see Fig. 10c).
Besides analyzing the energy variation, water flux, relative concentration profiles, orientation distribution and H-bonds, further analysis is here focused on describing motion trajectories for elucidating the flow mechanism, by depicting the Cartier coordinates (x, y, z) of O atoms in H 2 O molecules located at different position along Z -axis direction. Figure 11a gives their initial positions of O n (n = 1, 2, 3). The moving trajectories of O 1 , O 2 and O 3 in X-, Y-and Z-axis directions were extracted and exhibited in Fig. 11b-d. What is worth mentioning is that the tagged O atoms experienced a zigzag movement instead of linear motion as expected, roughly exhibited the same target direction, moving towards the exit.
The seepage characteristics in grotto relics experienced a complex evolution process and three types can be summarized as water seepage propagation behavior in grotto relics (see Fig. 12): for type 1, the pores interconnected with other pores to form micro-channels and water infiltrates through micro-channels of fracture zone or the small original cracks in the upper area, where the flow rate is low. As a result of washout of atmospheric rainfall, water flows through fracture with a relatively high flow rate (type 2). With the expanding and connecting of fracture or cracks, water seepage maybe turned into a kind of analogous pipe flow, resulting in water seepage hazards (type 3).

Conclusions
In this study, a new attempt was made to employ FE and MD simulations for elucidating the water infiltration and seepage characteristics through the macro-channel and micro-channel of fractured rock in grotto relics in consideration of the rainfall infiltration. The feldsparquartz slab was generated for MD simulation because it is mainly composed of 54.39% quartz and 31.20% feldspar in grotto rock, as verified by XRD patterns. The water occupancy increased monotonously with simulation time and there were 270 H 2 O molecules passing through channel from the initial position along the flow direction, from which the water flux and mass flow rate were calculated. It is found that the broader angular distribution and disordered configuration existed in confined H 2 O molecules was unfavorable for forming H-bonds among confined H 2 O molecules, which was verified by the H-bond analysis. The motion trajectories of H 2 O molecules along Z-direction were investigated and they showed a zigzag movement instead of linear motion as expected. Three types can be summarized as water seepage propagation behavior in grotto relics: for type 1, the pores interconnected with other pores to form microchannels and water infiltrates through micro-channels in the upper area with a low flow rate; for type 2, water flows through fracture with a relatively high flow rate; for type 3, it turned into a kind of analogous pipe flow as the expanding and connecting of fracture or cracks, resulting in water seepage hazard. The seepage rate k The permeability h The total head l The length of seepage path Q Flow rate w The cross-section area k x , k y , k z The permeability in the X-, Y-and Z-direction γ w The specific weight of water m w Storage curve slope SPC/E Simple point charge-extended NVT Isothermal-isometric vdW Van der Waal