In situ time-zero correction for a ground penetrating radar monitoring system with 3000 antennas

The time-zero correction is an essential step in the data pre-processing of ground penetrating radar (GPR) measurements to obtain an accurate signal propagation time between transmitting and receiving antennas. For a novel custom GPR monitoring system with about 3000 antennas and corresponding transceiver structures placed around a soil sample (lysimeter), an in situ approach for the time-zero correction is required. In particular, unknown material properties between any pair of transmitting and receiving antennas prevent a conventional time-zero correction. We present and compare two calibration approaches, namely a pairwise and a mesh calibration, both utilizing the ability of the monitoring system to conduct reciprocal measurements between any pair of antennas. The pairwise calibration enables an individual calibration for any antenna pair, whereas the mesh calibration reduces the influence of the soil between antenna pairs compared to the pairwise calibration. The developed approach is verified by utilizing a mathematical model. Experimental results from a simplified setup show that the lysimeter filling has a negligible impact onto the calibration approach based on adjacent measurements for the mesh calibration. In addition, it is shown that a state of the art time-zero calibration can be used to measure the signal delays within the analog circuit of the measurement system with an accuracy of ±4 ps. The simulation results indicate that by using the developed concept, no prior air calibration between every possible antenna combination is necessary. Thus, this work provides a crucial contribution towards an automated in situ time-zero correction for 3D GPR monitoring systems with many antennas.


Introduction
Ground penetrating radar (GPR) is a non-invasive electromagnetic (EM) geophysical technique to investigate the shallow subsurface of soils by emitting radio frequency EM waves into the soil and receiving a response signal via antennas (Jol 2008). GPR yields higher spatial resolutions than electromagnetic induction (e.g. Doolittle and Brevik 2014) and has a faster data acquisition than electrical resistivity surveys (e.g. Samouëlian et al 2005). The propagation of emitted EM waves depends on the relative dielectric permittivity ε r and the electrical conductivity σ of the subsurface. Typical permittivity values vary in the range of ε r ≈ 1 for air and ε r ≈ 80 for water, making the permittivity highly sensitive to soil water content (SWC) (Topp et al 1980, Dobson et al 1985, Roth et al 1990. This allows to observe SWC variations using appropriate petrophysical models (Huisman et al 2003, Paz et al 2017, Klotzsche et al 2018, Wu et al 2019. Measuring the SWC and, thus, flow and transport processes with a high accuracy assists agriculture (Bolten et al 2010, De Benedetto et al 2019, Yang et al 2021, land rearrangement (Xu et al 2014), and leak detection (Hyun et al 2007). In general, a better understanding of soil parameters helps to develop more sophisticated soil-plant-atmosphere models (Vereecken et al 2016). Besides applications related to the SWC, GPR is used for the noninvasive inspection of roads (Saarenketo and Scullion 2000, Hammond et al 2011, Ihamouten et al 2018, concrete bridges (Hugenschmidt andMastrangelo 2006, Dinh et al 2019) and railways (Ciampoli et al 2019), the assessment of contaminated areas (Catapano et al 2014, Wijewardana et al 2017, and the detection of landmines (Sai andLigthart 2004, Lameri et al 2017), among others.
One particular GPR measurement configuration is the crosshole setup, where transmitting and receiving antennas are placed inside boreholes (Slob et al 2010, Qin et al 2021. The crosshole setup in combination with acquiring multi-offset gathers enables the GPR operator to determine a tomographic image of ε r and σ and, therefore, the SWC between the boreholes with a higher spatial resolution in deeper soil areas compared to on-ground or off-ground surveys (Binley et al 2001, Liu et al 2017. Especially in the last decade, crosshole GPR full-waveform inversion (FWI, Ernst et al 2007, Meles et al 2010, Klotzsche et al 2019 has shown a high potential to derive higher decimeter scale images of the subsurface compared to standard approaches (Gueting et al 2017, Zhou et al 2020. In general, GPR systems derive the soil's permittivity by determining the propagation time and, hence, the velocity of EM waves in the soil. In order to determine the signal propagation time between transmitting (Tx) and receiving (Rx) antennas, a time reference for the measured signals has to be defined, termed as time-zero, t 0 (Cassidy 2009). The time-zero correction aims to find a reliable time reference and ensures that no additional signal propagation time from within the measuring instruments, e.g. system timing and cable length, is included in the data (Olhoeft 2000). Beyond that, the time-zero correction monitors changes in time-zero due to, e.g. thermal drifts, damaged cables or electronic instability (Huisman et al 2003) and is, therefore, time dependent. Furthermore, Klotzsche et al (2019) pointed out that the correct time-zero correction is a crucial step in the pre-processing to accurately apply an FWI to the data.
Time-zero is commonly defined as where t a is the arrival time of the direct air wave between Tx and Rx, d is the propagation distance of the direct air wave and c 0 is the speed of light (Annan 2003). Note that the air wave always arrives first at Rx under a line of sight condition, where the shortest path between Tx and Rx runs through air. Due to antenna near-field effects and an interference of the air wave and ground wave, the direct air wave can have a velocity that differs from the assumed c 0 (Gerhards et al 2008), resulting in an erroneous time-zero when using (1). Thus, it is advantageous to use a linear regression to deduce time-zero and the signal's velocity v from multiple measurements with different distances d between Tx and Rx, where The calibration data required to use (2) is acquired either by changing the distance between Tx and Rx for multiple measurements or by sequentially lifting Tx and Rx above a known reflector and measuring the two-way propagation time (Yelf 2004). Instead of measuring time-zero directly via (1) or (2), another approach is to determine time-zero via an inversion. Gerhards et al (2008) implemented a multipoint inversion that minimizes an objective function, which contains the air wave propagation time as a free parameter to calibrate a multichannel GPR system. Kaufmann et al (2020) recently presented a novel time-zero calibration method to calibrate the wide angle reflection and refraction (WARR) machine (Diamanti et al 2018), a GPR on-ground system that utilizes one transmitting antenna and seven receiving antennas. Each receiving antenna has its own data acquisition hardware, making a time-zero correction for each receiving antenna necessary. The final calibration scheme yields precise time-zero by finding the minimum of a defined objective function, while the calibration measurement itself is time consuming with measurement times of about 40 min. In this study, it was observed, while the setup was always exactly the same, that the time-zero difference of each Tx-Rx pair to the first pair was stable and only an overall time-zero shift needs to be determined. Radzevicius (2015) used a combined least-squares and grid search algorithm to determine time-zero along with further model parameters. Viriyametanont et al (2008) conducted an effective time-zero correction based on the direct wave propagating inside concrete instead of using the direct air wave. This method requires known reflector depths inside the concrete. For crosshole measurements, a direct air wave between Tx and Rx placed in separate boreholes is not available. Therefore, the antennas are repeatedly placed in air after several crosshole acquisitions, e.g. after ten measurements (Oberröhrmann et al 2012), to conduct time-zero correction measurements according to (2) and to capture the time variability of time-zero. Afterwards, time-zero values from successive time-zero correction measurements are interpolated to obtain an individual time-zero for each crosshole measurement. This approach is very time consuming as calibration measurements have to be performed repeatedly. Oberröhrmann et al (2012) introduced a cross-correlation method to enhance the time-zero accuracy for crosshole measurements by linking traces that are captured at comparable Tx-Rx positions. This is especially important due to the sensitivity of the FWI to small time-zero errors as emphasized by Klotzsche et al (2019).
Until now the time-domain crosshole GPR FWI was mainly applied to investigate aquifers and related processes within decimeter scale. Within the last years, a need to detect smaller structures and to map more detailed and faster processes like preferential flow paths or root zones becomes more and more important to enhance model predictions. To achieve this, higher frequency antennas that provide a higher resolution and faster measurement techniques are essential. Currently, we are investigating a GPR monitoring system that extends the crosshole setup. The monitoring system will allow to radiate EM waves into a soil column from all spatial directions, similar to, e.g. Stoffregen et al (2002) and Schmalholz et al (2004). This is achieved by permanently mounting up to 3000 antennas on the outside of a lysimeter (a cylindrical shell, used by, e.g. Hannes et al (2015) and Pütz et al (2016)). The aforementioned FWI is then used to analyze the measured GPR data. The interested reader is referred to, e.g. Klotzsche et al (2019) or Mozaffari et al (2020) for exemplary images of 2D and 3D inversion results and their applications. So far, we are not aware of any system that uses such a large number of antennas for tomographic images of the soil. Such a permanently installed monitoring setup prevents the use of existing calibration approaches that utilize the air wave for an accurate timezero correction since there is always unknown soil between the antennas. Additionally, the large amount of possible Tx/Rx combinations requires a novel automated in situ approach for the time-zero correction, because a manual time-zero correction for every possible Tx/Rx combination is not feasible. Finally, the method must be able to perform the time-zero correction periodically to capture possible changes. Applications that also use a large number of antennas are found, for example, in beamforming radar. Here, the calibration of the antennas is often performed via measurements within a known medium (Hu et al 2013) or with the help of target objects (Harter et al 2016). Such approaches are not applicable for our system as the signal propagation time is initially unknown due to the unknown medium between antennas.
In this paper, we address five related topics. After (a) introducing the setup of the monitoring system, (b) two novel timezero correction methods are presented, and (c) the accuracy of calibration measurements is assessed. Next, a (d) simulation is utilized to analyze the feasibility and in particular the accuracy of the methods. A mathematical model is used within the simulation to verify the correction approach. As the calibration method is an essential part of the monitoring system and as the system is currently in the design phase, it is important to develop and evaluate the calibration approach at an early stage. Thereby, the system can be designed in such a way that the presented approach can be implemented. Finally, (e) preliminary experimental measurements are evaluated to assess the influence of the lysimeter filling onto coupling signals as these can be used for the correction. This reduced measurement setup nevertheless contains all relevant aspects of the later system setup (size and spacing of the antennas and thickness of the lysimeter wall) and can therefore be used for verification. Our approach is able to time-zero calibrate the multichannel monitoring system in situ without the need of direct air signals between transceiver pairs. Our goal is to determine the final time-zero value with an accuracy of 25 ps. This work is intended as a feasibility study of the calibration concept in the design phase of the monitoring system.

Monitoring system
The planned monitoring system consists of around 3000 transceivers that are placed on the outside of a lysimeter, which is a cylindrical column filled with soil (see figure 1). Each transceiver is in principle composed of one antenna, one Tx circuit and one Rx circuit. The 3000 transceivers are, e.g. equally distributed into 120 transceivers per row around the lysimeter and 25 transceivers per column. The lysimeter is made of polyvinyl chloride (PVC), has a height of 1.5 m, an inner cross-sectional area of 1 m 2 , and a wall thickness of 30 mm. A master module is used to trigger the transceivers, where t i defines the time that the trigger needs to arrive at the corresponding transceiver i. For 3000 transceivers, t = (t 1 , . . ., t 3000 ) T marks the trigger times of all transceivers.
For an entire tomographic measurement, each transceiver is used once to transmit, while the other 2999 transceivers serve as receivers. Unknown soil properties and, therefore, signal velocities between any pair of antennas initially prevent the use of direct signals as in (1) for a time-zero correction. Note that in most cases the lysimeter is filled with soil by pushing it into the ground such that the antennas can only be installed after the filling.
Next, important time segments for a measurement between two transceivers i and j are discussed to derive a method for the time-zero correction (see figure 2). The time axis starts with the emission of the trigger at the master module. The signal generator of the transmitting transceiver i starts at time t i , while the data acquisition of the receiving transceiver j starts at time t j . The signal is emitted by the antenna of i into the lysimeter at t i + t Tx . In figure 2 the times are rearranged to illustrate the context between time-zero and the trigger times of the corresponding transceivers. Using the definition of t 0 = t a − t p with t p as the propagation time of the signal between the antennas, it is readily observed that where we define  t Rx are the analog circuit propagation times and describe the times that a pulse needs to propagate within the Tx and Rx circuit, respectively. t a i j describes the arrival time of the signal when i transmits and j receives. t p i j is the propagation time of the signal from the transmitting antenna of transceiver i to the receiving antenna of transceiver j. Note that the propagation time t p i j is initially unknown and cannot be determined without a proper calibration. The gap between t j and t Rx exemplary illustrates that the data acquisition of transceiver j may start before the signal is emitted by transceiver i. Note that changing the order of t Rx and t p i j does not change the arrival time. Finally, t 0 i j is time-zero when i transmits and j receives.
as the trigger offset of transceivers i and j. Therefore, the determination of time-zero for any transceiver pair requires the knowledge of the difference between the trigger times of the corresponding transceivers. The trigger times primarily differ due to different cable lengths to the transceivers. If t Tx and t Rx are equal, respectively, and stable for all 3000 transceivers, these times just have to be calibrated once. This is reasonable because the system design aims to have equally built transceivers with negligible differences in t Tx and t Rx compared to the trigger times. However, this calibration approach is also feasible for systems with different analog signal propagation times, i.e. t Tx i ̸ = t Tx j or t Rx i ̸ = t Rx j , to the antennas as long as these delays are known and stable over time. This allows to connect multiple antennas to the same data acquisition channel to keep hardware costs low.
In the following, for simplicity, we discuss the system with equal analog transceiver propagation times t Tx and t Rx and one Tx/Rx unit per antenna.

Pairwise calibration
The pairwise calibration aims at determining the propagation time of a signal between any two antennas. The approach is based on the ability to use any transceiver either as Tx or Rx. This allows to conduct reciprocal measurements between any two transceivers i and j, yielding The term reciprocal measurement always refers to a pair of two normal measurements as in (5) and (6). The advantage of reciprocal measurements between two antennas is that the propagation time t p of the signals is exactly the same in both directions. As a consequence, the two propagation times t p i j and t p j i to be determined originally can now be simplified into one propagation time t p , i.e. t p := t p i j = t p j i . This approach is a significant improvement as it reduces the number of unknowns in relation to the number of independent measurements, enabling to uniquely determine the remaining unknown t p . Reciprocity applies to all passive linear systems, such as the soils under study. Adding (5) and (6) and rearranging results in This demonstrates that the unknown propagation time of the signal between i and j through the unknown medium, which usually can only be determined via a calibrated time-zero, is directly determined by conducting a reciprocal measurement as long as the signal delay within the analog circuit, t Tx + t Rx , is known. The advantage of this method is that each individual transceiver pair can be calibrated regardless of their position, and that t Tx + t Rx only has to be calibrated once, assuming that it is a constant value for all transceiver pairs. In order to visualize fast soil changes, fastest possible tomographic measurements are required. In these cases, the pairwise calibration has the disadvantage that each signal trace needs to be included twice in the measurement data to obtain reciprocal measurements. If static soil conditions throughout the measurement are assumed, the reciprocal trace is redundant and increases the measurement time for a complete tomogram. In addition, dynamic irrigation processes might disturb the time invariance of the soil (t p i j ̸ = t p j i ) and would, therefore, lead to erroneous calculations when using (7).
The pairwise calibration is well suited for passive linear time invariant systems that show a reciprocal behavior. In these scenarios, the pairwise calibration should be the method of choice. However, if fast measurements are required, the calibration needs to be extended or replaced by the following mesh calibration approach.

Mesh calibration
The mesh calibration aims to determine t p via calibrated timezero values and is based on a superposition of reciprocal measurements between adjacent transceivers. Thereby, the influence of the nearby soil on the signal is minimized in our lysimeter setup. The reciprocal measurements ensure that t Tx , t Rx and t p are eliminated from the following equations and, as described in section 2.2, reduce the number of unknowns. The trigger offset t i j between any two transceivers is calculated via a superposition of the coupling signals between adjacent transceivers. To illustrate this, we start by subtracting (6) from (5), which yields so that we obtain the sought value in (3) by subtracting the known arrival times. Note that only the reciprocity of the medium between the antennas allows to eliminate the unknown signal propagation times between any two antennas. Equation (8) also holds for neighboring transceivers. If we want to determine, e.g. t 1 3 , then a first reciprocal measurement is conducted between transceivers 1 and 2, whereby 2 is a direct neighbor of 1. Using (8), it follows that Continuing the same way with a reciprocal measurement between transceivers 2 and 3 yields Calculating (9) + (10) results in This procedure is repeated until any delay time is calculated. Care should be taken that each additional transceiver that is used to determine the trigger offset for any two transceivers might lead to an error amplification due to inaccurate measurements. A more sophisticated approach is to define a matrix A such that where t are the trigger times of the 3000 transceivers, and t obs is a vector containing the observed trigger offsets of neighboring transceivers as in (9) and (10). For where B defines the sought transceiver combinations and t searched is the vector containing the searched trigger offsets (e.g. t 1 3 ), it follows that A + denotes the pseudoinverse of A, which can be calculated by a QR factorization or singular value decomposition. Using the matrices A and B, therefore, allows to determine any trigger offsets simply based on coupling measurements. This approach is similar to the approach of Xu and Noel (1993), who used the superposition of a primary data set to synthesize resistivity data for arbitrary measurement configurations. We solve (14) with a classical optimization method. The innovative part here, however, is that only the reciprocal measurements allow us to set up a solvable system of equations. Finally, the propagation time of a signal between any transceivers i and j is determined as Note that the usage of coupling signals for the calibration does not require additional measurements to be performed during a tomography, since every transceiver is once used as Tx for a complete tomography. The only additional effort is to use neighboring transceivers of a transmitter as receivers. The mesh calibration can be performed before irrigation processes to obtain a time-zero correction based on (3) for any two transceivers. When the soil then experiences fast changes, it is sufficient to use 3000/2 transceivers on one half of the lysimeter as Tx for the soil measurement and, therefore, the acquisition speed is increased. In contrast to the pairwise calibration method, calibration data for the mesh calibration, i.e. trigger offsets, only need to be determined once or after a specific period of time to re-calibrate the system. Lastly, this approach does not rely on coupling measurements between neighboring transceivers. The system of linear equations can be constructed using any combination of reciprocal measurements, as long as all transceivers are connected via a mesh.

Random timing errors
The treatment of timing errors is important for both, the pairwise and mesh calibration. We distinguish between two kinds of timing error sources, namely systematic and random differences between t Tx and t Rx for different transceivers. While the impact of random errors, namely jitter, is reduced by stacking the signals, systematic differences are unaffected by stacking. The ability to measure the trigger offsets of transceivers based on reciprocal measurements as in (8) relies on equal analog circuit propagation times t Tx and t Rx . Any inequalities of analog circuit propagation times would result in systematic errors for both, the pairwise and mesh calibration, which are the subject of a future investigation. In this work, the effects of random timing errors on the trigger offsets are analyzed by means of simulations. Three possible measurement topologies for reciprocal measurements are depicted in figure 3.
While the usage of Mesh2 is sufficient to solve (14), Mesh4 and Mesh8 are introduced to reduce the influence of jitter by utilizing a larger system of linear equations. The 3000 transceivers are separated into 120 columns around the lysimeter and 25 rows on top of each other.
Equations (16) and (17) are used to implement random error effects: and where ξ i and ξ j ∼ N(0, η 2 ) represent errors that are introduced by transceivers i and j either as Tx or as Rx. Subtracting the erroneous measurements (16)-(17) and rearranging yields where ξ i j is a new random variable with ξ i j ∼ N(0, η 2 ) following the properties of normal distributed random variables. Therefore, the true trigger offsets t obs in (14) are replaced by measurest Elements in the error vector ξ are not necessarily independent, because when one transceiver transmits and multiple transceivers receive, then an error in t Tx affects all receivers. We use (19) to implement the mathematical model of our simulation. The same calculation can be done for the summation of erroneous reciprocal measurements as it is required for the pairwise calibration, which again leads to an error term ξ i j ∼ N(0, η 2 ): It follows that the pairwise calibration exhibits the standard deviation of a single measurement (η), which can be reduced by repetitive measurements (stacking). A stacking factor s would finally reduce the random errors standard deviation by a factor of √ s which is also true for the mesh calibration. For the simulation, the matrix B in (14) is defined such that all possible 3000 · 2999 trigger offsets are calculated. Finally, we assess the results by calculating the normalized root mean square (nRMS) of the calculated erroneous trigger offsetst searched between any two transceivers i and j by utilizing the true trigger offsets t searched as where η is the predefined standard deviation of random errors and N is the number of simulations conducted with the same value for η. This allows us to investigate how the single standard deviation η transfers to the searched trigger offsets. Besides stacking the data by repeating the measurements to obtaint obs , random errors are furthermore reduced by performing reciprocal measurements between as many neighboring transceivers as possible.

Detection of systematic timing errors
Equal analog propagation times t Tx and t Rx , respectively, are a prerequisite for the calibration presented so far (or at least have to be known and stable). The aim of this section is to establish an in situ measurement approach to test this prerequisite. In order to detect systematic timing errors, (7) is used. This is possible under the condition that t p is equal for all coupling signals. We use the term coupling signals to refer to signals that propagate between adjacent transceivers, more specifically between their respective antennas. If it holds that t p is equal for all coupling signals, then t a i j + t a j i in (7), where i and j are now adjacent transceivers, is equal for all reciprocal measurements as long as t Tx + t Rx is equal for all transceivers. This means that t a i j + t a j i is a measure to control whether t Tx + t Rx is equal for all transceivers. If, e.g. coupling measurements between transceivers 1-2 and 2-3 are conducted, then as long as t p is the same for both measurements. A deviation of t a 1 2 + t a 2 1 and t a 2 3 + t a 3 2 indicates a violation of the calibration requirement that t Tx + t Rx is equal for the transceivers. This calculation is only valid if the coupling signal's propagation time t p does not depend on the lysimeter filling. This will be further analyzed in subsection 3.2. Note that we have to distinguish between three propagation times for coupling signals, namely for horizontal, vertical and diagonal coupling signals. Similar approaches are used in mutual coupling calibrations for phased-array antennas, where phase differences between active antenna elements are of interest (Aumann et al 1989, Shipley and Woods 2000, Bekers et al 2013.

Experimental setup
To measure the influence of the lysimeter filling on t p between adjacent antennas (horizontal, vertical, diagonal) as required for the detection of systematic errors, a test setup as depicted in figure 4 is used. In order to have a realistic estimation of coupling signals, an antenna array as planned for the GPR system is used ( figure 4(a)). This array consists of 16 circular bowtie antennas. To refer more easily to specific antennas, a number is assigned to each antenna. The size of a single bowtie antenna of 3 cm × 6 cm (Mester et al 2020) is derived from the requirement to place 3000 antennas around the lysimeter. The wall of the lysimeter is represented by a 1 m × 1 m PVC plate with a thickness of 30 mm. The antenna array is taped to the PVC plate. An additional box made of PVC can be filled with water and placed behind the PVC plate to simulate changing material properties inside the lysimeter and, thereby, test the influence of the lysimeter filling onto coupling signals.
The experimental setup represents a reduced setup, but it contains all relevant elements of the later system, i.e. antenna size, antenna spacing and thickness of PVC wall, and thus allows a realistic estimation for the true system behaviour. A Ricker wavelet (second derivative of a Gaussian function) with an amplitude of 1 V, a center frequency of f c = 750 MHz and a bandwidth of 620 MHz is used as the transmit signal and is generated by a Keysight arbitrary waveform generator (AWG) M8190A. The signals are sampled with 4 GSa/s and averaged by a factor of 64 using a Keysight DSAX91304A Infiniium high-performance oscilloscope. One coaxial cable is used to connect the AWG to the Tx and one to connect the Rx to the oscilloscope. Finally, the signals are interpolated using a 250 point sinc interpolation to get higher timing resolutions when picking time lags between signals via a cross-correlation. Figure 5 exemplifies the signal coupling between adjacent antennas. All three coupling directions (horizontal, vertical, diagonal) can in principle be used for the mesh calibration and will be analyzed for their dependence with respect to the lysimeter filling.

Determination of analog circuit propagation times
For a successful calibration, the sum of analog propagation times t Tx +t Rx needs to be determined (compare (3) and (7)). Equation (7) can be used to determine t Tx +t Rx by conducting, e.g. reciprocal WARR measurements in air between two transceivers, similar to a time-zero correction as described by (2). The intercept of the linear regression then yields the value for t Tx +t Rx . For this purpose, two transceivers can be placed in air to conduct reciprocal measurements with varying transceiver separations before mounting the transceivers to the lysimeter.
We utilize the AWG and oscilloscope to perform a classical time-zero correction to evaluate the measurement accuracy that we can achieve with an ideal measurement setup. Later, measuring t Tx +t Rx can be carried out in the same way. For this purpose, a WARR measurement in air is conducted between Tx6 and a single bowtie antenna, where the single  bowtie antenna is successively moved away from Tx6 (see figure 4(c)). The PVC plate is removed for this measurement and the antenna array is located freely in air using an additional bracket (not shown in figure 4). Both Tx6 and the single bowtie antenna are mounted at a height of about 0.55 m above the ground and have an initial distance of 0.15 m to each other. The single bowtie antenna is then moved away from Tx6 in 0.05 m steps. At each step, one measurement with a stacking of 256 is conducted.
The arrival times that are needed for the calibration are determined as follows: first, the arrival time of the raw signal with the smallest distance to Tx6, i.e. d = 0.15 m (see figure 6 inner box), is determined using a threshold voltage. The threshold value is set to 5γ, where γ is the standard deviation of the noise in the absence of a clear signal. Secondly, the cross-correlation of the signal measured at d = 0.15 m with signals measured for larger d is calculated to determine the arrival times of the corresponding signals (see figure 6, red crosses). The linear regression yields t a (d) = 30.604 + d · 0.2986 −1 . The 99% confidence interval for time-zero is 30.604 ± 0.004 ns. The R 2 value almost approaches one. This indicates a high accuracy of the performed time-zero correction via the WARR measurement. Moreover, the calculated velocity of the signals, v = 0.2986 m ns −1 , shows a deviation of 0.4% compared to the expected value of c = 0.2998 m ns −1 , giving a further indication for a successful time-zero correction measurement. Note that the initial arrival time of the signal measured at d = 0.15 can also be selected differently. This would lead to a new absolute time-zero, but the width of the confidence interval would still be the same. This measurement procedure ensures a highly accurate determination of the sum of analog propagation times t Tx +t Rx .

Simulation of calibration
The standard deviation of the mesh calibration is complex to analytically analyze and is, therefore, analyzed in this section using a simulation following the considerations of section 2.4. Figure 7 shows the nRMS values according to (21) for Mesh2, Mesh4, and Mesh8 in the range of transceivers 1200-1560. The corresponding transceivers 1200 to 1560 are placed in the middle of the lysimeter or mesh, respectively. Figure 7(a) exhibits a clear pattern that separates the nRMS values into 120 × 120 areas. Each area represents one row around the lysimeter. Remember that 120 transceivers are placed around the lysimeter in one row. The last transceiver per row performs a reciprocal measurement only with the left neighbor in the same row (see figure 3(a). It thereby follows that to calculate the trigger offsets between the first and last transceiver per row, 119 intermediate measurements are required compared to e.g. Mesh4 and Mesh8, where one additional measurement is performed between the first and last transceiver per row. This pattern becomes visible in figure 7(a), where the nRMS increases with the distance between first and last transceiver per row. Here, the term distance refers to the number of reciprocal measurements required to connect any two transceivers. After each 120 steps, the nRMS value again decreases due to a direct vertical measurement between the two first transceivers of neighboring rows. The diagonal is zero by definition (t i i = 0). Furthermore, the nRMS values are symmetric around the diagonal, because t i j = −t j i by definition. The overall mean nRMS value for Mesh2 yields 7.70.
For Mesh4 and Mesh8, the observed pattern for the nRMS of Mesh2 is not visible due to further measurements between additional neighboring transceivers (see figures 7(b) and (c)). The largest distance for two transceivers in one row now exists between transceivers on opposite sides of the lysimeter. The striking feature of figures 7(b) and (c) is that in addition to the main diagonal, further diagonal lines are visible. These lines result from further direct vertical measurements, e.g. between transceivers 2 and 122 or transceivers 3 and 123. Thereby, the error between the corresponding trigger offsets is reduced. The diagonal lines repeat every 120 transceivers because, again, this is the number of transceivers per row. For Mesh4, the mean nRMS equals 1.22, indicating that the larger system of linear equations reduces the error by a factor of about 6 compared to Mesh2. For Mesh8, the mean nRMS is 0.93. This result shows that the standard deviation of the trigger offsets for any transceivers can be lower than the initial standard deviation for a direct reciprocal measurement between two transceivers. In addition to the nRMS, we calculate the mean nRMS for a fix transceiver i as  8). For Mesh2, a sawtooth like pattern is observable. This results from the arrangement of reciprocal measurements between neighboring transceivers, because the first transceiver per row has a smaller distance to transceivers in other rows compared to the last transceiver per row, as discussed before for figure 7. For Mesh4 and Mesh8, the sawtooth pattern is removed due to additional reciprocal measurements between neighboring transceivers. The nRMS i values show that transceivers in the upper and lower rows exhibit larger mean nRMS values. This is because transceivers in the top and bottom row perform less reciprocal measurements, i.e. three or five instead of four or eight, respectively, and thereby have an influence on the following rows. Therefore, the most significant jump for the nRMS i values appears after transceiver index 120 and before transceiver index 2880 for Mesh4 and Mesh8, respectively. This should be considered within the FWI, because these In principle, the construction of the system of linear equations as in (12) is not limited to measurements between neighboring antennas as long as reciprocity is guaranteed. As a consequence, one can use the signals through the soil together with (8) to solve the system of linear equations by utilizing an even larger matrix A compared to that of Mesh8. This will be followed up in subsequent work. Nevertheless, the results of Mesh8 are already satisfactory and show that the standard deviation increases with less than √ k · η for increasing distance k between transceivers as it is the case for Mesh2. The usage of a larger linear equation system, therefore, offers significant advantages in terms of error propagation. This result states that expectable random error variances in the low picosecond area are negligible for the calibration approach. Note that using all eight neighboring transceivers of a Tx transceiver as additional Rx does not increase the tomographic measurement time significantly, because the calibration measurements can be included in the actual tomographic measurement.

Impact of lysimeter filling
Equal propagation times t p for coupling signals allow to verify and test the assumption that t Tx + t Rx is equal for all transceivers, as described in subsection 2.5. Therefore, different materials are placed behind the PVC and the impact on the propagation time is analyzed. Before measuring the material's impact on coupling signals, reciprocal measurements are conducted between antennas 6, 7, 10 and 11 (compare figure 4(a)) to assess the general measurement uncertainties. For this purpose, the four antennas are successively used as Tx and Rx by interchanging the cables, resulting in six reciprocal measurements. The time lags for all reciprocal measurements are 4 ± 4 ps, yielding an uncertainty of 4 ps which is in good agreement with the measurement accuracy obtained by the time-zero correction measurements in section 2.7. Note that the time lags can be different for the measurement system due to different hardware. However, the variation of ±4 ps is taken as a reference for an achievable measurement accuracy.
Next, antenna 6 is constantly used as Tx and the remaining 15 antennas are successively used as Rx. For each Tx/Rx combination, two measurements are conducted with either air or water behind the PVC (compare figure 4(c)). This experiment is intended to provide a worst case scenario regarding the permittivity of the material behind the PVC (ε air ≈ 1, ε water ≈ 80). The comparison between air and water coupling is evaluated by the time lags between the signals. The time lags are calculated via the temporal difference between the first local extrema of both signals. Figure 9 shows three examples of received signal pairs (Tx6 is fixed). Comparing the associated signals in figure 9 gives insight into the impact of changing permittivities behind the PVC on the signals. Figure 9(a) shows that when comparing the local minima and maxima of the two signals for Rx7, the time points of similar peaks arrive later for the air signal, which is unexpected, and tend to shift apart for later time points. For Rx10 (see figure 9(b), the peaks of the air signal appear earlier than the peaks of the water signal. For Rx11 (see figure 9(c), more pronounced time lags are visible before and after the maximum peak. In this case, the air signal arrives after the water signal for times earlier than the maximum peak and arrives before the water signal for the following negative peak. It follows that a changing material behind the PVC has an influence on the propagation of coupling signals. However, the influence varies in the course of the signal and is smaller for earlier arrival times for horizontal and vertical signals. To further illustrate this, the differences of first local extrema times for air and water signals are calculated in picoseconds for different Tx/Rx combinations (see figure 10). Here, the time lags indicate a smaller material impact for horizontal and vertical measurements. The impact of water seems to be systematically and should be further analyzed in future work. The horizontal measurements show a time lag between air and water measurements of up to 10 ps. It needs to be tested whether the position of the receivers (Rx5 is located at the edge of the array) has an influence on the signals by, e.g. using a larger antenna array. For vertical measurements, the time lags are already below the determined measurement accuracy of ±4 ps and, therefore, make vertical coupling signals well suited for testing the analog circuit propagation times. Diagonal coupling signals should not be used for the analysis as they show the strongest dependency on the lysimeter filling.
Comparing air and water signals represents a worst case scenario regarding the material permittivities. For the later soil measurements, permittivities vary in a range of ε r = 5 − 40. Therefore, it is reasonable to assume that comparing early signal times and having more realistic material variations behind the PVC reduces the influence of the lysimeter filling onto the coupling signals even more. Then, horizontal coupling signals might be used for testing the analog circuit propagation times in addition to vertical coupling signals. It follows that the assumption that t Tx + t Rx is equal for all transceivers can be tested, as the influence of the lysimeter filling onto coupling signals is negligible.

Conclusion
This work presents a time-zero calibration strategy for a monitoring system with about 3000 transceivers that are placed around a soil sample. Because the system is installed after the filling of the lysimeter, the calibration cannot be performed with known material, e.g. air, between any two antennas. The transceivers are assumed to be equal, while the trigger lines to control each transceiver can be of different length. Therefore, each possible Tx/Rx combination requires a separate timezero correction. Our approach utilizes the ability of the system to use every antenna as Tx and Rx, and, thereby, allows to conduct reciprocal measurements between any two antennas. This reciprocal measurement arrangement makes it possible to reduce the number of unknowns and enables the measurement of the propagation time of the signal between two antennas directly and uniquely regardless of the position of the antennas. We term this procedure pairwise calibration. Furthermore, reciprocal measurements between transceivers can be combined within a superposition to calculate time-zero for any Tx/Rx combination, termed as mesh calibration. Solving a sufficiently large system of linear equations leads to the fact that the standard deviation for time-zero increases significantly slower than with √ k for an increasing distance k between two transceivers. For realizable jitter values in the low picosecond range, these time inaccuracies can, therefore, be neglected. Also, the mesh calibration can be applied before measuring fast soil processes and, thereby, reduces the measurement time compared to the pairwise calibration.
The described methods require equal or known signal propagation times within the transceivers. In following works, it must be verified whether and with what reliability this condition is fulfilled. The pros and cons of the pairwise and mesh calibration are summarized in table 1.
In order to test the prerequisite that the signal propagation times within the transceivers are equal (or known and stable), the influence of the material behind the PVC needs to be negligible for coupling signals. We found that the influence of the material behind the PVC, or inside the lysimeter, is minimal for horizontal and vertical measurements. Therefore, the coupling signals between horizontally and vertically neighboring transceivers have the potential to control the assumption that t Tx + t Rx is equal for all transceivers. The diagonal neighbors are significantly influenced by the soil and should not be used for this purpose. This result enables the use of a flexible mesh for the mesh calibration, so whenever a violation of the assumption arises, the corresponding transceivers are excluded from the inversion. Further, simulations should be applied to gain a better understanding of the field distribution between neighboring antennas and especially between diagonal neighbors and why these are especially influenced by the soil.
We assume that the time delay within the identical transceivers t Tx + t Rx is constant over time and for all transceivers, and, therefore, can be calibrated by a single WARR measurement. The time-zero correction measurement yields an accuracy of 4 ps for the reference measurement. The same accuracy is expected for the trigger lines. The final time-zero error is about √ 4 2 + 4 2 ≈ 6 ps for the pairwise calibration, which is significantly better than the required 25 ps. Note that this analysis only addresses statistical variations. The magnitude of systematic shifts due to thermal effects and transceiver variations in the analog circuits will be investigated in a following study. However, the methods and results presented in this paper provide an important contribution towards an automated in situ time-zero correction of the designed monitoring system. The presented approaches are not limited to the lysimeter application. In principle, the methods can be used by all measurement systems that support reciprocal measurements such that no distinct time-zero measurement setup is needed for calibration. For example, the approach could also be used in multichannel surface measurements, eliminating the need for time-consuming repeated calibration measurements of the system. In addition, our approach enables tomographic measurements via radar signals for any application where manual recalibration is not feasible.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.