Gas hydrate saturations estimated from pore-and fracture-filling gas hydrate reservoirs in the Qilian Mountain permafrost, China

Accurate calculation of gas hydrate saturation is an important aspect of gas hydrate resource evaluation. The effective medium theory (EMT model), the velocity model based on two-phase medium theory (TPT model), and the two component laminated media model (TCLM model), are adopted to investigate the characteristics of acoustic velocity and gas hydrate saturation of pore- and fracture-filling reservoirs in the Qilian Mountain permafrost, China. The compressional wave (P-wave) velocity simulated by the EMT model is more consistent with actual log data than the TPT model in the pore-filling reservoir. The range of the gas hydrate saturation of the typical pore-filling reservoir in hole DKXX-13 is 13.0~85.0%, and the average value of the gas hydrate saturation is 61.9%, which is in accordance with the results by the standard Archie equation and actual core test. The P-wave phase velocity simulated by the TCLM model can be transformed directly into the P-wave transverse velocity in a fracture-filling reservoir. The range of the gas hydrate saturation of the typical fracture-filling reservoir in hole DKXX-19 is 14.1~89.9%, and the average value of the gas hydrate saturation is 69.4%, which is in accordance with actual core test results.

In the sample of gas hydrate from the QMP, there are two obvious different filling patterns: pore-filling and fracture-filling 3,4 . The gas hydrate displaces the fluid within the pores of the rock matrix in the pore space and is indistinguishable to the naked eye. The gas hydrate that is produced in the cracks does not occupy the pore space, but forces the formation rock to open and form cracks and fill them with the visible white ice layer 14,16 . For the different types of gas hydrate reservoir, adopting a velocity model of an isotropic porous medium for the inversion of gas hydrate saturation will lead to large errors in the resource estimations. Logging must usefully identify gas hydrate reservoirs and calculate gas hydrate saturation, and also to provide a reference for regional gas hydrate evaluation. To do this, an analysis of the acoustic velocity characteristics must be made for the gas hydrate saturation estimation within the different types of reservoirs present in the study area.
In this paper, the reservoir characteristics of gas hydrate in the study area are analyzed using ultrasonic imaging logging and drilling core data. The acoustic velocity characteristics of pore filling gas hydrate reservoir in the QMP are simulated using the forward modeling velocity model that is established by the effective medium theory (EMT) model, and the simulation results are compared with the elastic wave velocity model based on the two-phase medium theory (TPT) model. The acoustic velocity characteristics of fracture-filling gas hydrate reservoir in the QMP are also simulated using the two component laminated media model (TCLM). Finally, the acoustic velocities that are simulated by the pore-and fracture-filling reservoirs are used to estimate the gas hydrate saturation in the study area.

Geological characteristics of gas hydrate
The Qilian Mountain permafrost is located on the northern margin of the Tibetan Plateau, and the permafrost area is about 10 × 10 4 km 2 . In recent years, the gas hydrate samples are obtained firstly by drilling in the Juhugeng mining of the Muli coal field 3,4 . The drilling area is tectonically situated in the western Middle Qilian block formed during the Caledonian Movement, adjacent to the South Qilian structural zone. Except for the Quaternary deposits, the exposure strata mainly include Jiangcang Formation (J2j) and Muli Formation (J2m) of middle Jurassic age. The lithology of the drilling formation is mainly mudstone, siltstone, fine sandstone, medium sandstone, and coarse sandstone. Drilling test wells of DK-1, DK-2, Dk-3, DK-4, DK-5, DK-6, DK-7 and DK-8 were carried out, and the gas hydrate was acquired in multiple layers of the wells DK-1, DK-2, DK-3, DK-7 and DK-8 (Fig. 1). Gas hydrate was mainly found in the Jiangcang Formation of the middle Jurassic, and occurred within intervals of about 133-396 meters below the surface, below the permafrost in the drilling area.
The combined information from drilling and core experiment studies shows that gas hydrates mainly occur in the Jiangcang Formation of Middle Jurassic in the Muli permafrost. The geological and occurrence characteristics indicate that gas hydrates occur in separate reservoirs below the base of permafrost. In addition, gas hydrates are vertically distributed non-continuously in each borehole, and the nature of the gas hydrate distribution in the lateral areas between drill holes is not apparent due to the rock fracture system that plays an important role in gas hydrate distribution 5 .
Specific geological characteristics of the gas hydrate reservoirs are summarized in Table 1. Compared to gas hydrate samples from the Mackenzie Delta of Canada, gas hydrate in the QMP can be regarded as a new type gas hydrate characterized by a relatively shallow buried depth, thin thickness of permafrost, complex gas components, and coal-bed methane gas source, which is of scientific, economic, and environmental significance 4 .

Numerical simulation methods
When the lithology of the gas hydrate reservoir matrix is sandstone and formation fracture development is small, the gas hydrate will be disseminated to fill the pores of the sandstone 17,18 . In previous studies of the acoustic velocity characteristics of pore-filling gas hydrate reservoirs, marine gas hydrates have been studied in the early stages of formation. The effective medium theory, contact media theory, and relevant theory of porous media were used to establish the different forward modeling velocity models, included the EMT model and the TPT model [19][20][21][22][23] . In order to associate the microscopic distribution form with the physical parameters of gas hydrate formation, the elastic wave velocity model based on effective medium theory can be used 24 . The model is established based on the effective medium theory, considering the rock mechanics of elasticity and statistics; the combination of the macroscopic and microscopic aspects of wave propagation in the medium has a relatively wide range of use. In this paper, the EMT model is used to carry out the research work on the pore-filling gas hydrate reservoir in the QMP.
Because there are many parameters that must be set in the simulation of elastic wave velocity using the TPT model, there are some limitations in practical application 22,25 . In order to verify the accuracy of simulation results by the EMT model, the TPT model is selected for comparison study.
When the cracks of gas hydrate reservoir are more developed, the gas hydrates usually fill the cracks in the formation rocks with meshes, nodules, or veins 3,26 . Because the dips of the cracks in the gas hydrate reservoir are generally relatively steep, and the distributions of cracks are often controlled by the regional tectonic principal stresses, the cracks will in general be directionally aligned, leading to anisotropic characteristics in the gas hydrate occurrence zones 27 . If we were to assume that the gas hydrate filled the pores in isotropic rock, then the velocities of the forward modeling used for gas hydrate saturation inversion would yield large errors 26,28 .
For the stratigraphic anisotropy caused by directional alignment of the cracks, many previous studies have been carried out, and a variety of simplified models are proposed, including a laminated media model 29 , a media model where cracks are embedded in the pores 30,31 , periodic thin interlayers and expansion model 32 , and so on. For the fracture-filling gas hydrate reservoir, some scholars have used the TCLM model to carry out the acoustic velocity simulation and gas hydrate saturation estimation study in several gas hydrate potential areas, and have obtained the better application effect 26,33,34 . Therefore, this paper also selects the TCLM model to carry out related research work on the fracture-filling gas hydrate reservoir in the QMP.
Numerical simulation of pore filling reservoir. Acoustic velocity forward simulation. In the actual use of the model, the gas hydrate reservoir can be regarded as the equivalent homogeneous medium, where microscopic distribution patterns affect the elastic parameters, and we study the characteristics of the elastic wave. According to the effective medium theory, the P-wave and S-wave velocities of gas hydrate formations are as follows: where V p is the P-wave velocity, V s is the S-wave velocity, K sat is the bulk modulus of the fluid-saturated rock formation, μ sat is the shear modulus of fluid-saturated rock formation, and ρ b is the density of the rock formation. In order to determine the physical parameters (V p , V s , ρ b ) of the gas hydrate reservoir in the study area, we need to determine the K sat and μ sat . The calculation of the two modulus parameters is referred to as the EMT model by Helgerud et al. 24 . The model considers the impacts of the formation physical parameters changing by mineral constituent, porosity, and gas hydrate saturation of the formation rock, and can be used to simulate the change characteristics of the gas hydrate reservoir in the QMP.
In order to compare the forward simulation results with the TPT model, the TPT model by Domenico 35 is selected to carry out the relevant forward simulation work. The P-wave and S-wave velocities of the TPT model are as follows: where V p is the P-wave velocity, V s is the S-wave velocity, φ eff is effective porosity, μ is the average rigidity of the matrix, ρ m is the average density, ρ f is the density of the fluid phase, β is the proportionality coefficient, k is the coupling factor, and C b , C f , and C m are the compressibilities of the solid phase, the fluid phase, and the matrix, respectively. The use of equations (3) and (4) for computing the P-wave and S-wave velocities are as reported in Tinivella 22 .
Gas hydrate saturation inversion. When the acoustic velocity obtained from the forward modeling is used for inversion for the gas hydrate saturation, we need to compare the results with the actual borehole acoustic logging in the study area. The difference between the actual acoustic velocity of the gas hydrate reservoir and the acoustic velocity as simulated using the water-saturated pore space reflects the saturation of gas hydrate. The actual P velocity data of the gas hydrate reservoir can only be obtained relative to the normal acoustic travel time logging in the gas hydrate drilling borehole of QMP. So this research paper on gas hydrate saturation inversion of pore-filling reservoir will only be based on the P-wave velocity obtained from the EMT model used in the inversion to obtain the gas hydrate saturation. Figure 2 shows a flowchart for the inversion for the gas hydrate saturation using the P-wave velocity in the pore-filling reservoir. We set an initial gas hydrate saturation value, and then input the parameters of the formation rock matrix in the study area, such as porosity, density, bulk modulus, and shear modulus, and obtain the P-wave velocity under this saturation condition by use of equation (1). We compare the difference between the simulated P-wave velocity and the corresponding P-wave velocity observed in the acoustic log. If the difference between the values is within a preset permissible error range, then the gas hydrate saturation is considered to be the actual gas hydrate saturation of the reservoir. Otherwise, if the difference between the predicted and measured values is not within the permissible error range, then the modified saturation's initial value is repeated until the error is met.
Numerical simulation of fracture filling reservoir. Acoustic velocity forward simulation. Some drilling boreholes have implemented advanced ultrasonic imaging logging in scientific drilling project of QMP. This logging method can effectively identify formation fractures and analysis of fracture occurrence, and thus provides a basis for numerical simulation of the acoustic log in a fracture-filling gas hydrate reservoir. Figure 3(a) shows an ultrasonic imaging log of a gas hydrate sample in the study area. From the image, it can be seen that the cracks of this interval are well developed, mainly consist of high angle fractures, and the gas hydrates are mainly located in these high angle cracks. Figure 3(b) shows the two component laminated media model of the corresponding fractured reservoir, which is used to simulate the acoustic velocity characteristics of the gas hydrate in the cracks.
The TCLM model is composed of two component elements: I and II ( Fig. 3(b)). The component I is the anisotropic fracture medium, and the cracks are 100% filled with gas hydrate. The component II is isotropic pore medium with fully saturated water in the pore. The volume fraction of the fracture medium in the component I is η 1 , fracture porosity is φ 1 , and assumes φ 1 = 100. The volume fraction of the pore medium in the component II is η 1 , and saturated water porosity is φ 2 . The elastic parameters of the fractured gas hydrate reservoir can be expressed as follows: where M is the combination of arbitrary elastic parameters of the components I and II in Fig. 3(b). For the component I, the gas hydrate completely fills the cracks, and the elastic parameters can be replaced by the elastic parameters of pure gas hydrate. For component II, each elastic parameter can be calculated by the EMT model. The P-and S-wave phase velocities can be calculated from the following equations using the Lame constants λ and μ 36 : where ϕ is the propagation angle relative to vertical, V P is the P-wave phase velocity, V S h is the horizontally polarized S-wave phase velocity, and V S v is the vertically polarized S-wave phase velocity. The P-and S-wave velocities calculated from equations (14), (15), and (16) are the phase velocities for the fractured gas hydrate reservoir. When the P-and S-wave phase velocities are subsequently used to invert for the gas hydrate saturation, they need to be converted into the transverse velocities. For a vertical borehole, modeling using an incidence angle of 0° represents wave propagation for a horizontal fracture and using an incidence angle of 90° represents wave propagation for a vertical fracture. The incidence angle is related to the angle of fracture, and the angle of incident angle reflects the dip angle of the fracture in the formation 26 . When the dip angle of the fracture is horizontal or vertical, the transverse velocity and phase velocity are consistent. When the fracture dip angle is any other angle, the phase velocity can be converted into the transverse velocity by use of Thomsen's equation 37 .
Gas hydrate saturation inversion. Figure 4 shows a flowchart for the inversion of gas hydrate saturation using the P-wave velocity in a fracture-filling gas hydrate reservoir. We set an initial gas hydrate saturation value, and then input the elastic parameters of the components I and II for the study area, such as porosity, density, bulk modulus, and shear modulus. We obtain a P-wave phase velocity under this saturation condition by use of equation (14). According to the statistics of fracture occurrence in the gas hydrate reservoir based on the ultrasonic imaging well logging image, it is necessary to determine whether the simulated P-wave phase velocity can be converted into the P-wave transverse velocity. Finally, we compare the difference between the simulated P-wave transverse velocity and the corresponding actual P-wave velocity determined from the acoustic log. If the difference is within the permissible error range, then the gas hydrate saturation is considered to be the actual gas hydrate saturation of the reservoir. Otherwise, if the difference between above two is not within the permissible error range, then the modified saturation's initial value is repeated until the error is met.

Results and Discussion
Research on acoustic velocity and gas hydrate saturation from pore filling reservoir. The hole DKXX-13 has implemented an integrated conventional well logging of the whole borehole section. The logging projects include natural gamma, well diameter (caliper), density, acoustic travel time, resistivity logging, and so on. Figure 5 shows the conventional well logging curves of the gas hydrate in a sample drilling interval (X29.8~X32.2 m). Due to the presence of gas hydrate, this interval has obvious responses in the resistivity and acoustic travel time logging curves, and shows the response characteristics of high resistivity (RT) and low acoustic travel time (AC), which are consistent with the logging response characteristics of gas hydrates in other permafrost regions around the world 14 . The natural gamma and density curves are mainly the reflection of strata lithology, which is mainly composed of argillaceous sandstone. At the same time, the well diameter curve can also be used as the effective parameter to identify the gas hydrate reservoir; the interval of the well diameter curve is relatively smooth and unchanging, and the other well log responses are not related to the borehole condition.
In hole DKXX-13, we also carried out ultrasonic imaging logging during the drilling process. High quality borehole formation information was obtained, which can be used to analyze and study the development and occurrence distribution of fractures in the gas hydrate reservoir. In this paper, the cracks in the borehole are collated by means of the ultrasonic imaging logging image, and the characteristics of fracture occurrence with depth change are analyzed (Fig. 6). As shown for the depth range of X29.8~X32.2 m, 267 cracks are present, and the absolute frequency is 9.7 stripes/10 m, but we also see that the stratigraphic fracture in the whole borehole is not developed. In depth range from X29.8~X32.2 m of the reservoir, 32 cracks are collected, the absolute frequency is 13.7 stripe/10 m; the dip angles are mainly distributed over a range of 45°~75°, which indicated high angle fracturing 38 . The above analysis shows that the development of high angle fracture in gas hydrate reservoir of the hole   Acoustic velocity characteristics analysis. By use of the acoustic velocities simulated by the EMT and TPT models, the formation porosity, φ, needs to be determined. Formation porosity is a key parameter for gas hydrate estimation, so appropriate log data should be selected to estimate the formation porosity of hole DKXX-13 in the study area. The log data that can be applied to determine the formation porosity include density, acoustic travel time, and resistivity logging. When the acoustic log is used to determine formation porosity, the calculated porosity needs to be corrected using the regional core data because gas hydrate is buried shallow in the study area. However, the core data is always insufficient, and it is difficult to determine the compaction correction coefficient, so the acoustic log is not suitable. The Archie formula should be used to calculate porosity for the resistivity log, and the Archie constants and formation water resistivity need to be known for the Archie formula 39 . However, the above two parameters are generally determined by some empirical equations, which will lead to significant estimation error. The density log is less affected than the resistivity and acoustic logs in gas hydrate reservoir, and can generally reflect the situation of formation porosity, so we select the density log to estimate the formation porosity in this paper.
The intensity of scattering gamma ray is measured for the density log, which reflects the electron density of the formation and the volume density of rock (ρ b ). The porosity estimation by density log can be expressed as 40 : where ρ ma is the matrix density and ρ ma = 2.64 g/cm 3 by use of the mineral component of the rock skeleton, ρ f is the fluid density, and ρ f ≈ 1.00 g/cm 3 by use the density of pure water. Because of the high content of mud in the gas hydrate reservoir of hole DK XX-13, equation (17) needs to include a mud correction. So equation (17) can be corrected as follows 40 : where V sh is the volume content of the mud, SH is the content index of the mud, ρ sh is the density of the mud, and GR, GR min , and GR max are the values of the natural gamma log in the target layer, in the pure sandstone layer, and in pure mud layer, respectively, and GCUR is the Hilchie index 41 . Equation (18) is used to calculate the formation porosity in depths of X28.0~X34.0 m of hole DKXX-13 (Fig. 7). The mud content in the depth ranges of X29.3~X30.0 m and X30.7~X31.1 m are quite high, so the porosities of those intervals are 0, which indicates that the lithology is mudstone. Because of the effect of borehole enlargement in the depth ranges of X28.0~X29.3 m and X33.4~X34.0 m, the calculated porosities of those intervals are too large. In the depth range of X29.8~X32.2 m, the calculated porosities vary from 1.0~8.0%, the average value is 4.0%, and standard error is 0.1%. According to the porosity test results of different granulometric class sandstone, the test porosities vary from 2.2~5.1%, the average value is 3.5% 42 . The calculated porosities are basically the same as the core porosity test results. The result indicates that the formation porosity of the gas hydrate reservoir is quite low, which means that the lithology of the reservoir in hole DKXX-13 is more compact.
According to the natural gamma log data, the shale content of borehole formation in hole DKXX-13 is high, and large parts of the borehole intervals are mudstone. Therefore, only the non-mudstone formation is carried out acoustic velocity characteristics simulation by use of the EMT model and TPT model. If we assume that the borehole formation is water saturated (i.e., the gas hydrate saturation S h = 0), then the P-wave velocities in the  (Fig. 8). Meanwhile, the simulation results of the TPT model are comparable (Fig. 9). Some values of the parameters in above two velocity models are shown in Tables 2 and 3.
From Fig. 8, the P-wave velocity curve for the water-saturated formation obtained from the EMT model and the actual curve of the P-wave velocity log are almost coincident in the depth ranges of X28.0~X29.8 m and X32.2~X34.0 m (the non-gas hydrate interval), and the two curves coincide in the depth range of X33.5~X34.0 m. We can be seen that the EMT model and the parameters used can be used to reliably analyze the velocity characteristics of the pore-filling gas hydrate reservoir in the study area.
From Fig. 9, the results of the P-wave velocity for the water-saturated formation using the TPT model are consistently less than the values obtained by the EMT model. In the non-gas hydrate interval, the P-wave velocity curve for the water-saturated formation from the TPT model and the actual P-wave velocity log are almost coincident, but do not coincide, which can not reflect the actual saturation water formation situation. The P-wave velocity curve by the TPT model has some deviation from the actual P-wave velocity curve. Therefore, the EMT model can be better used to analyze the P-wave velocity and gas hydrate saturation characteristics of the pore filling gas hydrate reservoir in the study area.    The difference between the actual log P-wave velocity and saturated water P-wave velocity reflects gas hydrate or free gas saturation concentration, and can be used to qualitatively identify gas hydrate reservoir 22,43 . The method for the identification of gas hydrate or free gas in the formation by using the simulated saturated water P velocity curve is as follows: If the actual log P-wave velocity is higher than that for a water-saturated formation, then the formation interval may contain gas hydrate; If the actual log P-wave velocity is lower than that for a water-saturated formation, then the formation interval may contain free gas.
Gas hydrate saturation estimation. When we estimate the gas hydrate saturation by inversion using the EMT model, we can set a relatively large initial guess for the gas hydrate saturation, and use forward modeling of the P-wave velocity curve to determine the scope of the change of the gas hydrate saturation Then the scope of the search for a fit can be reduced to reduce the inversion workload. From Fig. 10, with the increase of gas hydrate saturation (S h ), the simulated P-wave velocity increases, and when = .
S 85 0% h , the whole curve of P-wave velocity by simulated coincides with the actual P-wave velocity curve from the acoustic login the depth range X28.0~X34.0 m of hole DKXX-13, which reflects the range of gas hydrate saturation of 0~85.0%.
The gas hydrate saturation for depths of X28.0~X34.0 m in hole DKXX-13 is inverted using the EMT model (Fig. 11). At the same time, in order to compare and analyze the inversion results, gas hydrate saturation is also estimated from the resistivity logging using the Archie equation 39 . From the blue curve, the estimated gas hydrate saturation for depths of X28.0~X34.0 m is relatively large. Because the lithology of the borehole formation is relatively dense, the porosity of the formation is small, and so the volume percentage of gas hydrate in the pore of the borehole formation is large enough to be detected. In the intervals of abnormal gas hydrate concentrations (X28.0~X29.3 m, X32.3~X32.7 m, X32.8~X33.4 m), the range of gas hydrate saturation is 3.0~86.0%, the average value is 52.3%, and standard error is 1.7%. The inversion results show existing gas hydrate occurrences for depths of X30.0~X30.2 m, X30.3~X30.4 m, X31.1~X31.6 m, X31.7~X31.9 m, and X32.0~X32.2 m. The range of gas hydrate saturation is 13.0~85.0%, the average value is 61.9%, and standard error is 1.9%, which reflects the relatively high gas hydrate saturation.
From the pink curve in Fig. 11, the gas hydrate saturation curve from the EMT model and the saturation curve by the density equation are comparable, although there are some differences in the peak values of the two  .5~X30.7 m, X31.1~X31.6 m, X31.7~X31.9 m, and X32.0~X32.2 m; the range of gas hydrate saturation is 18.0~77.0%, the average value is 61.6%, and standard error is 1.3%. The estimation results from the Archie equation show that while the peak value of gas hydrate saturation is lower than the EMT model, the average value of gas hydrate saturation is basically the same. What is more, the actual core test data for gas hydrate bearing sediments of pore-filling reservoir in study area show the gas hydrate saturations range from 52.9~61.8%, the average value is 58.8% 44 . The above study results are basically in accordance with the calculated result made by the EMT model. Thus the estimation of the gas hydrate saturation in pore-filling gas hydrate reservoir is reliable using the EMT model.
Research on acoustic velocity and gas hydrate saturation from fracture filling reservoir. The gas hydrate samples were acquired for depths of X10.9~X14.2 m in hole DKXX-19; the total thickness of the occurrence of gas hydrate is 33.7 m (Fig. 12). Where the gas hydrate occurs, conventional logging curves also show the response characteristics of high resistivity and low acoustic travel time. For depths of X11.4~X11.9 m, X12.6~X12.7 m, and X12.8~X12.9 m, conventional logging curves show the response characteristics of low resistivity and high acoustic travel time, which reflect a lack of gas hydrate in these borehole formations. For the caliper curve, multiple intervals display diameter fluctuation variation in the occurrence of gas hydrate, which reflect the high degree of fracture development in hole DKXX-19.
Ultrasonic imaging logging work was also carried out in hole DKXX-19 during the drilling process, and the characteristics of fracture occurrence with depth change are also analyzed (Fig. 13). As shown in all sections, 1986 cracks are present, and the absolute frequency is 33.1 stripes/10 m, which shows that the stratigraphic fracture in the whole borehole is quite developed. For depths of X10.9~X14.2 m in the gas hydrate reservoir, 157 cracks are collected, the absolute frequency is 46.6 stripes/10 m; the maximum absolute frequency is 64 stripes/10 m, and the dip angles are mainly distributed in a range of 50°~70°, which indicates high angle fractures. The foregoing analysis shows that the development of high angle fractures in the gas hydrate reservoir of the hole DKXX-19 is relatively predominant, which reflects gas hydrate filling in the cracks of the borehole formation. The results agree with the response characteristics of conventional logging curves and the gas hydrate sampling from field drilling.
Acoustic velocity characteristics analysis. Equation (18) is used to calculate the formation porosity for depths of X10.0~X15.0 m in hole DKXX-19 (Fig. 14). The calculated porosities for depths of X11.5~X11.9 m and X12.8~X12.9 m are 0, which indicate that the lithologies in those intervals are most likely mudstone. For depths of X10.9~X14.2 m, the calculated porosities vary in the range of 5.0~20.0%, the average value is11.9%, and standard error is 0.2%. According to the core test results for the reservoir properties in Qilian mountain permafrost, the porosities range from 4.2~13.3%, the average value is10.6%, the calculated porosities is basically in accordance with the core test results 45 . The result indicates that the formation porosity of the gas hydrate reservoir is quite high. Because the density of gas hydrate is similar to that of the formation pore water, the calculated porosity from the density log can be approximated as the total porosity of the formation, and will include the fracture porosity and the porosity of pores saturated with water 46,47 . For the fractured gas hydrate reservoir, the P-wave and S-wave velocities simulated by the TCLM model are affected by the changes in the gas hydrate content and the fracture dip angles in the cracks. Therefore, it is necessary to analyze the relations between these two parameters and the P-wave and S-wave velocities, so as to better carry out the relevant analysis.
In order to determine the effect of the gas hydrate content on the P-wave and S-wave velocities, we set the formation porosity φ = 30.0%, the fracture dip angle is set to vertical and horizontal, respectively, and the volume and parameter values of each component of the formation rock skeleton are shown in Table 3. The P-wave velocity (V P ) and horizontally polarized S-wave velocity (V S H ) versus volume fraction of gas hydrate are simulated using the TCLM model ( Fig. 15(a)). When = V 0 h , which indicates that there is no hydrate in the fracture, the simulated velocities represent the saturated water situation of the formation. When = V 1 h , which indicates that the model pore space is completely filled with gas hydrate in the fracture, then the simulated velocities represent the parameter values of pure gas hydrate. For the vertical and horizontal fractures, when the volume fraction of gas hydrate increases, V P and V S H increase rapidly. For the same volume fraction of gas hydrate, when the fracture dip angle changes from horizontal to vertical, the V S H increases rapidly, and the two S-wave velocity curves are clearly separated; the V S H increases slowly, and the two P-wave velocity curves are effectively coincident. To analyze the effect of fracture dip angle on the P-wave and S-wave velocities, we set the formation porosity φ = .
30 0%, and the volume fraction of gas hydrate = .   simulated using the TCLM model ( Fig. 15(b)). When the fracture dip angle is small (<40°), V P and V S H change slowly. When the fracture dip angle is greater than 40°, the V P and V S H increase slowly with increasing fracture dip angle, and the range of velocity increase is quite little.
Through the above analysis, we see that when φ < . 30 0%, the change of volume fraction of gas hydrate filling in the fracture has great influence on the P-wave and S-wave velocities, whereas the change in the fracture dip angle has little influence on the velocities. Given that the average value of the formation porosity for depths of X10.0~X15.0 m in hole DKXX-19 is 11.9%, and the fracture dip angles vary mainly from 50°~70°, the fracture dip angle of the borehole formation can be assumed to be vertical. Therefore, the simulated P-wave phase velocity of hole DKXX-19 can be converted directly to the P-wave transverse velocity. Assuming = S 0 h , the P-wave velocity is simulated by the TCLM model for depths of X10.0 ~X15. 0 m in hole DKXX-19 (Fig. 16). The P-wave velocity curve for a water-saturated formation using the TCLM model and the actual curve of the P velocity log are almost coincident for depths of X10.0~X10.9 m and X14.2~X15.0 m (i.e., the intervals with no gas hydrate), and the two curves coincide for depths of X10.6~X10.8 m and X14.6~X14.7 m. It can be seen that the TCLM model and parameter settings can be used to reliably analyze the velocity characteristics of the fracture filling gas hydrate reservoir in the study area.
Gas hydrate saturation estimation. When we obtain gas hydrate saturations by inversion using the TCLM model, the set inversion parameter is V h , and the conversion relation between V h and S h is S h = V h /φ, where φ is the total porosity by equation (18). Therefore, if we set a relatively large initial V h , and perform forward modeling of the P-wave velocity curve to determine the scope of the change of the gas hydrate saturation, then the search scope can be reduced to reduce the inversion workload. From Fig. 17, when V h increases, the simulated P-wave velocity increases, and when = V 20% h , the whole curve of P-wave velocity by simulated is located at the top of the actual P-wave velocity curve as determined from the acoustic log, which reflects the range of gas hydrate saturation is 0~20% in depths of X10.0~X15.0 m of hole DKXX-19.
The gas hydrate saturation for depths of X10.0~X15.0 m of hole DKXX-19 is detemined using the TCLM model (Fig. 18). In the intervals of abnormal gas hydrate (X10.1~X10.6 m, X14.3~X14.6 m, and X14.7~X15.0 m),  The actual core test data for gas hydrate bearing sediments of fracture-filling reservoir in study area show the gas hydrate saturations range from 72.3~75.8%, the average value is 73.5% 44 . What is more, the rang of simulated gas hydrate saturation is 0~80.0% based on porosity of 7.0% and 10.0% by Liu et al. 48 . The above previous study results are basically in accordance with the calculated result made by the TCLM model. This indicates that using the TCLM model to estimate gas hydrate saturation of fracture-filling reservoir is available.
Through the simulation results in this paper, in unfractured gas hydrate bearing sands in the QMP, China, an advantage of the EMT model by Helgerud et al. over the TPT model by Domenico is that whereas the EMT model accurately predicts P-wave velocities. In fracture-filling gas hydrate reservoirs, the method of combining ultrasonic log data with the TCLM model can be effectively used to evaluate gas hydrate saturations of gas hydrate formation. When gas hydrates usually fill the cracks in formation, the dip angle of the cracks has a certain influence on the estimation of gas hydrate saturations by the TCLM model. In this paper, fracture dip angles obtained by gas hydrate drilling borehole are steep and assuming the cracks are vertical, so a gas hydrate saturation error associated with the fracture dip angle error may exist. The fracture dip angle is related to the regional tectonic environment and presents the orientation characteristic, indicating that the fracture-filling gas hydrate is different from the pore-filling gas hydrate, and gas hydrate bearing sediments have anisotropic properties. For the gas hydrate saturation estimated by the EMT model, we have referred some empirical parameters from permafrost zone of other countries, so the inversion results may have some limitations. Although there are some empirical parameters by using above methods, the methods proposed in this paper are very necessary for the study of gas hydrate in this area. The research results could provide important reference for gas hydrate reservoir logging evaluation and seismic exploration in the QMP, China.

Conclusions
(1) The combination of ultrasonic imaging logging and drilling core data can be used to identify the filling type of gas hydrate reservoir in Qilian mountain permafrost. The hole DKXX-13 illustrates the typical pore filling gas hydrate reservoir, and the hole DKXX-19 represents the typical fracture filling gas hydrate reservoir. (2) Using the EMT and TPT models of pore filling gas hydrate reservoirs, the characteristics of the acoustic velocity for hole DKXX-13 are analyzed. The P-wave velocity simulated by the EMT model is more consistent with the actual log data than the TPT model.   range of gas hydrate saturation is 13.0~85.0%, and the average value of the gas hydrate saturation is 61.9%, which is largely in accordance with the results derived using the standard Archie equation. (4) Using the TCLM model of a fracture-filling gas hydrate reservoir, the characteristics of the acoustic velocity of hole DKXX-19 are analyzed. The change in the volume fraction of gas hydrate filling the fractures greatly influences the P-wave and S-wave velocities, but the change in fracture dip angle has much less effect on the velocities. The P-wave phase velocity simulated using the TCLM model can be transformed directly into the P-wave transverse velocity. (5) The occurrence of gas hydrate is estimated to occur at depths of X11.1~X11.3 m, X11.9~X12.2 m, X12.5~X12.6 m, X12.7~X12.8 m and X12.9~X14.2 m of hole DKXX-19.The inversion estimation results using the TCLM model estimate that the range of gas hydrate saturation is 14.1~89.9%, and the average value of the gas hydrate saturation is 69.4%, which is largely in accordance with the core test results.