Dynamic groundwater recharge simulations based on cosmic‐ray neutron sensing in a tropical wet experimental basin

Although cosmic‐ray neutron sensing (CRNS) is probably the most promising noninvasive proximal soil moisture measurement technique at the field scale, its application for hydrological simulations remains underexplored in the literature so far. This study assessed the use of CRNS to inversely calibrate soil hydraulic parameters at the intermediate field scale to simulate the groundwater recharge rates at a daily timescale. The study was conducted for two contrasting hydrological years at the Guaraíra experimental basin, Brazil, a 5.84‐km², a tropical wet and rather flat landscape covered by secondary Atlantic forest. As a consequence of the low altitude and proximity to the equator low neutron count rates could be expected, reducing the precision of CRNS while constituting unexplored and challenging conditions for CRNS applications. Inverse calibration for groundwater recharge rates was used based on CRNS or point‐scale soil moisture data. The CRNS‐derived retention curve and saturated hydraulic conductivity were consistent with the literature and locally performed slug tests. Simulated groundwater recharge rates ranged from 60 to 470 mm yr–1, corresponding to 5 and 29% of rainfall, and correlated well with estimates based on water table fluctuations. In contrast, the estimated results based on inversive point‐scale datasets were not in alignment with measured water table fluctuations. The better performance of CRNS‐based estimations of field‐scale hydrological variables, especially groundwater recharge, demonstrated its clear advantages over traditional invasive point‐scale techniques. Finally, the study proved the ability of CRNS as practicable in low altitude, tropical wet areas, thus encouraging its adoption for water resources monitoring and management.

as a natural or artificial mechanism of groundwater storage renewal (Gleeson et al., 2012). The groundwater recharge connects the atmospheric, surface, and subsurface components of the water balance (Mohan et al., 2018), being mainly controlled by the precipitation amount and intensity, boundary layer climatology, topography, water table level, watershed geomorphology, soil and vegetation characteristics, and irrigation return flow (Jasechko et al., 2014;Moeck et al., 2020). Understanding the seasonal controls upon recharge requires its accurate estimation (Jasechko et al., 2014), which can only be obtained by experimental methods or modelling since the groundwater recharge cannot be directly measured (Melo et al., 2015).
A considerable range of experimental methods and models is available for groundwater recharge estimation at various spatiotemporal scales, whose correct selection depends on data availability, desired spatiotemporal resolution, and result representations (Walker et al., 2019). Amongst the existing approaches, the following methods are frequently used to estimate groundwater recharge: (a) monitoring of groundwater level in unconfined aquifers, such as the water table fluctuation (WTF) method (Cai & Ofterdinger, 2016;Wendland et al., 2007); (b) tracer techniques, which estimate aquifer renewal through the concentrations of chemical elements and substances in the water, such as the chloride mass-balance method (Brunner et al., 2004;Hornero et al., 2016); (c) the water balance method, which uses the main components of the water cycle as inputs and outputs of the system (Hornero et al., 2016;Wendland et al., 2007); Darcy's law application, which uses the hydraulic gradient and conductivity to calculate the soil water percolation rates (Callahan et al., 2012;Yin et al., 2011); and (d) numerical modeling, which consists of a mathematical representation of the groundwater processes (Melo et al., 2015;Šimůnek, 2015).
The vertical percolation rate from the numerical simulation of the vadose zone is often considered as a proxy of groundwater recharge (Mathias et al., 2017). The simulation of these groundwater fluxes, however, come with the necessity of the proper calibration of some soil hydraulic parameters by means of soil moisture observations (Graham et al., 2018).
Three different basic approaches are currently available for retrieving soil moisture information: in situ measurements, remote sensing data, and proximal soil sensing. In situ measurements are generally the most accurate and direct methods for collecting soil moisture data but provide point-scale information only, whereas spatially extended information demands time-consuming and burdensome efforts, as the case for wireless sensor networks (Bogena et al., 2010) or spatial timedomain reflectometry (TDR) probes (Graeff et al., 2010). Satellite-borne remote sensing products feature a large-scale estimation of soil moisture (e.g., through active microwave sensors), but usually provide near-surface or root-zone information with coarse spatiotemporal resolutions (Famiglietti

Core Ideas
• Cosmic-ray neutron sensing (CRNS) can detect field-scale soil moisture also at low latitudes. • Soil hydraulic parameters were calibrated using CRNS and inverse simulations. • Daily groundwater recharge rates were simulated using CRNS-based soil hydraulic parameters. • CRNS-based recharge was matching independent area average recharge estimates. et al., 2015). Lastly, proximal soil sensing methods became a promising alternative method to monitor soil moisture noninvasively and at the field scale Robinson et al., 2008) based on the relationship from targeted relevant subsurface variables (e.g., groundwater table) to soil properties that could relatively easily be mapped from the surface. Amongst them, the cosmic-ray neutron sensing (CRNS) method (Zreda et al., 2008) is currently the most promising technique, closing the inconvenient spatial scale gap between invasive sensor data and satellite-based remote sensing products. The CRNS uses the inverse relationship between aboveground neutron count rate and all close-by hydrogen content pools (Desilets et al., 2010;Hendrick & Edge, 1966). It continuously measures average soil moisture from about 130to-240-m effective radius and from about 15-to-55-cm penetration depth for wet and dry conditions, respectively (Köhli et al., 2015;Schrön et al., 2017). Another advantage of CRNS includes the calibration dependence on a single site-specific parameter through a standard curve. During the last decade, several studies applied and progressed the CRNS method on stationary and mobile platforms up to the scale of square kilometers (Fersch et al., 2020;Schrön, Rosolem, et al., 2018) and by monitoring stations installed in a broad variety of climate conditions-namely, continental (Baatz et al., 2014), temperate (Evans et al., 2016), semiarid (Zreda et al., 2012), and tropical (Hawdon et al., 2014). The CRNS was further successfully applied in vegetated areas such as cropped fields (Franz, reda, Rosolem, et al., 2013;Rivera Villarreyes et al., 2011) and forests (Bogena et al., 2013;Heidbüchel et al., 2016). Most CNRS studies presented good performances when compared with soil moisture datasets from gravimetric experiments, invasive pointscale probes, and land surface models. The advantages of the CNRS method have boosted its application  in various scientific topics, like runoff dynamics (Dimitrova-Petrova et al., 2020), snow dynamics (Schattan et al., 2017), canopy water interception (Baroni & Oswald, 2015), vegetation biomass calculation , remotely sensed products validation (Montzka et al., 2017), land surface modeling (Baatz et al., 2017), agricultural support and irrigation (Finkenbiner et al., 2018;Li et al., 2019), and urban effects .
Although there is a variety of CRNS studies, the assessment of its performance and application remains poorly explored in regions that gather unfavorable environmental conditions for CRNS operation (Bogena et al., 2013). One of the more unfavorable conditions can be found in tropical regions that often exhibit: (a) a high geomagnetic cutoff rigidity, which reduces the secondary cosmic-ray intensities and thus incoming neutron intensities (Andreasen et al., 2017); (b) the tropical rainforest, whose higher atmospheric humidity and biomass content decrease the possibilities of aboveground neutron detection Jakobi et al., 2018;Köhli et al., 2021); and (c) high soil moisture values shrinking the CRNS support volume (penetration depth and horizontal footprint). Due to these unfavorable conditions, the CRNS data have been very rarely studied in tropical regions, an exemption being (Hawdon et al., 2014) in Northern Australia, whose information has been used for space-borne satellite product validation (Duygu & Akyürek, 2019;Montzka et al., 2017) and model calibration . The potential of CRNS for deriving also groundwater recharge rates has been suggested (Baroni et al., 2012), but not specifically targeted so far.
Few studies have used CRNS data to derive soil hydraulic parameters in combination with orthogonal functions (Finkenbiner et al., 2018;, in inverse modeling (Rivera Villarreyes et al., 2014), and also in combination with neutron transport models (Baatz et al., 2017). The COsmic-ray Soil Moisture Interaction Code (COSMIC) developed by Shuttleworth et al. (2013) provided a simplified physically based analytical approach that reduces the computational cost of neutron transport simulations for the inverse calibration of soil hydraulic parameters. Recently, it was efficiently coupled with the one-dimensional soil hydrological model HYDRUS-1D, enabling its use for vadose zone simulations at intermediate field scales (Brunetti et al., 2019).
Based on this knowledge, this study assesses the use of noninvasive CRNS data for two contrasting hydrological years to inversely calibrate the soil hydraulic parameters at intermediate field scale and, with this information, simulating the dynamic groundwater recharge rates through the physically based soil hydrological model HYDRUS-1D (Šimůnek et al., 2016). The hydrological simulation of quantities in the vadose zone based on the CRNS data (i.e., soil moisture, soil hydraulic properties, and groundwater recharge) were tested against additional measurements and compared with the results obtained using point-scale conventional time series (e.g., TDR probes). The experimental site of this study holds unfavorable conditions for CRNS operation with its location close to the equator, low-altitude tropical wet climate, and shallow groundwater table and thus fills a gap within CRNS studies worldwide. In this study, we hypothesized that the use of CRNS data for estimating hydrological variables, also in such regions, will result in a higher prediction accuracy than the traditionally applied invasive point-scale techniques.

Study area and monitoring network
This study was carried out in the Guaraíra experimental basin (GEB, Figure 1), located between the coordinates 07˚16′35"-07˚19′11″ S and 35˚01′18"-35˚02′27″ W, which is operated as part of the Brazilian Hydrology Network for the Semi-Arid Region (REHISA) since 2003. The GEB is located at ∼25 km from the shoreline from northeastern Brazil, situated on the surrounding rural area of João Pessoa city (Paraíba state's capital). The GEB comprises a total area of 5.84 km 2 with a mean elevation of 100 m asl and a rather flat landscape with a mean slope equal to 1.8%. The GEB lies within the Gramame river basin that is the topic of recent investigations, also in respect to larger-scale estimations of groundwater recharge from satellite remote sensing (Barbosa, 2020). The Gramame river basin has a mean potential groundwater availability of ∼1.07 × 10 8 m 3 yr −1 , whereas the cumulative withdrawals sum up already to 83.9%, much larger than the maximum recommended exploitation of 60% (PERH-PB, 2006). The groundwater recharge in the region is mainly promoted by vertical infiltration (i.e., diffuse recharge) over these sedimentary terrains.
The GEB features a tropical wet climate with dry summer and well-distinguished rainy season during the austral autumn and winter (As), according to the Köppen climate classification (Alvares et al., 2013). The rainfall in the region is mainly formed from shallow convective clouds, with a mean annual rainfall of 1,700 mm, ∼70% of which occurs from March to July. The air temperature and potential evapotranspiration are high in the experimental basin, with mean annual values of approximately 25.4˚C and 1,500 mm, respectively. The predominant land use and land cover in GEB are forest (native and secondary Atlantic remnants), cropland (sugarcane [Saccharum spp.] and pineapple [Ananas comosus (L.) Merr.]), and built-up areas (rural housing and unpaved roads). The soil in GEB originated from crystalline basement erosion during the Late Quaternary, which constitutes massive deposits of friable sands (Rossetti et al., 2012). The soil texture in the experimental basin is mainly composed of coarse sand (∼46%), fine sand (∼52%), and silt and clay (∼2%), belonging to the haplic acrisol and carbic podzol soil types encompassing the study area (Barbosa et al., 2019;IUSS Working Group WRB, 2015).
The study used data from two monitoring sites in GEB (S1 and S2, Figure 1), the vegetation at these sites being F I G U R E 1 Representation of (a) Brazil (BR), (b) Paraíba and Gramame river basins with João Pessoa city, (c) Guaraíra experimental basin, (d) location S1 with a cosmic-ray neutron sensing probe, and (e) location S2 with a meteorological station including a time-domain reflectometry soil moisture profile sparse tropical forest of 3-to-4-m height in secondary-growth stage with no substantial litter layer to observe. The sites were equipped with rain gauges, meteorological sensors (e.g., radiometer, anemometer, and thermo-hygrometer), invasive and noninvasive soil moisture probes, and groundwater monitoring wells. One noninvasive CRNS probe was installed at S1 for recording hourly soil moisture values. On its surrounding, three soil profile probes type PR2 (Delta-T Device) were vertically inserted into the soil at radial distances of 17, 66, and 153 m from the CRNS probe, for monitoring the hourly soil moisture in the depths of near-surface in 10, 20, 30, 40, 60, and 100 cm. The generalized polynomial calibration curve optimized by the PR2 manufacturer was used, as it covers a wide range of mineral soil types with organic content of ∼1% (Delta-T, 2016), which is the case at GEB. Hardware problems have made the measured data in depths >10 cm unreliable for two profile probes, whereas only the one at 66-m distance to the CRNS sensor delivered reliable values for all depths during the entire observation period. For testing the CRNS derived soil moisture, the soil moisture time series monitored in the near surface (10-cm depth) by the three PR2 profile probes were smoothed through a moving average and integrated to an area average soil moisture using the weighting function (Schrön et al., 2017). Three invasive soil moisture sensors (Campbell Scientific, CS616 water content reflectometer) were horizontally inserted in the soil at S2 to monitor a profile of soil moisture in the depths of 5, 20, and 40 cm, recorded also at hourly timescale. The measurement principle is based on a transmission line oscillator and works similar to TDR systems (Blonquist et al., 2005) and for simplicity the sensors will be called TDR probes within this study. Groundwater levels were recorded using fast-response pressure transducers. The groundwater table depth in GEB showed a mean value of ∼157 cm below the surface, influenced by the shallow occurrence of a Si-cemented impermeable layer in the mean depth of ∼200 cm below the surface. All the datasets were averaged in this study to values at daily timescale. We assume similar hydrologic behavior at the sites S1 and S2, based on the measured data on texture and bulk density, which is confirmed by measurements from Santos et al. (2012), who reported extensive measurements of soil characteristics and hydrological properties within the catchment (raw data can be obtained from the appendix of Silva (2009). Figure 2 shows the usage of the data within this study at each field site and comparisons drawn between datasets and field sites. Firstly, the CRNS calibration based on gravimetric soil sampling and a comparison with independent soil moisture profile measurements (PR2) is intended to test the usability of CRNS as soil moisture monitoring method in humid tropic conditions. Secondly, after this basic validation of CRNS applicability, the time series of corrected neutron count rates is used directly to inversely estimate the soil hydraulic F I G U R E 2 Workflow (black arrows) showing the different measurement setups at the two locations used within this study (S1 = top and S2 = bottom). Comparisons between data products are marked with a gray, two-sided arrow. A device for cosmic-ray neutron sensing (CRNS) was calibrated at S1 and its applicability in humid tropic conditions was validated by independent soil moisture profile measurements (PR2). Soil hydraulic parameters (SHPs) were optimized independently within HYDRUS at each field site (CRNS and Cosmic Operator at S1, using a time-domain reflectometry [TDR] soil moisture probes at S2), derived parameters were compared with values from literature and a slug test validated the saturated hydraulic conductivity (K s ) at S1. The CRNS-based simulated soil moisture time series at S1 were validated by independent measurements (PR2). Based on the simulated fluxes at S1 and at S2, at each site a time series of groundwater recharge estimates was obtained and compared to independent estimates derived from water table fluctuation (WTF) measured in groundwater wells parameters within HYDRUS 1D coupled with the COSMIC operator (Brunetti et al., 2019). The resulting fluxes due to the nature of the CRNS data represent effective field-scale estimates. Thirdly, a standard approach for inverse calibration of the hydraulic parameters using soil moisture time series (TDR) is applied at S2, representing point-scale parameters. Fourthly, the calibrated parameter sets of both locations were compared with data reported in literature (Santos et al., 2012) reflecting the variability observed within the catchment. Saturated hydraulic conductivity based on CRNS data was independently validated with slug tests at S1. Moreover, the simulated soil moisture time series resulting from the CRNSoptimized parameter set at S1 is compared with the independent profile soil moisture of the PR2. Finally, the simulated vertical percolation (constituting groundwater recharge) was derived at S1 based on the area average CRNS data and at S2 based on the point-scale TDR data. Both recharge estimates were compared with recharge estimates using fixed-interval WTF from groundwater wells at each site.

CRNS data
In this study, a CRNS probe equipped with a 3 He neutron counter (manufactured by Hydroinnova, type CRS 1000) was installed at a location seen as typical for GEB conditions to record the neutron count rate at hourly intervals as cumulative neutron counts. Standard correction procedures for atmospheric pressure, incoming neutron count rate and air humidity were applied to the raw neutron count rate (see Appendix A for details). The theoretical relationship between corrected neutron count rate N' (in counts per hour [cph]) and soil moisture depends on the dry neutron count rate parameter N 0 (in cph) (Desilets et al., 2010;Equation 1). This parameter was calibrated based on soil moisture measured by the gravimetric method on soil samples collected at 17 locations ( Figure 3) in three depth intervals (0-5, 5-10, and 15-20 cm), so that different moisture conditions in the CRNS footprint were covered. The weighted average soil moisture using the weighting Vadose Zone Journal F I G U R E 3 Calibration of the cosmic-ray neutron sensing (CRNS). Distribution of (a) the vertically weighted soil moisture around the CRNS probe retrieved from gravimetric sampling, the locations of the PR2 profile probes, and (b) the average water equivalent from soil organic matter and lattice water content per depth. A medium-size footprint of 180-m radius is indicated here for illustration, but it varies between dry and moist environmental conditions to some extent function proposed by Köhli et al. (2015) was then calculated and used to calibrate the CRNS standard calibration curve. For calculating N 0 , the corrected neutron count rate was averaged during the calibration period and used alongside with the weighted average soil moisture and the total bound water content. The latter component corresponds to the sum of soil organic matter and lattice water contents (θ SOM and θ LW , in cm 3 cm −3 ), which were calculated in this study according to Heidbüchel et al. (2016), by weighting soil samples after drying at 105˚C for 24 h (m 105 ), and burning 10 g of dry samples at 400˚C for 16 h (m 400 ) and 1,000˚C for 12 h (m 1,000 ) (Equations 2 and 3). The soil bulk density (ρ bd ) was obtained by averaging the gravimetric-based estimates, and the density of water (ρ w ) was adopted as 1 g cm −3 .
The effect of biomass water equivalent was neglected in this study, because forest biomass and litter layer are much smaller than in other CRNS studies (Bogena et al., 2013;Heistermann et al., 2021) and the effect of a tree canopy as an additional moderator and the creation of neutrons with detectable energy ranges are still under discussion (Baroni & Oswald, 2015;Li et al., 2019). The CRNS-based soil moisture determined for hourly integration periods was smoothed to a 24-h timescale (i.e., representing daily values), using the moving average calculation to reduce the stochastic measurement noise caused by low incoming neutron count rate (Bogena et al., 2013). This is longer than in other studies but related to the high cutoff rigidity close to the equator.

Soil hydrological modelling and estimation of soil hydraulic parameters
Understanding of how soil moisture can be affected by land surface water fluxes allows for comprehensively modeling the groundwater recharge (Sadeghi et al., 2018). Simulations of soil moisture and land surface water fluxes were carried out using the HYDRUS-1D model because of the large number of implemented processes, the flexibility for application in different soil water conditions, and the possibility of inverse calibration (Šimůnek et al., 2016). The recent integration of the COSMIC operator enables the proper handling of CRNS data in this context (Brunetti et al., 2019).
The HYDRUS-1D model is capable of solving the governing flow equations for partially-saturated porous media through Richards' nonlinear partial differential equation (Equation 4; Šimůnek et al., 2008).
where θ is the soil moisture, h is the water pressure head (or matric potential), t is time, S(h) is the sink term, K(h) is the unsaturated hydraulic conductivity function, and β is the angle between the flow direction and the vertical axis (= 0˚for vertical flow, in this study). The lateral flux was neglected in this study because of the high permeability expected from the thick loose sand of the study area, which impairs the modification of the vertical percolation to an inclined direction (Jankowski, 2014). The sink term of Richards' equation refers to the actual root water uptake (RWU a , for details see Appendix B). A rooting depth of 0.70 m was determined by in situ observations. Variable hydro-phenomenological effects were disregarded on the calculation of RWU in GEB due to the negligible temporal variations observed in situ on the rooting depth and root density. The soil water output is mainly triggered by potential evapotranspiration ET p , which was estimated using the Penman equation (Penman, 1948) modified by (Shuttleworth, 1993). Beer's equation was then used to split the ET p into potential soil surface evaporation E p and the transpiration T p by the Equations 5 and 6 (Sutanto et al., 2012).
These potential components were calculated based on the leaf area index (LAI) (Bastiaanssen et al., 1998;Huete, 1988), which was averaged over GEB using two Moderate Resolution Imaging Spectroradiometer (MODIS) reflectance products (MOD09Q1 and MYD09Q1). The LAI estimate is based on the main land cover of GEB, which was Atlantic secondary forest. The cloud cover effects within the remote sensing products were attenuated by carrying out monthly map compositions by substituting the pixels with disparate LAI values by pixels of other images in the same month (Barbosa, 2020). These equations use the light extinction coefficient k that, in this study, was set equal to 0.463 due to the arborealshrub physiognomies of the vegetation that features a generally rounded canopy (Chen et al., 2014;Raoufi & Beighley, 2017), whose average height of 330 cm was determined by in situ observations. The component T p was then used to estimate RWU p by subtracting the interception I, which was calculated by an empirical equation presented by Šimůnek et al. (2013).
The actual evaporation E a and the infiltration F were calculated by the numerical solution proposed by Neuman et al. (1974), which depends on the upper boundary condition. In this study, the upper boundary was set as the atmospheric condition with maximum surface water layer of 2.5 cm similar to (Batalha et al., 2018), which was observed in situ after ponding occurrences. This information was set to obtain the actual net surface water fluxes S a , calculated as the F minus E a . The resulting fluxes were accumulated at daily time scale for the hydrological study years and plotted alongside with their annual total proportions in relation to the potential net surface water fluxes S p , which refers to F minus E p .
The simulation of water fluxes depends on the proper calibration of soil hydraulic properties (retention and conductivity) and is crucial for consistently simulating the vertical percolation rates under shallow groundwater conditions (Wessolek et al., 2008), such as in GEB. A study comparing four parametrization methods showed that the inverse calibration based on field soil moisture data implemented in the HYDRUS-1D model provided the best water flux predictions (Graham et al., 2018). Therefore, HYDRUS-1D was used in this study for calculating the matric-potentialdependent soil water retention [i.e., θ(h)] and hydraulic conductivity [i.e., K(h)] through the Mualem-van Genuchten soil hydraulic model (van Genuchten, 1980), including the option of air-entry detailed by Vogel et al. (2000) (Equations 7-12): Vadose Zone Journal where in which = (1 − 1∕ ) , and n > 1. The α is the air entry suction factor, n is the pore size distribution parameter, l is the pore tortuosity and connectivity lumped parameter, θ r is the residual soil moisture, θ s is the saturated soil moisture, and K s is the saturated hydraulic conductivity. This model reduces to the conventional van Genuchten-Mualem equation (van Genuchten, 1980), in the case of h s = 0 (at site S2), leading to θ m = θ s . The air entry (h s ) of −2 cm related to factor α was considered (at site S1) for CRNS based calibration to reduce the nonlinearity effect of the hydraulic conductivity function close to saturated conditions (Ries et al., 2015;Vogel et al., 2000). The l factor was assumed as 0.50, as usually adopted in the literature, because of its underexplored effects on simulations (Mualem, 1976;Ventrella et al., 2019;Wang et al., 2009). Amongst the soil hydraulic parameters, the saturated hydraulic conductivity was specifically assessed in this study, since it affects the vertical percolation rates and, in turn, the groundwater recharge and was shown in previous studies to inhibit a high level of uncertainty (Brunetti et al., 2019;. The hydraulic conductivity was experimentally estimated from water-level displacement datasets resulting from the overdamped injection and withdrawal of a slug in a monitoring well at S1. Such a slug test is expressed by an empirical logarithmic equation that considers: (a) the well, perforation, and effective radii; (b) the well filter length and saturation thickness; and (c) the groundwater table dynamic (Bouwer & Rice, 1976). Six slug-in and slugout attempts were carried out since the adjusted water level showed an uneven recovering to the static level; therefore, the groundwater level was monitored using a fast-response pressure transducer.
The convergence of the soil hydraulic parameters also relies on the proper input of their initial values for the inverse calibration. Thus, the Rosetta Lite v.1.1 module (Schaap et al., 2001), incorporated in the HYDRUS-1D model (Šimůnek et al., 2008), was used to calculate the initial values of the soil hydraulic parameters from soil texture and hydrophysical data. The vertical soil texture down to 100-cm depth was calculated from the mean soil particle distributions determined by sieving (particle size > 0.063 mm) and sedimentation methods (particle size ≤ 0.063 mm). The hydrophysical data were consistent with soil water contents experimentally obtained in pressured chambers for the suctions at 33 and 1,500 kPa.
In this study, all soil hydraulic parameters except for l were calibrated using CRNS neutron count rate information at S1 and using a profile by TDR soil moisture at S2. For handling the CRNS data within the inverse calibration procedure, the physically based and analytic forward neutron operator COS-MIC was used . It simulateds a neutron count rate based on soil moisture information. Using measured CRNS neutron count rates thus enables to calibrate the effective soil hydraulic parameters at the intermediate scale. COSMIC is capable of representing the (a) exponential reduction of high energy neutron flux towards deeper soil layers, (b) fast neutron generation from high energy neutrons at different depths, and (c) additional depth-dependent reduction of lower energy neutrons before their detection at the soil surface. The COSMIC operator includes the siteindependent parameters L 1 = 162.0 g cm −2 , L 2 = 129.1 g cm −2 , and L 4 = 3.16 g cm −2 . The COSMIC operator also includes the site-dependent parameters L 3 (g cm −2 ) and α (cm −3 g −1 ), which are calculated based on the soil bulk density that was averaged in this study down to 20-cm depth and the total bound water content (θ SOM + θ LW ). The soil surface high-energy neutron count N C scales the simulated count rate to the observed count rate and was set in this study to 90 by a trial-and-error approach with forward simulations of the HYDRUS-COSMIC model (based on initial estimates of soil hydraulic properties) to derive lowest deviations to the observed neutron count rate. It was additionally validated by the relationship of the N 0 calibration parameter derived by Baatz et al. (2014) and showed reasonable agreement (data not shown). This leaves the same hydraulic parameters to calibrate for using CRNS data or TDR profiles.

Groundwater recharge fluxes
The dynamic vertical percolation rates were simulated for the two contrasting hydrological years 2017-2018 (dry) and 2018-2019 (rainy) by computing the cumulative water fluxes in the depths of 20, 30, 40, 60, and 100 cm, and the bottom boundary. These simulations were performed at an daily timescale to avoid the bias caused by the temporal averaging of the meteorological data (Batalha et al., 2018). The dynamic groundwater recharge rates correspond to the vertical percolation rates in the bottom boundary, which was set as the variable groundwater level due to the shallow occurrence of a Sicemented impermeable layer in 194 cm at S1 and in 215 cm at S2. In this study, a homogeneous soil sandy profile discretised to 389 nodes at S1 and 431 nodes at S2. Homogeneity can be assumed because of the low variation in texture (>96% sand) and organic matter with depth (>2%). The groundwater recharge (R) was calculated by integrating the water flux over time, applying a water budget equation (Equation 13) that includes precipitation (P), l, RWU a , E a , runoff (Q), and soil water storage changes (ΔS). The runoff was only observed in few stretches of GEB, wherein flat natural landscape was modified to install gas pipelines or at the more sloped areas closer to the outlet, leading to negligible results found in simulations that were disregarded in the analyses of this study.
The dynamic of groundwater recharge rates simulated by HYDRUS-1D model was validated by assessing the cumulative groundwater recharge rates obtained from experimental estimates calculated by the WTF method (Equation 14). The WTF method was applied for the two monitoring wells present at S1 and S2. For such purpose, an alternative solution was considered due to daily timescale of simulation, which consists of a fixed-interval WTF method that combines two approaches: the so-called RISE and the Master Recession Curve (MRC) procedures (Nimmo et al., 2015). In this study, the RISE and the MRC procedures were used to determine the contributions of the actual groundwater table rise and the potential groundwater decline, respectively. Moreover, the MRC was defined in this study by linear regression over the longer groundwater recession period, differently from Delin et al. (2007). These approaches involve less subjectivity by circumventing the lack of user's judgment expertise since it adopts a fixed-interval frequency in the accurate measurement of groundwater level (Nimmo et al., 2015).
where R WTF denotes the groundwater recharge estimated using the WTF method, S y represents the aquifer specific yield coefficient, Δt is the time from the beginning of the rise to the peak, and ΔH corresponds to the sum of the actual groundwater rise and the potential groundwater decline for the same period, the latter being obtained by extrapolating the recession curve (Jie et al., 2011). The recession curve is the behavior that the groundwater level would have if no elevation in the groundwater level occurred. The S y values of 0.10, 0.16, and 0.24 obtained by pumping tests over the region were tested and analyzed for comparison purposes since the S y can change sensibility according to the depth to the groundwater table (Dettmann & Bechtold, 2016). Values of S y ranging from 0.10 to 0.26 were similarly reported in other unconfined aquifers of northeastern Brazil (Vasconcelos et al., 2010).

Performance metrics
In this study, four statistic metrics were selected to assess the performance of the hydrological simulations. The correlation coefficient (CC,Equation 15) and Kling-Gupta efficiency (KGE,Equation 16) were used to describe the relationship and agreement between simulated and observed datasets, respectively. The mean absolute percentage error (MAPE, Equation 17) and the percent bias (PBIAS, Equation 18) were used to describe the error and bias between simulated and observed datasets, respectively. Percentage values of these two metrics were adopted so that a comparative discussion with the literature results can be performed.
where N is the length of a dataset,̄and̄are the mean values of the simulated (S i ) and observed (O i ) datasets, respectively, and σ S and σ O are the standard deviations of the S i and O i datasets, respectively.

Calibration of a CRNS sensor under humid tropic conditions
The first set of analyses concerned the calibration of the CRNS, which was operated at S1 to estimate the soil moisture at the intermediate field scale from July 2018 to June 2019 (Figure 3). Gravimetric sampling was performed on 1 Aug. 2018, resulting in the weighted average values of soil moisture and soil bulk density equal to 0.088 cm 3 cm −3 and 1.11 g cm −3 , respectively, following the weighting procedure outlined in Schrön et al. (2017). Surrounding the CRNS, the soil moisture varied from 0.056 to 0.296 cm 3 cm −3 . The spatial variation was likely influenced by the uneven vegetation cover, implying also an average soil bulk density smaller than that found at four points within the experimental catchment by Barbosa et al. (2019). Additionally, the total bound water content summed up 0.009 cm 3 cm −3 in the near-surface, which, alongside with soil moisture and soil bulk density, led the dry neutron count rate parameter N 0 to be equal to 550 cph.
Although the total bound water content was relatively small, it was also considered in the CRNS calibration curve for refining the soil moisture estimation, as demonstrated by Baroni et al. (2018).
The resulting daily soil moisture time series consistently responded to the rainfall events and was able to detect the influence of the groundwater level fluctuations during the recession periods (Figure 4). This information was then compared with the PR2 soil moisture time series (weighted area average value using measurements <10-cm depth), showing similar behaviors throughout most of the monitoring period except for the dry season where a slight shift was found. Other studies recommend to calibrate using more than one calibration data set, preferably at different soil moisture stages (Iwema et al., 2015), to improve the calibration result. However, they also adjusted the parameters of the calibration curve to match the independently observed soil moisture. When staying with the standard values of the parameters, deviations in the very wet and dry range of soil moisture were also found by other studies (Baatz et al., 2014;Heidbüchel et al., 2016). Calibration of the parameters and local adjustment of the F I G U R E 5 Comparison among the (a) retention curves and (b) hydraulic conductivity curves based on (c) the van Genuchten-Mualem parameters (values shown are rounded to the last digit), obtained from time-domain reflectometry (TDR), cosmic-ray neutron sensing (CRNS), and field sampling data (Silva, 2009;Santos et al., 2012). *Adoption of air-entry of −2 cm. K s , saturated hydraulic conductivity; θ r , residual soil moisture; θ s , saturated soil moisture; α, air-entry suction factor; n, pore size distribution parameter; l, pore tortuosity and connectivity lumped parameter curve can be a solution to this as demonstrated by Heidbüchel et al. (2016), but was not the goal of this study. Further, a scatter plot between CRNS and PR2 datasets revealed a decent correlation of 0.80 and a tolerable error and underestimation of 14 and −9%, respectively. Moreover, a nearly perfect fit of the linear regression was also found between CRNS and PR2 datasets (slope = 1.13 and y intersect = −0.011 cm 3 cm −3 ), even though only PR2 data in <10-cm depth were available to derive the weighted soil moisture. Using a more complete dataset of the soil moisture profile probes, the vertical sensitivity as inherent in CRNS-derived soil moisture could be reduced (Scheiffele et al., 2020). Still, the results in this study show that CRNS can consistently detect soil moisture when operating in humid tropic areas and highlight its potential as a noninvasive soil moisture monitoring method to derive area average soil moisture values.

Soil hydraulic parameters
The van Genuchten-Mualem parameters were inversely calibrated from the CRNS neutron count rate data at S1 and based on the TDR-soil moisture data at S2 ( Figure 5). Overall, the inverse calibrations provided similar saturated and residual soil moisture values but different shape parameters, showing relatively distinctive retention and hydraulic conductivity curves. The CRNS-based retention was greater than the values derived from the TDR, mainly because of its lower parameter n equal to ∼1.26, which reduced the soil hydraulic con-ductivity, as shown by the lower K s value of 1.26 × 10 −4 m s −1 found S1 in comparison with 6.47 × 10 −4 m s −1 found in S2. The K s values were expected to contain higher uncertainties as the soil saturation during the study site rarely reached saturation (Figure 4). Sales et al. (2014) however, found a mean value of K s equal to 1.20 × 10 −4 m s −1 in forest areas of a representative watershed that encompasses GEB, which is very similar to the CRNS result. Silva (2009) and Santos et al. (2012) conducted a 200-m-mesh experimental study in GEB to obtain soil hydraulic parameters. Their results built the baseline for the water retention and hydraulic conductivity curves and their interquartile range (based on values of |h|) as used in this study to compare with the CRNS results. The CRNS-based curves were closer to the interquartile range than the TDR-based curves. Therefore, it was presumed that the soil hydraulic parameters calibrated from the CRNS data were a better match to the area average hydrologic properties. The difference in hydraulic conductivity curves between the S1 and S2 were, to some extent, the consequence of the different nature of the datasets used to optimize the soil hydraulic parameters (field-scale with CRNS and point-scale with TDR at S1 and S2, respectively). The adoption of the air-entry value of −2 cm for the CRNS-based simulations affects the shape of the curve in the range near saturation (Section 2.3). As the soil moisture rarely approaches saturation (Figure 4) this was expected to be of minor importance on the simulated fluxes from the soil column. On the other hand, the differences between the curves could potentially also be caused by slightly different soil properties at the two locations. To rule Vadose Zone Journal F I G U R E 6 Saturated hydraulic conductivity estimated by the slug test. Solid line = mean; dashed line = mean ± standard deviation out the effect of the location, further validations within this study were thus drawn only between datasets measured at the same site.
Slug tests were carried out in the well at S1 to verify the calibrated values of K s (Figure 6). The resulting K s values estimated by the six slug-in and slug-out attempts ranged from 1.18 × 10 −4 to 6.96 × 10 −4 m s −1 , which encompasses the results obtained by inverse calibration from the CRNS as well as from the TDR data. The slug test results thus confirmed the consistency of the inversely calibrated K s values, matching the soil hydrological characteristics of high permeability in medium clean sand (Cooper et al., 2017;Ishida et al., 2014).

Simulated neutron count rate and soil moisture
Aboveground neutron count rate simulated within HYDRUS 1D by the COSMIC operator was plotted against the measured CRNS neutron time series to demonstrate the good fit following the calibration of the soil hydraulic parameters (Figure 7). Overall, COSMIC-based time series adequately responded to the rainfall events, being able to consistently simulate the neutron count rate variability throughout the available data period. This simulated time series tended to slightly underestimate the CRNS neutron count rate particularly during the dry season, likely influenced by the standard parametrization of the calibration curve. The scatter plot between COSMIC and CRNS datasets highlighted the excellent fit of the linear regression (slope = 0.90 and y inter-sect = 44 cph), as well as a good performance expressed by a high correlation of .89 and an irrelevant mean error of 4% (16 cph). These results confirmed the good performance of the COSMIC operator as found by other studiers, because of its ability to deal with the decreasing sensitivity of CRNS with depth (Baatz et al., 2017;Brunetti et al., 2019;Shuttleworth et al., 2013). The calibrated soil hydraulic parameters enabled the model to simulate the integrated signal from CRNS very well.
To validate the CRNS-based HYDRUS-1D simulation at S1, the simulated soil moisture time series at the depths of 2, 20, 30, 40, 60, and 100 cm is compared to independent soil moisture measurements of PR2 (Figure 8). The results showed that the performance decreases as depth increases, except at the depth of 100 cm, but still providing acceptable levels of correlation varying from .70 to .90 and satisfying linear regression values (slope = 0.58 to 1.23 and y intersect = 0.000 to 0.061 cm 3 cm −3 ). Conversely, the accuracy was particularly low at −40 cm, showing overestimation of 36% and mean error of 44%, where the greatest leaching through the soil was observed. The overall decreasing correlation of CRNS to soil moisture with increasing depth is in alignment with observations from previous studies (Brunetti et al., 2019;Rosolem et al., 2014) and could be caused by several factors. Most obvious reason is the higher sensitivity of CRNS recordings to the shallower soil layers (Franz, Zreda, Ferre, & Rosolem, 2013;Lv et al., 2014). When comparing CRNS-derived soil moisture to profile measurements, usually stronger dynamics are observed in CRNS soil moisture (compare also Figure 4). Recent studies suggested the application of a profile shape correction, to reduce the inherent vertical sensitivity of the CRNS (Scheiffele et al., 2020), and thus possibly the discrepancy to other measurements with increasing depth, but is still subject of ongoing research. The observed deviation in correlation might also be in parts the consequence of measurements uncertainties at this depth. Finally, since only one profile probe was used (section 2.1), the deviations could also be the result of differences in investigation volume, i.e. the footprint.
After optimizing the soil hydraulic properties based on CRNS data, the integrated signal was reproduced very well ( Figure 7) but also the simulated time series show a more immediate response to the rainfall events ( Figure 8). Overall, the results showed that the CRNS based simulations were able to consistently simulate the soil moisture responses to rainfall events throughout most of the soil profile.

Land surface water fluxes and groundwater recharge
Simulated soil water fluxes at S1 and S2 were analyzed by the cumulative water fluxes of the RWU, evaporation, and F I G U R E 7 Comparison between the (a) neutron count rate simulated by the COsmic-ray Soil Moisture Interaction Code (COSMIC) and the cosmic-ray neutron sensing (CRNS)-corrected neutron count rate, with (b) correlation of between the COSMIC simulated and CRNS corrected count rates. DOY, day of year; CC, correlation coefficient; MAPE, mean absolute percentage error; PBIAS, percentage bias net surface fluxes (Figure 9). Ratios of actual to potential annual fluxes revealed pertinent spatial changes over GEB, since such ratios of RWU and evaporation remained on ∼64 and 63% at S1 and were ∼54 and 44% at S2, respectively. These differences were likely caused by the wetter soil condition at S1, as result from its shallower groundwater table since the capillary fringe rise exert a stronger influence on plant transpiration (Han et al., 2015). The actual RWU values were slightly greater than actual evaporation values, likely due to the high net surface fluxes (infiltration minus evaporation), whose actual-to-potential ratio remained on ∼89% at S1, and ∼84% at S2. These net surface fluxes were similar between S1 and S2 despite a higher hydraulic conductivity found at S2 (see Figure 5). Annual actual evapotranspiration rates from the CRNS-based simulation at S1 are 1,117 and 1,039 mm for 2017-2018 and 2018-2019, respectively. This is close to the mean annual evapotranspiration rate for Atlantic Forest that were estimated using distributed data (MOD16, MODIS Global Evapotranspiration Project) ranging from 1,161 to 1,389 mm (Barbosa, 2020;Cristiano et al., 2015;de Oliveira et al., 2016). The TDR-based simulations at S2 resulted in much lower and less realistic values of 877 and 884 mm.
The simulated effects of the rain events on percolation and groundwater fluxes at S1 and S2 are shown in Figure 10. Due to the shallow groundwater table at GEB, the recharge can already occur at depths ≥100 cm. The vertical percolation fluxes were reduced by the RWU throughout the soil profile. The groundwater recharge varied temporally due to combination of rainfall, RWU, and evapotranspiration, resulting in annual rates ranging from 60 to 470 mm yr −1 at S1, corresponding to 5 and 29% of rainfall, respectively. These results were lower than the recharge rates at S2, which ranged from 230 to 620 mm yr −1 , corresponding to 18 and 38% of rainfall, respectively. Given the results of the previous sections, the most important factors distinguishing the recharge rates at S1 (upstream area) and S2 (downstream area) are the soil hydraulic properties and groundwater table depths.
The cumulative groundwater recharge flux was compared with the WTF estimates at S1 and S2 for S y equal to 0.10, 0.16, and 0.24 ( Figure 11). The specific yield of 0.16 found at S1 was similar to the mean estimate of ∼0.15 found in a management zone of water resources close to GEB (Braga et al., 2015), which acted as a reference value, whereas 0.10 and 0.24 are an estimated minimum and maximum limit, respectively. The simulations at S1 showed good performance in both hydrological years for S y equal to 0.16, underlined by KGE equal to 0.83 and 0.79. However, using the same S y value, estimations at S2 in the rainy hydrological year resulted in a lower accuracy, indicated by a KGE equal to just 0.36. On the other hand, the KGE was 0.84 when using S y equal to 0.24. This behavior can be explained in part by the sensitivity of S y values according to the depth to the groundwater table (Dettmann & Bechtold, 2016), as well as by the TDRbased point-scale calibration of soil hydraulic parameters used for the GWR estimation. The better match between the simulated recharge at S1 to the WTF estimates can be explained by the integrative nature and wider spatial range of influencing moisture as association of both data products, the neutron count rate being representative of the field-scale soil moisture and the groundwater table fluctuations also representing the area average behavior. Although the study period comprises two contrasting hydrological years, it needs to be emphasized that limitations to the implications may arise from the short timeframe. Strong seasonality and extreme hydrological conditions, which could be expected in this climate, were not sufficiently represented by this study (Wendland et al., 2007). Nevertheless, the basic conclusion is that using CRNS data compared with point-scale soil moisture data to inversely estimate soil hydraulic parameters to derive dynamic groundwater recharge rates at field scale is at least a suitable, probably more representative alternative. F I G U R E 8 Soil moisture simulation at S1, based on soil hydraulic properties calibrated with cosmic-ray neutron sensing (CRNS) data plotted for validation against the PR2 data in the depths of (a, b) 2 cm, (c, d) 20 cm, (e, f) 30 cm, (g, h) 40 cm, (i, j) 60 cm, and (k, l) 100 cm (only one probe delivered consistent data for the full measurement period). CC, correlation coefficient; MAPE, mean absolute percentage error; PBIAS, percentage bias F I G U R E 9 Cumulative water flux of evaporation, transpiration, root water uptake, infiltration, and net surface water flux at (a) S1 and (b) S2. Showing similar ratios of potential to actual net surface fluxes (infiltration minus evaporation) for both sites but distinct evapotranspiration rates. Especially the ratio of potential to actual evaporation is low at S2 (∼44%) compared with S1 (∼63%). CRNS, cosmic-ray neutron sensing; TDR, time-domain reflectometry; DOY, day of year F I G U R E 1 0 Cumulative vertical percolation flux at (a) S1 and (b) S2, in the depths of 20, 30, 40, 60, and 100 cm, and the variable bottom boundary for the dynamic groundwater recharge rate estimation. CRNS, cosmic-ray neutron sensing; TDR, time-domain reflectometry; DOY, day of year

CONCLUSION
The study assessed the capacity of CRNS for (a) monitoring soil moisture at the intermediate scale under the unfavorable, humid tropic conditions, and in (b) the use of CRNS to inversely calibrate soil hydraulic parameters and simulate dynamic groundwater recharge rates. The study used CRNS data from two contrasting hydrological years in an experimental basin and supported by a range of ground truth information. The results showed that the CRNS-based soil moisture time series responded consistently in alignment with the observed soil moisture profile measurements, despite its relatively scattered dataset resulting from the low incoming neutron count rate. Using the COSMIC operator coupled with HYDRUS-1D model enabled to inversely calibrate the soil hydraulic parameters based on CRNS data, and, in turn, to satisfactorily simulate the aboveground neutron count rate and soil moisture time series. The resulting soil hydraulic parameters were found to be in line with ranges reported in literature, whereas the saturated hydraulic conductivity was successfully validated F I G U R E 1 1 Comparison of the cumulative groundwater recharge flux estimated by the one-dimensional soil hydrological model HYDRUS-1D and water table fluctuation method for specific yield (S y ) equal to 0.10, 0.16, and 0.24, at (a, b) S1 and (c, d) S2. The negative sign for the flux indicates the downward direction. CRNS, cosmic-ray neutron sensing; TDR, time-domain reflectometry; KGE, Kling-Gupta efficiency; R WTF , groundwater recharge estimated using the water table fluctuation method by local experiments. Moreover, the HYDRUS-1D model enabled us to simulate the vertical water percolation fluxes of the root zone and the dynamic groundwater recharge rates at the variable bottom boundary. Compared with WTF-based groundwater recharge estimates, the CRNS-based simulated recharge performed better than simulated recharge rates based on point-scale data (used to optimize the soil hydraulic properties). Given the results presented here, we concluded that CRNS derived hydrological parameters could provide a more effective support to derive field-scale groundwater recharge rates, as it represents integral information over an intermediate scale, even under tropical conditions. Recent developments in CRNS detectors towards different designs and more sensitive devices (Stevanato et al., 2019;Weimar et al., 2020) could help to address unfavorable conditions with respect to low neutron count rates and might provide even more precise results in the future. This encourages the adoption of such technology in developing tropical countries, whose forest protection and aquifer groundwater level conservation are essential for preserving the terrestrial and estuarine biodiversity, as well as whose economy strongly depends on agricultural practices. Accordingly, it is expected that these findings will increasingly promote the CRNS hydrological applications to provide appropriate estimates of groundwater recharge rates, important for diligent and sustainable management of groundwater resources by decision-makers in regions with similar characteristics. However, further studies exploring the hydrological poten-tials of the CRNS-based simulations, related to dynamic groundwater recharge rates in other tropical wet and semiarid regions, are strongly recommended to expand the knowledge of water resources and to foster the good sustainable practices.

A C K N O W L E D G M E N T S
The authors would like to acknowledge the financial support granted by (a)  Baroni to this study was partly conceived during the IAEA Coordinated Research Project (CRP) "Enhancing agricultural resilience and water security using Cosmic-Ray Neutron Sensor (D1.20.14)." The author Suzana Montenegro appreciates the PQ scholarship granted by CNPq. The authors are finally grateful for the thoughtful inputs of Professor Carlos de