The TOC analysis was performed using a LECO CS–230HC analyzer. Approximately 250–500 mg of 100–120 mesh rock was required. To remove the inorganic carbon in the form of carbonates, it was necessary to acidify the sample before analysis. The samples were rinsed with distilled water to remove the HCl solution. After this, the samples were dried to eliminate moisture prior to analysis. The prepared sample was placed in the oven of the analyzer and heated at 1,100 °C, and the amount of carbon dioxide (CO₂) produced was measured using an infrared cell. The samples were subjected to quality assurance/quality control (QA/QC) analysis in duplicate to ensure that the analysis accuracy is better than ±0.2%.

A Rock–Eval–VI pyrolysis apparatus was used for both the rock pyrolysis and the open system thermal simulations. The rock pyrolysis allows for the detection of the hydrocarbon generation potential of the sample by heating the sample in an open system under non–isothermal conditions. The hydrocarbons released from the sample are monitored by a flame ionization detector (FID). The detection items include free hydrocarbon (S1), pyrolytic hydrocarbon (S2), organic carbon dioxide (S3), hydrogen index (HI), and maximum pyrolysis peak temperature (Tmax).

In the open system thermal simulation, 100 mg of the sample was used in a single experiment. During the heating process, first, the pyrolysis apparatus was rapidly heated to 200°C, and the free hydrocarbons were removed at constant temperature for 3 min. Then, under different heating rates (10°C/min, 20 °C/min, 30°C/min, 40°C/min, and 50°C/min), the sample was heated from 200 to 700°C, and the volume of the products was recorded in real time.

The PY–GC consisted of an SRA–TEPI pyrolyzer, which controlled the pyrolysis temperature, and an Agilent 6,890 gas chromatograph detection system (Xie et al., 2020). The sample was placed in the sample tube, and the tube was placed into the pyrolysis probe. The sample was heated from 200 to 630°C at heating rates of 10°C/min and 30°C/min. The pyrolysis products were collected at a temperature interval of 30°C, and gas chromatography analysis was performed to determine the gas chromatogram relative content of the heavy hydrocarbons (C₁₄₊), light hydrocarbons (C_6–14), and gaseous hydrocarbons (C_1–5) at each temperature. A DB–502 50 m × 0.2 mm × 0.5 um capillary column was used.

In order to compare the hydrocarbon generation characteristics of the open system and the closed system, a closed system gold tube thermal simulation experiment was performed on sample No. 9. The kerogen sample was sealed in a gold tube under an argon atmosphere. The gold tube was placed in an autoclave, and the autoclave was filled with water using a high–pressure pump. The high–pressure water caused the gold tube to deform flexibly, thereby exerting pressure on the sample. The samples were heated at heating rates of 20°C/h and 2°C/h, and the temperature difference of each autoclave was less than 1°C. The pressure was 5 MPa, the temperature range was 150–600°C, and the temperature fluctuation was less than 1. After the experiment, the three components, i.e., the gases (hydrocarbon gas and non–hydrocarbon gas), light hydrocarbons (C_6–14), and heavy hydrocarbons (C_14+`), were analyzed. The determination methods used for the different components have been described in previous studies (Ungerer et al., 1988; Hill et al., 2003).

Assuming that at a certain heating rate l, when a certain temperature j is reached, the hydrocarbon production rate measured in the experiment is X1_lj. Under the same conditions, assuming E_i, A_i, and X_i0, the hydrocarbon production rate calculated using the model is X_lj. If there is a certain group of E_i, A_i, and X_i0 values such that X1_lj − X_lj = 0 for all l and j, then this group of E_i, A_i, and X_i0 values represents the correct parameters. However, due to experimental errors and other reasons, this is actually impossible. Therefore, the values of E_i, A_i, and X_i0 that make X1_lj − X_lj as small as possible are the substitutes. Thus, the objective function can be constructed as Eq 8 to characterize the calculation error. Q ( A , X ) = ∑ l = 1 L 0 ∑ j = 1 J 0 ( X 1 l j − X l j ) 2 (8) where L0 is the number of experiments with different heating rates, and J0 is the number of sampling points from an experimental curve.

In addition, for each parameter in the formula, certain constraints listed in Eq 9 must be met. E_i can be solved by determining the distribution range of the activation energy of the parallel reactions and the activation energy interval of the adjacent parallel reactions. { A i > 0 0 ≤ X i 0 ≤ 1 ∑ i = 1 N X i 0 = 1   o r   | 1 − ∑ i = 1 N X i 0 | ≤ ε ( ε   i s   a   inf i n i t e l y   s m a l l   p o s i t i v e   n u m b e r ) (9)

Thus far, the problem of obtaining the kinetic parameters has been changed to the problem of finding the minimum point at which the non–negative objective function satisfies the constraint conditions.

The abovementioned problem of finding the minimum value that satisfies the constraint conditions is more complicated because in addition to the gradual decrease in the value of the objective function, attention must be paid to the feasibility of the solution, that is, to determine whether the solution is within the range defined by the constraint conditions. Thus, the penalty function method is adopted to turn the constrained extreme value problem into an unconstrained extreme value problem. The process is as follows.

For any constraint condition, a function can be constructed. When the obtained extreme point meets the condition, the function’s value is 0; otherwise, it is a positive number.

If A_i > 0, this constraint condition can be obtained as follows: G 1 ( A i ) = { 0       w h e n   A ≥ 0 Α i 2     w h e n   A ≤ 0 (10) That is,G ₁ (A _i ) = [min(0, A _i ]². For the other constraints, the method is the same as for A_i. Because of space limitations, the constraints on the remaining parameters are not repeated here. After the penalty function of each constraint condition is constructed, the penalty items can be obtained as follows: G 1 ( A i , X i 0 ) = G 1 + G 2 + G 3                     = [ min ( 0 , A i ) ] 2 + [ min ( 0 , X i 0 ) ] 2 + [ min ( 0,1 − X i 0 ) ] 2 + [ min ( 0 , ε − | 1 − ∑ i = 1 N X i 0 | ) ] (11)

Use a sufficiently large positive integer R to construct a penalty function F ( A , i X i 0 ) = Q ( A , i X i 0 ) + R ⋅ G ( A , i X i 0 ) (12)

If the obtained minimum point exceeds the constraint condition, the coefficient R is gradually increased. When R is sufficiently large, the minimum solution of the penalty term is the minimum solution of the objective function, thus turning the constrained extreme value problem into unconstrained extreme value problems, which are relatively easy to solve.

The necessary condition for the existence of the minimum is that the first–order partial derivative of the penalty function is 0.

First find the partial derivative of the objective function: ∂ Q ∂ A m = ∑ l = 1 L 0 ∑ j = 1 J 0 ( − 2 ( X 1 l j − X l j ) ⋅ ∂ X l j ∂ A m ) (13) where ∂ X l j ∂ A m = ∂ ∑ i = 1 N ( X i 0 ( 1 − exp ( − ∫ T 0 T l j A i D l exp ( − E i R T ) d T ) ) ) ∂ A m (14)

When i≠m, partial derivative is 0, so, ∂ X l j ∂ A m = X m 0 ⋅ exp ( − ∫ T 0 T l j A m D l exp ( − E m R T ) d T ) ⋅ ∫ T 0 T 1 D l exp ( − E m R T ) d T (15) ∂ Q ∂ X m 0 = ∑ l = 1 L 0 ∑ j = 1 J 0 ( − 2 ( X 1 l j − X l j ) ⋅ ∂ X l j ∂ X m 0 ) = ∑ l = 1 L 0 ∑ j = 1 J 0 ( − 2 ( X 1 l j − X l j ) ⋅ ( 1 − exp ( − ∫ T 0 T l j A m D l exp ( − E m R T ) d T ) ) ) (16) m = 1, 2, 3, … , N.

partial derivatives of the penalty term are as follows: ∂ G ∂ A m = 2 ⋅ min ( 0 , A m ) (17) ∂ G ∂ X m 0 = 2 ⋅ min ( 0 , X m 0 ) − 2 ⋅ min ( 0,1 − X m 0 ) − 2 ⋅ min ( 0 , ε − | 1 − ∑ i = 1 N X i 0 | ) ⋅ F N ( ∑ i = 1 N X i 0 − 1 ) (18)

FN means: F N ( ∑ i = 1 N X i 0 − 1 ) = { 1     w h e n   ∑ i = 1 N X i 0 − 1 > 0 − 1     w h e n   ∑ i = 1 N X i 0 − 1 < 0 (19) m = 1, 2, 3, … , N

Obtain the partial derivative of the objective function and the penalty term. After the partial derivative is obtained, theoretically speaking, the minimum point should meet the following conditions: { ∂ F ( A i , X i 0 ) ∂ A O m = ∂ Q ∂ A m + R 1 ⋅ ∂ G ∂ A m = 0 ∂ F ( A i , X i 0 ) ∂ X m 0 = ∂ Q ∂ X m 0 + R 1 ⋅ ∂ G ∂ X m 0 = 0 (20) m = 1, 2, 3, … , N

Therefore, if the equations can be obtained accurately, several possible minima can be obtained, which can be used to solve for the 2 × N undetermined kinetic parameters (A_i, X_i0) and to complete the model calibration. Although it is impossible to find an exact solution for such a complex non–polynomial function as the above equations, an approximate solution can be found.

A variety of optimization algorithms are available to solve the unconstrained value problem. In this study, the variable–scaling method with fast convergence speed and no need to calculate the cumbersome second derivative matrix and its inverse matrix was selected for the optimization calculation. The detailed derivation of the variable–scale optimization algorithm can be found in the literature (Powell, 1978).

When calibrating the overall reaction model, the difference from the parallel first–order reaction is the change in the order of the reaction (n) and the change in the constraint conditions.

Compared with the parallel first–order reaction model, the constraints of the overall reaction model are as follows: { E > 0 A > 0 n > 0 (21)

In addition, the corresponding penalties are as follows: G ( A , E , n ) = G 1 + G 2 + G 3 = [ min ( 0 , A ) ] 2 + [ min ( 0 , E ) ] 2 + [ min ( 0 , n ) ] 2 (22) F ( A , E , n ) = Q ( A , E , n ) + R 1 ⋅ G ( A , E , n ) (23)

The partial derivative of the objective function to the activation energy is as follows: ∂ Q ∂ E = ∑ l = 1 L 0 ∑ j = 1 J 0 ( − 2 ( X 1 l j − X l j ) ⋅ ∂ X l j ∂ E )             = ∑ l = 1 L 0 ∑ j = 1 J 0 2 ( X 1 l j − X l j ) ⋅ ( 1 − ( 1 − n ) ∫ T 0 T l j A D exp ( − E R T ) d T ) n 1 − n ⋅ ∫ T 0 T l j A D ⋅ 1 R T exp ( − E R T ) d T       ( n ≠ 1 ) (24) ∂ Q ∂ E = ∑ l = 1 L 0 ∑ j = 1 J 0 ( − 2 ( X 1 l j − X l j ) ⋅ ∂ X l j ∂ E )                   = ∑ l = 1 L 0 ∑ j = 1 J 0 ( 2 ( X 1 l j − X l j ) ⋅ exp ( − ∫ T 0 T l j A D exp ( − E R T ) d T ) ⋅ ∫ T 0 T l j A D ⋅ ( 1 R T ) ⋅ exp ( − E R T ) d T )         ( n = 1 ) (25)

The partial derivative of the objective function to the frequency factor as follows: ∂ Q ∂ A = ∑ l = 1 L 0 ∑ j = 1 J 0 ( − 2 ( X 1 l j − X l j ) ⋅ ∂ X l j ∂ A ) (26) where ∂ X l j ∂ A = ∂ ( 1 − ( 1 − ( 1 − n ) ∫ T 0 T l j A D exp ( − E R T ) d T ) 1 1 − n ) ∂ A O = ( 1 − ( 1 − n ) ∫ T 0 T l j A D exp ( − E R T ) d T ) n 1 − n ⋅ ∫ T 0 T l j 1 D exp ( − E R T ) d T           ( n ≠ 1 ) (27) ∂ X l j ∂ A = ∂ ( 1 − exp ( − ∫ T 0 T l j A D exp ( − E R T ) d T ) ) ∂ A O = exp ( − ∫ T 0 T l j A D exp ( − E R T ) d T ) ⋅ ( ∫ T 0 T l j 1 D exp ( − E R T ) d T )         ( n = 1 ) (28)

The partial derivative of the objective function to the order of the reaction is as follows: ∂ Q ∂ n = ∑ l = 1 L 0 ∑ j = 1 J 0 ( − 2 ( X 1 l j − X l j ) ⋅ ∂ X l j ∂ n ) (29) where ∂ X l j ∂ n = − [ 1 − ( 1 − n ) ∫ T 0 T l j A D exp ( − E R T ) d T ] 1 1 − n          ⋅ [ ln ( 1 − ( 1 − n ) ) ∫ T 0 T l j A D exp ( − E R T ) d T) ⋅ ( 1 1 − n ) 2 + ∫ T 0 T l j A D exp ( − E R T ) d T ⋅ ( 1 − ( 1 − n ) ∫ T 0 T l j A D exp ( − E R T ) d T ) − 1 ⋅ 1 1 − n ]           ( n ≠ 1 ) (30)

When n = 1, there is no need to calculate the partial derivative of the objective function to the order of the reaction.

The partial derivative of the penalty term is as follows: ∂ Q ∂ A = 2 ⋅ min ( 0 , A ) (31) ∂ Q ∂ E = 2 ⋅ min ( 0 , E ) (32) ∂ G ∂ n = 2 ⋅ min ( 0 , n ) (33)

After obtaining the partial derivatives of A, E, and n, the subsequent optimization of the activation energy at the minimum point of the objective function, the frequency factor, and the reaction order are the same as for the parallel first–order reaction model.