Finite-difference time domain (FDTD) modeling of ground penetrating radar pulse energy for locating burial sites

Analysis of the finite-difference time domain (FDTD) numerical simulation of ground penetrating radar (GPR) measurement for locating burial sites is described in this paper. Effective, efficient, and reliability interpretation of GPR field data obtained from clandestine sites is very crucial in forensic investigations. The main goal of the study is the prediction of the change in the interaction of the electromagnetic incident on changes in buried bodies with time. In order to achieve this, the research involves the modeling of the GPR electromagnetic pulse energy responses to simulated changes in buried body with time with a view to understand what the results of real field measurement will give. The field measurements were conducted with GPR system manufactured by Mala Geoscience with antennae frequency of 500 MHz, 250 MHz, and 100 MHz. Responses from both synthetic and field radargrams depict the target was intercepted at same time (approximately 25 ns). The results have demonstrated that FDTD modeling is an important tool for enhancing the reliability of GPR data interpretation particularly for forensic study.


Introduction
Organic bodies after burial undergo decomposition with time. The decomposition is largely dependent on so many factors which hamper the simple detection of the point of burial that may be crucial in some circumstances. Locating buried bodies is primarily related to forensic investigations. Forensic investigation is the gathering and analysis of all crime-related physical evidence in order to come to a conclusion about a suspect. Investigators will look at blood, fluid, fingerprints, residue, hard drives, computers, or other technology to establish how a crime took place.
The overall forensic investigation is vital to the mission of truth seeking, because the forensic evidence does not lie. Moreover, gathering and examining forensic evidence is the only way to ensure that false information is exposed and ultimately kept from the gates of justice (Turvey and Crowther 2017).
Clandestine graves contain buried human remains not intended to be found and are the focus of forensic or criminal investigations led by police, increasingly drawing upon the expertise of anthropologists and other professionals. Clandestine graves may be a result of murder case, war and natural disasters such as landslides, avalanches, etc.
Locating hidden and buried objects has long been a subject of interest to archeologists seeking buried sites, tombs of historical objects. Forensic investigations are concerned with locating, identifying, collecting, and cataloguing physical evidence for the purpose of presenting it in court. One aspect of forensic investigations concerns locating clandestine evidence which is often concealed in the subsurface. This task is typically guided by information provided by informants (pathological liars), witnesses, psychics, and sometimes suspects (Davenport 2001). The 1 3 various methods of locating secret buried bodies as found in literature include: human eye witness, use of cadaver dogs, chemical analysis method, remote sensing (LiDAR, aerial photography, satellite imagery, thermovision imaging, sonar method, and geoscientific methods).
Eye witness individuals account may be valuable to locating clandestine grave if given information is truthful. Cadaver dogs are dogs trained to use its sense of smell to locate remains, usually human. The science of cadaver dogs is a fascinating combination of biology and singleminded training, said Cat Warren, the author of "What the Dog Knows: Scent, Science, and the Amazing Ways Dogs Perceive the World" (Warren 2015). Soil and solid materials sampled from suspected scene and chemically analyzed may give useful clues to locating clandestine graves (Vass et al. 2004;David et al. 2008). Brilis et al. (2000) and Grip et al. (2000) provided a comprehensive overview of remote sensing methods in forensic investigations, which includes aerial photography, topographic mapping, satellite imagery, and global positioning systems (GNSS) applications. They detail the usefulness of early forensic balloon unbiased photographs, taken just after the 1906 San Francisco earthquake, to detail the location and extent of damage caused both by the earthquake itself and by non-earthquake insurance homeowners setting fire to their own homes. Ultraviolet (UV) photography is increasingly being used in forensic investigations to map the maturity of vegetation-on the premise that any vegetation growing over recently buried material will be younger than the surrounding plants ( Ruffell and McKinley 2008). Satellite imagery datasets are many and varied; hyperspectral imagery, for example, has proven useful to determine locations of individual clandestine graves and mass graves in areas of rapid vegetation growth, changing land-use and humid climates in Colombia (Ruffell and McKinley 2008;Equitas 2010). Light detection and ranging (LiDAR) airborne data capture systems can be used to very rapidly collect laser scan threedimensional points in most weather conditions. This has a significant advantage over aerial or ground photographs (Ruffell and McKinley 2008;Mallet and Bretar 2009).
Thermal imaging equipment may be aircraft mounted (Dickinson 1976), vehicle mounted or handheld (Davenport 2001;Statheropoulos et al. 2011). Handheld devices have been used for landmine detection, and objects can be detected even if temperature differences between the target and background levels are less than 1 °C (Deans et al. 2006). Clandestine graves are thought to be most easily detected with thermal techniques during the first few weeks of burial (Dickinson 1976). It has also been found that buried animal was detectable for up to 17 days post-burial, although this time frame will be influenced by the local geology, ground water level, and soil type(s); after this time, it gets progressively more difficult to pinpoint burial locations using this technique (Larson et al. 2011).
Knowledge of the local geology and soil types can be critical in forensic investigations. Williams and Aitkenhead (1991) discuss how the lack of knowledge of the local geology and the poor understanding of the overlying landfill geochemistry in Loscoe, Derbyshire (UK), culminated in a methane gas explosion and subsequent destruction of a domestic house in 1986 seriously injuring three occupants. Lee (2004) also described a house explosion in 2004, where on investigation it was found that rather than a local gas leak or malicious damage, methane from underlying coal seams were the cause of the explosion.
Although all the above range of techniques has been found successful under some circumstances, yet they are not perfect methods for the location of clandestine graves. This may be largely attributed to the fact that they are not reliable, difficult to be carried out, and destructive (jeopardize evidence), not 100% accurate. Moreover, resultant searches may be time-consuming and frustrating to understaffed departments. No doubt, any nondestructive method used to reduce the time spent on searches and excavations and to increase the probability of locating physical evidence are of prime interest to the law enforcement community. Essentially, geophysical method involves the application of physics in the study of the earth crust (Powell 2004). One of the geophysical methods used for near-surface investigation is ground penetrating radar-GPR (for instance, Giang et al. (2013) reported the efficiency application of GPR in near-surface pollution study.) GPR application for the study and monitoring of landslides in part of Carpathian flysch is also found in the literature (Pilecki et al. 2007). GPR, an electromagnetic geophysical method, which is nondestructive, is currently being considered with a view to circumvent the aforementioned limitations of the conventional methods of locating clandestine graves. Similar approach to the one used in this research study is recorded in the literature. Widodo et al. (2016) utilized GPR for locating buried bodies due to landslide avalanches triggered by frequent earthquake in part of Indonesia. The study distinguished detection of new buried bodies from the old buried bodies with variation in signal amplitude. Moreover, Solla et al. (2012) in their work experimented the characterization of simulated GPR response of crime scene buried items such as bone remains, gun, and drug catches. The method adopted was the FDTD method of numerical simulation, and the results showed that GPR method may be effective in actual forensic investigation. Similarly, Fernández-Álvarez et al. (2016) have demonstrated the efficiency of GPR application to locating mass buried site during Spanish civil war between 1936 and 1939. In most of these previous similar researches, no one considered the GPR response to changes in the buried organic bodies with time. It is on this premise and other considerations that this research attempt was made. In the work of Carcione et al. (2017), a pseudo-spectral method of modeling was used to locate organic body.
The focus of the research is the modeling of the GPR electromagnetic pulse energy responses to simulated buried body with a view to understand what the results of real field measurement will give. It is aimed at knowing the appropriate parameters to be taken into consideration for actual field measurement.

Method of study
GPR near-surface geophysical method was used for this study. The nondestructive nature of the method which is an advantage to forensic investigations has informed its choice. A comprehensive introduction to the GPR method can be found, among others in the works (Daniels 2004;Jol 2008;Karczewski et al. 2011). Basically, GPR is a geophysical method that uses radar pulses to image the subsurface. It is a nondestructive method that engages electromagnetic radiation in radio spectrum to detect the reflected signals from the subsurface. GPR uses high-frequency (usually polarized) electromagnetic waves, in the range 10 MHz to 2.6 GHz. A GPR transmitter emits electromagnetic energy into the ground. When the energy encounters a buried object or a boundary between materials having different permittivities (dielectric constants), it may be reflected, refracted, or scattered back to the surface. A receiving antenna can then record the variations in the return signal.
In this study, buried body (cadaver) within the near surface has been modeled and its responses to GPR signals were also simulated using an open source computer software GprMax2D developed by Antonis Giannopoulos, University of Edinburgh, UK (Giannopoulos 2005;Warren et al. 2016). The program was written in C language and is based on finite-difference time domain (FDTD) method of approximation. FDTD is based on Maxwell's equations which describe the behavior and effect of electromagnetism. It is the method for simulating computational electromagnetism and considered the simplest and most efficient way to model the effects of electromagnetism on a certain material or object. In the FDTD approach, both space and time are divided into discrete segments. Space is segmented into box-shaped cells, which are small compared to the wavelength (Luebbers et al. 1993;Taflove and Umashankar 1989;Taflove and Hagness 2000). The electric fields are located on the edges of the cube, and the magnetic fields are positioned on the faces as shown in Fig. 1. In a simple explanation, FDTD computation is performed in three steps: discretization of the whole area into small cells (definition Δx, Δy, and Δz) in the Yee cell (Yee 1966); definition electrical properties (permittivity and conductivity) of the model in mesh grid; and substitution of partial derivatives by differential quotient in Maxwell's equations.
Most modeling and simulation of GPR pulse energy which analyze material properties of subsurface geological media are defined by values of electrical conductivity (σ), dielectric permittivity (ε), and magnetic permeability (µ). For most soils, magnetic permeability is equal to 1. The choice of these parameters is based on standard range of values as found in literature (Jol 2008;Daniels 2004).
The hypothetical model is essentially four layers. It has zone of buried body at the center which is based on forward modeling techniques with numerical parameters (as shown in Table 1). However some assumptions were made. The zones represent the affected parts by body decomposition within the model. Model 1 represents the subsurface layers, while model 2 has the subsurface layer and the pit for the burial of the body. In models 3-6, the subsurface layers, the body, and zones depicting decomposition influence on the surrounding layers were represented.
The choice of dielectric permittivity of the zones (1, 2, 3) was assumed based on the knowledge of decomposition rate of organic body which occurs in five general stages: fresh, bloat, active decay, advance decay, and dry/remain. The rate of the decomposition is dependent on the climatic conditions which are function of temperature, humidity, etc. (Vass et al. 2004).
The geometry of the model has dimension of 5 m in length and depth of 3 m with point of the burial between 2.2 and 2.8 m at the center with a depth of 1.2 m. The cadaver (organic body) was given a radius of 0.2 m. Moreover, to represent the model in a more reality, some disturbances were introduced into the model layers. For example, limestone with a random dimension (0.01-0.04 m) was added to the basal wet clay layer. Also gravel of random dimension (0.005-0.01 m) was added to the upper layers. Essentially, these additions are to make the model display heterogeneity.  (Karczewski et al. 2011) The properties of the model components are given in Table 2.
The field measurement was made on a piece of ground where pig cadaver was buried (Fig. 2a) with GPR equipment (ProEx System) manufactured by MALA GeoScience Sweden (now: ABEM/MALA). Antennae with central frequency of 100 and 200 MHz (unshielded); 250, 500, and 800 MHz (shielded) were used (Fig. 3). The measurements were deployed in the constant offset mode along a profile of 5 m perpendicular to the strike of the buried spot (Fig. 2a).
Field measurement was deployed in the short offset mode after the configuration of the system. Subsequently, the entire arrangements were mounted and data were collected by pulling the wheel antenna along the profile at a walking speed. The system generates radar impulse energy at the given frequency of its antenna. The pulse energy propagates into the media and gets scattered back at electromagnetic discontinuities mainly due to contrast in dielectric constant between the media and the buried  Modroo and Olhoeft (2004) Model no Body Zone 1 Zone 2 Zone 3 3  50 a  8 a  ------4  45  7  40  6  ----5  40  6  35  5  30  4  --6  35  5  30  4  25  3 20 2  body. Prior to the measurement, the GPR system was configured with the basic parameters which is dependent on the target. Time window was set at a range of 70-180 ns depending on the antennae frequency. A trace stacking of 16 was selected to improve S/N ratio. Numbers of sample per trace were also pegged at 516. The choice of the range of the antennae central frequency was based on the resolution and depth of the target.

Results and discussion
In this section, the results of the study are presented with explanation of the findings and the inferences made.

Models
The numerical models of the buried body with variation of parameters of the material properties are as shown in Figs. 4 and 5. They are numerical representations of the dimension and physical properties variations of the soil and the buried body. Moreover, the models (Fig. 5) depicted the changes in the buried body with time along with changes in the surrounding soil properties. These variations are represented as zones that increase with time in the model.

GPR data simulation
Results of the modeling of GPR signal responses for different simulated central frequency antennae are presented as synthetic radargrams (as shown in Fig. 6). Essentially, the image of the radargrams represented a plot of GPR signal amplitudes against time along the scanned profile. What is important in the results is the delineation of the model's layers and hidden body via the anomalies responses as shown in Fig. 6. The wave velocity was estimated based on the dielectric permittivity of the layers of the models. The waged average value is v = 0.07 m/ns. Although the zones represented as variations due to changes in the buried body with time are not well elucidated, the responses pointed out contrasts in the properties of the model. The reflection coefficient was calculated using Eq. 1. Fig. 6 Results of modeling. a Six different synthetic radargrams for six models. b Comparison of synthetic radargrams of models 3 and 6 with antenna frequency of 250 and 100 MHz where ɛ 1 , ɛ 1 -electric permitivities of different zones.
In Table 3 are the calculated values of the coefficient of reflectivity for different zones.
Results of the simulation are a pointer to better field measurement and interpretation.

Field data results
The field data after editing and processing are presented as real radargrams (Fig. 7) and interpreted on the basis of available velocity through which the electromagnetic signal propagated. Data were processed by REFLEXW software (Sandmeier 2012) using following procedures: move start time, DC Shift, dewow, gain, background removal, frequency filtering and so on. On the knowledge of the geology of the area, it was assumed that the velocity of the electromagnetic wave was v = 0.07 m/ns. The results have given insights to the buried body via the hyperbola anomaly which indicated its delineation.

Correlations of synthetic and field radargams
In order to prove the efficiency of the modeling and the numerical simulations, the synthetic and field radargrams were compared to evaluate similarities (Figs. 8 and 9). The whole essence of the correlation is to have a deeper insight on the responses of the different approaches to detecting buried organic body. Hence, hidden burial sites may be better investigated.
It can be seen (Figs. 8 and 9) that both images indicated anomalies at the same distance and depth. Thus, the positions of the buried body (pig cadaver) are discernable. Real geological model is more complicated than numerical model. In the real model, there are variations in the dielectrical parameters which may be as a result of heterogeneity of the layer. The changes in humidity also influence the response to the electromagnetic energy. The success of the modeling process lies largely on the appropriate choice of modeling parameters and careful study of the material properties of the target.

Conclusions
The clutter of hyperbolic diffractions displayed in the synthetic models radargrams may be attributed to the reflection from the disturbance (limestone/gravel) introduced and not from the targeted objects. Resolution analysis of both synthetic and field results shows that antenna frequency 250 MHz may be the appropriate to be deployed for buried organic body at the depth of investigation as illustrated in this study. As shown in Fig. 9, results of simulation with 250 MHz antenna gives better information due to better resolution that was not affected by the clutter disturbances introduced in the model. In the 500 MHz antenna results, the effects of the disturbances of clutters (gravel, limestone) are pronounced while in the result of 100 MHz, the model layers could not be distinguished.
Although the field and synthetic radagrams are not perfectly matched, the hyperbola marks at both approaches occurred at similar time range (22-25 ns). This implies that the reflection may be assumed to be the buried body.
The synthetic response of the different models is similar contrary to the assumption that they should be different. This may also be attributed to the model parameters gradation such as conductivities that were assumed. Moreover, the discrepancy in the models' responses may be due to the low values of reflection coefficient which is less than 0.05 as shown in Table 2; hence, the different zones in the models (3, 4, 5, and 6) could not be visible in the results of the simulation (synthetic radargrams).
The study has demonstrated that numerical modeling can be used to predict GPR signal response. Moreover, subtle information that may be difficult to extract during GPR field data interpretation may be discerned using method of modeling such as FDTD. The approach has also proven useful in the selection of appropriate GPR antenna central frequency for the real field measurement. Fig. 9 Correlation of field (a) and synthetic (b) radargrams for antennae frequency 250 MHz with red circle indicating the buried body