Paleostress Analysis from Calcite Twins at the Longshan Dome (Central Hunan, South China): Mesozoic Mega-Fold Superimposition in the Reworked Continent

: It is generally accepted that during the Mesozoic NE − NNE-trending folds overprinted E − W-trending folds to form the Longshan dome in the central South China continent, although the interference map does not tell the relative ages of the fold sets. In an effort to deepen our understanding of the process of reworking the continent, paleostress analysis using calcite twins was carried out in this study to verify or falsify this model. Ten limestone samples were collected from Upper-Paleozoic limestones on the ﬂanks of the dome and were measured using the universal stage for calcite e-twins. E-twins in the samples are divisible into two kinds, thick ( ≥ 1 µ m) and thin (<1 µ m), indicative of relatively higher and lower deformation temperatures, respectively. Stress estimates obtained using the improved version of Shan et al.’s (2019) method were grouped into two layer-parallel shortening (LPS) subsets and three non-LPS subsets. These subsets comprise four tectonic regimes: NWW − SEE compression (LPS1 and non-LPS1), NNE − SSW compression (LPS2 and non-LPS2), NW − SE extension (non-LPS3a) and NNE − SSW extension (non-LPS3b). They were further arranged in a temperature-decreasing order to establish a complex deformation sequence of the study area. In the sequence NE − NNE-trending folds have an older age than E − W-trending folds, something different from the model. The approximately N − S regional compression responsible for the former folds should have a profound effect on the intensely deformed continent, something ignored in earlier work.


Introduction
Kilometer-scale superimposed folds are common in ancient accretionary [1][2][3][4] and collisional [5,6] orogens. They are generally diagnostic of a complex plate convergence and, hence, of practical significance in understanding the deformation history of the orogens. They are however not peculiar to the orogens, for example, the classic dome-basin structure or Ramsay's [7] type 1 cross fold formed in the central portion of the South China continent (Figure 1a), when the continent had undergone intense deformation and intense magmatism during the Mesozoic [8][9][10][11][12][13]. This interference was caused by two fold sets: NE−NNE-trending tight folds and nearly E−W-trending open folds (Figure 1b).
Within the dome-basin structure the oval-shaped Longshan tectonic dome (Figure 2a) has received the attention of many Chinese geologists ever since [14][15][16][17][18][19][20]. Although the map of fold interference does not tell the relative age of two fold sets, previous workers [14][15][16][17][18][19][20] have tried to correlate the two fold sets with the parallel fold sets on the continental margin. They reached a unanimous agreement that the NE−NNE fold set should have a younger age than the E fold set. What they argued about is when and how the fold sets workers [14][15][16][17][18][19][20] have tried to correlate the two fold sets with the parallel fold sets on the continental margin. They reached a unanimous agreement that the NE−NNE fold set should have a younger age than the E fold set. What they argued about is when and how the fold sets formed. The mainstream viewpoint relates the E and NE−NNE fold sets to the collision between the North China and South China continents in the Middle−Late Triassic and the westward subduction of the Izanagi oceanic plate in the time of Middle−Late Jurassic to Early Cretaceous [15,18,19], respectively. However, in the belief that this collision-related deformation is confined to the northern continental margin [21], Li et al. [17] ascribed the fold sets to the northeastward subduction of the Paleo-Tethys oceanic plate and subsequently the collision of the Indosinian block to the clockwise rotating South China continent, respectively, during the Late Permian to Triassic. The principal reason for this inconsistency is, as described below, the inhomogeneous, polyphase and even repetitious nature of regional deformation in the South China during the Mesozoic.    [22]). The dome flank is divisable into four parts, eastern, southern, western and northern, for the convenience of this study. Equal-area, lower hemisphere projection (d-o) of beddings, slickenlines on bedding-parallel faults, and cleavages measured in the Upper Paleozoic and Triassic in each part. Each slickenline is represented by a slickenline-parallel short line at the bedding pole, with or without a directed arrow toward the relative slip movement of the overlying bed. Shown in (a) are locations of limestone samples collected in this study.  [22]). The dome flank is divisable into four parts, eastern, southern, western and northern, for the convenience of this study. Equal-area, lower hemisphere projection (d-o) of beddings, slickenlines on bedding-parallel faults, and cleavages measured in the Upper Paleozoic and Triassic in each part. Each slickenline is represented by a slickenline-parallel short line at the bedding pole, with or without a directed arrow toward the relative slip movement of the overlying bed. Shown in (a) are locations of limestone samples collected in this study.
In South China, severe vegetation coverage makes it difficult to use secondary deformation structures to determine the relative ages of the fold sets that made the Longshan dome. In these circumstances, it is necessary to look for an appropriate kind of secondary deformation structure that we can correlate from site to site within the cross fold. Cal- cite e-twinning is perhaps the best kind of deformation mechanism, as it is sensitive to stress [23,24], and because limestone is abundant in the platformal cover that was involved in the cross fold. Paleostress analysis using calcite e-twin data has been widely used to reveal the structural evolution of the fold belts [25][26][27][28][29].
This communication applied paleostress analysis to the Longshan tectonic dome to better understand the mega-fold superimposition in the interior of the South China continent during the Mesozoic. For this purpose we conducted field observations of secondary deformation structures on the flanks of the dome, collected oriented samples of sparitic limestone of the platformal cover, and measured calcite e-twins under the universal stage. Meanwhile, Shan et al.'s [30] inversion method for polyphase twin data was improved and applied to the measured datasets.

History
In the Late Neoproterozoic the primary South China continent disintegrated into two principal portions: the Yangtze platform in the west and the less stable Cathaysian block in the east (Figure 1a). These portions merged as a result of something like intracontinental orogeny [31][32][33]. This orogeny affected the Yangtze platform by means of uplift and erosion but produced intense deformation, low-grade regional metamorphism and the emplacement of post-orogenic granites in the Cathaysian block [33,34]. A ca. 8 km thick succession of Late Neoproterozoic-Silurian sedimentary rocks deposited on the continental shelf and slope was thrusted and folded throughout the block. Regional folds have a generally E−W trend in the central Hunan [35]. Consequently the platform regime prevailed in the remerged continent and did not terminate until the Early Triassic [31].
The early Mesozoic marks the onset of regional shortening in a far-or near-field response to the multiple plate convergence around the South China continent. In the Triassic the continent collided with the Indosinian block to the southwest [28,31,36] and with the North China continent to the north [37,38] and had been subducted by the Izanagi or older oceanic plate to the east [10]. This caused the contraction of the continent, as a whole, which differed in pattern between the portions. The fault-and-thrust belts formed as narrow forelands on the platform margin and were hardly associated with magmatism [37,38]. In contrast, one or more phases of deformation happened in most, if not all, of the Cathaysian block, in a close relation with the emplacement of numerous granites or granodiorites [10] and hydrothermal mineralization, for example, many Au or Au-Sb deposits in the Hunan province [39].
A more intense and more extensive contraction occurred in the Middle or Late Jurassic−Early Cretaceous. The pre-existing foreland fold-and-thrust belts on the margins were reactivated to propagate towards the inland. This resulted from a combination of the continuous westward subduction of the Izanagi plate to the east [40,41], the southward collision of the Siberian plate with the Sino-Korea plate far to the north [11,42], and the collision between the Lhasa and the Qiangtang blocks to the west [43,44]. The first tectonic process played the most important role in reworking the South China continent by producing a 1300 km wide tectono-magmatic belt in the central and eastern portions. Apart from the intact Shichuan basin in the west, all the portions of the South China continent were involved in structural deformation to varying degrees.
In the Cretaceous, margin-perpendicular regional extension took place throughout the South China continent, due to possible tectonic processes such as the collapse of the Early-Mesozoic orogen [45], the high-angle subduction of the Pacific oceanic plate [9,40] and the roll-back of the subduction zone [46,47]. Numerous NE− or NNE−trending extensional basins were produced and infilled with a several-kilometers thick sequence of fluvial and lacustrine sediments and, to variable degrees, volcanic rocks in the central and eastern portions.

Macroscopic Deformation Structures
Two sets of regional folds are present in the central Hunan province: longitudinal and transverse with respect to the fold belt ( Figure 1b). The longitudinal folds trend toward the northeast or north-northeast and are upright, inclined or locally overturned. They have a straight or curved fold axial trace, with a wavelength of ca. 20−30 km. The anticlines are closer than the synclines at most places. They become more open at a few places, for example, where the basement-involved tectonic domes formed. These cases seem to represent thin-skinned and thick-skinned tectonics, respectively. In contrast, the transverse folds are less developed. They trend toward the east-southeast, and occasionally the northeast. They are gentle and nearly upright and have a wavelength of ac. 90 km, about three times larger than the longitudinal folds. There are two kinds of fold interference pattern in the region: dome-and-basin structure or Ramsay's [7] type 1 and bending of longitudinal folds as a transverse fold or Ramsay's [7] type 2. Theoretically speaking, in the case of such fold sets, the pattern kind is controlled by the pitch of the transverse fold axis [48,49]. The former pattern is gradually replaced by the latter pattern when the transverse fold has an increasing pitch from 0 • -90 • .

1.
Dome-and-basin structure. Tectonic domes occur in an E-W or NE-SW-trending chain. They are in a circular or elliptical shape, with a size of 20-30 km. They involve both the Devonian-Lower Triassic platformal cover and the Upper Proterozoic-Silurian folded basement. Between the neighboring domes is a NNE-or NE-trending syncline, close and long-axis. Tectonic basins are much less developed and are characterized by NNE-or NE-trending short-axis synclines, 7-15 km in width and 20-60 km in length. The oval-shaped Longshan dome to the southeast of Lianyuan County is the most typical tectonic dome. It is represented by a gentle anticline with a half wavelength of 35 km in the NE-SW direction ( Figure 2b) and a box anticline with a half wavelength of 20 km in the NW-SE direction (Figure 2c). It is generally accepted that the dome is a good example of Ramsay's [7] type 1 cross folds and formed during the Mesozoic [15,[17][18][19][20]. Theoretically speaking, the map of such a cross fold does not tell the relative chronology of the aforementioned fold sets. Even now this issue is not well resolved because, as described below, there is a lack of reliable field observations about the relationship between the sets.

2.
Bending of longitudinal folds as a transverse fold. In plan-view the longitudinal fold set is bent into a few open arcuate folds with a roughly E−W axial trace ( Figure 1b). Such an interference map is representative of Ramsay's [7] type 2 cross fold and, hence, indicative of a younger age of the transverse folds than the longitudinal folds. This is however ignored by existing models that tentatively ascribe in various ways the transverse folds to WNW−ESE or NW−SE regional compression responsible for the longitudinal folds. They include the local disturbance around some pre-existing granite intrusive bodies of the regional deformation field [50,51], the propagation of arcuate thrust faults [16] and the fault-bend folds of NE-striking sinistral strike-slip faults [52]. The fact that all of the intrusive bodies, thrusts or strike-slip faults have a much smaller size than the longitudinal folds ( Figure 1b) is strongly suggestive of no genetic relationship between them. These models have not been verified ever since, and they will not be discussed below to save the room in this paper.
In addition, NE−NNE-trending regional faults have a complex history of formation and reactivation [22]. They may be reverse, strike-slip, normal or in combination. The reverse faults are generally parallel to the longitudinal fold set (Figure 1b). Some of them had been subsequently reactivated as normal faults controlling Cretaceous half grabens. At places the mylonitization of felsic rocks was induced on the contacts of the footwalls of the boundary faults. 40 Ar/ 39 Ar dating of newly-formed muscovite, sericite, and biotite in these ductile shear zones yields a deformation age of 195−217 Ma for reverse faults [53] and 211.7 ± 1.5 Ma [18] for a dextrally strike-slip fault in the Xuefengshan structural belt in the west, and 120−140 Ma and 85−100 Ma for normal faults in the southeast of Xinning town [54].

Mesoscopic Deformation Structures
Present on the flank of the Longshan dome are a variety of secondary deformation structures including slaty cleavages (Figure 3a

Western Flank
Secondary folds, centimeters ( Figure 3c) to meters or a little more in wavelength, are common in Carboniferous and Permian mudstones or shales on this flank. They possess a horizontal NNE-trending fold axis, as is defined on the stereonet by bedding measurements (Figure 2d), in accordance with a close regional syncline on the west. Associated with the syncline is a major set of NNE-striking vertical or sub-vertical slaty cleavages in mudstones ( Figure 2e). Pre-existing bedding-parallel or -oblique calcite veins were strongly folded to accommodate the amount of shortening absorbed by slaty cleavages (Figure 3c).
However, a few slaty cleavages gently dipping towards the northeast (Figure 2e) are present at the place where the syncline axis was bent. They are probably attributed to the transverse folding. Unfortunately no superposition of cleavages is found in the field.
As a flexural slip, slickenlines on bedding-parallel faults are reliable recorders of multiple folding. They have two major sets, NWW−SEE and NNE−SSW (Figure 2f), which seem to correspond to the longitudinal and transverse folds, respectively. It is therefore crucial but difficult to determine the relative ages of two or more sets of slickenlines on a certain bedding-parallel fault at outcrops. Only at the location of sample H34 is there an evident relationship between the sets of slickenlines, dip-parallel and strike-parallel, on a certain bedding-parallel fault. The former set turns gradually into the latter set ( Figure 3d). Given the close relationship between the former set and the longitudinal folds, this strongly suggests that the longitudinal folds have an earlier age than the transverse folds.

Northern Flank
This flank was involved in complex deformation that accommodated the deformation incompatibility between the tectonic domes and basins. For this reason, two vertical great circles, north-northeast-striking and northwest-striking, are poorly defined by bedding measurements (Figure 2g). They represent west-northwest-and northeast-trending folds, respectively. Considering a vertical, northwest-striking great circle of slaty cleavages (Figure 2h), the absence of cleavages of west-northwest-trending folds implies that the folds have a younger age.
Locally exposed near Baima village are spaced cleavages, a set of nearly beddingperpendicular, relatively densely-spaced stylolites, in a gently dipping thick bed of Upper Devonian limestone (Figure 3e). They are associated with and perpendicular or subperpendicular to lenticular, V-shaped or tabular calcite veins or veinlets. All of them reflect a local ESE−WNW shortening, generally in correspondence with the longitudinally folding. These structures are offset by a nearly bedding-parallel fault bearing a set of east-pitching slickenlines. Calcite fibers and steps on the fault indicate an oblique reverse slip sense. As such the fault is most likely to form in a different stress state from the spaced cleavages. It is however unknown whether the fault is related to the longitudinally folding as a result of progressive deformation, the transversely folding or other possible tectonic processes.

Eastern Flank
On this flank there is a north-northeast-striking vertical great circle of bedding measurements (Figure 2j), as similar to that on the western flank. Slaty cleavages, although in a small number, do not cluster around a certain great circle (Figure 2k). Half of them strike towards the east or southeast and are probably associated with the transverse folds.

Southern Flank
There is a southwest-dipping great circle of bedding measurements (Figure 2m). Slaty cleavages have high dips and strikes generally toward the northeast, as similar to the ones on the northern flank ( Figure 2n).
On the ridge near Hetianpu village there are two principal sets of calcite veins in a subhorizontal thick bed of Lower Carboniferous limestone (Figure 3f). A north-striking set comprises tabular or, in an en echelon pattern, lenticular veins. Short stylolites are present in the overlaps between the veins and are generally perpendicular to the veins. This set is cut through by an east-striking set of tabular calcite veins. Veins of each set have a serrated contact in some portions, owing to pressure dissolution under the subsequent perpendicular or subperpendicular shortening. Both sets are cut through or abut by two minor sets of tabular veins or veinlets, northeast-and northwest-striking, respectively. These scenarios give a variation in direction of local compression from S−N to E−W to SE−NW and SW−NE with time.
However, at another location on the flank (Figure 3g) a different deformation history is established: the local compression direction varies from E−W to SE−NW to N−S with the time. The difference in deformation history between the neighboring sites means that the occurrence of parallel veining makes it difficult to unravel the complexity in time and space of deformation. The east-and north-trending sets involve two low-temperature deformation mechanisms: tensile fracturing and, insignificantly, dissolution (stylolites) − precipitation (veins). Considering that this combination works at a relatively high temperature [55], some if not all of the veins of each set formed at a relatively deep depth and therefore must have the oldest age.

Methodology
In a certain calcite crystal e-twining or mechanical twinning occurs by gliding along a certain e-plane toward the C-axis when the resolved shear stress (RSS) on the crystallographic plane equals or exceeds a certain critical value of about 10 MPa [56][57][58][59]. This inequality is the basis of inversion of calcite monophase e-twins for stress [60][61][62][63][64][65][66]. According to the inequality, a relative or absolute deviatoric stress tensor can be defined by a population of monophase twin data, if they have a relatively large number and are widespread in orientation space. The relative deviatoric stress tensor is also called the reduced stress, from which one knows only the principal stress axes and the stress ratio.
Polyphase e-twins are more common in nature due to the variation in time of the tectonic stress field during the geological history. How to invert them for polyphase stresses has received rigorous research during the decades. There are two categories of existing inversion methods that are strictly based on the inequality. The first category includes Etchecopar's [64] method and its modifications [67,68]. Etchecopar's [64] method sequentially applies his proposed method for monophase twin data to a certain population of polyphase twin data until the remaining twinned e-planes are insufficient to define a new stress. As such the driving stresses responsible for the data are required to have as much large differences as possible, which limits the wide application of the method. This inherent drawback has been tried to be overcome in several ways, for example, the re-separation of twinned e-planes with respect to the previously obtained stresses for a newly obtained stress [67], and the separation of all twinned e-planes before inversion [68].
The second category is the solution of an optimization problem that depends on the compatibility or incompatibility of twinned or untwinned e-planes and the shear stresses along the gliding directions on compatible twinned e-planes under the unknown stresses [30]. For example, Shan et al. [30] take into account two criteria: minimizing the number of incompatible e-planes (twinned and untwinned) shared by the unknown stresses, N SI , and maximizing the sum of shear stresses of compatible twinned e-planes. Their minimization problem has the following object function: where N T is the number of twinned e-planes, K is the number of tectonic phases, M i is the number of twinned e-planes for the i-th stress, σ (i) is the i-th modified reduced stress vector, t I(j) is the I(j)-th modified datum vector, and I(j) is the position in the e-plane sequence of the j-th compatible twinned e-plane for the i-th stress. In order to facilitate the search of a reduced stress tensor in five-dimensional stress space, a constant critical resolved shear stress value, 10 MPa, is used to rescale the tensor. The real deviatoric stress tensor can be restored with the knowledge of the real critical value. The evolutionary algorithm was adopted to solve the above optimization problem. Application of this method to a series of two-phase artificial twin datasets indicates the relatively high accuracy of the stress estimate, even in the case of datasets generated under the reasonably similar imposed stresses, but the relatively low accuracy for stress difference [30].
In this study, maximizing the number of compatible untwinned e-planes shared by the unknown stresses, N SC , is chosen as a third criterion about the unknown stresses. The above object function is improved in the following way: Some of the artificial two-phase e-twin datasets made by Shan et al. [30] were used to check the validity of this object function. The first-phase subsets, each comprising 20, 40 or 60 twinned calcite crystals, are generated under the Andersonian stress that has a stress ratio of 0.5 and a maximum differential stress of 60 MPa. They are rotated around the maximum, intermediate and minimum principal stress axes, respectively, with an increment of 5 • to have the second-phase subsets. Before combining two of the subsets into a two-phase equally-sized dataset we need to calculate whether each untwinned e-plane of one subset is activated or not under the imposed stress with respect to another subset (see [30] for more information). For each dataset, the estimated stresses obtained using this improved method are correlated with the imposed stresses. Two parameters are adopted here to characterize the difference between the estimated stresses and their corresponding imposed stresses [30]: The maximum stress difference (MSD). This is the maximum of differences between the estimated stresses and their corresponding imposed stresses. The stress difference (SD; [69]) describes the difference between the reduced stresses. It is equal to the Euclidean distance between the stress vectors in reduced and partially distorted stress space [70]. It ranges from 0-2, for the strictly similar to strictly different reduced stresses.

2.
The maximum relative error of the maximum differential stress (MREMDS). This is the maximum relative error of the maximum differential stresses of the estimated stress tensors to the imposed stress tensors. As previously described, the imposed stress tensors are initiated with a maximum differential stress of 60 MPa.
Shown in Figures 4 and 5 are results of applying the improved method and, for comparison, the primary method to the datasets. Both methods have a strong tendency to give an increasing accuracy of reduced stress with the ascent of rotation angle and crystal number (Figure 4). As a whole, the improved method has a smaller MSD value and hence a higher accuracy than the primary method, particularly in the case of rotation around the intermediate and minimum principal stress axes, respectively. The MREMDS is generally smaller in value for the improved method than for the primary method in all cases ( Figure 5). These illustrate the feasibility of the improved method.
sequence of the j-th compatible twinned e-plane for the i-th stress. In order to facilitate the search of a reduced stress tensor in five-dimensional stress space, a constant critical resolved shear stress value, 10 MPa, is used to rescale the tensor. The real deviatoric stress tensor can be restored with the knowledge of the real critical value. The evolutionary algorithm was adopted to solve the above optimization problem. Application of this method to a series of two-phase artificial twin datasets indicates the relatively high accuracy of the stress estimate, even in the case of datasets generated under the reasonably similar imposed stresses, but the relatively low accuracy for stress difference [30].  In this study, maximizing the number of compatible untwinned e-planes shared by the unknown stresses, NSC, is chosen as a third criterion about the unknown stresses. The above object function is improved in the following way: Some of the artificial two-phase e-twin datasets made by Shan et al. [30] were used to check the validity of this object function. The first-phase subsets, each comprising 20, 40 or 60 twinned calcite crystals, are generated under the Andersonian stress that has a stress ratio of 0.5 and a maximum differential stress of 60 MPa. They are rotated around the maximum, intermediate and minimum principal stress axes, respectively, with an increment of 5° to have the second-phase subsets. Before combining two of the subsets into a two-phase equally-sized dataset we need to calculate whether each untwinned e-plane of one subset is activated or not under the imposed stress with respect to another subset (see [30] for more information). For each dataset, the estimated stresses obtained using this improved method are correlated with the imposed stresses. Two parameters are adopted here to characterize the difference between the estimated stresses and their corresponding imposed stresses [30]: 1. The maximum stress difference (MSD). This is the maximum of differences between the estimated stresses and their corresponding imposed stresses. The stress difference (SD; [69]) describes the difference between the reduced stresses. It is equal to the Euclidean distance between the stress vectors in reduced and partially distorted stress space [70]. It ranges from 0-2, for the strictly similar to strictly different reduced stresses. 2. The maximum relative error of the maximum differential stress (MREMDS). This is the maximum relative error of the maximum differential stresses of the estimated stress tensors to the imposed stress tensors. As previously described, the imposed stress tensors are initiated with a maximum differential stress of 60 MPa. In addition, the improved method yields a MREMDS range of 0-1 for a MSD value of <0.3 and a crystal number of 60 ( Figure 6). That is to say, when the estimates of reduced stress are acceptable, there is still a very large range or uncertainty of their maximum differential stresses. Figures 4-5 are results of applying the improved method and, for comparison, the primary method to the datasets. Both methods have a strong tendency to give an increasing accuracy of reduced stress with the ascent of rotation angle and crystal number (Figure 4). As a whole, the improved method has a smaller MSD value and hence a higher accuracy than the primary method, particularly in the case of rotation around the intermediate and minimum principal stress axes, respectively. The MREMDS is generally smaller in value for the improved method than for the primary method in all cases ( Figure 5). These illustrate the feasibility of the improved method.

Shown in
In addition, the improved method yields a MREMDS range of 0-1 for a MSD value of < 0.3 and a crystal number of 60 ( Figure 6). That is to say, when the estimates of reduced stress are acceptable, there is still a very large range or uncertainty of their maximum differential stresses.

Samples and Measurements
In this study ten oriented samples (Table 1) were collected from the Middle Devonian and Carboniferous sparitic limestones on the flanks of the Longshan dome and, for comparison, on the eastern and southern limbs of the triangular tectonic basin abutting the northern flank (Figure 2a). All sampling locations are away from outcrop-scale interlayer fold(s) and fault(s). Each sample was cut into three mutually perpendicular thin sections, perpendicular to the reference plane along the dip and strike, respectively, and parallel to the reference plane.
Under the microscope, nearly all calcite grains are intensely involved in twinning, with an average density of ca. 50 twins per mm. They host one or, more significantly, two sets of twins. Twins in each set are thin, thick or hybrid (Figure 7a-b). Thin twins are < 1 μm thick and usually straight. Thick twins are generally 1−5 μm thick and rarely 10 μm or more. They have straight (Figure 7b), lenticular and curved (Figure 7c-d) shapes. Thin twins are parallel to thick twins (Figure 7b), crosscut them (Figures 7a-b) or stop at them (Figures 7b-c). These relationships indicate that thin twins have a similar age to thick twins, as a whole. Nowadays there are two different mechanisms about the formation of thick twins: the increase in strain [71] and the relatively high temperature [72,73]. Considering that the curved shape of some thick twins (Figure 7b-c) is diagnostic of the relatively high temperature [72], the measured twins in each sample are separated into two groups, thin-twin and thick-twin, for estimation of stress.

Samples and Measurements
In this study ten oriented samples (Table 1) were collected from the Middle Devonian and Carboniferous sparitic limestones on the flanks of the Longshan dome and, for comparison, on the eastern and southern limbs of the triangular tectonic basin abutting the northern flank (Figure 2a). All sampling locations are away from outcrop-scale interlayer fold(s) and fault(s). Each sample was cut into three mutually perpendicular thin sections, perpendicular to the reference plane along the dip and strike, respectively, and parallel to the reference plane. Under the microscope, nearly all calcite grains are intensely involved in twinning, with an average density of ca. 50 twins per mm. They host one or, more significantly, two sets of twins. Twins in each set are thin, thick or hybrid (Figure 7a-b). Thin twins are < 1 µm thick and usually straight. Thick twins are generally 1−5 µm thick and rarely 10 µm or more. They have straight (Figure 7b), lenticular and curved (Figure 7c-d) shapes. Thin twins are parallel to thick twins (Figure 7b), crosscut them (Figure 7a-b) or stop at them (Figure 7b-c). These relationships indicate that thin twins have a similar age to thick twins, as a whole. Nowadays there are two different mechanisms about the formation of thick twins: the increase in strain [71] and the relatively high temperature [72,73]. Considering that the curved shape of some thick twins (Figure 7b-c) is diagnostic of the relatively high temperature [72], the measured twins in each sample are separated into two groups, thin-twin and thick-twin, for estimation of stress.   In each thin section, the crystal C-axes (Figure 8a,d) and twinned e-planes (Figure 8b,e) of more than 30 twinned calcite grains were measured under the universal stage. Each measurement was calibrated according to the crystallographic relationship and then used to calculate the untwinned e-plane(s). The calculated untwinned e-planes (Figure 8c,f) at an angle of <30 • with the thin section are in fact uncertain under the universal stage and are, for safety, discarded. The measured C-axes for each sample tend to have an even distribution (Figure 8a,d), as required in stress inversion.

Tables 2 and 3 list the results obtained by applying the improved version of Shan et al.'s [30]
inversion method to the measured datasets. The optimal phase number of each dataset is determined under the following restraints: a relatively small value of the minimized object function, a very small number of the common incompatible twinned e-planes, and the absence of a new different stress for the next phase number. As such thin-twin and thick-twin datasets are separated into two or three subsets. An exception is the thick-twin dataset of sample H27 that has an optimal number of 4. By comparing the best estimates in the current and untilted states we classified them into two types: layer-parallel shortening (LPS) and non-layer-parallel shortening (non-LPS). For LPS, the maximum principal stress axes (PSA) are more horizontal or less upright in the untilted state than in the current state. Afterwards, according to the similarity of the maximum PSAs, the stress estimates are further grouped into five subsets, two LPS and three non-LPS (Figure 9). There is no great difference in subsets between the thick-twin estimates (TKE) and the thin-twin estimates (TNE).     The first LPS subset (LPS1) represents a WNW-ESE compression and has an average maximum PSA of 101 • (Figure 10a). Most of TKEs in the untilted state have a sub-vertical minimum PSA and a relatively wide range of the maximum PSA from NE to SE. The TNEs in this state have a high-angle intermediate PSA, and the maximum PSAs concentrate in the nearly east-west direction. A possible reason for this difference is the slight change of the tectonic stress field during the twinning.
The second LPS subset (LPS2) corresponds to a NNE-SSW compression with an average maximum PSA of 199 • (Figure 10a). All TKEs and TNEs in the untilted state have a nearly horizontal maximum PSA and a high-angle or sub-vertical minimum PSA. This subset has a greater concentration of the maximum PSAs than the previous subset.
The first non-LPS subset (Non-LPS1) is characterized by a NW-SE compression, with an average maximum PSA of 120 • (Figure 10b). An exception to this tectonic regime is the strike-slip regime of the TNE of sample H24. The presence of this subset in all samples suggests that the most intense deformation occurred in this region since the early Mesozoic. Besides, a few seemingly spurious estimates, such as TKEs for samples J90 and H30 and TNE for sample H27, have a medium-angle maximum PSA. It would be better to ascribe these to the syn-folding compression with respect to the non-LPS1. This was not endertaken for simplicity.
for thick twins (Figure 11). There is a larger mean corrected MDS for thin twins than for thick twins in the LPS and Non-LPS1 subsets, and for thick twins in the Non-LPS2 and Non-LPS3 subsets. As a whole the former subsets, in particular the LPS subsets, have a narrower corrected-MDS range than the latter subsets. For thin twins in the LPS1 subset ( Figure 12a) and for both thin twins and thick twins in LPS2 subset (Figure 12b), the corrected MDSs have a good positively-linear correlation to the buried depths, with slope values of 17.7-19.7 MPa/km. This illustrates the validity of the improved method, although it yields, as shown in Figure 6, a generally low accuracy of the MDS.   The second non-LPS subset (Non-LPS2) is dominated by a N-S or NE-SW compression that has an average maximum PSA of 194 • (Figure 10b). Only the TNE of sample H30 has a strike-slip stress regime. There are however two clusters, NNW-SSE and NE-SW, of the maximum PSAs. This may be indicative of two tectonic phases, but we would like to relate it to the spatial variation of the tectonic stress field since this subset has a relatively small number of estimates, 7.
The third non-LPS subset reflects the extensional regime. It has two sub-subsets, NNE-SSW extension (non-LPS3a) and NW-SE extension (non-LPS3b), with an average minimum PSA of 20 • and 116 • , respectively (Figure 10c). Similar to the other subsets, this subset is present on the four flanks of the dome. This suggests a regional extent of the extension.
Interestingly, both the TKEs and the TNEs even belong to a similar subset for some samples, for example, the LPS subsets of sample H20 in the northern flank ( Figure 9). This corresponds with the aforementioned observations about the similarity in age between thin twins and thick twins. In this manner thin and thick twins seem not diagnostic of relatively low and relatively high temperature, respectively, as found in the axial shortening experiment of single calcite crystal [58]. However, as previously described, slaty cleavages in the cover attests to the occurrence of low-grade regional metamorphism during the formation. Such a relatively high temperature is in accordance with the curved shape of some thick twins (Figure 7c-d; [72]). As a compromise, thin twins are considered to occur in a wider range of deformation temperature.
In addition, the maximum differential stresses (MDS) obtained from the samples were corrected (Figure 11), according to Rocher et al.'s [67] empirical correlation between the mean calcite size and the critical resolved shear stress. As this correlation is based on samples subjected to very small strain, ca. 1% [58,67], it may be inapplicable to many of our samples that had undergone relatively large strain in terms of densely-spaced twins. In light of an apparently smaller strain for thin twins than for thick twins, the corrected MDSs would be more accurate for thin twins. This is in agreement with the fact that the mean corrected MDSs have a narrower range for thin twins than for thick twins ( Figure 11). There is a larger mean corrected MDS for thin twins than for thick twins in the LPS and Non-LPS1 subsets, and for thick twins in the Non-LPS2 and Non-LPS3 subsets. As a whole the former subsets, in particular the LPS subsets, have a narrower corrected-MDS range than the latter subsets. For thin twins in the LPS1 subset ( Figure 12a) and for both thin twins and thick twins in LPS2 subset (Figure 12b), the corrected MDSs have a good positively-linear correlation to the buried depths, with slope values of 17.7-19.7 MPa/km. This illustrates the validity of the improved method, although it yields, as shown in Figure 6, a generally low accuracy of the MDS.
Geosciences 2021, 11, x FOR PEER REVIEW 17 of 2 Figure 10. Lower-hemisphere, equal-area projection of the maximum principal stress axes of each compression subset (a-b) and the minimum principal stress axes of the extension subset (c). For each subset, the mean direction and angular dispersion of these axes are displayed as a plus symbol and a small circle centered at it. Figure 11. The range and mean value of the corrected (see [67]) maximum differential stresse (MDS) of each subset.

Interpretation
Our stress results (Figure 9) reveal the complex evolution of the tectonic stress fiel in the study area since the Mesozoic. This complexity is represented by four tectonic re gimes: NWW−SEE compression (LPS1 and non-LPS1), NNE−SSW compression (LPS and non-LPS2), NW−SE extension (non-LPS3a) and NNE−SSW extension (non-LPS3b The compressional regimes are nearly perpendicular to the transverse and longitudina Figure 11. The range and mean value of the corrected (see [67]) maximum differential stresses (MDS) of each subset. Regional compression with respect to the second compressional regime took place in the Middle-Late Triassic and again in the Late Jurassic-Early Cretaceous, as a result of the long-standing northwestward or westward subduction of the Izanagi and subsequently Pacific oceanic plates. A majority of the South China continent was involved in deformation and magmatism, presumably due to the flat-slab subduction [10]. In the study area, this Triassic compression is evidenced by the intrusion of a few NW-striking felsic dykes with a Zircon U-Pb age of 217−220 Ma on the eastern and northern flanks [75], and an angular unconformity between the Jurassic and Triassic (Figure 1).
The first extensional regime (non-LPS3a) belongs to NW−SE regional extension on the eastern margin of the Eurasian continent in the Cretaceous. This regional extension is characterized by a spread of NE-trending grabens or half-grabens in the central and eastern portions of the South China continent [76,77]. It is attributed to a variety of possible tectonic processes such as the collapse of the Early-Mesozoic marginal orogen [45], the high-angle subduction of the Pacific oceanic plate [9,40] and the roll-back of the subduction zone [46,47].
The second extensional regime (non-LPS3b) seems not to have a counterpart of regional extension. More about this will be discussed below.
Furthermore, in the ideal upper crust, where deformation is accommodated by the reactivation of pre-existing faults, the minimum differential stress gradient (MDSG) under the thrust regime is calculated in the below expression [78]: where and are the maximum and minimum principal stresses, z is the buried depth, is the minimum ⁄ , is the rock density, g is the gravitational acceleration, and is the pore-fluid factor. is the hydrostatic pressure gradient at the depth (HPG). Under the nearly hydrostatic condition the friction coefficient reaches ca. 0.75, yielding a value of 4. Then for a typical value of 0.4 [78,79] and the normal HPG,

Interpretation
Our stress results (Figure 9) reveal the complex evolution of the tectonic stress field in the study area since the Mesozoic. This complexity is represented by four tectonic regimes: NWW−SEE compression (LPS1 and non-LPS1), NNE−SSW compression (LPS2 and non-LPS2), NW−SE extension (non-LPS3a) and NNE−SSW extension (non-LPS3b). The compressional regimes are nearly perpendicular to the transverse and longitudinal fold sets in the region (Figure 1b), respectively. In this manner both regimes have a regional extent. Their presence in the pre-folding stage indicates a much complex deformation history of the Longshan dome, in contrast with the generally accepted model, where transverse folds are simply overprinted by longitudinal folds [15][16][17][18][19].
Regional compression with respect to the first compressional regime happened in the Middle-Late Triassic, due to the collision between the South China and North China continents [37,38], and re-occurred in the Late Jurassic-Early Cretaceous, principally owing to the closure of the Mongolia-Okhotsk Ocean in the long distance to the north [11]. As a result a fold-and-thrust belt formed [37,38,74] and was subsequently reactivated [11,42] on the northern margin of the South China continent.
Regional compression with respect to the second compressional regime took place in the Middle-Late Triassic and again in the Late Jurassic-Early Cretaceous, as a result of the long-standing northwestward or westward subduction of the Izanagi and subsequently Pacific oceanic plates. A majority of the South China continent was involved in deformation and magmatism, presumably due to the flat-slab subduction [10]. In the study area, this Triassic compression is evidenced by the intrusion of a few NW-striking felsic dykes with a Zircon U-Pb age of 217−220 Ma on the eastern and northern flanks [75], and an angular unconformity between the Jurassic and Triassic (Figure 1).
The first extensional regime (non-LPS3a) belongs to NW−SE regional extension on the eastern margin of the Eurasian continent in the Cretaceous. This regional extension is characterized by a spread of NE-trending grabens or half-grabens in the central and eastern portions of the South China continent [76,77]. It is attributed to a variety of possible tectonic processes such as the collapse of the Early-Mesozoic marginal orogen [45], the high-angle subduction of the Pacific oceanic plate [9,40] and the roll-back of the subduction zone [46,47].
The second extensional regime (non-LPS3b) seems not to have a counterpart of regional extension. More about this will be discussed below.
Furthermore, in the ideal upper crust, where deformation is accommodated by the reactivation of pre-existing faults, the minimum differential stress gradient (MDSG) under the thrust regime is calculated in the below expression [78]: where σ 1 and σ 3 are the maximum and minimum principal stresses, z is the buried depth, R is the minimum σ 1 /σ 3 , ρ is the rock density, g is the gravitational acceleration, and λ is the pore-fluid factor. ρg is the hydrostatic pressure gradient at the depth (HPG  (Figure 12), are much smaller than this MDSG value but slightly smaller than the normal HPG. This indicates a difference in stress profile between the LPS state and the typical thrusting state, as it is.

Chronology of Stress Subsets
As previously described, field observations do not provide direct evidence on the relative ages of the two fold sets, longitudinal and transverse, that made the Longshan dome. They are however informative in understanding the deformation history of the dome in the following respects: (1) Secondary deformation structures pertaining to the longitudinal folds, such as bedding-parallel striations ( Figure 3d) and spaced cleavages (Figure 3e), have an older age than the ones most probably pertaining to the transverse folds.
(2) Cleavages associated with the longitudinal folds are much more common than the ones associated with the transverse folds ( Figure 2). This suggests a greater deformation temperature for the longitudinal folds.
(3) A combination of wedge-shaped calcite veins and stylolites (Figure 3f), diagnostic of the relatively high deformation temperature, has a generally older age than planar calcite veins. It formed under the local shortenings, roughly N-S and roughly E-W.
(4) The difference in vein sequence between the neighboring locations on the southern flank (Figure 3f,g) seems to characterize a complex deformation history of the dome.
These observations are helpful but insufficient to determine the relative and even absolute ages of the stress subsets. For example, the occurrence of parallel veining does not permit, as utilized [23], the use of the crosscutting relationships between the vein sets in timing the stress subsets favorable to form the sets. For this purpose, it is necessary to take into consideration more relevant geological observations and inferences as follows: 1.
In the east of Dongkou county there exists an angular unconformity where the folded platformal cover is overlain by Lower-Jurassic lacustrine sandstones (Figure 1b; [22]). The overlain rocks are much less strongly folded than the underlying rocks, where the longitudinal folds predominate. In this light the NW-SE compression should have made a more significant contribution to the deformation of the cover in the Middle-Late Triassic than in the Late Jurassic-Early Cretaceous. 2.
Exposed in Hunan province are numerous granitic intrusive bodies with zircon U-Pb ages of 210−225 Ma [81][82][83][84][85][86]. This widespread intrusion is presumably attributed to the mechanically thickening of the crust caused by intense compression at the time [87]. Meanwhile, Mesozoic granitic intrusion on the continent has a strong tendency to migrate towards the east with time [10]. The crust of the region would have been in a generally cooling state since the Middle-Late Triassic compression.

3.
As previously described, thick twins abundant in samples are indicative of a relatively high deformation temperature (HDT), in correspondence with slaty or spaced cleavages in mudstones and muddy limestones, associated principally with the longitudinal folds. An issue thus arises about the recognition of the HDT stress subsets.
In order to resolve this issue, the ratio of stress estimates for thick or thin twins to samples (RES) for a certain stress subset is introduced ( Figure 13). It is more likely to have a HDT stress subset when there is a larger RES value for thick twins. In Figure 11, most of the subsets are separated into two groups: HDT (LPS1, LPS2 and non-LPS1) and non-HDT (non-LPS2 and non-LPS3b). The former group has a larger RES value for thick twins than for thin twins, and the latter group possesses a larger RES value for thin twins. Between the groups is the non-LPS3a subset that has a larger RES value for thick twins.
Based on these considerations the stress subsets are rearranged temperature-decreasing order and further grouped into three stages in the Middle-Late Triassic (LPS1, LPS2, non-LPS1 and non-LPS3a), Jurassic-Early Cretaceous (non-LPS2) and extension in the Creta However, there are some uncertainties regarding this arrangement. older or younger than the LPS2 is yet unknown and requires addi The non-LPS3a is placed immediately after the non-LPS1, since th approximately parallel to the longitudinal fold axes, is considered a ation of strain accumulated by the latter intense deformation. Interes principal stress axis is nearly perpendicular to the strikes of the mo portant Au-Sb-bearing quartz veins in the basement core. These ve WNW-striking and have a 40 Ar- 39 Ar age of 161.1±1.2 Ma [88]. Figure 13. The ratio of stress estimates for thick and thin twins to samples fo set.

Tectonic evolution
According to the established sequence of the stress subsets (F area has the following scenarios of tectonic deformation during the M 1. In the Middle-Late Triassic, both the westward subduction of oceanic basin(s) and the collision between the South China and nents were responsible for the tectonic stress field in the central China continent (Figures 14a-b). The subduction process was su to deform the crust of the portion by means of folds and thru strongly implied by the HDT state of the LPS subsets ( Figure 1 have been heated and hence weakened before the onset of inten line of probable evidence on this heating is the widespread oc granitic intrusive bodies in the region. Given the common phenomenon of fold-axis parallel extensio [89][90][91][92][93], we would like to relate the non-LPS3a extension ( non-LPS1 compression (Figure 14c) that generated the NNE-tren however a lack of field observation about the relationship betw mechanisms about such a syn-folding extension include the diffe Figure 13. The ratio of stress estimates for thick and thin twins to samples for a certain stress subset.
Based on these considerations the stress subsets are rearranged in a younging and temperature-decreasing order and further grouped into three stages ( Figure 14): folding in the Middle-Late Triassic (LPS1, LPS2, non-LPS1 and non-LPS3a), refolding in the Late Jurassic-Early Cretaceous (non-LPS2) and extension in the Cretaceous (non-LPS3b). However, there are some uncertainties regarding this arrangement. Whether the LPS1 is older or younger than the LPS2 is yet unknown and requires additional observations. The non-LPS3a is placed immediately after the non-LPS1, since the former extension, approximately parallel to the longitudinal fold axes, is considered as a product of relaxation of strain accumulated by the latter intense deformation. Interestingly, its minimum principal stress axis is nearly perpendicular to the strikes of the most economically important Au-Sb-bearing quartz veins in the basement core. These veins are vertical and WNW-striking and have a 40 Ar-39 Ar age of 161.1 ± 1.2 Ma [88].
Geosciences 2021, 11, x FOR PEER REVIEW 22 of 26 type 2 cross fold, whose interference map tells the relative ages of the fold sets. What is more, the approximately N−S regional compression responsible for transverse folds happened in Late Jurassic-Early Cretaceous and should have a profound effect on the intensely deformed continent. This has been ignored in previous works.

Tectonic Evolution
According to the established sequence of the stress subsets (Figure 14), the study area has the following scenarios of tectonic deformation during the Mesozoic:

1.
In the Middle-Late Triassic, both the westward subduction of the western Pacific oceanic basin(s) and the collision between the South China and North China continents were responsible for the tectonic stress field in the central portion of the South China continent (Figure 14a-b). The subduction process was sufficiently enhanced to deform the crust of the portion by means of folds and thrusts ( Figure 14c). As strongly implied by the HDT state of the LPS subsets (Figure 13), the crust should have been heated and hence weakened before the onset of intense deformation. One line of probable evidence on this heating is the widespread occurrence of Triassic granitic intrusive bodies in the region. Given the common phenomenon of fold-axis parallel extension in orogenic belts [89][90][91][92][93], we would like to relate the non-LPS3a extension (Figure 14d) to the non-LPS1 compression (Figure 14c) that generated the NNE-trending folds. There is however a lack of field observation about the relationship between them. Possible mechanisms about such a syn-folding extension include the difference in rate of amplification along an individual fold [94], and flattening on the fold axial plane [94].

2.
In the Late Jurassic-Early Cretaceous, the pre-existing NE-or NNE-trending fold set was refolded under the N-S or NNE-SSW compression (Figure 14e). The regional refolds comprise the dome-and-basin structure and the bending of longitudinal folds, as previously described. They belong to Ramsay's [7] type 1 and type 2 superimposed folds, respectively. Both are common in nature for the interference of buckles. According to Ghosh et al.'s [95] experimental study, the type of fold interference depends on the tightness of pre-existing folds. The former type prefers relatively gentle preexisting folds and the latter type would produce relatively tight pre-existing folds. The non-HDT state of the non-LPS2 subset indicates the cooling of the study area caused by either uplift and erosion or the descent of geothermal flow. Under either condition, there is an increase in intensity of tectonic deformation from the LPS1 to non-LPS2. That is the reason why at this time the ancient fold-and-thrust belt on the northern margin was reactivated and propagated inland toward the south [96,97]. These deformation structures are the products of far-field compression that originated principally from the closure of the Mongolia-Okhotsk Ocean to the north [11,42]. Meanwhile, the continuous subduction of the western Pacific oceanic basin(s), for example the Izanagi basin, resulted in a more intense NW-SE compression that generated a remarkable fold belt on the eastern side of the Yangtze platform, about 400 km in width [98]. How this compression affected the pre-existing folds and faults in the region is poorly understood. 3.
In the Cretaceous after intense compression, NW-SE regional extension prevailed in the South China continent. It is characterized by grabens or half grabens bounded by NE-striking high-angle normal faults and infilled by red beds. In the uplift mountains this extension is well recorded by calcite e-twins in multi-deformed limestones (Figure 14f; [29]).
In the scenarios, NE−NNE-trending folds have an older age than E-trending folds. This is completely different from the generally accepted model for the formation of the dome [15][16][17][18][19][20]. The most important line of field evidence on this is the bending of longitudinal folds to the south of the study area (Figure 1b). It is typical of Ramsay's [7] type 2 superimposed folds, and its interference map tells the relative age of two fold sets: the curved fold set has an older age. Unfortunately this phenomenon has been ignored or, as above mentioned, mis-explained ever since.

Conclusions
The oval-shaped Longshan dome is located in the central portion of the South China continent. It is the product of mega-fold superimposition during the Mesozoic. It is generally accepted that NE−NNE-trending folds overprinted E−W-trending folds to form the dome. Unfortunately, this model has not been confirmed ever since, because the crossfold interference map does not tell the relative ages of the fold sets, and because there is no reliable field observation about the relationship regarding secondary deformation structures with respect to the sets.
In this paper we tried to resolve this issue by applying paleostress analysis of calcite e-twins to the dome. Ten limestone samples were collected from the Upper Paleozoic platformal cover on the flank of the dome, and measured for e-twins under the universal stage. E-twins in samples are thin (<1 µm) or thick (≥1 µm), and thick twins are diagnostic of relatively high temperature. Shan et al.'s [30] inversion method for polyphase twin data was improved and then applied to the measured datasets. The obtained stress estimates were classified into two layer-parallel shortening (LPS) subsets and three non-LPS subsets. These subsets comprise four tectonic regimes: NWW−SEE compression (LPS1 and non-LPS1), NNE−SSW compression (LPS2 and non-LPS2), NW−SE extension (non-LPS3a) and NNE−SSW extension (non-LPS3b). They are further arranged in a temperature-decreasing order to establish a complex deformation sequence of the study area.
There are three distinct stages of the sequence: folding in the Middle-Late Triassic (LPS1, LPS2, non-LPS1 and non-LPS3a), refolding in the Late Jurassic-Early Cretaceous (non-LPS2) and extension in the Cretaceous (non-LPS3b). In this manner E−W-trending folds overprinted NE−NNE-trending folds to form the dome. This is supported by the presence of E-W-trending transverse folds to the south. These folds are the bends of NE−NNE-trending longitudinal folds. Such a transverse fold is a typical Ramsay's [7] type 2 cross fold, whose interference map tells the relative ages of the fold sets. What is more, the approximately N−S regional compression responsible for transverse folds happened in Late Jurassic-Early Cretaceous and should have a profound effect on the intensely deformed continent. This has been ignored in previous works.