## 2 Materials and Methods
### GPR Data Processing and Solid Modeling  
A Ground-Penetrating Radar (GPR) output typically comprises a radargram, which consists of time-series reflections that provide cross-sectional profiles of subsurface depth over time [43]. The collected radargrams were processed using Reflexw software version 4.0 (Sandmeier Scientific Software, Karlsruhe, Germany). The processing protocol included the following: (i) trace interpolation to ensure equidistant spacing between markers at survey line intersections, with an increment of 0.2 m, (ii) static time shift adjustment to synchronize the arrival of direct ground waves to 0 ns, (iii) bandpass filtering within the temporal domain, (iv) application of an exponential gain function to counteract GPR signal attenuation with depth; as per [PERSON] et al. [10], the same gain curve, established from pre-wetting data, was uniformly applied across all time-lapse data, (v) implementation of a background removal filter to mitigate horizontal noise, and (vi) compression in both temporal and spatial dimensions to minimize computational time required for subsequent processing steps, such as 3D interpolation.  
Amplitude values, G, from all radargrams within the grid were aggregated to form a comprehensive XYZG dataset, where X, Y, and Z denote spatial coordinates (easting, northing, and elevation). Elevation, Z, was derived from wave travel time recorded during GPR acquisition, with an assumed wave velocity of 0.11 m ns\({}^{-1}\), as determined from the detection of subsurface objects at known depths. Consequently, both pre- and post-wetting XYZG datasets were obtained. Subsequently, additional datasets were created based on the absolute differences between pre- and post-wetting G values. The movement of water within the unsaturated zone induced variations in the dielectric contrast between layers, thereby altering the reflection coefficient [10]. These alterations were manifested as changes in amplitude. Therefore, the differenced datasets were constructed to emphasize amplitude fluctuations between repeated GPR radargrams acquired along the same survey lines before and after the infiltration test, e.g., [23; 44]. It was hypothesized that the absolute difference between pre- and post-wetting amplitude values was inversely correlated with the soil water pressure head [45]. Consequently, the most significant reflection differences were anticipated at the location of the saturated soil bulb extending outward from the infiltration surface, with decreasing reflection differences in the partially saturated zone surrounding the wetting bulb. In this zone, the water pressure head diminishes as the wetting front advances due to soil capillarity from the saturated bulb [37].  
The RockWorks 17 software [46] was employed in conjunction with the inverse-distance anisotropic modeling technique to execute a 3D interpolation of the differenced GPR datasets. This anisotropic search methodology enhances the interpolation accuracy of voxel (volumetric picture element) values located between clusters of data points and has been effectively utilized in numerous geophysical studies, e.g., [47; 48; 49]. The interpolation algorithm assigned values to voxel based on a weighted average of adjacent data points within each 90-degree sector around the node. The weighting of each data point's amplitude value was inversely proportional to the square of its distance from the voxel node (inverse-distance squared).