Shallow Crustal Thermal Structures of Central Taiwan Foothills Region

Crustal thermal structures are closely related to metamorphism, rock rheology, exhumation processes, hydrocarbon maturation levels, frictional faulting and other processes. Drilling is the most direct way to access the temperature fields in the shallow crust. However, a regional drilling program for geological investigation is usually very expensive. Recently, a large-scale in-situ investigation program in the Western Foothills of Central Taiwan was carried out, providing a rare opportunity to conduct heat flow measurements in this region where there are debates as to whether previous measured heat flows are representative of the thermal state in this region. We successfully collected 28 geothermal gradients from these wells and converted them into heat flows. The new heat flow dataset is consistent with previous heat flows, which shows that the thermal structures of Central Taiwan are different from that of other subduction accretionary prisms. We then combine all the available heat flow information to analyze the frictional parameters of the Chelungpu fault zone that ruptured during the 1999, Chi-Chi, Taiwan, earthquake. The heat flow dataset gave consistent results compared with the frictional parameters derived from another independent study that used cores recovered from the Chelungpu fault zone at depth. This study also shows that it is suitable for using heat-flow data obtained from shallow subsurface to constrain thrusting faulting parameters, similar to what had been done for the strike-slip San Andreas Fault in California. Additional fieldworks are planned to study heat flows in other mountainous regions of Taiwan for more advanced geodynamic modeling efforts.


INTRODUCTION
The thermal structure of the crust may affect geological and geophysical processes, e.g., metamorphism, changes in crustal rheology, seismogenic process, seismic attenuation and the maturity of hydrocarbons, to name a few. Previously, a heat flow dataset (Lee and Cheng 1986) was used to help formulate the critical taper wedge model using Taiwan as an example (Barr and Dahlen 1989). Such study has dramatic impacts on how scientists look at the thermal structures of mountain building processes in an accretionary prism. However, there are still debates on whether the heat flow dataset used in that study is representative of the thermal structures of Taiwan. As a result, it is advantageous if we can collect additional heat flow data and compare the new results with previously published datasets. In addition, after two decades, it is necessary to update the onshore heat flow dataset for Taiwan.
One particular controversy is that in Central Taiwan the heat flows from previous studies (Lee and Cheng 1986;Barr and Dahlen 1989;Hwang and Wang 1993) show an increase towards the hinterland, while most of the subduction zone accretionary prisms have decreasing heat flows towards the hinterland (e.g., Hyndman et al. 1993). Whether such a heat flow pattern in Taiwan is true is still unknown. Recently, we have started to conduct field experiments to measure shallow crustal heat flows in the foothills of Central Taiwan to address this question. Dozens of wells are drilled each year in this high cost investigation project. Such wells provide a chance for us to measure temperature profiles. As a result, there will be more heat flow data available for researchers to better understand the heat flow patterns.
These additional heat flow data can be used to constrain seismogenic processes. Repeated slips on a large thrust generate frictional heating and can perturb crustal thermal structures (e.g., Lachenbruch and Sass 1980). Previously, there have been attempts to use heat flow dataset to constrain the strength of the San Andreas Fault in California, USA (Brune et al. 1969;Zoback et al. 1987;Scholz 2000). Similarly, we can use heat flow data to constrain the frictional parameters on faults in Taiwan. In 1999, a M w 7.6 earthquake occurred in Chi-Chi, Taiwan. This earthquake inflicted huge damage, but at the same time, provided a rare opportunity for scientists to study the earthquake rupture parameters. A large-scale drilling experiment designed to penetrate the rupture fault plane at depth was conducted (Ma et al. 2006). Using the core recovered from the deep hole, they studied many earthquake rupture processes during the Chi-Chi main shock, including frictional heating-related parameters (Wang et al. 2002;Chen et al. 2004;Yabe et al. 2008). In this study, we apply these derived parameters to heat flow modeling to see if they are consistent with the heat flow pattern measured in this region. The Foothills of Central Taiwan provide a unique opportunity to test these parameters and the heat flows in an independent way.

REGIONAL TECTONICS OF TAIWAN
Taiwan ( Fig. 1) is located along the boundary between the Eurasian plate and Philippine Sea plate, where to the south, the oceanic lithosphere of the South China Sea is subducting eastward beneath the Philippine Sea plate (Bowin et al. 1978;Ho 1986a) with a high rate of 7 to 9 cm yr -1 (e.g., Yu et al. 1997). Near Taiwan, subduction changes into an arc continent collision where the Chinese continental margin enters the subduction zone (Suppe 1981). The initial collision zone is probably near latitude 21°N where the continent-ocean boundary (COB) enters the trench. As  the Luzon arc approached Taiwan and the Chinese passive margin, volcanism started to cease. The forearc basin is becoming consumed until the Luzon volcanic arc is juxtaposed next to the fold and thrust belt of Taiwan onshore in the mature collision zone. Although the timing of the initial collision is still under debate, Teng (1990) proposed that it probably began in the late middle Miocene and is currently propagating southward toward the modern subduction zone offshore Taiwan (Suppe 1984).
Further to the north in the mature collision zone on land Taiwan, the geologic provinces from west to east include: Coastal Plains, Western Foothills, Hsuehshan Range, Central Range and Coastal Range. The Coastal Plains are composed mostly of recent alluvium deposits overlying continental shelf materials of the Chinese passive margin. The Western Foothills are composed of Oligocene to Pleistocene shallow marine materials. The Hsuehshan Range consists of sandstone, argillites and shales. The western side of the Central Range consists of slaty materials while the eastern part is mostly underlain by the pre-Tertiary metamorphic rocks, including gneiss, marbles and schists (e.g., Ho 1986b). Overall the metamorphic grade increases eastward in the Central Range. Some of the metamorphic rocks have been interpreted as the basement rocks of the Eurasian continental crust (Lan et al. 1990). Further east is the Coastal Range, containing rocks of Luzon arc affinity, both volcanic and sedimentary rocks.
In the northern part of Taiwan, the polarity of the subduction flipped and the Philippine Sea plate is subducting northward underneath the Eurasia plate along the Ryukyu trench. The Benioff zone is well-imaged underneath the eastern part of northern Taiwan at least north of 24°N (Kao et al. 1998).
Recently, Chi and Reed (2008) systematically studied the shallow crustal thermal structures from the subduction zone to the collision zone. They found that ( Fig. 2), in the subduction zone, the geothermal gradient ranges mostly within 30 to 80°C km -1 , and decreases toward the hinterland due to intensive dewatering at the toe, sediment blanketing and topographic effects, cooling from the subducting plate and other processes. The geothermal gradient ranges mostly from 30 to 90°C km -1 in the collision zone and increases, instead of decreasing, toward the ridge, possibly caused by exhumation, erosion, topographically-induced ground water circulation, and some upper mantle processes related to collision. Heat flow in the collision zone ranges from 80 to 250 mW m -2 . In this study we focus mainly on the thermal structures of the Western Foothills of Central Taiwan.

FIELD SURVEYS
For this study we used the compact high-resolution temperature logger (CHTL) developed by the Institute of Oceanography, National Taiwan University (Shyu and Chang 2005;Chang and Shyu 2011). Figure 3a shows the instruments including a 24-bit, low noise A/D converter, capable of resolving 0.0001°C in the -1 to 25°C range. The small size of the thermal sensor at the tip of the CHTL makes Fig. 2. Heat flow patterns from the subduction zone to the collision zone in the Taiwan region, modified from Chi and Reed (2008). Large stars indicate Q1 class, or highest confidence in Bottom-Simulating Reflector (BSR) identification. Medium stars indicate Q2 class, or probable BSR identification. Small stars indicate Q3 class, or possible BSR identification. For details on how to use the BSR to derive heat flow, please refer to Chi and Reed (2008). The large circles are from previously published thermal probe measurements (Shyu et al. 1998(Shyu et al. , 2006. The different sizes and shapes of the marks represent different geothermal gradient data sets on land in the collision zone, which are discussed in Lee and Cheng (1986). Note higher geothermal gradients along the trench in the subduction zone and along the topographic high in the collision zone. BHT: bottom-hole temperature; CPC: Chinese Petroleum Corporation. it quickly attain temperature equilibrium with the water in the well. The compact size of the logger makes it possible to collect temperature profiles using only one field worker, although we usually have a field party of 2 or more persons when conducting such measurements. More than 30 days of fieldwork were required for all measurements.
Starting in 2010, the Central Geological Survey and Sinotech Engineering Consultants Inc. of Taiwan conducted a large scale, high cost hydrogeological investigation in the Western Foothills of Central Taiwan to measure the long term ground water level changes and study the hydrogeo-logical properties of rocks in the higher elevation mountains. This provides a rare opportunity for us to use these wells to measure temperature profiles. There are also some other opportunities when the Central Weather Bureau and Water Resource Bureau of Taiwan has wells available for us to measure geothermal profiles. Forty-four temperature profiles were measured from 42 wells. The length of the wells ranged from 100 to 500 m, mostly between 100 to 200 m.
Here we describe the procedure for temperature profile measurement: as a standard procedure, Sinotech engineers will conduct well logging measurements right after a well is drilled reaching the total depth. The well logging includes acoustic and optical borehole televiewer, caliper, electrical logging, sonic logging, fluid temperature and conductivity, flowmeter measurement. All of these are done using commercial logging tools. Before using the CHTL for temperature measurement, we will wait at least two months until the borehole has reached thermal equilibrium. In this study we only used the data collected by CHTL. The flow chart shown in Fig. 3b describes the measurement procedure. In the field, we first measure the groundwater levels in the borehole; then decide the initial logging depth due to the higher equilibrium rate in water than in air. After that we set up the CHTL so it will start recording temperature every 1 sec. We then lower the CHTL into the well and stop at every 5 or 10 meters. The CHTL will stay at that depth for 45 sec to make sure that it is in equilibrium with the ambient water temperature at that particular depth. A similar process is repeated until the CHTL reaches the bottom of the well. We then retrieve it back to the surface and download the data to a PC for further processing (Chi et al. 2011). After that we make a plot of temperature vs. time, which then will be converted into a temperature vs. depth plot, i.e., a temperature profile.

RESULTS OF FIELD SURVEY
To verify the reliability of the instruments, two different instruments were used for in-situ measurements and they both produced consistent results (Chi et al. 2011). We also measured one single well two different times, with the re-sults nearly identical except at very shallow depth. Shallow depth measurements are impacted by the daily temperature fluctuation. The skin depth of the daily temperature fluctuation is very shallow, about 0.17 m based on Turcotte and Schubert (2002), thus should not affect our measurements at depth.
To derive the preferred geothermal gradient for each temperature profile, we first manually discarded the upper part of the temperature profile that shows ground surfacerelated anomaly. We then applied the least-squares method to derive the geothermal gradient. Overall, the lower part of the temperature profiles are usually quite linear (see Fig. 4) if there are no clear ground water perturbations.
A total of 44 temperature profiles were measured successfully in Central Taiwan. Some of the temperature profiles in the ChoShuiHsi alluvial fan gave abnormal geothermal gradients, and thus were discarded due to possible groundwater circulations. A couple of temperature profiles in the foothills region were very close to rivers and gave very low geothermal gradients, which we also discarded. In the end, a total of 28 geothermal gradients were derived. Figure 5 shows the geothermal gradient measurement results in Taiwan. Table 1 lists the data collected in 2011 for this study. The general trend is that the geothermal gradients increase toward the mountainous regions. We then performed a comprehensive literature search to compile the available thermal conductivities of the rocks based on different sedimentary rock formations (Lee and Cheng 1986). The exhumed metamorphic rocks in the mountainous regions also have higher thermal conductivities. Using a Fig. 4. A typical temperature profile measured by the CHTL. Left bracket shows which part of the curve we used to derive the linear geothermal gradient. WB-OWD site located in the Aowanda National Forest Recreation Area. Other information is listed in Table 1. The maps that combine data from this study and from Lee and Cheng (1986). The triangles are from this study, and other symbols are from Lee and Cheng (1986).
(a)  published geological map (Ho 1986a), we can assign the major rock type, and thus the thermal conductivity, to each well. We then multiplied the geothermal gradient with the thermal conductivity to derive the heat flow. We found some very high heat flow regions, mostly close to hot springs or along the proposed Lishan fault, Coastal Range fault, and Chaochou fault. From the foreland to hinterland, we found increased heat flow, consistent with the previous study (Fig. 6) of Lee and Cheng (1986). Again, such pattern contradicts the typical decreasing heat flow towards the arc in subduction zones elsewhere around the world. In other words, the heat flow pattern in Taiwan should not be used to explain the critical wedge model of convergent boundaries.

IMPLICATIONS FOR THE CHELUNGPU FAULT PROPERTIES
Despite the overall trend in heat flow increasing eastward, we also found some small-scale variations. For example, the heat flow actually decreases, instead of increases, across the Chelungpu fault that ruptured in the 1999 Chi-Chi earthquake in the Western Foothills of Central Taiwan. This pattern has been well-documented (Molnar and England 1990) as a result of down-going footwall that cools the over-riding hanging wall.
Here we studied the cooling effect due to the rapid thrusting of the Chelungpu fault based on the model proposed by Molnar and England (1990) and modified by Von Herzen et al. (2001), the expression is as follows: where Q 0 is the reference heat flow in the region, V is the rate of orthogonal thrusting, and τ is the shear stress associated with slip on the fault. The term τV is for frictional heating on the fault. The z f is depth of the fault below the heat flow measurement, δ is the dip of the fault separating the upper and lower plates, κ is the average thermal diffusivity of the crust. We then used the following parameters to calculate the east-west trending profile showing heat flow values across the frontal thrust sheet of the Chelungpu fault: Q 0 = 69 mW m -2 is the average heat flow at the earth surface for the incoming plate (Lee and Cheng 1986); V = 8.7 mm yr -1 × cos 50° which is the average slip rate on the Chelungpu fault times cosine of the angle between relative motion and the strike of the fault plane (Chen et al. 2004). Based on a crustal model derived from a seismic profile (Wang et al. 2002), we used an average dip of 23° to calculate the depth of the fault at different distances from the surface rupture (Johnson et al. 2001). This is mainly because we are only modeling the ramp of the folded thrust sheet so we do not need to consider the flat part of the thrust sheet in the hinterland in our 30-km length model. We used a specific heat of 775 J kg-K -1 , an average thermal conductivity of 2.9 W mK -1 , and a density of 2400 kg m -3 . The diffusivity of 1.42 × 10 -6 m 2 s -1 was derived using the equation which states that diffusivity equals to the thermal conductivity divided by the product of density and specific heat. We used 1 for the dimensionless factor b, as Molnar and England (1990) found this a reasonable number for most of the cases. The results are plotted in Fig. 7, which shows a first-order match to the heat flows in this region. Fig. 6. From the contour map before and after the adding of our 2011 survey, it demonstrates the detailed thermal distribution in the Western Foothills of Central Taiwan, which is consistent with the previous study by Lee and Cheng (1986).  One important objective of this study was to examine the stress and frictional heating properties of the Chelungpu fault using the heat flow dataset. We tested a range of parameters to see which one gives a better fit. The preferred model has a stress gradient of 14 MPa km -1 . This stress gradient is the average of that estimated by Yabe et al. (2008) using the Chelungpu drill cores near the fault trace, which is from 13.7 to 16.3 MPa km -1 at a depth of 1082 m.
Using the published fault parameters, we are able to calculate the theoretical heat flows as a function of the distance from the Chelungpu fault surface rupture. The theoretical heat flows derived from our preferred model are within the 25% error bounds of the data collected in this study. The reasonable fits suggest that these parameters derived from seismic reflection and deep drilling data are consistent with our heat flow dataset, even though these datasets are totally independent from each other.

CONCLUSIONS
A recent large-scale drilling experiment in the Western Foothills of Central Taiwan provided a rare opportunity to study the shallow thermal structures in that region. We successfully used high-resolution temperature loggers to measure 44 temperature profiles and derived the geothermal gradients. Using the published thermal conductivity of different rock types in Taiwan, we derived a dense heat flow dataset in this region.
The new dataset gave similar heat flow patterns as the Lee and Chen's 1986 dataset, suggesting that there are increasing heat flows to the hinterland, contradicting the typical heat flow patterns found in other normal accretionary prisms around the world. However, there are some smallscale variations in this general trend, including a slight decreasing heat flow across from a thrust that ruptured in the 1999, Chi-Chi, Taiwan, earthquake. Such features are welldocumented as a result of the over-riding block cooled by the underthrusting block. Using an analytical solution, we were able to fit the new measured heat flows using fault parameters recently derived from reflection seismic data and drilling data. This study shows that the stress gradient, thus the fault frictional parameters derived from the drilling are consistent with the new heat flow dataset collected from this study.
We are planning to collect additional heat flow data in the mountainous regions of Taiwan and hope the results will help us better understand the regional shallow thermal structures of Taiwan in an arc-continent collision setting.