# Physical analysis and inverse methods applied to archaeological FIRE replications on ATACAMA desert soils, northern Chile

Édouard Canot, Renaud Delannay, Calogero M Santoro

► To cite this version:

Édouard Canot, Renaud Delannay, Calogero M Santoro. Physical analysis and inverse methods applied to archaeological FIRE replications on ATACAMA desert soils, northern Chile. Archaeometry, 2025, 67 (3), pp.625-641. 10.1111/arcm.13030. hal-04717063

HAL Id: hal-04717063

https://hal.science/hal-04717063v1

Submitted on 25 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ORIGINAL ARTICLE

# Physical analysis and inverse methods applied to archaeological FIRE replications on ATACAMA desert soils, northern Chile

Édouard Canot<fnref n="1" /> | Renaud Delannay<fnref n="1" /> | Calogero M. Santoro<fnref n="2" />

**Correspondence**

Édouard Canot, IPR, UMR 6251, CNRS, Université de Rennes, F-35000 Rennes, France.

Email: edouard.canot@univ-rennes1.fr

**Funding information**

Field and lab work in Chile and France were funded by ANID-FONDECYT project 1201786 and FONDECYT 1120454 and 1160744, and by a Coopération Scientifique ECOS-CONICYT (grant number C13H02). Lab experimental research in France was funded by CNRS and OSU-Rennes.

**Abstract**

Physical analysis of *in situ* fire experiments on soils are useful for the estimation of subsurface thermal diffusivity, which is affected by factors such as water, heterogeneity and heating conditions. To address the uncertainties due to these factors, a new data-processing procedure based on inverse methods was developed and experimentally applied to soils from an archaeological site in the Atacama Desert, Chile. By combining experimental data and numerical simulations, we determined the dominant physical processes arising during the heating. The analysis succeeded in defining practical procedures to obtain a more accurate estimation of the diffusivities, thus reducing the above-mentioned uncertainties.

**KEYWORDS**

Atacama Desert northern Chile, *in situ* hearth experiments, inverse problem, late Pleistocene hunter-gatherers, numerical simulations, physical modelling, soil thermophysical properties

INTRODUCTION

Hearth structures are of prime interest in archaeology in order to understand past fire and hearth management, which in turn illuminates the permanence and intensity of settlement occupation (Strauss, 1989; Brace, 2000; Fernandez Peris et al., 2012; March et al., 2012; Aldeias et al., 2016). In particular, the a posteriori estimated duration of the hearth offers a window onto the usage of the fire, which may range from a simple cooking fire to complex metal-firing workshop furnaces. Fire causes soil alterations, which occur at precise temperatures and remain

This is an open access article under the terms of the [Creative Commons Attribution-NonCommercial-NoDerivs License](http://creativecommons.org/licenses/by-nc-nd/4.0/), which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

© 2024 The Author(s). *Archaeometry* published by John Wiley & Sons Ltd on behalf of University of Oxford.

<footnote n="1">¹IPR, UMR 6251, CNRS, Université de Rennes, Rennes, France</footnote>

<footnote n="2">²Laboratorio de Arqueología y Paleoambiente, Instituto de Alta Investigación, Universidad de Tarapacá, Arica, Chile</footnote>
for a long time (Canti & Linford, 2000; Ulery & Graham, 1993). The duration of heating can be determined by computational methods from the depth of reddening (Laloy & Massard, 1984) if the thermal properties of the soil are known (Brodard et al., 2016). Other kinds of alteration may occur due to a strong and rapid rise of temperature, including changes of Fe oxide, with its concomitant change of colour, changes of some magnetic properties (at about 800°C; Tsatskin & Gendler, 2016) or the formation of anhydrite concretion on superficial layers.

Determining the thermal diffusivity of a soil, however, is not an easy task because soils are generally heterogeneous. Therefore, extracting a few samples for delayed analysis in the laboratory is usually not relevant. *In-situ* measures are preferable, and a simple method is to heat the soil and then to deduce its thermal diffusivity from the temperature recorded at several points (Laloy & Massard, 1984; March et al., 2012). Although inverse methods have some constraints and drawbacks – discussed in the paper – they appear the only practical way to measure the thermal properties in a global way (Mansour et al., 2016). For this study, we chose to apply inverse methods only to dry soils, since the treatment of wet soils requires a greater amount of information and has a greater number of unknown variables. Treating humidity requires measuring its amount during the experiment, but the data we have available only include temperature measurements.

The archaeological site Quebrada Maní 12 (QM12; 21°S, 1240 masl, placed on top of a fluvial Miocene terrace) is one of the few well-dated hunter-gatherer open-air camps (Figure 1, a). Hunter-gatherers colonized the core of the northern Chile Atacama Desert by the end of the Pleistocene when its current extreme hyper-arid conditions were ameliorated by increasing precipitation over the western slope of the Andean Cordillera and the Sierra Moreno, a lower mountain range, 80 and 30 km from the site, respectively. Consequently, this territory was not a barrier for early human dispersal throughout South America, and a suitable ecosystem for niche formation c.13,000–12,000 cal years BP (Dillehay, 2014; Gayo et al., 2012; Joly et al., 2017; Latorre et al., 2013; Quade et al., 2008; Santoro et al., 2017). At QM12open camp a prepared fireplace was dated to 12,200 and 11,900 cal years BP (Latorre et al., 2013). Despite these archaeological features, we still do not know certain aspects about the activities of these people. For instance, radiocarbon dates and the stratigraphic refuse accumulation of this site show that people revisited this camp for around 800 years, yet the frequency and time of permanence are still unknown.

In this paper, we briefly describe the Laloy and Massard (1984) method, which is a simple way to determine the thermophysical properties of a soil, but which has too many constraints to be applied in real situations, whereas a numerical global inverse problem (GIP) method based on computer simulations is more useful even though it requires homogeneous soil (Wan et al., 2022). Consequently, for the heterogeneous soils of Quebrada Maní, we developed a numerical variant of the inverse method, called a local inverse problem (LIP), which we validated with laboratory experiments. The results of these methods applied to hearth replications are then presented, which in turn are discussed for the conclusions of this study.

# MATERIAL AND METHODS

## Determination of the thermophysical properties of the soil using inverse methods

All soils can be considered as composite media. They are made of a solid matrix (grains, for example) and empty pores which can be filled with air and/or water.

Heat transfer by conduction in a granular medium is complex. Even if the pores contain only dry air, estimating the effective thermal conductivity is a problem. This property strongly depends on the structure of the grain pile and the quality of the solid contacts between the
<figure><img src="image_2.png" /><figcaption>**FIGURE 1** (a) Location of the archaeological site quebrada Maní 12 (QM). (b) Stratigraphic view of the soil (scale given by the colour rule; image by Renaud Delannay). Layer 1 is very loose (mainly sand); layer 2 contains a mixture of loose sediments and salt blocks; layer 3 is the most rigid layer with the highest content of compacted salts, difficult to excavate; and layer 4 was not studied in this work. The most frequent minerals reported in each layer correspond to: Plag = plagioclase, Quar = quartz, gyps = gypsum, Anhy = anhydrite, Musc = muscovite, then = thenardite, and Glaub = glauberite.</figcaption></figure>

grains. Another problem is the gas expansion during heating because gas thermal properties change significantly with temperature (e.g., the thermal diffusivity of dry air is multiplied by 4 between 25 and 325°C at atmospheric pressure). Humid granular media are even more complex, especially when the temperature is high enough to induce evaporation. Liquid water entrapped between the grains may absorb an important part of heat energy when changing from liquid to vapour. Vapour in the gas changes its thermophysical properties. Moreover, liquid bridges between the grains can increase the thermal conductivity dramatically (Canot et al., 2016; Smits et al., 2010).

Since thermal properties fluctuate greatly at the pore scale, one must work with averaged thermal values, and heat transfer through porous material must be considered at a larger scale called the representative elemental volume (REV). All soil components have their own thermal properties, and it is worth noting that averaged properties at the REV scale are not easy to predict theoretically, mainly because some of them depend on the topology of the pores (Kunii & Smith, 1960; Canot et al., 2016). Hereinafter, homogeneity and heterogeneity will always refer to the REV scale.
An inverse method needs two ingredients. The first is called the forward problem: it is a method (usually a numerical one) used to compute the temperature field in a given domain, knowing the thermal properties of the material and appropriate boundary conditions. The temperature is usually time-dependent, though it can also be stationary. The inverse problem consists of finding the thermal diffusivity, knowing only the temperature at some finite locations and times if applicable. It is formulated as a non-linear least-square problem (Mansour et al., 2016; Muhieddine et al., 2012): after choosing an initial guess for the diffusivity, we iterate it by minimizing the discrepancy between the known temperatures and those provided by the forward problem. In our case, the validation of this procedure has been achieved using synthetic data. For details, see the supplemental data online.

## The Laloy and Massard method

The method proposed by Laloy and Massard (1984) is well adapted to *in situ* estimation of the soil diffusivity if it is homogeneous. This simple method determines the duration of a hearth, knowing the thickness of some altered zone in the upper part of the heated soil, and using a well-known mathematical solution (the ‘Gauss Error Function’) of the diffusion equation in a semi-infinite medium of constant thermophysical properties (Carslaw & Jaeger, 1959). This analytical solution is the forward problem. Their method includes two steps: (1) the determination of the thermal diffusivity of the soil by heating it under some conditions and measuring the temperatures as a function of time at some location inside the medium; and (2) the recovery of the duration of a real or ancient hearth using the estimation of the thermal diffusivity, the thickness of a thermo-altered layer and the temperature threshold at which this change occurs. The main advantage of this method is its simplicity. The computation can be easily programmed using an Excel-type spreadsheet thanks to a linearization of the equations. The assumptions of this method are, however, numerous. Among them, we must emphasize the following which have caused some issues in this study:

* Heat diffusion must be one-dimensional (1D). The thermal sensors for measuring the temperature under the heating device must be centred and very close to it. Such a constrained configuration tends to increase the relative error in sensor positioning.

* The initial heating condition is assumed to be a sudden jump in temperature. Reproducing this condition is difficult; it is not compatible with a wood fire.

* The soil must have constant (in both space and time) thermophysical properties. This means that these properties must be spatially uniform and should not depend on temperature.

These constraints are severe, and when some of them are not verified, the method can yield erroneous results. For example, the soil must not be compacted too much by the insertion of the thermal probe, otherwise an artificial heterogeneity could appear (Grott et al., 2010). For the determination of the thermal diffusivity, we could apply another global inverse method, which has fewer constraints than the method of Laloy and Massard. For heterogeneous soils, however, global inverse methods must be discarded, and we had to develop a local method.

## Two-dimensional (2D) axisymmetric model and inverse methods

For several years the ArPhyMat project (http://arphymat.univ-rennes1.fr/) was dedicated to numerical modelling of heat and mass transfers in soils below fires. In this framework, appropriate numerical tools have been developed and validated by both laboratory and field experiments (Muhieddine et al., 2011). The numerical model and tools presented here concern only
dry soil without chemical reaction or phase change. The choice of a forward model involving a 2D axisymmetric geometry seems to be a good compromise between performance and reliability: a 1D model also has constraints, and a fully three-dimensional (3D) model would lead to costly simulations, not to mention that the ideal number of thermal sensors would be very great in this last case. The constraints of this 2D axisymmetric numerical model are much less severe, and are summarized below:

* The geometrical configuration and any possible heterogeneity, as well as the initial temperature of the soil and the heating, must be axisymmetric (invariant by rotation around the vertical axis).

* The heating device must be circular (due to the first constraint) and its temperature profile must be radial at each time.

* The soil surface outside the heating device should be insulated. This is generally well verified since air has a lower thermal conductivity than the soil.

* The medium temperature does not need to be initially uniform but the axi-symmetry constraint remains.

The GIP uses all sensors for the whole measurement domain to determine a global estimation of the thermal diffusivity, except that one sensor (the highest one) is used to obtain the upper boundary condition along the time. This is preferable: the result is more reliable because it is less sensitive to the uncertainties of the sensor positions. This is only relevant, however, when the soil is homogeneous, so that the diffusivity is constant throughout the whole domain.

When we suspect that the soil is not homogeneous, we must develop a local method. The LIP uses only a subset of three sensors, generally close to each other, to determine a local value (more precisely, a mean local value in a small volume containing the three sensors) of the thermal diffusivity, with a 1D-underlying model since the control volume used is small in comparison with the whole domain. Though this second approach may appear more interesting, it uses a 1D equation and assumes that the diffusivity is constant in the domain that embeds the three sensors. Specifically, two of these sensors are used to give the temperature boundary conditions and only the values of the middle sensor drive the inverse problem iterations. Outside the subdomain defined by the three sensors, the soil may have any heterogeneous thermal properties. The LIP is more sensitive to errors because only one sensor's data are used to find the local diffusivity, but it is complementary to GIP since it allows local measurements.

These new inverse methods are more versatile than the Laloy and Massard method, and can be used to estimate the fire's duration from the diffusivity, based on the thickness of the thermo-altered layer at the soil surface. Such numerical recovery is relatively easy to perform when the soil structure is homogeneous. Irregular thickness of the thermo-altered zone shows how heterogeneous are the thermophysical properties of the soil. In this kind of soil, accurate measurements of diffusivity are difficult to obtain. A layered medium, where the depth of the thermo-altered layer is uniform, will not prevent highly erroneous measurements. Appendix S1 in the supplemental data online discusses a special case where a deeper layer, unknown after the excavation, has thermal diffusivities that significantly differ from the diffusivity of the upper layer. These results show the difficulties in achieving accurate estimation of the time duration of a fire when the thermal properties of soil are not well known.

## LABORATORY EXPERIMENTS

Laboratory experiments are necessary to understand how heat transfers occur in granular media, and how the thermal properties are affected by humidity. They also validate
numerical codes based on simplified physical models, and provide a way to improve field experiments, as well as give some insight into the information that should be prioritized in those experiments.

## Experimental set-up

The experiments were conducted in our laboratory at the Institut de Physique de Rennes (IPR), France. A cylindrical box (30 cm in diameter, 30 cm high) was filled with Fontainebleau white sand whose grains are approximately spherical, 200 µm in diameter. The surface was heated by an infrared electric bulb (power 100 W) located inside a cylindrical aluminium foil (diameter 12 cm) to get a uniform infrared flux on the sand surface after multiple reflections. All was done in a way to ensure that the configuration was axisymmetric (Figure 2, a). Heating duration ranged from 1 h to a few hours. A few thermal probes were located at different depths (ranging from zero to 5 cm deep) (Figure 2, a). Note that the rod diameter (2 mm) is small enough in comparison with the box radius (40 cm) to avoid a perturbation of the temperature field by the presence of these metal pieces. Even using 10 sensors does not affect very much the heat flow across the medium. The vertical position of the thermocouples rod leads to an accurate location of their heads. Initial temperature was about 25°C and its time evolution at each probe was recorded. The heating of the porous media took around 1 h to obtain a substantial temperature increase for the deeper probe.

## Role of the humidity in the heat transfer

To demonstrate the role of humidity, we performed two different experiments. Humidity in porous medium or hydration, even in low quantity, may have a significant influence on the thermal properties and strongly affects the thermal energy balance. Min and Emmons (1972) experimentally studied such a situation and emphasized the consequences, but for relatively large humidity (the liquid fraction was 20%). We claim that the influence of humidity is noticeable even for a liquid fraction lower than 10% (Figure 2, c). At a small liquid fraction, a plateau temperature can be observed as low as 35°C.

In a first experiment (dry case), sand was stored for few weeks in a room with a standard relative humidity of about 40-50%. This state of equilibrium implies that some water vapour condenses inside the porous medium (Kelvin equation effect), but the quantity of liquid is extremely small. There were no apparent internal cohesive forces (the sand easily flows through the fingers when we handled it), therefore we consider it as dry. The aluminium reflector bottom was placed in contact with the sand (Figure 2, a). In Figure 2 (b) shows that temperature curves were strictly increasing, with smooth variations because the heating was constant. Particularly, the slope of each curve is continuous, without any break (it will be shown below that the presence of humidity leads to a slope break).

For the second experiment, we filled a bowl with dry sand, spread the appropriate amount of water on its surface and mixed it with the sand. The humid sand was placed in the cylindrical box and packed layer by layer using a flat plate with a constant vertical force. This procedure aims to remove most of air pockets naturally contained in humid sand. The initial theoretical volume ratio of liquid to pore (called the water saturation) can be predicted by means of simple formulas, knowing the initial volume ratio of water to sand, on one hand, and the measured solid fraction of our dry sand (0.62), on the other. To measure the humidity of some samples, we weighed the humid sand, dried it in an oven and weighed it again. We found a water saturation of 8% instead of the predicted 10%, due mainly to the variation of the pore volume during the process. During heating, the aluminium reflector was moved up from the sand surface by
<figure><img src="image_3.png" /><figcaption>FIG URE 2 (a) Schematic of the laboratory heating device. Multiple reflections of infrared rays inside the aluminium cylinder produce a uniform heat flux on the sand surface. Both the sand box and the heating are axisymmetric. (b) Laboratory experiment 1 (LE#1) in dry sand; smooth curves stand out. (c) Laboratory experiment 2 (LE#2) in wet sand. The curves may present an inflection point, related to the evaporation of a low water content in the sample. Temperature is as a function of time, at different positions (cm) of the sample, with saturation $S = 8\%$, heated by a constant flow infrared device.</figcaption></figure>
1 cm to allow water vapour to escape. Figure 2 (c) shows some curves (green, light blue, orange) presenting a slope decrease followed by a slope increase. This is typical of phase change and is due to energy absorption (latent heat) to change liquid water into vapour. Usually, a sudden slope jump is observed during a phase change in the saturated case (Muhieddine et al., 2009). In our case study, given that the water content is always small, the jump is weak and smoothed by the fact that water liquid bridges are disconnected, allowing the vapour to move and diffuse easily in the porous medium (Min & Emmons, 1972). The upper curve (red) records the temperature at the surface of the sand; the perturbation of the values comes from air displacement in the experimental room.

It should be emphasized that, as opposed to the saturated case ($S = 100\%$) or the classical Stefan problem (https://en.wikipedia.org/wiki/StefanProblem), the slope jump observed in the temperature curves did not occur at $100^{\circ}$C, which is the boiling point of water at atmospheric pressure. Here, the strong shift in boiling temperature (in Figure 2, c, the slope break occurs from 35 to $50^{\circ}$C) is related to the very small size of the concave liquid bridges between the grains and to Laplace's law. Anyway, a true plateau at $100^{\circ}$C, or a pseudo-plateau at lower temperatures, can be related to water inside the porous medium (Min & Emmons, 1972).

## Validation of the LIP

To validate the LIP, the experimental set-up of Figure 2 (a) was used to carry two different experiments. For the first, we used Fontainebleau sand, and for the second, we used fine cooking sea salt La Baleine™ (grains diameter around 500 µm). The thermal properties for our two granular media are affected by the compacity of the pile of grains, the shape of the grains and humidity, which can increase greatly the thermal conductivity, especially for salt. General data found in the literature cannot be applied directly to validate our method. An estimation of the thermal diffusivity for each granular medium was obtained by the global inverse method (Mansour et al., 2016), via a numerical processing of the recorded temperatures. The use of all sensors prevents an erroneous diffusivity value due, for example, to non-uniform heating or to an imperfectly flat surface of the material. The estimation for the thermal diffusivity of the sand is $2.0 \times 10^{-7} m^2/s$, whereas it is $2.7 \times 10^{-7} m^2/s$ for that of the salt (data obtained by GIP from Figure S4 in supplemental data online).

The LIP needs temperature evolution at three probes. As there are many probes in the porous medium, we could select different triplets in the same experiment. In some cases, especially when the assumption of 1D heat transfer is not well verified, the method can lead to odd results, with out-of-range numerical results. For that reason, such aberrant results were removed (but for the whole data of results, see in supplemental data online, Table S2). Moreover, there are two ways of interpolating the initial temperature profile: a linear interpolation and a quadratic one. The comparison between the two results gives an idea of this error due to the linearity of the initial temperature.

As shown in Figure 3, the resulting estimated thermal diffusivity may vary depending on the selected triplet. First, because all sensors are not located on a vertical line (see Figure 2a); and second, because some sensors are more sensitive to a non-horizontal heating or non-plane isotherms (violation of the 1D model assumption).

In sum, the LIP gives correct numerical results if the local temperature isotherms are flat and horizontal, and the location of the sensors is known with precision. For heterogeneous materials, the three sensors must be located within a domain whose length is less than the length of the heterogeneity.
<figure><img src="image_4.png" /><figcaption>**FIGURE 3** Results of the application of local inverse problem (LIP) to the sand (blue bars) and salt (red bars) samples. Each bar corresponds to the thermal diffusivity measurement obtained with a set of three thermocouples. Dashed lines represent the expected values, which were obtained from all sensors.</figcaption></figure>

# APPLICATION OF THE INVERSE METHODS TO HEARTH REPLICATIONS IN THE ATACAMA DESERT

## Field experiments: hearths replication

Three experimental hearths (FE-L1–L3) were conducted in a manner like that described by March et al. (2012); they were located in a natural soil of one of the Miocene alluvial remnants, about 400 m south of QM12. The soil stratigraphy consists of three layers (Figure 1, b): layer 1 beneath the desert pavement contains mainly silt, plus gypsum and anhydrite; layers 2 and 3 were composed of sand, very fine white silt and caliche, respectively (Bastos, 2018; see also Latorre et al., 2013: fig. 5, for a detailed description of an exposed section of site QM12). Caliche, strong salt blocks (Figure 4, d), results from the transport of the salt by some internal flow of soil moisture (Ugalde et al., 2020). Such blocks have been reproduced in our laboratory by drying a mix of sand, salt and water. Due to the strong compactness of this type of structure, their thermal conductivity can be up to four times those of dry granular salt. For a brief description of the experiments and the raw numerical data, see https://arphymat.univ-rennes1.fr/data/august_2015.html. For resampled temperature data adapted to our inverse method, see the supplemental data online.

The boundary condition temperature used in the numerical model comes from the upper thermal sensor located just below the surface. In Figure 4 (a–c), the dashed black line is the temperature of this upper sensor, and we take its position as the reference depth $z = 0$. The coloured solid curves correspond to the other sensors, whose relative depth $z$ is indicated on the right side of the graph. The continuous curves correspond to the measured temperature,
<figure><img src="image_5.png" /><figcaption>FIGURE 4 (a–c) Field experiments FE-L1, FE-L2 and FE-L3. The dashed black line is the fire temperature; solid coloured curves are the experimental temperatures; and dotted curves are the temperatures reconstructed using the global inverse problem (GIP). (d) The white fragments found during the hearth excavation correspond to the presence of some salts (caliche). The data shown in the grey rectangle are not used for technical reasons; this does not affect the application of the inverse methods.</figcaption></figure>
whereas the dotted ones are constructed from the thermal diffusivity estimated by GIP. Note that not all the sensors are located along the same vertical line, but are arranged in a relatively small volume surrounding the central axis (see the supplemental data online, sheet ‘Other Data’).

A sound practice is to use both field and laboratory experiments and numerical simulations (based on simplified models). It is essential first to examine the experimental data to identify the cases where the physical constraints of the model are satisfied. If we treat them blindly, we run the risk of making serious errors of interpretation. Therefore, we performed qualitative analyses to obtain useful feedback that could constrain the application of the LIP.

## Applying LIP

First, it should be noted that the soil of the experiments was not homogeneous (Figure 1, b), so GIP is not appropriate. Its use, however, allows us to obtain reference temperature curves and to justify the choice of the three sensors for each LIP application by comparing the measured temperature curves with the reference ones. This comparison provides useful information about spatial location (where the soil may be heterogeneous) and time (when moisture may disappear completely). Considering two sensors located at different depths, if we visualize the two measured temperature curves with the two simulated ones, we conclude that GIP finds an average value between two different diffusivities, since the inverse method minimizes the residue, that is, the cumulative difference between the experimental and synthetic curves. In this case, we can safely conclude that each of the two sensors experienced a different thermal diffusivity. On the other hand, if the experimental curves show a plateau that disappears after a certain time while the simulated curve does not show such a plateau, we can conclude with certainty that there was some moisture at the beginning of the heating, but it disappeared after a given time. This leads us to use only the experimental data related to the dry time interval. All the thermal diffusivity estimations obtained from GIP and LIP are reported in Table 1.

Figure 4 (a) shows the temperature curves of the FE-L1 field experiment. The red, purple and orange experimental curves show discrepancies with their reconstruction by GIP. This discrepancy reveals the presence of water, as (1) the temperature slope decreases before increasing again (reminiscent of the phase change described in the Laboratory Experiments section) and (2) there is a weak slope jump between 70°C (orange curve) and 85°C (red curve). This moisture may come from free water within the pores and/or hydrated minerals. Dense fog, known locally as *camanchaca* or *garúa*, is a common phenomenon originating on the Pacific coast that advances several kilometres inland fully reaching the locality of QM12. Although this fog does not produce rain, it does contain moisture droplets measuring between 1 and 40 µm. During the night and early morning hours, some of this moisture condenses on the cold surface of the ground. Another source of moisture, although less frequent, comes from organic matter contained in the soil (Aldeias et al., 2016). The stratigraphy of hearth FE-L1 (Figure 1, b) shows that some gypsum was present between 6 and 13 cm depth. Gypsum (dihydrate) starts to lose water at 60°C and produces bassanite (hemihydrate) through an endothermic reaction. Consequently, given these potential energy sinks, the application of GIP to this experiment cannot provide a good approximation of diffusivity. Figure 4 (a) shows, however, that the soil at $z = -6.0$ cm is dry for more than 2.2 h. Therefore, LIP is applicable for this time interval by using the first three sensors (black, green and red dashed) that are in the dry zone.

The soil in the third field experiment (FE-L3) (Figure 4, b) was dry. The temperature data, shown in the red and orange continuous curves, cross the dashed ones at some point, which may be due to a non-uniform initial temperature condition. The orange and purple sensors show inconsistent temperature curves. Although the purple sensor is located deeper than the orange one, the purple curve starts at a higher temperature than the orange one. Moreover, the
<table><caption>TABLE 1 Summary of soil thermal diffusivities estimated from three different field experiments at quebrada Maní.</caption><thead><tr><th>Experiment and date</th><th>Thermal diffusivity (m²/s, one GIP estimation)</th><th>Sensors used for LIP</th><th>Thermal diffusivity (m²/s), three LIP estimations</th><th>Remarks</th></tr></thead><tbody><tr><td>FE-L1<br>13 August<br>2015</td><td>2.0 × 10<sup>-7</sup></td><td>Dashed<br>black +<br>green + red</td><td>α<sub>1</sub> = 2.5 × 10<sup>-7</sup><br>α<sub>2</sub> = 2.7 × 10<sup>-7</sup><br>α<sub>3</sub> = 1.6 × 10<sup>-7</sup></td><td>Upper dry layer of the soil (about 6 cm thickness) containing gypsum (as found in layer 1 of Figure 1, b)</td></tr><tr><td>FE-L2<br>17 August<br>2015</td><td>2.7 × 10<sup>-7</sup></td><td>Dashed<br>black +<br>orange +<br>red</td><td>α<sub>1</sub> = 0.33 × 10<sup>-7</sup><br>α<sub>2</sub> = 0.43 × 10<sup>-7</sup><br>α<sub>3</sub> = 0.12 × 10<sup>-7</sup></td><td>Two sensors in the upper layer of the soil (about 4 cm thickness)</td></tr><tr><td>FE-L3<br>15 August<br>2015</td><td>2.8 × 10<sup>-7</sup></td><td>Orange +<br>red +<br>purple</td><td>α<sub>1</sub> = 44.1 × 10<sup>-7</sup><br>α<sub>2</sub> = 44.1 × 10<sup>-7</sup><br>α<sub>3</sub> = 17. × 10<sup>-7</sup></td><td>Two sensors below 1 cm depth in a layer of soil containing salt blocks (as found in layer 2 of Figure 1, b)</td></tr><tr><td></td><td></td><td>Dashed<br>black +<br>blue + red</td><td>α<sub>1</sub> = 4.3 × 10<sup>-7</sup><br>α<sub>2</sub> = 4.3 × 10<sup>-7</sup><br>α<sub>3</sub> = 3.5 × 10<sup>-7</sup></td><td>Upper layer of the soil (about 4 cm thickness) (probably again layer 2 of Figure 1, b)</td></tr></tbody></table>

GIP = global inverse problem, using all sensors. LIP = local inverse problem, using only three sensors, where selection is detailed in the third column. For LIPs, α₁ is the raw value, using a linear temperature profile for the initial conditions; the same for α₂, but using a parabolic temperature profile instead a linear one. α₃ is obtained by perturbing the middle sensor position up to the point where the residue reaches a minimum value. Given layer identification is coherent with soil examination conducted during excavation.

two experimental curves cross after 2 h. This could be explained by the fact that as the sensors were not in a vertical line, the horizontal displacement of the heating centre above the ground causes a change in the inclination of the isothermal surfaces, which in turn causes an apparent change in the depth of the sensors. Consequently, the LIP was applied using the top three sensors.

In the second field experiment (FE-L2) the soil was also dry. The orange and green curves in Figure 4 (c) correspond to sensors located at $z = -1.1$ and $-1.3$ cm depth, respectively. The dashed black line shows a continuous increase (except for the noise) in fire temperature, up to 5.5 h. Unexpectedly, the orange and green curves, on the other hand, are flat after 4 h, and even decrease well before 5.5 h. Theoretically, the cooling should first appear on the surface and then propagate downward to the other sensors. But theodolite space measurements show that the orange and green sensors are horizontally offset from black by 3 cm (green) and 5 cm (orange). The unexpected drop would be explained by the heterogeneity of the surface temperature: the green and orange sensors are under a part of the fire that started to cool down much earlier than the part containing the black sensor. Nevertheless, the behaviour of the temperature curves seems reasonable for all sensors during the first half of the time interval. A closer inspection of the curves, however, shows that the soil was not homogeneous. The fact that two curves (orange and green) are below the dotted curves of the same colour, while two other curves (red and purple) are above the corresponding dotted curves, indicates different thermal diffusivities for depths $z = 0$ to $-1.3$ cm, and $-4.3$ cm. This is consistent with the stratigraphic conditions of the soil, in which compact and hard salt blocks stand out, strongly affecting its diffusivity (Figure 4, d). To apply LIP, we have to discard the green curve; it is too close to the orange one. We apply LIP twice as follows: one simulation for the first volume, using the first three sensors: black, orange and red (Figure 4, c); and one for the second volume (using the last three sensors: orange, red and purple) (Figure 4, c). Moreover, to avoid the above-mentioned
problem due to the displacement of the heat centre, we restrict the duration to the first half-time. It can be observed that these two volumes overlap and that each contains some heterogeneity. This last choice of triplets was deliberately made to test LIP when a jump in diffusivity occurs. It would also be interesting to apply the same analysis to synthetic data to see how the method performs in the presence of such a jump.

## Results

Table 1 presents a summary of all measured values. For the four presented cases, we computed the soil thermal diffusivity. Two main uncertainties affect their values. The first is difficult to quantify: it comes from the fact that the initial temperature profile between the three points of LIP is not known; instead we have only three temperature points. We then made two different computations, using a piecewise linear temperature interpolation with respect to the depth, giving the value $α_1$, and a parabolic interpolation, giving the value $α_2$. The second source of uncertainty is the sensor location, which is of the order of 1 mm. To see its influence, we perturb the location of the middle sensor up to the point where the residue in the LIP reaches a minimum value. This gives the value $α_3$. Comparing the three values then gives an idea of the overall uncertainties and also shows by coherence if the homogeneity hypothesis is valid. The results of the FE-L2 experiment confirm, as predicted by the behaviour of the GIP curves, that the means of the two selected layers are very different: the upper layer has a very low thermal conductivity, while the lower layer has a much higher value, which corresponds to a very heterogeneous soil. For the two cases, $α_3$ is always very different from the two first values, showing a high sensitivity to the sensor position. Given these conditions, the large discrepancy in the calculated diffusivities for the second layer confirms that they are unreliable; this can be caused by heterogeneity within the LIP volume.

For each of the FE-L1 and FE-L3 experiments, by comparing the diffusivities reported in Table 1, we can observe similar values.

## CONCLUSIONS AND PERSPECTIVES

We present the local inverse method (LIP), a new method for determining the thermal diffusivity of a soil. We validated it in the laboratory and applied it to numerical data available from three field experiments in the Atacama Desert, to demonstrate its applicability to the study of campfire soil replicas. These numerical data, consisting of temperature inside the soil during a heating process, came from experiments carried out on soil of the Quebrada Mani 12 archaeological site. A first stage of our analysis consisted of applying, blindly, a global inverse problem (GIP) for each experiment, keeping as many sensors as possible to reveal the defects of each experiment and to identify the cases that satisfied the physical constraints of LIP (dry and homogeneous part of the soil containing at least three sensors). In a second stage, LIP was applied to sets of three sensors selected among the available sensors. This approach revealed the heterogeneity of the subsoil in terms of its thermal behaviour. The scattered diffusivities allows us to point out, in general, that the estimation of the subsurface thermal diffusivity field in the soil of the Atacama Desert is a complicated problem.

Consequently, we suggest considering the following factors for future analysis: (1) the presence of water, either free or in hydrated minerals, since even a small amount of moisture can affect the temperature curves, leading to erroneous thermal diffusivities (cf. section Laboratory Experiments); (2) heterogeneous thermal properties of the soil due to the influence of different structures, such as salt, weakly cemented silt with gypsum and anhydrite known locally as chusca; (3) heterogeneity restricts the use of the GIP and requires many sensors to achieve
meaningful results using LIP; (4) non-uniform fire temperature, due to the position of the wood pieces in the hearth, or partial activation of the embers by the wind, may violate the axisymmetric condition of our GIP; and (5) the use of flexible rod thermocouples in heterogeneous or too hard soils is inappropriate, given the uncertainty about the positions of their head, which in turn hinders the application of LIP.

Improvements in field experimental methods for complex soils such as those of the Atacama Desert should considered the following:

* Precise control of the heating process during field experiments.

* Sufficient thermocouples should be well placed according to the heterogeneity of the soil. But note that the heterogeneity scale cannot be too small, because having too many sensors can also disrupt the temperature field of the heterogeneous soil to be measured.

* During the heating process, avoid possible horizontal displacements of the heater on the soil. This ensures the generation of more accurate information and data.

* If water is present in the soil and evaporates or condenses, the inverse method must be extended to reach a more complete physical model (Mansour et al., 2016), including, in addition to calculating the temperature, the evolution of humidity that could be measured with capacitance technique (Louge et al., 2013).

We suggest the use of an experimental set-up such as that of Laloy and Massard (1984), in which the positions of the thermocouples are precisely known, since they are fixed to a rigid, vertical, cylindrical iron bar. Moreover, a heating plate with a constant, easily controlled energy source, such as electricity or gas, should be used. In all cases, quantitative relationships between the relaxation of some constraints and the accuracy of the results can be calculated by further numerical simulations. Ideally, future technological developments could create devices with special thermal sensors that combine, in a small volume, several thermocouples to measure temperature, heat flux and their 3D orientation.

In summary, physical analysis and inverse methods applied to archaeological fire replicates in the Atacama Desert show that the combination of data obtained from field and laboratory experiments and numerical simulations is a reliable way to identify and understand the dominant physical processes involved in archaeological hearth fires (thermal conduction, phase change, salt transport). The characterization of the thermal behaviour of Atacama Desert soils clearly shows, despite the inevitable uncertainties, that thermal diffusivities could be significantly higher than those of dry but non-saline, clayey or sandy soils (the latter soils usually have a thermal diffusivity close to $2 \times 10^{-7} \text{ m}^2/\text{s}$). Future laboratory experiments should test whether some salts, in association with other factors such as transient moisture, for example, can influence soil thermal behaviour. Additional laboratory experiments on the heating of hydrated minerals (such as gypsum) should also be considered. Future work in our laboratory will validate the proposed LIP with synthetic data perturbed in different ways.

## ACKNOWLEDGEMENTS

We thank Ramiro March, Cecilia Rodríguez-Loredo, Wilfredo Faundes and Felipe March for carrying out the field experiments at Quebrada Maní in the Southern Hemisphere winter of 2015. We also thank Sean McNamara for his help and for fruitful discussions.

## DATA AVAILABILITY STATEMENT

The raw data that support the findings of this study are available at https://arphymat.universnes1.fr/data/august_2015.html.

## ORCID

Renaud Delannay https://orcid.org/0000-0002-8361-9604
REFERENCES

Aldeias, V., Dibble, H. L., Sandgathe, D., Goldberg, P., & McPherron, S. J. P. (2016). How heat alters underlying deposits and implications for archaeological fire features: A controlled experiment. *Journal of Archaeological Science*, 67, 64–79. https://doi.org/10.1016/j.jas.2016.01.016

Bastos, G. (2018). *Caractérisations physico-chimiques des sols archéologiques du Nord du Chili*. Unpublished master thesis. Université de Rennes 1. Available at: https://arphymat.univ-rennes1.fr/publis/2018BastosM2Internship.pdf

Brace, C. L. (2000). *Evolution in an anthropological view*. AltaMira Press.

Brodard, A., Lacanette-Puyo, D., Guibert, P., Lévêque, F., Burens, A., & Carozza, L. (2016). A new process of reconstructing archaeological fires from their impact on sediment: A coupled experimental and numerical approach based on the case study of hearths from the cave of les Fraux (Dordogne, France). *Archaeological and Anthropological Sciences*, 8, 673–687. https://doi.org/10.1007/s12520-015-0250-7

Canot, É., Delannay, R., Mansour, S., Muhieddine, M., & March, R. J. (2016). Effective thermal conductivity of a wet porous medium – Presence of hysteresis when modeling the spatial water distribution for the pendular regime. *Journal of Heat Transfer*, 138, 091011. https://doi.org/10.1115/1.4032950

Canti, M., & Linford, N. (2000). The effects of fire on archaeological soils and sediments: Temperature and colour relationships. *Proceedings of the Prehistoric Society*, 66, 385–395.

Carslaw, H. S., & Jaeger, J. C. (1959). *Conduction of heat in solids*. Oxford University Press.

Dillehay, T. D. (2014). Entangled knowledge: Old trends and new thoughts in first South American studies in *Paleoamerican Odyssey*. In K. E. Graf, C. V. Ketron, & M. R. Waters (Eds.), *Center for the Study of the first Americans* (pp. 377–396). Texas A&M University press. College Station

Fernandez Peris, J., Barciela Gonzalez, V., Blasco, R., Cuartero, F., Fluck, H., Sanudo, P., & Verdasco, C. (2012). The earliest evidence of hearths in southern Europe: The case of Bolomor cave (Valencia, Spain). *Quaternary International*, 247, 267–277. https://doi.org/10.1016/j.quaint.2010.10.014

Gayo, E. M., Latorre, C., Jordan, T. E., Nester, P. L., Estay, S. A., Ojeda, K. F., & Santoro, C. M. (2012). Late quaternary hydrological and ecological change in the hyperarid core of the northern Atacama Desert (~21°S). *Earth Science Reviews*, 113, 120–140. https://doi.org/10.1016/j.earscirev.2012.04.003

Grott, M., Knollenberg, J., & Krause, C. (2010). Apollo lunar heat flow experiment revisited: A critical reassessment of the *in situ* thermal conductivity determination. *Journal of Geophysical Research*, 155, E11005. https://doi.org/10.1029/2010JE003612

Joly, D., Santoro, C. M., Gayo, E. M., Ugalde, P. C., March, R. J., Carmona, R., Marguerie, D., & Latorre, C. (2017). Fuel management and human colonization of the Atacama Desert, northern Chile, during the Pleistocene–Holocene transition. *Latin American Antiquity*, 28, 144–160. https://doi.org/10.1017/laq.2016.8

Kunii, D., & Smith, J. M. (1960). Heat transfer characteristics of porous rocks. *AICHE Journal*, 6, 71–78. https://doi.org/10.1002/aic.690060115

Laloy, J., & Massard, P. (1984). Nouvelle méthode thermique d'étude des foyers préhistoriques. *Revue d'Archéométrie*, 8, 33–40. https://doi.org/10.3406/arsci.1984.1209

Latorre, C., Santoro, C. M., Ugalde, P. C., Gayó, E. M., Osorio, D., Salas-Egaña, C., De Pol-Holz, R., Joly, D., & Rech, J. A. (2013). Late Pleistocene human occupation of the hyperarid core in the Atacama Desert, northern Chile. *Quaternary Science Reviews*, 77, 19–30. https://doi.org/10.1016/j.quascirev.2013.06.008

Louge, M. Y., Valance, A., Ould el-Moctar, A., Xu, J., Hay, A. G., & Richter, R. (2013). Temperature and humidity within a mobile barchan sand dune, implications for microbial survival. *Journal of Geophysical Research*, 188, 2392–2405. https://doi.org/10.1002/2013JF002839

Mansour, S., Canot, É., & Muhieddine, M. (2016). Identification of the thermophysical properties of the soil by inverse problem. *Journal of Heat Transfer*, 138, 091010. https://doi.org/10.1115/1.4032947

March, R. J., Lucquin, A., Joly, D., Ferreri, J. C., & Muhieddine, M. (2012). Processes of formation and alteration of archaeological fire structures: Complexity viewed in the light of experimental approaches. *Journal of Archaeological Method and Theory*, 21, 1–45. https://doi.org/10.1007/s10816-012-9134-7

Min, K., & Emmons, H. W. (1972). *The drying of porous media*. Stanford University Press.

Muhieddine, M., Canot, É., & March, R. J. (2009). Various approaches for solving problems in heat conduction with phase change. *International Journal on Finite Volumes, AMU*, 6, N. 1 <hal-01120384>

Muhieddine, M., Canot, É., & March, R. J. (2012). Heat transfer modeling in saturated porous media and identification of the thermophysical properties of the soil by inverse problem. *Applied Numerical Mathematics*, 62, 1026–1040. https://doi.org/10.1016/j.apnum.2012.02.008

Muhieddine, M., Canot, É., March, R. J., & Delannay, R. (2011). Coupling heat conduction and water–steam flow in a saturated porous medium. *International Journal for Numerical Methods in Engineering*, 85, 1390–1414. https://doi.org/10.1002/nme.3022

Quade, J., Rech, J. A., Betancourt, J. L., Latorre, C., Quade, J., Rylander, K. A., & Fisher, T. (2008). Paleowetlands and regional climate change in the Central Atacama Desert, northern Chile. *Quaternary Research*, 69, 343–360. https://doi.org/10.1016/j.yqres.2008.01.003
Santoro, C. M., Gayo, E. M., Capriles, J. M., de Porras, M. E., Maldonado, A., Standen, V. G., Latorre, C., Castro, V., Angelo, D., McRostie, V., Uribe, M., Valenzuela, D., Ugalde, P. C., & Marquet, P. A. (2017). Continuities and discontinuities in the socio-environmental systems of the Atacama Desert during the last 13,000 years. *Journal of Anthropological Archaeology*, *46*, 28–39. https://doi.org/10.1016/j.jaa.2016.08.006

Smits, K. M., Sakaki, T., Limsuwat, A., & Illangasekare, T. H. (2010). Thermal conductivity of sands under varying moisture and porosity in drainage-wetting cycles. *Vadose Zone Journal*, *9*, 1–9. https://doi.org/10.2136/vzj2009.0095

Strauss, L. G. (1989). On early hominid use of fire. *Current Anthropology*, *30*, 488–491. https://doi.org/10.1086/203772

Tsatskin, A., & Gendler, T. S. (2016). Identification of “red ochre” in soil at Kfar HaHoresh Neolithic site, Israel: Magnetic measurements coupled with materials characterization. *Journal of Archaeological Science: Reports*, *6*, 284–292. https://doi.org/10.1016/j.jasrep.2016.02.027

Ugalde, P. C., Quade, J., Santoro, C. M., & Holliday, V. T. (2020). Processes of Paleoindian site and desert pavement formation in the Atacama Desert, Chile. *Quaternary Research*, *98*, 58–80. https://doi.org/10.1017/qua.2020.39

Ulery, A. L., & Graham, R. C. (1993). Forest fire effects on soil color and texture. *Soil Science Society of America Journal*, *57*, 135–140. https://doi.org/10.2136/sssaj1993.03615995005700010026x

Wan, S., Wang, K., Xu, P., & Huang, Y. (2022). Numerical and experimental verification of the single neural adaptive PID real-time inverse method for solving inverse heat conduction problems. *International Journal of Heat and Mass Transfer*, *189*, 122657. https://doi.org/10.1016/j.ijheatmasstransfer.2022.122657

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Canot, É., Delannay, R., & Santoro, C. M. (2025). Physical analysis and inverse methods applied to archaeological FIRE replications on ATACAMA desert soils, northern Chile. *Archaeometry*, *67*(3), 625–641. https://doi.org/10.1111/arcm.13030

## APPENDIX: THE EFFECT OF A LOWER LAYER OF HIGHER THERMAL CONDUCTIVITY

We have realized from the field experiment FE-L2 that strong differences in local values of the medium conductivity may occur. For example, a deep layer, located below the sensors, with higher thermal conductivity, can yield dramatic behaviour. We applied here the typical Laloy and Massard (1984) method to estimate the time duration of heating, based on the depth of an altered layer observable by its colour difference. We suppose that the soil contains a certain component for which the temperature alteration is 360°C. We also suppose the existence of two homogeneous layers: the upper one of depth −2 cm, mainly constituted by sand, with a thermal diffusivity $\alpha_u = 2 \times 10^{-7} \text{ m}^2/\text{s}$, whereas the deeper one ($z < -2 \text{ cm}$) has a greater thermal conductivity, noted $\alpha_d$. In all the following virtual experiments, the soil surface is heated at a constant temperature of 500°C.

Figure 5 (a) shows the results of 6 h heating: the different curves present the temperature profiles for different deep layer diffusivities, $\alpha_d$, without a change in the upper layer diffusivity, $\alpha_u$. It shows that the temperatures in the upper layer depend strongly on $\alpha_d$. It also shows the modification of the depth of the altered layer, along the vertical dashed red line.

Figure 5 (b) shows that if we ignore the presence of the more conductive deep layer, large errors for the time duration estimation can be obtained. For example, here the estimated duration would be 2 h instead of a real duration of 6 h, when the diffusivity of the unknown layer $\alpha_d$ is $10 \times \alpha_u$.
<figure><img src="image_6.png" /><figcaption>FIGURE 5 (a) Influence of the presence of a deep layer more conductive than the upper one. Temperature profiles were simulated for different diffusivity ratios. (b) Error in the estimation of fire duration in the case of a deeper, more conductive layer (case of ratio = 10). The reference curve used to estimate the heating duration was calculated for a semi-infinite homogeneous soil of usual diffusivity.</figcaption></figure>