Performance of three reanalyses in simulating the water table elevation in different shallow unconfined aquifers in Central Italy

Water table elevation is a key feature for identifying the groundwater behaviour. Accordingly, appropriate measurements—in terms of both frequency and spatial distribution—play a crucial role for capturing the aquifer response to recharge and withdrawals. However, numerical models simulating the main features of the behaviour of the water table elevation may help groundwater management, as an additional tool. In this article, soil moisture data from three well‐established global reanalyses (ERA5, CFS, and JRA‐55) are used for evaluating the flux in the vadose zone towards shallow unconfined aquifers, Fg$$ {F}_g $$ , in the Umbria region (central Italy). Then, according to the methodology proposed in Bongioannini Cerlini et al. (2021), where for the considered aquifers most of the recharge derives from the unsaturated zone, Fg$$ {F}_g $$ is used for simulating the water table evolution in time. With the aim of assessing which reanalysis is the most appropriate in simulating the evolution of groundwater levels, the properties of the correspondent land surface models (LSM) are examined, as they provide Fg$$ {F}_g $$ . For the considered aquifers, the analysis of the performance of the selected reanalyses confirms the validity of the proposed approach. Moreover, it points out the crucial role of the spreading of the water table elevation with respect to its mean value, as a significant parameter for selecting the most adequate reanalysis to use. In addition, the role in the LSM of the explored soil depth, hydraulic conductivity curve, and spatial resolution is highlighted. These results, in line with recent literature on the performance of the reanalyses, suggest to extend future work to other regions of the world.


| INTRODUCTION
In most countries, the socio-economic development and population growth cause an increase of water consumption that may lead to exploiting indiscriminately groundwater resources. A possible effect of the groundwater mismanagement is the depletion of the aquifer storage and then the decline of its water table elevation (Banerjee et al., 2009). The drop in the groundwater levels has remarkable effects for both the environment and human activities. In fact, it reflects in the decrease of the discharge from groundwater to streams, with potentially severe effects on aquatic ecosystems (de Graaf et al., 2019), and then it may put at risk ecological maintainment, as an example, of surface water (e.g., Caschetto et al., 2014). Moreover, since groundwater is one of the largest freshwater resources, its depletion may restrict irrigation and drinkable water supply (Giordano, 2009). As a consequence, information about the temporal and spatial variability of the aquifer recharge is of noticeable importance for water managers (Li et al., 2021). Groundwater recharge can be evaluated by means of several techniques (e.g., the water budget model, the waterbalance method) (Scanlon et al., 2002). However, in the water-balance equation, precipitation and runoff are relatively easy to measure, whereas recharge remains an elusive process to quantify (Lee et al., 2006). As a matter of fact, it depends not only on several parameters (e.g., soil type, vegetation cover, slope) but it also derives from understanding the hydrological processes that are the result of the water-heat balance between the land and atmosphere that are in a continuous, dynamic interaction (Assefa & Woodbury, 2013). In particular, shallow unconfined and semi-unconfined aquifers are in direct interaction with the atmosphere and their behaviour is strongly linked to both precipitation and evapotranspiration (Li et al., 2021). For example, a study by Shrestha (2021) analyses the effect of the spatial distribution of clouds on evapotranspiration. Hence, the soil moisture state (i.e., the amount of water stored in the unsaturated zone) is the final result of the mentioned interaction phenomena, reflecting changes in the meteorological conditions. Particularly, in shallow unconfined aquifers, on which attention is focused on this article, soil moisture may contribute significantly to the aquifer recharge.
The use of general circulation models (GCMs), which include land surface models (LSMs), is the most appropriate approach for evaluating soil moisture as a result of joint investigation of soil moisture-climate interactions by various disciplines (Seneviratne et al., 2010). However, such an approach requires suitable field observations for model calibration, with those concerning the soil moisture state being very demanding. Such an integrated approach is followed in Assefa and Woodbury (2013) and Zhang et al. (1999) where climate data, in situ observations of soil moisture (measured at a depth ranging from 10 cm to 1.0 m), water table elevation, and temperature are combined with a Richards equation-based LSM to provide the groundwater recharge. Beyond the evident large cost of field measurements and model tuning, a limit of such an approach is that the procedures are strongly linked to the available datasets and used models. This makes quite arduous transferring the refined procedure to other catchments. Very recently, a suite of recharge estimates is provided in Li et al. (2021), where the LSM simulates evapotranspiration, soil moisture state, and runoff. Then these quantities are used for evaluating the monthly and annual recharge.
The approach based on the use of GCMs and related LSMs requires advanced multidisciplinary skills. This makes it substantially impractical for most of water supply companies. As a consequence, several approximate methods have been used in the literature for evaluating the groundwater recharge. In this context, two main families of methods can be identified: the data-driven models and the physically based ones.
A review of the numerous data-driven models for simulating the water table elevation is offered in Rajaee et al. (2019); the one using the artificial neural network (Banerjee et al., 2009) and discrete space-state model  are cited here as noticeable examples. Such an approach is attractive because of its ability to capture non-linear trends without prior knowledge on the underlying physical processes due to data limitations. Moreover, at present, no particular type of data-driven model can be recommended for a given problem (Rajaee et al., 2019).
Physically based models correlate water table fluctuations to groundwater recharge, bypassing the evaluation of the soil moisture state. This approach produces recharge estimates that compare well other methods when the water table elevation varies smoothly and aquifer properties are relatively well known (Cuthbert, 2010). The water table fluctuation method requires the specific yield as the only parameter (Gumuła-Kawęcka et al., 2022). However, such a quantity is difficult to determine, as it depends on the type of soil, depth to the groundwater table, and water dynamics in the vadose zone (Healy & Scanlon, 2010).
Within an alternative approach, Cerlini et al. (2017) and Bongioannini Cerlini et al. (2021) highlight that, for shallow unconfined aquifers in Umbria region, a good correlation exists between water table observations and fluxes from the vadose zone, as given from ERA5, the global atmospheric dataset (reanalysis) provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). Such a result allows assuming that, for the considered aquifers, the water flux from the vadose zone is a relevant portion of the aquifer recharge. Later on, in Li et al. (2021), it is shown that the link between the maximum of the measured seasonal soil moisture and groundwater storage stands also for aquifers located in the eastern United States and Pacific Northwest.
This article aims to investigate how the simulated flux towards groundwater from the vadose zone and its correlation with water table observations depends on the LSM features. For this purpose, three well-established global reanalyses with a different LSM have been compared: ERA5 reanalysis, the Climate Forecast System Reanalysis (CFSR/CFSv2), provided by the National Centers for Environmental Prediction (NCEP), and the JRA-55 reanalysis, provided by the Japan Meteorological Agency (JMA). The analysis focuses on the main LSM features that are directly linked to the evaluation of the flux in the vadose zone: soil moisture state, hydraulic conductivity curve, and soil column discretization. The performance of the selected reanalyses is evaluated by correlating such a flux with the observations at 10 piezometers from different shallow unconfined aquifers in the Umbria region (central Italy). Accordingly, this article is structured as it follows. The characteristics of the piezometer network in the Umbria region are reported in Section 2. The main features of the LSMs used in the considered reanalyses are discussed in Section 3. The methodology used for correlating the flux towards the aquifer and water table elevation is described in Section 4. The obtained results with a comparison between the reanalysis products are illustrated in Section 5. A discussion and the conclusions are finally provided.

| THE STUDY AREA AND THE PIEZOMETER NETWORK
The measurements of the water table elevation, which play an irreplaceable role in the proposed approach, are provided by the regional piezometric network of the Umbria region (central Italy), managed by the Regional Environmental Protection Agency [Agenzia Regionale per la Protezione Ambientale (ARPA)]. Daily water table elevation measurements, h d w (in m a.s.l.) have been collected from early 2000s to present (the subscript d indicates daily quantities). Daily data have been subjected to preliminary quality control procedures, before being considered in the following analysis (see Bongioannini Cerlini et al., 2021, for a detailed description of the quality controls).
Within the available piezometers, 10 have been selected according to the following main criteria (Bongioannini Cerlini et al., 2021): 1. the aquifer is alluvial and unconfined, for a direct interaction between land and atmosphere to happen; 2. no pumping is active around the piezometer, to avoid artificial changes in the water table elevation; 3. the mean depth of the water table, D w (in m from the ground level), ranges between 4 and 10 m, for the flux in the vadose zone to be a possible, significant recharge mechanism (Seibert et al., 2003).
The selected piezometers ( Figure 1) monitor five different alluvial and unconfined aquifers; their main characteristics are reported in Table 1 for the period 2001-2012, the one used for the model calibration (the remaining period 2013-2018 is considered for evaluating the model performance, see Section 4). In Table 1 σ w , the mean absolute deviation of h d w is a measure of its variability in time.
According to the value of σ w , two groups of piezometers can be identified: (i) the piezometers P4, P5, P6, and P10, with a highly variable in time water table (σ w > 1.5 m), and (ii) the piezometers P1, P2, P3, P7, P8, and P9, with a more stable behaviour (σ w ≤ 1.5 m). For these groups of piezometers, the behaviour in time of the fluctuation of the water table elevation around the mean value, is the mean daily water table elevation), is shown in Figure 2a, whereas their power spectrum, PS, is depicted in Figure 2b. These plots clarify the different dynamical behaviour of the two mentioned groups of piezometers: high-variable water tables have a smoother behaviour with more energy concentrated in larger scales and a larger seasonal variability, whereas low-variable water tables have more energy concentrated in smaller fluctuations and a smaller seasonal variability.

| SOIL MOISTURE BUDGET MODELLING IN LAND SURFACE MODELS (LSMS)
As mentioned, in this performance analysis, the following three reanalyses have been selected: (i) the Climate Forecast System Reanalysis (CFSR), for the period 1979-2011, and its upgraded version 2 (CFSv2), from 2011 to present, provided by the National Centers for Environmental Prediction (NCEP)-hereafter referred to simply as CFS; (ii) the ERA5 reanalysis provided by the European Centre for Medium-Range Weather Forecasts (ECMWF); and (iii) the JRA-55 reanalysis provided by the Japan Meteorological Agency (JMA). A brief description of the performance of the "third" generation reanalysis is offered in Supplementary Material S1.
A reanalysis system is made by different components: an atmospheric model, an LSM, an ocean model, a sea-ice model. The reanalysis performance depends not only on how each model simulates the physical processes, but also on the assimilated observations and the way the system components interact with each other. A complete description of the components of the reanalyses-which is beyond the scope of this article-can be found in Saha et al. (2010Saha et al. ( , 2014 for CFS, Hersbach et al. (2020) for ERA5, and Kobayashi et al. (2015) for JRA-55. A comprehensive analysis of the philosophy, perspectives, and numerical tools implemented in the LSMs is offered in Fisher and Koven (2020). Instead, this article focuses on a specific aspect of the LSM: the soil moisture budget and the computation of the water flux in the unsaturated zone.   LSMs are responsible for the parametrization of the land surface processes in the GCMs: CFS uses the NOAH LSM (Ek et al., 2003), ERA5 the H-TESSEL LSM (Balsamo et al., 2009), and JRA-55 the SiB LSM (Sellers et al., 1986). ERA5 has also a specific component for the analysis of the behaviour of water bodies (e.g. FLake, Bongioannini Cerlini et al. (2022)).
In the LSMs, two main tasks can be identified: the land surface energy balance and surface hydrology. While the former provides fluxes at the surface and couples them with the overlying atmosphere model, the latter computes the water and heat budget in the investigated soil column. A detailed description of the land surface energy balance is given in Chen and Dudhia (2001); Ek et al. (2003); Pan and Mahrt (1987) for NOAH; ECMWF (2016) and Balsamo et al. (2009) for H-TESSEL, and Sellers et al. (1986) for SiB.
The boundary condition at the top of the soil column is given by the surface energy budget that is different for each LSM. On the contrary, at the bottom all LSMs assume free drainage, with the flux towards groundwater, F g , given by: where γ bot = hydraulic conductivity of the bottom (i.e., the deepest) soil layer. In ERA5, the soil hydraulic properties, λ and γ, are evaluated by means of the Van Genuchten (1980) analytical formulation: where p is the pressure head, which is linked to the soil moisture using the following relationship: where θ r = residual soil moisture, and parameters α and n depend on the soil texture type. CFS and JRA-55 adopt the Clapp and Hornberger (1978) where parameters θ s , γ s , and b depend on the soil texture type, as specified in Cosby et al. (1984).
In the considered reanalyses, the soil texture type does not change in time and is constant through the investigated soil column, but it is different for each grid point. ERA5 adopts seven types of soil texture (i.e., sea, coarse, medium, medium fine, fine, very fine, and organic), and the related parameters can be found in Table 1 of Balsamo et al. (2009). CFS adopts 19 different types of soil texture (e.g., sandy, loamy, silty, clayey) and the related parameters can be found in Table 2 of Chen and Dudhia (2001). In JRA-55, the soil texture type depends on the vegetation type. This approach differs from the one followed in CFS and ERA5 where the soil texture and vegetation type are two separated invariant parameters. JRA-55 soil parameter tables can be found in Japan Meteorological Agency (2016) and Table 3c of Hirai et al. (2007), where only θ s depends on the vegetation type, whereas the other parameters are kept constant. For the considered reanalyses, only one soil texture type is present in the study domain, with the related parameters reported in Table 2.
The main characteristics of the LSMs used in the considered reanalyses are reported in Table 3. A more indepth discussion of both temporal and spatial resolution and soil column discretization is given in Supplementary Material S2.

| THE MODEL FOR SIMULATING THE WATER TABLE ELEVATION
In the used approach, the flux from the deepest soil layer, integrated on a daily and monthly scale, provides the daily, F d g , and monthly, F m g , fluxes towards groundwater from the vadose zone (the subscript m indicates monthly quantities). To analyse the correlation between such integral fluxes and the observed water table elevations, a cross-correlation function based on the Spearman rank correlation coefficient is used. Such a parameter allows evaluating the time lag, τ, between the computed flux and observed change in the water table elevation. Then, the water table elevation is obtained from the fluxes by T A B L E 2 Soil parameters of the Clapp and Hornberger (1978) and Van Genuchten (1980) formulations for the soil texture type present in the study domain and considered reanalyses.

Clapp and Hornberger
Van Genuchten where h max w is a reference maximum value of the daily water table elevation (calculated as the 99th percentile), and F max g is a reference maximum value of the daily integral flux (evaluated as the maximum of the reanalysis fluxes). Equation (7) can be reshaped in a more compact non-dimensional form: =σ w is the standardized relative water table elevation, and F dÃ g ¼ F max g =F d g is the relative flux towards the aquifer. Coefficients k d log and k d lin are obtained by means of the non-linear least square interpolation method.
Equations (7) and (8) refer to daily quantities, but they can be applied to monthly quantities in the below equivalent form: with the monthly relative fluxes, F mÃ g , and standardized relative elevation, h mÃ r , evaluated by considering the monthly mean of the daily fluxes and water table elevation, respectively.
The above dimensionless equations allow comparing piezometers with a different dynamics. In particular, the value of the two dimensionless daily (monthly) parameters k d log (k m log ) and k d lin (k m lin ) regulates the relative magnitude of the non-linear and linear components in the flux-water table elevation relationship. With respect to Bongioannini Cerlini et al. (2021), in this article, the calibration of the k log and k lin parameters has been done by using the observations from 2001 to 2012. This period has been chosen in order to ensure at least 6 years of observations for the time series starting in 2006 (Table 1). Different values of the k log and k lin parameters are expected for each location since they account for the local variability of the flux towards the aquifer as well as its dynamical properties. For checking the model performance, the Kling-Gupta efficiency parameter, KGE (Gupta et al., 2009): has been evaluated in the remaining 6-year observation period (i.e., from 2013 to 2018). In Equation (10), r is the linear coefficient of correlation, η is the ratio between the simulated and observed standard deviation, and μ is the ratio between the simulated and observed mean.

| PERFORMANCE OF THE LAND SURFACE MODELS
According to Equation (2) the flux in the vadose zone towards groundwater depends on the soil texture type and value of the soil moisture at the bottom layer; the latter can be extracted from the related time series at the model grid point nearest to the selected piezometer. Accordingly, in the below two subsections, the considered reanalyses are compared with regard to the soil moisture state and water flux towards the aquifer, both affecting the aquifer behaviour. Then, the correlation between F g and water table elevation is discussed. Finally, the performance of the used model for simulating the water table elevation is evaluated.

| Comparison of the soil moisture state
In Figure 3, the soil moisture time series at the available depths for piezometers P1 and P8 are reported, as representative examples, for a qualitative comparison. In this figure, the curves at z = 1.95 m for ERA5, z = 1.07 m for JRA-55, and z = 1.5 m for CFS, represent the series of the soil moisture at the deepest layer.
The very different dynamics and average values of such time series are the inevitable result of the differences between the used LSMs. In particular, the different soil parameters and hydraulic conductivity equation, as well as the different depths of the soil layers, make comparisons quite arduous. To get out of this impasse without deviating from the goal of this article, the selected reanalyses are compared by considering the soil moisture integral, SMI: expressed in m 3 m À2 , where z c is the smallest soil column depth considered in the LSMs. This is the case of JRA-55 for which it is z c = 1.57 m. In other words, according to Equation (11), SMI takes into account the soil depth in common between ERA5, CFS, and JRA-55. According to Bongioannini Cerlini et al. (2021), the value of SMI does not depend significantly on the assumed soil moisture profile within each layer (i.e., linear or constant). As a representative example, Figure 4a shows a qualitative comparison of the SMI values for the grid point close to piezometer P8 from 2009 to 2011. The seasonal cycle of the soil moisture is clear in all reanalyses, with ERA5 and CFS oscillating with a larger amplitude than JRA-55. Such a difference is confirmed by the monthly statistics summarized in Figure 5 where the JRA-55 median values exhibit a smaller variability over months. On average, ERA5 has a larger SMI value over all months when compared with CFS and JRA-55, which show similar average values of SMI, except for the period between June and September, where CFS exhibits the smallest SMI value among all reanalyses. The SMI variability over space (i.e., at the selected grid points) and time (2001present) depends also on the month: ERA5 and CFS show a larger variability from October to June and a smaller variability from July to September. In this period, there is a large difference in the SMI distribution between ERA5 and CFS, with the latter showing very dry soils with very small F I G U R E 5 Monthly box plots of the soil moisture integral, SMI, for the considered reanalyses and all piezometers.
F I G U R E 6 Hydraulic conductivities curves of the considered reanalyses, as obtained from parameters in Table 2. variability. Such a feature is also evident in the behaviour of the soil moisture in time and with depth, as shown in Figure 4b-d for ERA5, CFS, and JRA-55, respectively. During the driest period (i.e., from June to September), CFS shows the smallest values of θ and the fastest propagation of the dryness from the surface to the bottom soil layer. As shown in Figure 6, in this range of θ values, the hydraulic conductivity is very small and the curves show a very large gradient, meaning that a small variation of the soil moisture implies a large variation of the conductivity. As a consequence, the smaller soil conductivities of CFS in the first three soil layers cause a persistent dry period with a very small variability of SMI extending from July to September. This feature is significant in the analysis of the flux towards groundwater (Section 5.2), especially from July to September. Even if ERA5 shows the wettest soil in all periods, it exhibits a slower propagation of the soil moisture from the surface to the bottom layer. This is probably due to the large difference in the formulation of the hydraulic conductivity with respect to CFS and JRA-55, as discussed in the next subsection.

| Comparison of the flux in the vadose zone towards groundwater
According to Equation (2), the flux in the vadose zone towards groundwater depends on the bottom soil moisture, θ bot , and hydraulic conductivity curves ( Figure 6). The values of θ bot are very different among the considered reanalyses as an unavoidable result of the used LSM and soil depth, z bot . As shown in Figure 7a,b, the bottom layer of JRA-55 exhibits large variations of the soil moisture, especially in summer, where it frequently becomes very dry, reaching values smaller than 0.25 m 3 m À3 . The same oscillations are present in ERA5, but with a much larger mean soil moisture, which rarely is smaller than 0.25 m 3 m À3 . CFS shows the smallest variation of θ bot , especially during the period between June and October. Such a behaviour is strictly connected to the SMI dynamics highlighted in the previous subsection. Figure 7c qualitatively points out the difference between the monthly flux towards groundwater given by the considered reanalyses. The first feature to notice is the discrepancy between the magnitude of CFS and JRA-55 fluxes with respect to the one of ERA5, which is almost two orders of magnitude smaller. The monthly statistics in Figure 7d highlights the different range of values between ERA5 and the other reanalyses in the whole recharge period (i.e., from December to May). Such a large difference can be ascribed to the different hydraulic conductivity curve used in the respective LSM. According to the curves reported in Figure 6, all reanalyses show a great sensitivity to small variation of the soil moisture, especially in relatively dry soil (Chen & Dudhia, 2001): the larger the soil moisture, the smaller the gradient of hydraulic conductivity. The significant differences between ERA5, CFS, and JRA-55 in terms of the hydraulic conductivity curve are highlighted in Figure 6 by considering the following three soil moisture ranges: 1. for θ ≤ 0:125, ERA5 and CFS have similar conductivites at very small values of θ, but for larger values of θ they start to diverge. Precisely, for a given θ, JRA-55 provides smaller values of γ with respect to CFS, but the same gradient. Similar values of the conductivities for JRA-55 and ERA5 can be found only for θ ≈ 0:125; 2. for 0:125 < θ ≤ 0:225, the JRA-55 conductivity curve has a larger gradient than ERA5 and CFS. In JRA-55 and CFS, similar values of γ happen for θ ≈ 0:225. ERA5 has a gradient close to the one of CFS, but one order of magnitude smaller conductivites; 3. for θ > 0:225 up to the saturation value θ s (Table 2), ERA5 maintains a positive gradient very similar to the one of CFS, but almost two order of magnitude smaller conductivities. In this range, JRA-55 exhibits conductivities larger than CFS. In ERA5 the Van Genuchten formulation implies a sharp gradient close to the saturation value. This implies that for values θ ≈ θ s , ERA5 has a large increment of the flux for a small increment of the soil moisture compared with the other reanalyses.
Accordingly, the large differences in monthly fluxes between ERA5 and other reanalysis (Figure 7d) are due to the hydraulic conductivity curve of ERA5 that, for values of θ larger than 0.2 m 3 m À3 , is always two orders of magnitude smaller than the one of CFS and JRA-55, for a given soil moisture content. In other words, the difference induced by the hydraulic conductivity curve prevails on the one due to the value of θ bot . The situation is different in the period between June and November when the very small values of θ bot of JRA-55 ( Figure 7a) cause a strong decrease of the hydraulic conductivity whose values are comparable with those of ERA5. This is not true for CFS that shows fluxes larger than those given by JRA-55 and ERA5. Such a feature could be ascribed to the effect of evapotranspiration that is more significant in JRA-55 than in ERA5 and CFS, since the JRA-55 bottom soil depth is smaller.

| Correlation between flux in the vadose zone and water table elevation
As mentioned, the proposed model assumes a correlation between the fluxes calculated from the reanalysis data and measured water table elevations. Accordingly, the correlation coefficient has been evaluated for different time lags, τ, in order to investigate the time shift between the flux in the vadose zone and water table elevation. In Table 4, where the piezometers are ordered by increasing σ w , the maximum value of the Spearman rank correlation coefficient, R max , is reported along with the time lag, τ max , that maximizes it. As expected, τ max is proportional to the soil column depth investigated in the reanalysis: the deeper the soil column the smaller the time lag between the flux and water table elevation. In a short soil column, as that of JRA-55, the maximum value of the correlation can be reached even after 40 days. For all the selected reanalyses, the larger (smaller) values of τ max concern piezometers with larger (smaller) σ w . The only exception to such a behaviour is piezometer P9 for CFS and JRA-55 reanalysis. Table 4 shows also that, for a given reanalysis, the largest correlation coefficients happen at the same F I G U R E 8 Scatter plot of daily standardized relative elevation, h dÃ r , against daily relative fluxes, F dÃ g , from 2001 to 2012. Red curves represent logarithmic-linear relationship of Equation (8), with parameters reported in Table 5. Piezometers are in order, from top to bottom, by increasing σ w and, from left to right, by reanalysis (ERA5, CFS, and JRA-55). locations (i.e., P4 and P5). In these locations, all reanalyses show a time lag proportional to the distance between the bottom soil layer and mean water table elevation. On the whole, the correlation of the ERA5 and CFS fluxes with the water table elevation is generally larger than for the JRA-55 ones.
The correlation between the daily fluxes and water table observations can be further analysed by considering the scatter plots reported in Figure 8. The analysed data cover the time period from 2001 to 2012, the same used later to calibrate the parameters considered for simulating the water table elevation. As discussed in Section 4, dimensionless flux and water table elevation are considered for all piezometers, sorted by increasing σ w . The range of relative fluxes of JRA-55 is two order of magnitude larger than ERA5 and CFS. The reason for such a difference can be ascribed to the large variability of the JRA-55 fluxes, as shown in Figure 7d. The ratio between the maximum and minimum daily recharge throughout the year can reach values of 10 5 , in contrast with the smaller values of ERA5 (10 3 ) and CFS (10 2 ). The relationship between the standardized relative elevation, h dÃ r , and the relative flux, F dÃ g , is very different among piezometers and reanalyses.
Only in few cases, such as P3 and P4, it is clear whether the relationship can be assumed purely linear (P3) or purely logarithmic (P4). This is the reason why a mixed linear-logarithmic model has been introduced in Bongioannini Cerlini et al. (2021). Red curves in Figure 8 show the relationship obtained by using such a model, that is, Equation (8). As discussed in Section 4, the model needs some a priori information describing the main features of each location: σ w and h max w , concerning the water table elevation and obtained from the measurements, and F max g , evaluated by data from reanalysis. For the considered locations, the values of such parameters are reported in Table 5.
The ratio between the two model regression coefficients, k lin =k log , as obtained from non-linear least square interpolation, is also reported in Table 5 for both daily, k d lin =k d log , and monthy, k m lin =k m log , values. Such a ratio indicates whether the linear component prevails on the logarithmic (positive values) or viceversa (negative values). Values of k lin =k log much larger than unity mean that the relationship between the relative fluxes, F dÃ g , and standardized relative water table elevation, h dÃ r , is almost linear. The magnitude and the sign of this ratio are different for the considered reanalyses. Concordant signs are mostly found for ERA5 and CFS, but their magnitude differs for most of the piezometers. On the contrary, JRA-55 shows very large positive values for all piezometers, meaning that a pure linear relation could be used for simulating the water table elevation on the basis of its data. The dependence of k lin =k log on the considered reanalysis is a clear indication that this ratio is not an intrinsic property of the considered aquifer. However, the monthly values of the ratio k lin =k log are very close to those of the daily quantities, indicating that the above considerations about the linear or logarithmic behaviour of the piezometers are still valid also for monthly quantities. Moreover, a link between the mean relative standardized elevation, h dÃ r : values of k lin =k log ) exhibit smaller oscillations (i.e., with small values of σ w and large values of h dÃ r ) than those characterized by a logarithmic behaviour. Such a feature could be due to different recharge mechanisms. Precisely, as discussed in Section 2, piezometers P4, P5, P6, and P10 show highly variable water table elevations that could indicate that a diffuse recharge happens, with a more direct connection between the recharge and atmospheric processes. Such piezometers show a predominant logarithmic behaviour. On the contrary, the other piezometers, with the only exception of P9, are characterized by a small variability and linear behaviour, which could depend on a focused recharge.

| Simulation of the water table elevation
The water table elevation is simulated by substituting the computed fluxes and model parameters (Table 5) into Equation (7) or in its dimensionless form Equation (8). As anticipated, the performance of the model is evaluated for a period of time, that is, the 6-year period from 2013 to 2018, different from the one considered for the parameter calibration (i.e., 2001-2012). The model has been tested on a daily and monthly scale, since both are relevant for groundwater management purposes.
The parameters of the dimensionless flux and water table elevation (h max w , σ w , and F max g ) are kept equal to those used for the daily quantities (Table 5), whereas different monthly linear and logarithmic coefficients (k m lin , k m log ) have been obtained by least square interpolation. The comparison between the observed and simulated water table elevations is shown in Figure 9 for all the considered reanalyses and piezometers. The resulting KGE coefficients for daily, KGE d , and monthly, KGE m , simulations are reported in Table 6.
The analysis of the curves in Figure 9 indicates that the simulations based on ERA5 generally perform better for smooth and largely oscillating water tables (see the KGE coefficients for piezometers P4, P5, P9, and P10 in Table 6). On the contrary, the CFS simulations perform similarly to ERA5 for small fluctuating water tables (see the KGE coefficients for piezometers P1, P2, P3, P8 in Table 6). The driest periods and water table elevation extremes remain the most critical periods to be simulated, especially when the CFS fluxes are used. However, such a worse performance in the driest periods can be ascribed to the fact that in such condition the relevance of the flux in the vadose zone, on which the proposed model is based, can be negligible with respect to the flow in the aquifer. The water table simulation concerning piezometers P8, P9, and P7 by CFS shows a very poor representation of the seasonal cycle. On the other hand, ERA5 cannot capture all the variability of the water table in P8 and P7. In terms of the KGE coefficient, this means that CFS reanalysis performs better than ERA5 over piezometers P7 and P8, because of the larger α coefficient in Equation (10).
ERA5 shows the largest KGE coefficient for piezometer P10, altough the simulation in Figure 9 shows a larger bias of ERA5 if compared with JRA-55 and CFS. However, the adopted formulation of KGE in Equation (10) assigns an equal weight to all the three KGE components (correlation, variability, and bias). Assigning a different weight to the bias component, μ, would have raised KGE coefficient from CFS and JRA-55 (Knoben et al., 2019).
JRA-55 performs better than the other reanalyses only for piezometer P7, while it shows very small KGE coefficients for piezometers P9. The reason for such a negative performance could be imputed to the smaller investigated soil depth and coarser spatial resolution of JRA-55.
On average, the simulations performed using ERA5 and CFS show the largest KGE efficiency, which is never smaller than the limit value (= 0.4), both on a daily and monthly scale. From a qualitative point of view, Figure 9 highlights a better agreement between the measured values and those based on the ERA5 and CFS reanalyses with respect to those given by JRA-55. Drought periods remain challenging for the model and all the considered reanalyses: larger errors characterize all piezometers with regard to the simulation of the extreme values of the water table elevation.

| DISCUSSION AND CONCLUSIONS
This article analyses how the flux in the vadose zone towards the aquifer, estimated from soil moisture data given by reanalyses, is influenced by the properties of the land surface models (LSMs) of each reanalysis. It also evaluates the connection of such a flux with the observed water table elevation. With this aim, three different global reanalyses have been compared: ERA5, provided by ECMWF, CFS, provided by the NCEP, and JRA-55, provided by the JMA. For the considered aquifers in the Umbria region, the analysis of the performance of the selected reanalyses has confirmed the reliability of the use of the estimates of the flux in the vadose zone to simulate the water table elevation. The executed numerical simulations show the crucial role played by the spreading of the water table elevation with respect to its mean value. In particular, on the basis of the measured water table elevations, two types of aquifers have been identified (Figure 2): the ones characterized by a fast and small-spreading water table elevation (σ w ≤ 1.5 m), and the ones with a slow and large-spreading water table elevation (σ w > 1.5 m). This distinction is a significant factor for defining which reanalysis performs better. Precisely, ERA5 shows a better performance for the aquifers with a slow and largespreading water table, whereas CFS performs better for fast and small-spreading water tables. To a considerable extent, this behaviour can be ascribed to the different depths of the explored soil column. In fact, ERA5, which reaches a larger soil depth (= 2.89 m), produces smoother water fluxes with a slower dynamics. On the contrary, CFS, which reaches a smaller soil depth (= 2 m), produces water fluxes with a faster dynamics as an effect of the larger influence of the surface processes (e.g., precipitation and evaporation). A worse performance characterizes JRA-55 that reaches the smallest soil depth (=1.57 m). Accordingly, it can be affirmed that the explored soil column depth is a significant factor for selecting the most appropriate reanalysis for simulating the underlying groundwater, but, as mentioned, such a feature must be considered in conjunction with the physical characteristics of the water table, that is, the value of σ w .
Besides water table properties and soil column depth, the spatial resolution of the reanalysis is another significant factor to be considered. In the reanalyses, a fine spatial resolution implies not only the availability of grid points closest to the piezometers where water table measurements are executed, but, most significant, it implies a better representation of the land-atmosphere processes. As an example, a finer spatial resolution allows representing properly the topography and, as a consequence, the surface and sub-surface runoff. This is also the reason why global atmospheric reanalyses are usually accompanied by the development of higher resolution land reanalyses such as ERA5-Land for ERA5 (Muñoz Sabater, 2019). The importance of the spatial resolution in the considered case study is confirmed by the fact that JRA-55, the reanalysis with the coarser resolution, shows the worst performance in almost all the considered piezometers. Such a feature could be exalted by the fact that the territory of Umbria is hilly. However, for any type of territory, a fine spatial resolution implies also a better representation of land vegetation, soil texture types, and, therefore, soil hydraulic properties. Finally, for a given soil moisture state, the formulation of the hydraulic conductivity is a key factor for evaluating the flux towards groundwater. Precisely, the Van Genuchten analytic formulation-used by ERA5-gives rise to fluxes towards groundwater of two order of magnitude smaller than those given by CFS and JRA-55, which use the Clapp and Hornberger formulation. With this regard, the executed analysis indicates a better performance of the Van Genuchten relationship on average. Soil hydraulic properties represent a key issue for modelling soil water transport processes within a Richard equation-based approach (Brunone et al., 2003;Romano et al., 2011). The subject remains challenging and under debate even if innovative approaches-for example, the use of bimodal lognormal function (Romano et al., 2011)-are emerging. However, the comparisons executed in Braun and Schadler (2005) indicate that the Van Genuchten model gives the best agreement between observed and simulated soil moisture state on average. Future work will include more observations to investigate whether the above findings can be extended to other areas, with, as an example, a different topography. Indeed, reanalysis performances over a certain area are greatly influenced by the local observations available for being assimilated in the global system and climate peculiarities.