Climate Changes of the Recent Past in the South American Continent: Inferences Based on Analysis of Borehole Temperature Profiles

1. Introduction

A detailed understanding of the nature of past climate changes is important in assessment of the effects of global warming trends identified in meteorological records (Hansen and Lebedeff, 1987). Nevertheless, there are large uncertainties in the reconstruction of the climate history of times prior to the period of instrumental records, there being considerable difficulties in experimental determination of past climate changes. In this context geothermal methods, based on results of temperature logs in boreholes stands out as one of the few methods that allow direct measurement of thermal signals in the subsurface induced by climate changes of the past.

Geothermal methods have been employed during the last few decades in extracting information on climate changes of the recent past for several regions of the northern hemisphere (e.g. Cermak, 1971; Lachenbruch et al, 1982; Beltrami et al, 1992; Bodri and Cermak, 1995; Duchkov and Sokolova, 1998; Harris and Chapman, 2001; Majorowicz and Safanda, 2001; Golovanova et al, 2001; Beltrami and Burlon, 2004). Nevertheless, very few attempts have been made in using geothermal data for examining climate variations in low latitudes of the southern hemisphere. Among the studies in this category are the works carried out in Australia (Cull, 1979, 1980; Torok and Nicholls, 1996; Taniguchi et al, 1999a, 1999b; Appleyard, 2005), Brazil (Hamza, 1991; Hamza,1998; Hamza et al, 1991; Cavalcanti and Hamza, 2001; Cerrone and Hamza, 2003) and South Africa (Tyson et al, 1998; Jones et al, 1999). Much of the work carried out in Brazil remain as publications of limited access, such as internal reports (Hamza et al, 1978; Eston et al, 1982; Hamza et al, 1987), academic theses (Vitorello, 1978; Araújo, 1978; Santos, 1986; Ribeiro, 1988; Del Rey, 1989; Cavalcanti, 2003) and meeting proceedings (Ribeiro, 1991; Souza et al, 1991; Hamza and Cavalcanti, 2001; Cerrone and Hamza, 2003; and Conceição and Hamza, 2006). In the present work we provide a synthesis of these earlier works, with emphasis on progress obtained during the period 2006 – 2011.

As prelude to the discussion of results presented in this work we provide first a brief description of the sources and characteristics of the geothermal data employed and the criteria used for data selection. Details of the methods used for extracting information on past climate are set out in the next section. The results obtained in model simulations of temperature-depth profiles are classified into groups, representative of the major geographic zones. Finally, the climate history of South America, deduced from geothermal data, is compared with results of geothermal climate reconstructions from other continental areas.

2. Characteristics of the data base

According to the recent compilations carried out by the National Observatory (Observatório Nacional – ON/MCT) in Brazil geothermal measurements have been carried out in over 5000 localities in South America (Hamza et al, 2010; Vieira and Hamza, 2010). Most of the earlier data were acquired as parts of basic research projects for heat flow determinations and also as parts of applied research projects for oil exploration and geothermal energy assessments. The focus of data acquisition in the earlier works has been on determining temperature gradients in the deeper parts of the boreholes.

The characteristics of these data sets are variable, depending on the methods used for primary data acquisition. Of these, only the ones acquired using the so-called conventional (CVL) method provide direct information on the vertical distribution of temperatures at shallow depths and hence are potentially suitable for climate related investigations. The conventional method has been employed for geothermal studies in 134 localities, which is slightly more than 10% of the overall data set. It includes mainly temperature logs in bore holes and wells and thermal property measurements on samples representative of local geologic formations. In some cases estimation of radiogenic heat production was also carried out. The details of the experimental techniques employed for temperature and thermal conductivity measurements have been discussed in academic theses (e.g. Vitorello, 1978; Araújo, 1978; Del Rey, 1989) and in publications dealing with heat flow measurements (Hamza et al, 1987; Hamza and Muñoz, 1996; Gomes and Hamza, 2005). A direct evaluation of the quality of data acquired in the earlier works is a difficult task since the experimental techniques used for temperature and thermal conductivity measurements have undergone substantial changes over the last few decades.

The sources of conventional data sets employed in the present work may be considered as falling into five main groups:

• New data acquired in ten different localities in the Amazon region, during the period of 2006 to 2008 (Hamza, 2006;, Hamza 2007);

• Results of recent geothermal measurements in the cordilleran region of Colombia (Hamza et al, 2009; Alfaro et al, 2009).

• Data acquired during the period of 2000 - 2005, as part of geothermal projects for mapping heat flow variations in the coastal area of southeast Brazil (Gomes, 2003; Gomes and Hamza, 2005; Hamza et al, 2005). Some of these have been employed in studies of climate change the state of Rio de Janeiro (Cerrone and Hamza, 2003; Hamza et al, 2003);

• Data acquired during the decade of 1980, mainly in the state of São Paulo, as parts of hydrocarbon and geothermal energy exploration programs. Some of these data are reported as parts of academic theses of the 1980s (Santos, 1986; Ribeiro, 1988; Del Rey, 1989) while some are parts of related publications (Santos et al, 1986; Del Rey and Hamza, 1989; Hamza et al,1987);

• Data acquired during the decade of 1970 in the southern and eastern parts of Brazil. Most of the results have been published in the Brazilian Geothermal Data Collection – Volume 1 (Hamza et al, 1978). Some are part of academic theses of the late 1970s (Vitorello, 1978 and Araújo, 1978) and related publications (Vitorello et al, 1978; Hamza, 1982). Data acquired for four boreholes in the cordilleran region of Peru is reported in the IHFC data collection by Huang and Pollack (1998). The geographic distribution of the overall data sets is illustrated in the map of Figure (1).

Figure 1.

Localities of geothermal measurements in South America (Vieira and Hamza, 2010). Letters in the legend refer to Koppen climate classification: A – Tropical; B – Dry; C – Temperate; D – Continental; E - Polar (Koppen, 1936; Pidwirny, 2006).

2.1. Selection criteria for climate studies

The characteristics of the conventional (CVL) data set were examined carefully to screen out records with indications of possible perturbations arising from non-climatic effects. Also, it was necessary to eliminate those which do not provide fairly reliable determinations of both the steady and the transient components of the subsurface thermal field. In an attempt to guarantee the reliability of the data set the following quality assurance conditions were imposed:

1. The depth of borehole is sufficiently large that the lower section of the thermal profile allows a reliable determination of the geothermal gradient, free of the effects of recent climate changes. Order of magnitude calculations indicate that surface temperature changes of the last centuries would penetrate to depths of nearly 150meters, in a medium with a thermal diffusivity of 10^-6 m²/s. Thus boreholes of at least 200 m deep are necessary for a reliable determination of the local geothermal gradient. The choice of this depth limit is rather arbitrary since the possibility that low amplitude climate signals of earlier periods are present at larger depths of up to several hundreds of meters cannot entirely be ruled out. However it is a reasonable compromise for examining subsurface thermal effects of ground surface temperature (GST) variations of the last few centuries;

2. The temperature-depth profile is free from the presence of any significant non-linear features in the bottom parts of the borehole, usually indicative of advection heat transfer by fluid movements, either in the surrounding formation or in the borehole itself;

3. The time elapsed between cessation of drilling and measurements in boreholes is at least an order of magnitude large compared to the duration of drilling, minimizing thereby the influence of eventual thermal perturbations generated during the drilling activity;

4. The lithologic sequences encountered in the borehole have relatively uniform thermal properties and are of sufficiently large thickness that the gradient changes related to variations in thermal properties does not lead to systematic errors in the procedure employed for extracting the climate related signal; and

5. The elevation changes at the site and in the vicinity of the borehole are relatively small so that the topographic perturbation of the subsurface temperature field at shallow depths is not significant.

The sites of these selected boreholes are distributed in the eastern parts of Brazil (in the states of Santa Catarina, Paraná, São Paulo, Minas Gerais, Rio de Janeiro and Bahia), in the Amazon region, central cordillera in Colombia and eastern cordillera in Peru. The majority of the selected temperature logs are from boreholes with depths greater than 200 meters. Some log data for depths less than 200 m were also considered, as these are found to provide complementary information on subsurface temperature fields at shallow depths which may be compared with those encountered in areas where deeper boreholes are situated. Most of the data sets acquired during the decades of 1970 and 1980 have temperature measurements at depth intervals of ten to twenty meters. In more recent logs measurements have been made at intervals of two meters. Such recent logs may be considered as capable of providing more robust estimates of the background temperature gradients. Data from boreholes with depths shallower than 150 meters were not considered in the present work in view of the potential uncertainties in the determination of the local undisturbed gradient and consequent difficulties in extraction of the climate signal. On the other hand, information on climate changes of the recent past available in such logs may be used in obtaining qualitative estimates of GST changes. Data acquired at shallow depths of less than 20 meters were excluded from analyses for climate changes, avoiding thereby eventual perturbing effects of diurnal and seasonal variations in the reconstruction of surface temperature history.

In an earlier study Hamza et al (2007) reported results of geothermal measurements for a number of sites in the Brazilian territory. In the present we have included results of additional studies carried out during the period of 2006 to 2011. The data sets have been classified into groups, designated as subtropical highlands, subtropical humid zones (of the interior and of the coastal areas), tropical Amazon region, semi-arid zones and cordilleran regions in western parts of the continent. The area extents of these geographic sectors are in large part similar to the prevailing climate zones indicated in Figure (1). Typical examples of temperature profiles encountered in boreholes in Brazil, Peru and Colombia are illustrated in Figure (2). In this figure some of the temperature-depth profiles have been shifted laterally to convenient positions along the temperature axis, to avoid overlap and to allow for easy visualization. Consequently, the temperature axis in Figure (2) displays only relative values.

Figure 2.

Vertical distributions of temperatures in selected localities in Brazil, Colombia and Peru. The climate changes of the recent past are often considered responsible for the non-linear features in temperature logs.

The vertical distributions of temperatures in boreholes, some examples of which are illustrated in Figure (2) reveal several remarkable features. For example, the temperature depth profiles at shallow depths are consistently concave towards the depth axis, which is indicative of surface warming events of relatively recent times. The widespread occurrence of such temperature-depth profiles in almost all of the major geographic zones, irrespective of the local geological complexities and changes in soil type, is considered a clear indication that the observed features are generated by surface warming events of large spatial dimensions. Temperature profiles that are convex towards the depth axis, and hence characteristic of cooling events, were not encountered.

3. Methods employed in data analysis

Three different methods have been employed in data analysis: forward models, inversion methods and signal backstripping methods. Brief descriptions of these methods are provided in the following sections.

3.1. Forward models

The basics of forward modeling approach has been discussed extensively in the literature (Birch, 1948; Cermak, 1971; Vasseur et al, 1983; Lachenbruch et al, 1986, 1988). For the case where surface temperature variation can be represented by a power law relation, analytic solutions are readily available. Thus, for a linear (or ramp type) change in surface temperature the relation between the amplitude of the climate signal (ΔT) and the time elapsed (t) at any depth (z) is given by the relation (Carslaw and Jaeger, 1959):

$$\Delta T(z)=4\Delta T{i}^{2}erfc\left(z/\sqrt{4\kappa t}\right)$$E1

where i² erfc is the second integral of the complementary error function and κ the thermal diffusivity of the medium. The best fitting ramp function is obtained by inverting the above relation using iterative procedures such as linearized Newton’s method. The iterative procedure for this model, referred to as Ramp Inversion (Chisholm and Chapman, 1992; Golovanova et al, 2001; Roy et al, 2002), allows simultaneous determination of the magnitude of surface temperature change and the period of its occurrence. The main limitation of the forward model approach is that it resolves mainly the first-order features in the GST history. This is a consequence of the implicit assumption that the bottom parts of the log data, employed in determination of background temperature gradients, are free of transient perturbations.

The basic steps in forward model approach include identification and separation of the steady and transient components present in temperature profiles at shallow depths. The steady state component is determined by the flow of heat from the interior of the Earth while the transient component is induced by downward propagation of a climate related thermal signal induced at the surface. Usually a linear fit to the deeper portion of the log data, where the climate perturbation is practically absent, allows determination of the local temperature gradient. However some care is necessary in selecting the depth interval for determination of the gradient. If the temperature gradient is calculated using a small subset of data from the lowermost part of the borehole its standard deviation (σ_G) is likely to be relatively large, a consequence of the large root mean square (rms) deviation associated with the small number of data points. Progressive inclusion of data from the overlying parts in least square analysis leads to a steady initial decrease in σ_G, as the estimation of gradient becomes more robust. However, as more data are included from shallower depths (where non-linear features are present) this tendency is reversed and σ_G increases. In the present work, the depth corresponding to the minimum value of σ_G is considered as indicative of the top of the unperturbed zone. The background gradient determined for the depth interval below this zone is used for calculating the steady component of the temperature field. Subtracting it from the observed temperatures allows determination of the transient component. A similar procedure has also been employed by Roy et al (2002) in the separation of steady and transient components of temperature profiles in the Indian subcontinent.

3.2. Inverse models

In the inverse problem approach (Tarantola and Valette, 1982) a priori information is explicitly incorporated in constraining the solution. The functional space inversion (FSI) method discussed by Shen and Beck (1991, 1992) and Shen et al (1992) makes use of the non-linear least squares theory in solving the one dimensional heat conduction equation in a layered half space. The algorithm employed finds the model that minimizes the misfit function:

$$S(m)=\frac{1}{2}\left\{\left[{\left(d-d_0\right)}^t{C}_d^{-1}\left(d-d_0\right)\right]+\left[{\left(m-m_0\right)}^t{C}_m^{-1}\left(m-m_0\right)\right]\right\}$$E2

where d and d₀ are respectively the calculated and observed temperatures, m and m₀ the calculated and a priori model parameters and C_d and C_m the covariance matrices of d₀ and m₀. The term C_d indicates the uncertainty in the observed temperature-depth data while C_m indicates uncertainty in the a priori model. The selection of appropriate values of a priori standard deviations for the temperature (σ_d0) and thermal conductivity (σ_k0) data are important in determining the solutions.

The main advantage of FSI formulation is that it does not predetermine the steady state temperature profile. Instead, both the steady state and transient profiles are estimated simultaneously. In addition, it allows consideration of the vertical distribution of thermo-physical properties and their uncertainties as model parameters, allowing thereby determination of a more detailed GST history where it is possible to identify second order features. Also, FSI formulation includes as model parameters all variables that govern the conductive thermal regime (background heat flow density, the GST history as a function of time and thermal properties as functions of depth).

As prelude to the presentation of the results obtained by the functional space inversion method we provide brief descriptions of the steps taken in data processing and analysis. These include specifying the depth intervals and time periods used for inversion, setting a priori standard deviations of temperature and thermal conductivity data sets and measures taken to minimize the undesirable consequences of the null hypothesis in the inversion scheme.

The estimates of the depth at which the thermal regime is supposedly untouched by the GST variations and the time limit beyond which GST variations cannot be recovered from the given temperature log data, were chosen in accordance with the depth extent of the available temperature log data. In particular, the depth estimate is set to be greater than the deepest data point because the calculated data are projected (interpolated) from the finite element solution. As for the time limit there is no harm in setting a value compatible with or greater than the depth extent of the borehole (Shen and Beck, 1992). On the other hand, use of a shorter-than-necessary time span would end up in “telescoping” the GST history. Unless there are independent evidences indicating a rapid return to unperturbed conditions it seems prudent to assume that this return take place gradually. Use of shorter time spans leads to slight reductions in the magnitude and duration of the cooling events of the earlier periods. As pointed out by Shen and Beck (1991; 1992) the results of GST history, determined by functional space inversion, is sensitive to a priori standard deviations of thermal conductivity. The preferred values for the standard deviations are based on considerations of the trade-off between consistency of the solution and data resolution. In the present work we have used 50mK for standard deviation of the temperature data (σ_d0) and 1W/m/K for standard deviation of the thermal conductivity (σ_k0). These values are of the same order of magnitude as those adopted by Safanda and Rajver (2001) and Golovanova et al (2001).

FSI inversion makes use of an a priori null hypothesis in obtaining robust estimates of the prior steady state. However, this built-in feature can potentially lead to undesirable results when inversions are attempted for determining GST history of periods not correlated with the subsurface temperature data. Thus, in the absence of suitable temperature data for shallow depths the inversion scheme generates artificial values for the late part of GST history. It is clear that acquisition of reliable temperature data for shallow depths is important in determining GST history of the last few decades. On the other hand, the results also show that occurrence of artificial cooling trends for the decades prior to 1970 are possible only in cases where temperature measurements are restricted to depths greater than 100m.

3.3. Signal back stripping approach

Both the conventional and the Bayesian inversion methods have inherent difficulties in identification of individual thermal signals originating from climate variations that are episodic. In such cases methods based on signal back stripping approach are more convenient. Following the standard practice we also assume that the residual temperature profile represents is a superposition of individual perturbations. The characteristics of such perturbations vary, mainly as a result of differences in the magnitudes and time periods of individual GST episodes. The essence of the procedure adopted in the present work can be understood by considering magnitudes of temperature perturbations at two conveniently selected depths z₁ and z₂ in the residual profile:

$$d{T}_1=Aerfc\left(z_1/\sqrt{4\kappa t\text{'}}\right)$$E3a

$$d{T}_2=Aerfc\left(z_2/\sqrt{4\kappa t\text{'}}\right)$$E3b

Dividing (3a) by (3b) and designating the ratio dT₁/dT₂ by δ we have:

$$\frac{d{T}_1}{d{T}_2}=\delta =\frac{erfc\left(z_1/\sqrt{4\kappa t\text{'}}\right)}{erfc\left(z_2/\sqrt{4\kappa t\text{'}}\right)}$$E4

Note that equation (4) does not depend on the magnitude of the temperature perturbation but only on the selected value of the time period. Iterative methods may now be employed for determining the appropriate value of δ for the depth interval. The magnitude of this perturbation can be obtained from the relation:

$$A=\frac{d{T}_1}{erfc\left(z_1/\sqrt{4\kappa t\text{'}}\right)}$$E5

In applying this procedure it is necessary to start with results of the lowermost section of the residual temperature profile. If the borehole is sufficiently deep it is fairly reasonable to assume that the residual temperatures of the deeper parts retain the effects of only the earliest perturbation. The magnitude and period of this earliest perturbation can be determined through the use of equations (3) and (5). It also opens up the possibility of removing the effects of this particular perturbation from the original residual temperature profile. The result is a back stripped profile free of the effects of the earliest perturbation. The procedure is repeated successively for identifying and removing the perturbations arising from later climate episodes.

4. Estimates of surface temperature changes

4.1. Forward model results

A summary of the results obtained in fitting forward models to the observational data discussed in this work, is presented in Table (1). It includes magnitudes of the GST change (ΔT) and their duration (t) as well as the values of the root mean square (rms) misfit between the model and the observational data. For reasons of brevity, we present here only the results for the ramp function model. Climate changes inferred on the basis of this model are grouped together for the five major geographic zones: subtropical highlands (elevations >400 meters) subtropical lowlands (elevations <400m), humid regions, semi arid zones, tropical rain forest and cordilleran regions.

The data in table (1) reveal differences in GST between the major geographic zones. For example, the highland areas seem to be characterized by changes in GST in the range of 2 to 3.8^oC, the mean value being 3.3^oC. The lowlands of the interior also seem to have nearly similar GST changes. On the other hand, the humid zones have relatively low GST values, in the range of 0.4 to 1.8^oC, with a mean of 1.1^oC. Low GST values are also found in the semiarid zones in the northeast, in the interval of 1 to 2^oC. Also, the times of duration of climate changes are different for the geographic zones. Thus, the subtropical highlands, subtropical lowlands and humid zones are characterized by climate episodes with durations varying from 40 to 120yrs, the mean value being 80 yrs. The semi arid zones on the other hand have climate variations with relatively larger duration, the mean value being 150yrs.

The vertical distributions of transient components derived from forward model results (which are considered as indicative of temperature signals resulting from climate changes) are illustrated in Figure (3). Note that the magnitudes of the transient components decrease rapidly with depth. At depths greater than 150 meters these values fall below the experimental detection limits of temperature changes in boreholes. The period of GST changes fall in the range of 130 to 260 years. Hence the beginning of the climate change seems to have taken place during the period of approximately 1750 to 1850. A closer examination of the results illustrated in Figure (3) reveals some marked differences in the GST values within individual geographic zones. The primary reason for the occurrence of such intra-zonal variations is unknown at the moment, but it is likely that they are related to microclimatic histories of individual sites.

Also, the depth distribution of the transient components indicates that the magnitude of warming event is relatively smaller (in the range of 1.4 to 2.2⁰C) for the semi-arid zone in north eastern parts of Brazil. The close agreement between the results for the localities in Table (1) is considered as indication that local changes in vegetation cover and soil types have only a minor influence on the surface thermal budget in semi-arid zones. The GST change in this region seems to have had its beginning during the time period of 1850 to 1900, significantly earlier than the corresponding periods for other geographic zones.

Climate Zone	Locality	Coordinates	ΔT; (^oC)	Climate variation		rms; (mK)
Duration	Year of Onset
Subtropical Highlands	Á. Lindóia	22^0 29’/46^0 38’	3.2	105	1877	4.4
Amparo	22° 43’/46° 46’	2.8	115	1852	5.5
Araras	22° 21’/47° 22’	3.4	40	1942	6.4
Cosmópolis	23^0 43’/47^0 12’	3.8	75	1907	5.5
Itapira	22^0 28’/46^0 43’	3.8	40	1942	10.2
Rafard	23° 00’/47° 31’	3.8	60	1922	6.8
Jacarei	23^0 18’/45^0 57’	2.4	105	1880	3.5
São Paulo	23° 34’/46° 44’	2.1	105	1905	7.2
Serra Negra	22° 36’/46° 32’	2.0	105	1877	3.9
P. Caldas	21^0 55’/46^0 25’	2.0	50	1926	13.7
Teresópolis	22° 26’/ 42° 57’	2.0	55	1945	5.5
L. Muller	28º 40´/49º 30´	2.6	80	1895	10.0
Subtropical Lowlands	Seropédica	22º 46´/43º 39´	3.8	145	1860	23.1
Miracema	22º 01'/41º 06'	3.6	30	1969	7.3
Campos	21º 46'/41º 17'	3.0	20	1980	6.7
Itapemirim	19° 06’/41^0 04’	3.6	65	1910	14.1
S. Sebastião	23° 48’/45° 25’	2.7	110	1880	7.8
Subtropical Humid	Itu	23^0 15’/47^0 19’	1.2	60	1922	3.2
Jundiaí	23^0 10’/46^0 52’	1.2	50	1932	7.2
Papanduva	26^0 23’/50^0 08’	1.2	90	1885	15.0
Cach. Sul	30º 00´/52º 55´	1.8	70	1905	12.0
Maricá	22° 54’/42° 45’	0.4	15	1985	2.9
Rio Bonito	21° 25’/42° 12’	1.8	40	1960	7.7
Semi-Arid	Arraial	12° 32’/42° 50’	1.4	85	1890	8.4
Caraiba	09° 28’/39° 50’	1.4	110	1865	9.3
Jacobina	11° 11’/40° 31’	1.9	150	1828	8.7
Poço Fora	09^0 41’/39^0 51’	1.8	105	1870	4.9
Tropical Rain Forest	Manaus	02º 57´/60º 01´	1.3	40	1960	15.0
Belém	01º 27´/48º 27´	2.2	55	1953	9.4
Salinópolis	00º 38´/47º 20´	2.1	65	1943	8.6
Dom Eliseu	04º 17´/47º 34´	1.9	60	1948	9.2
Cordilleran Regions	Pen 742	81º 10´/04º 17´	4.0	200	1779	7.4
Lobitos	81º 16´/04º 27´	2.5	200	1779	8.3
Lomitos	70º 39´/17º 16´	2.0	200	1779	7.2
Bogotá	74º 05´/04º 38´	2.8	150	1860	9.5

Table 1.

Results of Ramp Inversions of GST changes, for selected localities in South America. ΔT is the magnitude of ramp change, and rms the root mean square misfit.

Figure 3.

Magnitudes of climate changes and its vertical distributions, deduced from forward model fits to borehole temperature profiles at selected localities in South America.

A careful examination of figure (3) reveals that some of the transient temperature profiles do have small negative values for depth intervals corresponding to the lower parts of the transient sections. The negative residuals are usually considered as arising from cooling events prior to recent warming episodes. However, the magnitudes of such events appear as subdued features, because of the implicit assumption in forward model approach that the bottom parts of the borehole are free of transient perturbations. Results of numerical simulations indicate that occurrence of negative residuals is quite sensitive to subtle changes in the value adopted for the background temperature gradient. This is the main limitation of the forward models, which resolves mainly for the first-order features in the GST history. We conclude by noting that residual temperature profiles similar to those obtained in the present work were also reported by Golovanova et al (2001) for the Urals region and Roy et al (2002) for the Indian Peninsula. A complementary analysis of this problem is provided in section (4.2) below, where we discuss the residual temperature profiles derived using the inversion method.

4.2. Results of functional space inversion

The method of Functional space inversion (FSI) was employed in analysis of temperature-depth profiles for 14 localities distributed over three main geographic zones of Brazil. The criteria used in the selection of profiles included availability of both thermal property data of subsurface layers and supplementary information on the history of changes in the vegetation cover. In discussing the results it is important to point out that the FSI method provides more detailed information on the GST history than that provided by the forward model approach. However, in comparing the GST histories of several localities it is convenient to work with deviations from the site specific mean rather than the absolute value. In the present work, GST deviations are calculated by subtracting the model results from the a posteriori estimate of the site specific mean.

The results of FSI method, illustrated in figure (4), reveal several characteristic features in the GST history of the study area. Foremost among these are the indications that surface temperatures have increased by as much as 1 to 4^oC, over the last century. This observation is in reasonable agreement with the results of the ramp function model discussed in the previous section. However, in all localities the warming events seem to be preceded by cooling episodes occurring over the time period of approximately 1700 to 1900. The amplitudes of the cooling events are much less, falling generally in the range of 0.5 to 1^oC. For time periods prior to the 17th century the resolving power of the inversion method is poor, a consequence of the limitations in sensitivity and precision of sensors used for temperature measurements in boreholes. Hence FSI model calculations for periods prior to 1700 may not necessarily be representative of true climate history.

Even though the data set is poor the results of Figure (4) seem to indicate that the warming trends are less pronounced in semi-arid regions relative to the tropical humid regions. The results also reveal a small time shift in the occurrence of climate warming in highlands regions compared with that for semiarid regions. For example, the warming event in the highlands region had its beginning during the time period of 1850 to 1900 while that for semi-arid regions seems to have had its beginning during the period of 1670 to 1860. Also the durations of both warming and cooling episodes appear to be relatively smaller in highland regions when compared with those for the semi-arid regions. Such differences have important implications for understanding evolutionary trends in climate history of eastern Brazil.

The vertical distributions of transient components derived from FSI inversions are similar to those found for the forward model results, illustrated in Figure (3). The magnitudes of the transient components decrease rapidly with depth. However, the residual temperatures become significantly negative at depth intervals corresponding to the lower parts of the transient sections. Such negative residuals are indicative of the occurrence of cooling events prior to recent warming episodes. We conclude that the vertical distributions of residual temperatures in FSI method are better representations of subsurface transient thermal regimes than those derived from the forward model approach.

Another conspicuous feature of the results by the FSI method in Figure (4) is the presence of brief cooling episodes for the recent decades, since 1970. There is a possibility that short period cooling episodes are spurious, a consequence of the null hypothesis employed in the computational process of the FSI method (Shen and Beck, 1992). The potential undesirable effects of the null hypothesis may be minimized by limiting the GST history to periods compatible with reliable subsurface temperature data. In the present work temperature data are available for depths less than 50 meters in all of the data sets employed in inversion schemes. Consequently the cooling trends observed in our results, for the decades prior to 1970, cannot be attributed to the spurious effects of the null hypothesis. Examples of cooling trends similar to those found in the present work can also be seen in the results reported for several localities in Europe and North America (Bodri and Cermak, 1995; Rajver et al, 1998; Stulc et al, 1998; Golovanova et al, 2001).

Figure 4.

History of ground surface temperature variations associated with climate changes of the recent past, deduced on the basis of functional space inversion models fits to borehole temperature profiles, for selected localities in South America (Adapted with modifications from Hamza et al, 2007).

The uncertainties in the estimates of GST in the inversion method can somewhat be improved by carrying out simultaneous inversion of temperature profiles of several sites in the same geographical province. In the present work we have carried out simultaneous inversions only for repeat measurements at Seropédica, in the state of Rio de Janeiro. In view of notable differences in the local soil conditions and vegetation cover, no attempt has been made for carrying out simultaneous inversions of temperature profiles from different locations discussed in the present work. Also, most of these sites are separated by large distances and fall within areas with distinctly different microclimate conditions and geographic characteristics.

A summary of the results of inverse modelling is presented in Table (2). The summary includes maximum and minimum values of ground surface temperatures and their respective times of occurrence. Also given in this table are the difference between the maximum and minimum values of GST, the time elapsed between the maximum and minimum and aposteriori estimates of undisturbed GST. The values of the differences in magnitudes are in the range of 1.3 to 3⁰C, which is slightly lower than the range indicated by the results of the forward model approach. On the other hand, the values of the differences in the time periods are in the range of 62 to 147 years, which is comparable to the ranges indicated by the results of the forward model approach.

Locality	Maximum		Minimum		Difference		T0 (°C)	Log Year
°C	Age	°C	Age	Mag.	Yr
Amparo	21,1	1949	18.1	1802	3.0	147	18.9	1982
Itu	20.7	1982	19.4	1920	1.3	62	19.9	1982
Araras	22.8	1982	19.8	1910	3.0	72	20.9	1982
Campos	25.7	1999	22.7	1936	3.0	64	23.9	2000
C. Itapemirim	23.7	1971	20.0	1883	3.7	88	21.2	1976
Arraial	28.8	1967	27.2	1863	1.6	104	27.7	1976
Caraíbas	30.1	1946	28.6	1812	1.6	134	29.0	1976
Poço de Fora	30.4	1976	28.2	1822	2.1	154	28.7	1976

Table 2.

Results of functional space inversions for selected data sets. Maximum and minimum values of GST variations and their respective years of occurrences are given in columns 2 - 5. Also given are model results for GST deviations and the respective periods. T₀ is the aposteriori steady state temperature.

4.3. Results of signal back-stripping approach

Techniques of signal back-stripping were used in analysis of temperature log data from 15 localities in south east Brazil. In the present work we limit the discussion to the results obtained for the locality of Seropédica (Rio de Janeiro). The vertical distributions of temperatures obtained during the various stages of the back stripping process are illustrated in the set of upper and lower panels of Figure (5).

The segment on the left side of the upper panel in this figure refers to results of the first stage of the back stripping process. Here the blue curve indicates the initial reduced temperatures and the red curve the first signal extracted by the back stripping method. The extracted signal points to an episode of climate change that took place 115 years back in time and had a magnitude of 2.9⁰C. Similarly, the segments on the right side of the upper panel and that on the left side of the lower panel of this figure illustrate, respectively, the results of the second and third stages of the back stripping process. According to the results obtained the signal extracted in the second stage points to an episode of climate change that took place 93 years back in time and had a magnitude of 0.2⁰C. On the other hand, the signal extracted in the third stage points to an episode of climate change that took place 21 years back in time and had a magnitude of 0.8⁰C. The segment on the right side of the lower panel illustrates the residual temperatures, after extracting the above mentioned three climate signals. The back stripping process is terminated at this third stage since the magnitudes of the residual temperatures are less than 0.05⁰C, below the sensitivity limit of experimental system used in data acquisition. The overall history of climate change derived from the back stripping method is illustrated in Figure (6).

Figure 5.

Results of signal back stripping method employed in disentangling the climate history of Seropédica (Rio de Janeiro).

Figure 6.

Climate history of Seropedica (Rio de Janeiro) derived from results of the signal back stripping method.

5. Discussion

5.1. Climate changes in South America

The results obtained in the present work have contributed to substantial improvements in data base for climate change in South America. The area covered in the data base includes several geographic zones in Brazil, Colombia and Peru. Classical inverse models were employed in the analysis of temperature logs from over 30 localities and, in addition, Bayesian inverse modelling was carried out for data from 20 selected sites. The model results have allowed determination of the magnitude as well as the duration of ground surface temperature changes in the major geographic zones of South America. The map of recent climate changes derived from geothermal data is presented in Figure (7).

According to the results obtained the magnitude of GST changes are in the range of 2 to 3.5^oC but have had their beginning during the early decades of the 20th century. Nearly similar trends are seen in temperature-depth profiles of bore holes in tropical as well as subtropical zones of the interior and coastal areas. The data from semi arid zones also indicate occurrence of surface warming events but the magnitudes are in the range of 1.4 to 2.2^oC while the duration of the warming event is larger, extending back into the last decades of the 19th century. The magnitudes of GST variations are relatively large in localities which have undergone recent changes in vegetation cover. Also there are indications that GST changes are practically insignificant in areas of tropical rain forest.

5.2. Global climate changes

The improved data base on GST changes in the South American continent has contributed to a better understanding of climate changes in the southern hemisphere. In this context, it is convenient to examine also its implications for global variations. It is clear that improved assessments of global changes require an integrated analysis of both the geothermal data sets as well as the meteorological data sets (such as HadCRUT3) reported for oceanic and continental regions (Brohan et al, 2006; Jones and Moberg, 2003; Rayner et al, 2003). Major obstacles in such a venture springs from the wide disparities in data density, and it is

Figure 7.

Map of recent changes in ground surface temperatures derived from geothermal data sets, and supplemented with ground data for selected localities.

convenient to adopt procedures that minimize problems arising from the non-homogeneous distribution. In the present case, the surface area of the globe is divided into a regular grid system composed of 5 x 5 degree grid cells and average values of surface temperature variation in the grid elements calculated. Experimental data are available for almost all of the grid elements covering the continental regions and a significant part of the grid system for the oceanic regions. Following the common practice employed in deriving maps of global representations interpolated values were used for grid elements without data.

Techniques of spherical harmonic representation (Hamza et al, 2008) were employed in analysis of global variations. The harmonic representation of surface temperature variations (q) may be represented as:

$$q\left(\theta ,\varphi \right)={\displaystyle \sum _{n=0}^N}{\displaystyle \sum _{m=0}^n\left[A_{nm}\cos \left(m\varphi \right){P^{\prime }}_{nm}\left(\cos \theta \right)+B_{nm}\sin \left(m\varphi \right){P^{\prime }}_{nm}\left(\cos \theta \right)\right]}$$E6

where is the longitude = 90 – ψ, is the colatitude, P´nm(cos ) is the associated Legendre function that is fully normalized and A_nm and B_nm the coefficients of the harmonic expansion. The expression for evaluation of P´_nm is:

$$P_{nm}^{´}=P_{nm}\sqrt{K_n^m}$$E7

where P_nm is the associated Legendre function given by:

$$P_{nm}(\cos \theta )=\frac{se{n}^m\theta }{2^n}{\displaystyle \sum _{t=0}^{Int\left(\frac{n-m}{2}\right)}\frac{{(-1)}^t(2n-2t)!}{t!(n-t)!(n-m-2t)!}{\cos }^{(n-m-2t)}\theta }$$E8

and

$$K_n^m=\frac{H(2n+1)(n-m)!}{(n+m)!}\stackrel{}{},\stackrel{}{}\{\begin{array}{c}if\stackrel{}{}m=0\Rightarrow H=0\\ if\stackrel{}{}m\ne 0\Rightarrow H=2\end{array}$$E9

In equation (9) $$Int\left(n-m/2\right)$$is the largest integer that is lower than$$\left(n-m\right)/2$$.

Full normalization of associated Legendre functions (P_nm) requires that the following equations be satisfied:

$${\displaystyle \underset{0}{\overset{2\pi }{\int }}}{\displaystyle \underset{0}{\overset{\pi }{\int }}{\left[{P^{\prime }}_{nm}(\cos \theta )sen(m\varphi )\right]}^{2}sen\theta \stackrel{}{}d\theta \stackrel{}{}d\varphi }=4\pi$$E10a

$${\displaystyle \underset{0}{\overset{2\pi }{\int }}}{\displaystyle \underset{0}{\overset{\pi }{\int }}{\left[{P^{\prime }}_{nm}(\cos \theta )\cos (m\varphi )\right]}^{2}sen\theta \stackrel{}{}d\theta \stackrel{}{}d\varphi }=4\pi$$E10b

The coefficients A_nm and B_nm are evaluated by fitting the harmonic expansion to the set of experimental data, which are the values of surface temperature changes (q) and their respective geographic coordinates ( and ).

The results obtained on the basis of harmonic representation of the mixed data set are illustrated in the global map of surface temperature changes, in Figure (8). According to features discernible in this figure the thermal signals of climate change in continental areas of North America, Europe, Asia, West Africa and Eastern South America have magnitudes in excess of one degree centigrade, while those in most of the oceanic regions is less than one degree centigrade. On the other hand most of the oceanic regions are characterized by surface temperature variations of less than 0.8⁰C.

6. Conclusions

Much of the information on surface temperature changes of the recent past in South America are derived from historical records of local meteorological stations. Some of the

Figure 8.

Spherical harmonic representation of surface temperature variations associated with climate changes of the recent past. Note the systematic differences in magnitudes of climate change between continental and oceanic areas.

main problems in the use of such records include non-uniformity in the accuracy and precision of sensors employed in data acquisition, gaps in the time series of observations and influence of urban warming effects (most of the meteorological observatories are located in or close to towns or cities). The results reported in the present work, based on geothermal methods, are in general free of such problems.

An important result emerging from studies based on geothermal methods is that the climate was relatively cooler during the 17th and 18th centuries. The climate histories, deduced from geothermal data, are found to be consistent with the results of available meteorological records in Brazil, Colombia and Peru. Comparative studies indicate that the magnitudes and duration of recent climate changes in South America are similar to those found in other continental areas such as North America, Asia and Europe.

Analysis of surface temperature changes on global scale reveal that climate warming effects in the northern hemisphere are more pronounced because of the relatively large areas of continental blocks in these regions. On the other hand, the climate warming effects are relatively subdued in the southern hemisphere, where relatively larger proportions of oceans dominate. It is possible that systematic difference in magnitudes of global warming between continental and oceanic regions is a consequence of large scale mixing of near surface waters with those of deep ocean circulation systems. This mixing process is considered to be responsible for pronounced mitigation of the effects of global warming in the southern hemisphere.

Acknowledgments

This work was carried out as part of a research project with funding from Conselho Nacional de Desenvolvimento Científico – CNPq (Project No. 301865/2008-6; Produtividade de Pesquisa - PQ). We thank the topic editor and reviewers for critical comments and suggestions. Dr. Andres Papa of the Geophysics Department (ON-MCT, Rio de Janeiro) provided institutional support. The source code for functional space inversion program was provided by Dr. Paul Shen, University of Western Ontario, Canada