Post-collisional exhumation and geotherm pattern in northern Tananao Complex, northeastern Taiwan

The northeast region of Taiwan has experienced a stress transition from the compressive collision to the extension induced from the Ryukyu subduction system in a post-collisional orogenic belt. In order to clarify when the extensional stage began and how the cooling and exhumation history varied with the transition of tectonic stress, apatite and zircon fission-track ages on the granitic gneiss in northern Tananao Complex were collected from a ~2700 m high sub-vertical profile. The exhumation rate was estimated as ~4 mm yr -1 since 4.1 - 3.0 Ma using previous biotite 40 Ar- 39 Ar ages with an assumed uniform geothermal gradient. During the period of 1.5 - 0.2 Ma, the exhumation rate remained at 2.3 - 1.6 mm yr -1 , evidenced by the similar slopes for apatite and zircon in the age-elevation relationship (AER). The decrease in exhumation occurred in between 3.0 and 1.5 Ma is interpreted as the onset of the extensional im-pact into NE Taiwan from the Ryukyu subduction zone. After ~0.2 Ma, the intercepts of apatite AER suggest that exhumation accelerated to ~5 mm yr -1 at the mountain range but was steady at the valley bottom. The offset in ages of the vertical profile is supportive of the normal faulting at the shallow crust. The accelerating cooling with the steady exhumation from 1.5 - 0.2 Ma implies the increase in paleogeothermal gradient with time. The paleogeothermal gradients are estimated to be over 100°C km -1 after 1.5 Ma, suggesting that the Hoping area possesses heat sources induced from the complex tectonic configuration. However, the modern geotherm is relatively low and may be influenced by the water circulation near surface.


INTRODUCTION
The well-known Penglai Orogeny of Taiwan mountain building belt is a result of the oblique collision between the northern Luzon volcanic arc of the Philippine Sea Plate and the continental margin of the Eurasian Plate since Late Miocene (Ho 1986;Teng 1990). Taiwan Island is distinctively located above two subduction systems with a flipping of polarity: to the south offshore, the Eurasian Plate subducts eastwards underneath the Philippine Sea Plate along the Manila Trench while to the eastern offshore, the Philippine Sea Plate subducts northwards below the Eurasian Plate along the Ryukyu Trench (Fig. 1a). At present, the Philippine Sea Plate is moving northwestwards with a rate of ca. 8 cm yr -1 (Yu et al. 1997). The development of Tai-wan orogen has been proposed showing the progressively southward propagation of the arc-continent collision at rates ranging from 50 -90 mm yr -1 (Suppe 1981;Liu et al. 2001). As a result, different geographic domains of Taiwan display asynchronous mountain building stages across the island: the post-collisional stage in the north, the steady-state collision in the central, and the incipient stage of collision in conjunction with advanced subduction in the south (Teng 1996;Willett et al. 2003;Stolar et al. 2007).
The Ryukyu subduction zone follows its westward direction and terminates against the eastern Taiwan mountain belt. To the west of the subduction zone, the western edge of subducted Philippine Sea plate (WEP) runs beneath northern Taiwan (Teng 1996;Wu et al. 2009). To the northeast of WEP, the northern part of the orogen has been subject to crustal subsidence and extensional stress resulting from Terr. Atmos. Ocean. Sci., Vol. 31, No. 4, 369-381, August 2020 the post-collisional collapse since ~2 Ma (Teng 1996;Teng et al. 2000). A series of post-collisional volcanoes in Pliocene-Quaternary have been proposed to be related to the southwestward migration of the Ryukyu Trench (Wang et al. 1999). In addition, to the northeast offshore of Taiwan, the north-dipping subduction of the Philippine Sea Plate proceeds and is accompanied by the rifting of the Okinawa Trough due to the southward rollback of the Philippine Sea slab (Teng 1996). The extensional phase of the southern Okinawa Trough has been proposed to occur since Late Pleistocene (~0.1 Ma) based on seismic refraction data (Sibuet et al. 1995(Sibuet et al. , 1998, characterized by normal faults with E-W strike around the Ilan Plain. The Ilan Plain is located at the western tip of the N-S opening Okinawa Trough and is presently subsiding at rates between 2 and 18 mm yr -1 (Ching et al. 2011). However, the seismogenic patterns indicate that the Philippine Sea slab has remained dominated by NW-SE compression towards the Eurasian Plate at a greater depth (Kao and Jian 2001;Wu et al. 2009).
The tectonic stress of northeastern Taiwan has been regarded as a spatial stress transition from the compressional state into the extensional stage (Shyu et al. 2005;Huang et al. 2012). The transition pattern can be illustrated by the regional seismicity ( Fig. 1b) with earthquakes whose magnitudes are 3 -6 and depths are shallower than 100 km from 2010 -2015 (Institute of Earth Sciences, Academia Sinica, Taiwan 1996). The distributions of focal mechanisms show the change of fault types from south to north, indicating the transition of the regional stress field. To the south where the Coastal Range locates, a cluster of thrust faults parallel to the subducting interface between the Philippine Sea Plate and the Eurasian Plate occurs. In Suao-Hoping area, it is dominated by left-lateral strike-slip faults with a NNW-SSE minimum principal stress axis, consistent with the rifting of the Okinawa Trough (Huang et al. 2012). To the north, the Ilan Plain consists of a series of normal faults.
Although the geological structures, geodetic measurements, and seismogenic patterns all indicate a transition of tectonic stress in northeastern Taiwan, the link between the Ryukyu subduction system and the exhumed Taiwan mountain belt has not been well constrained by the thermochronolgic data so far. In order to clarify the influence superimposed by the extensional stress induced from the Ryukyu subduction system on the exhumation history of northeastern Taiwan, we collected and analyzed 15 rock samples from a sub-vertical profile of the Tachoshui (a) (b) Fig. 1. (a) The regional tectonic map of Taiwan and the adjacent regimes. The tectonic configurations follow Teng (1996) and the current velocity vector of the Philippine Sea Plate is relative to the Penghu Islands (Yu et al. 1997). The western edge of subducted Philippine Sea plate (WEP) beneath northeastern Taiwan is recognized as the boundary between the Eurasia plate and the Philippine Sea plate (Wu et al. 2009 Ho 1986). (b) The regional stress tensor patterns proposed by Huang et al. (2012) and the distribution of focal mechanism around the hinterland and offshore of northeastern Taiwan granitic gneiss body in the northern Tananao Complex. The average exhumation rate is derived from the fission-track age-elevation relationship (AER) without the presumption of a fixed geothermal gradient. The comparisons between the modern geothermal gradient and the available thermochronologic data further enable us to detect the perturbation of paleo-isotherms and the accelerated cooling rates in the northern Tananao Complex.

GEOLOGICAL BACKGROUND
From west to east, the geological units of Taiwan can be divided as the Coastal Plain, mostly composed of Holocene clastic sediments; the Western Foothills, the sedimentary unit mostly comprising Miocene to Pleistocene sedimentary strata; the Hsueshan Range, a low grade metamorphic belt consisting of Eocene to Miocene metamorphic rocks; the Backbone Range, a slate belt mainly composed of Eocene and Miocene strata; the Tananao Complex, a severely metamorphosed complex consisting of blackschist, greenschist, meta-chert marble and gneiss; and the Coastal Range, the remnant of the Luzon Arc and the fore-arc basins (Fig. 1a).
The target Hoping area is geographically positioned along the tip of WEP and under a post-collisional environment with complex stress derived from the deep slabcontinent collision, the subduction of the Ryukyu Trench, and also the rifting of the Okinawa Trough (Fig. 1a). When crossing the Chingshui Cliff, the Hoping area faces a sudden drop in altitude and enters into one of the adjacent fore-arc basins, the Hoping Basin. The Hoping Basin is a triangular graben surrounded by normal faults and dextral strike-slip faults resulted from the extensional stress above the Ryukyu subduction zone (Shyu et al. 2005). The significantly negative free-air gravity anomaly over the Hoping Basin probably results from the presence of thick sediments or subsidence of the oceanic crust, both indicating a downwarping of the WEP against the eastern side of the island (Yen et al. 1995).
The Tananao Complex has experienced multiple episodes of metamorphism and deformation as summarised by Simoes et al. (2012). The studies on K-mica crystallinities have suggested that it has explicitly archived the last metamorphism in greenschist facies with the calculated P-T condition of ~5 kb and 350 -475°C during the Penglai Orogeny (Chen and Wang 1995). The Hoping area sits at the northern tip of the Tananao Complex and is mainly composed of Cretaceous granitic gneiss bodies, quartz schist, graphite schist, and marble. The granitic gneiss bodies (e.g., the Chipan and the Tachoshui granitic gneiss body) have been considered subduction-related plutonic rocks in origin (Lan et al. 1996). The Tachoshui granitic gneiss body is affected by a series of anticlines and synclines which are parallel to the orientation of the Taiwan mountain belt (Fig. 2a). The mylonitic deformation in the Tachoshui granitic gneiss body was found on the valley of the Hoping River and suggested to be associated with the shortening and thickening of the mountain building during the Penglai Orogeny (Wang et al. 1998). Prominent foliations occur in the studied granitic gneiss body with a strike of N 66 -115° and a moderate dip around 35° toward north. The foliation features were later truncated by pseudotachylyte veins (Chu et al. 2012). These pseudotachylyte veins were recognized with normal faulting kinematics, indicating a stage of exhumation-related extension in Pliocene-Pleistocene (Korren et al. 2017).
Geochemical and thermochronological data suggest that the Taiwan mountain belt has undergone the most critical thermal event (i.e., Penglai Orogeny) since the late Cenozoic. The Cretaceous granitic gneiss body in the Chipan area yields a flat biotite 40 Ar-39 Ar plateau age of ~7.9 Ma (Lo and Onstott 1995). On the other hand, the mylonitic gneiss at the bottom of the Hoping River ( Fig. 2a) shows the biotite 40 Ar-39 Ar ages ranging from 4.1 -3.0 Ma (Wang et al. 1998). These datasets suggest that the temperature of the last metamorphism could reach 350°C, consistent with the occurrence of greenschist facies in the northern Tananao Complex. In addition, abundant results of fission-track dating in the Tananao Complex show that the apatite and zircon grains have been totally annealed by the last thermal event with the cooling ages of ca. 0.6 -0.3 and 2.0 -0.8 Ma, respectively (Liu 1982;Tsao 1996;Willett et al. 2003;Fuller et al. 2006;Lee et al. 2015;Hsu et al. 2016). The (U-Th)/He dating for apatite recorded a very young age of 0.25 Ma in the southernmost tip of the Tananao Complex (Hsu et al. 2016).

METHODOLOGY
Temporal variations in exhumation rates can be estimated by constraining the cooling rates of different systems with high to low closure temperatures in the same rock under the assumption of a thermal steady state (Reiners and Brandon 2006). Theoretically, the average exhumation rate is inferred as dividing the cooling rate by an assumed geothermal gradient. However, if there is insufficient time to achieve the thermal steady state, the exhumation would cause the upward migration of the isotherms and increase the thermal gradient. The evolving of the thermal field through time can be complicated and mislead our interpretations of temperature distribution within the crust to transform from a linear relationship into an exponential increase of temperature with depth as Braun et al. (2006) suggested. The variation of geothermal gradient can therefore result in a misunderstanding about the exhumation history of a rock.
To avoid the uncertainty and complexity of an assumed geothermal gradient, alternatively, we can examine the AER of the bedrock to simply discuss the vertical velocity (i.e., exhumation) of the bedrock relative to each closure isotherm at the time of closure. The slope changes of AER then consequently represent the changes of exhumation rates (Farley et al. 2001;Reiners and Brandon 2006). The approach is based on the one dimensional model with three assumptions: (1) a flat isotherm at the time of closure; (2) spatially uniform erosion rates; and (3) a constant closure depth of the closure isotherm during the span of the discussed AER. These assumptions are to assure that the slope of the AER can reflect the actual movement of the bedrock to the surface through a stationary isotherm.
In addition, the value of the zero-age intercept calculated from the AER trend can predict the isotherm depth of a specific closure temperature with the aid of modern geothermal gradient data and can provide more insights on the recent exhumation history. When samples of a vertical transect pass through a flat closure isotherm at a depth below the mean elevation with a constant exhumation rate, the predicted depth of its closure isotherm from AER will accordingly archive the response to the post-closure system enhanced or reduced exhumation. If the modern isotherm is below (above) the predicted paleo-isotherm, then it suggests that the exhumation rate has increased (decreased) after the age of the youngest sample (Reiners and Brandon 2006).
The sampling strategy in this study is to retrieve sam-ples from a single rock body in a vertical profile, it is able to extract average exhumation rate and paleogeothermal gradient simply derived from the AER without assumption on the past crustal thermal structure (Gallagher et al. 1998). A total of 15 samples were collected from the Tachoshui granitic gneiss body along a former logging trail and from a borehole near the valley floor in the Hoping area. The samples display a narrowly horizontal extent of ca. 8 km and a vertical span of ca. 2700 m. Eleven granitic gneiss samples were retrieved from a sub-vertical transect at different elevations: nine of them were sampled from 2100 -840 m above sea level (masl) and two from the valley bottom (100 masl). Four samples were collected from borehole HCBH01 at depths ranging from 48 -587 m below the surface (Fig. 2b).
The experimental procedure of fission-track dating followed Liu et al. (2000) and was processed for both zircon and apatite. Most of the analyzed gains were polished to expose internal grain surfaces to give 4π geometry for fissiontrack density counting while three zircon samples (i.e., 45K, 23K, and 0K) were prepared with 2π geometry. Grain-bygrain and external detector methods were adopted to obtain individual grain ages (Galbraith 1984 1985). The age standard samples of the Fish Canyon tuff zircon with standard glass NBS SRM-612 and Durango apatite with NBS SRM-610 (Hurford and Gleadow 1977) were used to determine the zeta (ζ) values (Hurford and Green 1983) as 28 ± 5 and 345 ± 29 (1σ) yr cm -2 , respectively.

RESULTS
Apatite and zircon fission-track ages are listed with 1σ (68%) confidence interval in Table 1. Both apatite and zircon fission-track ages are significantly younger than the formation age (ca. 90 Ma) of the Tachoshui granitic gneiss body (Jahn et al. 1986). All results passed the χ 2 test, indicating that the fission-tracks have been totally annealed by the last thermal event during the Penglai Orogeny.
Both the apatite and zircon fission-track ages display the good linear correlation with respect to elevations and depths (Fig. 3). The sample acquired from the bottom valley, however, exhibits an unexpected age that is older than the other surface samples in higher elevation. We had collected duplicate samples from the valley bottom (i.e., HPV and 0K) to re-examine the fission-track ages at the same elevation and they both display older ages (0.6 -0.5 Ma for apatite and 1.2 -1.1 Ma for zircon). When compared with the other surface samples, the two results from the bottom of the valley lie outside of the 68% confidence interval and show no correlation coefficient (r 2 < 0.1) in the linear regression analysis. Meanwhile, when compared with borehole samples, they fit into the age-depth trend and obtain a strong linear relationship for both apatite and zircon fissiontrack results. In the following context, we will present and discuss the results as two major groups: the surface group of samples collected from higher elevation of 2100 -840 masl; and the lower group retrieved below 100 masl, including valley bottom and borehole.
The uncertainties for the slopes and intercepts of AER are also shown with ±1σ confidence interval. The apatite fission-track results of nine gneiss outcrops from the subvertical transect range from 0.7 ± 0.1 (at 2100 masl) to 0.3 ± 0.1 Ma (at 840 masl) with the AER slope of 2.3 ± 0.8 mm yr -1 (r 2 = 0.6) and the intercept of 231 ± 406 m, as determined by the least squares regression. Samples from the valley bottom and the borehole range from 0.6 ± 0.1 (at 100 masl) to 0.2 ± 0.1 Ma (at the depth of 402 m) with an age-depth slope of 1.6 ± 0.6 mm yr -1 (r 2 = 0.7) and the intercept of -774 ± 242 m.
Zircon fission-track ages display a very similar trend to the apatite fission-track results. Zircon fission-track ages from the exposed gneiss outcrops range from 1.5 ± 0.1 to 1.0 ± 0.1 Ma at the elevation of 2100 -840 masl with the AER slope of 2.1 ± 0.6 mm yr -1 (r 2 = 0.7) and the intercept of -1467 ± 821 m. Zircon fission-track ages from the valley bottom and borehole range from 1.2 ± 0.1 (at 100 masl) to 0.8 ± 0.1 Ma (at the depth of 587 m) with an age-depth slope of 1.8 ± 0.3 mm yr -1 (r 2 = 0.9) and the intercept of -1979 ± 359 m.

Exhumation History Between 4.1 and 0.2 Ma and Offset in Ages of Vertical Profile
Since 4.1 -3.0 Ma, the exhumation rate was estimated as 3.9 -2.6 mm yr -1 on the basis of the biotite 40 Ar-39 Ar ages from the mylonitic gneiss found at the valley bottom of the Hoping River (Fig. 2a) (Wang et al. 1998), whereas a surface temperature of 15°C, a closure temperature of 370 -330°C, and a geothermal gradient of 30°C km -1 based on metamorphic P-T conditions of the greenschist facies in the Tananao Complex are assumed (Chen and Wang 1995). In this study, the AER slopes of the apatite and zircon fissiontrack results display almost identical trends within the regression error bars for both surface and borehole samples from ca. 1.5 -0.2 Ma, suggesting that the Tachoshui granitic gneiss body was exhumed in a steady pace between 2.3 and 1.6 mm yr -1 during this period so that the two minerals could pass through their closure temperatures with a similar exhumation rate (Fig. 3). Comparing the aforementioned exhumation rates derived from the high-to low-temperature thermochronometers, the exhumation rate in the Hoping area seems to decline in between 3.0 and 1.5 Ma.
The exhumation rates have been estimated as 6 -4 mm yr -1 since ca. 5 Ma for the steady-state section of central Taiwan , consistent with the rates of 6 -3 mm yr -1 within the Central Range, using the one-dimensional thermal modelling (Dadson et al. 2003). A rapid exhumation was recorded as ~6.3 mm yr -1 in the Tananao Complex, estimated by the zircon (U-Th)/He data and thermokinematic model ), while Liu (1982) estimated the exhumation rates as 9 -7 mm yr -1 in the Central Range from the apatite and zircon fission-track ages. In addition, Hsu et al. (2016) have also reported several age-elevation profiles in the central-southern Central Range suggest the exhumation rates from 4 -2 mm yr -1 during ca. 2 -0.5 Ma and a further increase of 8 -4 mm yr -1 from 0.5 Ma to present. The significantly rapid exhumation in the eastern part of the Taiwan mountain belt is in consequence of the intense compression from the arc-continent collision. However, the exhumation rates in the Hoping area of 2.3 -1.6 mm yr -1 are temporally and spatially 2 to 3 times slower than those estimated in the central-southern part of Taiwan mountain belt. The relatively slow rate in the Hoping area probably coincides with the initiative extension stress above the Ryukyu subduction system to the east offshore of Taiwan.
Despite that the linear regression error bars of AER are statistically large, the surface samples (collected from 2100 -840 masl) and borehole samples (collected below 100 masl) did exhibit the positive correlations between fission-track age as a function of elevation and depth. The apatite and zircon fission-track ages from the valley floor are older than the ages projected from the surface AER trend and the old -test with probability P(χ 2 ) value > 5%, indicating a single thermal event (Galbraith 1981). ages are confirmed by the reduplicated samples (HPV and 0K). This exceptional and unique characteristic of the offset in the fission-track age occurs between sample HP08 (840 masl) and HPV/0K (100 masl). These offset in ages could be induced by the tectonic or topographic perturbation of the thermal field. The topographic effect on the fission-track closure isotherm can lead to a steeper AER slope of apatite than that of zircon and give an overestimate of the real exhumation rate (Stüwe et al. 1994). However, our data show similar AER trends for both apatite and zircon fission-track results, implying that the topographic effect on shallower isotherms may not significantly disturb the surface samples and borehole samples. A simple kinematic solution is to propose a fault or landslide boundary in between the two terranes. One of the possibilities appears to be a SW-dipping thrust fault system. In this case, the samples higher than 840 masl are considered to be located on the hanging-wall and being thrust up towards northeast (Fig. 4a). This thrust mechanism probably attributes to the NW-SE compression induced from the western margin of the Philippines Sea Plate slab colliding with the continental margin. Beneath the Central Range in northeastern Taiwan, the Eurasian Plate is likely to vertically contact against the Philippines Sea Plate slab. Offshore of northeastern Taiwan, the plate interface of the subducted slab is characterized by the NW-trending shear zones with sinistral component (Lallemand et al. 2013). Because the differential motion between the plates on the two sides, the seismic stress around the Hoping area shows NW-SE compression at the deeper depth and transits to NW-SE tension toward the shallower crust at < 35 km (Wu et al. 2009). As the thrust faulting is proposed to be placed after ~0.2 Ma, it only considered a few kilometer depths. The thrust faulting seems less likely to form within the very shallow crust given independent evidence of stress tensor.
The other option is a NE-dipping normal fault system along the pre-existing foliation planes on the Tachoshui granitic gneiss body. In this case, all the samples below 100 masl are considered to be located on the hanging-wall and hence may preserve older fission-track ages even though they are lower than the foot-wall in elevation (Fig. 4b). The possibilities of normal faulting activities in the surrounding area of northeast Taiwan have been reported. For instance, the occurrence of pseudotachylyte veins is believed to be associated with the exhumation-related extension during ~1.5 Ma (Korren et al. 2017;Chen et al. 2018). As Dewey (1988) elaborated, gravitational sliding or gravitational spreading at shallow crust often occurs in the post-collisional orogens when the lithosphere fails to maintain stability. This post-collisional extension process could be behind the exhumation of the mountain range and alternatively expose the high-grade metamorphic basement of Taiwan, e.g., the amphibolites and gneiss. Continuous detection of levelling survey (Ching et al. 2011) and sea-level variation ) demonstrate modern subsidence along the east coast of Taiwan, illustrating the ongoing influence of extensional stress. Considering the stress pattern of regional tectonics around northeastern Taiwan, the possibility of forming the normal fault system in the extension-dominated shallower crust seems to be more appealing. Although there is insufficient field evidence to date that provides a robust trace of the discussed fault systems, we expect more detailed identification of fault activities in the future to make furthermore reasoning on the kinematics of the structure.

The Present and Past Closure Depth: Implications for Recent Exhumation
To further constrain the post-0.2 Ma exhumation history, we take advantage of comparison between the present-day closure depth below the mean elevation derived from the modern geothermal gradient and the zero-age intercept of apatite fission-track AER trend. The topography of the Hoping and its nearby area shows a sudden decrease in elevation towards north and east, and hence the mean elevation of the first-order topography is adopted when considering the modern closure depth of apatite fission-track. A mean elevation of 1250 m according to a 50 km long topographic profile along the coastline of northeast Taiwan and a modern geothermal gradient of 56°C km -1 (Hsieh et al. 2014) are applied to describe the modern relief. The modern isotherm of the closure temperature for apatite fission-track of 130°C should be situated at the depth of 800 m below the sea level (surface temperature is 15°C). The extrapolated isotherm positions corresponding to the closure temperature of apatite for the surface and borehole samples can be estimated with their AER lines at 231 and -774 m, respectively (Fig. 4c). For the surface samples, the discrepancy on the position of the predicted isotherms (from 231 to -800 m) suggests an exhumation of the upper terrane about 1 km after the closure of the youngest sample at 0.2 Ma, equivalent to an increased exhumation rate of ca. 5 mm yr -1 .
On the other hand, the predicted isotherm of the lower terrane is close to the one estimated from the mean elevation (-774 and -800 m), suggesting that the borehole samples did not experience severe perturbations of enhanced exhumation and have maintained a relatively stable rate since 1.5 Ma.
Given the short horizontal distance between the surface and borehole samples, the enhanced exhumation observed is likely the result of the proposed faulting, and the measured 1 km total exhumation might be equivalent to the vertical offset of the fault. This is also supported by the vertical shift between the topmost samples of the bottom profile (at 100 masl) and its age-equivalent (at 1400 -1500 masl) in the upper profile. However, the contribution of the fault slip would require a better understanding of the fault geometry to quantify in the future. In order to constrain the exhumation behaviour of the proposed faulting, the timetemperature histories of the upper terrane and lower terrane are analyzed by the inversion model of HeFTy v.1.8.0 (Ketcham 2005). The depositional age of the northern Tananao Complex is set to Late Eocene (38 -34 Ma) as suggested by Hsu et al. (2016). Two thermal constraints are superimposed in the modelling of thermal history: a peak metamorphic temperature is set to 550 -450°C based on RSCM dataset ) at 6.5 -6.0 Ma and a thermal event of 370 -330°C at 4.1 -3.0 Ma according to the biotite 40 Ar-39 Ar dating (Wang et al. 1998). (c) Fig.4. (a) The schematic diagram showing the difference in paleo-isotherm depth between foot-wall and hanging wall due to the inferred SWdipping thrust faulting. The samples from a vertical profile passed through a flat closure isotherm at a depth (Zc) below the mean elevation with a constant exhumation rate. (b) The diagram represents the the inferred NE-dipping normal faulting. (c) The AER plot of the apatite fission-track in the Tachoshui granitic gneiss body and the estimated closure depth of 130°C isotherm at 2050 m beneath the mean elevation with the assumptions that the surface temperature is 15°C and the modern thermal gradient is 56°C km -1 as reported by Hsieh et al. (2014). The closure depth of modern isotherm is below the predicted paleo-isotherm, suggesting an increased exhumation in the upper terrane after ~0.2 Ma.
The inverse modelling shows that all weighted mean cooling paths are overlapped within a region and thereby we are unable to distinguish the differences in cooling histories between the upper terrane and the lower terrane (Fig. 5). The similar cooling paths for all samples suggest that both apatite and zircon grains passed through their fission-track closure temperature at a steady exhumation in 1.5 -0.2 Ma, evidenced by the slopes AER (Fig. 3). However, the timetemperature paths for the Tachoshui granitic gneiss body presented an accelerated cooling history after ca. 1.5 Ma. The temperature-time path could appear to be an exponential decrease in temperature with time even with a constant exhumation rate. The conflict may be due to the non-linear correlation between the geothermal gradient and the depth because the isotherms tend to be more compressed toward the shallower crust as a result of heat advection during exhumation. The geothermal gradient is critical in the subsurface thermal structure and is a function of thermochronologic age (Mancktelow and Grasemann 1997;Braun et al. 2006), thus the age is used to estimate the rates of cooling and exhumation. The intercepts of apatite AER show a stage of enhanced exhumation after ~0.2 Ma (Fig. 4c), possibly related to the proposed faulting. However, the rapid acceleration in cooling rate after ~0.2 Ma is not predicted in the HeFTy results (Fig. 5). The modelling results appear to be strongly controlled by the input of older fission-track ages, such as 0K on the lower terrane and 45K on the upper terrane, which generate relatively slow cooling histories. Acquisition of apatite fission-track length and He-dating would be greatly helpful to clarify the more recent exhumation process in the future.

Calculation of Paleogeothermal Gradient in the Hoping Area
Thermochronologic data with low to high closure temperature systems allow us to estimate the cooling rates in different time spans and the concurrent paleogeothermal gradients when samples passed through their paleo-closure depth. The closure temperature of biotite 40 Ar/ 39 Ar system is set to 370 -330°C from 4.1 -3.0 Ma as previously reported by Wang et al. (1998). The inversion modelling demonstrates that the range of closure temperature for zircon fission-track system is 250 -225°C from 1.5 -0.8 Ma. Later on, the range of closure temperature for apatite fission-track system is 140 -110°C from 0.7 -0.2 Ma (Fig. 5). In consequence, the cooling rates from biotite 40 Ar/ 39 Ar to zircon fission-track system can be estimated as 55 -40°C m.y. -1 (50°C/my on average). Likewise, the cooling rates from zircon to apatite fissiontrack system can be estimated as 190 -140°C m.y. -1 (160°C/ my on average), and for the range of apatite fission-track system to the surface (15°C), it is equivalent to the cooling rates of 620 -140°C m.y. -1 (350°C/my on average). The corresponding geothermal gradients within the temperature ranges for each time span can therefore be discussed and inferred. At 4.1 -3.0 Ma, the initial paleogeothermal gradient is inferred as ~30°C km -1 based on the P-T range of the deeply buried biotite as discussed in the geologic background section. From 1.5 -0.8 Ma, the paleogeothermal gradient can be estimated to be 85°C km -1 on average when dividing the cooling rates by the observed exhumation rates (i.e., 2.1 -1.8 mm yr -1 ). From 0.7 -0.2 Ma, the paleogeothermal gradient can be estimated at 180°C km -1 on average as the exhumation rates are 2.3 -1.6 mm yr -1 (Fig. 5).
In order to better incorporate the effect of the variable geothermal gradient, we applied the analytical solution proposed by Willett and Brandon (2013). The exhuming half-space modelling considers the dependence of closure temperature on cooling rate and the advection of heat upwards thereby change the geothermal gradient. It therefore provides the exhumation rates under no steady-state (transient) condition. In principle, when all ages are lying in the linear sequence with respect to elevations with a constant exhumation rate from the time of closure to the present-day, all of the resultant curves would intersect at a point in the plot of exhumation rate versus geothermal gradient. Given a mean elevation of 1250 m in the Hoping area, a surface temperature of 15°C, and the time for initiation of exhumation of 6.5 Ma, the analytical solution can accordingly predict the associated exhumation rates and geothermal gradients. Because the solution requires the surface boundary condition of topography in elevations, only nine surface samples from 2100 -840 masl are used in the modelling.
The modelled results of zircon fission-track ages (Fig. 6a) show that all age-curves generally concentrate to an area without crossing at a point even though the samples are expected to fit in the linear regression of AER. The preferred range of geothermal gradients of 120 -100°C km -1 corresponds to the exhumation rates of 2.7 -1.6 mm yr -1 (~2.2 mm yr -1 on average). On the other hand, the resultant curves of fission-track age for apatite (Fig. 6b) are more scattered than those for zircon especially for sample HP07, whose apatite fission-track age lies outside of the linear regression analysis of AER (Fig. 3). After excluding the outlier, the geothermal gradients can be estimated in a range of 150 -120°C km -1 and the exhumation rates between 3.2 and 1.6 mm yr -1 (~2.4 mm yr -1 on average). In general, the modelled results of apatite and zircon fission-track both suggest that the geothermal gradients can be over ~100°C km -1 consistent with the direct estimates derived from the cooling/ exhumation histories.
Although such a high paleogeothermal gradient is rare, we argue that the intricate tectonic evolution of northeast Taiwan may have perturbed the regional thermal structure as the area has undergone a transition from the arc-continent collision to the subduction and extension since the past 1.5 Ma. Previous seismologic and geophysical researches have proposed the potential of having the high geothermal gradient in the Tananao Complex presently. For instance, the seismicity around the Hoping area is mostly distributed on the lithologic formation boundary and its offshore area and is related to the northward-subduction of the Philippine Sea plate. An aseismic zone is observed under the rigid Tachoshui granitic gneiss within the depth of ca. 15 km. This signature is compared by Lin (2000) with the silent zone located in the eastern Central Range where the highest exhumation rate and geothermal gradient are found in Taiwan. Hsieh et al. (2014) have reported the depth of Curies point in the northern Tananao Complex, located at ca. 10 -11 km deep with an av-erage geothermal gradient of 58 -53°C km -1 derived from magnetic anomaly data. Chi and Reed (2008) have also suggested that the high heat flow near the surface observed in the Tananao Complex could be the result from the higher thermal conductivity of the exhumed basement rocks.
However, the drill hole HCBH01 on the valley bottom of the Hoping River has only recorded core temperature to be around 30°C at the depth of 600 m, suggesting a geothermal gradient of ca. 25°C km -1 as the surface temperature is 15°C. It is obviously much lower than the estimated paleogeothermal gradients for the earlier periods. One of Fig. 5. Modelled cooling paths and estimated paleogeothermal gradients of the Tachoushui granitic gneiss body in the Hoping area, using the HeFTy program (Ketcham 2005). The two black boxes represented additional thermal constraints during the exhumation: one is the peak metamorphic temperature (P) based on RSCM ) at 6.5 -6.0 Ma; the other is biotite Ar-Ar date (B) proposed by Wang et al. (1998). The weighted mean time-temperature paths are indistinguishable between the two terranes. The closure temperature of fission-track system for apatite (green box) and zircon (orange box) are constrained with its corresponding ages. Cooling rates and paleogeothermal gradients are shown an average value in different time span and gradually increasing to the present. the possible factors for the heat loss in this area could be the transport and circulation of liquid within the aquifers. The adjacent Hoping Basin has developed since ca. 1 Ma (Lallemand et al. 1997) and brought up the possibility of seawater intrusion through the fractures in the basement rocks. Groundwater from the mountain ranges in the vicinity would be another important source of the cooling media since the drill core from the borehole has reported the occurrence of dislocated fractures, suggesting possible circulation of groundwater within the permeable units at the shallower crust.

CONCLUSION
The AERs derived from apatite and zircon fission-track ages of the Tachoshui granitic gneiss body shed light on the cooling and exhumation histories of the post-collisional orogen in northeastern Taiwan. On the basis of the previously published biotite 40 Ar/ 39 Ar ages and a geothermal gradient of 30°C km -1 is assumed, the exhumation rates can be estimated as 3.4 -2.5 mm yr -1 since 4.0 -3.1 Ma under the compressive stress conditions induced from the arc-continent collision. During a period of ~1.5 -0.2 Ma, the exhumation rates remained steady to be ca. 2.3 -1.6 mm yr -1 based on the similar slopes of AER for apatite and zircon fission-track datasets. Therefore, the declined in exhumation rates happened sometime in between 3.1 -1.5 Ma, suggesting that the extension above the Ryukyu subduction system might have started to influence the Hoping area. After ~0.2 Ma, the exhumation rates have increased to ca. 5 mm yr -1 at the higher elevation of the mountain region but undisturbed at the valley bottom, evidenced by the projection of the apatite AER to a zero age. The offsets in ages for the apatite and zircon fission-track data is supportive of the proposed normal faulting, resulting in the difference in exhumation rates between two terranes after ~0.2 Ma.
The paleogeothermal gradients at specific periods of time and temperature ranges were estimated as dividing the cooling rates by the exhumation rates. The modelled cooling paths of the granitic gneiss show an accelerated cooling history for the recent 1.5 Ma, while the exhumation rates have kept constant. Therefore, the estimated paleogeothermal gradients in the Hoping area might have temporally increased from ca. 75 -150°C km -1 during 1.5 -0.2 Ma. The calculated geothermal gradients are in good agreement with the analytical solutions and can reach as high as 100 -150°C km -1 , suggesting the occurrence of high geothermal gradients in the Hoping area. The existence of high geothermal gradients agrees with the aseismic pattern observed from the distribution of earthquakes beneath the granitic gneiss body where the stress is spatially transiting from collision to the extension.