## Abstract

The early postmortem interval (PMI), i.e., the time shortly after death, can aid in the temporal reconstruction of a suspected crime and therefore provides crucial information in forensic investigations. Currently, this information is often derived from an empirical model (Henssge’s nomogram) describing posthumous body cooling under standard conditions. However, nonstandard conditions necessitate the use of subjective correction factors or preclude the use of Henssge’s nomogram altogether. To address this, we developed a powerful method for early PMI reconstruction using skin thermometry in conjunction with a comprehensive thermodynamic finite-difference model, which we validated using deceased human bodies. PMIs reconstructed using this approach, on average, deviated no more than ±38 minutes from their corresponding true PMIs (which ranged from 5 to 50 hours), significantly improving on the ±3 to ±7 hours uncertainty of the gold standard. Together, these aspects render this approach a widely applicable, i.e., forensically relevant, method for thermometric early PMI reconstruction.

## INTRODUCTION

The early postmortem interval (PMI), i.e., the time shortly after death, plays a key role in forensic investigations, as it aids in the temporal reconstruction of events. Consequently, the development of a method to determine the PMI remains one of the most important challenges in forensic medicine to date (*1*). In its pursuit, many pathophysiological changes have been investigated as potential measures of the PMI. These measures can be divided into two groups. The first relies on sampling of tissue or bodily fluids (*2*), followed by laboratory examinations to quantify, e.g., nucleic acid degradation (*3*), changes in the ocular potassium concentration (*4*, *5*), and microbial (*6*) and metabolomic (*7*) changes. In contrast, the second group involves probing optical, mechanical, or thermal changes in human tissue (*8*–*12*). This latter group does not require any sample extraction or laboratory examinations; these pathophysiological changes can therefore be quantified directly at the crime scene.

Of these measures, the change in body temperature is most frequently probed to determine the PMI in the early postmortem period. Thermometric PMI determination was first introduced in the 19th century (*13*), and since then, many models have been developed aiming to map postmortem body temperature to time since death (*13*–*16*). The current gold standard in forensic practice is a model by Henssge (*17*–*19*), relating rectal (core) temperature to PMI. Usually presented in the form of a nomogram, this model is based on a limited set of measurements and the underlying assumption that the postmortem rectal temperature of any given human body follows a typical cooling curve. In Henssge’s model, this typical cooling curve is described by a double exponential decay, the exponents of which are derived empirically and related to the victim’s body weight, its coverage, and surface contact. While this approach is widely used, it is subject to substantial limitations. First, the underlying dataset was acquired under standardized conditions; deviations from these standard conditions necessitate the use of subjective correction factors or preclude the use of Henssge’s nomogram altogether. Moreover, use of this model requires an invasive measurement of the victim’s rectal temperature, risking contamination, and destruction of other traces. Last, the model classifies human bodies only by weight, thereby introducing a consequential thermodynamic degree of freedom: under identical circumstances, two bodies of equal weight but different stature or body composition (body fat percentage) will cool at different rates. As a result, the uncertainties of PMIs determined using Henssge’s nomogram vary broadly from ±3 to ±7 hours on a 20-hour time scale.

More rigorous, i.e., nonsubjective, approaches to thermometric early PMI determination have been developed (*20*–*24*). In these approaches, the thermodynamic processes governing body cooling are modeled using numerical methods, e.g., finite elements. While these efforts expand the applicability of thermometric PMI determination in theory, in practice, they are subject to considerable limitations. First, some require computed tomography (CT) data, which are not commonly available in forensic case work. Second, their computational implementation necessitates highly specialized technical expertise. Last, and perhaps most importantly, none of these approaches have been validated using human bodies. Consequently, there is a clear need to develop and validate a numerical approach to early PMI reconstruction straddling the divide between model complexity and usability in forensic practice.

To address this, we developed an approach that overcomes the above limitations by using a simplified but versatile numerical (finite difference) model. Body posture, stature, and composition as well as (time-dependent) environmental variables such as contact surface, (partial) submersion in water, and (partial) coverage by clothes are all readily integrated in our model (including environmental changes before and after discovery of the body), rendering it applicable in a wide variety of forensic cases. Furthermore, it allows computation of the body temperature at external body locations, in turn, enabling a noninvasive experimental protocol. We validated and benchmarked our approach using deceased human bodies. To this end, we recorded time-resolved skin temperature curves, at several body locations, on deceased human subjects and then evaluated the performance (accuracy and precision) of PMI reconstruction using our model. To determine the accuracy, we compared measured abdominal skin temperatures of four deceased subjects to their corresponding numerical predictions to reconstruct a wide range of PMIs, ranging from 5 to 50 hours since death. Next, we determined the variability in these PMI reconstructions resulting from uncertainty in the model input parameters. The outcome of this evaluation has important practical implications, as in forensic practice, some of these parameters will only be known within parameter-specific margins of error.

## RESULTS

### Experimental design

To reconstruct the PMI from a measured body temperature, we developed a finite difference model (see Materials and Methods and section S1) allowing time- and spatially resolved simulation of the body temperature following death. To this end, the model uses a discretized three-dimensional representation of the body and its surroundings to calculate the heat exchange between the involved materials. By repeating this calculation for consecutive time intervals, the change in body temperature can be simulated. The sum of the necessary computational repetitions required to reach a location-specific measured reference temperature then corresponds to the numerically reconstructed PMI.

### Theoretical model validation

To determine the model accuracy, we first tested our finite difference model on a simple system for which an analytical solution exists. To this end, we solved the problem of a solid sphere cooling in air using our finite difference algorithm. The resulting numerical solution was then compared to the analytical solution (*25*) at different spatial positions on the sphere as well as different points in time (Fig. 1). The numerical simulations are in close agreement with the analytical solution demonstrating the accuracy of the computational implementation.

### Experimental model validation using deceased human bodies

Next, we evaluated the capacity of our model to describe the change in skin temperature for recently deceased human subjects. To this end, we recorded skin cooling curves for four deceased subjects lying on a medical dissection table and stored at 2° to 4°C. These measurements were performed at specific body locations (abdomen, chest, forehead, and thigh) and compared with the corresponding model predictions. These results are presented as a function of PMI in Fig. 2. The change in environmental temperature following storage in the mortuary refrigeration unit is included in the computation and manifests as a change in cooling rate, i.e., as an inflection point, in all simulated temperature curves.

Figure 2A shows numerical and experimental data of a 79-year-old male weighing 64.5 kg, with a body length of 169 cm and with a calculated fat percentage of 29% wearing a shirt and a diaper. The body was additionally covered with a sheet, and the subject’s head was resting on a pillow. This case-specific coverage and surface contact were included in the simulations by assigning the thermal conductivities of the clothes and the pillow to the appropriate locations in the grid. Skin temperatures were recorded between 6 and 30 hours postmortem. All numerically derived temperature curves are in excellent agreement with the measured cooling curves. Figure 2B depicts simulated and measured data of a 60-year-old male with a body weight of 99 kg, a body length of 181 cm, and a calculated body fat percentage of 29%. Besides a long-sleeved shirt and a diaper, the subject also had long and facial hair; all coverage was included in the simulations by adjusting the thermal properties of the corresponding grid elements accordingly. Skin temperatures were measured between 21 and 43 hours after death. While simulation and experiment are in close agreement for the abdomen and the forehead, the measured temperatures of the chest and the thigh exceed the simulated ones at these locations. The subject of the experiment summarized in Fig. 2C was a 94-year-old female and weighed 39 kg with a body length of 159 cm and a computed fat percentage of 21%. Skin temperatures were measured between 24 and 45 hours postmortem during which the subject was naked. Despite the progressed cooling state upon arrival, the simulated temperatures are still in close agreement with the measured temperatures. Only the simulated temperatures of the chest and the thigh exhibit moderate deviations from the measured data. Last, Fig. 2D depicts the measured and simulated location-specific skin cooling curves of a 61-year-old female with a body weight of 87 kg, a body length of 157 cm, and an estimated fat percentage of 34%. This subject also wore a shirt and a diaper; however, no sheet or a pillow was present. Again, the clothes were incorporated by assigning the respective thermal conductivities to the appropriate locations within the grid. Thermal measurements took place at PMIs ranging from 26 to 50 hours. The simulated cooling curves are generally in good agreement with the measured curves. For three of the four measurement locations, however, small deviations are visible: While the time-dependent temperatures are underestimated for the forehead, they are overestimated for the thigh and the chest. Notwithstanding, all simulated decay rates accurately model those of the measured temperatures.

### PMI reconstruction and parameter sensitivity

Using both the measured and simulated abdominal data shown in Fig. 2, we reconstructed PMIs for each measurement time point of all subjects. By comparing these reconstructed PMIs to their corresponding true values (which ranged from 5 to 50 hours), we determined the accuracy of our method for PMI reconstruction. Figure 3A shows the numerically reconstructed PMIs as a function of the true PMI. The absolute errors of the PMI reconstructions (ΔPMI) are shown in Fig. 3B. All reconstructed PMIs lie within at least 3.2 hours of the true PMI, while the average ΔPMI is ±38 min. Moreover, 83.3% of the reconstructed PMIs deviate no more than ±1 hour from their corresponding true PMI. Next, we investigated the extent to which this error depends on the uncertainty in the model input parameters. To this end, we chose a specific set of input parameters as a starting point, yielding a reference PMI. We then systematically and sequentially varied these input parameters and compared the resulting reconstructed PMIs to this reference PMI. The results of this parameter sensitivity analysis are summarized in Fig. 4 and fig. S2. Both figures show the maximum error of the reconstructed PMI as a function of the variation in five model input parameters: initial body temperature (Fig. 4A), body fat percentage (Fig. 4B), thermal conductivity of the clothes (Fig. 4C), thermal conductivity of adipose tissue (fig. S2A), and thermal conductivity of nonadipose tissue (fig. S2B). The black stars denote the values for the reference dataset. Of these parameters, variation in the thermal conductivity of the clothes induces the largest variation in reconstructed PMI. On average, the reconstructed PMIs deviate no more than 2 hours from the reference PMI. Moreover, the deviation in reconstructed PMI induced by uncertainty in the thermal conductivity of the clothes, adipose tissue, and nonadipose tissue or the body fat percentage is independent of the environmental temperature. In contrast, the extent to which uncertainty in the initial body temperature induces deviation in the reconstructed PMI depends on the environmental temperature. Together, these results demonstrate that the accuracy of our method surpasses that of the gold standard even in cases where input parameters are known only within a margin of uncertainty.

## DISCUSSION

In this study, we developed and validated an advanced approach to describe postmortem body cooling. Computational robustness and numerical accuracy of this approach were established by solving the heat exchange problem for a simple geometry both analytically and numerically and demonstrating close agreement of these solutions. Overall, we found the predicted cooling curves of (partially) clothed human bodies to be in close agreement with the measured temperatures (mean deviation of 1°C) of subjects spanning large ranges in age (60 to 94 years), weight (39 to 99 kg), body length (157 to 181 cm), and fat percentage (21 to 34%) showcasing the broad applicability of the model. To demonstrate the feasibility of PMI reconstruction, we used our model in conjunction with measured abdominal temperature curves of four subjects to reconstruct PMIs ranging from 5 to 50 hours. We found the maximum and average error of these reconstructed PMIs to be as low as ±3.2 hours and ±38 min, respectively. Moreover, 83.3% of the reconstructed PMIs deviate no more than ±1 hour from their corresponding true PMIs. Parameter sensitivity analysis revealed the extent to which these errors depend on uncertainties in the studied model input parameters. Errors of the reconstructed PMIs induced by parameter-specific uncertainties remain within ±2.5 hours for true PMIs ranging from 6 to 40 hours. While the errors induced by uncertainty in the thermal conductivity of the clothes and the body fat percentage are independent of the environmental temperature, the effect of uncertainty in the initial body temperature clearly is not. This latter finding is in agreement with the results of a similar study (*26*). The initial body temperature determines the temperature gradient between the body and its surroundings; this gradient, in turn, is the fundamental property underlying and driving the heat exchange process. Therefore, uncertainty in the initial body temperature plays a bigger role at higher environmental temperatures. Notwithstanding, these results represent a notable improvement over the current gold standard (Henssge’s nomogram) where uncertainties range from ±3 to ±7 hours.

In some of the subjects, there is an observable discrepancy between the simulated and measured temperatures of the chest. The deviation of the temperature of the chest may be a result of placing the temperature sensor on the sternum. Our model does not include the location of bones within the body. This, in turn, may lead to measured temperatures exceeding simulated ones: The higher thermal conductivity of bone will transport core body heat to the skin more efficiently than accounted for in the model (see fig. S1). This deviation can therefore easily be avoided in the future by placing the sensor next to the sternum. While we successfully validated our approach on cooling human bodies, all experimental work was conducted in a standardized environment (hospital morgue), significantly reducing the variability of many environmental parameters. Consequently, performance evaluation of the approach in forensic fieldwork is paramount in determining its added value to forensic practice. The parameter most likely to vary considerably in most cases will be the posture of the victim. Currently, our approach implements the body in a straight configuration (see Fig. 5). The required individualization of the virtual body posture could be addressed using photogrammetric image processing techniques, e.g., Structure from Motion (SfM). SfM allows noncontact measurement of the three-dimensional shapes of objects from two-dimensional images (*27*); this information (the dimensions and the posture of the body) could then serve to render the virtual body in a straightforward way. Furthermore, spatial coregistration of measured and simulated skin temperatures can be achieved through the addition of coded imaging targets, increasing modeling accuracy and potentially allowing the integration of thermal imaging. Similarly useful to the individualization of the virtual body could be the inclusion of postmortem CT scans (*22*): Besides improving the assignment of the different tissue types (*28*), this tomographic information would also allow both the inclusion of cavities filled with air (or other materials) and the locations of bones in the model, the latter addressing issues such as that of the sternum mentioned earlier. While integration of these data could potentially increase model accuracy, it may hamper practical applicability of the approach by the same token: The increased model complexity would manifest as a higher computational workload, larger datasets, and decreased ease of use. However, postmortem CT scans are increasingly performed in forensic practice, rendering the required information available for a larger number of cases.

Improving the error margin in our PMI estimates necessitates the reduction of the uncertainty in the model input parameters. It would therefore also be desirable to measure the thermal properties, e.g., thermal conductivity, of the materials, which are in contact with the body, and hence relevant to the heat exchange problem, directly at the crime scene.

Last, another varying environmental parameter is the ambient temperature. An unknown ambient temperature significantly increases the uncertainty in the time of death estimate (*23*). However, in many situations, this information may be available as thermostat or meteorological data, in which case they could be easily included in the thermodynamic computations. Moreover, a general strength of our model is that multiple scenarios can be simulated to reconstruct a range of possible PMIs, the maximum and minimum of which can be reported as the most likely time frame for the PMI. While other numerical descriptions of the postmortem body cooling process exist, they lack either geometric accuracy (*20*) or experimental validation using deceased human bodies (*20*, *24*). Moreover, currently, none of these approaches include modeling of (partial) clothing, (partial) contact with surfaces, or submersion in water. To our knowledge, we are the first to validate a numerical description of postmortem skin cooling with realistic and noninvasive measurements on deceased human bodies. Last, multiple (simultaneous) skin temperature measurements may reduce uncertainty in the PMI calculation by increasing the number of independent measurements. Together, these aspects render this study a considerable advance in the pursuit of a widely applicable, and hence forensically relevant, method for temperature-based PMI estimation.

## MATERIALS AND METHODS

### Modeling the body and its environment

Thermometric PMI reconstruction using the numerical approach outlined in the theory section (see section S1) requires modeling (discretization) of the body and its environment. Here, a discretized three-dimensional representation of the body and its surroundings is obtained as follows. First, the individual body parts are approximated as cones (e.g., arms and legs), ellipsoids (e.g., head), and cylinders (e.g., neck and torso), the proportions of which are dictated by standardized anatomical measurements (e.g., length of the arms and legs and the circumference of the head, torso, upper arms, and wrists). Second, these individual body parts are then assembled to form the entire body. Third, this simplified model of the body is placed within an isotropic cubic mesh (which is generated automatically based on the body dimensions) of desired cube size (in this study, 1 cm^{3}) where it then serves to determine material (i.e., thermal properties) assignments on a cube-by-cube basis (see Fig. 5A). Similarly, to obtain spatially resolved temperature data for different time points, it is necessary to discretize the time before equilibrium into finite time intervals Δ*t*, in our case 60 s (matching the sampling period of our temperature measurements). This, in turn, determines the number of computational iterations needed to simulate the heat exchange process: e.g., for a period of 24 hours, a Δ*t* of 60 s corresponds to 1440 iterations.

Besides modeling the geometry of the body and its surroundings, it is crucial that the material composition of the computational representation (i.e., the cube-wise assignment of the thermal properties) closely resembles reality to ensure accurate simulation of the body temperature change. Important factors therefore include body fat percentage, coverage by clothes, and contact area with other surfaces, such as the floor or water. The body’s adipose tissue is modeled as an outer layer surrounding the individual body parts. To accurately model the layer thickness, we estimate the fat percentage using the U.S. Navy circumference method (see Fig. 5A) (*29*).

The thermal conductivities of the nonadipose tissue, adipose tissue, and the clothes (in our case cotton) were set to 0.55, 0.2, and 0.03 Wm^{−1} K^{−1}, respectively (*30*, *31*). For the emissivity ϵ, we used the emissivity of human skin, namely, 0.96 (*32*), while for the characteristic length of the convective heat transfer, we chose the characteristic length of a cylinder approximated by the width and depth of the torso of the virtual body (*33*). Last, the specific heat capacities of nonadipose tissue, adipose tissue, skin, and clothing were chosen as 4.5, 1.96, 3.77, and 4 JK^{−1} cm^{−3}, respectively (*30*, *34*–*36*). All simulations and data analyses were carried out using custom-made scripts written in MATLAB (The MathWorks Inc., Natick, Massachusetts, USA) and executed on a laptop with 8 GB RAM and an Intel Core i5-8250U CPU operating at 1.6 GHz. Typically encountered runtimes of our simulations are on the order of 10 s to 20 s.

### Theoretical model validation

The accuracy of the finite difference model was validated by numerically solving the problem of a solid sphere (radius = 30 cm) of nonadipose tissue (thermal conductivity = 0.55 Wm^{−1} K^{−1}) in air over a period of 50 hours. The spatial and temporal changes of the temperature distribution of the sphere were then compared to the corresponding analytical solutions (*25*).

### Experimental model validation using deceased human bodies

Our finite difference approach yields spatially and temporally resolved computations of the postmortem body temperature (see Fig. 5B). This, in turn, allows for the temperature measurements to be performed on the surface (i.e., the skin) of the body rather than rectally as required for the use of Henssge’s nomogram. To validate this approach, we conducted postmortem skin temperature measurements on deceased subjects available through the body donation program of the Department of Medical Biology, Section Clinical Anatomy and Embryology, of the Amsterdam University Medical Centers [location Academic Medical Center (AMC)] in The Netherlands. The donation of these bodies to science occurred in accordance with Dutch legislation and the regulations of the medical ethical committee of the Amsterdam UMC, location AMC. Bodies were included if they were suitably warm to ensure a measurable change in skin temperature. All measurements were performed in the AMC morgue using small (⌀ = 16 mm) contact temperature sensors (DS1922L iButtons, Thermochron, USA) attached to the bodies using adhesive tape. Using this approach, skin cooling datasets were gathered from four deceased subjects, two males and two females, aged between 60 and 94 years. The body weights, body lengths, and calculated fat percentages of the subjects ranged from 39 to 99 kg, 157 and 181 cm, and from 21 to 34%, respectively. The amount of clothing and surface contact with insulating materials as well as the PMI at the beginning of the measurements varied between experiments. As all of these are input parameters in our model, their variation between experiments was accounted for by adjusting them accordingly in the respective calculations.

### PMI reconstruction

Reconstruction of the time since death, using our model, comprises computation of the heat exchange until a measured location-specific reference temperature is reached. The sum of the necessary computational time steps, i.e., *n**Δ*t*, required to reach the measured temperature then corresponds to the numerically reconstructed PMI (see section S1). The accuracy of this numerical PMI reconstruction was evaluated using abdominal skin temperatures of four deceased subjects comprising 96 measurement time points, i.e., true PMIs, between 5 and 50 hours with the corresponding measured abdominal reference temperatures ranging from 28° to 5°C.

### Parameter sensitivity analysis

Our computational approach requires the values of several physical quantities (model parameters): the initial body temperature, the environmental temperature, the body fat percentage, and, if applicable, the thermal conductivity of the clothes and the floor, the floor temperature, as well as the flow speed of the surrounding medium (e.g., air). The extent to which uncertainty in these input parameters generates variation in our PMI predictions was evaluated by simulating the cooling of a body for 13 different sets of model input parameters at two different environmental temperatures (10° and 20°C) yielding a total of 26 set-specific simulations. These 26 parameter sets were generated as follows: For both environmental temperatures, an original set of input parameters was chosen (see Table 1). Subsequently, we assigned physiologically and physically reasonable parameter ranges for three of the seven parameters, namely, initial body temperature, thermal conductivity of the clothes, and fat percentage of the body (see Table 1). This particular choice of parameters to investigate was motivated by their likelihood of being unknown and/or exhibit most variation in forensic casework. Next, for each of these three parameters, four (equidistant) values were chosen within their respective parameter ranges. Together with the original set of input parameters, this yields 13 separate sets of model parameters (see table S1) for each environmental temperature, i.e., 26 separate cooling simulations in total. For each of the two environmental temperatures, skin temperatures (of the abdomen, chest, upper arm, and thigh) simulated using the original parameter set served as a reference in the calculation of the deviation in predicted PMIs. The reference temperature ranged from 34.7° to 10.6°C (environmental temperature = 10°C) and 35.6° to 20.4°C (environmental temperature = 20°C) corresponding to true PMIs between 6 and 40 hours. Deviations from these true PMIs were then calculated for each of the four body locations for all 12 parameter sets for both environmental temperatures. Last, we determined the maximum deviation in reconstructed PMI (over all reference temperatures, i.e., “reference PMIs”, and body locations) for every parameter set for both environmental temperatures yielding a total of 24 data points. In addition, we investigated the parameter sensitivity of our model with respect to two more model parameters: the thermal conductivities of adipose and nonadipose tissue. Here, we chose two equidistant values on either side of the original parameter value within their respective parameter ranges (see Table 1) and calculated the maximum deviation in PMI at two environmental temperatures (10° and 20°C). The results are summarized in fig. S2.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/22/eaba4243/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank M. Clerkx, R.-J. Oostra, H. Boesveld, J. Woertman, M. van den Born, and B. F. L. Oude Grotebevelsburg for help in carrying out measurements and improving the algorithm.

**Funding:**Ministerie van Veiligheid en Justitie, Innovatieproject ronde 2019 “Therminus”.

**Author contributions:**L.S.W. developed software, performed measurements, and analyzed data. R.J.M.H., G.J.E., and H.V. performed measurements. H.J.J.H. partly developed the theoretical framework. S.S. developed software. M.C.G.A. performed measurements, partly developed the theoretical framework, and conceived and supervised the project. All authors contributed to planning, design of experiments, discussion of results, and writing of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).