GMODL: An Indonesian MATLAB-based ground-penetrating radar data modeling and processing software

Ground-penetrating radar (GPR) data modeling and processing is crucial to near-surface geophysics. Its use in geohazard mitigation, construction, shallow hydrocarbon contamination, and other shallow subsurface detection is undeniable due to its high-resolution imaging. GPR software (MATGPR) that we used currently requires access to MATLAB which not everybody can use because of its licensing prices. Thus, we developed this program (GMODL) that uses finite-difference time-domain (FDTD) and split-step Fourier algorithms. To test the software, we created synthetic models. The synthetic model used for the testing is a river model for flood mitigation that consists of a layer of freshwater with ρ = 20 Ωm, k = 81 and μ r = 1 of depth 5 m, two layers of sandstone with ρ = 850 Ωm, k = 2.5 and μ r = 1 of total depth 4 m, and a layer of claystone with ρ = 120 Ωm, k = 11 and μ r = 1 of depth 1 m. The GPR antenna frequency is set to 250 MHz. The testing algorithm of GMODL shows that it can create a synthetic radargram for the river model. The boundary layers are obvious to identify. Also, the freshwater thickness can be determined from the radargram.


Introduction
GPR is a non-intrusive, near-surface geophysics method that uses electromagnetic (EM) waves to image the subsurface. GPR used two antennas, a transmitter that transmits EM waves and is reflected by a reflective surface into the receiver. The two-way travel time (TWTT) is then plotted against the scan-axis, creating a radargram. The radargram's TWTT is then converted to depth estimates using 1D EM wave propagation velocity structure [1].
To produce a radargram from synthetic modeling or field data acquisition, there needs to be software that encompasses the algorithms necessary for GPR modeling. GMODL provides the tools necessary to model radargrams such as FDTD and split-step. The program is executable and doesn't need to be accessed from MATLAB. For testing purposes, we have decided to make a synthetic radargram of the river model to determine freshwater thickness because river monitoring using GPR is used for flood mitigation.

Finite-Difference Time-Domain (FDTD)
The FDTD modeling in GMODL uses a concept from Irving and Knight 4 . It is used by applying the TM (Transverse Magnetic) concept to an antenna perpendicular to the x-z plane. To prevent reflection from the edges of the model, we used the coordinates for the absorption boundary at a perfectly matched layer (PML) [2].
For TM cases on the x-z plane, the magnetic waves are in the x-direction ( ) and the z-direction ( ). Then the electric waves in the y-direction ( ) are as follows: σ, σ, ϵ, and μ are each physical parameters in the form of electrical conductivity, dielectric constant or permittivity, and magnetic permeability. k is stretching the coordinates in the x,y, and zdirection. The equation above is modeled with finite-difference using leapfrog. With this approach, the electric field and the magnetic field become balanced as in the partial derivative approach on the same grid using a fourth-order finite-difference for space and second-order for time.
FDTD uses an algorithm proposed by Kane Yee in 1966. This algorithm uses second-order on central difference [3]. The algorithm is as follows: 1. Ampere's Law and Faraday's Law in Maxwell's equation are modified, the differential form is converted into the finite difference form. Discretization of space and time so that the electric field and the magnetic field are both in space and time. 2. Then the problem of the modified equation is solved to obtain the "update equation" as the shape of the plane of the object that is not known from the shape of the plane of the object that is known. 3. Evaluate the magnetic field at one time on an unknown object so that we find that unknown object (effectively this is obtained from a known object). 4. Do the same step as step 3 but this time we evaluate the electric field. 5. Then, repeat the previous two steps until the field obtains the desired results.
Within GMODL, there exists 2 different data set to store data. The first one is called Old Data (OD). OD dataset is used for storing imported data and processed data that will be processed further. The second one is called New Data (ND). This ND data set is used to store data resulted from the process within GMODL. In GMODL field data processing, imported data from field acquisition gets saved to the OD. When field data is processed, the results will be stored temporarily in the ND dataset. The user then can review the data within ND and choose whether they want to keep the data within ND or to discard it. If the user keeps the data in ND, then the data within OD will be replaced by the new data from ND. Then the ND data set will be emptied and ready to store another data resulted from the process within GMODL. The same process happened for modeling. The results from modeling will be stored in the ND data set and can be discarded or be stored in OD for further processing.

Split-Step Fourier Method
In GMODL, there are two main concepts used for Split Step Fourier Method, constant reference slowness and perturbation term [4]. Fourier transformation methods are used to solve the acoustic wave equation with constant lateral velocity. Within the split-step Fourier method, the reciprocal of IOP Publishing doi:10.1088/1755-1315/1031/1/012026 3 interval velocity is divided into two parts: a constant reference slowness and a perturbation term. Reference slowness is used to migrate across intervals for the wave-frequency domain as a form of the constant-phase value for phase-displacement.
Consider a wave equation with the propagation of compressional waves within an acoustic constant density medium: With p = p(x,y,z,t) is the pressure dan u = u(x,y,z) is the medium slowness, which is defined as the inverse of half the propagation velocity u(x,y,z) = 2/ν(x,y,z) with exploding reflector model. After equation (4) has been transformed, the domain frequency can be obtained as: with, Where r is the horizontal position vector which is defined as the following: With slowness field value u(r,z) into 2 components: ( , ) = 0 ( ) + Δ ( , ) By defining to(z) as reference slowness value, equation (7) is then substituted to equation (5): With, ( , , ) = 2 (2 0 Δ ( , ) + Δ 2 ( , )) ( , , ) We then obtain, Therefore, homogeneous acoustic wave equation (5) is transformed into the inhomogeneous wave equation (8) because there is a source S(r,z,ω) during the varied slowness process. The solution of equation (8) is then used for the split-step Fourier method by ignoring Δ 2 as a contributing factor. We can conclude that in the previous Fourier transformation, the wave is reflected as a field wave at depth , − ( , , ), from r to k. Then, Then, we use the second phase shift for perturbation inside slowness, Δu (r,z) = u(r,z) -u0(z), within interval Δz : − ( , +1 , ) is then integrated within the frequency, from 1 to 2 to obtain the displacement data for every depth +1 : We do the first identical phase shift for the constant-velocity phase shift migration value. The secondary phase shift functions as a correction factor provided that it acts as a time shift based on the difference between the actual value and the reference slowness of each spatial position.

Results
GMODL can make synthetic models for forward modeling and process field data. The synthetic model is set with an antenna frequency of 250 MHz. First, a homogeneous model is made to check the fitness of the grid and radargram size. The model is based on the standard river stratigraphy. The model's parameter values are reference values [5,6] . The model consists of a layer of freshwater with ρ = 20 Ωm, k = 81 and μr = 1 of depth 5 m, two layers of sandstone with ρ = 850 Ωm, k = 2.5 and μr = 1 of total depth 4 m, and a layer of claystone with ρ = 120 Ωm, k = 11 and μr = 1 of depth 1 m. The model's radargram is made using FDTD and split-step. The models and their radargrams can be seen in figures 2 and 3. The split-step Fourier method's radargram is closer to an actual radargram than its FDTD counterpart. The boundary layers of the split-step radargram are distinguishable ( figure 1b and 2b). And the layering of the homogenous and river model (figure 1b and 2b) is perfectly aligned with the geological model ( figure 1a and 2a) Based on the TWTT (top TWTT = 0 ns, bottom TWTT = 300 ns) and velocity structure [1] of the freshwater layer (calculated to be 0.033302 m/ns) and using the layer thickness equation [7] the freshwater thickness of the split-step layer is calculated to be 4.9953 m. The FDTD radargram's TWTT range of 0 -15 ns (figure 1c and 2c) produces a layer of fresh water with a calculated thickness of 0.249765 m. The distinct difference of the FDTD radargram's TWTT is likely caused by the velocity pull-up caused by the contrast claystone layer. The split-step radargram is closer to the geological model than the FDTD radargram. Despite the differences, the FDTD radargram of the homogeneous and river model ( figure 1c and 2c) shows correctly the layering of the geological model ( figure 1a and 2a)

Conclusion
The program that used in GMODL is executable and doesn't need to be accessed from MATLAB. GMODL provides the tools necessary to model radargrams such as FDTD and split-step. The actual radargram result is closer to an actual radargram result when using split-step Fourier method's than using FDTD. For the boundary layers of the split-step radargram are distinguishable and obvious to