Deep Rock Behaviour in Engineering EnvironmentsView this Special Issue
Research Article | Open Access
Guangdi Deng, Heping Xie, Mingzhong Gao, Cong Li, Zhiqiang He, "Numerical Simulation on the Evolution of Mining-Induced Fracture Network in a Coal Seam and Its Overburden under the Top Coal Caving Method", Advances in Civil Engineering, vol. 2020, Article ID 8833193, 14 pages, 2020. https://doi.org/10.1155/2020/8833193
Numerical Simulation on the Evolution of Mining-Induced Fracture Network in a Coal Seam and Its Overburden under the Top Coal Caving Method
The evolution of the fracture network induced by mining has an important influence on the mechanical behavior of the rock and the safety of the mine. In this study, a new 3D modeling approach based on discrete fracture network (DFN) and fractal theory was developed and applied to the Tashan coal mine. The model results of the evolution of fracture in the coal seam are consistent with field observations. The evolution law of fracture network in the overburden strata was studied, and the results show that the fractal dimension of the fractures in the lower strata increased linearly and stabilized quickly within 50 m behind the mining face. In the higher strata, most of the fractures were generated behind the mining face and continued to develop farther than 100 m. More fractures were generated in the lower strata than in the higher strata, and the fractures more easily developed and expanded in the soft rock than in the hard rock. The evolution of the fractures of the main thick hard roof in the lower strata had a great impact on the generation of fractures in the higher strata.
The longwall top coal caving (LTCC) method is a new type of mining technique that is popular for thick coal seams in China due to its high efficiency and output . However, these large-scale excavations seriously disturb the original equilibrium state underground, and the mining-induced redistribution of stress leads to the initiation and propagation of fractures, creating fractured zones with high connectivity around the mining excavations. These fractured zones weaken the strength of the rock mass and provide a pathway for confined water and gas inflow, leading to the movement of overburden strata and water and gas outbursts, which severely affect mine safety [2–6]. It is an essential issue in mining science to recognize the evolution of fracture networks during mining processes. Many researchers have studied this topic and achieved fruitful results.
Xu and Qian [7, 8] studied mining-induced fractures in the overlying rock layer of a coal mine and revealed their “o”-shaped distribution in the longwall face. Li et al.  discussed the effect of key strata on the distribution of fractures and concluded that an elliptic zone would be formed in the overlying strata. Kang et al. investigated the mechanism of the occurrence of mining-induced fractures and analyzed the morphological characteristics of a fracture zone [10–12].
Nevertheless, it is difficult to develop potential prevention mechanisms and the associated spatiotemporal effects throughout the mining process from a quantitative perspective due to the extreme complexity and disorder of mining-induced fractures . The fractal geometry proposed by Mandelbrot  can quantitatively describe complex objects and irregular phenomena. Xie et al.  introduced the fractal geometry of fractures and discovered the self-similarity behavior of the spatial distribution of fracture networks. Since then, the evolution of fracture networks has been studied with fractal theory based on similar material models [16, 17], on-site field observations [18, 19], and 2D discrete element methods [20, 21]. However, field observations are limited by local conditions and observation tools. In similar material methods, the mining-induced fractures need to be determined manually from the model with subjectivity. 2D discrete element methods present the fracture evolution in only a single plane.
In this study, a new 3D modeling approach for the calculation of the fractal dimension of a fracture network was developed. This approach makes use of the built-in discrete fracture network (DFN) module in the commercial software FLAC3D to create stochastically distributed fractures. After each longwall face excavation, deterministic global stress analysis is first performed, followed by an evaluation of all the discrete fractures to identify fracture slip. A fracture is deemed to “generate” if it slips. Finally, the fractal dimension of the generated fracture network was calculated based on specified section planes. With the advance of mining, the evolution of the mining-induced fracture network can be determined.
This modeling approach was applied to the Tashan coal mine in Shanxi Province, China, where very thick coal seams were extracted using the LTCC method. Section 2 presents the detailed methodology of the numerical model. The establishment and running of the model based on the 8212 working panel are introduced in Section 3. The evolution of the mining-induced fracture network in the coal seam and overburden strata is discussed in Section 4.
2. Description of the Modeling Procedure and Method
The modeling procedure involves three main steps. First, a 3D model representing the working longwall panel and the strata is built with a group of discrete fractures (a DFN), and then mechanical properties and initial conditions are applied to zone elements to reach an initial equilibrium state. Next, progressive longwall mining is simulated by assigning null mechanical properties to zone elements. At the next excavation step, together with the simulation of coal extraction, the previously “removed” elements are assigned goaf material to simulate the recompaction of the goaf. In the final modeling step, a fracture stress calculation is performed to determine whether fracture slip occurs based on the prevailing stress state. Specific section planes are set up, and the fractal dimensions of the slipped fracture network on the section planes are calculated. Steps 2 and 3 are repeated in each coal excavation step until the completion of the longwall mining simulation in the study area. Notably, the model progress from Step 2 to Step 3 is a one-way process, which means that the state of the fractures is determined by the zone element stress states, while it is assumed that the stress states of the rock mass are not affected by fracture slippage.
2.1. DFN in FLAC3D
A randomly distributed network of disc-shaped discrete planar fractures can be created in the FLAC3D model by setting a series of attributes. These attributes include the fracture stopping criterion, fracture size distribution, and fracture orientation.
The total number of fractures was used as the stopping criterion to generate the fracture network. In consideration of the model size and computing ability, the total number of fractures was set to twice the total number of zone elements. The power law relationship [22, 23] has been adopted to describe the distribution of fracture sizes:where is the number of fractures with size per unit volume, is the density term, and a (>1) is the power law scaling exponent and usually ranges between 2 and 4.
The total number of fractures within the range [lmin, lmax] can be obtained by integrating equation (1):
The cumulative distribution function is then given by
In view of the meshed element size and the scale of research, the minimum fracture size was set to the shortest zone element length, and the maximum fracture size was set to 15 m.
The fracture position was assumed to follow a uniform spatial distribution throughout the model domain. Due to the lack of geological measurements, the dip angles and dip directions that determine the fracture orientations also follow a uniform distribution.
2.2. Fracture Slip Criterion
The Mohr–Coulomb slip criterion is widely used to evaluate fracture slippage. The fracture is considered to have slipped if the shear stress along the fracture plane exceeds its resistance :where is determined by the properties of the fracture surfaces aswhere σ is the normal stress on the fracture surface and and are the static friction coefficient and the cohesion along the fracture, respectively.
Alternatively, the fracture shear strength in equation (5) can be estimated by applying a strength reduction factor to the intact rock strength given by the following equation:
The shear and normal stresses used in equations (5) and (6) can be calculated using the stresses on the zone element closest to the centroid of a fracture. Salamon  used this method for fracture slip evaluation. This approach is valid if the fracture is small compared to a zone element. As mentioned above, the fracture is larger than a zone element, so following Board  and Cao et al. , this study adopted nine test points for slip, one testing point at the center, and eight testing points uniformly distributed around the periphery of the fracture surface. The testing point is deemed to slip if the zone element closest to it is in the failure state or if equation (6) is satisfied, and the fracture is thought to slip when all nine testing points have slipped. The workflow of the evaluation of fracture slip is presented in Figure 1.
2.3. Calculation of the Normal and Shear Stresses at a Test Point
The test points on a fracture surface are specified in local coordinates with their origin at (x0, y0, and z0) in the global coordinate system (Figure 2). The direction cosines of the axis OX′ are l1, m1, and n1, and those of axes OY′ and OZ′ are l2, m2, and n2, and l3, m3, and n3, respectively. The transformation of the ith test point from its local coordinates to the global coordinates is given by
The normal stresses along the local coordinate axis directions at the test point can be written aswhere , , and and , , and are the normal and shear stresses of the zone element closest to the test point, respectively. The corresponding normal and shear stresses on the fracture surface associated with the test point are then given by
2.4. Calculation of the Fractal Dimension of the Fracture Network
The calculation of fractal dimension on the section plane is done by the “counting the boxes” method , and the diagram of the calculation is shown in Figure 3. A square mesh of scale R is used to cover the target area, and the number of grid cells in which the fracture length is equal to or greater than the corresponding grid size is counted and denoted by N (R). An array of corresponding N (R) values can be obtained with a series of R values to study the changes in log (R) and log (1/R). The fractal dimension D is defined as follows:
3. Establishment of the Numerical Model
3.1. Introduction of 8212 Working Facing
The Tashan coal mine is one of the largest underground coal mines in the world and produced more than 20 Mt of coal in 2019. The 8212 working face is mined using the LTCC method and is located in the northwest region of the second mining district 2614 m along strike and 230 m along dip. The working face mines the very thick #3–5 coal seam, and the average thickness of this coal seam is 11.6 m. The dip angle of the coal seam is 4°; the coal seam can be viewed as an approximately horizontal structure. The average mining height and caving height are 3.5 m and 8.1 m, respectively. The main roof strata in 8212 working face are sandstone which would accumulate large elastic energy before failure [27, 28], and the breaking of roof strata can cause dynamic disasters such as rockburst [29–31] and coal-gas outburst [32, 33]. The detailed rock types and thicknesses of the 8212 working panel strata are presented in Table 1.
3.2. Model Setup and DFN Generation in FLAC3D
A 3D numerical model based on the 8212 working face was constructed in FLAC3D. The size of the model was 600 m × 300 m × 500 m in the x, y, and z directions, respectively, and the mining area was 300 m × 230 m. The model comprised 620400 zone elements and 643720 grid-points. As discussed in Section 2.1, the number of fractures was set to twice the number of zone elements, i.e., 1240800. As the element size of the coal and roof strata was approximately 5 m × 5 m × 3 m, the fracture sizes, following a power law distribution, ranged from 3 m to 15 m, and the numerical model and the size of the generated fracture are shown in Figure 4. Notably, this paper focuses on the evolution of the fractures in only the coal seam and overburden strata, so discrete fractures were not generated in the floor strata.
3.3. Property and Initial Conditions
The Mohr-Coulomb model was used to represent shear failure in the coal and strata. The detailed physical and mechanical properties of the strata at the 8212 working panel are presented in Table 2 . The falling waste in the gob was considered an elastic material, and as the working face advances, the waste was progressively compressed by the overburden. Experimental results show that the density, Young’s modulus, and Poisson’s ratio can be characterized as follows :
In the model, the horizontal displacement of the two sides was limited, and the vertical displacement of the bottom was limited. A vertical load () was applied on the top to simulate the overburden weight.
3.4. Mining Excavation Procedure
After the initial equilibrium state, model mining started with the excavation of the intake airway and return airway, and thirty excavation steps were conducted to simulate progressive longwall mining. The excavation step is 10 m in the mining direction, representing two days of face advance. Each excavation also included the “backfill” of the previously removed elements by assigning waste material properties. A special subroutine was invoked after each step to determine and record the slippage of the fractures, and a fractal dimension calculation was implemented.
4. Model Results
4.1. Section Plane Setup
To describe the evolution of the fracture network with the mining process, eleven section planes are considered in the model, as shown in Figure 5. The first four vertical section planes are along the X direction (face advance direction) with coordinates of 250, 300, 350, and 400. Note that the coordinate in the X direction for the mining start line is 150, so as the mining face advances, and these vertical sections can record the evolution of the fracture network in front of and behind the mining face. Another seven horizontal section planes are placed at the center of each stratum to present the evolution of the fracture network in each stratum during the mining process.
4.2. Evolution of the Coal Seam Fracture Network
As the mining face advances, the redistribution of the stress in the coal seam leads to the initiation and propagation of fractures. Figure 6 shows the evolution of fractures in the coal seam in section plane 3. Fractures continue to grow with the approach of the mining face. To quantitatively describe the evolution of the fractures in the coal seam, the fractal dimensions of the fractures in the section plane were calculated based on the method mentioned above, as shown in Figure 7. The fractures exhibited a two-stage evolution as the mining face progresses: a stable growth stage (at a distance of 170∼60 m from the mining face) and a rapid growth stage (at a distance of 0∼60 m from the mining face). In the first stage, fractures were generated under the influence of dynamic mining, and the number of fractures gradually increased as the mining face progressed. In the second stage, as the mining face continued to approach, the impact of mining became stronger, and the number of fractures increased rapidly. At approximately 10 m in front of the mining face, most of the fractures developed and interconnected due to failure of the coal body.
Field experiments, including borehole observation and abutment pressure monitoring, were conducted by Gao et al.  at the 8212 working panel, and the results are shown in Figure 8. The on-site field data present great fluctuations due to the changing geological environment. However, these data reveal a two-stage fracture evolution, and the abutment pressure also shows that the coal rock would break and fail at approximately 10 m in front of the mining face, leading to a decrease in abutment pressure. In general, the model results of the fracture evolution in the coal seam are consistent with the field observations.
4.3. Evolution of the Fracture Network in the Overburden Strata during the Mining Process
Figure 9 shows the relationship between the fracture network and distance from the mining face in section plane 2. The different colors represent different strata. The fractures continue to develop and grow steadily with the advance of the mining face. The fractal dimension of the fractures in the vertical section plane with distance from the mining face is obtained to quantitatively analyze the fracture evolution, as shown in Figure 10. The development of the fractal dimension of fractures with distance from the mining face can be divided into two stages: a linear growth stage (in front of the mining face) and a slow growth stage (behind the mining face). In the first stage, as the mining face moved forward, the number of fractures increased steadily due to the influence of mining, and the fractal dimension increased linearly with the mining advance. Notably, these fractures were generated mainly in the coal seam and lower strata. In the second stage, the fractures in the lower strata continued to develop, leading to the deformation and movement of the lower strata. Meanwhile, fractures in the higher strata continued to develop, presenting a saddle shape in the vertical section plane. Therefore, the overall fractal dimension of the fractures in the second stage slowly increases and gradually stabilizes behind the mining face.
4.4. Evolution of the Fracture Network in the Overburden Strata in the Vertical Section Planes
Section 4.3 discusses the overall fracture network evolution with distance from the mining face. In this section, a detailed study of each overburden stratum is conducted. Taken from the vertical section planes, the development of the fracture network in each stratum with distance from the mining face can be obtained, and the fractal dimensions were calculated, as shown in Figure 11. The evolution of the fractures in each stratum is related to its location.
In the lower strata such as the pebbly sandstone and siltstone, the evolution of the fractures is directly influenced by the mining excavation. The fractal dimensions of the fractures in the different section planes present similar characteristics. In front of the mining face, the number of fractures increased as the mining face advanced, and the fractal dimension of the fractures exhibited a corresponding linear increase. Behind the mining face, the strata collapsed rapidly and fill the goaf. Most of the fractures in the strata were well developed around the mining face, and the fractal dimension of the fractures stabilized quickly within 50 m behind the mining face.
The higher strata were not affected at the beginning of the mining process until the impact of the mining extended upward. Therefore, the evolutions of the fractures in each section plane are different. In addition, in contrast to the lower strata, in the higher strata, only a few fractures were generated in front of the mining face, and most of the fractures continued to develop more than 100 m behind the mining face.
4.5. Evolution of the Fracture Network in the Overburden Strata in the Horizontal Section Planes
Figures 12 and 13 show the evolution of the fracture network in the typical lower rock strata (pebbly sandstone) and higher rock strata (medium sandstone) with mining advance based on the horizontal section planes. Fractures were generated around the mining excavations as the mining advanced, exhibiting an elliptical distribution in the strata. To quantitatively compare the evolution of the fracture network in each stratum, the fractal dimension of the fracture network was calculated based on the horizontal section planes, as shown in Figure 14. It should be note that the height “h” in the figure means the distance from the section planes to the coal seam. Two distinct features were presented from the figure:(1)In terms of the number of fractures, more fractures were generated in the lower strata than in the higher strata. In addition, a fracture could more easily develop and expand in the soft rock (mudstone) than in the hard rock (fine sandstone).(2)The higher strata (medium sandstone) were almost unaffected in the first 100 m of mining until massive fractures were generated in the main thick hard roof (siltstone), which shows that the evolution of fractures of the main thick hard roof in the lower strata had an important influence on the generation of fractures in the higher strata.
It is an essential issue in mining science to recognize the evolution of fracture networks during mining. In this paper, a new 3D modeling approach based on a DFN in FLAC3D and fractal dimension was developed and applied to the Tashan coal mine. The model results of the fracture evolution in the coal seam are consistent with field experiments. Finally, the evolution of a fracture network in overburden strata was studied, and the following conclusions were drawn:(1)The fractures in the coal seam exhibited a two-stage evolution as the mining face progresses: a stable growth stage and a rapid growth stage.(2)In the lower strata, the fractal dimension of the fractures linearly increased as the mining face moved forward and stabilized quickly within 50 m behind the mining face.(3)In the higher strata, only a few fractures were generated in front of the mining face, and most of the fractures in the higher strata continued to develop further than 100 m behind the mining face.(4)More fractures were generated in the lower strata than in the higher strata, and the fractures were more easily developed and expanded in the soft rock than in the hard rock. In addition, the evolution of the fractures of the main thick hard roof in the lower strata had an important influence on the generation of fractures in the higher strata.
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
The authors would like to acknowledge the support by the National Natural Science Foundation of China Youth Fund (grant no. 51822403).
- H. Xie, Z. Chen, and J. Wang, “Three-dimensional numerical analysis of deformation and failure during top coal caving,” International Journal of Rock Mechanics and Mining Sciences, vol. 36, no. 5, pp. 651–658, 1999.
- S. G. Li, M. G. Qian, and P. W. Shi, “Study on bed-separated fissures of overlying stratum and interstice permeability in fully-mechanized top coal caving,” Chinese Journal of Rock Mechanics and Engineering, vol. 19, no. 5, pp. 604–607, 2000.
- H. W. Zhou, T. Zhang, D. J. Xue et al., “Evolution of mining-crack network in overburden strata of longwall face,” Journal of China Coal Society, vol. 36, no. 6, pp. 1957–1962, 2012.
- M. Z. Gao, W. C. Jin, Z. X. Dai et al., “Relevance between abutment pressure and fractal dimension of crack network induced by mining,” International Journal of Mining Science and Technology, vol. 6, pp. 148–153, 2013.
- M. Z. Gao, S. Zhang, J. Li et al., “The dynamic failure mechanism of coal and gas outbursts and response mechanism of support structure,” Thermal Science, vol. 23, p. 122, 2019.
- M. Z. Gao, Z. L. Zhang, X. G. Yin et al., “The location optimum and permeability-enhancing effect of a low-level shield rock roadway,” Rock Mechanics & Rock Engineering, vol. 51, 2018.
- J. L. Xu, “Study of features of distribution of overlaying stratum fractures caused by mining,” Ground Pressure and Strata Control, vol. 14, no. 3/4, pp. 210–212, 1997.
- M. G. Qian, “Study on the “o shape” circle distribution characteristics of mining induced fractures in the overlying strata,” Journal of China Coal Society, vol. 23, no. 5, pp. 466–469, 1998.
- S. G. Li, M. G. Qian, and P. W. Shi, “Study on the dynamic distribution characteristic of fissures elliptic paraboloid zone in the mining overlying strata,” Mine Pressure and Roof Management, vol. 3-4, pp. 44–46, 1999.
- Y. H. Kang and G. X. Zhao, “The effect of overburden strata characteristics on heights of the fractured and caved zones,” Mine Coal Mine Technology, vol. 29, no. 1, pp. 52–54, 1998.
- J. L. Xu, M. G. Qian, and H. W. Jin, “Study on coal and coal-bed methane simultaneous extraction technique on the basis of strata movement,” Journal of China Coal Society, vol. 29, no. 2, pp. 129–132, 2004.
- G. M. Cheng, K. Huang, and X. J. Chen, “The stress distribution characteristics of top coal and its influence on permeability,” Coal Engineer, vol. 10, pp. 35–37, 2002.
- M. Z. Gao, R. Zhang, J. Xie, G. Y. Peng, B. Yu, and P. G. Ranjith, “Field experiments on fracture evolution and correlations between connectivity and abutment pressure under top coal caving conditions,” International Journal of Rock Mechanics and Mining Sciences, vol. 111, pp. 84–93, 2018.
- B. B. Mandelbrot, The Fractal Geometry of Nature, W. H. Freeman, San Francisco, CA, USA, 1982.
- H. P. Xie, G. M. Yu, L. Yang et al., “Research on the fractal effects of crack network in overburden rock stratum,” Chinese Journal of Rock Mechanics and Engineering, vol. 18, no. 2, pp. 147–151, 1999.
- L. Tao, X. L. Liu, S. J. Wang, E. Wang, and Q. Li, “Study on the fractal characteristics of fracture network evolution induced by mining,” Advances in Civil Engineering, vol. 2018, Article ID 9589364, 13 pages, 2018.
- G. Yu, H. Xie, J. Zhoa, and L. Yang, “Fractal evolution of a crack network in overburden rock rtratum,” Discrete Dynamics in Nature and Society, vol. 5, no. 1, pp. 47–52, 2000.
- Z. G. Wang, H. W. Zhou, and H. P. Xie, “Research on fractal characterization of mined crack network evolution in overburden rock stratum under deep mining,” Rock & Soil Mechanics, vol. 30, no. 8, pp. 2403–2408, 2009.
- M. Z. Gao, W. C. Jin, C. J. Zheng et al., “Real time evolution and connectivity of mined crack network,” Journal of China Coal Society, vol. 37, no. 9, pp. 1535–1540, 2012.
- X. D. Zhang, Z. R. Xu, Z. J. Su et al., “Research of the imitation about the fractal crack in overburden rock stratum by computer,” Chinese Journal of Rock Mechanics and Engineering, vol. 20, no. 6, pp. 809–812, 2001.
- J. Y. Feng, X. H. Liu, and Z. Q. Yu, “Numerical simulation study on the mining-induced fracture evolution of steep coal seam,” Journal of China Coal Society, vol. 42, no. 8, pp. 1971–1978, 2017.
- W. Cao, J.-Q. Shi, G. Si, S. Durucan, and A. Korre, “Numerical modelling of microseismicity associated with longwall coal mining,” International Journal of Coal Geology, vol. 193, pp. 30–45, 2018.
- E. Bonnet, O. Bour, N. E. Odling et al., “Scaling of fracture systems in geological media,” Reviews of Geophysics, vol. 39, no. 3, p. 347, 2001.
- M. D. G. Salamon, “Some applications of geomechanical modelling in rockburst and related research,” in Proceedings of the 3rd International Symposium on Rockbursts and Seismicity in Mines, pp. 297–309, Kingston, OT, USA, August 1993.
- M. Board, “Numerical examination of mining-induced seismicity,” International Society for Rock Mechanics and Rock Engineering, 1996.
- G. M. Yu, H. P. Xie, H. W. Zhou et al., “Distribution law of cracks in structured rock mass and its experimental investigation,” Journal of Experimental Mechanics, vol. 13, no. 2, pp. 145–154, 1998.
- H. Xie, L. Li, R. Peng, and Y. Ju, “Energy analysis and criteria for structural failure of rocks,” Journal of Rock Mechanics and Geotechnical Engineering, vol. 1, no. 1, pp. 11–20, 2009.
- F. Q. Gong, J. Y. Yan, S. Luo et al., “Investigation on the linear energy storage and dissipation laws of rock materials under uniaxial compression,” Rock Mechanics and Rock Engineering, vol. 52, no. 12, 2019.
- X. Si and F. Gong, “Strength-weakening effect and shear-tension failure mode transformation mechanism of rockburst for fine-grained granite under triaxial unloading compression,” International Journal of Rock Mechanics and Mining Sciences, vol. 131, Article ID 104347, 2020.
- Q. S. Wu, J. Q. Jiang, Q. L. Wu et al., “Study on the fracture of hard and thick sandstone and the distribution characteristics of microseismic activity,” Geotechnical and Geological Engineering, vol. 36, 2018.
- J.-A. Wang and H. D. Park, “Comprehensive prediction of rockburst based on analysis of strain energy in rocks,” Tunnelling and Underground Space Technology, vol. 16, no. 1, pp. 49–57, 2001.
- T. Liu, B. Q. Lin, W. Xiao et al., “A safe mining approach for deep outburst coal seam groups with hard-thick sandstone roof: stepwise risk control based on gas diversion and extraction,” Energy Science and Engineering, vol. 8, 2020.
- J. G. Fu, W. Wu, and G. Y. Hua, “Gas control technology for fully mechanized face of outburst coal seam with hard roof during initial caving,” Advanced Materials Research, vol. 347–353, pp. 974–978, 2012.
- W. Wang, Y.-p. Cheng, H.-f. Wang, W. Li, and L. Wang, “Coupled disaster-causing mechanisms of strata pressure behavior and abnormal gas emissions in underground coal extraction,” Environmental Earth Sciences, vol. 74, no. 9, pp. 6717–6735, 2015.
Copyright © 2020 Guangdi Deng et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.