Stochastic Generation of Subsurface Profiles for Realistic Simulation of Through-the-Earth Communication Systems

Channel modeling for Through-the-Earth (TTE) communication has mainly been based on deterministic and empirical formulations. Deterministic models are limited to homogeneous rock structures while empirical models are only valid for the same configuration as the measurements that were taken. Given the difficulties in implementing underground measurements taking years to achieve enough database for a variety of mines and depths, we propose, as a more versatile solution, realistic simulations of a variety of underground structures to statistically characterize magnetic field attenuation. An algorithm based on the Monte Carlo method is presented and uses the Finite Element Method (FEM) to compute field propagation through each random structure. An empirical model is produced through simulations, and some metrics are evaluated, such as the median of magnetic field intensity, and the propagation loss, specially at the optimum operation frequency, in which channel variability is less spread when compared to all possible frequencies. This model can be used to estimate link budget for TTE communication in coal mines.


Introduction
Through-the-Earth (TTE) communication is a practical solution to link underground environments to the surface.Such systems are more resistant to accidents in the interior of mines that could affect primary wired networks [1,2].In underground mines, TTE communication is constrained by range, operating frequency, and soil apparent conductivity, which is an artificial conductivity calculated through simplistic models and magnetic field measurements.Statistical metrics, such as the median of propagation loss and optimal frequency, can be useful to adjust communication parameters between coaxially aligned antennas in coal mines.
Between 1970 and 1980, the American Bureau of Mines measured the intensity of magnetic fields for Throughthe-Earth communication for frequencies between 630 and 3030 Hz in 94 coal mines all over the United States [3].Using part of this data (27 mines), Durkin [4] estimated the apparent conductivity based on the homogeneous half-space model [5].The results led to apparent conductivity decreasing with increasing frequency and depth.Trying to analytically model this behavior, Hill and Wait then proposed the thin-sheet model to emulate the presence of conducting materials, such as metal sheets, cables, pipes, or higher conductivity layers near the surface [6][7][8].By setting appropriate parameters of the conducting sheet, the magnitude of the decrease in apparent conductivity could be predicted.Several years later, Yan et al. [9] analyzed data from all 94 mines to determine the conductivity behavior of the overburden based on the homogeneous half-space model, the thin-sheet model, and a combination of both models using the so-called Q-factor model to predict depth and frequency dependency of apparent conductivity.
All models mentioned above tried to investigate the channel behavior of a complex environment composed of multiple layers of different rocks using models based on 2 or 3 homogeneous layers.By varying the conductivity of the upper layer, these simple models try to fit the channel's behavior.Yan also proposed a more complex and laborious analytic model that characterizes the propagation medium as a stratified soil with different electrical conductivities for each horizontal layer in [10].Nevertheless, due to the complexity of such a model, it was applied in [10] for five layers only.All these models are site-dependent, but [4,9] also present some statistics on apparent conductivity; the authors, however, could not analyze the channel's variability at any optimal frequency due to the small number of frequencies considered in [3].
A feasible way to understand the propagation of electromagnetic waves in complex environments, where analytical formulas become very difficult, is through simulations of the propagation environment in order to generate a database for parameter estimation and modeling of the propagation channel.Often, the simulated propagation environment depends on statistical distributions computed from measurement data of the channel, as in [11].The derivation of these physical characteristics from measurements can lead to simulation results that extrapolate the original measurement results.In this sense, we propose realistic simulations of a variety of underground structures to statistically characterize magnetic field attenuation.These simulations are based on the Monte Carlo method and the Finite Element Method (FEM) [12].The environment parameters are set according to a proposed algorithm that stochastically generates subsurface profiles for TTE simulation using statistics extracted from the measurements in [3].No assumption of higher conductivity for upper layers is considered.Soil is set as stratified to a random number of horizontal layers, with uniformly distributed layer thickness, and each homogeneous layer with log-normally distributed conductivity.To obtain a statistical analysis of field behavior, 150 trials were conducted to emulate the empirical model.With this method, we were able to characterize the channel in a way that was not possible in [4,9] due to their limitation of working only with the four frequencies used in the measurements, and consequently, they were unable to perform analysis at optimal frequencies.
In the remainder of this article, we present the best-known models for TTE propagation in Section 2. Section 3 presents the coaxial model of the magnetic field transmission, along with the simulator's settings and the earth profile, followed by the proposed multilayer statistical model generated through electromagnetic field computation.The results are presented and discussed in Section 4, and the conclusions end this paper in Section 5.

Reference Models for TTE Communication
A uniform plane wave decays exponentially with distance traveled through the material medium, and such decay can be expressed as a function of skin depth δ.For a uniform conductive medium, the skin depth is the distance at which the intensity of a propagating plane wave decays to e −1 of its value and can be expressed as 2/μσω for low frequencies ω, where μ is the magnetic permeability and σ is the electrical conductivity.
The simplest approximation for a magnetic field created by an electrically small loop antenna is made by considering vacuum as a homogeneous medium, neglecting any boundary condition [13] and assuming a uniformly distributed current I on a loop of surface S and number of turns N, resulting in magnetic moment m d = NIS.The standard derivation of the magnetic field is realized using the method of retarded potential [14] and described in [15] in spherical coordinates.The adaptation of the vacuum infinite plane model in [15] to an Infinite Conductive Medium (IC) is realized by altering the wave number k = 1 − j /δ.This field approximation, in spherical coordinates, is described in its phasorial fashion by H = m d 4πr 3 e −jr/δ e −r/δ 2 cos θ 1 + 1 + j r δ r where r and θ are vectors in radial and elevation spherical coordinates and r is the distance between the center of the loop and the studied point.At a very close distance to the source (r ≪ λ/2π), the intensity of a time-varying magnetic field is similar to a static field calculated with the Biot-Savart law [16].
A more sophisticated model proposed by Wait and Spies in [5] considers the conductive medium as a homogeneous half-space (HHS), the transmitting antenna not being a point source.It considers a small wire loop antenna parallel to and above the ground at some height h 0 , and the H-field is measured at a distance h under the earth surface.The cylindrical coordinate system (ρ, ϕ, z) is used, and the earth is considered homogeneous.All displacement currents are neglected due to the assumption that all distances involved are smaller than the free-space wavelength.Hence, the magnetic field can be expressed in cylindrical coordinates in both radial ρ (horizontal field) and depth ẑ (vertical field) directions as follows:

2
where and A = a/h, D = ρ/h, Z = h 0 /h, δ = 2/ωμσ,a is the loop radius, and J 0 and J 1 are Bessel functions of the first kind.Figure 1 shows this scenario.All models cited above are in the frequency domain, since δ varies with ω.As mentioned at the Introduction, Hill and Wait [8] proposed the thin-sheet (TS) model to consider high conductivity in the upper layers due to the usual presence of cables, equipment, saline soil, and other conductive materials on the surface or in lower depth positions.Such modeling was motivated by the estimations of apparent conductivity that indicated higher values for shorter depths.

2
International Journal of Antennas and Propagation The models previously presented are valid approximations for a few homogeneous layers of soil and air.Still, they may be valid for estimating the propagation loss at a mine site with definite depth and soil conductivity at a specific frequency.In this study, a more simplified model of the multilayered channel is proposed.Using several random scenarios, based on characteristics of coal mines, it was possible to develop an empirical model for the calculation of the H-field inside such sites.This model provides a more practical way to estimate the median of magnetic field and propagation loss variability in the interior of an underground environment.
Models in [4,9] also present regression curves for calculating the channel behavior from a variety of measurements present in [3], but they preferred to estimate the apparent conductivity.These models are limited to frequencies of 630, 1050, 1950, and 3030 Hz, which were not necessarily the optimal frequencies for each measurement.

Channel and Earth Profile for Data Generation
In this section, we present the normalized magnetic field for coaxial configuration and the characterization of electric conductivity distribution based on the measured data for setting the simulator, followed by other parameters that model the earth profile.

Normalized H-Field and Propagation Loss. Considering only the vertical component of the H-field H z in transmis-
sions with antennas coaxially aligned, the resulting magnetic field normalized per unit of magnetic moment is in which I TX , S TX , and N TX are, respectively, the current, the surface area, and the number of turns N on the transmission loop.The channel transfer function may thus be calculated in [17] through The propagation loss in coaxial transmission can be found using L ω = 1/∥F z ω ∥ 2 , which is given in dB: The channel behavior tends to be similar to a bandpass filter, which consequently establishes an optimal transmission frequency where propagation loss is minimal.

Electric Conductivity Characterization through
Measurements.The data presented in [3] are the most comprehensive regarding magnetic fields in underground mines that can be found in literature so far, considering the complexity of measurement and the large number of sites.In this report, the authors investigated the detection of electromagnetic signals above coal mines to assist and rescue trapped miners.This study was requested by the American Bureau of Mines, and they classified 1222 active coal mines across the United States.
Their selection criteria were based on two considerations: the group of mines had to be sufficiently large to represent their physical uniqueness and operation characteristics, and the number of workers inside, due to the correlation between each mine's size and the quantity of miners.To do this, they used the technique known as optimum allocation, which states that for strata with great variability, it is necessary to perform a larger number of samples.
After a thorough screening, they selected 94 sites, whose depth ranged from 20 to 400 meters, located in ten American states: Pennsylvania, Ohio, West Virginia, Virginia, Tennessee, Kentucky, Illinois, Alabama, Utah, and Colorado.International Journal of Antennas and Propagation Due to the complexity and extensiveness of this study, there is no reason to believe that if similar measurements were carried out nowadays, their statistic behavior would be different.

H-field database
The apparent conductivity estimated through the HHS model applied to this measured data is the key to set the distribution of the electric conductivity in the Monte Carlo simulations.These measurements were carried out in coal mines operating at frequencies of 630, 1050, 1950, and 3030 Hz.
The raw data of field intensity is normalized by the magnetic moment for each measurement as in (4).Taking into account the complexity of the integral in (2), the apparent conductivity estimation cannot be obtained in a direct and manner.Field curves are computed from these equations for a variety of conductivity values.Such curves, or lookup tables, are used as reference maps for selecting the most adequate apparent conductivity and must consider the field intensity, depth, and operating frequency.Through a data fit comparison for several distributions using three tests, Root-Mean-Square Error (RMSE), Kolmogorov-Smirnov, and Cramer-von Mises, it was observed that log-normal was the most suitable distribution.From the distribution parameters for each frequency, a cubic spline interpolation and subsequent extrapolation with a 3rd degree polynomial was carried out to estimate the log-normal parameters μ lg = −3 8 and σ lg = 1 2.This was necessary to predict the distribution parameters which the proposed algorithm should apply to the Monte Carlo trials.

Multilayer Statistical Model through Electromagnetic
Field Simulations.There are very few reports or papers on electromagnetic simulations in TTE communications [6,10].Durkin in [6] recommended the Finite Element Method for TTE simulations; thus, it is employed in this study with the use of the commercial software CST Studio Suite®.To emulate soil variability, we modeled the ground as homogeneous layers with random physical characteristics.Still, we used random variables to describe the number of layers, their height, and their conductivity.For ensuring the variation of each scenario, we implemented an algorithm, shown in Figure 2, that uses electromagnetic simulation software for building the physical structures and calculating the H-field intensity throughout the downlink transmitting antenna axis.The Monte Carlo method is employed through a supporting VBA macro that generates all random values, sets the geometrical structures along with their soil electromagnetic characteristics, sets each trial boundary conditions, and restarts the simulation software.
With the parameters estimated in Section 3.2, we created a set of trials to study the H-field behavior in scenarios that mimic real-world conditions.The scenario used for this simulation is described as follows and shown in Figure 3.
The transmitting antenna is a one-turn square loop with sides of π meters.Such size is small enough to approximate a punctual antenna with constant electric current.It is located 50 cm above the ground and surrounded by air.Magnetic field probes are placed underground to measure the H-field.They are evenly distributed along the loop axis with a step of 10 m, in order to obtain a smooth variation of field.
Horizontal layers are modeled in such a way so as to represent the sedimentary sequences of rocks found in coal mines [7].They are modeled by square-base parallelepipeds with 50 m sides, whose large width was considered to avoid boundary problems for conductive media and determined after simulation tests.Their heights were uniformly distributed from 1 to 300 m, and the sum of all heights was limited to 300 m.Such a large span of heights tries to factor in contributions from more stratified profiles along with those of lesser stratification.The latter configuration was considered since the conductivity values used in the simulations were based on the HSS model.A random number of rock layers were selected from a uniform distribution, with a minimum of 3 and a maximum of 25 Regarding the electromagnetic characteristics of rocks, the magnetic permeability was set as a vacuum condition μ 0 = 4π10 −7 (H/m), which is common in the overburden of coal mines [15].Electric conductivity was programmed to be a random variable with log-normal distribution according to the extrapolated parameters given in the last section by μ lg = −3 8 and σ lg = 1 2. None of the structure parameters was tuned to achieve any expected result.
Using a FEM solver with tetrahedral meshing and absorptive boundary walls, 150 independent and random trials were carried out resulting in more than 200 thousand measurements of field intensity.Such quantity was selected to surpass the number of scenarios considered in [3] (94 mines) and used in [4,10] to achieve their models.In fact, not many trials are necessary to achieve acceptable coefficients for equation ( 4), shown later, that fits well to the median of measured data in [3].The model in (4) and its coefficients computed using 150 trials attained a RMSE of 3.7 dB, and when any combination of only 10 trials was randomly chosen, the RMSE was always below 5 dB.However, when fewer trials are used to generate the coefficients, the discreet CDF of propagation loss may be compromised.
Since the optimum frequency in TTE transmission is usually below 10 kHz for moderate and great depths (>100 meters), the frequencies were chosen in logarithmic steps from 0.1 to 10 kHz, and the H-field was normalized according to (4).For each trial, a different seed was generated, by CST®, to guarantee the randomness of each scenario.The postprocessing data analysis was realized by the use of MATLAB®.

H-Field Median Fitting.
A useful way to describe the average behavior of H-field intensity is through its median.An equation to model the median of H-field intensity as a function of frequency and depth by fitting the median curve of synthetic data is proposed as follows: where ω = 2πf , T = h ωμσ/2, μ = μ 0, and σ = e −3 8 = 0 022 is the equivalent value of median conductivity used in all simulations.To develop the proposed model in (7), we tested several equations based on regression analysis.Every protomodel took into account a structure with only two independent variables (depth and frequency) and was adjusted through a mathematical tool considering the RMS error and R-squared.By using such an equation for the regression, we introduced another condition to the model, which is the complexity of the mathematical expression.Thus, a trade-off was established between the complexity and accuracy of the model.With the inclusion of another independent parameter to the model, complexity increases but the accuracy only improves slightly.Conversely, when the number of parameters decreases, the model loses precision.For example, by excluding the term associated to a 4 , the RMS error increases to 8 dB and R -squared drops from 98% to 91%.The coefficients found are a 1 = −2 834, a 2 = 0 222, a 3 = −2 316, and a 4 = 1 115, whose fitting has R-squared of 0.998 and RMSE of 1.144 dB.

H-Field
Comparisons.In order to compare the model with measured data, Figure 4 shows the median of fields calculated from (7) (surface) and the normalized field from [3], only for the four frequencies employed in the measurements.It is possible to infer that simulated data resemble the measured data up to 250 m.The RMSE between (7) and the median of the normalized measured H-field is 3.7 dB and R-squared is 0.98, considering the four frequencies for all depths.It was not necessary to adapt or adjust the parameters of any specific layer to achieve such results.We repeated these calculations for the models presented above, obtaining a RMSE of 4.77 dB and a R-squared of 0.98 for the HHS model, and 4.63 dB and 0.98 for the model, respectively.The theoretical models were set for an apparent conductivity of σ = 0 022 S/m, which is the median value of the log-normal distribution extracted from the measured data available in [3].
Therefore, it is possible to suggest that the simulated model is valid, including the layer choice and conductivity distributions and their parameters.
To supplement the results presented here, a comparison was made between the deterministic models HHS and IC, the median of the simulated data, and the empirical model presented in this article.Figure 5 shows the density of the simulation trial data expressed as a function of number of trials (N) per dB per kHz, along with the theoretical models set for an apparent conductivity σ =  7 International Journal of Antennas and Propagation 0 022 S/m.Also, the estimated median of measured magnetic field for the frequencies of 630, 1050, 1950, and 3030 Hz is presented.Observing the field data, it can be noted that for small depths, there is almost no variation in the field as a function of frequency, because the observing point is in a quasistatic zone with respect to the transmitting antenna, while larger attenuation for higher frequencies is observed for medium and larger depths.Such field behavior in the quasistatic zone is confirmed by low-range TTE measurements performed by the authors and published in [17].Those measurements were performed in a different scale with smaller depths and higher apparent conductivities; however, they follow the same propagation principles.
It is noticeable in Figure 5 that deterministic models HHS and IC for σ = σ are more optimistic when compared to the median of simulated data and the empirical model.This that multilayer scenarios offer an extra loss for higher frequencies compared to models with one or even when the median conductivity of the multilayer is equal to the conductivity of the homogeneous overburden.Higher losses for higher frequencies for models with more layers are observed in [10] and confirmed here with a different multilayer model.
Observing the shades of trial density, a larger field variability for higher frequencies and depths is observed.This is because the variability of sedimentary layer dimensions is more sensitive to smaller wavelengths, and also the larger number of independent layers for greater depths contributes to the field variation.To better understand channel variability, a 90% prediction interval was computed for data around the median of the H-field, at the depth of 150 m for the frequencies of 5 kHz and 10 kHz.The results were 31 dB and 38.5 dB, respectively.However, this behavior may not be clearly observed in the results seen in Figure 6, where only fields at the optimal frequency are selected.
This figure shows the evaluation of channel variability observing exclusively simulated field values of minimum attenuation for each trial.The discrete cumulative probability of the propagation loss at the optimal frequency is computed for depths of 150, 220, and 300 meters.The L opt for 5%, 50%, and 95% are indicated along with the median values in round brackets, computed from the empirical model from ( 6) and (7).For the worst case (300 m), it was possible to observe a maximum additional loss of 16 dB for the 95% best scenarios compared to the scenario with 5% of cumulative probability.Therefore, a 16 dB change can cover most of the channel variation.However, this variation is as high as 20 dB for a depth of 150 m.An even greater variation would be observed if the CDF of a 60 m scenario was measured, but the optimal frequencies would be much greater (>100 kHz) than those of interest to this paper.Therefore, a greater variability of H-field that could be 8 International Journal of Antennas and Propagation visually observed for frequencies above 5 kHz and depths above 150 m in Figure 6 may be of less importance; from a practical standpoint, transmissions are set to operate at the optimal frequency, which tends to be lower the higher the depth, producing less field variation.

Conclusions
The proposed stochastic generation of a variety of underground structures for electromagnetic simulations provides more adequate conditions for TTE channel modeling than simplified deterministic models and the empirical models in [4] and [9], especially in terms of frequency choice.
Our empirical model of the H-field median may provide improved orientation for link budget calculation in TTE communication in coal mines.When compared to deterministic models, the empirical model in ( 7) is better fitted to the median of simulated and measured data.Compared to models based on statistics of apparent conductivity, our model is valid for a wider range of frequencies and not only for the four frequencies in [3].Besides, it considers the channel variability for the optimal frequencies, which was not possible in [4,9].Regarding the variation in magnetic field, for smaller depths, the variance among trials is greater despite the lower overall attenuation.For deeper positions, the attenuation becomes more severe at higher frequencies as the field moves from the transition zone to the far field.
With regard to the CDF of propagation loss, it is observed that for a 300 m depth, which is the worst-case condition, there is an additional loss of 16 dB for the scenario with a cumulative probability of 95% compared to the scenario with 5% of cumulative probability.The variance, which may be seen from data collected at the optimal frequency for each trial, is smaller than that seen for all possible field values at frequencies higher than the optimum.Such smaller variation augments the predictability of the link budget in a TTE communication even though, in praxis, it also depends on atmospheric and manmade noise conditions at each site.The study of electromagnetic noise levels at the surface and underground as in [7,8,18,19] may be used in conjunction with this paper to compute the median and the limits of the link budget used for communications in coal mines through vertical-axis transmission.
The model, presented in this paper, used 150 independent data points to achieve a RMSE of 3.7 dB (while HHS and IC models presented 4.77 dB and 4.63 dB, respectively), and the RMSE was always below 5 dB for any randomly chosen combination of 10 trials.Although the estimation of the median of field intensity would converge with much fewer than 150 trials, the variability analysis from the discrete CDF would be impaired with fewer propagation loss samples.The modeling process used here can be applied not only to coal mines but also to any sedimentary structure.In such cases, the statistical characteristics of soil extracted from magnetic field measurements at different sites may be used to tune the simulations to new conditions.

Figure 1 :
Figure1: Geometry used for the calculation of the magnetic field at the point P h, ρ for a circular loop antenna at h 0 above ground.

Figure 2 :
Figure 2: Simulation process algorithm.The red boxes represent MATLAB processes, the orange boxes represent VBA processes, and the green boxes represent CST processes.

Figure 4 :
Figure 4: H-field median from the empirical model in (7) (surface) compared to the measured data (circles).

Figure 5 :
Figure 5: Comparison of different models and measures of normalized vertical magnetic field for depths of 60, 150, and 210 meters.

Figure 6 :
Figure 6: CDF of propagation loss at optimum frequencies for depths of 150, 220, and 300 meters.