Probability of Non-Exceedance of Arsenic Concentration in Groundwater Estimated Using Stochastic Multicomponent Reactive Transport Modeling

: Stochastic multicomponent reactive transport modeling is a powerful approach to quantify the probability of non-exceedance (PNE) of arsenic (As) critical concentration thresholds in groundwater. The approach is applied to a well-characterized shallow alluvial aquifer near Venice, Italy. Here, As mobility depends primarily on rainfall-controlled redox-dependent precipitation-dissolution of iron hydroxides. A Monte-Carlo analysis based on a calibrated three-dimensional ﬂow and transport model targeted the geochemical initial conditions as the main source of uncertainty of As concentrations in the studied aquifer. It was found that, during 115 simulated days, the fraction of the entire aquifer volume with As > 10 µ gL − 1 decreased on average from ~43% to ~39% and the average As concentration from ~32 µ gL − 1 to ~27 µ gL − 1 . Meanwhile, PNE increased from 55% to 60% when 10 µ gL − 1 was set as target threshold, and from 71% to 78% for 50 µ gL − 1 . The time dependence of As attenuation can be ascribed to the increase of oxidizing conditions during rainfall-dependent aquifer recharge, which causes As sorption on precipitating iron hydroxides. When computing the same statistics for the shallowest 6 m, As attenuation was even more evident. The volume fraction of aquifer with As > 10 µ gL − 1 dropped from 40% to 28% and the average As concentration from 31 µ gL − 1 to 20 µ gL − 1 , whereas PNE increased from 58% to 70% for As < 10 µ gL − 1 and from 71% to 86% for As < 50 µ gL − 1 . Thus, the wells screen depth in the aquifer can be a critical aspect when estimating As risk, owing to the depth-dependent relative change in redox conditions during rainfall events.


Introduction
Geogenic aquifer contamination by arsenic (As) is a common issue in several countries of the world [1][2][3]. An intake of high As concentrations can cause critical diseases such as skin lesions, ulcers or cancer. The World Health Organization (WHO) recommends As concentration in drinking water not to exceed a threshold of 10 µgL −1 (i.e., 10 ppb) [4]. However, up to 220 million people may be potentially exposed to As > 10 µgL −1 [3], the vast majority being in Asia, where groundwater in contact with shallow (<50 m), commonly Holocene sediments can exceed the WHO threshold by a factor of 100 [5]. In many countries, the As concentration threshold for drinking is set to 50 µgL −1 [6]. Yet, a lifetime exposure to just 50 µgL −1 of As in drinking water may kill one in a hundred people prematurely [7].
Most groundwater consumers are unaware of As-related health risks. In alluvial aquifers, shallow private tube wells, sometimes reaching just a few meters below the ground surface, can be easily drilled and ignored by authorities. Therefore, mapping the spatial and temporal (i.e., spatiotemporal) variations of elevated As concentrations in groundwater is of the uttermost urgence for public administrations and regulators. In the European Union, monitoring of groundwater health has been a mandatory task for local administrations since the 2000s, according to the Water Framework Directive (2000/60/EC) and the Groundwater Directive (2006/118/EC).
Reactive transport modeling has been increasingly adopted in recent years to assist decision makers to quantify the potential extension of As contamination in aquifers [8,9]. Specifically, multicomponent reactive transport modeling (MRTM) has been deployed as a tool that enables elucidating and simulating the controls of As mobility in aquifers and predicting its concentrations in space and time [10][11][12][13][14]. These models can mechanistically simulate the effects of As-bearing rocks or sediments, which are naturally found in aquifers.
The development of predictive MRTM to estimate the spatiotemporal variability of biogeochemical reactive species in aquifers remains a complicated task. A key reason is that aquifers are texturally and compositionally heterogeneous across multiple scales. This makes the complete mapping of model parameters incomplete, rendering the predictive model outputs epistemically uncertain [15]. Although averaging (i.e., homogenizing) aquifer parameters is a common practice among modelers, embedding parameter heterogeneities remains fundamental to properly mimic the transport dynamics of reactive dissolved species, such as As [16,17], particularly when performing risk assessment.
When making model-based predictions under uncertainty, stochastic modeling is widely used in hydrogeology to circumvent the limitations of deterministic modeling [18][19][20].
In stochastic analysis, model input parameters are treated as random spatial functions, whereas the model outputs are expressed in terms of probability density functions. Statistical indicators derived from these functions are used as metrics to quantify one or more desired target variables (e.g., the concentration of a polluting aqueous species) [18].
Stochastic MRTM are useful tools that can be used to estimate the probability of nonexceedance (PNE) of an aquifer toxic compound under uncertainty [21]. Quantifying PNE according to different concentration thresholds is a key input for As risk assessment [22,23]. Nonetheless, stochastic reactive transport modeling has been only occasionally applied to the analysis of As transport [24]. We have no knowledge of stochastic MRTM targeting regional scale As mobility in heterogeneous aquifers, as addressed in the present work.
The epistemic uncertainty related to the incomplete mapping of the geochemical initial conditions (GICs) represents a crucial limitation for the reliability of the MRTM predictions. When a system is far from chemical equilibrium, its initial geochemical status is perturbed with time as the flow, transport and geochemical transformations occur. Setting the correct GICs at each model cell is fundamental to properly compute the PNE of a desired toxic compound. Unfortunately, GICs are typically measured at a few locations (i.e., the boreholes in which geochemical analyses were performed) and repeated at time intervals which may exceed the frequency at which contaminants concentrations oscillate and temporarily exceed the regulatory limits. The resulting epistemic uncertainty of the GICs translates into uncertainty in the model outputs. In statistical analysis, the problem is known as error propagation or propagation of uncertainty. In fluid dynamics, error propagation as a function of the uncertain initial conditions has been addressed since the 1960s [25]. However, our literature review revealed that uncertainty in GICs has not been yet widely explored for MRTM-based forecasts, especially for As-related problems.
In this work, we developed a new Monte-Carlo-based stochastic analysis based on three-dimensional (3-D) MRTM. The model, never presented before, was used to compute the time dependent PNE of As concentrations in the shallow Holocene aquifers of the Western Agricultural Areas (WAA). This site is a well-characterized subset of the Venetian Alluvial Plain in Northern Italy ( Figure 1). Previous works on the area [26][27][28][29] have generally concluded that As concentrations can be considered a random spatiotemporal variable, amenable to be modelled using probabilistic and geostatistical analysis. This makes the WAA an appealing case study to evaluate the actual usefulness of the stochastic model to predict the PNE of arsenic concentrations above the 10 µgL −1 and 50 µgL −1 critical thresholds.
Water 2021, 13, x FOR PEER REVIEW 3 of 22 the time dependent PNE of As concentrations in the shallow Holocene aquifers of the Western Agricultural Areas (WAA). This site is a well-characterized subset of the Venetian Alluvial Plain in Northern Italy ( Figure 1). Previous works on the area [26][27][28][29] have generally concluded that As concentrations can be considered a random spatiotemporal variable, amenable to be modelled using probabilistic and geostatistical analysis. This makes the WAA an appealing case study to evaluate the actual usefulness of the stochastic model to predict the PNE of arsenic concentrations above the 10 μgL −1 and 50 μgL −1 critical thresholds. The trace of this section is parallel to the groundwater flow direction depicted in frame (a). Note that the amount of finegrained sediments increases towards the Adriatic sea, forming the multilevel aquifer around the WAA. The figure has been modified after [26].
Our primary goal is to demonstrate the use of a stochastic multicomponent reactive transport model to assess the probability of As non-exceedance, treating the geochemical initial conditions of the MRTMs as the main stochastic parameter. Indeed, for this study a limited and incomplete amount of geochemical data is available to parametrize the initial condition of the MRTMs. More broadly, we aim to showcase the benefit of stochastic MRTM as a decisional tool for administrations and regulators when dealing with complex hydrogeological systems at regional scales. Although our study targets a specific contaminant in a specific aquifer, the developed methodology is general and can be applied to virtually any kind of aquifer and toxic compound.

Background
The Venetian Alluvial Plain (VAP) is a multilevel stratified aquifer underlying the Drainage Basin of Venice Lagoon (in Italian, "Bacino Scolante della Laguna di Venezia") ( Figure 1a). This densely populated region near Venice (Italy) is composed of 108 municipalities, covering an area of about 2000 km 2 and with a population close to 1 million inhabitants (Italian National Institute of Statistics, 2001). Arsenic concentrations in the VAP shallow groundwater can exceed the Italian maximum concentration limit and the WHOrecommended threshold of 10 μgL −1 [26][27][28][29]. Following the EU Water and Groundwater Directive, implemented in Italy mainly through the National Decree D.Lgs 152/2006, an excess of As concentration in groundwater implies that local administration should properly monitor and characterize the extension of the contamination, the natural background levels and thresholds values that can potentially entail a risk for the groundwater consumers. Although the water distribution is supplied to virtually all households in the Drainage Basin of Venice Lagoon, groundwater is mainly used for irrigation purposes and there could be a potential for direct human consumption through private wells. Acronyms refer to the major cities of the Veneto regions. The arrow represents the regional-scale groundwater flow direction. (b) Conceptual not-to-scale vertical sketch of the stratified distribution of the main hydrofacies in the VAP. The trace of this section is parallel to the groundwater flow direction depicted in frame (a). Note that the amount of fine-grained sediments increases towards the Adriatic sea, forming the multilevel aquifer around the WAA. The figure has been modified after [26].
Our primary goal is to demonstrate the use of a stochastic multicomponent reactive transport model to assess the probability of As non-exceedance, treating the geochemical initial conditions of the MRTMs as the main stochastic parameter. Indeed, for this study a limited and incomplete amount of geochemical data is available to parametrize the initial condition of the MRTMs. More broadly, we aim to showcase the benefit of stochastic MRTM as a decisional tool for administrations and regulators when dealing with complex hydrogeological systems at regional scales. Although our study targets a specific contaminant in a specific aquifer, the developed methodology is general and can be applied to virtually any kind of aquifer and toxic compound.

Background
The Venetian Alluvial Plain (VAP) is a multilevel stratified aquifer underlying the Drainage Basin of Venice Lagoon (in Italian, "Bacino Scolante della Laguna di Venezia") ( Figure 1a). This densely populated region near Venice (Italy) is composed of 108 municipalities, covering an area of about 2000 km 2 and with a population close to 1 million inhabitants (Italian National Institute of Statistics, 2001). Arsenic concentrations in the VAP shallow groundwater can exceed the Italian maximum concentration limit and the WHO-recommended threshold of 10 µgL −1 [26][27][28][29]. Following the EU Water and Groundwater Directive, implemented in Italy mainly through the National Decree D.Lgs 152/2006, an excess of As concentration in groundwater implies that local administration should properly monitor and characterize the extension of the contamination, the natural background levels and thresholds values that can potentially entail a risk for the groundwater consumers. Although the water distribution is supplied to virtually all households in the Drainage Basin of Venice Lagoon, groundwater is mainly used for irrigation purposes and there could be a potential for direct human consumption through private wells.
The controls of elevated As concentrations in the VAP have remained unknown for years. Hot spots of dissolved As in groundwater, sometimes exceeding 200 µgL −1 , are variably distributed in the shallow VAP aquifer, without an apparent relationship with high concentrations of As-bearing minerals [26]. Moreover, concentrations measured at a few piezometers were found to fluctuate with time above or below the 10 µgL −1 concentration thresholds. This observation was described during different surveys [26], including the 2009 and 2017 surveys reported in Supplementary Information (Section SI-1) and later analyzed in this work. The erratic occurrence of As excess in groundwater urged environmental regulators to perform systematic analyses and establish collaborations with universities and other research centers to unravel the factors controlling As mobility in the VAP. A target of these activities was to quantify the probability of non-exceedance (PNE) of As concentrations above 10 µgL −1 , to be later used for measuring the risk to which the population could be exposed if As-contaminated groundwater is consumed unwarily.
In a recent study, Dalla Libera et al. [29] (hereafter, "DL2020") proposed a conceptual model that mechanistically explains the spatiotemporal As variations in the VAP aquifer. The study focused on the Western Agricultural Areas (WAA), a well characterized area of the distal portion of the VAP towards the Adriatic Sea. According to DL2020, the aquifer is stratified and composed of multiple sandy layers embedded into a fine-grained matrix ( Figure 1b). Between 6 December 2017 and 26 March 2018 (i.e., over 115 days), weekly geochemical measurements and continuous groundwater head logs were collected from three fully penetrating observation wells (piezometers PZP36, PZP40 and PZP41). The piezometers are drilled to an average depth of about 10 m from the ground surface. DL2020 found that rainfall events were well correlated to the decrease of concentrations of aqueous As and Fe and to the increase of oxidizing-reducing potential (ORP). The shift towards higher ORP values and associated attenuation of As and Fe as a consequence of rainfall events was more evident when a limited thickness (<1 m) of fine-grained sediments was observed above the sandy aquifer, according to the boreholes' stratigraphic logs. When such fine-grained sediments were thicker, the ORP positive shift was more limited, and no significant variations of As and Fe were observed.
Based on these observations, DL2020 hypothesized that rainfall-driven redox-controlled precipitation and dissolution of ferric hydroxides (HFOs) is a main control of As in the shallow VAP aquifers. Similar to the conceptual model proposed for other alluvial aquifers of the world [1,10], under circumneutral pH As can sorb on HFOs surfaces, when the groundwater is sufficiently oxidized to precipitate HFOs. Under reducing conditions, HFOs tend to dissolve, releasing sorbed As to the groundwater. In the VAP, reducing conditions are principally induced by organic matter degradation. Previous studies highlighted that As tends to precipitate as sulfide when groundwater tends to be very reducing (Eh < −200 mV) [27]. Such a situation is however never found in the studied shallow aquifers, where Eh remains between 0 and −200 mV.
A salient aspect of DL2020 s conceptual model was that rainfall could have a strong control on the redox status of the shallow aquifers in the VAP. The key hypothesis was that rainfall could bring oxidants to the shallow aquifers, therefore helping precipitate the HFOs and in turn attenuate As concentrations. Such conceptual model justifies the link between the magnitude of As, Fe and ORP variations and the thickness of the fine-grained sediments, acting as aquitards that partially confine the underlying aquifers. The number of oxidants entering the aquifer could be therefore connected to the effectiveness of aquifer recharge, which depends on the aquitards' hydraulic conductivity (K) and thickness (i.e., the aquitards' leakance).
The conceptual model proposed by DL2020 was validated using a geochemical batch model implemented with the code PHREEQC [30]. The simulated variations of ORP, As and Fe were consistent with the observed values at the boreholes with high-frequency monitoring. The PHREEQC model provided a quantitative and mechanistic way to simulate the impact of rainfall events on As attenuation. Yet, the main limitation of DL2020 model was that a batch (i.e., zero-dimensional) reactor cannot be directly used to estimate the spatio-temporal variations of As concentrations at a continuum, multidimensional scale. Thus, the DL2020 model cannot be used to quantify the probability of not-exceedance (PNE) of As concentration at the scale of the WAA.
The goal of the present work was the development and application of a MRTM to quantify the probability of As exceeding the 10 µgL −1 threshold as well as another relevant thresholds for As risk, the 50 µgL −1 , at the scale of the WAA, bypassing the limitations of DL2020 model. Leveraging the knowledge gained by DL2020, the new model aims at mechanistically predicting As mobility in the multidimensional heterogeneous WAA, testing the capability of DL2020 conceptual model to quantify the PNE of As concentrations at a wider, regional scale.

Model Setup
We utilized PHAST [31,32], an operator-splitting MRTM code that combines the flow and conservative transport code HST3D [33] and PHREEQC. The selection of PHAST was appropriate for this study for multiple reasons. PHAST shares the same structure as PHREEQC; thus, it is able to reproduce mechanistically the full range of equilibrium and kinetic geochemical reactions describing As mobility in a multidimensional dynamic system. As the input files are similar in PHAST and PHREEQC, the same geochemical setup developed by DL2020 could be used for the PHAST models. Additionally, PHAST includes a free-surface boundary condition that can simulate groundwater table oscillation in confined-unconfined aquifers. This was important for our study area where the rainfalldriven areal recharge and a large drainage system overlap to change the confining status of the model's layer in space and time. PHAST allows for cells rewetting when the hydraulic head of a dry cell raises above the cell's bottom elevation, and the input flux (both water and dissolved solutes) is applied to the highest active cells. Ultimately, PHAST is equipped with a parallel computation algorithm allowing for runs on a multiprocessor computer or on a collection of computers that are networked.
A main aim of this study was to generate a stochastic MRTM analysis with PHAST, treating the geochemical initial conditions (GICs) as uncertain parameters. Therefore, we proceeded through the following sequential steps: • A 3-D flow model embedding a heterogeneous distribution of lithological materials was setup to simulate groundwater flow in the WAA. • The flow model was used for the advective and dispersive transport of a multidimensional MRTM model. This model inherited the geochemical reactions proposed by DL2020 to describe As mobility in the WAA. • Within a Monte-Carlo framework, indicator-based geostatistical modeling was used to build spatially correlated random realizations of GICs. Such GICs were employed to run independent MRTM simulations.

•
The ensemble of results was then collected and examined to estimate the PNE and useful statistics that allow gaining insight into the expected variation of As and related key variables in the WAA.

Flow Model
The permost active cells in order to simulate the areal recharge due to rainfall events. Two recharge rate zones (RCH1 and RCH2) were considered, due to the different amounts of fine-grained materials on top of the aquifer. A third-type DRAIN (DRN) BC was set in the central zone of the WAA to simulate an existing drainage system that lowers the groundwater table in the area (as part of the agricultural land conversion initiated in the early 1900s). All BCs were set as time-dependent functions, except for the DRN which has no time-dependent parameters. To map the spatially-variable hydraulic properties of the WAA, we followed the methodology presented by Fabbri et al. [34] who developed a 3-D mapping of subsoil heterogeneity in a different, northern sector of the VAP. The methodology was based on the code "spMC" [35,36], a transition probability algorithm that simulates the juxtapositions of categorized variables (here, the "lithological materials") through embedded Markov chains. Stochastic lithological maps conditioned to the stratigraphic data available in the site are generated by spMC. Details regarding the construction of the lithological map are provided in the Supplementary Information (SI-2). In short, stratigraphic data are classified and categorized into four lithological materials: "sand", "silt", "clay" and "peat". Although more complexity could be used to classify soil materials, the adopted codification is based on the average grain size distribution of the sediments, which ultimately control the soil hydraulic conductivity (K). Simplified ways to classify soil type for transitional probability based geostatistical modeling have been commonly adopted in the past [37]. The code spMC computed internally the transiograms, which are used to produce a stochastic realization of the four materials in the studied area. Such realization, depicted in Figure 3, was then re-mapped by assigning individual hydraulic conductivity (K), specific storage (Ss) and specific yield (Sy) to each of the four materials. The resulting K, Ss and Sy maps were then used to populate the PHAST model. The porosity (ϕ) was assumed as homogeneously distributed, varying less than the other relevant hydraulic parameters The flow boundary conditions were set as follows. Along the lateral edges, Cauchytype General Head Boundaries (GHB) were assigned. Along the Western and Northern edges, the GHB1 and GHB2 simulated the lateral groundwater recharge from upgradient aquifer bodies. The GHB4 assigned along the Southern edge simulated the interaction between the study aquifer and the Naviglio Brenta River flowing in that area. The GHB3 located on the Eastern edge simulated the hydraulic interaction between the studied aquifer and the Venice Lagoon. A Neuman-type Recharge (RCH) BC was applied to the uppermost active cells in order to simulate the areal recharge due to rainfall events. Two recharge rate zones (RCH1 and RCH2) were considered, due to the different amounts of fine-grained materials on top of the aquifer. A third-type DRAIN (DRN) BC was set in the central zone of the WAA to simulate an existing drainage system that lowers the groundwater table in the area (as part of the agricultural land conversion initiated in the early 1900s). All BCs were set as time-dependent functions, except for the DRN which has no time-dependent parameters.
To map the spatially-variable hydraulic properties of the WAA, we followed the methodology presented by Fabbri et al. [34] who developed a 3-D mapping of subsoil heterogeneity in a different, northern sector of the VAP. The methodology was based on the code "spMC" [35,36], a transition probability algorithm that simulates the juxtapositions of categorized variables (here, the "lithological materials") through embedded Markov chains. Stochastic lithological maps conditioned to the stratigraphic data available in the site are generated by spMC. Details regarding the construction of the lithological map are provided in the Supplementary Information (SI-2). In short, stratigraphic data are classified and categorized into four lithological materials: "sand", "silt", "clay" and "peat". Although more complexity could be used to classify soil materials, the adopted codification is based on the average grain size distribution of the sediments, which ultimately control the soil hydraulic conductivity (K). Simplified ways to classify soil type for transitional probability based geostatistical modeling have been commonly adopted in the past [37]. The code spMC computed internally the transiograms, which are used to produce a stochastic realization of the four materials in the studied area. Such realization, depicted in Figure 3, was then re-mapped by assigning individual hydraulic conductivity (K), specific storage (S s ) and specific yield (S y ) to each of the four materials. The resulting K, S s and S y maps were

Deterministic Transport Model
According to DL2020, As mobility in the WAA is primarily controlled by redox-de pendent precipitation and dissolution of HFOs, on which As can sorb. Arsenic sorptio was simulated through a surface complex model (SCM) [10] that provides a mechanisti

Deterministic Transport Model
According to DL2020, As mobility in the WAA is primarily controlled by redoxdependent precipitation and dissolution of HFOs, on which As can sorb. Arsenic sorption was simulated through a surface complex model (SCM) [10] that provides a mechanistic quantification of As adsorption behavior as a function of the dissolved As and Fe concentrations. The amount of reactive surfaces (i.e., the HFOs) was assumed to be a proportion of the amorphous Fe-hydroxides (Fe(OH) 3 ). The SCM considers the influence of competing ions such as phosphate or bicarbonate as well as the density of adsorption sites on the host minerals. The thermodynamic database WATEQ4F [38] provided the reference stability constants for HFOs equilibrium based on Dzombak and Morel [39].
The stability of HFOs depends directly on the change in pH and redox conditions. The system pH was observed to remain stable around circumneutral values due to the buffering by calcite, which is present in great amounts in the studied area and included in the reaction network as an equilibrium phase. The change in redox conditions affect the stability of Fe(OH) 3 by altering the Fe 2+ /Fe 3+ redox couple. When the aquifer becomes more oxidized (i.e., the ORP increases), Fe is more prone to pass to a higher oxidation status (Fe 2+ → Fe 3+ ) and precipitate as amorphous Fe(OH) 3 . Since HFOs are directly linked to the amount of Fe(OH) 3 , a change in ORP intrinsically leads to a change in the concentration of surfaces on which As can sorb. Fe(OH) 3 was introduced to the reaction network as an equilibrium phase.
Rainfall events alter the redox status of the groundwater by recharging the aquifer with oxidized water. The oxidant mass flux entering the aquifer is computed internally by PHAST from the recharge rates imposed to the RCH boundary and from the concentrations assigned to that boundary, which is in direct contact with the atmosphere. DL2020 noted that the ORP variations could not be explained considering solely the ingress of O 2(aq) in equilibrium with the atmosphere. It was concluded that infiltrating water could be enriched in other dissolved oxidizing species (e.g., nitrates or sulfates generated by agriculture practices and leached through the topsoil). Alternatively, gas-phase oxygen could be pushed through the vadose zone, reaching the aquifer. We lack specific information regarding the dominant mechanism controlling the excess of oxidants in the infiltrating water. Since the goal of this work is to present a general tool to compute the excess of As at the regional scale using a probabilistic approach, we adopted a simplified approach and equilibrated the infiltrating water with a much higher O 2 concentration than the atmospheric gas. Hence, the infiltrating boundary brings a much higher amount of O 2 molecules than in real life, while not adding other species to the aquifer. Note that this assumption is not a limitation of the study, since we still preserve the very essence of the study, i.e., evaluating the probability of As non-exceedance as a function of the uncertain input factors (the geochemical initial conditions), as carefully explained in the next section.
The dispersivity coefficient (α) was treated as a homogeneous parameter. Lacking specific tracer tests, we followed empirical compilations of results [40] to set the longitudinal dispersivity (α x ) to 100 m, the horizontal transverse dispersivity (α y ) to 10 m and the vertical transverse dispersivity to 1 m, The chosen dispersivities satisfy the stability condition by which the Peclet number (Pe), such that Pe = ∆x i /α i ≤ 2, where ∆x defines the cell size and i the direction (i = x, y, z). The effective diffusion coefficient (D) was set to D = 10 −9 m 2 /s.

Stochastic Transport Model
PHAST requires a "solution" to be defined at each model cell. In this context, a "solution" is a PHREEQC modeling term [30] that refers to a list of concentrations that defines the geochemical chemical composition of the water at the beginning of the simulation. Since the geochemical composition of the aquifer is heterogeneous, a spatial distribution (i.e., a map) of solutions is required to parameterize each PHAST model cell. Our simulations assumed December 2017 as the initial time; therefore, the map of solutions should refer to December 2017. The December 2017 geochemical database contains measured concentrations at the eight boreholes shown in Figure 2a and reported in the Supplementary Information SI-1. The database was obtained from the Veneto region environmental agency (ARPAV) during the routine monitoring activities to monitor groundwater quality. Although a larger number of boreholes exist in the area, most of them were either dismissed with time or not included in the recent monitoring activities. The December 2017 dataset is broad in terms of variety of analyzed species. This allows setting a complete PHREEQC solution at least at the eight model cells that spatially correspond to the position of eight boreholes sampled by ARPAV. Sampling and analytical methods followed the standard protocols, as suggested by the Italian National Environmental Agency. Detail regarding these methods can be found in the Supplementary Information SI-1.
The limited number of sampled boreholes contained in the December 2017 prevented us from obtaining a "deterministic" distribution of solutions at each model cell. Geostatistical analysis could be invoked for the statistical interpolation of concentration maps [23,26,41]. Such maps can be conditioned to the existing geochemical measurements (in geostatistical terms, the "hard data"). Yet, a sound geostatistical inference cannot be performed using the December 2017 dataset. Too few spatially distributed data impede calculating the spatial structure of the geochemical data (e.g., the concentrations variograms), which is the crucial step to perform geostatistical analysis.
To circumvent this limitation, we considered an older geochemical characterization carried out by ARPAV in 2009. The 2009 characterization targeted about 50 piezometers, of which 34 were screened in the shallow portions of the VAP aquifer. The corresponding database is reported in the Supplementary Information SI-1. The 34 sampling points are enough to perform geostatistical inference and obtain information regarding the correlation of concentrations in the shallow aquifer [26]. However, in 2009 only specific contaminants were sampled, the majority being heavy metal(loid)s. The reason is that the 2009 sampling campaign was mainly focused on evaluating the extension of the groundwater contamination linked to the Porto Marghera chemical factory, an Italian site of major interest located in the proximity of the WAA [42]. While the 2009 database included key elements such as As and Fe, several other major, minor and trace elements were not measured. This issue implies that a complete list of solutions cannot be obtained from the 2009 dataset. Although a MRTM can be run using only a limited number of species, not including key geochemical species that could directly or indirectly control the fate of As limits the reliability of the stochastic MRTM outputs. We also considered that, even if we had a sufficiently complete set of geochemical species in the 2009 dataset, we could not directly use this dataset to simulate transport starting in December 2017. Indeed, according to our conceptual model, As mobility is controlled by aquifer recharge and our calibrated flow model refers specifically to the 2017-2018 period.
We combined pros and cons of the 2009 and 2017 databases to obtain stochastic geochemical initial conditions (GICs) for the PHAST models, as follows: 1.
The spatial structure (i.e., the variogram) of As was initially computed from the 34 boreholes sampled in 2009. The As variogram was used to generate N mc = 100 sequential indicator simulations (SISs) [43,44] of As. The simulations were created on a 2D grid having the same size as the PHAST models (43 × 46 cells).

2.
These simulations (or "realizations") were conditioned to the 2017 measurements of As and obtained from the eight sampled piezometers. Examples of maps generated using the SIS algorithm are shown in the Supplementary Information SI-4.

3.
The procedure at steps 1 and 2 was repeated for Fe, to obtain N mc = 100 Fe realizations.

4.
We treated the eight hydrogeochemical surveys performed in 2017 as eight PHREEQC "solutions". We then created N mc = 100 empty grids of the same size of the stochastic realizations (43 × 46 cells) and populated them with a unique solution in each grid cell. Such solution was chosen by minimizing, for each cell of each j-th pair of As and Fe realizations (j = 1, . . . , N mc ), the pairwise Euclidean distance between the simulated As and Fe concentrations and the observed concentrations from the 8 PHREEQC solutions.
In this way, we obtained N mc = 100 random maps of solutions (i.e., the geochemical initial conditions, GICs), which were used to populate the PHAST models. The mapping algorithm was coded through an in-house MATLAB script, which included the native pdist function to calculate pairwise Euclidean distance. The variograms and SIS computation was performed using SGEMS [45]. The resulting variograms are reported in the Supplementary Information SI-4.
An example of the approach is graphically depicted in Figure 4a. The map on the left represents a realization of As concentrations, while the map on the right represents a realization of Fe concentrations. In the cell marked by the green square, the SIS algorithm generated a concentration of As = 0.023 mgL −1 and of Fe = 2.214 mgL −1 , shown respectively in the map on the left and on the right. The closest pairwise distance between measured and computed As and Fe concentrations in this specific cell correspond to the concentrations represented by the "Solution 4". In the cell marked by the magenta square, the concentrations of As = 0.012 mgL −1 and Fe = 1.354 mgL −1 are instead closer to the "Solution 6". For each cell (i.e., x-y position), we assigned the same solution in all layers, since the composition of the groundwater which is sampled by fully penetrating observation wells is an average aquifer composition of the entire borehole.
The process was repeated for each domain cell until the map of solution was completed. The whole procedure was then looped N mc = 100 times to obtain N mc = 100 solution maps. Examples of two solution maps obtained using this approach are shown in Figure 4b. Each of the N mc = 100 solutions maps were used to populate and run independent PHAST simulations. After each run, the entire spatial and temporal distribution of As and other key variables was stored. The ensemble of results was analyzed in a probabilistic fashion. For a generic variable X, the ensemble mean X and standard deviation σ X were computed to evaluate the change in time of the expected behavior of geochemical system under uncertainty. Empirical cumulative distribution functions (ecdfs) were computed and compared at different time scales to calculate the PNE above key thresholds of As concentrations.
We finally recalled that groundwater could be withdrawn through "fully penetrating" pumping wells, i.e., wells which are screened along the entire 12-m-long vertical section of the aquifer, or through partially penetrating well, and in particular through wells are screened at shallower depths. We postulated that a different pumping well configuration could result in a different PNE and become a useful decisional criterion for local administration, for instance to compute well-type-specific As risk for groundwater consumers. To demonstrate this, we conducted two distinct analyses:

Flow and Deterministic Transport Model
The subsurface lithological model created using spMC highlighted a large variability of lithological units within the WAA, both in terms of relative proportion and in terms of spatial continuity of such units. Volumetrically, the statistical analysis of available borehole logs (Figure 5a) indicates that about one third of the aquifer is composed of sandy materials (35% of the total volume), while the remaining two thirds are composed of finegrained, low-permeable sediments (46% silt, 18% clay and 1% peat). The mean length of

Flow and Deterministic Transport Model
The subsurface lithological model created using spMC highlighted a large variability of lithological units within the WAA, both in terms of relative proportion and in terms of spatial continuity of such units. Volumetrically, the statistical analysis of available borehole logs (Figure 5a) indicates that about one third of the aquifer is composed of sandy materials (35% of the total volume), while the remaining two thirds are composed of fine-grained, low-permeable sediments (46% silt, 18% clay and 1% peat). The mean length of the materials along the x (i.e., E-W) direction shows (Figure 5b) similar lengths for sand and silt (respectively, 284 m and 297 m), while clay (122 m) and particularly peat (16 m) materials show shorter lengths. Comparing the lengths in the x direction with those computed for the y (N-S) direction (Figure 5c) reveals that sandy materials are quite isotropic in the horizontal direction, while silty material is slightly more anisotropic and elongated towards the x direction. Peat and particularly clay materials seem much more continuous along the y direction, when compared to the corresponding values in the x direction. For all materials, however, the mean lengths along the horizontal plane are much longer than in the vertical direction z (Figure 5d). Lengths of about 2 m were found for the sandy materials (the main aquifers of the WAA), with a ratio between the horizontal and vertical lengths close to 150.  (Figure 5c) reveals that sandy materials are quite isotropic in the horizontal direction, while silty material is slightly more anisotropic and elongated towards the x direction. Peat and particularly clay materials seem much more continuous along the y direction, when compared to the corresponding values in the x direction. For all materials, however, the mean lengths along the horizontal plane are much longer than in the vertical direction z (Figure 5d). Lengths of about 2 m were found for the sandy materials (the main aquifers of the WAA), with a ratio between the horizontal and vertical lengths close to 150. The flow model indicates a mean groundwater velocity of about 0.2 m/d. The stratigraphic variation determines a spatial variability of the hydraulic parameters, which translates into a variability in the local flow velocity. Table 1 reports the resulting hydraulic parameters at the end of the calibration process. The sandy materials show a larger K than the finer-grained materials, suggesting that sandy horizons could be acting as preferential flow zones, enhancing the mixing between infiltrating waters and existing water. In the northern part of WAA, the calibrated recharge rates associated to the transient RCH BCs linearly decreased from an initial value of 10% of the total rainfall to 1% after 97 days and remained at 1% until the end of the simulation. In the southern and in the middle zones, the calibrated recharge rate linearly decreased from 10% to 3% after 60 days and remained at 3% until the end of the simulation. This range of values agrees with the recharge rates reported for the shallow aquifers of the adjacent Porto Marghera site [42]. The flow model indicates a mean groundwater velocity of about 0.2 m/d. The stratigraphic variation determines a spatial variability of the hydraulic parameters, which translates into a variability in the local flow velocity. Table 1 reports the resulting hydraulic parameters at the end of the calibration process. The sandy materials show a larger K than the finer-grained materials, suggesting that sandy horizons could be acting as preferential flow zones, enhancing the mixing between infiltrating waters and existing water. In the northern part of WAA, the calibrated recharge rates associated to the transient RCH BCs linearly decreased from an initial value of 10% of the total rainfall to 1% after 97 days and remained at 1% until the end of the simulation. In the southern and in the middle zones, the calibrated recharge rate linearly decreased from 10% to 3% after 60 days and remained at 3% until the end of the simulation. This range of values agrees with the recharge rates reported for the shallow aquifers of the adjacent Porto Marghera site [42].  Figure 6 shows the interpolated map of resulting head-level distribution at the end of the simulated period. It can be observed that the head levels are strongly affected by the presence of the drainage systems, which locally lower the water table below the sea level. At the stage depicted by this specific snapshot, the minimum head level is at −1.37 m above sea level (a.s.l.). Considering that the model elevation is 3 m a.s.l. and the layer discretization is 1 m, at the end of the simulated time the top four layers in the proximity of the drainage system are dry. As mentioned earlier, however, this issue does not imply a limitation for the reactive transport model, since the input mass flux associated to the recharge (a key mechanism for the control of As in our study) still applies to the top active layers.   Figure 6 shows the interpolated map of resulting head-level distribution at the end of the simulated period. It can be observed that the head levels are strongly affected by the presence of the drainage systems, which locally lower the water table below the sea level. At the stage depicted by this specific snapshot, the minimum head level is at −1.37 m above sea level (a.s.l.). Considering that the model elevation is 3 m a.s.l. and the layer discretization is 1 m, at the end of the simulated time the top four layers in the proximity of the drainage system are dry. As mentioned earlier, however, this issue does not imply a limitation for the reactive transport model, since the input mass flux associated to the recharge (a key mechanism for the control of As in our study) still applies to the top active layers.

Whole Aquifer Depth
We initially analyzed the statistical change in relevant geochemical species and variables at the scale of the entire model, i.e., considering the entire vertical extension of the aquifer. Key results of this analysis are shown in Figure 7. For each variable, the blue line represents the ensemble mean obtained from averaging over the 100 stochastic simulations, while the red dotted lines represent the 95% (±2σ) uncertainty envelope. The grey bars represent the recharge rates (ms −1 ), which are directly correlated to the rainfall events.
Our first significant result is the relative volumetric proportion (the fraction) of aquifer where As concentrations exceeding 10 μgL −1 decrease with time. Initially, ~43% of the

Whole Aquifer Depth
We initially analyzed the statistical change in relevant geochemical species and variables at the scale of the entire model, i.e., considering the entire vertical extension of the aquifer. Key results of this analysis are shown in Figure 7. For each variable, the blue line represents the ensemble mean obtained from averaging over the 100 stochastic simulations, while the red dotted lines represent the 95% (±2σ) uncertainty envelope. The grey bars represent the recharge rates (ms −1 ), which are directly correlated to the rainfall events. The ensemble means of the ecdfs are reported in Figure 8c, along with the calculated mean As concentrations already reported in the previous panels. It was found that the expected PNE for the WHO-recommended threshold of 10 μgL −1 was 55% at t = 0, and increased to 60% at t = 115 d. This means that, on average, the fraction of aquifer that can be considered contaminated by arsenic at the end of the simulation was lower by 5% compared to the beginning of the simulation. For 50 μgL −1 , another common threshold which is sometimes a maximum concentration in drinking water in countries such as Bangladesh, the probability of non-exceedance has increased from 71% at t = 0 to 78% at t = 115 d.
The time dependance of the PNE is consistent with the rest of the analysis and linked to the effect of As attenuation as controlled by the adsorption of HFOs. Attenuation increases as time elapses, as a consequence of the hypothesized effects of rainfall-controlled HFOs in the studied aquifer. Since PNE is directly connected to risk indicators, such as to compute the chronic exposure of population to an excess of a toxic waterborne compound, our results stress the importance of rainfall events as a critical aspect controlling the risk of As contamination in the studied area. This is important information for decision makers when quantifying the risk of As contamination in the studied area. Since the PNE is nonstationary in time, administrations should regularly update their risk calculations depending on new concentration levels. Our first significant result is the relative volumetric proportion (the fraction) of aquifer where As concentrations exceeding 10 µgL −1 decrease with time. Initially,~43% of the entire aquifer is characterized by As concentrations above 10 µgL −1 . As time elapses, the fraction drops toward a final value of~39%. This behavior is linked to a decrease in average As concentration which tends to drop from~32 µgL −1 to~27 µgL −1 . Although the model was not calibrated using geochemical information, the decreasing trend in As concentration simulated by the stochastic model is consistent with the previous results by DL2020. Since As and Fe are strongly correlated, dissolved Fe concentrations also drop with time, following a similar pattern as As concentrations. The drop in dissolved Fe is opposed to an increase in reactive surfaces (HFOs), whose amount depends directly on the amount of amorphous iron Fe(OH) 3 existing in the system. The concentration of amorphous iron is inversely correlated to the amount of dissolved Fe. This is explained considering that, as the time proceeds, the ensemble mean of saturation index of Fe(OH) 3 tends towards higher values, i.e., SI → 0. In other words, as time elapses the system tends towards equilibrium conditions for Fe(OH) 3 .
The main driver of such reduction is the aquifer recharge. This is demonstrated by the clear correlation between rainfall events and the oxidation reduction potential (ORP), here expressed by the pe = − log(e+). Initially, ORP tends to higher values; then, a sudden decrease in ORP is found at~40 days, corresponding to the beginning of a dry period lasting several consecutive days. After that, new precipitation events occurred and ORP increased again. This is explained considering that precipitation triggers aquifer recharge, which adds oxidants to the aquifers. All the other redox-sensitive species displayed in Figure 7, including As, show a similar "bumped" behavior. Note instead that pH remains close to circumneutral conditions throughout the entire simulation time. This is consistent with the buffering of the calcite minerals. Figure 8 shows the empirical cumulative distribution functions (ecdfs) of As concentrations obtained from two different time steps of the MRTM simulations. In Figure 8a, the ensemble of ecdf calculated for each realization reflects the distribution of As at the beginning of the simulation (t = 0). The ecdf indicates that the distribution of As in the whole aquifer is very heterogeneous. The use of the log scale helps to better capture the variations of As concentrations which span over more than three orders of magnitude and deviate largely from the mean value of~32 µgL −1 (blue dotted vertical line). The shape of the ecdfs is alike among the different realizations. It is characterized by a stepwise increase in the cumulative probability, as a result of the SIS-based stochastic realizations used as GICs in the reactive transport models Figure 4. The SIS realizations include eight categories (i.e., the "solutions"), resulting in an equivalent number of steps, since each solution occupies a different proportion of the aquifer volume.

Shallower Aquifer Portion
We now reanalyze the results from the stochastic models using a similar approach, this time focusing on the top six meters of the aquifers. Figure 9 shows the changes in the same geochemical species and variables evaluated in Figure 7. Compared to variation occurring when integrating the results of the entire aquifer thickness, we found a more In Figure 8a, the ecdf reflects the distribution of As at the end of the simulation (t = 115 d). Note that the ensemble of ecdfs is now different than the ensemble of ecdfs calculated for t = 0. This is due to a combination of factors. First, transport determines a mixing among adjacent cells, such that the well-defined clustered As distribution at t = 0 is now more continuous, as reflected by smoother ecdf obtained at t = 115 d. Secondly, as time elapses, As is attenuated by the rainfall-driven redox-controlled adsorption on HFOs. Indeed, the ensemble of ecdfs at t = 115 d is shifted towards lower values, being consistent with the decrease in the mean As concentration to~27 µgL −1 .
These ecdf can be used to quantify the PNE of arsenic concentrations according to targeted thresholds. We note however that such estimation is hardly obtained from Figure 8a,b, where the ensemble of ecfds is quite spread. This is particularly true for values exceeding a concentration of 3 µgL −1 . The ecfds spreading is explained considering that the relative proportion of As in the GICs is random outcome of the SIS algorithm. While on average all realizations embed the same statistics, each individual realization can produce a slightly different proportion of a specific solution. As such, we calculated the ensemble mean of the ecdfs to obtain a unique reference PNEs for the two analyzed times reported in Figure 8a,b.
The ensemble means of the ecdfs are reported in Figure 8c, along with the calculated mean As concentrations already reported in the previous panels. It was found that the expected PNE for the WHO-recommended threshold of 10 µgL −1 was 55% at t = 0, and increased to 60% at t = 115 d. This means that, on average, the fraction of aquifer that can be considered contaminated by arsenic at the end of the simulation was lower by 5% compared to the beginning of the simulation. For 50 µgL −1 , another common threshold which is sometimes a maximum concentration in drinking water in countries such as Bangladesh, the probability of non-exceedance has increased from 71% at t = 0 to 78% at t = 115 d.
The time dependance of the PNE is consistent with the rest of the analysis and linked to the effect of As attenuation as controlled by the adsorption of HFOs. Attenuation increases as time elapses, as a consequence of the hypothesized effects of rainfall-controlled HFOs in the studied aquifer. Since PNE is directly connected to risk indicators, such as to compute the chronic exposure of population to an excess of a toxic waterborne compound, our results stress the importance of rainfall events as a critical aspect controlling the risk of As contamination in the studied area. This is important information for decision makers when quantifying the risk of As contamination in the studied area. Since the PNE is non-stationary in time, administrations should regularly update their risk calculations depending on new concentration levels.

Shallower Aquifer Portion
We now reanalyze the results from the stochastic models using a similar approach, this time focusing on the top six meters of the aquifers. Figure 9 shows the changes in the same geochemical species and variables evaluated in Figure 7. Compared to variation occurring when integrating the results of the entire aquifer thickness, we found a more marked decrease in As concentrations with time in the shallow portions of the aquifer. The relative aquifer fraction with As >10 µgL −1 dropped from 40% to 28%, which is linked to a decrease in average As concentrations from 31 µgL −1 at t = 0 to 20 µgL −1 at t = 115 d. A higher relative change was also observed for the other key parameters controlling As mobility. of these curves. While the distribution of As values at t = 0 is similar to the distribution calculated for the full-depth analysis (Figure 8c), at t = 115 d the ensemble-mean ecdf is smoother than in the previous case. The probability computed for different thresholds also changed. For As < 10 μgL −1 , the probability increased from 58% at t = 0 to 70% at t = 115 d, while for As < 50 μgL −1 , the PNE increased from 71%% at t = 0 to 86% at t = 115 d. Thus, the relative risk due to As excess in drinking water is more reduced if groundwater is withdrawn at shallower depths, owing to the better attenuating effects of rainfall-driven recharge on As. Figure 9. Key results of the stochastic transport analysis. This figure refers to the "shallow depth" analysis. For each variable, the blue lines reproduce the ensemble mean obtained from averaging over the 100 stochastic simulations, while the red dotted lines represent the 95% (±2σ) uncertainty envelop. The gray bars represent the recharge rates (ms −1 ).

Figure 9.
Key results of the stochastic transport analysis. This figure refers to the "shallow depth" analysis. For each variable, the blue lines reproduce the ensemble mean obtained from averaging over the 100 stochastic simulations, while the red dotted lines represent the 95% (±2σ) uncertainty envelop. The gray bars represent the recharge rates (ms −1 ).
The difference between the results of the shallower portions of the aquifer and those of the full aquifer is explained considering that the shallow aquifer is closer to the recharge boundary. Therefore, the recharge-driven redox changes are enhanced in the shallower parts of the aquifer than in the deeper part of it where groundwater remains more isolated from the ingress of oxidants. The relative control of recharge on As attenuation explains why the probability that As does not exceed a specific threshold also depends on the depth of the analyzed aquifer. Figure 10a,b show the ensemble of ecdfs for the top six meters of the aquifer calculated at t = 0 and t = 115 d, respectively, and Figure 10c the ensemble mean of these curves. While the distribution of As values at t = 0 is similar to the distribution calculated for the full-depth analysis (Figure 8c), at t = 115 d the ensemble-mean ecdf is smoother than in the previous case. The probability computed for different thresholds also changed. For As < 10 µgL −1 , the probability increased from 58% at t = 0 to 70% at t = 115 d, while for As < 50 µgL −1 , the PNE increased from 71%% at t = 0 to 86% at t = 115 d. Thus, the relative risk due to As excess in drinking water is more reduced if groundwater is withdrawn at shallower depths, owing to the better attenuating effects of rainfall-driven recharge on As.

Discussion
Stochastic modeling is a powerful method to compute the probability of non-exceedance (PNE) of toxic compounds in complex aquifers. While several models have been developed to compute the probability of As contamination in different environments [46], for groundwater such probability has been traditionally computed using data-driven statistical analysis [22]. For instance, the Bayesian Maximum Entropy theory was applied by Serre et al. [47] to study arsenic exposure across Bangladesh by the local population using a database that combined about 4000 tube-well and borehole samples from four different surveys and a total of 3373 As concentration measurements of the shallow aquifers. When such datasets do not exist, data-driven statistical analysis may be insufficient to perform probabilistic analysis.
The present study is among the first documented attempts to evaluate PNE for As critical concentrations in groundwater using a mechanistic MRTM formulation. Leveraging previous studies on the same area, the model accounted for key hydraulic and geochemical controls of As mobility in a shallow aquifer near Venice, Italy. Although stochastic mechanistic modeling coupled to geostatistical inference is a demonstrated tool that can generate probabilistic maps of time-dependent concentrations of toxic compounds

Discussion
Stochastic modeling is a powerful method to compute the probability of non-exceedance (PNE) of toxic compounds in complex aquifers. While several models have been developed to compute the probability of As contamination in different environments [46], for groundwater such probability has been traditionally computed using data-driven statistical analysis [22]. For instance, the Bayesian Maximum Entropy theory was applied by Serre et al. [47] to study arsenic exposure across Bangladesh by the local population using a database that combined about 4000 tube-well and borehole samples from four different surveys and a total of 3373 As concentration measurements of the shallow aquifers. When such datasets do not exist, data-driven statistical analysis may be insufficient to perform probabilistic analysis.
The present study is among the first documented attempts to evaluate PNE for As critical concentrations in groundwater using a mechanistic MRTM formulation. Leveraging previous studies on the same area, the model accounted for key hydraulic and geochemical controls of As mobility in a shallow aquifer near Venice, Italy. Although stochastic mechanistic modeling coupled to geostatistical inference is a demonstrated tool that can generate probabilistic maps of time-dependent concentrations of toxic compounds [18,21], it is not clear why stochastic MRTM have not been widely documented so far in the literature. One reason may owe to the well-known computational burden required to resolve the nonlinear equations characterizing this type of model, which can determine very long computational times. A second reason may be linked to the number of unknowns and input parameters required to run MRTM. However, the present study showed that stochastic MRTM of As mobility can be computed with reduced computational burden (all simulations run overnight on a modern workstation) if the model embeds a simplified geochemical system that captures the essential mechanisms controlling As mobility.
Our stochastic analysis revealed that the shallow aquifers in the WAA are more sensitive to rainfall-controlled reactions than deeper parts of the aquifers. On the one hand, this corroborates and extends the validity of the conclusions by Dalla Libera et al. [29]. On the other hand, the results are well aligned with the conclusions from other studies that already identified a link between As mobility and aquifer recharge. Rodriguez et al. [48] explained the variation in As concentration in terms of the precipitation regime and the oxidation conditions induced by local rainfall incorporation. Our model results could be interpreted as if relatively shallow wells are less prone to generate As risk than relatively deeper wells. We stress the term "relative", as this assessment refers uniquely to wells which are drilled up to a maximum depth of 10-15 m. In the Venice area, other aquifers located at much deeper depths exist and are utilized for drinking purposes.
It is also important to note that rainfall is not always linked to an attenuation of As concentrations. For instance, Duan et al. [49] found that rainfall implied an increase in As concentrations in Jianghan Plain in China. The authors suggested that rainfall could generate increasingly reducing conditions by increasing the groundwater table, while the dry season favors the ingress of oxidants and the adsorption of As of HFOs. Therefore, a sound evaluation of the site-specific mechanisms controlling As mobility is required prior to performing any stochastic MRTM.
Ultimately, we note that despite being based on a 3-D mapping of the subsoil, flow boundaries, recharge conditions, and the variations of the hydraulic parameters, the flow model was calibrated using the values of hydraulic heads measured in three aligned points, and the results were extended to the studied area. Future developments of the model should include more spatially distributed observation points, when available, in order to increase the validity of the model results.

Conclusions
We conducted a stochastic analysis based on multidimensional multicomponent reactive transport modeling (MRTM) to evaluate the probability that As does not exceed critical thresholds in the Western Agricultural Area (WAA). This area is a subset of the Venetian Alluvial Plain (VAP), near Venice (Italy), a densely populated region where As concentration in shallow Holocene aquifers (up to a depth of 10-15 m) exceeds the WHOrecommended threshold of 10 µgL −1 . The complex spatial and temporal variations of As concentrations complicate the decisions regarding control and mitigation of As risk in the area.
The analysis suggested that the proposed modeling approach is able to quantify such probability. When the entire modeled aquifer is analyzed in a Monte-Carlo-based stochastic fashion, the relative proportion of the aquifer where As > 10 µg/L tends to decrease with time from an initial~43% to a final~39% after 115 simulated days. The probability of As non-exceedance increased from 55% to 60% with time using the 10 µgL −1 as target threshold, and from 71% to 78% when the target was 50 µgL −1 . The result corroborates the conclusions of a previous study [29] who proposed that rainfall could add oxidants to the shallow aquifer and promote the formation of sorption surfaces, on which As can be temporarily removed from the aqueous environment.
When the same statistics are computed for the shallowest (top 6 m) portion of the aquifer, the relative aquifer fraction with As >10 µgL −1 dropped from 40% to 28% during the 115 simulated days, while the average As concentration dropped from 31 µgL −1 to 20 µgL −1 . The PNE computed for different thresholds also changed when shallower parts of the aquifer were sampled compared to the full-aquifer sampling. For As < 10 µgL −1 , the PNE increased from 58% to 70%, while for As < 50 µgL −1 the probability increased from 71% to 86%. Thus, the PNE was higher if groundwater was withdrawn at shallower depths. This is explained considering that the impact of rainfall on As attenuation is exacerbated in the shallower parts of the aquifer, since rainfall-controlled recharge has a more effective redox control in the shallower portions of the aquifer compared to the deeper parts.
We finally remark that the adopted modeling approach is general and can be applied to assess the risk of any toxic compound in complex systems. It can be easily extended to a different aquifer and/or to a different compound, provided that a mechanistic description of the processes controlling the mobility of the target variable are known.