An Improved Method for Calculating Paleoheat Flow from Vitrinite Reflectance Profiles

INTRODUCTION

It is well established that the vitrinite reflectance (R_o) depends on the temperature history that a basin has undergone. In general. R_o is a function of maximum paleotemperature and effective heating time. Thus, the paleotemperature of sedimenu and the value of vitrinite reflectance are controlled by the flow of heat Q (t) into the basin and by the thermal parameters of the sediments. Therefore, in principle, it is possible to calculate the thermal vitrinite reflectance of the basin if the thermal history is known and vice versa to invert the thermal history of the basin once the vitrinite reflectance has been determined through experiments. Over the last two decades, various techniques have been developed for calculating paleoheat flow using the measured vitrinite reflectance (Pang et al., 1993; He and Lerche, 1989; Lerche, 1988; Lerche et al., 1984; Waples, 1980). The method by Lerche et al. (1984) permits the heat flux to increase or decrease linearly in time. However, the methods by Lerche (1988) and He and Lerche (1989) allow the simultaneous determination of not only geological but also geochemical and geothermal parameters.

The objective of this paper is to determine the history of heat flow at the base of basins using measured vitrinite refleclance. Pangs model (Pang et al., 1993) for the calculation of the maturity index is adopted. The maturity index can be determined either through measured vitrinite reflectance (R_o) or through the temperature history T (t). As the temperature history T (t) depends, to a large degree, on the history of heat flow Q (t) into the basin, the maturity index calculated will match the index determined from measurement only if the heat flow Q (t) is accurately identified. The problem of identifying the thermal history of the basin can be cast into an optimization problem: Q(t) is located and such systems thermal parameters are obtained that the calculated maturity index determined from the temperature field based on Q (t) is as close as possible to the maturity index determined from measurement. The rest of the paper is organized as follows. Firstly. the optimization model for inverting the heat flow history is constructed. Then a genetic algorithm for solving the optimization problem is presented. Finally, the robustness of the method thus dcveloped is illustrated through an application example.

MATHEMATICAL MODEI

Lerche et al. (1984) found that the vitrinite reflectance, R_jk, in a well can be converted to a normalized maturity index, I_jk, which is a direct measure of the maturity of the sample. namely

where I_jk and R_jk are the maturity index and the measured vitrinite reflectance at the depth of the kth sample of the jth strato; R_jm is the measured vitrinite reflectance of the sample at the maximum depth of all the samples of the jth strata; R、the vitrinite reflectance of the strata at the time when they were deposited (generally around 0.20), although this reflE: ctance seems to he higher in some basins (Lerche et al., 1984). Based on Lerches model, Pang et al. (1993) developed a theoretical expression directly comparable to the index, I_jk, that is

where is the maturity index calculated from the temperature history; T_c and T_v are paleo-surface temperatures (K); T_d isa constant whose value ranges from 160 to 170 K; the subseript. k denotes the kth sample; j denotes the jth strata; i denotes the ith geological age; T_ijk, T_jjk, T_ijm and T_jjm are paleoternperatures (K) along the burial path, which are a function of both time and depth. Two different formulus for calculating temperature have been proposed, respectively by Lerche et al. (1984) and Pang et al. (1993), namely

where T_surface is temperature at ground surface (K); GT is temperature gradient (K/100 m); Z is depth (m). For certain special cases, the temperature can be calculated by using these formulas with reasonable accuracy. In sedimentary hasins, however, the distribution of temperature is usually affected by many parameters. such as thermal conductivity, specific heat, basal heat flow, and radial generation. It is thus necessary to use a more sophisticated model for calculating temperature. In this paper, we use the thermal conduction equation to (alculaW the temperature. The basic assumptions used in the model arc shown as follois: (1) The heat transport in sedimentary basins Is conducted mainly by heat conduction. (2) There is no significant heat transfer by convection. (3) The basins are not sediment-starved and are filled close to the sea level. Based on these assumptions, the governing equation of the heat transfer is

where T denotes the temperature [K]; c is the specific heat [J·kg^-1·K^-1]; ρ is the density of rock [kg·m^-3]; k is the thermal conductivity [W·m^-2·K^-1]; Q (t) is the basal heat flow[W·m^-2]; T_s is surface temperature [K], and T₀ is initial temperature [K].

In the porous media system. the parameters c, ρ and k are the function of porosity, and can be expressed as follows.

The specific heat c can be calculated by

where c_r and c_f are the specific heat of rock and fluid, respeclively [J·kg^-1·K^-1].

The density of formation is expressed as

where ρ_r and ρ_f are the density of rock and fluid, respectively [kg·m^-3]

The thermal conductivity k is determined by

where k_r and k_f are the thermal conductivity of rock and fluid. respectively [W·m^-2·K^-1].

Based on the discussion mentioned above, there are two ways for calculating the maturity indexes, I_jk and .Theoretically, both I_jk and should be the same. Hence, if the parameters of equation (2) are correctly determined, the theoretical maturity index, , will be close to the index, I_jk, calculated by the measured vitrinite reflectance. From this consideranon, the identification of paleoheat flow can be included in the following nonlinear optimum formula

equation (9) is subject to determined by equation (2) with T governed by equations (5)-(8) and

To calculate the value of the objective function, the partial differential equation (5) subject to conditions (6)-(8) should be solved.

It is almost impossible to use enumerative methods or gradint-based methods to solve this kind of problem. Hence, in this paper, a genetic algorithm is used to solve this complex and highly non-linear optimization problem.

PROCEDURE OF GENETIC ALGORITHMS

Genetic algorithm is a solution-search technique developed by Holland (1975). The technique uses the mechanisms of nat ural selection to search for optimal solutions through decision space (Goldberg, 1989).Genetic algorithm offers a distinct advantage over conventional optimization techniques because it is not affected by the complexity, convexity, or linearity of the objective function. For complicated problems, especially discontinuous or highly nonlinear and non-convex problems, genetic algorithm is usually computationally superior to traditional methods, such as gradient-based methods.

Genetic algorithm uses a random search procedure inspired by biological evolution and cross-breeding trial designs, and allows only the "fittest" design to survive and propagate to successive generations. In a genetic algorithm, first of all, the decision variable values are encoded as substrings of binary digits. These substrings of decision variables are concatenated to form longer strings representing a design or solution. In this paper, the traditional binary-coded substrings are used. For example, if the decision variable refers to x_i∈[a_i, b_i], this variable can be encoded as a binary substring with length m, where the length depends on the required precision of the representation a_i, the relationship between the number of design length and precision is

The entire of "population" of such designs represents a "generation".The procedure includes the following major steps.

(1) Generate randomly an initial population of binary strings representing possible designs.

(2) Decode the binary strings to obtain decision variables for each design in the population. Use the simulation model to calculate the fitness of each member of the current generation based on the objective function value.

(3) Start with the current population of designs, produce subsequent generations in the following three steps using the genetic operations of reproduction and breeding applied to the current generation. (a) Reproduction. Enhanced reproduction employing "elitism" and "tournament selection" here has been used to improve the genetic algorithm efficiency and performance. In an elitist strategies, the least fittest individual of the current generation is replaced by the fittest individual of previous generations (Davis, 1987). Tournament selection guarantees that the best design has probability 1. 0 of prevailing in the competition while the worst and middle-ranking designs have probabilities of 0 and 0.5, respectively. (6) Breeding. Produce new individuals by using genetic operators on the individuals chosen in the selection step. There are two main kinds of oper atOIS; ①Recombination. A new individual is produced by recombining features of a pair of'parent' solutions. The better performing individuals produce more offsprings.②Mutation. A new individual is produced by slightly altering an existing one occasionally changing some of the bits of the newly-formed point. A bit value 0 will become 1 and vice versa, with each bit only at a very small probability p_m, .

(4) Population update. The set is altered, typically by choosing to remove some or all of the individuals in the existing generation (usually beginning with the least fitness) and replacing these with individuals produced in the breeding step. The new population thus produced becomes the current generation.

(5) Repeat steps (2)-(4). The best point found is always recorded. Termination of searching can be conducted by specifying a total number of function evaluations.

The success and performance of genetic algorithm depends on several parameters, such as population size. number of generations, and the probabilities of crossover and mutation. 1t has been suggested (Goldberg, 1989) that good genetic algorithm performance requires the choice of high-crossover and low-mutation probabilities and a moderate population size.

Based on the genetic algorithm, a numerical scheme has been developed for calculating paleoheat flow from the measured vitrinite reflectance. which can be summarized as follows.

Step 1:determine the sublinear segments based on the measured vitrinite reflectance.

Step 2:calculate the index I_jk by equation (1).

Step 3:assign a feasible value to Q (i) by genetic algorithm.

Step 4:calculate the index by equation (2), with the temperature T determined by solving equation (5) subject to conditions (6)-(8).

Step 5:calculate the objective function by equation (9).

Step 6:identify the fitness and if the fitness satisfies the criterion, go to the next sublinear segment, otherwise, go back to step 2.

APPLICATION TO PEARL RIVER MOUTH BASIN, SOUTH CHINA SEA

To illustrate the application of the method presented in this paper, the R_o data from well Yj35-1-1 of the Pearl River Mouth basin are used to calculate the paleoheat flow.

Geological Setting and Parameters

The Pearl River Mouth is located at the northern continental shelf of the South China Sea (Fig. 1). The regional tectonic background of the continental shelf falls in a passive continental margin. The development of the Tertiary basin in the South China Sea shows a series of characteristic rifting events (Chen et al., 1993; Yang et al., 1992) and several faulting and subsidence phases (john and Pigo, 1987).Two angular unconformities and one parallel unconformity occur in the sequence of the Pearl River Mouth basin, the Early Paleocene, the Oligocene to Miocene and the Late Miocene to Pliocene. Each unconfortuity could be the reflection of the end of a rifting event re fated to lithosphere uplift after mantle intrusion. The depositional sequences in this basin in age dimension are shown in Table 1.

Figure  1.  Map showing location of Pearl River Mouth basin on northern continental shelf of South China Sea.

DownLoad: Full-Size Img PowerPoint

Table  1.  DEPOSITION SEQUENCE, FORMATIONS AND AGES FOR THREE BASINS

 | Show Table

DownLoad: CSV

Well Yj35-1-1 is located in the Pearl River Mouth basin (Fig. 1).The oil-immersion shows that the random vitrinite reflectance data (R_o) of the well were plotted, using a iogarithmic axis for R_o and an arithmetic axis for depth, and that the best-fitting lines were regressed by exponential relationship between R_o and the depth for each segment of a R_o profile, ax shown in Fig. 2. Figure 2 also shows four piece-wise linear lg R_o segments with different slopes. Steep segments of the lg R_o profiles correspond to the stages with high depositional rates. while the gentle segments of the 1g R_o profiles are also dated with periods with low depositional rates.

Figure  2.  Diagrams of vitrinite reflectance profile of well Yj35-1-1.

DownLoad: Full-Size Img PowerPoint

The decision variables include the heat flow Q_i (i= 1.2, 3, 4) at the bottom of the basin. According to the piec。-wine linear characteristics of vitrinite reflectance. the heat flow issumed to be a constant within the period in which the curve of lg R_o vs depth is a straight line as shown in Fig. 2.

The parameters used in the model are shown as follows.

The porosity of sandstone and mudstone are calculated, respectively, 6Y

and

the porosity of matrix hulk, φ, is calculated by

where φ_s, and φ_m are the porosity of andstone and mudstonc, respectively: φ denotes the porosity of matrix bulk; p_s is the percentage of sand content.

The initial parameters interval is

genetic algorithm parameters are number of population chosen to be 10. the probability of mutation is P_m =0.1. and the. probability of cross is taken to be p_c=0.5.

RESULTS AND DISCUSSIONS

The paleoheat flow at the well Yj35-1-1 was calculated by using the method discussed in previous section. Figure 3 shows. that the value of the objective function varies with the increase of the generation of genetic algorithm, indicating that after 200 generations. the value of the objective function changed very little. Figure 4 shows that the decision variables change with the genetic generations, indicating that the values of the decision variables converged after being iterated about 280 generations. Figure 5 suggests the paleoheat flow against time in comparison with rifting events. It is noted that the paleoheat flow changed with rifting phases during geological time, indicating that the pattern for Yj35-1-l displays two cycles of heat flow ranging from about 60 Ma to 10 Ma, and these cycles correspond to rifting I and rifting I based on the study of rifting phases by Yang et al. (1992). Based on the work of Yang et al. (1992), rifting ] occurred at about 60.5 Ma, rifting I began at about 23. 3 Ma, and rifting III took place at about 10.6 Ma. (onpared with riftings I and I, rifting [0 can be considered as a post-rifting of rifting I, so the heat flow during rifting ll is lower than that of the other two riftings. In general, the average of paleoheat flow is 8095 mW/m² in the rifting stage. In the postrifting stage with thermal subsidence and thermal decay, the average paleoheat flow is about 70 mW/m². Although it is possible that the gradient of a lg R_o segment may vary due to the thermal accumulative effect of geothermal field variations, the obtained results following the optimal method 1 mentioned above suggest that a an entire rifting event consists of a gentle and a steep lg R_o segment, reflecting different geothermal backgrounds.

Figure  3.  Relation of error and number of iterative generations of genetic algorithm.

DownLoad: Full-Size Img PowerPoint

Figure  4.  Paleoheat flow varies with number of iterative generations of genetic algorithm.

DownLoad: Full-Size Img PowerPoint

Figure  5.  Paleoheat flow against time in comparison with rifting events (Yang et al. (1992)) for Yj35-1-1.

DownLoad: Full-Size Img PowerPoint

CONCLUSIONS

(1) The method proposed in this paper is an improvement of Pangs model (1993). The major advantage of this model is that a heat transfer equation is used to calculate temperature instead of the formula used by Lerche et al. (1984) and Pang et al. (1993).

(2) To solve the models, a highly nonlinear optimization should be removed with a set of decision variables. Genetic algorithm is used to remove the nonlinear optimi2ation. The numerical experiment shows that the genetic algorithm is an effective method for this removal.

(3) This method is applied to the calculation of the paleoheat flow of the Pearl River Mouth basin. The results indicate that a rifting episode is accompanied by a thermal cycle. There are three periods of the thermal cycle for Yj35-1-1. The paleo heat flow values are about 80-95 mW/m² in the rifting stages and about 65 mW/m² in the post-rifting stages.

ACKNOWLEDGMENTS

Acknowledgments are due to China Offshore Petroleum Exploration and Development Corporation for providing support; to Yang Jiaming for his helpful comments and encouragement; to Wu Jinfu, Cui Hanyun and Wu Peikang for their support; to Mike Middleton and Lindsay Collons for their support. The genetic algorithm code was modified from a version provided by David Carroll.