Background information on numerically solving partial differential equations and on FE simulation can be found e.g. in [16].

For one exemplary case of human corpse with a given CT-scan, two different methods were used to generate FE models. In one of the methods, the FE Model (CTM) was developed from the segmentation of a CT scanned human body and in the other, the human body was manually approximated by hexahedral elements (see e.g. [2, 17]). Two variants of this manually generated model were constructed: One model (CM) is floating freely in air and a copy (CMS) of CM is laying firmly on a wet soil substrate. Detailed simulations of physical heat transfer processes generated the rectal temperature curves T(t). The corpse cooling was computed in all of the FE models based on the well known heat transfer Eq. (1), a partial differential equation where c is the specific heat capacity, ρ the mass density, and κ the heat conductivity of the tissue:

$$c\rho \fracT = \nabla \cdot (\kappa \nabla T)$$

(1)

Heat transfer from the body to the environment across the skin due to convection and surface-to-ambient radiation was captured by a Robin boundary condition with effective heat transfer coefficient γ as given in Eqs. (2), (3) (see e.g. [7]):

$$\gamma = h + 4 \epsilon \sigma _}^$$

(2)

$$^\kappa \nabla T = \gamma (_-T)$$

(3)

Here, TA is the ambient temperature, h is the heat transfer coefficient, ε is the emissivity and σ is the Stefan-Boltzmann constant. In the CTM developed from the CT scan of a human body, the initial temperature field T0 defined at time t = 0 satisfies Pennes’ Bio-Heat-Transfer-Equation (BHTE) [18]

$$-div (\kappa \nabla _)+__w(_-_)=0,$$

(4)

where ρbcb is the heat capacity of blood, w is tissue perfusion, and Tcore is the body core temperature.

In the manually developed FE model CM, the initial temperature field T0 is defined with a gradient between core and outside elements as in [2] referring to physiology literature.

For reliable comparison between the different models, the equivalent effective heat transfer coefficient was applied to the CTM corresponding to the convection and radiation terms applied on CMS and CM. The model curves T(t) were sampled at the designated nodes C and SPRk (k = 1,…6) corresponding to the anatomical RTM locations (RTML). The point C was placed at the intended measurement location in the respective FE model, whereas the six points SPRk (k = 1,…,6) lay on the vertices of an octahedron with radius R. The SPRk were defined to represent possible locations of a misplaced temperature sensor. If no confusion can arise about the value of R the SPRk are abbreviated SPi.

Our results consist of maximum deviations DMAXM comparing two cooling curves computed at two locations of RTM’s in a specified FE model M = CTM, CM, CMS. The distance DMAXM is defined between the cooling curves TC(t) and TSPi(t) at the center point C of a measurement point octahedron (see Fig. 4) of radius R and at its vertices SPRk (k = 1, …, 6). We also computed differences DMAXM1,M2: = DMAXM1 – DMAXM2 of those maximum deviations for two pairs (CM-CMS and CM-CTM) of different FE models (M1, M2) for all C – SPRi pairs. All computations were performed for different constant ambient temperatures TA and for different octahedra radii R. We also performed all of the computations for each of four Q-ranges, where the first three are Henßges Q-ranges (see e.g. [15]).

All of the modeling work done on the CT meshed FE model CTM was performed at the Zuse Institute Berlin (ZIB). The CT scan of a male corpse of length L = 1.74 m and weight M = 62 kg was segmented in the software AMIRA and the segmented data was converted into an FE model, comprising of 961,234 tetrahedral elements, via the FE-program Kaskade (see e.g. [19]), an in-house code of ZIB. As depicted in Fig. 1, the feet of the corpse were not included in the model and the weight, size and volume of the whole model were scaled according to the scanned corpse dimensions. The segmentation of different body parts was carried out using differences in density and location of different tissue types like bone, fat, bladder, kidney, abdomen, liver, heart, lungs and muscle as seen in Table 1 in the Appendix.. In CTM, skin tissue was not segmented separately from the subcutaneous fat tissue since skin tissue was assumed to have a negligible effect on the TDE. Hence, the skin material properties were not considered in CTM.

a: CT Meshed male corpse model with magnified mesh detail. b: Right half view of CTM model with magnified measurement location showing measurement points C, SPk (k = 1,…,6) in three different octahedra of radii R = 0.5 cm, 1 cm, 2 cm. For the octahedra see Fig. 4. Only the contour lines of the octahedra are shown for clarity

Our coarse meshed FE model CM was already described in several publications [2, 3, 20, 21]. Its mesh was manually constructed using the FE preprocessor Mentat. It consists of approx. 12,000 nodes in approx. 8300 hexahedral elements with trilinear shape functions. The model contains numerous compartments (see Fig. 2) standing for distinct anatomical structures of the body.

a: Coarse model lying on a support structure (CMS), b: Coarse model without support structure (CM), c: Detailed view of measurement points C, SP1,….,SP6 along different radii R = 0.5 cm, 1 cm, 2 cm. For the octahedra see Fig. 4. Only the contour lines of the octahedra are shown for clarity

Our standard FE model of length L = 1.64 m and mass M = 64 kg was scaled geometrically to a length of L’ = 1.74 m and weight of M’ = 62 kg, corresponding to the male reference body of our study. The scaling factors are defined as in the Eqs. (5) and (6) where k1 is linear 1-dim scaling along the body length axis and k2 is linear 2-dim dilation in the transverse plane [2].

$$}_=\sqrt^}\cdot L}^}}}$$

(6)

Generation and scaling of FE mesh and post-processing were performed in MSC’s pre- and post-processing tool Mentat.

Further, to investigate the TDE sensitivity on more complex boundary conditions favoring thermal energy transfer from the body core in dorsal direction, a supporting structure or floor is modeled with an hexagonal FE mesh as shown in Fig. 2a. The FE model lies on its back on the support structure defined, presumably causing large differences between C-temperature and temperature at SP4 or SP3, i.e. along y direction. The support structure is defined with the thermal properties of wet soil such as thermal conductivity c = 2 W/m2K, specific heat κ = 2200 J/kgK, density ρ = 1900 kg/m3 [22] and emissivity ε = 0.95. The coarse model laying on its back on the support structure is abbreviated CMS, whereas the coarse model without support structure is named CM (see Fig. 2c). CMS and CM were scaled to the same height and weight as CTM. This will eliminate the influence of weight on the body cooling.

TDE estimations are usually based on a single RTM. Hence, the measurement position is subject to variations due to different anatomies, thermometer angles and insertion depth. An approximate ideal measurement position can be specified from anatomy including forensic knowledge on measurement locations used in practice. The central position C for our measurements was determined in the pars ampullaris of the rectum near the incisura transversalis, which is the passage from pars ampullaris to the pars sacralis of the rectum (Fig. 3).

Anatomical sketch depicting ideal measurement location [23]

Additional measurement positions SPRk (k = 1, …,6) for our sensitivity studies were chosen in spatial x, y and z direction at distances (octahedral radii) R = 0.5 cm, 1 cm, 2 cm from the central measurement position C. For each fixed distance R six additional measurement locations SPRk were established. The six additional measurement positions SPRk form an octahedron with the central point C in its middle. The schematic sketch representing the additional positions in the model was illustrated in Fig. 4.

Central position C and six additional measurement positions (SPRk =) SPk (k = 1,…,6) in X, Y and Z direction in CTM, CMS, CM model at different radii R = 0.5 cm, 1 cm, 2 cm from center position C

The anatomical position of additional measurement points in CTM are described as follows: SP1 lies in the center of abdomen surrounded by abdomen tissues, SP3 lies nearby bone, SP4 lies in the center of abdomen, SP6 lies very near to abdomen. Figure 1b shows the anatomical positions of measurement points in CTM. Similarly, the anatomical position of additional measurement points in CMS and CM are depicted in Fig. 2c, which helps to identify the position of additional measurement points. SP1 lies lateral from C in positive X direction. Due to CM’s left–right symmetry the temperature value at SP2 is equal to the temperature at SP1. The points SP2, SP3, SP5 lie in the gastrointestinal tract. SP6 and SP4 lie on the border between gastrointestinal tract and pelvis.

The CT-generated model CTM was imported from the FE research code Kaskade 7 [19] into the commercial MSC-Marc Mentat FE system and the material parameters are defined as given in Table 1 of [2] (see Table 1 in the Appendix). The same material properties were used also for CMS and CM models. The initial temperature field T0(r) at the time of death on the body was calculated via the Pennes’ Bio-Heat-Transfer-Equation BHTE on the original CTM model [18]. The field T0(r) was approximated by 10 discrete node sets corresponding to 10 discrete temperatures in the converted version of CTM. In the CMS and CM models, the initial temperature was defined according to [2] where it was taken from physiological textbooks. The simulation is carried out for three different constant ambient temperatures TA = 5° C, 15° C, 25° C. The convection and radiation parameters in CTM were defined by the effective heat transfer coefficient γ = 4∙ε∙σ TA3 [7] and the calculated values in increasing order of TA = 5° C, 15° C, 25° C are γ = 7.93 W/m2K, 8.45 W/m2K, 9 W/m2K. In CMS and CM, a convection coefficient of h = 3.3 W/m2 K is applied [2]. Radiation view factors in the CMS and CM models were computed by Monte-Carlo-simulation [2]. Although the method of defining boundary conditions differs between CTM and CM models, it was found that the linearization of the radiation term applied in CTM has negligible effect on the cooling curve [7]. A sparse direct algorithm in MSC-Marc 2020 and adaptive time stepping with maximum allowed temperature change of 1°K and 0.1 s initial time step were used.

TDE deviation from the central position C to any position SPRk is quantified by the maximum distance DMAX between the respective cooling curves, which is defined in this section. Let the time interval [a, b] be given during body cooling which starts at 0 h and let a = 1 h be the left point and b = 45 h be the right point of our interval. Let further T1(t) and T2(t) be two cooling curves defined on the interval [0, b]. Evaluation of the cooling curve distances started at a > 0 h for reasons of numerical stability. The ambient temperature TA is assumed constant in space and time outside the body. The reference curve TR(t) for two cooling curves T1(t) and T2(t) is constructed as the pointwise mean of T1(t) and T2(t):

$$\forall t\in \left[a,b\right]: _\left(t\right):=(_\left(t\right)+_(t))/2$$

(7)

Let T0: = TR(0) be the initial temperature at the center point C. Now the reference curve TR(t) can be normalized to the function Q(t) taking values in the real number’s interval [0, 1]:

$$\forall t\in \left[a,b\right]: Q\left(t\right):=(_\left(t\right)-_)/(_-_)$$

(8)

This makes it possible to define the Q-intervals Q1 to Q4 (where Q1, Q2, Q3 correspond to Henßge’s [15] normed temperature ranges for tolerance radii) and their left and right boundaries q0, …, q4 by:

$$_:=1, _:=0.5, _:=0.3, _:=0.2, _:=0.1$$

(9)

$$\forall 1\le n\le 4: _:=]_,_]$$

(10)

Assuming TR(t) to be strictly monotonically decreasing, we can uniquely map the Q-interval boundaries qn to the respective temperature values TR(4-n).

$$\forall 0\le n\le 4: _:=_\cdot \left(_-_\right)+_$$

(11)

Since the two cooling curves T1 and T2 are monotonically decreasing it is possible to define their inverse curves t1, t2, called the (absolute) time since death (TSD) estimation result curves w.r.t. the cooling curves T1, T2 on the domain [TMIN,TMAX] which is the range of the reference curve TR on the temperature axis by:

$$_:[_,_]\to IR:(T\to _\left(T\right):=_^\left(T\right))$$

(12)

$$_:[_,_]\to IR:(T\to _\left(T\right):=_^\left(T\right))$$

(13)

The time functions t1, t2 are frequently concatenated with the reference curve TR, thus giving t1(TR(t)), t2(TR(t)) in the following definitions. For short those concatenations will be abbreviated t1(t), t2(t). If the point t in time is clear, we will even write t1, t2. We will now quantify the distance between the cooling curves T1 and T2, which are actually constructed as distances between the inverse curves t1 and t2. The reason for this construction lies in our primal interest in time differences because our research’s number one target is time since death. All of the temperature curve distance measures are constructed using the (time oriented!) maximum distance DMAX(T1,T2) between two cooling curves T1 and T2, which is shown here for continuous functions first. The maximum domain of the reference curve TR on the time axis is the interval [tMIN,tMAX] which is taken therefore as the domain for maximum-finding in our definition:

$$_\left(_,_\right):=ma__,_] }\left|_\left(_\left(t\right)\right)- _(_(t))\right|$$

(14)

We will now consider the equivalent of the definition (14) in terms of real measured temperature curves Ti(t). Each curve Ti(t) is represented as a finite series of real measurement values (Ti1, …, TiN) on a finite (regular) grid of time values (t1, …, tN) which is the same for both of the curves. The first point t1 of the time grid is identical to the starting time of cooling computation, while the last point tN marks the end of the cooling computation interval. To provide a joint domain of definition for both our inverted temperature curves TR−1 we define:

$$_:=min\left\_^,_^\right\}$$

(15)

$$_:=max\left\_^,_^\right\}$$

(16)

Now we have constructed the range of our reference curve TR(t). The limits of its domain of definition [tMIN,tMAX] were computed:

Let K be a natural number and let (t1, …, tK) be a regular time grid on the interval [tMIN,tMAX] with t1: = tMIN and tK: = tMAX. This provides the equidistant sampling points for the final definition of the functions ti(TR(t)). Let further be (T1, …, TK) the corresponding temperature grid with T1 ≤ … ≤ TK and Tk: = TR(tk) for all k = 1, …, K and TK: = TMAX. The time grid’s width Δt is:

$$\Delta t:=\frac_-_}$$

(19)

Let further be tqn the inverse image of TRn under the reference function TR for all n = 0, …, 4:

$$\forall 0\le n\le 4: _:=_^\left(_\right)$$

(20)

The five points tqn on the time scale constitute the endpoints of four time intervals tQ1, …, tQ4:

$$\forall 1\le n\le 4: _:=]_,_]$$

(21)

We will redefine our measure (14) quantifying the distance between the cooling curves T1 and T2 in terms of the real samples (t1, …, tK) and (Ti1, …, TiK). The distance measure is the distance DMAX (T1, T2). It is defined in a global version on the whole time interval [tMIN, tMAX] and in four local ones residing on one of the intersections tQj ∩ [tMIN, tMAX] on the time axis each. For all 1 ≤ i ≤ 4 the number of tk lying in tQj ∩ [tMIN, tMAX] is denoted by Ki.

The (absolute global) maximum distance DMAX(T1, T2) is defined as:

$$_\left(_,_\right):=_|_\left(_\left(^\right)\right)-_\left(_\left(^\right)\right)|$$

(22)

while the (absolute) Qi-local maximum distance DMAX,Qi(T1, T2) is defined for i = 1, …, 4:

$$__}\left(_,_\right):=__, ^\in _ \cap [_, _]}|_\left(_\left(^\right)\right)-_\left(_\left(^\right)\right)|$$

(23)