A Numerical Simulation Method Of Karst Fracture Network Evolution In Laminar-turbulent Coexistence Flow Field

Karst is widely distributed in China,and its types are complex and varied,accounting for more than one third of the land area in China.The karst forms such as fissure,pipelines and caves developed in karst geologic bodies not only provide space for groundwater and mineral resources,but also are the places where karst collapse and water inrush occur.Although many domestic scholars have studied the karst aquifer by numerical simulation,there are few studies on the evolution of karst fracture network in laminar-turbulent coexistence flow field.With the continuous dissolution process,the width of the fracture is widened continuously,and the flow pattern of groundwater flow gradually transits from laminar flow to turbulent flow.When the groundwater flow in the fracture is turbulent,the water head distribution,flow rate and the rate of fracture widening in the fracture network are quite different from those in the laminar flow.In this study,a numerical simulation method of karst fracture network evolution in laminar-turbulent coexistence flow field is proposed by using finite difference numerical simulation technology.The software FractureToKarst which was developed by our research group is improved.A turbulent module is added to the water head calculation part by using C++ language,and the function of the software to simulate the evolution of karst fracture network in laminar-turbulent coexistence flow field is realized.Based on the assumption of two-dimensional smooth parallel plate model,a two-dimensional fracture network seepage model of karst aquifer system is constructed by discontinuous medium analysis method.Fracture network is discretized and seepage model is coupled with dissolution model.Firstly,water balance principle is used to establish the equations of the nodal water head of the fracture network.When the flow in a fracture is laminar,the flow rate is expressed by the cubic law.When the flow in a fracture is turbulent,the flow rate is expressed by the turbulence formula proposed by Lomize G(1951),and the water head value of each node and the flow rate in each fracture are calculated by iteration method.Then,the dissolution rate of each fracture is calculated according to the empirical formula of carbonate dissolution rate proposed by Dreybrodt W(1990,1996).According to the principle of mass conservation,the fracture dissolution enlargement in this time step is calculated,and finally the fracture width after dissolution is obtained.After calculating the new fracture width,the above steps can be repeated,the water head distribution can be recalculated,and the fracture width after the next time step can be obtained,thus the evolution process of karst fracture network can be obtained.In this study,the turbulence models of confined karst aquifer and phreatic karst aquifer were designed to verify the correctness of software simulation,and the evolution process of fracture network in two models of karst aquifer overflow spring and paleokarst basin was simulated.It is found that neither the overflow spring model nor the paleokarst basin model,which only receives rainfall recharge,nor turbulence occurs until the watershed disappears.The overflow spring model receiving river and rainfall recharge first occurred turbulence at the discharge point of the spring in 1200 years,and the turbulence gradually developed to the internal system with the evolution time.In the whole process of karst system evolution,karst development is the strongest near the phreatic water level.