Models of the rapid post‐rift subsidence in the eastern Qiongdongnan Basin, South China Sea: implications for the development of the deep thermal anomaly

The Qiongdongnan Basin is one of the largest Cenozoic rifted basins on the northern passive margin of the South China Sea. It is well known that since the Late Miocene, approximately 10 Ma after the end of the syn‐rift phase, this basin has exhibited rapid thermal subsidence. However, detailed analysis reveals a two‐stage anomalous subsidence feature of the syn‐rift subsidence deficit and the well‐known rapid post‐rift subsidence after 10.5 Ma. Heat‐flow data show that heat flow in the central depression zone is 70–105 mW m−2, considerably higher than the heat flow (<70 mW m−2) on the northern shelf. In particular, there is a NE‐trending high heat‐flow zone of >85 mW m−2 in the eastern basin. We used a numerical model of coupled geothermal processes, lithosphere thinning and depositional processes to analyse the origin of the anomalous subsidence pattern. Numerical analysis of different cases shows that the stretching factor βs based on syn‐rift sequences is less than the observed crustal stretching factor βc, and if the lithosphere is thinned with βc during the syn‐rift phase (before 21 Ma), the present basement depth can be predicted fairly accurately. Further analysis does not support crustal thinning after 21 Ma, which indicates that the syn‐rift subsidence is in deficit compared with the predicted subsidence with the crustal stretching factor βc. The observed high heat flow in the central depression zone is caused by the heating of magmatic injection equivalently at approximately 3–5 Ma, which affected the eastern basin more than the western basin, and the Neogene magmatism might be fed by the deep thermal anomaly. Our results suggest that the causes of the syn‐rift subsidence deficit and rapid post‐rift subsidence might be related. The syn‐rift subsidence deficit might be caused by the dynamic support of the influx of warmer asthenosphere material and a small‐scale thermal upwelling beneath the study area, which might have been persisting for about 10 Ma during the early post‐rift phase, and the post‐rift rapid subsidence might be the result of losing the dynamic support with the decaying or moving away of the deep thermal source, and the rapid cooling of the asthenosphere. We concluded that the excess post‐rift subsidence occurs to compensate for the syn‐rift subsidence deficit, and the deep thermal anomaly might have affected the eastern Qiongdongnan Basin since the Late Oligocene.


INTRODUCTION
The widely used uniform lithosphere-stretching models (McKenzie, 1978;Jarvis & McKenzie, 1980) provide a classic quantitative framework for predicting the theoretical tectonic subsidence history of a rifted basin or margin. However, tectonic subsidence anomalies (significant deviations from the theoretical tectonic subsidence pattern) are observed in many rifted basins and passive margins (e.g. Royden & Keen, 1980;Buck, 2004;Davis & Kusznir, 2004;Ziegler & Cloetingh, 2004). Particularly, some of these rifted basins and passive margins show synrift tectonic subsidence deficit and extra post-rift thermal subsidence. In the distal domains of Angolan and Brazilian conjugated margins (Moulin et al., 2005;Contreras et al., 2010), the NW Australian margin (Karner & Driscoll, 1999) and the Iberia-Newfoundland margin (P eron-Pinvidic & Manatschal, 2008), shallow marine or even sub-aerial conditions developed during extreme crustal thinning (P eron-Pinvidic & Manatschal, 2008;Franke et al., 2014), and the present deep water environments suggest that the subsidence patterns in these areas were in significant subsidence deficit at the onset of seafloor spreading and rapid syn-seafloor spreading subsidence. The continental margins of the South China Sea also show anomalous subsidence patterns. For example, extra post-rift subsidence in the Late Cenozoic is revealed in the Yinggehai and Qiongdongnan Basins in the NW part of the South China Sea (e.g. Gong et al., 1997;Xie et al., 2006) and the Pattani and Malay Basins in the Gulf of Thailand (Morley & Westaway, 2006) (Fig. 1a). Recently, using the indication of shallow-water platform carbonates, Franke et al. (2014) and Ding et al. (2013) argued for a syn-rift tectonic subsidence deficit on the continental margins of the South China Sea.
Various syn-to post-rift tectonic mechanisms have been proposed to account for tectonic subsidence anomalies. These mechanisms include depth-dependent lithospheric stretching (Royden & Keen, 1980;Davis & Kusznir, 2004), simple shear model (Driscoll & Karner, 1998), phase transformation (Podladchikov et al., 1994), injection and underplate of mantle-derived melts , mantle serpentinization (Davis & Kusznir, 2004), lower crustal flow (Morley & Westaway, 2006), mantle plume (Nadin et al., 1997;Clift & Turner, 1998), influx of a warmer asthenosphere (Reston, 2009), dynamic topography (Xie et al., 2006), vigorous asthenospheric convection (Calv es et al., 2008), lateral flow of low-density depleted asthenosphere (Buck, 2004), variation in intraplate stress  and inelastic yielding caused by loading (Burov & Diament, 1995). Some of these mechanisms could generate a permanent subsidence anomaly by changing the crustal or mantle materials, whereas other mechanisms related to thermal anomalies or mechanics simply generate temporary subsidence anomalies. When these mechanisms generating temporary subsidence anomalies lose their supporting conditions, the tectonic subsidence would be reversed to compensate for the earlier subsidence excess or deficit (e.g. Nadin et al., 1997). Thus, the obtained crustal or  Fig. 1b show two seismic images of exposed igneous bodies. lithospheric stretching factor might be over-or underestimated based on just one phase's anomalous subsidence pattern (Ziegler & Cloetingh, 2004;Morley & Westaway, 2006). Therefore, to find out the most possible mechanisms accounting for the anomalous subsidence of a rifted basin and margin, it is necessary to investigate its entire subsidence history and episodic subsidence features.
The Qiongdongnan Basin is one of the major Cenozoic rift basins on the northern margin of the South China Sea (Fig. 1). Over more than 30 years of hydrocarbon exploration high-quality well and seismic reflection data have been acquired in the Qiongdongnan Basin. Based on these data, all of the studies on tectonic subsidence of the basin reveal a similar rapid post-rift subsidence phenomenon (e.g. Gong et al., 1997;Xie et al., 2006;Yuan et al., 2008;Li et al., 2012;Zhao et al., 2013). Notably, the rapid post-rift subsidence did not immediately follow the end of the syn-rift phase at 21 Ma; it began after a 10 my period of relative tectonic quiescence. Although the phenomenon is clearly observed and well accepted, its tectonic mechanism is still open to debate. Different tectonic mechanisms have been proposed to account for the rapid post-rift subsidence event, such as further crustal thinning due to a new rifting episode and polarity change in its western boundary fault (Yuan et al., 2008), lower crustal flow (Lei et al., 2013;Zhao et al., 2013), dynamic topography (Xie et al., 2006), thermal anomalies and magmatic injection (Li et al., 2012). All of these proposed mechanisms focus only on the rapid post-rift subsidence episode since the Late Miocene. However, our detailed analyses reveal a two-stage anomalous subsidence feature of the syn-rift subsidence deficit and the well-accepted rapid post-rift subsidence. Here, we first present the entire tectonic subsidence and thermal features, and then, we use a coupled rifting lithosphere-scale and basin-scale numerical model to analyse the origin of a recognized high heat-flow zone and the two-stage anomalous subsidence pattern. Our analyses suggest that the two-stage anomalous subsidence might be related to the development of the deep thermal anomaly, which might have affected the eastern Qiongdongnan Basin since the Late Oligocene, and the rapid post-rift subsidence occurs to compensate for the syn-rift subsidence deficit.

GEOLOGICAL BACKGROUND
The northern continental margin of the South China Sea is a rifted continental margin, which has experienced multi-episodic rifting since the Late Cretaceous (Ru & Pigott, 1986). It has evolved into a small-scale passive continental margin with seafloor spreading commencing at 30 Ma and terminating at 16 Ma (Taylor & Hayes, 1983;Briais et al., 1993;Li et al., 2014). Several Cenozoic rift basins, such as the Beibuwan, Pearl River Mouth, Qiongdongnan and Yinggehai Basins, developed on the northern margin of the South China Sea (Fig. 1). The first three basins, formed under a NW-SE extensional stress field, are NE-trending, showing a southward younging of the initial rifting from the Late Cretaceous in the north to 45 Ma in the south (Sun et al., 2009), whereas the NW-trending Yinggehai Basin, located along the offshore extension of the Red River Fault Zone (Fig. 1), is thought to be a transform-extensional basin related to the activities along the Red River Fault Zone (Sun et al., 2003;Clift & Sun, 2006;Zhu et al., 2009).
The Qiongdongnan Basin is located between the Hainan Island and Xisha Block (Fig. 1). To the west, it is separated from the Yinggehai Basin by the No. 1 Fault, which is an offshore extension of one main branch of the Red River Fault Zone (Clift & Sun, 2006;Zhu et al., 2009). To the east, the Qiongdongnan Basin is connected through the Xisha Trough to the oceanic NW sub-basin of the South China Sea, whose seafloor spreading began at ca. 30 Ma, and ceased at 28.5 Ma (Briais et al.,1993;Barckhausen et al., 2014). From north to south, the Qiongdongnan Basin is divided into a northern depression zone, northern rise zone, central depression zone and southern rise zone. From west to east, the central depression zone consists of the Yanan, Ledong, Lingshui, Beijiao, Songnan, Baodao and Changchang Sub-basins. The syn-rift phase began in the Eocene (45 Ma or earlier) and ended 21 Ma (e.g. Xie et al., 2006). A clear unconformity, developed 21 Ma, divides the Cenozoic into syn-rift and post-rift sequences. After 21 Ma, the majority of the faults became inert (Gong et al., 1997;Xie et al., 2008), and water depth gradually increased due to the post-rift thermal subsidence. During the late Early Miocene through the Middle Miocene, extensive shallow-water carbonate platforms and reefs began to develop on the uplifted fault blocks of the south rise zone and the Songnan Low Rise (e.g. Chen et al., 2011;Huang et al., 2012). At the end of the Middle Miocene, a clear unconformity was developed in the shallow northwestern basin, and in the eastern central depression zone, there was an abrupt depositional environment transformation from an onshore, neritic environment to a shelf-slope, bathyalabyssal environment (Xie et al., 2008;Su et al., 2011), suggesting that the water depth has increased rapidly since 10.5 Ma (Xie et al., 2008;Su et al., 2011). The rapidly increased water depth also resulted in the drowning of the shallow-water carbonates and reefs. The thickest Cenozoic strata appears in the Ledong Sub-basin and Lingshui Sub-basin, where seismic investigations indicate that the largest Cenozoic thickness is larger than 10 km, due to the greatly increased Red River-derived sediment since the Late Miocene (Xie et al., 2008;Zhao et al., 2013). In the eastern basin, where the water depth is >2000 m, the Cenozoic thickness is only 6 km because most of the Red River-derived sediment has been trapped in the western basin.
Since the cessation of the seafloor spreading at 16 Ma (post-spreading phase), basaltic magmatism has become considerably more widespread and active, particularly in South China (e.g. Leiqiong Depression of the Beibuwan Basin), Indochina and the South China Sea (Fig. 1). After analysing the Cenozoic volcanism age data, Xu et al. (2012) suggested that the greatest amount of igneous activity appears to have occurred from the Pliocene to recent times in the South China Sea and its adjacent regions. Regional and global tomographic studies suggest that there is a SE dipping low-velocity zone beneath Hainan region feeding the widespread late Cenozoic magmatism (e.g. Huang & Zhao, 2006;Lei et al., 2009). Based on seismic reflection profiles, magnetic data, geomorphology and drill samples, magmatic injections during the post-spreading phase have also been identified in the Xisha Trough and Qiongdongnan Basin (Liu & Wu, 2006;Lu et al., 2011;Yang et al., 2011; (Fig. 1). For example, a tholeiitic basalt sill with a thickness of 115 m was found in the Pliocene of Drill Y32-1-1, offshore of SW Hainan (Li et al., 1998). Igneous injections could be identified on seismic profiles in the eastern Qiongdongnan Basin (e.g. Ji et al., 2014), some of which are exposed at the seafloor ( Fig. 1). Two examples of these igneous bodies are shown in Fig. 1b. Because they are in deep water, the exact ages of these igneous injections are unclear, but the magmatism in the eastern basin typically occurred in the Miocene or more recently.

TECTONIC SUBSIDENCE AND GEOTHERMAL FEATURES Tectonic subsidence
There are numerous studies on the tectonic subsidence history of the Qiongdongnan Basin (e.g. Gong et al., 1997;Xie et al., 2006;Yuan et al., 2008;Li et al., 2012;Zhao et al., 2013). Although they employ different methods and are based on somewhat different data due to uncertainty in the paleobathymetry and stratigraphic data, all of these studies concluded that there was a rapid subsidence event after 10.5 Ma, approximately 10 Ma after the end of the syn-rift phase, suggesting that the Late Miocene-Quaternary rapid tectonic subsidence phenomenon is robust. To further analyse the tectonic subsidence pattern, we selected 41 representative points as synthetic drills that cross the main deep-water sub-basins of the Qiongdongnan Basin (Fig. 2), and then with the stratigraphic data derived from the commercial seismic reflection data, we calculated three types of tectonic subsidence histories as follows: 1 Calculate the observed tectonic subsidence history (TSH o , solid curves in Fig. 2) with the traditional backstripping and decompaction method from stratigraphic data (Sclater & Christie, 1980;Shi et al., 2005); 2 Calculate the predicted tectonic subsidence history (TSH p1 , dashed curves in Fig. 2). We first calculated the strain rate history by fitting the observed tectonic subsidence history (TSH o ) with the 1-D strain rate inversion method of White (1994) and Li et al. (2012), and then calculated the stretching factor (b s ) and the predicted tectonic subsidence history (TSH p1 ) assuming that the lithosphere stretches with the obtained strain rate history during the syn-rift phase (45-21 Ma) and then experiences thermal relaxation after 21 Ma with the forward method of White (1994) and Li et al. (2012); 3 Calculate the predicted tectonic subsidence history (TSH p2 , dotted curves in Fig. 2) assuming that the lithosphere stretches with a constant strain rate In b c / Δt during the syn-rift phase (45-21 Ma) with the same forward method as (2). b c is the observed crustal stretching factor, obtained by dividing the initial crust thickness by the present crust thickness, and Dt is the duration of the syn-rift phase.
In order to derive observed tectonic subsidence, paleobathymetry, the effect of compaction, sedimentary loading, and sea level variations during deposition must be treated with great precautions (Sclater & Christie, 1980). We determined the palaeobathymetry and lithology of each synthetic drill mainly based on the sedimentary facies and palaeogeography environment of Zhu et al. (2007), commercial drill data, the developments of the northern shelf-slope system (Xie et al., 2008) and the late Miocene reef-banks, the sedimentation facies variation at ca. 10.5 Ma (Su et al., 2011), time-dependent distance from the shelf/slope break and present water depth. Since the sedimentary facies during the syn-rift phase are lacustrine in the Eocene through the Early Oligocene, marsh and coastal plain in the late Early Oligocene, and littoral and neritic facies in the Late Oligocene, the uncertainty in estimating the syn-rift palaeobathymetry is generally no more than 100 m. However, after 21 Ma, with the depositional environment transformation from dominantly neritic environment during the Early Miocene through the Middle Miocene, to a shelf-slope, bathyal-abyssal environment (Xie et al., 2008;Su et al., 2011), the uncertainty in palaeobathymetry estimation might be up to 200 m at 15.5 Ma, 300 m at 10.5 Ma, and 500 m at 5.5 Ma for those synthetic drills lying in the central depression zone. These palaeobathymetry uncertainties, which are less than the amount of the anomalous subsidence ( Fig. 2), do not significantly affect the following subsidence features.
For decompaction, the exponential law /(z) = / 0 exp (Àcz) was used to calculate the porosity / at depth z, where / 0 is the depositional porosity and c is the compaction factor depending on the lithology (Sclater & Christie, 1980). We used common physical parameters (grain density, initial porosity and compaction factor) of different sediment types (Sclater & Christie, 1980) (Table 1). We found that the exponential law /(z) = / 0 exp(Àcz) could predict the neutron porosity-depth data of eight drills (including five deep-water drills). The third order sea level reconstruction of Haq et al. (1987) was used to estimate the global eustatic sea level.
The present Moho depth map (Fig. 2) presented by Zhang et al. (2008) is used to calculate the observed crustal stretching factor b c . The Moho depth map was derived mainly based on gravity data. Table 2 lists the gravity-derived Moho depths and the Moho depths derived from expanding spread profile data along Profile EF (Nissen et al., 1995b), wide-angle reflection seismic data along Profiles GH (Qiu et al., 2001) and IJ (Huang et al., 2011), and earthquake seismology (Huang et al., 2011) (Fig. 2), showing that the difference of the Moho depths is generally lower than 3 km. Gravity-derived results have a smoothing trend by increasing Moho depths in the centre depression zone and the Xisha Trough where the basement crust has been greatly extended and thinned, and decreasing the Moho depths in the shelf and the Xisha Block. The initial crust thickness is assumed to be 32 km, which is the thickness of the crust close to the coast (Nissen et al., 1995a;Clift et al., 2002). Thus, the calculated b c might be slightly smaller in the centre depression zone and the Xisha Trough. Figure 2 shows the tectonic subsidence histories of TSH o (solid curves), TSH p1 (dashed curves) and TSH p2  . The profiles of EF, GH and IJ are the deep crustal structure sections of Nissen et al.(1995b), Qiu et al.(2001) and Huang et al.(2011). The map also shows other representative synthetic drills (open circles) whose results are not described. The solid, dashed and dotted black lines represent observed tectonic subsidence (TSHo) and the predicted tectonic subsidence with stretching factors of b s (TSHp1) and b c (TSHp2) respectively. The basement subsidence history of Drill PA is shown in a red curve. For legibility, the observed subsidence just present the error bars of palaeobathymetry which are larger than AE100 m.
(dotted curves) of 10 synthetic drills and one real drill (PA), which suggest rapid post-rift subsidence occurred approximately 10 Ma after the end of the syn-rift phase, particularly in the central depression zone of the Qiongdongnan Basin. Since the calculation of the observed tectonic subsidence involves various parameters, it is very hard to estimate the error bars of the observed tectonic subsidence. Figure 2 only shows the most important error due to the uncertainty in palaeobathymetry estimation. Because TSH p1 is calculated assuming the syn-rift strain rate history obtained by fitting TSH o , TSH p1 is nearly the same as TSH o during the syn-rift phase. However, after entering the post-rift phase, the TSH o curve begins to deviate from the predicted thermal subsidence curve of TSH p1 , and the deviation becomes considerably greater, particularly after 10.5 Ma. The final anomalous subsidence, which is equal to the difference between TSH o and TSH p1 , could reach 1000-2500 m in the central depression, with a decreasing trend from the Ledong Sub-basin (P5) to the Songnan Sub-basin (P6) and an increasing trend from the Songnan Sub-basin (P6) to the Changchang Sub-basin (P8) (Fig. 2). Along the NW-trending profile across P1-P4, the total tectonic subsidence and post-rift anomalous subsidence increase from the shelf to the central depression zone. Figure 2 demonstrates the differential subsidence between the northern depression zone, the central depression zone and the Xisha Block. Well Xiyong-1 (PA in Fig. 2), lying in the Xisha Block, drilled the basement of Precambrian granite gneiss at the depth of 1251 m, over which are all shallow-water platform carbonates and reefs developed since the Miocene (Wei et al., 2005;Wu et al., 2014b). The observed subsidence history of the synthetic drill PB shows a rapid subsidence feature since the Late Miocene with an average thermal subsidence rate of 250 m Ma À1 , which is much larger than the rate at PA (43-63 m Ma À1 ), and another synthetic drill PC (50 m Ma À1 ) in the Songdong Sub-basin. The difference in the thermal subsidence rate suggests that after 10.5 Ma, there is significant differential subsidence between the central and the northern depression zones and the Xisha Block, and the rapid post-rift subsidence is mainly confined in the central depression zone, where the curst has been highly extended and thinned.
However, a comparison of TSH o with TSH p2 in Fig. 2 indicates that the observed tectonic subsidence TSH o is considerably less than the predicted TSH p2 at the end of the syn-rift phase. The difference would be slightly larger if the gravity-derived Moho depth overestimates the real depth in the central depression zone. Another impressive feature is that the difference between the predicted today's TSH p2 and the observed today's TSH o (0 Ma) is lower than 0.5 km at most synthetic drills except P5, P6, P2 and PA. The tectonic subsidence difference at these three synthetic drills is about 1000 m. The larger difference at P5, P6 and P2, might be due to the errors in estimating crust thickness and b c , particularly in the western basin, where the sediment cover is thicker than 10 km (Xie et al., 2008;Zhao et al., 2013). The larger difference at PA might suggest that we could not ignore the lithospheric strength of the Xisha Block, and its real tectonic subsidence history should range from the calculated TSH o and the total basement subsidence (Fig. 2). Thus, if the lithosphere was extended with stretching factor b c during the syn-rift phase, any explanation must accommodate both the observed syn-rift tectonic subsidence deficit and rapid post-rift subsidence that occurs 10 Ma after the end of the syn-rift phase.

Geothermal distribution features
Previous thermal regime analyses of the Qiongdongnan Basin and Xisha Trough were mainly based on limited heat-flow data from commercial drills and seafloor probe measurements (Nissen et al., 1995a;He et al., 2001;Shi et al., 2003;Yuan et al., 2009;Shan et al., 2011).
Recently, a set of valuable drill-derived and seafloor probe heat-flow data has been obtained in the Qiongdongnan Basin with increasing deep-water hydrocarbon exploration and drilling (Xu et al., 2006;Wang et al., 2014;Shi et al., 2015), which greatly improve the understanding of the thermal regime in the deep-water area. To further reveal the thermal regime of the eastern Qiongdongnan Basin, in September 2011, we carried out a new heat-flow survey onboard the R/V Shiyan 3 at nine sites (SY3-1109) along the 18N survey line in the eastern Qiongdongnan Basin using a 4-m-long Ewing-type heat-flow probe. We calculated environment temperatures and geothermal gradients following the method of Pfender & Villinger (2002), and we measured the thermal conductivity of the collected sediment samples with a TK04 thermal conductivity meter from TeKa (Berlin, Germany). Figure 3 presents the corresponding environmental temperature at different tilt-corrected depths (referring to the depth of the uppermost effective temperature logger), and their tilt-corrected temperature gradient lines obtained by the linear least square method. Figure 3 shows that temperature gradients in the sediment are mostly linear, suggesting that hydrothermal activity has little effect on these temperature gradients. We now first present the six heat-flow values obtained by multiplying the geothermal gradient by the thermal conductivity (Table 3). Thus, a data set of 154 heat flows (Fig. 4a), including 44 drill-derived heat flows and 110 seafloor heat flows from the deep-water area, could be used to constrain the heat-flow distribution.
Because sediment thermal blanketing and heat generation have different influence on the drill-derived heat flow    (average heat flow in the upper sediment column of several kilometres) and the seafloor probe heat flow (average heat flow in the uppermost sediment column of several metres), particularly in a region with recent high sedimentation rate (Rupke et al., 2008), generally they are not equal to each other. However, based on the results of Wang et al. (2014), who studied the thermal histories of the Qiongdongnan Basin, the largest difference ( Fig. 4a) between sediment heat flow (averaged from basin basement to seafloor) and seafloor heat flow is typically <8 mW m À2 . Because a commercial drill is typically 3000-4000 m in depth, the difference between the drillderived heat flow and seafloor heat flow is typically <3 mW m À2 in the Qiongdongnan Basin (Shi et al., 2015). This suggests that we can use all of the heat flows together to analyse the heat-flow distribution features of the Qiongdongnan Basin. Figure 4b shows that the study area could be divided into three regions according to the heat-flow values. Heat flow on the northern shelf ranges from 50 to 70 mW m À2 (green area), the range of the central area is 70-85 mW m À2 , and in the eastern basin, there is a NE-trending high heat-flow zone of 85-105 mW m À2 . The eastward extension of this zone might join the high heat-flow zone in the lower slope of the northern margin (Shi et al., 2003).

NUMERICAL EXPERIMENTS Numerical model
We use a structural and thermal coupled kinematic numerical model modified from the numerical model of Shi et al. (2008) to explain the area of anomalous tectonic subsidence and high heat flow, considering the close relationship between tectonic subsidence and the thermal process. Shi et al. (2008) used a model in which the thermal evolution of basin-and lithosphere-scale processes were decoupled, and the temperature at the upper boundary of the lithospheric thermal field was set to be equal to the bottom temperature of the sediment column, which was approximated by assuming an appropriate constant geothermal gradient. This numerical model could be resolved simultaneously for basin-scale processes (e.g. sedimentation, compaction, stretching, heat transfer and maturation) and lithosphere-scale processes (e.g. rifting, flexure, heat transfer), similar to the forward model of Rupke et al. (2008). It allows for multiple rifting events of finite duration (Jarvis & McKenzie, 1980) controlled by a given strain rate history (White & Bellingham, 2002). To conserve mass, during the rifting phase, the lithosphere and sediment are both stretched by assuming a sediment column coupled with basement deformation, which is similar to the models of Theissen & R€ upke (2010) and Wangen et al. (2008). This model consists of a forward algorithm and an inverse algorithm (Fig. 5). The forward algorithm has four main subroutines. First, the velocity field is calculated with a given horizontal strain rate distribution to determine the original model length and subsequent lithospheric deformation. Second, the velocity field is used to reconstruct sedimentary history with input observed strata data. Third, the thermal evolution of the lithosphere and basin is calculated with the new model geometry, and the deposited sediment is incorporated into the thermal solver to include the effects of sediment blanketing and radioactive heat generation (Lucazeau & Le Douaran, 1985;Theissen & R€ upke, 2010). Fourth, the loads over the depth of compensation are calculated with the density structure constrained by the thermal structure and used to determine the new model geometry and basement depth through the flexural equation. The inverse algorithm is similar to the one used by Rupke et al. (2008) to obtain the average strain rate of each rifting episode by minimizing the difference between the expected and predicted basement depth at key times with the given palaeobathymetry.

Velocity field
In the model, x is the horizontal coordinate and z is the vertical coordinate, which increases from bottom to top (Fig. 6). During stretching, the left-hand boundary is fixed, and the right-hand boundary moves away with a horizontal velocity u(x,t). At the same time, the upper surface and base of the lithosphere move towards the necking plane, and the asthenospheric mantle flows upwards across the lower plane z = 0 km. The horizontal velocity is determined from the given horizontal strain rate G(x,t) (White & Bellingham, 2002): The vertical velocity is given as defined previously   Z neck is the necking depth ( Fig. 6) . In the Lagrangian mesh, the node coordinates change with the deformation. The stretching factor is defined as follows: Thus, given a b(x,t) of a rifting episode, the strain rate G(x,t) can be calculated by assuming that G(x,t) is constant during the episode, or with an inverse algorithm, the b(x,t) and strain rate G(x,t) of each rifting episode can be determined by minimizing the difference between the observed/expected and predicted basement depth at the end of each rifting episode.
With the velocity field, we can determine the original model length of the studied profile by automatic trial and error and rebuild its lithospheric deformation history.

Sedimentation
Given a strain rate G(x,t) and the observed strata data, the sediment history can be determined as follows. First, assuming that the porosity / varies with depth z and obeys the compaction law /(z) = / 0 exp(Àcz), where / 0 is the depositional porosity and c is the compaction factor depending on the lithology (Sclater & Christie, 1980), the matrix thickness Gth i for layer i from depth  . Triangles are seafloor heat-flow stations, and solid triangles represent new observations presented in this work. The black thick isolines with values in (a) show the difference (mW m À2 ) between the calculated drill-derived heat flow (average heat flow from basement to seafloor) and the calculated seafloor heat flow. The absolute value of the difference in the region outside the isolines of 2 and À2 mW m À2 is <2 mW m À2 . The solid black curve EF is the heat-flow profile of Nissen et al. (1995a,b). EM, Igneous body exposed at the seafloor. z 1 to z 2 can be obtained by calculating R z 2 z 1 ð1 À /ðzÞÞdz using Gth i as the initial matrix thickness Gth i0 . Second, assuming that the deposition rate is constant during t 1 through t 2 , the original deposited matrix thickness in a time step dt is Gth i0 *dt/(t 2 À t 1 ). Third, to conserve mass, the matrix thickness, once deposited, would be coupled with the basement and subsequently thinned based on the given G(x,t); thus, the matrix thickness evolving history Gth(x,t) can be obtained. Fourth, the detailed sedimentary history can be calculated by expanding the compaction law /(z) = / 0 exp(Àcz). Fifth, the initial matrix thickness Gth i0 of layer i is adjusted, and the process begins again at the second step until the difference between the calculated and input observed strata data is sufficiently small. If the given strain rate G(x,t) changes, the sedimentary history would be needed to be recalculated.

Temperature structure
The thermal field is derived from a 2D unsteady thermal diffusion equation: where q is the density, c is the specific heat, T is the temperature, t is time, k is the thermal conductivity, Q is the internal heat sources and z is the depth. The boundary conditions are shown in Fig. 6b. The upper boundary condition T(x, z 0 (x,t)) is set to be the seafloor temperature, which is related with water depth (Shi et al., 2015). If the temperature is higher than the local annual mean temperature 22°C or there is no water, T(x, z 0 (x,t)) is set to be 22°C. z 0 (x,t) represents the upper surface of the model. The lower boundary temperature T(x,0) can vary with time, and here, it is set to be T 1 = 1390°C according to the 3D S-velocity model of the South China Sea (Tang & Zheng, 2013) and the method of Priestley & McKenzie (2006). The initial thermal field is assumed to be steady with T(x, a) = 22°C. The code for solving Eqn (4) is generated with a commercial finite element package FEPG 3.0 (www.fegensoft.com/english/index.htm). The thermophysical parameters in Eqn (4) for the Cenozoic sediment are functions of porosity, temperature and lithology. The sedimentary thermal conductivity k is the geometric average between pore fluid and matrix conductivities, which are related to temperature according to the method of Makhous & Galushkin (2005). The other sediment thermophysical parameters (such as density q, and volumetric heat capacity qc) are the weighted arithmetic averages of the density or volumetric heat capacity of the matrix and that of the pore fluid. These parameter values are listed in Table 1.
Thermal conductivity for the crust and mantle can be expressed as the sum of the lattice conductivity k l , and radiative conductivity k r , which are temperature dependent. For the upper crust, lower crust and lithosphere mantle, k l can be defined, respectively, as follows (Buntebarth, 1984): k À1 l ¼ 0:4 þ 0:33 Â 10 À3 Tð CÞ k À1 l ¼ 0:41 þ 0:29 Â 10 À3 Tð CÞ k À1 l ¼ 0:21 þ 0:5 Â 10 À3 Tð CÞ For the upper crust, the original value of the righthand side of the equation is 0.33. We now adjust it to be 0.4 because we found that the original value 0.33 would make thermal conductivity overly large in the basement rise region with thin sediment. As a result, the heat-flow refraction effect would produce considerably higher heat flow on the shelf than the observed value. For the crust, k r can be ignored, but for the upper mantle, where the temperature is >230°C, the radiative conductivity k r for upper mantle can be defined as follows (Buntebarth, 1984): k r ¼ À0:52 þ 2:3 Â 10 À3 Tð CÞ ð 6Þ

Loading and Flexure
Based on Eqns (2)-(3), the upper surface initial depression S due to necking can be defined as follows: where a is the lithospheric thickness (Fig. 6). This initial basement subsidence is progressively modified by the flexure due to vertical loading and lateral force. The flexural deflection w(x) is computed from the plate equilibrium equation: Here, x is the horizontal coordinate, w is the vertical deflection, P is the fibre force, g is the acceleration due to gravity, D is the flexural rigidity, where D(x) = ETe(x) 3 / 12(1 À v 2 ), q sw is the normal load, including sediment and seawater load above the actual basement, q cm is the load due to thinned lithosphere, q r is the load due to the asthenosphere above the compensation depth and q 0 is the weight of an isostatically compensated column referred to the depth of compensation for the original undeformed area. The lithosphere is adjusted to account for thermal expansion based on the linear relation q(x,z, t) = q 0 (1 À aT), where a is the lithospheric thermal expansion coefficient. In this study, we ignore the horizontal force P(x), and we assume that E = 8 9 10 10 Pa and m = 0.25 (Burov & Diament, 1995). The basement subsidence St(x,t) is the sum of the deflection w(x,t) and the initial subsidence S(x,t) due to necking.
The above subroutines could be assembled into different codes for the following cases to find the potential origins of the anomalous subsidence and thermal feature in the central depression zone of the basin. To eliminate the possible effect of the western boundary fault (the No. 1 fault) on subsidence and thermal anomalies (e.g. Xie et al., 2006;Yuan et al., 2008;Hu et al., 2013), thus simplifying the origin analysis of the anomalous subsidence and thermal feature, profiles AB and CD, which are relatively far from the No. 1 fault, were selected for numerical simulations (Figs 2 and 7). Furthermore, there are some observed heat flows in close proximity to these two profiles that are highly valuable for constraining the calculated results (Figs 2 and 3).
The values of the model parameters are listed in Table 1. The 3D S-velocity tomography model of Tang & Zheng (2013) suggests that the seismic lithosphere thickness of the northern coastal area was 90 km, and we would predict a thermal lithosphere thickness of 100 km according to the temperature-S-velocity relationship of Priestley & McKenzie (2006). Based on the thermal regime of the northern margin of the South China Sea, Shi et al. (2000) Fig. 6. Illustrations of the numerical model. (a) The model is divided into four layers: Layer 1 (magnified vertically) contains asthenosphere material whose initial thickness is 1 km, and Layers 2-4 represent the lithospheric mantle, lower crust and upper crust (including basin sediment) respectively. The horizontal u(x,t) and vertical v(x,z,t) velocities for Layers 2-4 are defined in the text. The lower surface of Layer 1 is fixed vertically when its upper surface moves with the lithosphere. (b) Thermal model with boundary conditions. The elements used in this model are quadrangular, and their original widths are variable. The final width after stretching is <5 km, and the original heights are 0.6, 1.6, 6.8 and 0.1 km for upper crust, lower crust, mantle lithosphere and Layer 1 respectively. In the model, a = the original lithosphere thickness, T is temperature, t = time, Z neck = necking depth.
suggested that the thermal lithospheric thickness in the coastal area was 90-100 km. Therefore, we set the lithosphere thickness of our model (Fig. 6) to 100 km. We have tried different thicknesses for the lithosphere: a thicker lithosphere (e.g. 125 km) yields a slightly shallower basin basement than a thinner lithosphere, but it does not significantly affect the following results. The present crustal thickness, palaeobathymetry, lithology and the physical parameters of different sediment types for calculating sediment porosity are determined in the same way as the above 1D tectonic subsidence analyses. The density and heat generation of the upper and lower crusts are determined based on the gravity modelling and geothermal study along a crustal structure profile across the Xisha trough (Qiu et al., 2001;Shi et al., 2002). The sediment heat generation is given based on Makhous & Galushkin (2005). The calculated average sediment heat generation is close to the average sediment heat generation, 1.34 lW m À3 (Shi et al., 2015), which was determined based on gamma logging data. Previous studies (e.g. Clift et al., 2002;Shi et al., 2005;Wu et al., 2014a) suggested that in the hyper extension slope area of the northern margin of the South China Sea, the lithosphere was sufficiently weak that it was nearly in local isostasy. Considering that the basement crust in the central depression zone has been greatly thinned and that the heat flow is extremely high (70-105 mW m À2 ), local isostasy was assumed in the following cases. The results of the following Cases 2, 4 and 5 could predict satisfactorily the present basement depth, suggesting that the assumption of local isostasy is reasonable.

Experimental results
Before analysing the following cases, we first calculated the stretching factors b s1 (45-36 Ma), b s2 (36-30 Ma) and b s3 (30-21 Ma) from the syn-rift sequences by minimizing the referred ('observed') basement depth and predicted ('calculated') basement depth at the end of each rifting episode (i.e. 36, 30 and 21 Ma). Figure 8 shows the variations of b s1 , b s2 , b s3 , their product b s (=b s1 9 b s2 9 b s3 ), and the observed crustal stretching factor b c for profiles AB and CD (Fig. 2). Although b s and b c have similar variations along these profiles, b s is considerably smaller than b c , which indicates that b s is not sufficiently large to generate the present thinned crust thickness.
Because there is uncertainty in the palaeobathymetry, the comparison between the predicted basement depth (total subsidence) and referred basement depth is performed at certain specified times (21, 15.5, 10.5 and 0 Ma) when palaeobathymetry can be constrained reasonably. The referred ('observed') basement depth at time t is the sum of the water depth and Cenozoic sediment thickness at that time. However, similar to what is noted in the subsection of Sedimentation, to maintain material conservation, our model assumes that the previously deposited sediment would be stretched and thinned during a new rifting. Thus, sedimentary history is related to the strain rate G(x,t) (rifting history) in our 2D model. This approach is different from the traditional 1D backstripping technique (Sclater & Christie, 1980) (Fig. 9). The latter does not consider the thinning of older strata during later rifting. If the older strata are extremely thick and if the thickness of the newly deposited sediment from the new rifting combined with the water depth is small, the predicted basement would rise rather than subside (Figs 9 and 10c). Therefore, in the following cases, the referred basement depth curves are different if their rifting histories are different. Though the rifting histories of Cases 1 and 2 are different, the referred basement depths after 21 Ma are the same for Cases 1 and 2 because there is no rifting event after 21 Ma.
Case 1: Assuming the lithosphere was thinned uniformly by b s1 , b s2 and b s3 sequentially during the syn-rift phase, followed by normal lithospheric thermal relaxation and contraction.
The results show that the predicted basement at the end of the syn-rift phase (21 Ma) corresponds well with the referred basement (which is used to derive b s ), whereas the predicted present basement is considerably shallower than the observed basement (Figs 10a, d, 11b, and 12b). The predicted heat flow is lower than the observed one, particularly in the central depression (Figs 11a and 12a).
Case 2: Assuming the lithosphere was thinned uniformly by b c during the syn-rift phase, followed by normal thermal relaxation.   and Moho depth (Fig. 2). See Fig. 2 for their locations. The black solid circles labelled as P1, P2 and P3 are the selected points for comparison.
The predicted basement at 21 Ma is considerably deeper than the referred basement, but the predicted present basement is close to the observed one (Figs 10b, e, 11b  and 12b). The minor difference (<700 m) between the observed and predicted present basement might be caused by the error of gravity-derived Moho depth or the observed basement depth converted from the time profile. Although the predicted heat flow is somewhat higher than in Case 1, it is still lower than the observed heat flow, particularly in the central depression (Figs 11a and 12a).
The two cases above indicate that the stretching factor b c might be more suitable than b s to predict the observed present basement depth and generate the present crust thickness, which is similar to that suggested in the above 1D subsidence analyses.
Case 3: Assuming that there are four rifting episodes: the first three occurred as Case 1, and the last episode began at 5 Ma with a stretch factor b Figures 10c, f, 11b and 12b show that the difference between the observed present basement and the predicted present basement is <900 m, suggesting that the present basement could be reasonably predicted if assuming a new rifting event. However, the predicted basement subsidence curves are quite different from the referred curves after 21 Ma (Fig. 10c, f), which indicates that a new rifting might not be feasible. Figure 11a shows that although the predicted heat flow is slightly higher than that of Case 2, it is still considerably lower than the observed heat flow in the central depression. Furthermore, the clear geological fact that faulting has been weak since 21 Ma also does not support a new rifting episode.
In any case, the employed model could be reasonable only when the predicted basement subsidence path is in accordance with the referred one and when the observed heat flow can be predicted. The results show that the above three cases could not reproduce simultaneously the referred basement depth (Fig. 10), and the observed heat flow (Figs 11a and 12a). However, these analyses might suggest that the lithosphere was stretched and thinned during the syn-rift phase (before 21 Ma) with a stretching factor b c rather than b s , and the observed heat flow might result from a thermal event, such as magmatism, rather than further crustal thinning due to a new rifting.
Case 4: To account for the observed high heat flow, assuming that the lithosphere was thinned uniformly by b c is the observed crustal stretching factor calculated directly by dividing the prerift crust thickness by the present crust thickness. The present crust thickness is obtained by subtracting the Moho depth from the basement depth ( Fig. 7) along the two profiles. Fig. 9. Basement depth variation comparison between the traditional method (Sclater & Christie, 1980) and the new method during a rifting episode. The new method, which considers the thinning and compaction of the older strata (black column) after a rifting episode, would produce a deeper basement if the thickness of the new deposited sediment (white column) during the rifting episode is larger than the thickness variation due to thinning and compaction of the older strata (a), and it would produce a shallower basement if the new sediment is not thick enough (b). However, the traditional method, which only considers compaction, would produce a deeper basement regardless of the thickness of the new sediment (c and d).
b c during the syn-rift phase, followed by a normal thermal relaxation, and in the Pliocene, there were equivalent heating events due to magmatic injection into the crust (Fig. 13). When calculating the heat effect of a dike, Bialas et al. (2010) assumed that the newly accreted material entered at 1200°C and added the latent heat of dike solidification to the system. However, our study area might be affected by multiple magmatisms, and it is difficult to determine the size and injection time of each dike. Therefore, to fit the observed heat flow, our model simplifies the magmasim as an instantaneous equivalent injection of 1000°C material, but it affects a wider region by replacing lighter crustal material with basaltic material of 2950 kg m À3 at 0°C (Turcotte & Schubert, 1982) and increasing the lower temperature to 1000°C (Fig. 13). Although the simplified treatment of instantaneous temperature and density change could not generate any clear

Rifting phase βc
Rifting phase βs1 βs2 βs3 βs4 Rifting phase βs1 βs2 βs3 βs4 Rifting phase βs1 βs2 βs3 Time/Ma uplift (Fig. 14) (slight thermal uplift has been balanced by an increase in material density), as expected for the thermal event due to magmatism, geological observations do indicate that at the end of the Middle Miocene, a clear unconformity was developed in the shallow NW basin (Su et al., 2011).
The results show that the integrated thermal effect of 1000°C magmatic injection in the Pliocene can account for the observed high heat flow (Figs 11a and 12a), and the predicted present basement depth could correspond well with the observed depth (Figs 11b, 12b, and 14a, c). However, our calculations suggest that earlier injection events in the Early Miocene through the Middle Miocene have negligible effect on the present heat flow due to rapid thermal diffusion, and younger heating event has yet to be observed at the surface. However, as in Case 2 (Fig. 10b), the predicted basement at 21 Ma is considerably deeper than the reference basement (Fig. 14a, c).  Temperature in the corresponding red area (T < 1000°C) and material density above Moho are changed respectively to 1000°C and 2950 kg m À3 at the equivalent injection time t e .
Since the timing and geometry of magmatism is unknown, Fig. 13 is a weakly constrained model. We have tried to improve the match between the predicted and observed heat flows by modifying such model parameters as the heat generation rate, thermal conductivity and initial lithosphere thickness. All of these efforts, except magmatism injection, failed to predict the heat flow as well as the observed data.
Case 5: We assume that the lithosphere was thinned uniformly by b c during the syn-rift phase. Warmer asthenosphere rises to replace cold lithosphere, buffering some of the syn-rift subsidence (Fig. 15). As in Case 4, the recent equivalent magmatism (Fig. 13) was added to match the observed high heat flow (Figs 11a and 12a).
Because our model is a kinematic numerical model, it is unable to implement the above complex processes without simplification. To reflect the influx of warmer asthenosphere and the deep thermal anomaly during the syn-rift phase, the temperature of Layer 1 at the beginning of each time step was increased with a constant warming rate W T (°C Ma À1 ), followed by thermal diffusion of a time step. W T was equal to the temperature increment DT divided by the duration of syn-rift phase (=45-21 Ma). Thus, the temperature T N of the model bottom increases from initial T o to T o +DT at the end of the syn-rift phase (21 Ma). To match the referred basement depth at 21 Ma, Layer 1 should be 360°C (DT) warmer than the initial temperature of 1300°C (T o ). However, the influx of warmer asthenosphere might cause thermal-mechanical erosion (similar to depth-dependent lithospheric thinning) in the lowermost Layer 2 (Burov et al., 2007;Brune et al., 2013). Thus, if the assumed thermal erosion zone is heated synchronously with Layer 1 with a warming rate W ET that decreases upward from W T to zero across the erosion zone (whose thickness h e is 20 km in profile AB, 15 km in profile CD), a smaller DT, 240°C, could generate the observed syn-rift deficit (Fig. 14b, d). The subsequent normal thermal relaxation could not produce rapid  post-rift subsidence. However, if we consider the cooling of the asthenosphere and thermal erosion region due to the decay of deep thermal sources and rapid thermal diffusion by magmatism, the rapid post-rift subsidence could be recovered (Fig. 14b, d) by decreasing the temperature linearly in Layer 1 and the erosion zone. The cooling rate C T in Layer 1 is set to be equal to the temperature difference between T N and the present deep temperature 1390°C divided by the cooling duration Dt, which is equal to (T N -1390)/Dt (°C Ma À1 ). Here, Dt is set to be 15 Ma by assuming that the cooling began at 15 Ma. The cooling rate C ET in the erosion zone decreases upwards from C T to zero across the thermal erosion zone. The subsidence matching would be improved if the cooling processes could be completed more elegantly. Our calculations show that the model could not predict the observed present high heat flow only when the recent equivalent magmatisms of Case 4 are included (Figs 11a and 12a).
The results of the above experiments could be summarized as follows: The stretching factor b s is considerably less than b c , and the present basement depth could be predicted fairly accurately by assuming that the lithosphere was thinned with b c ; A new rifting after 21 Ma is not supported, which indicates that the syn-rift subsidence might be in deficit; The observed high heat flow in the central depression might be related to magmatic injections in the Pliocene; The influx of warmer material, combined with thermal erosion, could generate the syn-rift subsidence deficit, and the rapid post-rift subsidence might be due to rapid cooling of the asthenosphere and thermal erosion caused by the decay of deep thermal sources and rapid thermal diffusion by magmatism.

DISCUSSION Origin of the high heat flow and magmatic activities
The above analyses suggest that the high heat flow in the central depression zone, particularly in the eastern high heat-flow zone, might be affected by local thermal events, such as magmatism occurring during the Miocene and Pliocene. The observed heat flow across the Xisha Trough (Nissen et al., 1995a) (Fig. 16) show that the seafloor heat flow in the trough is not as high as expected, indicating that the high heat flow in the eastern Qiongdongnan Basin is not the result of the lateral heat transfer of the NW sub-basin. The fact that recent magmatism could greatly increase heat flow has been observed at the ocean-continent transition of the Eastern Gulf of Aden. Lucazeau et al. (2008) reported a heat flow of approximately 900 mW m À2 observed over a volcanic edifice, and they proposed that the high heat flow was the result of the latest volcanic activity occurring at 0.1 Ma. In the South China Sea and its adjacent area, magmatic activities have become extensive since the cessation of seafloor spreading at 16 Ma, particularly since the Late Miocene . In the Qiongdongnan Basin, particularly in the eastern part, some igneous injections have been identified on seismic profiles (e.g. Ji et al., 2014;Fig. 1b), some of which could be identified directly from seafloor topography (inset seismic images in Figs 1b, 7a). The numerical analyses (Cases 4 and 5) show that equivalent magmatic injection in the Pliocene could predict the observed high heat flow quite well (Figs 11  and 12). Nissen et al. (1995a) also preferred recent (5-1.5 Ma) igneous intrusions, rather than hydrothermal activities, to interpret the local high heat flow on the northern region of profile EF (Figs 4a, 16). The numerical analyses of Wang et al. (2014) and Cases 4 and 5 (Figs 11 and 12) suggested that the heat contribution of magmatic thermal event could be up to 25 mW m À2 in the eastern Qiongdongnan Basin and 10 mW m À2 in the west, which indicates that the magmatism in the east might be more active.
The extensive magmatism might be fed by a deep thermal source. Geochemical and petrological evidence, such as the existence of extensive synchronous OIB-type basalts (Tu et al., 1991;Hoang & Flower, 1998;Zou & Fan, 2010;Wang et al., 2013), high mantle potential temperature (Hoang & Flower, 1998;Wang et al., 2013) and geochemical signatures of the basalts from a lower mantle plume (Zou & Fan, 2010;Wang et al., 2013), suggest that the extensive magmatism in the South China Sea and its adjacent area since the Miocene is closely related with a deep thermal source. Tomography studies suggest that there is a low-velocity structure in the lower mantle beneath Hainan Island, called the Hainan mantle plume (Lebedev & Nolet, 2003;Huang & Zhao, 2006;Zhao, 2007;Lei et al., 2009). By receiver function analyses, Wang & Huang (2012) showed that there is a 25-km thinner mantle transition zone than the global average thickness within an area approximately 200 km in diameter beneath NE Hainan, indicating an excess temperature of 180°C. Extensive Neogene magmatism and high heat flow suggest that there might be a failed thermal upwelling which linked the deep thermal source with the eastern Qiongdongnan Basin. Recent observations and studies suggest that the emplacement of large quantities of basalt  Fig. 16. Observed heat-flow (solid circles) variation with latitude along heat-flow profile EF (see Fig. 4a for location). The black curve shows the seafloor below the sea surface. The thick grey line represents a heat flow of 93 mW m À2 , corresponding to the heat flow of 30 Ma oceanic crust based on the GDH1 model of Stein & Stein (1992).
can accommodate extension with little or no crustal thinning (e.g. Buck, 2004;Bialas et al., 2010;Ebinger et al., 2013). The magmatic injections suggest that during the Neogene, there might be a failed extension event with little or no crustal thinning in the eastern Qiongdongnan Basin.

Original analyses of anomalous tectonic subsidence
The above analysis of tectonic subsidence reveals a two-stage anomalous subsidence feature of the syn-rift subsidence deficit and the well-known rapid post-rift subsidence after 10.5 Ma. Particularly, the present basement depth could be satisfactorily predicted by stretching factor b c (Figs 2, 10b, e, 11b, and 12b), suggesting that the excess post-rift subsidence occurs to compensate for the syn-rift subsidence deficit. One more clear feature is that the rapid post-rift subsidence occurred after 10.5 Ma, approximately 10 Ma after the end of the synrift phase (Fig. 2), indicating that tectonic subsidence was still in deficit during the early post-rift phase. This is different from the lithospheric depth-dependent thinning (DDT) model (Royden & Keen, 1980). Although the DDT model could reproduce the syn-rift subsidence deficit, it predicts a rapid post-rift subsidence immediately following the end of rifting. Further crustal thinning due to a new rifting episode or lower crustal flow after the end of the Middle Miocene was used to account for the rapid post-rift subsidence (Yuan et al., 2009;Lei et al., 2013;Zhao et al., 2013). However, the results of Case 3 and the geological observation that faulting became weak after 21 Ma suggest that further crustal thinning due to a new rifting after 21 Ma is unlikely. The mechanism of the lower crustal flow model (Morley & Westaway, 2006) was thought to be against the gravitational buoyancy forces and unlikely to work (Allen et al., 2004) because it requires lower crustal material flowing from the basin area (lower pressure) to the sediment source area (higher pressure). Furthermore, based on the high-quality crustal structure sections obtained from seismic refraction data across the study area (Qiu et al., 2001;Huang et al., 2011), the thinning factors of the whole crust and the upper crust (Fig. 17) show that the crust was stretched nearly uniformly, which indicates that the lower crustal flow or crustal depth-dependent thinning was unlikely to occur in the study area. Therefore, the rapid post-rift subsidence in the Qiongdongnan Basin might be caused by other mechanisms. Because the present basement depth is satisfactorily predicted by stretching factor b c (Figs 2, 10b, e, 11b, and 12b), the mechanisms accounting for syn-rift subsidence deficit might be interrelated with those resulting in rapid post-rift subsidence. These mechanisms generate temporary subsidence anomalies. We employ a lithospheric depth-dependent thinning, as in earlier kinematic models (Royden & Keen, 1980). The weak syn-rift magmatism on the northern margin of the South China Sea implies that during the syn-rift phase, the continental sub-lithospheric upper mantle was relatively cold, possibly colder than the potential temperature T p of 1300°C (Reston & Morgan, 2004;Yan et al., 2006). Rifting on such a cold sub-lithosphere mantle could lead to the influx of warmer materials, which would cause thermal uplift buffering syn-rift subsidence (Reston, 2009). Small-scale convection may have been induced by the steep gradients at the transition between thinned and unthinned lithosphere during the syn-rift stage (e.g. Ebinger et al., 2013) (Fig. 7). Geochemical analyses on the late Cenozoic basalts in the Hainan hotspot suggest that the mantle potential temperature T p is 1440-1550°C . The numerical analyses of Case 5 indicate that the predicted syn-rift subsidence at 21 Ma would match the observed if the asthenosphere is 240°C warmer than at the onset of rifting (Figs 12 and 14). In other words, an increase in potential mantle temperature T p of 1440-1550°C could cause the syn-rift subsidence deficit. According to the low-temperature thermochronological data, Shi et al. (2011) suggested that the eastern Hainan Island also went through a rapid cooling episode during the Late Oligocene even though it is far from large active faults, indicating that the Late Oligocene rapid uplift and cooling episode could be caused by a mantle upwelling. These studies suggest that during the late syn-rift phase, Fig. 17. (a) Crustal structure sections based on seismic refraction data in the study area (see Fig. 2 for locations) and (b) Upper crust (UC) thinning factors & whole crust (WC) thinning factors of the two sections, thinning factor = 1.0 -Hp/Ho, where Hp is the present thickness of the upper crust or the whole crust given by the crustal sections in (a), and Ho is the initial thickness of the upper crust or the whole crust. The upper crust Ho is set to be 16.0 km, and the whole crust, Ho, is 32 km. The structure of H1H6 is from Huang et al. (2011), the other one of Q1Q10 is from Qiu et al. (2001). the observed deep thermal anomaly may have thus buffered syn-rift tectonic subsidence and caused the syn-rift subsidence deficit in the study area (Fig. 18). Unlike the lithospheric depth-dependent thinning (DDT) model (Royden & Keen, 1980), the rapid post-rift subsidence in the Qiongdongnan Basin occurred approximately 10 Ma after the end of the syn-rift phase (Fig. 2), suggesting that the dynamic support generated by the deep thermal anomaly below the study area might persist during the early post-rift phase (Fig. 18). Lucazeau et al. (2008Lucazeau et al. ( , 2009 argued that the persistence of thermal anomalies after the continental break up in the eastern Gulf of Aden was likely a common and long-term feature of continental margins. When the deep thermal anomaly decays or moves away and can no longer maintain sufficient dynamic support to buffer the subsidence, rapid subsidence occurs to compensate for the earlier subsidence deficit (Fig. 18). The rapid post-rift subsidence in the eastern Qiongdongnan Basin might be the result of decay of a deep thermal source and rapid cooling of the asthenosphere. The extensive magmatism and high heat flow suggest that there might be a failed small-scale thermal upwelling which linked the deep thermal source with the eastern Qiongdongnan Basin. However, this small-scale thermal upwelling is usually short-lived (Davies & Bunge, 2006), and it may wane with the decay or moving away of the deep thermal source. Geochemical analyses suggested that the recent (<0.35 Ma) upwelling rate is <1 cm yr À1 for the Hainan plume (Zou & Fan, 2010), whereas a longer average long-term upwelling rate is >1 cm yr À1 , which indicates that the deep thermal anomaly beneath Hainan may have been decaying. Furthermore, denser magmatic intrusions into the crust could increase the load, increasing subsidence (Shi et al., 2005). The simplified numerical analyses of Case 5 show that the combination of these processes could satisfactorily predict the preferred subsidence path and could generate rapid post-rift subsidence (Figs 11,12,and 14) to compensate for the earlier subsidence deficit (Fig. 18). Besides the eastern Qiongdongnan Basin, the Yinggehai Basin and the western Qiongdongnan Basin also show very similar rapid post-rift subsidence phenomenon, supporting our regional dynamic mantle model.

CONCLUSIONS
1 Tectonic subsidence analyses show that in the Qiongdongnan Basin, the stretching factor b s based on synrift sequences is considerably less than the observed crustal stretching factor b c , and the present basement depth could be predicted fairly accurately, assuming that the lithosphere was thinned with b c . Syn-rift subsidence at 21 Ma is in deficit, and post-rift subsidence after 10.5 Ma is anomalously rapid. 2 Heat-flow data show that heat flow in the central depression is 65-105 mW m À2 , which is considerably higher than that on the northern shelf. The observed high heat flow in the central depressions was caused by the heating of a recent magmatic injection, which affected the eastern basin more than the western basin. The extensive magmatism might be fed by a deep thermal source. 3 The observed deep thermal anomaly might have affected the study area since the Late Oligocene. The syn-rift subsidence deficit might be mainly caused by the dynamic support of the deep thermal anomaly beneath the study area, and the rapid post-rift subsidence might result from the decay or moving away of the deep thermal source and the rapid cooling of the asthenosphere. The rapid post-rift subsidence occurred to compensate for the syn-rift subsidence deficit.

ACKNOWLEDGMENTS
We thank the crew of R/V Shiyan 3 and the scientists who participated in the fieldwork. We are grateful to Dr. Ruiqing He and Dr. Bin Zhu for improving the readability of the original manuscript, and Prof. Cynthia Ebinger and two anonymous reviewers for their constructive and valuable comments that improve greatly the quality of the manuscript. Figures 1, 2 and 4 are plotted using GMT (Wessel & Smith, 1995). The geothermal features of the deep-water area were analysed mostly based on the heatflow data of Prof. Xing Xu of the Guangzhou Marine Geological Survey. This work was jointly supported by the National Natural Science Foundation of China (Grant 41176050, 41576036, 41376059). The solid line represents the variation in the dynamic support from the deep thermal anomaly, and the dark grey region represents the magmatism based on the age histogram of Cenozoic volcanism .