University of Southern California Dissertations and Theses

Groundwater contaminant transport predictions in a sustainable water management scenario

(USC Thesis Other)

Groundwater contaminant transport predictions in a sustainable water management scenario

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Groundwater contaminant transport predictions in a sustainable water management scenario by Arianna Libera A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Environmental Engineering) August 423; Copyright 423; Arianna Libera Dedication This dissertation is dedicated to my dear grandmother, Chiara Cornaggia, who always profoundly loved me and fulfilled my vivid curiosity with all her wonderful stories. Questa tesi è dedicata alla mia cara nonna, Chiara Cornaggia, che mi ha sempre amato pro- fondamente e ha sempre soddisfatto la mia viva curiosità con tutte le sue fantastiche storie. ii Acknowledgements At the end of this journey that occupied the last five years of my life I would like to express my gratitude to a few people that played a fundamental role. First of all, I would like to thank my advisor, Professor Felipe de Barros, who believed in me from day one and never stopped. Thank you Felipe for everything you did for me, for all your help, for all your knowledge, for all your positivity. Thank you for always cheering and supporting my work. Thank you for helping me grow from both a professional and a personal point of view and always reminding me to enjoy the journey, while we strive to reach our ambitious goals. I profoundly admire you and I will always consider myself extremely lucky to have shared this journey with you. I will never forget these beautiful years, and you played a major role in making them so beautiful, so thank you, from the bottom of my heart! I would also like to thank Felipe for introducing me to wonderful collaborators. Between those, IamparticularlyfondtoHarukuWainwrightwhohostedmeattheLawrenceBerkeleyLaboratory during summer 4239. Thank you for everything you thought me Haruku and for always being so nice and caring. I would also like to thank Haruku for all the connections she provided by introducing me to people in the field and inviting me to conferences. I would like to thank Professor JJ Lee and the financial support provided by the Foundation for Cross Connection Control and Hydraulic Research at the University of Southern California for supporting my work in the first years. I would also like to thank Professor Kelly Sanders for always being an inspiring figure and for being part of all the important milestones of my PhD. iii I also express my gratitude to Professor Alberto Guadagnini, who connected me with Felipe and who provided wonderful support during my first two publications. A special thank you goes to my research group members, old and new. I was lucky to be part of such a beautiful group and I am particularly fond to Maria, Alessandra and Gerry for also being great friends who were there during some though moments and are part of my personal life. Finally, I would like to thank my parents, which have given me wings to fly, have supported me through everything and have always been on my side since the day I was born. Thank you mom and dad, I couldn’t have better parents and I would have never come so far without you, you are my true strength! Thank you for accepting and cheering my though choices and for always believing in me. A very special thank you goes to my little sisters, Ileana and Letizia, also always on my side, cheering and admiring me, and to my ants, uncles, cousins, I am so lucky to have you all! Infine, vorrei ringraziare i miei genitori, che mi hanno dato ali per volare, mi hanno supportato in tutto e son sempre stati dalla mia parte dal giorno in cui sono nata. Grazie mamma e papà, non avrei potuto avere genitori migliori e non sarei mai arrivata così lontano senza voi, siete la mia vera forza! Grazie per aver accettato e supportato le mie scelte difficili e per aver sempre creduto in me. Un grazie speciale va alle mie sorelline, Ileana e Letizia, che son sempre state dalla mia parte e mi han sempre incoraggiato e ammirato, e alle mie zie, i miei zii e tutti i miei cugini, sono cosi fortunata ad avere tutti voi! iv Table of Contents Dedication ii Acknowledgements iii List Of Tables viii List Of Figures ix Chapter 3 Motivation 3 3.3 Challenges in modeling groundwater contaminant transport in a sustainable water management scenario . . . . . . . . . . . . . . . . 3 3.4 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 3.5 Brief overview of the current state-of-the-art . . . . . . . . . . . . . . . . . . . . . . 8 3.6 Outline of the document . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . : Chapter 4 Hydraulic conductivity and porosity heterogeneity controls on environmental performance metrics: Implications in probabilistic risk analysis 33 4.3 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.4 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.4.3 K−φ model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.4.4 Flow and transport model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3: 4.4.5 Health risk model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.5 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.6 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4: 4.6.3 Analysis for non-reactive contaminants . . . . . . . . . . . . . . . . . . . . . 4: 4.6.3.3 CDF of first arrival times and peak mass fluxes . . . . . . . . . . . 4: 4.6.3.4 Relative plume dispersion . . . . . . . . . . . . . . . . . . . . . . . 54 4.6.3.5 Impact of connectivity on risk . . . . . . . . . . . . . . . . . . . . 56 4.6.4 Analysis for reactive contaminants: implications in increased lifetime cancer risk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.6.4.3 Low order moments of R T (x) . . . . . . . . . . . . . . . . . . . . . 59 4.6.4.4 PDF of R T . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Chapter 5 Influence of pumping operational schedule on solute concentra- tions at a well in randomly heterogeneous aquifers 68 5.3 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.4 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.5 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 v 5.6 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.6.3 Concentration breakthrough curves observed at the pumping well for differ- ent pumping strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7: 5.6.4 Mean and variance of solute concentration at the pumping well . . . . . . . 7; 5.6.5 Probability of exceedance of concentration limit value . . . . . . . . . . . . 85 5.6.6 Local scale dispersion and hydrogeological heterogeneity . . . . . . . . . . . 88 5.6.6.3 Impact of local scale dispersion on the concentration statistics . . 88 5.6.6.4 Impact of aquifer heterogeneity on the concentration statistics . . 8: 5.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Chapter 6 Soluteconcentrationatawellinnon-Gaussianaquifersundercon- stant and time-varying pumping schedule 97 6.3 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.4 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9; 6.5 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . :3 6.5.3 Domain configuration and numerical implementation . . . . . . . . . . . . . :3 6.5.4 Random Y Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . :5 6.6 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . :7 6.6.3 Temporal evolution of low-order moments of solute concentration at the well :8 6.6.4 Statistical analysis of the peak concentration . . . . . . . . . . . . . . . . . :; 6.6.5 Outliers in the peak concentration distribution . . . . . . . . . . . . . . . . :; 6.6.6 Probability density function of the maximum concentration . . . . . . . . . ;3 6.6.7 Average of the maximum concentration . . . . . . . . . . . . . . . . . . . . ;5 6.7 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ;6 Chapter 7 Climate change impact on residual contaminants under sustainable remediation ;: 7.3 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ;: 7.4 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322 7.4.3 Conceptual model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 322 7.4.4 F-Area site description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324 7.5 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325 7.5.3 Flow and transport simulations . . . . . . . . . . . . . . . . . . . . . . . . . 325 7.5.4 Modeling scenarios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 328 7.6 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 332 7.6.3 Spatiotemporal Dynamics of the Contaminant Plume . . . . . . . . . . . . 332 7.6.4 Contaminant concentration and export BTCs . . . . . . . . . . . . . . . . . 334 7.6.5 Impact of recharge perturbations on key Environmental Performance Metrics342 7.6.5.3 Maximum difference between baseline and perturbed BTCs . . . . 342 7.6.5.4 Difference of time above MCL between baseline and perturbed BTCs346 7.7 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 348 Chapter 8 Quantification of solute arrival time uncertainty in the planning of Managed Aquifer Recharge 355 8.3 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 355 8.4 Model formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 358 8.5 Numerical implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35: 8.6 Computational results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363 8.6.3 Residence time reliability . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363 8.6.4 Impact of areal recharge and Y variance on arrival times uncertainty . . . . 364 8.6.5 Impact of well pumping rate on arrival time uncertainty . . . . . . . . . . . 366 vi 8.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 366 Chapter 9 Conclusions 369 Reference List 374 Appendix A Low-order moments of the flux-averaged concentration . . . . . . . . . . . . . . . . . . . 389 Appendix B Supplementary material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38; vii List Of Tables 4.3 Risk model parameters used in the numerical simulations. . . . . . . . . . . . . . . 45 4.4 Parameter set used for the numerical flow and transport simulations. . . . . . . . . 48 4.5 Reaction parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.6 Mass transfer parameters of the source zone selected for the reactive transport simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4: 5.3 Parameter set employed in the numerical simulations. . . . . . . . . . . . . . . . . 77 6.3 Main parameters employed in the study. . . . . . . . . . . . . . . . . . . . . . . . . :3 6.4 Parameters characterizing the transmissivity field. . . . . . . . . . . . . . . . . . . :7 7.3 Parameters adopted in the numerical model presented in this study. The table reports values for porosity (φ [−]), permeability (k [m 2 ]), water retention-curve parameters (α [kg −1 m s 2 ],n [−],m [−]), residual liquid pore saturation (S rl [−]), Mualem (1976) [364] parameter (p [−]) (for details, see [34]) . . . . . . . . . . . . . 325 8.3 Model discretization and physical parameters employed in the study. . . . . . . . . 35; viii List Of Figures 3.3 Schematic representation of the main factors influencing groundwater contaminant transport: climatic conditions, natural recharge regime, pumping schedule, popu- lation water demand and subsurface heterogeneity. R denotes aquifer recharge rate and Q w corresponds to the pumping rate of groundwater extraction. . . . . . . . . 5 3.4 The main factors impacting subsurface contaminant transport are enclosed in three categories: climatic variability, engineering variability and natural variability of porous formations. These three areas are extremely interconnected. . . . . . . . . . 6 4.3 Log-logplotofpermeability(k)asafunctionofporosity(φ)fordifferentformations showing a positive correlation trend: Fontainebleau sandstone data from the study of [75] in purple, Upper Jurassic Brae Formation, East Brae field Offshore United Kingdom, North Sea [:7,365] in blue light, Permian-Triassic sandstones, Ivishak Formation, Sadlerochit Group, Prudhoe Bay Field, Alaska [:,365] in red, core-plug samples of Cenozoic sediments from geothermal well AST-02 in the Roer Valley Graben, southern Netherlands [353] in green. . . . . . . . . . . . . . . . . . . . . . 39 4.4 Schematic illustration of the problem under investigation: the contaminant is re- leased within the areal source zone,A s . A uniform-in-the average natural base flow, q 0 , takes place along the x−direction, contaminant mass fluxes or concentrations are measured at control planes (CPs). . . . . . . . . . . . . . . . . . . . . . . . . . 3; 4.5 Histogram of porosity values obtained by applying the KC equation (3) to one random realization of the conductivity field characterized by σ 2 Y = 1. The average porosity is 23.9%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.6 Cumulative distribution function (CDF) of the first arrival times, t 5% , of the con- taminantmassattheCPsforσ 2 Y = 1(a)andσ 2 Y = 3(b). TheCDFofdimensionless t 5% is indicated with dashed lines for heterogeneous φ fields and with solid lines for homogeneous φ fields. The green color indicates results at the CP located at dimensionless distance ζ = 4.29 whereas the red color pictures results at the CP located at dimensionless distance ζ = 15.71. . . . . . . . . . . . . . . . . . . . . . . 52 ix 4.7 Cumulative distribution function (CDF) of the contaminant mass peak, ˙ m p , at the CPs for σ 2 Y = 1 (a) and σ 2 Y = 3 (b). The CDF of ˙ m p is indicated with dashed lines for heterogeneous φ fields and with solid lines for homogeneous φ fields. The green color indicates results at the CP located at dimensionless distance ζ = 4.29 whereas the red color pictures results at the CP located at dimensionless distance ζ = 15.71. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.8 Box plots of Δτ (4.37) for homogeneous and heterogeneous φ and for σ 2 Y = 1 (a) and σ 2 Y = 3 (b). The green and red colors respectively indicate the results at the CP located at dimensionless distance ζ = 4.29 and ζ = 15.71. . . . . . . . . . . . . 55 4.9 Peak mass flux, ˙ m p , versus connectivity metric,CI (4.38), forσ 2 Y = 1 (a) andσ 2 Y = 3 (b). The red color indicates results at the CP located at dimensionless distance ζ = 1 whereas the green color pictures results at the CP located at dimensionless distance ζ = 4.29. Full circles indicate homogeneous φ and empty circles indicate heterogeneous φ. The regression line is plotted in gray for homogeneous φ and in black for heterogeneousφ for both CPs locations. The coefficient of determination (R 2 ) of the regression line for σ 2 Y = 1 at the CP located at dimensionless distance ζ = 1 are R 2 = 0.778 for homogeneous φ, R 2 = 0.599 for heterogeneous φ and R 2 = 0.656 for homogeneous φ, R 2 = 0.444 for heterogeneous φ at the CP located at dimensionless distance ζ = 4.29. The coefficient of determination (R 2 ) of the regression line for σ 2 Y = 3 at the CP located at dimensionless distance ζ = 1 are R 2 = 0.871 for homogeneousφ,R 2 = 0.816 for heterogeneousφ andR 2 = 0.845 for homogeneousφ,R 2 = 0.728 for heterogeneousφ at the CP located at dimensionless distance ζ = 4.29. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.: Longitudinal spatial distribution of the mean total increased lifetime cancer risk, hR T i, when the maximum concentration (a) or the maximum running averaged concentration over the exposure duration (ED) time (b) are considered to compute R T . The results forσ 2 Y = 1 are represented in magenta and the results forσ 2 Y = 3 are pictured in blue. Dashed curves correspond to heterogeneousφ fields and solid curves picture homogeneous φ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5: 4.; Longitudinal spatial distribution of the coefficient of variation of the total increased lifetime cancer risk,CV R T , when the maximum concentration (a) or the maximum running averaged concentration over the exposure duration (ED) time (b) are considered to compute R T . The results for σ 2 Y = 1 are represented in magenta and the results for σ 2 Y = 3 are pictured in blue. Dashed curves correspond to heterogeneous φ fields and solid curves picture homogeneous φ. . . . . . . . . . . 62 4.32 Probability distribution function (PDF) of the total increased lifetime cancer risk, R T , at the CP located atζ = 10 (a) andζ = 15.71 (b). The results forσ 2 Y = 1 are represented in magenta and the results for σ 2 Y = 3 are pictured in blue. Dashed curves correspond to heterogeneous φ fields and solid curves picture homogeneous φ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 x 5.3 Schematic representation of the problem studied. A two-dimensional flow field characterized by the presence of a pumping well and a uniform (in the mean) base flow from left to right of the domain is considered. The contaminant source is rectangular with longitudinal dimension L 1 and transverse dimension L 2 , located at distance L from the well positioned at x w = (x w ,y w ). The contaminant source is aligned with respect to the pumping well center. . . . . . . . . . . . . . . . . . . 73 5.4 Pumping scenarios analyzed in the study: Scenario (a) I; (b) II; (c) III. . . . . . 78 5.5 Pumping scenarios analyzed in the study: Scenario (a) I; (b) II; (c) III. . . . . . 82 5.6 Mean concentrationhCi observed at the pumping well for the uniform pumping strategy of ScenarioIII (a) and the transient pumping strategy of ScenarioIII (b). 84 5.7 Variance of the concentration Var[C] observed at the pumping well for the uni- form pumping strategy of Scenario III (a) and the transient pumping strategy of Scenario III (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 5.8 Coefficient of variation CV of the concentration observed at the pumping well for the uniform pumping strategy of Scenario III. . . . . . . . . . . . . . . . . . . . . 86 5.9 Coefficient of variation CV of the concentration observed at the pumping well for the transient pumping strategy of Scenario III. . . . . . . . . . . . . . . . . . . . . 86 5.: ProbabilityofexceedanceoftheconcentrationthresholdP (C≥C ∗ )fortheuniform pumping strategy of Scenario III. . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 5.; Probability of exceedance of the concentration threshold P (C≥C ∗ ) for the tran- sient pumping strategy of Scenario III. . . . . . . . . . . . . . . . . . . . . . . . . 87 5.32 Impact of the Péclet number (Pe) on the concentration mean for (a) constant and (b) time-varying pumping rate. Impact of the Péclet number (Pe) on the concentration variance for (c) constant and (d) time-varying pumping rate. Purple curves refer to Pe = 800 and blue curves refer to Pe = 200. . . . . . . . . . . . . . 8: 5.33 Mean concentration for (a) constant and (b) time-varying pumping rates. Variance behavior is also depicted for (c) constant and (d) time-varying pumping schemes. Blue curves refer to low heterogeneity and purple curves to high heterogeneity. . . 8; 5.34 Probability of exceedance of the concentration threshold P (C≥ C ∗ ) for (a) con- stant and (b) time-varying pumping rates. Blue curves refer to low heterogeneity and purple curves to high heterogeneity. . . . . . . . . . . . . . . . . . . . . . . . . 92 5.35 MeanconcentrationhCiobservedatthepumpingwellfor(a)constantand(b)time- varying pumping strategy of ScenarioIII for different levels of aquifer heterogeneity. 92 6.3 Sketch of the problem analyzed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . :2 6.4 Pumping strategies (a) constant flow rate S I and (b) variable flow rate S II . . . . . :4 xi 6.5 Temporal evolution of the normalized mean concentrationhCi observed at the pumpingwellforthreevaluesofα = 1.2, 1.5, 1.8andforaGaussianY field(α→ 2). Results for pumping strategy (a) S I and (b) S II . . . . . . . . . . . . . . . . . . . . :8 6.6 Temporal evolution of the normalized concentration variance Var[C] observed at the pumping well for three values of α = 1.2, 1, 5, 1.8 and for a Gaussian Y field (α→ 2). Results for pumping strategy (a) S I and (b) S II . . . . . . . . . . . . . . . :: 6.7 Pumping strategy S I : Box plots of C p for three values of α = 1.2, 1, 5, 1.8 and for a Gaussian Y field (α→ 2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ;2 6.8 Sketch of the problem analyzed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ;2 6.9 Pumping strategy S I : Peak concentration PDF p(C p ) for three values of α = 1.2, 1, 5, 1.8 and for a Gaussian Y field (α→ 2). See inset for peak concentration CDF P (C p ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ;4 6.: Pumping strategy S II : Peak concentration PDF p(C p ) for three values of α = 1.2, 1, 5, 1.8 and for a Gaussian Y field (α→ 2). See inset for peak concentration CDF P (C p ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ;4 6.; Relative impact of non-Gaussianity on the mean peak concentration C p measured through η, see equation 6.8, for three values of α = 1.2, 1, 5, 1.8 and pumping strategy S I (light grey) and S II (dark grey). . . . . . . . . . . . . . . . . . . . . . ;6 7.3 Schematic illustration of the hydrological conceptual model under investigation, representing a vertical two-dimensional cross-section driven along the middle line of the contaminant source zone. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323 7.4 Two-dimensional model domain adopted in our study. . . . . . . . . . . . . . . . . 326 7.5 Schematic illustration of the perturbed recharge scenarios analyzed in the study: (1) the constant positive recharge shift starting in 2020 (green line), (2) the con- stant negative recharge shift starting in 2020 (blue line), (3) the one-year extreme recharge in 2020 (magenta line) and (4) the hypothetical cap failure and constant positive recharge shift starting in 2020 (green line). The baseline recharge scenario is denoted by the horizontal black line. . . . . . . . . . . . . . . . . . . . . . . . . . 329 7.6 Schematic illustration of the impact of increased recharge on the concentration breakthrough curve (BTC) at an observation well located downgradient from the source zone. The recharge rate increase causes a slight dilution (Phase I), followed by a “rebound” effect in the BTC (Phase II) due to contaminants’ mobilization. The continuous blu curve indicates the BTC under baseline recharge whereas the dashed red curve corresponds to the concentration BTC under a change in the recharge conditions (see inset figure). . . . . . . . . . . . . . . . . . . . . . . . . . . 332 7.7 Illustration of the simulated tritium plume in the SRS. The initial condition is displayed in panel (a). The plume snapshots are shown for years (b) 1961, (b) 1988, (c) 2000 and (d) 2033. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333 xii 7.8 Temporal evolution of 3 H concentration at the: (a) source-zone well, (b) down- gradient well and (c) 3 H export at the CP for the baseline scenario, indicated by = 0, and the constant positive recharge shift scenario, characterized by = 0.1, 0.2, 0.3, 0.4, 0.5. The results of different recharge perturbations are por- trayed in different colors. The thin horizontal dashed black line represents the MCL for 3 H. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 335 7.9 Temporal evolution of 3 H concentration at the: (a) source-zone well, (b) down- gradient well and (c) 3 H export at the CP for the baseline scenario, indicated by = 0, and the constant negative recharge shift scenario, characterized by = 0.1, 0.2, 0.3, 0.4, 0.5. The results of different recharge perturbations are por- trayed in different colors. The thin horizontal dashed black line represents the MCL for 3 H. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 336 7.: Temporal evolution of 3 H concentration at the: (a) source-zone well, (b) downgra- dient well and (c) 3 H export at the CP for the baseline scenario, indicated by = 0, and the one-year extreme recharge scenario, characterized by = 0.5, 1, 4, 9. The results of different recharge perturbations are portrayed in different colors. The thin horizontal dashed black line represents the MCL for 3 H. . . . . . . . . . . . . 337 7.; Temporal evolution of 3 H concentration at the: (a) source-zone well, (b) down- gradient well and (c) 3 H export at the CP for the baseline scenario, indicated by = 0, the cap failure under the baseline scenario, indicated by = 0 (c) , and the cap failure under the constant positive recharge shift scenario, indicated by = 0.1 (c) , 0.2 (c) , 0.3 (c) , 0.4 (c) , 0.5 (c) , with (c) indicating the capping failure. The results of different recharge perturbations are portrayed in different colors. The thin horizontal dashed black line represents the MCL for 3 H. The inset shows the log-scale plot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 338 7.32 Maximum difference between contaminant concentrations and export BTCs of the baseline scenario and the perturbed recharge scenarios (γ c ) versus perturbation () for the constant positive recharge shift scenario (a), the constant negative recharge shift scenario (b), the one-year extreme recharge scenario (c), the cap failure and constant positive recharge shift scenario (d), where the superscript (c) indicates the capping failure. Results at the source-zone well are illustrated in red, at the downgradient well are pictured in blue and at the CP are indicated in green. . . . 344 7.33 Difference of time of exceedance of the MCL for 3 H between the baseline scenario and the perturbed recharge scenarios (γ t ) versus perturbation () for the constant positive recharge shift scenario (a), the constant negative recharge shift scenario (b), theone-yearextremerechargescenario(c), thecapfailureandconstantpositive recharge shift scenario (d), where the superscript (c) indicates the capping failure. Resultsatthesource-zonewellareillustratedinredandresultsatthedowngradient well are pictured in blue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 346 7.34 Temporal evolution of 3 H concentration at the: (a) source-zone well, (b) down- gradient well and (c) 3 H export at the CP for the baseline scenario, indicated by = 0, and the positive three months recharge shift scenario, characterized by = 0.1, 0.2, 0.3, 0.4, 0.5. The results of different recharge perturbations are por- trayed in different colors. The thin horizontal dashed black line represents the MCL for 3 H. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352 xiii 8.3 Conceptual model of the problem under investigation: a three-dimensional uncon- fined aquifer subject to a distributed areal recharge and managed aquifer recharge (MAR) through a perculating basin. Five extraction wells are operating. Emerg- ing contaminants move along the flow direction and reach a control plane located downgradient. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 359 8.4 Mean temporal evolution of the fraction of particles recovered at the control plane over 500 MC simulations for time step (δt) of 1 day and 10 5 inserted particles (blue curve) and for time step of 100 days and 10 4 inserted particles (orange circles). For this comparison we considered K fields characterized by σ 2 Y = 1, recharge R, and Q w =1000 m 3 d −1 applied to the five extraction wells. . . . . . . . . . . . . . . . . 362 8.5 Cumulative distribution function (CDF) of dimensionlesst 2% forσ 2 Y = 1, recharge R, and Q w = 1000 m 3 d −1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364 8.6 Probability density function (PDF) of t 5% (a), t 50% (b), t 95% (c) for higher Y heterogeneity (σ 2 Y = 1, continuous curves) and lower Y heterogeneity (σ 2 Y = 0.25, dashed curves) and recharge conditions R. PDF of t 5% (d), t 50% (e), t 95% (f) for higher Y heterogeneity (σ 2 Y = 1, continuous curves) and lower Y heterogeneity (σ 2 Y = 0.25, dashed curves) and lower recharge conditions. . . . . . . . . . . . . . . 365 8.7 Probability density function (PDF) of t 5% , t 50% and t 95% for Q w =500 m 3 d −1 (a) and Q w =1000 m 3 d −1 (b) and σ 2 Y = 1. . . . . . . . . . . . . . . . . . . . . . . . . . 366 A.3 Schematic representation of the problem analyzed using the Lagrangian framework. 389 SM3 Mean concentrationhCi observed at the pumping well for the uniform pumping strategy of Scenario I. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38; SM4 Variance of the concentrationVar[C] observed at the pumping well for the uniform pumping strategy of Scenario I. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 392 SM5 Mean concentrationhCi observed at the pumping well for the transient pumping strategy of Scenario I. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 392 SM6 VarianceoftheconcentrationVar[C]observedatthepumpingwellforthetransient pumping strategy of Scenario I. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393 SM7 Mean concentrationhCi observed at the pumping well for the uniform pumping strategy of Scenario II. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393 SM8 Variance of the concentrationVar[C] observed at the pumping well for the uniform pumping strategy of Scenario II. . . . . . . . . . . . . . . . . . . . . . . . . . . . 394 SM9 Mean concentrationhCi observed at the pumping well for the transient pumping strategy of Scenario II. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394 SM: VarianceoftheconcentrationVar[C]observedatthepumpingwellforthetransient pumping strategy of Scenario II. . . . . . . . . . . . . . . . . . . . . . . . . . . . 395 xiv SM; Coefficient of variation CV of the concentration observed at the pumping well for the uniform pumping strategy of Scenario I. . . . . . . . . . . . . . . . . . . . . . 395 SM32Coefficient of variation CV of the concentration observed at the pumping well for the transient pumping strategy of Scenario I. . . . . . . . . . . . . . . . . . . . . 396 SM33Coefficient of variation CV of the concentration observed at the pumping well for the uniform pumping strategy of Scenario II. . . . . . . . . . . . . . . . . . . . . 396 SM34Coefficient of variation CV of the concentration observed at the pumping well for the transient pumping strategy of Scenario II. . . . . . . . . . . . . . . . . . . . . 397 SM35ProbabilityofexceedanceoftheconcentrationthresholdP (C >C ∗ )fortheuniform pumping strategy of Scenario I. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 397 SM36Probability of exceedance of the concentration threshold P (C >C ∗ ) for the tran- sient pumping strategy of Scenario I. . . . . . . . . . . . . . . . . . . . . . . . . . 398 SM37ProbabilityofexceedanceoftheconcentrationthresholdP (C >C ∗ )fortheuniform pumping strategy of Scenario II. . . . . . . . . . . . . . . . . . . . . . . . . . . . 398 SM38Probability of exceedance of the concentration threshold P (C >C ∗ ) for the tran- sient pumping strategy of Scenario II. . . . . . . . . . . . . . . . . . . . . . . . . 399 SM39Transmissivity field T 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39: SM3:Transmissivity field T 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39; SM3;Transmissivity field T 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3:2 SM42 Mean concentrationhCi observed at the pumping well for the constant pumping strategy of Scenario III and a homogeneous aquifer. . . . . . . . . . . . . . . . . 3:3 SM43Mean concentrationhCi observed at the pumping well for the transient pumping strategy of Scenario III and a homogeneous aquifer. . . . . . . . . . . . . . . . . 3:3 xv Chapter 3 Motivation 3.3 Challenges in modeling groundwater contaminant transport in a sustainable water management scenario Predicting fate and transport of contaminants in groundwater systems represents a non-trivial task. Aquifers are heterogeneous by nature and capturing the variability of key hydrogeological properties, such as hydraulic conductivity (K) and porosity (φ) is unfeasible due to the high costs of data acquisition and the restricted availability of financial resources. Given limited site characterization data sets, model predictions regarding contaminant transport are subject to un- certainty. To overcome these difficulties, geostatistical models are combined with phisically driven probabilistic mathematical frameworks to forecast flow and contaminant transport as well as cor- responding uncertainty bounds. Geostatistics is employed to estimate the spatial distribution of hydrogeological properties, whereas probabilistic methods enable to estimate the uncertainty in transport model predictions associated with the lack of information regarding the structure of the geological properties. On the other hand, groundwater contamination is influenced by the different stresses associated with the water cycle and the societal demand. Amongst these we highlight external factors such as precipitation, evapotranspiration and natural recharge. All of these external factors are directly or indirectly linked to climatic variability. Simultaneously, 3 anthropogenic stresses are associated with engineered groundwater extraction or injection rates which depend on the available resources and on society needs. Therefore, reliable predictions of groundwater contamination rely on an accurate characterization of subsurface’s hydrogeological properties and on a realistic consideration of the external stresses to the aquifer system (e.g., precipitation, recharge, groundwater extraction, etc.). Figure 3.3 pictures an example of the connection between multiple elements influencing sub- surface contaminant transport. Factors connected to climatic conditions (e.g, precipitation, evap- otranspiration, natural recharge) and engineering stresses (e.g, groundwater pumping) are ex- tremely interdependent in the water management scenario. Climatic conditions dictate the time variability of precipitations, that in turn affect the natural infiltration rates to the subsurface en- vironment (Figure 3.3, segment 1), which impacts the availability of groundwater to satisfy society needs (Figure 3.3, segment 2). Groundwater extraction is, at the same time, scheduled to satisfy drinking, irrigation and industrial purposes, i.e., population demand (Figure 3.3, segment 3). The time variability of natural aquifers recharge and engineered extraction are therefore connected in a non-linear and complex manner (Figure 3.3, segment 2). The mentioned stresses (i.e., climatic and anthropogenic) to the subsurface system impact the transport of hazardous substances throughout aquifers (Figure 3.3, segment 4). Figure 3.3 shows a planar sketch of a contaminant plume moving through an heterogeneous aquifer, where red areas inthecontourplotindicatehighK channelswhereasbluezonesarecharacterizedbylowK values. All factors depicted in Figure 3.3 lead to the need to develop a holistic framework that captures the complex interaction of these factors on groundwater management. For example groundwater pumping enhances the mixing of the contaminated plume with clean water, therefore possibly leading to dilution at observation locations [344]. At the same time, well pumping operations can increase the likelihood that the contaminant plume will be captured by the well [344]. High precipitations increase the recharge rates and dilute contaminants but also lead to increased mobility of the contaminant plume. On the other hand, drought conditions decrease the plume 4 Figure 3.3: Schematic representation of the main factors influencing groundwater contaminant transport: climatic conditions, natural recharge regime, pumping schedule, population water demand and subsurface heterogeneity. R denotes aquifer recharge rate and Q w corresponds to the pumping rate of groundwater extraction. 5 mobility but can increase its persistence in the subsurface. Aside from climatic and anthropogenic stresses, solute transport phenomena are significantly affected by groundwater flow variability induced by spatial fluctuations of K. Variability in K leads to variability in the flow field which distorts the plume erratically as it migrates through the aquifers [63]. As a consequence, natural heterogeneity has a clear impact on solute dilution and mixing [6;,333,345]. All these elements (aquifer properties and external stresses) are keyed together and their simultaneous consideration leads to non-trivial trade-offs in subsurface solute transport studies. Figure 3.4: The main factors impacting subsurface contaminant transport are enclosed in three categories: climaticvariability,engineeringvariabilityandnaturalvariabilityofporousformations. These three areas are extremely interconnected. 6 In the overall, the main elements dictating subsurface contaminant transport in modern day society can be enclosed in three main categories, represented through three circles in Figure 3.4: climatic variability, engineering variability and natural variability of the subsurface formation. Climatic and engineering variability embed the “external factor” acting on the groundwater sys- tem. 3.4 Objectives Our research study focuses on better understanding the effect of these interconnected factors on subsurface contamination. It is indeed important to systematically quantify the spatio-temporal dynamic behaviors of the variables under investigation under different conditions, e.g., how peak concentration at an environmentally sensitive location changes with engineering variability (e.g., selection of the pumping scheme) or how solute travel time statistics from a contaminant source to different observation locations vary with the natural variability of the porous media. Through the use of stochastic computational modeling and fundamental flow and mass transport physics we evaluate contaminant transport and the uncertainty associated with risk predictions in multiple scenarios, characterized by different heterogeneity configurations of the porous media and/or by the presence of external factors stressing the system. We are interested in evaluating the trade- offs emerging from the simultaneous consideration of two or more of the three areas displayed in Figure 3.4, i.e., the red contoured area of Figure 3.4. These trade-offs can potentially provide important guidelines for water and risk managers to direct decisions and risk analysis. Some outcomes of our work help identifying conditions in which one or multiple factors prevail over the others in subsurface contamination predictions. For example, there can be subsurface contamination problems whose best management requires a proper and realistic consideration of external factors (e.g., variable pumping), but a detailed characterization of aquifer’s proper- ties, notably associated with significant financial efforts, is not required. On the other hand, 7 the situation could be opposite, leading to a justified decision of directing financial resources to characterize aquifer’s heterogeneity. Moreover, the outcomes of our work show that predictions of contaminant transport under realistic conditions (e.g., temporally variable pumping) should be incorporated in real life water management practices. For instance, knowing that the vari- ability of the groundwater extraction rate controls the temporal evolution of the contaminant concentration recovered at the operating well, could lead to more accurate predictions of water treatment costs and provide guidelines to schedule optimal pumping sequences associated with reduced water treatment and/or energy expenses (e.g., the expenses associated with an increase of the pumping rate could be justified by a decrease of water treatment costs). Therefore, in- corporating realistic subsurface contamination forecasts that simulate the presence of external factors acting upon the heterogeneous hydrogeological contaminated system into water manage- ment could lead to significant improvements. Moreover, considering the uncertainty connected with aquifers’ heterogeneous properties together with the impact of external factors on subsurface contaminant transport could help improve the design of engineered water management operations, e.g., managed artificial recharge (MAR) systems. Our research line is motivated by the research gaps presented in the upcoming chapter. 3.5 Brief overview of the current state-of-the-art The effect of spatial variability of aquifers’ hydraulic properties, such as K, on contaminant transport, has been broadly studied [55,:3,3:4]. However, the impact of including the variability of other hydrogeological properties, as porosity (φ), on contaminant transport and risk analysis has been only minimally tackled. There have been studies, e.g., [393,396], which showed that the description of transport processes highly profits from the incorporation of the spatial variability of theφ field. A deeper understanding of the impact of including the spatial variability ofφ, together with the heterogeneity of K, on risk analysis is important to achieve more realistic predictions. 8 Moreover, the impact of temporally variable “external factors” related to climate (e.g., precipi- tation, recharge, etc.) or to realistic engineering decisions (e.g., employment of transient pumping schedules), on subsurface contamination, has been only marginally studied within the hydroge- ology community. Indeed, groundwater pumping has been mainly investigated by considering constant extraction schedules [49,95,436,437]. In reality, water management agencies typically schedule groundwater extraction through a well-defined temporal sequence of pumping rates to maximize benefits to anthropogenic activities and minimize the environmental footprint of the withdrawal process. The temporal dynamics of well pumping operations induce changes in the subsurface flow field which in turn can affect spreading and mixing of contaminant plumes. In order to achieve truthful predictions of contaminant transport and direct risk analysis it is fun- damental to evaluate the joint effect of the temporal dynamics of an operating pumping well and the (typically uncertain) spatial variability of the hydrogeological properties of the subsurface reservoir within which a potentially harmful solute is residing and migrating. Furthermore, a common procedure adopted in the stochastic hydrogeology community is to model the spatial variability of hydrogeological properties through multivariate Gaussian geo- statistics [332,3:4]. In general, the log-conductivity field Y = ln(K) is assumed to be multi- Gaussian, thus we rely on two points statistics. Spatial variability of hydrogeological attributes is, however, known to be more complex than described by a Gaussian model [9:,:5,;6,324,355, 445,452]. It is therefore relevant to analyze the feedback between transient pumping operations and the non-Gaussian nature of Y in influencing contaminant transport. Such analysis would allow us to fundamentally understand the relative importance of these non-Gaussian features in the Y-field in transport when engineering factors are accounted for. We additionally concur that hydrogeological studies on the climate variability effects on contaminant transport and environmental remediation have been extremely limited. Although the effect of climate change has been investigated extensively from the water resource perspec- tive [:4,:8,35:,383], a limited number of studies have addressed water quality issues [439]. Most 9 literature focuses on surface-water systems [9;,3;2,434,435,447] because, as Green et al. [:8] noted, surface water is more visible and accessible. Current literature urgently lacks site-scale hy- drological studies that can guide environmental remediation under extreme climatic events within risk and performance assessments and engineering designs, as well as regulatory frameworks. Finally, only a limited number of studies, e.g., [48,427], have addressed the management of artificial recharge operations, which is employed in California to increase the sustainability of water resources. Given that in different cases reclaimed water is utilized to recharge the aquifers, a minimum residence time of water injected underground is established in order to ensure that natural degradation processes occur before groundwater extraction. At the current stage, this estimation is based on simple deterministic models and tracer experiments. Those models are do not take into account the uncertainty connected with the spatial variability of hydrogeological properties or do not properly account for the presence of external factors (e.g., varying precipitation, well pumping) therefore there is a need of modeling studies that quantify the uncertainty connected with reclaimed water arrival time to specific locations and estimate connected risks by providing probability bounds. 3.6 Outline of the document Our work is described in the following chapters: the upcoming five chapters correspond to work in progress or published, according to the gaps in the literature identified above, whereas the last chapter summarizes the importance of our work and illustrates the key conclusions. Chapter 4 shows the importance of considering the spatial variability of subsurface properties (i.e.,K andφ) in stochastic contaminant transport studies in the absence of external factors. This first chapter aims at deepening our understanding of the sole effect of aquifers’ heterogeneity on transport and decision making, i.e. investigates the fundamentals of flow variability on transport : phenomena. The chapter also provides a systematic analysis on the impact of the aquifer disorder level on risk analysis and corresponding uncertainties. Following chapters 5-6 consider “external factors” stressing the heterogeneous subsurface sys- tem. In particular, chapter 3 shows the importance of accounting for time-varying pumping schedules when analyzing contaminant breakthrough curves (BTCs) at environmentally sensitive locations, in the presence of heterogeneous multi-Gaussian hydraulic conductivity fields. Chapter 4 expands the results of chapter 3 by including non-Gaussian models to characterize the hetero- geneous subsurface system. Chapter 7 illustrates a real case study, applied at the Savannah River Site, South Carolina, USA. Chapter 7 investigates the effect of extreme events in a climate change setting on ground- water contamination and environmental remediation. Chapter 8 presents ongoing research on the importance of stochastic modeling to quantify the uncertainty connected with solutes’ arrival times within the planning of aquifers’ artificial recharge, which can lead to a better design of the operation. The conclusions of our work are finally presented in the last chapter. The full references to the published and ongoing research are listed below. • Chapter 4: Libera, A., Henri, C.V., de Barros, F.P.J. (2019). “Hydraulic conductivity and porosity heterogeneity controls on environmental performance metrics: Implications in probabilistic risk analysis.” Advances in Water Resources, 27 : 1− 12. • Chapter 5: Libera, A., de Barros, F.P.J., Guadagnini, A. (2017). “Influence of pumping operational schedule on solute concentrations at a well in randomly heterogeneous aquifers.” Journal of Hydrology, 546 : 490− 502. • Chapter 6: Libera A., de Barros, F.P.J., Riva, M., Guadagnini, A. (2017). “Solute con- centration at a well in non-Gaussian aquifers under constant and time-varying pumping schedule.” Journal of Contaminant Hydrology, 205 : 37− 46. ; • Chapter 7: Libera A., de Barros, F.P.J., Faybishenko, B., Eddy-Dilek, C., Denham, M., Lipnikov, K., Moulton, D., Maco, B., Wainwright, H. (2019). “Climate change impact on residual contaminants under sustainable remediation.” Journal of Contaminant Hydrology, 226 : 103518. • Chapter 8: Libera A., de Barros, F.P.J. (2019). “Quantification of solute arrival time uncertainty in the planning of Managed Aquifer Recharge.” On going research. 32 Chapter 4 Hydraulic conductivity and porosity heterogeneity controls on environmental performance metrics: Implications in probabilistic risk analysis The content of this chapter is published in Advances in Water Resources [346]. 4.3 Introduction It is well established that heterogeneities in natural porous formations largely control subsur- face groundwater flow and contaminant transport. These heterogeneities are mainly manifested through the hydraulic conductivity (K) and, to a lesser degree, the porosity (φ). The influ- ence of the heterogeneousK-field on flow and solute transport processes has been widely studied (see [3:4] and references therein), whereas less attention is dedicated to the effect of heteroge- neous φ. Hydraulic conductivity is commonly modeled as a random space function (RSF) in the stochastic hydrogeology community because of its erratic spatial variability and the large uncertainty associated with incomplete site characterization [53,:2]. φ variability is commonly regarded as a secondary factor when compared to K heterogeneity. Indeed, most studies in the field of stochastic hydrogeology assume that aquifers are characterized by spatially heterogeneous K and homogeneous φ. This is partially justified because K can vary in space by 3− 4 orders of 33 magnitude within small distances while the range of variability of φ in unconsolidated granular aquifers is generally between 0.1 and 0.55 [9,97,323]. Because of the limited availability of K and φ measurements, it is difficult to establish a clear correlation between the two geological properties. However, in specific formations, a level of correlation is likely to be present [;8]. The literature presents different studies exploring the relationship between K and φ. [6] and [75] found positive correlation between the variables. Additionally, different authors investigating the estimation of K from φ and other measurable parameters (e.g. grain size, pore surface area, pore dimension) established positive correlation between K and φ for different types of soil (e.g., [3,7;,96,366,376,39;,3;4]). Among them, [366] analyzed different models to predict permeability in sedimentary rocks and concluded that, in most cases, the permeability is related to a power of φ and to the square of a measure of surface area or a characteristic length. It is typically assumed that K and φ are positively correlated in unconsolidated aquifers [323], however some studies observed negative correlation between the two properties under specific packing arrangements and grain size distributions [362]. Onlyasmallnumberofstudiesinvestigatedtheeffectsofspatiallyvariableφfieldsontransport predictions. Among them, [443] did not consider a correlation between K and φ and concluded that the impact of φ variations on macroscopic dispersion is minor with respect to the effect of K variability. [347] analyzed various porosity functions to represent φ spatial variability in porous media and indicated that the solute concentration predictions are within 15% of the concentration resulting from an averageφ when a linear or slightly curved equation that expresses porosity as a function of the spatial coordinate, is employed. The author pointed out that the discrepancies of solute concentrations become very substantial as the deviation of the porosity functionfromlinearityincreases. Afewnumericalmodelingstudies(e.g.,[;7,;8,323])showedthat the consideration of correlated heterogeneous K and φ fields significantly impacts contaminant transport for two-dimensional domains. [;8] indicated that a positive correlation betweenφ andK decreases dispersion whereas a negative correlation enhances plume spreading, particularly along 34 the longitudinal direction. [;8] also considered reactive transport and indicated that reaction rates impact the relative importance of φ variability on concentration estimation, which appears to be amplified in the presence of slow reactions. Along the same line of work, [;7] highlighted that φ spatial variability, correlated to heterogeneousK, impacts the longitudinal velocity covariance as well as its cross covariance with the head and the log-conductivity (Y ≡ ln[K]) fields. [396] studied theinfluenceofK heterogeneitiesonsolutetransportwithconstanteffectiveφ, andlaterexpanded theiranalysisin[393], wherebothK andφwereconsideredrandomvariables.[393]highlightedthat considering the spatial distribution ofφ is fundamental to mimic the early arrival of breakthrough curves (BTCs). The studies of [393,396] showed that the site description of transport processes highly profits from the incorporation of the spatial variability of the φ field. In light of the mentioned findings, the influence of heterogeneous K coupled with spatially variableφ fields, especially when the two properties are correlated, as shown in most cases, should be further investigated to increase the reliability of contaminant transport and associated hu- man health risk predictions. In this context, our work aims to study the significance of spatial variability in the porous medium’s φ to conservative and reactive solute transport predictions in the context of probabilistic risk analysis. Through the use of stochastic numerical simulations we systematically illustrate the impact of K-φ spatial variability on key Environmental Perfor- mance Metrics (EPMs) such as solute arrival times, peak mass fluxes, increased lifetime cancer risk and corresponding uncertainties at environmentally sensitive locations. For our illustrations, we consider three-dimensional (3D) heterogeneous aquifer systems. In particular, the objectives of this work are twofold. First, we systematically investigate the coupled influence of φ and K heterogeneity on the statistical characterization of non-reactive transport. Second, within an application-oriented context, we show the implications of these coupled effects when estimating adverse human health effects (e.g. cancer risk) due to exposure to chlorinated solvents in ground- water. Through extensive stochastic computational analysis, we explore different scenarios that, to the best of our knowledge, have not been investigated, with the goal of identifying conditions 35 where the role of spatially variable φ, correlated to heterogeneous K, is relevant. Our analy- sis simulates contaminant transport through an unconsolidated sandy aquifer where, given the likelihood [;7,323,353, ], positive K-φ correlation is adopted, together with different levels of heterogeneity of theK system. Our study addresses fundamental questions such as: is the uncer- tainty associated with transport in the subsurface increased or decreased when variable φ fields are considered? Does the relative importance of variable φ depend on the level of heterogeneity of the K-field? What is the impact of this joint heterogeneity on the flow connectivity structure of the aquifer? How are the effects of heterogeneity propagated to decision making metrics such as human health risk? The paper is structured as follows. We present the problem under investigation in Section 4.4. In Section 4.5, we illustrate the methodology employed in our study and provide details on the computational implementation. Next, Section 4.6 is devoted to the analysis of the results for conservative tracers (Section 4.6.3) and reactive contaminants (Section 4.6.4). We finally summarize our study and outline the key findings in Section 4.7. 4.4 Problem Statement Inthisworkweconsideracontaminated,fullysaturated 3DheterogeneousaquiferwhoseCartesian coordinate system is indicated by x = (x,y,z). We are interested in understanding the impact of heterogeneity, stemming from both K and φ, on a series of EPMs [5:]. Let Ω denote a generic EPM. For example, Ω can represent solute arrival times or solute peak concentrations at an environmentally sensitive location. If the contaminant is toxic, then Ω can denote, for instance, an adverse health effect. Due to incomplete characterization of the subsurface, hydrogeological properties such as K and φ are modeled as spatially random. As a consequence, the quantity Ω is subject to uncertainty and quantified through its statistical moments, such as its meanhΩi (where the angled brackets denote ensemble expectation), its 36 variance σ 2 Ω , and its probability density function (PDF) p(Ω) or cumulative density function (CDF) P (Ω). WeareprimarilyinterestedinanalyzingtheeffectofK-φvariabilityondifferentEPMswiththe goal of: (a) improving our fundamental understanding of the significance of these variabilities to contaminant transport and (b) exploring how the latter variabilities propagate to decision making metrics relevant to human health risk analysis. In this study we perform both conservative and reactive contaminant transport simulations and select the specific EPMs indicated below: • We statistically analyze the first arrival times and peak mass fluxes through control planes (CPs) resulting from conservative tracer simulations. The first arrival time is defined as the time of arrival of 5% of the initial injected mass (i.e., t 5% ) at the CP of interest. An accurate detection oft 5% is of fundamental importance for risk and water entities in charge of accurately planning aquifer remediation, risk analysis and managing water resources. Moreover, peak mass fluxes at an environmentally sensitive location can be viewed as an indicator of the risks associated with groundwater contamination and can also serve as a proxy for dilution [89]. • Wethenexploretheincreasedlifetimecancerriskestimatesassociatedwithchronicexposure to chlorinated solvents that are the reactive contaminants of interest in our study. The latter metric is of fundamental importance for public health authorities and its quantification provides useful insights for risk and remediation experts in order to characterize risk of real-life operations. 4.4.3 K−φ model The interdependence between the heterogeneousK andφ fields is largely uncertain as explained in Section 4.3. However a positive correlation is likely to be present in unconsolidated sandy aquifers (e.g., [323,353]). With the goal of systematically exploring the influence of spatially variable φ 37 andK on the EPMs described before, we model our aquifer system as being characterized by (a) spatially heterogeneous K and homogeneous φ, and (b) spatially heterogeneous K correlated to spatially heterogeneous φ. In order to illustrate the possible relationship between K and φ data, we report in Figure 4.3 permeability (k) and φ data for different formations. Data presented in purple refer to ex- perimental measurements of a sequence of seven samples of Fountainebleau sandstone [75], in blue light we show samples from a turbidite reservoir of Late Jurassic age located in the North Sea [:7,365]. Furthermore, data from Permian-Triassic sandstones from Prudhoe Bay Field, Alaska, USA, are pictured in red [:,365], whereas core-plug samples of Cenozoic sediments from the southern Netherlands [353] are shown in green in Figure 4.3. We recall that permeability,k, is directly proportional toK throughK =kρg/μ [35], whereρ indicates water density (1000 g/m 3 ), g is the gravity acceleration (9.8 m/s 2 ) and μ is the dynamic viscosity of water. We adopted a value of μ equal to 8.84× 10 −4 Pa s, i.e., a temperature of 25 ◦ C. Note that k in Figure 4.3 is expressed in md, 1 md being equivalent to 9.87× 10 −16 m 2 . Interpolation of φ and k data of Figure 4.3 suggests a positive correlation for all the formations presented. In our work we consider a porous formation composed by unconsolidated sandy material and we assume positive correlation between the heterogeneous K and φ fields. For the purpose of illustration we estimate that the relationship amongφ andK can be described by the well-known Kozeny-Carman (KC) empirical equation [3;,42,336]: K = ρg μ φ 3 (1−φ) 2 d 2 e 180 , (4.3) where d e indicates the representative diameter of grains and the remaining parameters have been previously defined. The literature presents different interpretations ford e (e.g., [335]) and we assume that the latter corresponds to d 10 , i.e., the grain diameter at which 32% of the particles are smaller, in our study. 38 10 20 30 40 φ [%] 1 10 100 1000 10000 k [md] Fontainebleau sandstone Upper Jurassic Brae Formation Permian-Triassic sandstones Cenozoic sediments Figure 4.3: Log-log plot of permeability (k) as a function of porosity (φ) for different formations showing a positive correlation trend: Fontainebleau sandstone data from the study of [75] in purple, Upper Jurassic Brae Formation, East Brae field Offshore United Kingdom, North Sea [:7, 365] in blue light, Permian-Triassic sandstones, Ivishak Formation, Sadlerochit Group, Prudhoe Bay Field, Alaska [:,365] in red, core-plug samples of Cenozoic sediments from geothermal well AST-02 in the Roer Valley Graben, southern Netherlands [353] in green. TheKCequationhasbeenproveneffectiveinpredictingthepermeabilityofsandfromporosity data and generally gives good predictions of hydraulic conductivity when the value of d 10 is within the range 0.1− 3 mm [;]. The KC formula has been employed in different studies (e.g., [;,39,353,372,395]). [353] concluded that permeability of pure sand can be predicted with high confidence (R 2 ≥ 0.9) using the KC equation. [39] used the relationship to first estimate the K values of lithofacies during three-dimensional stochastic flow and transport modeling at the MADE site. [372] tested the reliability of different empirical formulas to estimate the hydraulic conductivity of unconsolidated aquifers by particle size analysis and identified the KC equation as the overall best estimator. Other authors, as [33] and [43], also encouraged the use of the KC formula. We point out that other empirical relationships betweenφ andK, which are site-specific, can be utilized within the methodological framework adopted in this study. 39 4.4.4 Flow and transport model Groundwater flow A steady state uniform-in-the mean base flow (q 0 ) takes place along the longitudinalx-direction within the 3D heterogeneous aquifer in the absence of sinks and sources. Constant head conditions are established at the west and east boundaries and no-flow conditions are set at the remaining boundaries. The steady-state groundwater flow is governed by the following partial differential equation: ∇·q(x) = 0, (4.4) with q(x) = −K(x)∇h(x), (4.5) where q [m/d] is the specific discharge vector andh [m] is the hydraulic head. Note that (4.5) is written in terms of q, therefore the spatial variability of the φ field does not play a role in the solution of the groundwater flow problem. After solving (4.5), groundwater velocity is obtained from v(x) =−K(x)∇h(x)/φ(x), and used as input to the contaminant transport model. The contaminant transport equation is then including theφ field variability for both conservative and reactive chemical species. As previously mentioned, to improve our understanding of the joint K-φ heterogeneity on transport and its role in contaminated groundwater management, we simulate contaminant trans- port for non-reactive and reactive species. The latter case will be linked to a dose-response model to evaluate the adverse health effects associated with contaminant exposure. Reactive transport A contaminant is released along a vertical 2D plane of areal extent A s . Contaminant mass fluxes or concentrations are measured at multiple control planes (CPs) located 3: Control planes, CPs Areal source zone, A s q 0 x y z Figure 4.4: Schematic illustration of the problem under investigation: the contaminant is released within the areal source zone,A s . A uniform-in-the average natural base flow,q 0 , takes place along thex−direction, contaminant mass fluxes or concentrations are measured at control planes (CPs). at different longitudinal distances from the source zone. Figure 4.4 illustrates a sketch of the problem under investigation. Following the work of [;;], the aquifer is contaminated by tetra- chloroethylene (PCE), a DNAPL product often identified in groundwater [83]. We assume that PCE is entrapped in the source area and slowly dissolves in time, originating a long-term contam- ination that moves downgradient. In the presence of an anaerobic environment, PCE undergoes chemical degradation that generates trichloroethylene (TCE), a byproduct that in turn reacts and degrades into cis-Dichloroethylene (DCE). DCE degrates then into vinyl chloride (VC), which will finally produce ethene. The US Environmental Protection Agency (EPA) classifies PCE, TCE and DCE as probable human carcinogenic products, whereas categorizes VC as a human carcino- gen [429]. The reactive transport of PCE and its daughter products can be described by the following system of equations, where subscript i indicates a specific contaminant within the toxic PCE degradation chain (i.e., i = 1: PCE; i = 2: TCE; i = 3: DCE and i = 4: VC): φR i ∂C i ∂t −∇· (φD∇C i ) +∇· (qC i ) =y i k i−1 φC i−1 −k i φC i +s(x,t)δ i1 ,∀i = 1,..., 4, (4.6) 3; whereC i [g/m 3 ] is the resident contaminant concentration of thei th species, D [m 2 /d] represents the hydrodynamic dispersion tensor whose components along the x-, y- and z- direction are respectivelyD xx =α x v,D yy =α y v andD zz =α z v, withv [m/d] indicating groundwater velocity of the grid block. Here, α x , α y and α z are respectively the longitudinal (along x), transverse horizontal (along y) and transverse vertical (along z) dispersivity coefficients, assumed to be constant for all the contaminant species. Furthermore, R i represents the retardation factor, k i [d −1 ] indicates the first-order contaminant degradation rate andy i [g g −1 ] symbolizes the effective yield coefficient for any reactant or product pair. The latter represents the ratio between the generated mass of speciesi and the degraded mass of speciesi-1. We assume that biodegradation does not occur in the sorbed phase and that sorption reactions follow a linear sorption isotherm (see [;;] and references therein). δ i1 indicates the Kronecker delta function. We consider a time-dependent dissolution rate s(x,t) of PCE leaving the areal source, A s , starting from no contamination at initial time t = 0. Following the approach described in [;;], we employ the mass-depletion model proposed in the literature [379,387]: c s (t) c 0 = m(t) m 0 Γ , (4.7) wherec s is the flux-averaged concentration of the dissolved DNAPL species, i.e., PCE, leavingA s , c 0 is the initial contaminant concentration at the source zone,m is the mass of DNAPL remaining in the source zone and Γ represents the mass-depletion constant accounting for changes of the interfacial surface area while the mass at the source diminishes. The flux concentration of PCE leaving the areal source is then expressed as in [379]: c s (t) = c 0 m Γ 0 − Q s c 0 λ s m Γ 0 + m 0 1−Γ + Q s c 0 λ s m Γ 0 e (Γ−1)λst Γ 1−Γ , (4.8) wherem 0 is the initial contaminant mass at the source,Q s represents the groundwater volumetric discharge rate through the source zone while λ s is the first-order degradation constant of PCE 42 at the source [;;]. The function s(x,t), i.e., the time-dependent dissolution rate at the source, is given by: s(x,t) =q s c s (t)δ(x−x inj )Ω(x∈A s ), (4.9) where q s = Q s /A s , and Ω(x∈ A s ) is an indicator function equal to one if x∈ A s and equal to zero elsewhere. In Section 4.4.5, we present the dose-response model employed to compute the adverse health risk caused by the reactive chemicals degradation chain (i.e., PCE→ TCE→ DCE→ VC). Non-reactive transport With the goal of isolating the influence of heterogeneity in the aquifer geological propertiesK andφ on transport from the chemical reactions component, we also simu- late transport of a non-reactive tracer. In this case, similarly to the reactive transport simulations, groundwater flow is also governed by equations (4.4) and (4.5), whereas contaminant transport through the porous media is expressed by (4.6) with chemical reactions set to zero. During the conservative contaminant transport simulations the tracer is instantaneously and homogeneously released within A s and we investigate the first arrival times and peak mass fluxes through the control planes and corresponding uncertainties. 4.4.5 Health risk model Being the human health risk the EPM of concern for the reactive transport case, we follow the work of [356], [;;], and [454] and adopt a Poisson model [428] to evaluate the human incremental lifetime cancer risk (ILCR), R i (x) associated to an exposure to the contaminant i at the CP location x: R i (x) = 1− exp[−ADD i (x)×CPF i ], (4.:) 43 where CPF i [kg d/mg] indicates the metabolized cancer potency factor related to contaminant i and ADD i [mg/(kg d)] is the average daily dose. In this work, risk is due to human exposure by direct ingestion and: ADD i (x) = ¯ c i (x) IR BW ED×EF AT . (4.;) Inequation(4.;),IRisthewateringestionrate[L/d],BW isthebodyweight[kg],ED symbolizes the exposure duration [y],EF is the daily exposure frequency [d/y] andAT indicates the expected lifetime [d]. In our study, ¯ c i (x) [g/m 3 ] in (4.;) represents the flux-averaged concentration [337] at a CP situated at a given longitudinal distance x. This variable is identified both as: (a) the peak flux- averaged concentration and (b) the maximum exposure duration (ED) averaged (flux-averaged) concentration. The manner in which ¯ c i (x) is evaluated can lead to different groundwater man- agement strategies [66]. It is therefore important to evaluate both cases (a) and (b) above, the first case being the worst scenario that can provide a conservative risk estimate. In the second case (b) ¯ c i (x) is equal to the maximum running averaged concentration, over theED interval, of the contaminant concentration BTC at the CP [;;,356] and is mathematically given by: ¯ c i (x) = max t>0 ( 1 ED Z t+ED t c i (τ;x)dτ ) , (4.32) where c i (τ;x) is the flux-averaged concentration observed at the CP location. Given that we consider the degradation of PCE into sub-products, parent and daughter species are likely to be present simultaneously in the aquifer. The total risk of the chemical mixture (denoted byR T ), in the presence of low concentrations (less than 300 ppm), can then be identified as the summation of the individual risks caused by each chemical product [86] as follows: 44 Parameter Symbol Value Ingestion rate IR 1.4 l/d Body weight BW 70 kg Exposure duration ED 30 y Exposure frequency EF 350 d/y Average time of expected lifetime AT 25550 d PCE TCE DCE VC Cancer potency factor CPF i 0.0021 0.011 0.6 1.5 Table 4.3: Risk model parameters used in the numerical simulations. R T (x) = 4 X i=1 R i (x). (4.33) The health risk model parameters employed to compute the total ILCR, i.e.,R T (4.33), at the CPs, are listed in Table 4.3. 4.5 Methodology General set-up In order to analyze contaminant transport and corresponding uncertainties throughout the 3D aquifer we employ a stochastic Monte Carlo (MC) framework that accounts for uncertain K and φ distributions. Five hundred equally likely K and φ fields are generated and, for each field, the groundwater flow and contaminant transport equations are numerically solved. Results are then statistically analyzed over the ensemble of MC realizations. The 3D computational domain is characterized by dimensions L x ×L y ×L z = 700 m× 400 m× 100 m and discretized into 350× 200× 100 cells (i.e., cells dimensions of Δ x = Δ y = 2 m and Δ z = 1 m). A large number of particles (10 4 ) is uniformly injected within the source area A s = L sy ×L sz = 70 m× 20 m, located at x = x inj = 50 m and transversely centered within the aquifer domain. A set of 7 CPs are positioned at different distances along the x- direction (henceforth denoted by x CP ): 85 m, 100 m, 120 m, 200 m, 300 m, 400 m and 600 m. 45 The log-conductivity field Y ≡ ln[K] is spatially heterogeneous and represented by a Multi- Gaussian RSF, characterized by a Gaussian covariance function of ranges λ x , λ y and λ z respec- tively along thex-,y- andz-direction, meanhYi and varianceσ 2 Y . For the computational analysis, we set λ x = λ y = 35 m and λ z = 10 m. A sequential multi-Gaussian simulator [389] is employed to produce the ensemble of Y field realizations. Based on the data reported in Section 4.4.3, we assume that the heterogeneousK field is posi- tively correlated to the heterogeneousφ field and that the relationship between the two geological attributes can be described through the KC relationship (4.3). From each heterogeneous K field realization we then obtain the corresponding spatially variable φ field by applying the KC equa- tion (4.3). We analyze the influence of heterogeneity on transport in the presence of (a) constant φ and heterogeneous K and (b) spatially variable φ positively correlated to heterogeneous K. The outputs of scenarios (a) and (b) can be compared since the constant value of φ in case (a) coincides with the average value of the heterogeneous φ field realizations (i.e.,hφi). The average porosity,hφi, is roughly constant for all the MC realizations and can, therefore, be obtained by applying the geometric mean ofK (indicated asK G ) in the KC equation above (4.3). This leads to a meanφ value of 0.23. In order to give a better idea of the range ofφ values of our simulations, we present in Figure 4.5 the histogram ofφ for one random realization of theK field. The latter is characterized by σ 2 Y = 1 and the mean value of φ obtained is 0.239 (i.e., 23.9%). Furthermore, different levels of heterogeneity of the Y field are explored: we consider both a mildly heterogeneous aquifer, represented by σ 2 Y = 1 and a more heterogeneous system, charac- terized by σ 2 Y = 3 for both the inert tracer and the reactive transport simulations. In each case, the variance of the ln(φ) field, i.e., σ 2 ln(φ) , is one order of magnitude less that σ 2 Y , which is in agreement with the work of [;8]. Groundwater fluxes are computed by means of the widely tested numerical model (finite- difference) MODFLOW [;4]. Conservative and multispecies reactive transport are simulated through the RW3D code [85,;:] using the random walk particle tracking method. We select 46 Figure 4.5: Histogram of porosity values obtained by applying the KC equation (3) to one random realization of the conductivity field characterized by σ 2 Y = 1. The average porosity is 23.9%. appropriate dispersivity values in order to account for sub-grid heterogeneity. The main input parameters used in the numerical solution of the flow and transport models are reported in Table 4.4. Reactivetransport Asfarasthechemicalparameterschosenforthereactivetransportsimula- tions, we employ constant values of biodegradation rates, which we select accordingly to the range of observed values listed by the US EPA [42:]. Additionally, the choice of the retardation factors is based on the expected differences in mobility between the chlorinated solvents species [;;,34;]. Table 4.5 shows the chemical reaction parameters adopted in our study. During the reactive transport simulations particles are instantaneously and uniformly released within A s . The first arrival time of the particles at the CPs is recorded and cumulative concen- tration BTCs are computed to obtain c i h (t;x). The flux-averaged concentration of any reacting speciesi,c i (t;x), caused by the time-dependent injection described in (4.9), is computed through the principle of superposition as follows: c i (t;x) = t Z 0 c s (τ)c δ i (t−τ;x)dτ, (4.34) 47 Parameter Symbol Value Aquifer domain (L x , L y , L z ) (700 m, 400 m, 100 m) Grid discretization (Δ x , Δ y , Δ z ) (2 m, 2 m, 1 m) Y ranges (λ x , λ y , λ z ) (35 m, 35 m, 10 m) Variance of Y σ 2 Y 1, 3 Average porosity hφi 0.23 Average hydraulic head gradi- ent J 0.07 Geometric mean of K K G 4.48 m/d Water density ρ 1000 kg/m 3 Gravity acceleration g 9.81 m/s 2 Dynamic viscosity μ 8.84 x 10 −4 Pa× s Mean particles diameter d e 2× 10 −4 m Source area dimensions L sy x L sz 70 m× 20 m Source area location along x x inj 50 m Longitudinal dispersivity α x 0.2 m Transversal horizontal disper- sivity α y 0.02 m Transversal vertical dispersiv- ity α z 0.02 m Controlplaneslocationalongx x CP 85 m, 100 m, 120 m, 200 m, 300 m, 400 m, 600 m Table 4.4: Parameter set used for the numerical flow and transport simulations. Value Parameter Symbol PCE TCE DCE VC First-order decay k i 0.005 d −1 0.004 d −1 0.003 d −1 0.001 d −1 Yield coeffi- cient y i 0.79 g g −1 0.74 g g −1 0.64 g g −1 Retardation factor R i 7.1 2.9 2.8 1.4 Table 4.5: Reaction parameters. 48 wherec δ i (τ;x) represents the Dirac-input solution (i.e., instantaneous injection), of the concentra- tion associated with the i th species. For numerical purposes, the source concentration c s (t) can be discretized in step functions according to: c s (t) =c s,0 H(t) + X j=1 Δc s,j H(t−t j ), (4.35) with H(t) indicating the Heaviside step function and Δc s,j = c s,j −c s,j−1 . Equation (4.34) is then given by: c i (t;x) =c s,0 c h i (t;x) + tj<t X j=1 Δc s,j c h i (t−t j ;x). (4.36) In (4.36), c h i (t;x) indicates the cumulative BTC of the i th species computed from an unitary mass source [;;]. Table 4.6 presents the mass transfer parameters adopted in the reactive case study. The contaminant concentration BTCs, recorded from the MC ensemble at the observation CPs, are then employed to compute ¯ c i (x) in the two ways previously illustrated (see Section 4.4.5) and to compute the associated increased lifetime cancer risk. Non-reactive transport In the conservative tracer simulations we instantaneously inject 1 gr of contaminant mass uniformly within A s and observe the contaminant mass flux, indicated as ˙ m [g/d] at the CPs, from which the EPMs under investigation are computed. We recall that contaminant concentrations can then be obtained by the ratio between the contaminant mass flux (measured at the control plane) and the volumetric water flux through the CP. The results of the conservative and reactive transport MC simulations are post-processed to evaluate the impact of φ heterogeneity on risk statistics and illustrated in Section 4.6. 49 Parameter Symbol Value Mass depletion constant Γ 0.1 Initial contaminant mass m 0 900 kg Initial contaminant concentration C 0 0.02 g/m 3 Source first-order degradation constant λ s 5× 10 −5 d −1 Table 4.6: Mass transfer parameters of the source zone selected for the reactive transport simu- lations. 4.6 Results 4.6.3 Analysis for non-reactive contaminants In this Section we show the effect of K-φ heterogeneity on the transport of an inert tracer over the ensemble of MC realizations. We first present the Cumulative Distribution Function (CDF) of early arrival times of contaminants and of peak contaminant mass fluxes at the CPs. We then quantify the relative difference between first and late arrival times of the tracer at the observation locations. We finally show the impact of aquifer connectivity on peak contaminant mass fluxes. In order to understand the importance of realistically accounting for heterogeneous φ, we consider different scenarios: heterogeneousφ fields correlated to spatially variableK characterized by σ 2 Y = 1 and σ 2 Y = 3, and corresponding spatially averaged homogeneous φ coupled with heterogeneous K with σ 2 Y = 1 and σ 2 Y = 3. 4.6.3.3 CDF of first arrival times and peak mass fluxes We present in Figure 4.6 the CDF of dimensionless first arrival times, normalized by the mean advective time required to travel a distance of 3 range along the longitudinal direction, i.e., t 5% U/λ x , where U is the longitudinal mean velocity, given by U = K G J/hφi, with J being the hydraulic head gradient. For simplicity, we only show the CDFs at two CPs located, respectively, at x = 200 m, corresponding to dimensionless distance ζ = 4.29 (in green) and x = 600 m, corresponding to ζ = 15.71 (in red), with ζ = (x CP −x inj )/λ x . Results at the other CPs are alike. Figure 4.6a and 4.6b refers toσ 2 Y = 1 andσ 2 Y = 3, respectively. Dashed curves correspond 4: to heterogeneous correlated K-φ fields and solid curves indicate simplified homogeneous φ and heterogeneous K. Figure 4.6 clearly illustrates that spatially variableφ has a strong impact on the conservative tracer first arrival time statistics. Accounting for φ heterogeneity increases the probability of detecting larger t 5% at the CPs, for both levels of heterogeneity of the Y field. In particular, the importance of including the spatial variability of φ increases as the CPs’ distance from the contaminant injection source increases (compare red curves versus green curves). Indeed, as the distance between the source and CP grows, the plume samples more heterogeneity of the system and the effect of considering variable φ strengthens. We also notice that increasing the heterogeneity of Y (see Figure 4.6b) generally decreases t 5% . Indeed larger σ 2 Y increases the probability of occurrence of fast flow channels which leads to a decrease of first arrival times (compare Figures 4.6a and 4.6b). Subsequent Figure 4.7 pictures the CDF of dimensionless peak mass fluxes ( ˙ m p ), normalized by the initial injected mass (M 0 ) over the time required to travel one range alongx, indicated as ˙ m p λ x /(M 0 U), atthetwoCPslocatedat dimensionless distancesζ = 4.29(in green)andζ = 15.71 (in red). Figure 4.7a refers to σ 2 Y = 1 whereas Figure 4.7b to σ 2 Y = 3. As before, dashed curves correspondtoheterogeneousφwhilesolidcurvesindicatetheoutcomeofconsideringhomogeneous φ. Results at the remaining CPs show comparable outcomes. From Figure 4.7 we observe that the significance of accounting for coupled heterogeneous K-φ fields on transport predictions decreases when the EPM of interest is ˙ m p , as compared to t 5% (Figure 4.6). This effect is more evident when the heterogeneity of the K-field is higher (compare Figure 4.7a with Figure 4.7b). We generally observe that the presence of heterogeneous φ fields increases the chance of observing higher peak mass fluxes as compared to homogeneous φ realizations for both observation locations and levels of heterogeneity of the Y field. In fact the positive correlation between theK-φ fields decreases the macrodispersion of the plume along the longitudinal direction, which produces higher peak mass fluxes at the observation locations. 4; 0 2 4 6 8 10 12 14 t 5% U λ x [−] 0 0.2 0.4 0.6 0.8 1 CDF σ 2 Y = 1 ζ = 4.29 ζ = 15.71 (a) homogeneous φ heterogeneous φ homogeneous φ heterogeneous φ 0 2 4 6 8 10 12 14 t 5% U λ x [−] 0 0.2 0.4 0.6 0.8 1 CDF σ 2 Y = 3 ζ = 4.29 ζ = 15.71 (b) homogeneous φ heterogeneous φ homogeneous φ heterogeneous φ Figure 4.6: Cumulative distribution function (CDF) of the first arrival times,t 5% , of the contam- inant mass at the CPs forσ 2 Y = 1 (a) andσ 2 Y = 3 (b). The CDF of dimensionlesst 5% is indicated with dashed lines for heterogeneous φ fields and with solid lines for homogeneous φ fields. The green color indicates results at the CP located at dimensionless distanceζ = 4.29 whereas the red color pictures results at the CP located at dimensionless distance ζ = 15.71. 52 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ˙ m p λ x M O U [−] 0 0.2 0.4 0.6 0.8 1 CDF σ 2 Y = 1 ζ = 4.29 ζ = 15.71 (a) homogeneous φ heterogeneous φ homogeneous φ heterogeneous φ 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ˙ m p λ x M O U [−] 0 0.2 0.4 0.6 0.8 1 CDF σ 2 Y = 3 ζ = 4.29 ζ = 15.71 (b) homogeneous φ heterogeneous φ homogeneous φ heterogeneous φ Figure4.7: Cumulativedistributionfunction(CDF)ofthecontaminantmasspeak, ˙ m p ,attheCPs for σ 2 Y = 1 (a) and σ 2 Y = 3 (b). The CDF of ˙ m p is indicated with dashed lines for heterogeneous φ fields and with solid lines for homogeneous φ fields. The green color indicates results at the CP located at dimensionless distance ζ = 4.29 whereas the red color pictures results at the CP located at dimensionless distance ζ = 15.71. 53 The decrease of plume dispersion is a result of a diminished variance of the groundwater velocity field in the presence of a positive correlation between φ and K. In fact lower velocity variance implies smaller velocity fluctuations which are the main factors driving dispersion. Lower plume dispersionalongthelongitudinaldirectioninthepresenceofpositiveK-φcorrelationisinlinewith thefindingsof[;8]. Reducedlongitudinaldispersionisalsoinagreementwiththedetectionoflater first arrival times whenφ is modeled as heterogeneous and positively correlated toK (see dashed curves in Figure 4.6) [;8]. Moreover the effect of porosity variability on ˙ m p increases for mildly heterogeneous aquifers (i.e.,σ 2 Y = 1), as compared to an aquifer with higherY heterogeneity (i.e., σ 2 Y = 3). This is also in agreement with the work of [;8]. We also notice lower peak mass fluxes at the CP located further away from A s (i.e., at ζ = 15.71, in red) according to the fact that plume dispersion increases as the plume travels longer and samples more aquifer heterogeneity. The results presented in Figures 4.6 and 4.7 show the importance of incorporating the spatial variability in φ when predicting early arrival times and peak mass fluxes. Since the heterogeneity of the φ field impacts the plume macrodispersion, we analyze the relative plume dispersion for the different scenarios investigated in the following Section 4.6.3.4. 4.6.3.4 Relative plume dispersion The macro-scale dispersive behavior of the solute plume under uniform or spatially random φ fields, coupled with heterogeneous K realizations, is additionally explored in Figure 4.8. In the latter, we present box plots of Δτ, defined as: Δτ = t 95% −t 5% t 95% . (4.37) t 95% is identified as the time of arrival of 95% of the initial injected mass, M 0 . Δτ can represent a measure of the relative plume dispersion since it quantifies the relative difference between the first (i.e.,t 5% ) and late arrival times (i.e.,t 95% ) of the tracer at the CPs. In Figure 4.8 the green 54 Figure 4.8: Box plots of Δτ (4.37) for homogeneous and heterogeneous φ and for σ 2 Y = 1 (a) and σ 2 Y = 3 (b). The green and red colors respectively indicate the results at the CP located at dimensionless distance ζ = 4.29 and ζ = 15.71. boxes indicate the results at the CP located at dimensionless distance ζ = 4.29 whereas the red boxes refer to the outcomes at the CP situated at ζ = 15.71, results at the additional CPs being similar. The thickness of the box plots equals the lag between the first and third quartiles of the probability distribution of Δτ. Figure 4.8a refers to σ 2 Y = 1 whereas Figure 4.8b to σ 2 Y = 3. When we compare the box plots of Δτ for homogeneous and heterogeneous φ, for both CPs locations andσ 2 Y values, we observe that Δτ values are lower whenφ is modeled as heterogeneous. This is because, as explained before, the longitudinal dispersion of the plume decreases when φ is positively correlated to K, leading to a smaller difference between first (t 5% ) and late (t 95% ) arrival times of the plume at the CPs. Accordingly, the range of Δτ presents higher values for 55 all scenarios (homogeneous and heterogeneous φ and different CPs locations) when σ 2 Y = 3, i.e., under higher heterogeneity of theY field we observe a higher chance of fast flow channels (which decrease t 5% ) as well as increased longitudinal dispersion. Finally, in agreement with previous results, the difference between the range of Δτ values between the constantφ and heterogeneous φ cases increases for the CP located further away from the injection area (red box plots) for both levels of heterogeneity ofY. This confirms that the importance of modelingφ variability increases as the tracer travels longer through the heterogeneous aquifer domain. 4.6.3.5 Impact of connectivity on risk Our analysis on the importance of accounting forφ spatial variability on the inert tracer transport statisticsendswithanalyzingtheeffectofaquiferconnectivityonpeakmassfluxesforthedifferent scenarios considered. Previous results (see Figure 4.6 and Figure 4.8) showed how heterogeneity inφ andK affected solute arrival times. As a consequence, these heterogeneities can be expected to impact flow connectivity. We analyze here how K-φ variability affects connectivity and its predictive capabilities on the peak mass flux at a given CP. Connectivity can be described by different metrics [334,3:2] and, in line with the work of [;;], we employ a dynamic connectivity metric, which depends on both groundwater flow and contaminant transport parameters [334]. The selected connectivity metric, indicated as CI [-], is given by the ratio between the effective hydraulic conductivity, K eff , and the geometric mean of K, K G : CI = K eff K G ' 1 t 50% (x cp −x inj )hφi K G J , (4.38) with t 50% being the time of arrival of 50% of the initial injected mass at the CP of interest. The remaining parameters have been defined in Table 4.4. Note that the average φ value (i.e.,hφi) has been employed to compute CI for both the homogeneous and heterogeneous φ simulations. 56 In this analysis, we identify the peak contaminant mass flux ( ˙ m p ) at the observation location (i.e., CP) as a measure of risk. Figure 4.9 shows scatter plots of dimensionless ˙ m p versus CI (4.38). Every circle in Figure 4.9 corresponds to a different MC realization and for simplicity only results at CPs located at dimensionless distance ζ = 1 (in red) and ζ = 4.29 (in green) are reported as the observations at the remaining CPs are similar. Full circles correspond to the outcome of considering homogeneous φ, whereas empty circles indicate the output of more realistic heterogeneousφ conditions. Figure 4.9a refers to σ 2 Y = 1 while Figure 4.9b reports the results of σ 2 Y = 3. A linear regression between ˙ m p and CI is identified for both homogeneous (full circles) and heterogeneous (empty circles)φ fields, CP location andσ 2 Y value. Indeed, under higherCI values, corresponding to the presence of preferential flow channels, ergo lowert 50% (4.38), a higher chance of detecting higher ˙ m p values at the CPs under investigation is identified. We observe that the slope of the linear regression line, whose equation is reported in gray for homogeneous φ and in black for heterogeneous φ, is higher when φ is heterogeneous and positively correlated to K as compared to simplified homogeneous φ for both CPs and for both levels of heterogeneity of Y. Therefore, under realistic heterogeneous φ conditions and under the assumption that the K-φ relationship can be described through the KC equation (4.3), an increased chance of observing higher ˙ m p , i.e. having higher risk, emerges as compared to constantφ conditions, for the sameCI value. This could be attributed to the lower macro-scale spreading/dispersion of the contaminant plume under positive K-φ correlation, which results in higher ˙ m p at the sensitive location. We also notice that a larger σ 2 Y (Figure 4.8b) favors the formation of preferential channels, hence higher CI values and corresponding ˙ m p values (compare Figures 4.9a and 4.9b). Results at the remaining CPs locations are similar and we generally observe a decrease of the coefficient of determination (R 2 ) of the regression line between CI and ˙ m p values as the CP distance from A s increases. 57 0 1 2 3 4 5 CI [−] 0 1 2 3 4 5 ˙ m p λ x M O U [−] σ 2 Y = 1 y = 1.030x - 0.023 y = 0.347x - 0.050 y = 1.295x - 0.094 ζ = 1 y = 0.472x - 0.074 ζ = 4.29 (a) homogeneous φ homogeneous φ heterogeneous φ heterogeneous φ 0 2 4 6 8 10 12 CI [−] 0 2 4 6 8 10 12 ˙ m p λ x M O U [−] σ 2 Y = 3 y = 0.827x + 0.057 y = 0.242x - 0.029 y = 1.012x - 0.034 ζ = 1 y = 0.327x - 0.066 ζ = 4.29 (b) homogeneous φ homogeneous φ heterogeneous φ heterogeneous φ Figure 4.9: Peak mass flux, ˙ m p , versus connectivity metric, CI (4.38), for σ 2 Y = 1 (a) and σ 2 Y = 3 (b). The red color indicates results at the CP located at dimensionless distance ζ = 1 whereas the green color pictures results at the CP located at dimensionless distance ζ = 4.29. Full circles indicate homogeneous φ and empty circles indicate heterogeneous φ. The regression line is plotted in gray for homogeneousφ and in black for heterogeneousφ for both CPs locations. The coefficient of determination (R 2 ) of the regression line for σ 2 Y = 1 at the CP located at dimensionless distance ζ = 1 are R 2 = 0.778 for homogeneous φ, R 2 = 0.599 for heterogeneous φ and R 2 = 0.656 for homogeneous φ, R 2 = 0.444 for heterogeneous φ at the CP located at dimensionless distance ζ = 4.29. The coefficient of determination (R 2 ) of the regression line for σ 2 Y = 3 at the CP located at dimensionless distance ζ = 1 are R 2 = 0.871 for homogeneous φ, R 2 = 0.816 for heterogeneousφ andR 2 = 0.845 for homogeneousφ,R 2 = 0.728 for heterogeneous φ at the CP located at dimensionless distance ζ = 4.29. 58 Our results suggest that realistic heterogeneousφ conditions are associated with higher values of ˙ m p , corresponding to higher health risk at sensitive locations, therefore considering simplified constant φ conditions could be misleading, i.e. lead to underestimated risks. We also foresee that the effects of incorporating porosity variability on transport observables could potentially be emphasized in the presence of non-GaussianK fields which present higher probability of well- connected zones of high K values [:5,;;]. For example, the work of [345] showed that higher solute peak concentrations are detected at extraction wells as the departure of the non-Gaussian Y fields from a Gaussian structure increases. 4.6.4 Analysisforreactivecontaminants: implicationsinincreasedlifetime cancer risk Results in this Section refer to the reactive contaminant transport simulations, i.e. to the degra- dation chain of chlorinated solvents (PCE→ TCE→ DCE→ VC). We analyze the statistics of the total increased lifetime cancer risk, R T (x) (4.33), for the scenarios under investigation. We first present the spatial distribution of the low-order moments, i.e. the ensemble mean and the coefficient of variation, of R T , and we then show the complete statistical characterization of the PDF of R T . The outcomes of accounting for heterogeneous φ are represented with dashed lines while the results of constantφ are illustrated with solid lines in Figures 4.:-4.32. The results of considering σ 2 Y = 1 are pictured in magenta and the outcomes for σ 2 Y = 3 are illustrated in blue. 4.6.4.3 Low order moments of R T (x) Figure 4.: shows the spatial evolution of the mean value of R T (4.33), denoted ashR T i, over the ensemble of MC simulations, when the peak concentration of the contaminanti at the CP (Figure 4.:a) or the maximum running averaged concentration over the ED time (4.32) (Figure 4.:b) is adopted to compute R T (x).hR T (x)i quantifies the expected threat to the exposed community. 59 0 2 4 6 8 10 12 14 16 0 1 2 3 10 -5 0 2 4 6 8 10 12 14 16 0 1 2 3 10 -5 Figure4.:: Longitudinalspatialdistributionofthemeantotalincreasedlifetimecancerrisk,hR T i, when the maximum concentration (a) or the maximum running averaged concentration over the exposure duration (ED) time (b) are considered to compute R T . The results for σ 2 Y = 1 are represented in magenta and the results forσ 2 Y = 3 are pictured in blue. Dashed curves correspond to heterogeneous φ fields and solid curves picture homogeneous φ. 5: We first notice that, for the current numerical set-up, the formula adopted to compute the flux-averaged contaminant concentration in equation (4.;) has only a very minor influence on the spatial evolution ofhR T i (compare Figure 4.:a and 4.:b). We indeed notice thathR T i is just slightly lower when we employ the maximum running averaged concentration to compute the risk (Figure 4.:b), according to the fact that adopting the peak concentration represents the most conservative approach. We observe thathR T i increases with CP distance since more toxic species are produced with time along the chemical degradation chain of PCE [;;] in Figures 4.:a and 4.:b. In fact, cancer potency factors (CPF i ) increase during the degradation of PCE until VC is produced (see Table 4.3). In general, the influence of heterogeneity on the mean ILCR is almost unnoticeable in Figure 4.: (see dashed lines versus solid lines). The fact that accounting for φ heterogeneity has only a small impact onhR T (x)i confirms that, also in the presence of reactive contaminants, the influence of incorporating φ variability depends on the metric of interest. Increasing the heterogeneity of the Y field (i.e., σ 2 Y = 3, blue lines) decreases the magnitude of the mean ILCR along the longitudinal direction because of enhanced dispersion of the plume and subsequent dilution at CPs [;;]. Figure 4.; shows the spatial evolution of the coefficient of variation of the total ILCR (i.e., CV R T (x)). We recall that CV R T (x) = σ R T (x)/hR T (x)i, with σ R T (x) indicating the standard deviation of the total ILCR. Figure 4.;a refers to the computation of R T (x) by considering the maximum contaminant concentration whereas Figure 4.;b shows the same results when the maximum ED flux-averaged concentration (4.32) is adopted to compute R T (x). The coefficient of variation, CV R T (x), measures the uncertainty associated with the risk estimate. In agreement with what observed in Figure4.;, we also point out that the approach adopted to compute the flux-averaged concentration in equation (4.;) does not influence the risk uncertainty and the differences between Figure 4.;a and Figure 4.;b are imperceptible. We note that this remark is applicable to the current numerical set-up. We generally observe that the spatial evolution ofCV R T presents an inverse correlation withhR T (x)i, in agreement with [;;]. Therefore 5; 0 2 4 6 8 10 12 14 16 0 0.2 0.4 0.6 0.8 0 2 4 6 8 10 12 14 16 0 0.2 0.4 0.6 0.8 Figure 4.;: Longitudinal spatial distribution of the coefficient of variation of the total increased lifetime cancer risk, CV R T , when the maximum concentration (a) or the maximum running av- eraged concentration over the exposure duration (ED) time (b) are considered to compute R T . The results forσ 2 Y = 1 are represented in magenta and the results forσ 2 Y = 3 are pictured in blue. Dashed curves correspond to heterogeneous φ fields and solid curves picture homogeneous φ. weidentifylowermeantotalILCRandmoreuncertainriskpredictionsclosetothesourcewhilewe observe higher mean risk and lower risk uncertainty further away for all the simulated scenarios. Larger risk uncertainty close to the source is expected since the plume is mainly controlled by advection at that stage [67]. Our results indicate that accounting for the heterogeneity inφ affects the risk uncertainty. Indeed, the positive correlation betweenφ andK leads to a reduction of the sample-to-sample fluctuation of the risk, and therefore to a reduction of CV R T (x) for both levels of heterogeneity of theY field (see dashed lines versus solid lines in both Figures 4.;a and 4.;b). This effect is emphasized when σ 2 Y = 3. 62 0 0.5 1 1.5 2 2.5 3 3.5 4 10 -5 0 0.5 1 1.5 2 10 5 0 0.5 1 1.5 2 2.5 3 3.5 4 10 -5 0 0.5 1 1.5 2 10 5 Figure 4.32: Probability distribution function (PDF) of the total increased lifetime cancer risk, R T , at the CP located at ζ = 10 (a) and ζ = 15.71 (b). The results for σ 2 Y = 1 are repre- sented in magenta and the results for σ 2 Y = 3 are pictured in blue. Dashed curves correspond to heterogeneous φ fields and solid curves picture homogeneous φ. TheresultspresentedinFigures4.:and4.;showthattheinfluenceofφheterogeneityishigher on risk uncertainty, i.e., CV R T (x) (Figure 4.;) than on the mean total ILCR,hR T (x)i, (Figure 4.:). This can be justified by the fact that mean mass fluxes at large CPs are not significantly affected by local fluctuations of the flow field, as explained in [8:]. On the other hand, the first order analysis provided in [8:] and [5] shows that the variance of the mass flux over a control plane is sensitive to local heterogeneity. 63 4.6.4.4 PDF of R T We conclude the discussion on the influence of φ spatial heterogeneity on reactive transport by presentingthePDFoftheILCRinFigure4.32. TheresultsinFigure4.32refertothecomputation of R T (x) by employing the maximum running averaged concentration over ED (4.32). Results at the CP located at ζ = 10 are reported in Figure 4.32a, while the outcomes at ζ = 15.71 are presented in Figure 4.32b. The PDFs of R T present a similar “Gaussian”-like shape under constant or random φ. These results are consistent with the findings of [;;,322]. We observe that the impact ofφ heterogeneity on the PDF of the ILCR, described below, is more pronounced when the CP is situated closer to the source for the case of σ 2 Y = 3 (compare Figure 4.32a with Figure 4.32b). In general accounting for φ heterogeneity (dashed PDFs) only slightly influences the mean value of R T but significantly affects the uncertainty ofR T . We indeed notice that theR T PDF is less spread when the heterogeneity ofφ is modeled as compared to assuming constantφ (dashed PDFs versus solid PDFs). In particular the probability of observing values of R T aroundhR T i increases of more than 30% when the variability of φ is included in models’ predictions for all scenarios considered except for the case ofσ 2 Y = 3 at the CP located atζ = 15.71, where this probability is around 20% higher under heterogeneous φ. Moreover, the full risk PDFs show the probability of observing extreme values of R T . We observe that modeling the φ heterogeneity has influence on both the highest and lowest values of R T and the probability of detecting extreme values of R T decreases in the presence of spatially variable φ. 4.7 Conclusions In this study we analyze the impact of the coupled spatial variability ofφ andK on contaminant transport and probabilistic risk analysis. We employ a Monte Carlo framework to statistically analyze key EPMs, detected at observation locations. We analyze the impact of φ heterogeneity 64 on contaminant arrival times and peak mass fluxes of tracers at control planes as well as on risk and on connectivity. We then show the implications of accounting forφ variability in probabilistic human health cancer risk analysis by considering human exposure to chlorinated solvents through the ingestion pathway. Our work leads to the following main key conclusions: • The positive correlation betweenφ andK leads to a decrease of the plume macrodispersion with consequent higher peak mass fluxes, later first arrival times and earlier late arrivals at control planes. Higher peak mass fluxes correspond to higher health risk therefore including the spatial variability of φ in models’ predictions is essential for risk assessment associated with aquifer contamination and remediation. Moreover, a more accurate quantification of early and late arrival times is fundamental to respectively detect the potential hazard of a contaminant and to control the end-point of remediation [4,66]. • The relevance of incorporating the spatial variability of theφ field can depend on the metric of interest, for example our study shows that this relevance is emphasized when investigating contaminants’ arrival times over peak mass fluxes at the observation control planes. This hasimportantimplicationsboth toset-up modeling studiesand in decisionmaking processes (e.g., inriskassessmentoraquiferremediationactivities)inthecontextofabetterallocation of resources. Indeed, available resources need to be prioritized and assigned in a defensible manner to ensure that unacceptable identified risks are reduced to acceptable levels. • Modeling the φ field spatial variability as positively correlated to the heterogeneous K field significantly reduces the cancer risk uncertainty but only minorly influences the mean value of the total cancer risk at control planes. Our findings constitute then an important contribution for cancer risk uncertainty quantification however, adopting more simplified and conservative models which assume homogeneous porosity does not highly impact the mean value of the total cancer risk. 65 We finally noticed that the effect of considering φ heterogeneity can be affected by the level of heterogeneity of the K field and by the source-to-CP distance. We recall that the results of our work are limited to positive correlation between φ and K expressed through the KC equation (4.3). Other empirical models can be adopted to compute theK-φ relationship and could potentially lead to different conclusions. Furthermore, other factors such as pumping operations [344] and solute mass release rates at the source zone [58] can impact the relative importance of hydrogeological heterogeneity. On a final note, we emphasize that our work represents a first step towards a more comprehensive global sensitivity analysis (GSA) that could lead to more general conclusions. Usingpolynomialchaosexpansion(PCE)withinaGSAcanreducethecomputational burden associated with Monte Carlo simulations. Utilizing PCE to perform GSA would allow to better understand the relevance of heterogeneous porosity with respect to heterogeneous hydraulic conductivity. Under the context of risk, PCE-based GSA was employed by different authors [46,47,373]. [46] employed GSA to identify the influence of uncertain hydrogeological parameters on the first two statistical moments of the peak concentration. [47] applied a GSA through PCE to three transport models to identify the relative importance of model-dependent parameters, the space-time locations where the models are more sensitive to these parameters and to assist with parameters estimation. [373] illustrated their GSA approach, based on arbitrary polynomial chaos expansion (aPC), on a 3D groundwater quality and health risk problem in a heterogeneous aquifer. Incorporating a GSA through PCE into the problem investigated in this work is topic of future research with the goal of reducing the computational burden normally associated with Monte Carlo simulations and to investigate the impact of considering different soil types on risk. After illustrating the importance of characterizing aquifers’ hydrogeological properties in con- taminant transport predictions, the next chapters consider the presence of realistic external stresses to the heterogeneous groundwater system. The impact of external factors on ground- water contaminant studies is presented in the following chapters and we start by showing the 66 significance of incorporating realistic pumping schedules when measuring contaminant BTCs at the extraction well in chapter 5. 67 Chapter 5 Influence of pumping operational schedule on solute concentrations at a well in randomly heterogeneous aquifers The content of this chapter is published in Journal of Hydrology [344]. 5.3 Introduction Pumping wells are widely employed for groundwater supply in the context of urban, agricultural, and industrial activities. Water management agencies typically schedule groundwater extraction through a well-defined temporal sequence of pumping rates to maximize benefits to anthropogenic activities and minimize the environmental footprint of the withdrawal process. The temporal dynamics of well pumping operations induce changes in the subsurface flow field which in turn can affect spreading and mixing of contaminant plumes. A key challenge is to provide estimates of contaminant concentration at an environmentally sensitive location while accounting for the joint effect of the temporal dynamics of an operating pumping well and the (typically uncertain) spatial variability of the hydrogeological properties of the subsurface reservoir within which a potentially harmful solute is residing and migrating. 68 The effect of the spatial variability of hydraulic properties of aquifers, such as hydraulic con- ductivity, on contaminant transport has been broadly studied [55,:3,3:4]. Due to the high cost of data acquisition campaigns and limited financial resources, a highly detailed characterization of the spatial distribution of hydrogeological attributes (e.g., hydraulic conductivity) of an aquifer system at thefield scale ofinterest is virtually unfeasible. As a consequence, this lack of knowledge leads to uncertainty associated with predicted values of contaminant concentration at environ- mentally sensitive and strategically valuable locations, such as pumping wells. As noted in the following, while a considerable amount of published works address the analysis of contaminant transport under steady-state pumping, the way a given time-varying pumping schedule affects contaminant plume behavior in heterogeneous aquifers is tackled only marginally. Uncertainty in the quantification of the extent of pumping well protection regions (both in terms of well catchment and time-related capture zones) in the presence of random spatial dis- tributions of aquifer parameters and under the assumption of steady-state background ground- water flow to which one or more wells extracting a constant flow rate are superimposed has been addressed [49,88,95,:9,338,339,38;,436,437] through numerical Monte Carlo simulations (un- conditional or conditional on data). Stauffer et al. [424] performed a theoretical analysis based on first-order approximation and a Lagrangian approach to estimate the uncertainty associated with the delineation of well catchments in the presence of a spatially variable random hydraulic conductivity field. Lu and Zhang [352] used a Lagrangian framework for the study of time-related capture zones in spatially randomly heterogeneous aquifers. Riva et al. [392] used stochastic moment equations to analyze the uncertainty associated with well catchments. The work of de Barros et al. [64] addressed the value of hydrogeological information to assess the risk associated with pollution of an operating pumping well. Molson and Frind [35;] applied the concept of mean life expectancy (i.e., the time between plume detection and plume capture by a well) and capture probability to delineate well capture zones. The life expectancy approach constitutes an alter- native to delineate time-of-travel capture zones for wellhead protection. Pedretti and Fiori [382] 69 developed an analytical solution of contaminant transport in heterogeneous porous media under purely convergent flow to analyze solute breakthrough curves (BTCs) at the well. Indelman and Dagan [325] analyzed the spreading of a contaminant plume injected in a heterogeneous formation through a well and advected by a steady-state groundwater velocity. Onlyalimitednumberofworksconsidertheimpactoftransientflowregimesonthedelineation of capture zones or wellhead protection areas (WHPA). Amongst these, [386] used numerical modelstodelineatetheWHPAoftwomunicipalwellsundertheinfluenceofnearbyirrigationwells. Their conclusion was that seasonal variation in pumping rates affected the water table depth and the extent of a target time-related well capture zone. In this sense, significant seasonal variations inpumpingoperationsshouldbeexplicitlyembeddedinaprocedureforWHPAdelineation. Reilly and Pollock [388] employed a numerical flow model and a particle tracking method to investigate the effects of external stresses, i.e. recharge (that is typically cyclic over time), on the area contributing recharge to the wells. They observed that the ratio between the mean travel time of the solute to the well and the period of cyclic stress can provide a quantification of the importance ofconsideringthetransienteffectsofacyclicrecharge. Theauthorsnotedthat: (a)approximating the system behavior through an average steady-state condition yields results which are similar to those obtained by considering the actual transient pattern of recharge when the above ratio is much larger than one; and (b) considering transient flow conditions is important to characterize the proper behavior of solutes released in the proximity of the boundary of the well. Jacobson et al. [326] incorporated the effects of uncertainty in the transmissivity and the magnitude and direction of the hydraulic head gradient in the steady-state analytical solution obtained by [36] for capture zones delineation. These authors observed that uncertainty in the magnitude of the mean regional flow propagates to the uncertainty in the length (along the direction of the regional flow) of the time-dependent capture zone. Otherwise, uncertainty in the mean direction of the background flow influences the uncertainty in the maximum width (along a direction normal to the average regional flow) of the capture zone. Festger and Walter [87] analyzed the same problem 6: through a semi-analytical approach. They concluded that temporal variations in the direction of the background natural hydraulic gradient typically cause an expansion of a given capture zone when compared against its steady-state counterpart, the opposite being observed in some specific cases (e.g., in the presence of moderate directional variations and large gradient magnitudes). Neupaueretal.[369]investigatedtheeffectsofaquiferheterogeneityonthespreadingofaplumeof conservative solute which is advected within a temporally variable flow field induced by engineered injectionandextractionsequences. Otherauthors,e.g. Assoulineetal.[7],illustratethattransient effects associated with high-frequency irrigation affect solute and water regimes and residence times under partially saturated conditions. Additionalexamplesofcontaminanttransportstudieswhichincorporatetheeffectsoftransient pumping strategies include the works of [44], who simulated unsteady flow to optimize time- varying pumping rates for groundwater remediation; Vesselinov [438], who studied the influence of transient flow and of uncertainty in longitudinal and transverse dispersivities on the delineation of well capture zones; Chen et al. [45], who evaluated solute transport in divergent radial flow fields created by multistep pumping; and Leray et al. [33:], who detected the occurrence of non- negligible solute concentrations at a pumping well following induction of a transient flow field by the instantaneous activation of the well. In contrast to our work, Leray et al. [33:] did not investigate the uncertainty of the contaminant concentrations recovered at the well. Notably, in this broad context, the effects of a temporally varying pumping regime on contaminant solute BTCs detected at the location of the pumping well operating in a heterogeneous aquifer has not been systematically addressed. Our work is then geared towards a first exploration of the joint effects of (a) random spatial variability of hydraulic conductivity and (b) transient flow regime, as induced by a temporally varying pumping schedule, on this key feature of contaminant transport. These results will then be employed in the context of risk quantification at the well, which is considered as a sensitive strategic location with social and economic value. Risk is here quantified in terms of the probability that an established threshold for the concentration of a 6; target chemical species is exceeded. As such, the objective of our study is directly related to the assessment of risks associated with contamination of drinking water extracted at the well (see. e.g., the works of [67,79,99,;5]). 5.4 Problem Formulation We consider fully saturated groundwater flow taking place in a heterogeneous porous medium characterized by a spatially variable isotropic random transmissivity field T (x) [L 2 /T ], uniform porosityφ [−], constant thicknessb [L] and a uniform-in-the-average natural base flowQ 0 [L 3 /T ]. The location vector is here expressed as x = (x,y), in a Cartesian coordinate system. In Figure 5.3 depicts a sketch of the domain and problem we study. The log transmissivityY (x) = ln(T (x)) is assumed to be a normally distributed and a weakly stationarity correlated random field, which is fully described by its mean, m Y , and covariance function, C Y (e.g., [3:4]). The covariance function is statistically isotropic and characterized by the variance of Y (x), denoted by σ 2 Y , and its (isotropic) correlation scale λ [L]. A pumping well located at x w = (x w ,y w ) operates with pumping rate Q w (t) [L 3 /T ], with t [T ] denoting time. Therefore, the flow field is described by the depth-averaged saturated groundwater flow equation: S ∂h(x,t) ∂t =∇· [T (x)∇h(x,t)] +Q w (t)δ(x−x w ) (5.3) whereh [L]isthehydraulichead,S [−]correspondstostorativityandδ istheDiracdeltafunction. Fixed hydraulic head values are imposed on the west (left) and east (right) boundaries and no- flow boundary conditions are set on the north and south boundaries of the domain, resulting in a uniform mean base flow along the x-direction (Figure 5.3). A tracer with initial uniform concentration c 0 [M/L 3 ] is instantaneously injected at time t 0 [T ] within a region of areaA 0 [L 2 ]. Transport of the tracer plume is described by the standard advection-dispersion equation: 72 L x 2 x 3 x 1 L 2 L 1 Pumping well Capture zone Contaminant Source L y L x q 0 =Q 0 /bL y y x x 2 Figure 5.3: Schematic representation of the problem studied. A two-dimensional flow field charac- terized by the presence of a pumping well and a uniform (in the mean) base flow from left to right of the domain is considered. The contaminant source is rectangular with longitudinal dimension L 1 and transverse dimensionL 2 , located at distanceL from the well positioned atx w = (x w ,y w ). The contaminant source is aligned with respect to the pumping well center. 73 ∂c(x,t) ∂t =−v(x,t)·∇c(x,t) +∇· [D∇c(x,t)]; c(x,t 0 ) =c 0 for x ∈A 0 . (5.4) Here, c [M/L 3 ] is solute resident concentration, v = q/φ [L/T ] is effective velocity (and q [L/T ]isthespecificdischarge)andD [L 2 /T ]isthelocal-scaledispersiontensorwithlongitudinal component D x and transverse component D y . We consider that solute transport takes place sufficiently far away from the boundaries (Figure 5.3). Once we obtain the space-time distribution of the concentration we evaluate the flux-averaged concentration C [M/L 3 ] at the pumping well [337] and obtain the solute BTC. Details related to the computational method used to evaluate equations (5.3) and (5.4) are provided in chapter 5.5. Equations (5.3) and (5.4) are expressed in dimensional form. In this work, we will analyze our results in terms of relevant dimensionless groups. We start by the following functional relationship between system parameters and state variables of interest (see also, e.g., [95]): f 1 (t,x,y,T,b,q 0 ,C,Q w ,φ,D) = 0; (5.5) withq 0 [L/T ] being the imposed uniform Darcy’s background groundwater specific discharge. We consider q 0 , Q w , b and φ as deterministic in our analysis. The following relationship between dimensionless quantities can then be written under the assumption that φ and D are constant: f 2 tq 2 0 b Q w , xq 0 b Q w , yq 0 b Q w , T q 0 b , C C∗ = 0; (5.6) Note that the dependence on a typical Péclet number would appear in equation (5.6) in case one would consider analyzing the interplay between advective and dispersive processes in the system (see chapters 5.5 and 5.6.6 for a discussion on this aspect). The reference concentration C ∗ in equation (5.6) is taken as the contaminant concentration limit value for a target chemical species. The value of the latter is established, for example, by 74 regulatory bodies/agencies. In the regulatory context, this contaminant concentration limit can represent the Maximum Contaminant Level for drinking water. Note that Q w in equation (5.6) indicates the constant pumping well extraction rate. In this work, the constant pumping rate, which is selected in the steady-state setting, coincides with the average (over the temporal window of well operation) of the pumped flow rate in the transient scheme (see chapter 5.5). Adopting the constant (in time) extraction rate is then a natural choice to normalize the various terms in equation (5.6) for both the steady-state and transient settings analyzed. Due to the randomness of Y (x), the feedback between a transient pumping regime and con- taminant breakthrough curves must be studied in a stochastic framework. Chapter 5.5 provides details of the computation of the flux-averaged concentration at the well. 5.5 Methodology We consider a rectangular solute source, with longitudinal (along direction x) dimension L 1 and transverse (along directiony) dimensionL 2 , within which solute is instantaneously released in the system (i.e. the area of the source isA 0 =L 1 ×L 2 , see Figure 5.3). A pumping well is located at a distanceL (measured along thex-direction) from the contaminant source and aligned with the source centroid location. The computational domain has size Ω = [(x,y)|x∈ (0,L x );y∈ (0,L y )] with L x = 170 m and L y = 150 m. We considered a unit thickness aquifer in all simulations. In order to quantify the statistics ofC, we employ a numerical Monte Carlo (MC) framework. A collection of realizations of the Y (x) field are generated using a sequential Gaussian simulator (SGeMS; [389]). The spatial covariance C Y is modeled through a stationary Exponential model, characterized by isotropic correlation scale λ = 8 m and variance σ Y 2 = 3. Consistent with pre- vious studies (e.g., [64,363]), a uniform generation grid, with elements of side Δx = Δy = 1/8 λ, is employed. Transient groundwater flow and contaminant transport for each realization are respectively solved through the widely tested codes MODFLOW [;4] and MT3DMS [456]. A 75 Lagrangian-Eulerian method, the Method of Characteristics (MOC), is applied to solve contam- inant transport with MT3DMS. Further details pertaining the computational methods (and the associated numerical errors) adopted are documented in the literature [455]. A comparison be- tween the numerical methods utilized in this work with other existing codes can be found in the literature [;9]. We explore three diverse pumping scenarios, indicated as scenarios I, II and III. Each scenario is characterized by a fixed total volume of water extracted (V w ) and a given operational scheme according to which the pumped flow rate is varied in time. Each of the transient scenarios is compared against a corresponding setting associated with steady-state pumping characterized by the same total volume of water extracted from the system during the simulated transient regime. Table 5.3 lists the main parameters used in the study. Note that the longitudinal Péclet number is set to 800. The Péclet number is given by P e = λv 0 /D x , where v 0 = q 0 /φ denotes the imposed uniform-in-the-mean longitudinal velocity and D x is the local scale longitudinal dispersion coefficient, which we set at the constant value listed in Table 5.3. We perform MC numerical simulations of flow and contaminant transport according to the three scenarios displayed in Figure 5.4 (Figure 5.4a, b, and c respectively referring to scenariosI, II and III). The selected well pumping operation schedules are inspired by patterns employed in practical settings by groundwater management agencies to extract water from underground resources for diverse purposes (including, e.g., drinking and irrigation). As seen in Figure 5.4, scenarioIII ischaracterizedbythehighestvariabilityofthedynamicpumpingstrategy. Following a preliminary convergence study (not shown), we base our results on a set of 1000 MC simulations for each scenario and pumping strategy adopted. The results of the MC simulations are post- processed with the aim of capturing the temporal evolution of the mean and variance of the flux- averaged contaminant concentration recovered at the well. We note that the effects of variable well pumping operations on the solute breakthrough curve at the well can also depend on the ratio between the natural base flow Q 0 and the well pumping rate Q w . 76 (L x ,L y ) Flow domain (170 m, 150 m) (Δ x , Δ y ) Grid size 1/8 λ m, 1/8 λ m λ Correlation length of Y 8 m σ 2 Y Variance of Y 3 φ Porosity 0.2 J Mean head gradient 0.005882 T Average Transmissivity 1 m 2 /d D x Longitudinal dispersion coefficient 2.94 x 10 −4 m 2 /d D y Transverse dispersion co- efficient 2.94 x 10 −6 m 2 /d (x w ,y w ) Location of the pumping well (114.5 m, 74.5 m) Q w Temporal sequence of pumping rate for scenario I 0.6 (uniform pumping rate) vs 0.8/0.4/0.8/0.4/0.8/0.4/0.8/0.4/0.8/0.4/0.8/0.4/ 0.8/0.4/0.8/0.4 m 3 /d Q w Temporal sequence of pumping rate for scenario II 0.45 (uniform pumping rate) vs 0.8/0.7/0.6/0.5/0.4/0.3/0.2/0.1/0.8/0.7/0.6/0.5 0.4/0.3/0.2/0.1 m 3 /d Q w Temporal sequence of pumping rate for scenario III 0.3 (uniform pumping rate) vs 0.8/0/0.4/0/0.8/0/0.4/0/0.8/0/0.4/0/0.8/0/0.4/ 0 m 3 /d q o Groundwater specific dis- charge 5.882 x 10 −3 m/d V w Volume of extracted water for scenario I 2880 m 3 V w Volume of extracted water for scenario II 2160 m 3 V w Volume of extracted water for scenario III 1440 m 3 C ∗ Contaminant concentra- tion threshold 10 g/m 3 Table 5.3: Parameter set employed in the numerical simulations. 77 Figure 5.4: Pumping scenarios analyzed in the study: Scenario (a) I; (b) II; (c) III. 78 For each scenario we quantify the risk, defined here as the probability of failure, e.g., the prob- ability that the contaminant concentration exceeds a given critical value at the well. According to this, we define risk (ξ) in this work in terms of exceeding probabilities [62]: ξ =P (C(x w ,t)≥C crit ), (5.7) where C crit ≡ C ∗ , as introduced in chapter 5.4, and P denotes probability. Note that the risk definitioninequation(5.7)restsonafullstatisticaldistributionoftheflux-averagedconcentration at the pumping well as a function of time. As shown by several authors [62,8;,329], the mean concentration is insufficient to characterize risk and higher order moments are required. The mean concentration only becomes a meaningful measure of risk under ergodic conditions, i.e., large travel distances or large contaminant source dimensions with respect to the log-transmissivity correlation scale. If these conditions are satisfied, macrodispersion theory is expected to hold ( [52,32;]; chapter 10 of [3:4] and references therein). In particular, [62] provides a detailed discussion on the need to use of a full probabilistic distribution to quantify risk by taking into account also the probability of extreme events (e.g., probability of concentration exceeding a very high or low value). The exceedance probability 5.7 is computed directly from the numerical MC results and informs decision makers on the probability that C(x w ,t)≥C crit . 5.6 Results and discussion We start by illustrating the relevance of the problem tackled in this study through a discussion of some key differences between contaminant BTCs observed at the pumping well operating with a steady-statepumpingstrategyandatransientpumpingschemeforselectedexemplaryrealizations of theY (x) field. Second, we present the temporal evolution of the mean and variance (as well as theresultingcoefficientofvariation)oftheflux-averagedconcentrationatthewellforthepumping strategies of scenarioIII. The latter is representative of the main conclusions of our analysis, the 79 results of scenarios I and II being fully documented in Appendix B (Supplementary Material). Third, we display the temporal behavior of risk, as defined in equation (5.7), for scenario III, under steady-state and transient pumping conditions. Corresponding outputs for scenariosI and II are included in the SM. Finally, we demonstrate the effects of local scale dispersion and aquifer heterogeneity on the uncertainty of the concentration BTCs for scenario III. We show how the interplay between pumping operations, aquifer heterogeneity and local scale dispersion affect the temporal patterns of the flux-averaged concentration uncertainty. The upcoming results also emphasize the importance of pumping in diluting the concentration at the well. 5.6.3 Concentration breakthrough curves observed at the pumping well for different pumping strategies In this chapter, we aim to provide a qualitative illustration of the way the concentration BTC at the well can be affected by the pumping operation and the heterogeneous transmissivity field. A graphical depiction of the effect of a time-dependent withdrawal schedule on the concentration BTC is provided in Figure 5.5. The latter depicts the concentration BTCs observed at the pumping well for three diverse realizations of Y (x) under the constant and temporally varying pumping strategy of each scenario. Each plot within Figure 5.5 includes the contaminant BTCs resulting from a constant and a transient extraction strategy for a given realization of the Y (x) field. Figures 5.5a, 5.5b, 5.5c refer to scenario I, Figures 5.5d, e, f refer to Scenario II and finally, Figures 5.5g, 5.5h, 5.5i refer to Scenario III. The three distinct realizations of the transmissivity are denoted byT 1,T 2 andT 3 (see Figures SM39-SM3; of the SM) in Figure 5.5 and are selected solely for the purpose of illustration. These three realizations of the transmissivity fields (T 1,T 2 andT 3) were obtained through the software SGeMS [389] and are characterized by the geostatistical parameters described in chapter 5.5. As shown in Figure5.5, significant differences are detected between the concentration BTCs obtained 7: with uniform and variable in time pumping rates. From a qualitative standpoint, it can be noted that the concentration BTCs produced by a dynamic pumping regime are characterized by peaks of different strength and occurring at diverse times when compared to the corresponding BTCs obtained with uniform pumping. This result has strong implications in risk analysis, inasmuch as public health authorities are concerned with peak concentrations and the time of their occurrence [66,373,3;9]. TheresultsdisplayedinFigure5.5alsoillustratethecomplexinteractionbetweenthe hydrogeological heterogeneity of the aquifer and the transient flow regime induced by a temporally varying pumping schedule. These illustrative results point out the importance of considering temporally varying pumping settings, as employed in practical applications, when evaluating contaminant transport in a well field, with direct implication on risk caused by contamination of water extracted at the well. 5.6.4 Mean and variance of solute concentration at the pumping well We compute and discuss here the temporal evolution of the first two statistical moments of the flux-averaged solute concentration detected at the pumping well. To streamline the presentation, we only show the results related to scenario III, which encapsulates the key observations of the analysis. Results associated with scenarios I and II are documented in the SM. Figures 5.6 and 5.7 depict the mean and the variance of the flux-averaged concentration (respectively denoted ashCi andVar[C]) detected at the well following the adoption of a uniform (Figures 5.6a and 5.7a) and dynamic (Figures 5.6b and 5.7b) pumping strategy. These results indicatethatthetransientoperationalregimeappliedatthewellhasmarkedinfluenceontheshape of the concentration BTCs. A striking observation is that the local minimum and maximum of the mean and variance of the flux-averaged concentration at the well occur at the same dimensionless time for both types of pumping regimes. This observed behavior is supported by theoretical developments illustrated in Appendix A. 7; Figure 5.5: Pumping scenarios analyzed in the study: Scenario (a) I; (b) II; (c) III. 82 Figure5.6b shows that the mean flux-averaged concentration starts increasing when the plume approaches the well, until a maximum value is reached between dimensionless time values of 0.2 and 0.3. Figure 5.6b also reveals that an increase in the pumping rate can lead to a sudden decrease of the mean flux-averaged concentration. Due to the intermittent pumping schedule, this effect is seen to extend in time causing an oscillatory pattern of the mean BTC. It can be noted that activating the pumping well induces contaminant dilution at the well due to mixing of contaminated and clean water. This leads to a decrease of the contaminant concentrations in the water extracted at the well with respect to concentration levels in the contaminant plume [67,77]. A pronounced oscillatory pattern in time is detected for Var[C] (Figure 5.7b). We note that an increase of the pumping rate corresponds to a decrease of Var[C] across the various stress periods. This behavior is consistent with the observation that an increase of the extraction rate enhances the likelihood that the contaminant is captured by the well. This consequently tends to dampen the variability of contaminant concentrations observed at the well, with a corresponding reduction of Var[C]. The opposite tends to happen when the pumping rate is reduced (and eventually goes to zero), thus leading to increased values of Var[C] at the well. By comparing the variance of the flux-averaged concentration detected at the pumping well operating under a temporally variable strategy for the three investigated scenarios (Figure 5.7b addresses scenarioIII and Figures SM6 and SM: of the SM respectively refer to scenariosI and II), we observe that the maximum values in Figure 5.7b (scenario III) are much higher than their counterparts appearing in Figure SM6 (scenario I) and Figure SM: (scenario II) of the SM. We recall that scenario III is associated with the highest temporal variability of the well operational scheme. Our results indicate that the uncertainty associated with the flux-averaged concentration is higher for pumping strategy III than for pumping regimes I and II. This behavior is also supported by viewing the results of Figures 5.8 and 5.9, respectively depicting the coefficient of variation (CV) of the flux-averaged concentration for the constant and transient pumping strategy of scenarioIII and comparing them against the corresponding results depicted 83 in the SM (see Figures SM;, SM32 for scenario I and Figures SM33 and SM34 for scenario II). When an extraction schedule with high temporal variability (see, e.g., scenario III) is Figure 5.6: Mean concentrationhCi observed at the pumping well for the uniform pumping strategy of Scenario III (a) and the transient pumping strategy of Scenario III (b). Figure 5.7: Variance of the concentration Var[C] observed at the pumping well for the uniform pumping strategy of Scenario III (a) and the transient pumping strategy of Scenario III (b). 84 applied at the well, the values of CV tend to increase, in particular during time intervals when pumping is stopped. We conclude that considering a high temporal variability of the pumping rate has the effect of increasing the uncertainty associated with the concentration observed at the pumping well. We observe that while it is well known that aquifer heterogeneity leads to uncertain contaminant transport prediction, our results indicate that the engineered operational pumping schedule, manifested throughQ w (t), may play a major role in controlling the variability of the moments of flux-averaged concentration at the well. Results in Figures 5.8, 5.9, SM;, SM32, SM33, SM34 also show that uncertainty in the solute concentration at the well is large at both early and late times. This implies that the (ensemble) mean is not sufficient to characterize the transport behavior and quantification of higher order moments is needed. 5.6.5 Probability of exceedance of concentration limit value Wecomputeheretheprobabilityoffailure, i.e. theprobabilitythatthecontaminantconcentration exceeds an established limit value,C ∗ . This probability represents a measure of risk, as defined in equation (5.7) and is computed on the basis of the MC realizations performed. The exceedance probability provided in equation (5.7) enables one to evaluate the probability of extreme events since it contains information about the mean, variance and higher order moments of the flux- averaged concentration. Figures 5.: and 5.; depict the temporal evolution of P (C ≥ C ∗ ) for the uniform and time-varying pumping strategies of scenario III. The choice of the operational scheme (constant or variable in time pumping) to withdraw a fixed water volume has a clear influence on the temporal evolution of risk (see equation 5.7) at the pumping well. The results depicted in Figure 5.: are linked to a temporal pattern which is similar to that of the mean flux-averaged concentration in Figure 5.6a. For the cases analyzed, the results suggest that higher mean concentration values are associated with increased probability to exceed the concentration threshold C ∗ . When considering transient pumping, we note that the mean contaminant concentration (Figure 5.6b) and the probability of exceedance (Figure 5.;) display 85 Figure 5.8: Coefficient of variationCV of the concentration observed at the pumping well for the uniform pumping strategy of Scenario III. Figure 5.9: Coefficient of variationCV of the concentration observed at the pumping well for the transient pumping strategy of Scenario III. 86 Figure 5.:: Probability of exceedance of the concentration threshold P (C≥C ∗ ) for the uniform pumping strategy of Scenario III. Figure 5.;: Probability of exceedance of the concentration thresholdP (C≥C ∗ ) for the transient pumping strategy of Scenario III. 87 similarpatternsforthelowestvaluesofthevarianceoftheflux-averagedconcentrationobservedat the pumping well (Figure 5.7b). This correspondence is not as evident for increased concentration variance (Figure 5.7b). As noted in chapter 5.6.4, highly variable transient pumping strategies favor an increase of the uncertainty of solute concentrations at the pumping well. This, in turn, affects risk uncertainty, as quantified by equation (5.7). In general, we observe that the transient operational scheme selected in scenarioIII induces a temporal oscillating behavior of risk at the well (Figure 5.;), according to which risk tends to decrease with the temporal reduction of the mean concentration signal. 5.6.6 Local scale dispersion and hydrogeological heterogeneity The previous chapters considered a fixed level of heterogeneity and local scale dispersion (see list of parameters in Table 5.3). We now assess the impact of aquifer heterogeneity and local scale dispersion on the statistics of the flux-averaged concentration in the presence of variable pumping rates. 5.6.6.3 Impact of local scale dispersion on the concentration statistics Figure 5.32 illustrates the influence of increasing local scale dispersion on the mean and on the variance of the flux-averaged concentration for the diverse pumping operations of scenario III. The mean concentration for constant and variable pumping rates are shown in Figures 5.32a and 5.32b, respectively, whereas Figures 5.32c and 5.32d depict the concentration variance for constant and variable pumping schedules, respectively. The effect of local scale dispersion on the concentration variance is manifested through the longitudinal Péclet number, defined as P e = λv 0 /D x . Figure 32 shows the results for both longitudinal P e = 800 (corresponding to the parameters listed in Table 3) and P e = 200 for a fixed level of heterogeneity (as rendered by σ Y 2 = 3). The Péclet number of 422 corresponds to increased local scale dispersion (longitudinal and transverse dispersivities are set to 0.04 m and 0.0004 m, respectively). The remaining input 88 parameters for the simulations related toP e = 200 are the same as those reported in Table 5.3 of the manuscript. As depicted in Figure 5.32, increasing the value of local scale dispersion (i.e., reducing P e ) leads to decreased values of the concentration variance (Figures 5.32c and 5.32d). This is related to the observation that dispersion smooths out the concentration gradients, hence reducing the sample-to-sample variability of transport observables at the well [6:,8;,329]. Note that the tran- sient operational scheme applied at the well has marked influence on the concentration variance for both Péclet numbers considered, as opposed to what can be observed for the constant in time pumping scheme. This is mainly attributed to the following reasons. First, when the pumping rate is temporally variable, the solute plume residence time in the aquifer increases, thus allowing for the effects of local scale dispersion to become more pronounced. Second, temporally variable flows, induced by the action of dynamic pumping rates, tend to enhance solute spreading, which in turn augments dilution (see also [69]). Furthermore, this dilution enhancement (driven by the action of variable pumping rates) leads to an increase of the plume spatial extent (in the longitu- dinal and transverse directions) and to an increase in the probability of the plume being captured by the well, thus contributing to a reduction of the variance. For the set of parameters used in this numerical investigation, the significance of local scale dispersion on the reduction of the con- centration variance is diminished for a constant pumping rate (compare Figure 5.32c and 5.32d). It is important to note that the joint effects of local scale dispersion and pumping operations on the concentration variance will also depend on the distance between the contaminant source and the receptor (i.e., the pumping well). With reference to the effect of local scale dispersion on the concentration mean at the extracting well (Figures 5.32a and 5.32b), we notice that the differences between the two values ofP e considered are smaller as compared to the corresponding differences detected in the temporal evolution of the concentration variance displayed in Figures 5.32c and5.32d. When we increase local scale dispersion, the mean concentration does not change much for a constant pumping schedule (Figure 5.32a) because the distance between the source 89 and the operating well is not sufficiently large in our scenario (approximately equal to 10 λ). Therefore, the effects of local scale dispersion are not that evident on the temporal evolution of the mean concentration. We recall that the advective time scale is the same for both values of P e . On the other hand, for the case of transient pumping (Figure 5.32b), the residence time of the plume in the aquifer increases, thus allowing for the effects ofP e to become slightly visible. A further increase of the travel distance of the contaminated plume will enhance the effects of local scale dispersion on the mean concentration for both pumping operations (constant and variable in time). Figure 5.32: Impact of the Péclet number (Pe) on the concentration mean for (a) constant and (b) time-varying pumping rate. Impact of the Péclet number (Pe) on the concentration variance for (c) constant and (d) time-varying pumping rate. Purple curves refer to Pe = 800 and blue curves refer to Pe = 200. 5.6.6.4 Impact of aquifer heterogeneity on the concentration statistics In the following, we investigate the sensitivity of the concentration statistics to the degree of aquifer heterogeneity in the presence of both constant and time-varying pumping schedule. The 8: level of aquifer heterogeneity is epitomized by the log-conductivity variance σ Y 2 . Figure 5.33 provides a comparison between the mean (Figures 5.33a and 5.33b) and variance (Figures 5.33c and 5.33d) of the flux-averaged concentration at the well for σ Y 2 = 0.5 and 3 under a constant and temporally dynamic pumping operations. The larger log-transmissivity variance (σ Y 2 = 3) Figure 5.33: Mean concentration for (a) constant and (b) time-varying pumping rates. Variance behavior is also depicted for (c) constant and (d) time-varying pumping schemes. Blue curves refer to low heterogeneity and purple curves to high heterogeneity. yields an increased likelihood of the occurrence of preferential flow paths within the aquifer which contribute to faster solute arrival. Solute migration through these high velocity paths tends to be locally dominated by advection and the solute reaches the pumping well at earlier times (Figures 5.33a and5.33b). The temporal evolution of the mean concentration displayed in Figures5.33a and 5.33b for both log-transmissivity variances also shows that an increased strength of σ 2 Y can also augmenttheprobabilityofoccurrenceoflowconductivityzonesthatcanentrapsolute, resultingin late arrival times. This contrast between low and high transmissivity values for largerσ Y 2 causes the plume to spread erratically and facilitates the exchange of solute mass between streamtubes 8; (i.e., dilution enhancement). In order words, higher permeability variability tends to reduce the mean solute concentration and smooths out concentration gradients, thus reducing concentration variance (Figures 5.33c and 5.33d). The numerical results in Figures 5.33c and 5.33d illustrate that the differences between the concentration variance obtained forσ Y 2 = 0.5 and 3 are higher when the pumping operation varies in time. The dynamic pumping scheme also increases the magnitude of the concentration mean and variance whenσ Y 2 is lower (Figures 5.33b and 5.33d). We remind Figure 5.34: Probability of exceedance of the concentration thresholdP (C≥C ∗ ) for (a) constant and (b) time-varying pumping rates. Blue curves refer to low heterogeneity and purple curves to high heterogeneity. Figure 5.35: Mean concentrationhCi observed at the pumping well for (a) constant and (b) time-varying pumping strategy of Scenario III for different levels of aquifer heterogeneity. that the results displayed in Figure 5.33 are not completely general because other factors can play a role in the uncertainty of the concentration statistics. These additional factors include the solute travel distance, source release conditions, location of the solute source, the conceptualization of the heterogeneity model and the presence of chemical reactions. Nevertheless, it is evident that 92 the interplay between aquifer heterogeneity and pumping operation plays fundamental role in probabilistic risk assessment. Figure 5.34 shows the temporal variability of the risk (probability of exceedance), defined in equation (5.7), for constant and dynamic pumping operations and different values of σ Y 2 . With the exception of the early and late time regimes, the results in Figure 5.34 show that the probability of observing a concentration value larger than a critical value is higher for low than for high heterogeneity. The reason for this behavior is similar to the discussionrelatedtoFigure5.33cand5.33d; i.e.,largervaluesofσ 2 Y facilitatemasstransferbetween neighboring streamtubes, a phenomenon which enhances dilution and, as a consequence, results in a larger probability of having lower concentration values. Therefore, the effects of the interplay between advection and local-scale dispersion become more evident for larger σ Y 2 , resulting in a tendency to lower the risks at the operating well (see Figure 5.34 for both pumping schemes). Otherwise, larger heterogeneity leads to higher risks at the early and late times. We highlight that the risk computed for an aquifer with high heterogeneity persists for longer periods of times (compare results obtained for σ Y 2 = 0.5 and σ Y 2 = 3). This is a consequence of the increased tailing effects of the concentration BTC observed in each MC realization when heterogeneity is large. For the numerical set-up adopted in this work, we show that the pumping operation scheme hasasignificantcontrollingroleonthetemporalevolutionoftherisk. AsdepictedinFigure5.34b, the temporal variability of the risk and its multiple peaks are controlled by the temporal patterns of the pumping rate. Albeit the focus of the current work lies on the interplay between aquifer heterogeneity and pumping operation (i.e., engineering factors) on the uncertainty of the concentration BTC, we include results in the SM for a homogeneous aquifer (i.e., constant transmissivity field). Figure SM43showshowthepumpingoperationcanaffecttheasymmetryoftheconcentrationBTCunder uniform flow fields. Furthermore, as shown in Figures SM42 and SM43, the peak concentration in the deterministic homogeneous scenario is higher for both pumping operations when compared to the heterogeneous aquifers. Figure 5.35 depicts the temporal evolution of the concentration mean 93 for different values of aquifer heterogeneity (σ Y 2 = 0.5, σ Y 2 = 3 and for the homogeneous case). The increase of the mean concentration with lower heterogeneity is linked to the selected location of the solute source zone (recall that, in the settings we analyze, the centroid of the source zone is aligned with the location of the well). Plume meandering is less pronounced for lower heterogeneous aquifers (when compared to higher heterogeneity in Y), therefore the probability of the plume being captured by the well is higher for low values of σ Y 2 in our settings. Figures SM42 and SM43 depict the suite of results for the case of a deterministically homogeneous aquifer (σ Y 2 = 0). 5.7 Conclusions We study and compare the effect of temporally variable and uniform pumping regimes on key featuresofcontaminanttransportinarandomlyheterogeneousaquifer. Thischapterconsidersthe joint effects of spatially heterogeneous hydraulic conductivity (or transmissivity) and temporally varying well pumping rates and offers some insights on a realistic approach for the evaluation of the risk associated with contamination of groundwater extracted at the pumping well location. The analysis is performed within a stochastic framework upon relying on the numerical study of three distinct pumping scenarios. Thetwoleadingstatisticalmomentsoftheflux-averagedcontaminantconcentrationsrecovered at the well are computed. We document the way the temporally dynamic pumping rate augments the uncertainty in the flux-averaged concentration at the well. We then provide an appraisal of risk at the well by computing the probability that the contaminant concentration exceeds a defined threshold value (i.e., probability of failure). Our work leads to the following major conclusions. 3. In addition to aquifer heterogeneity, the use of a transient pumping strategy at the well can markedly affect the temporal evolution of contaminant concentration BTCs and statistical 94 moments. We show that dynamic pumping strategies induce multiple peaks in the mean and variance of the flux-averaged concentration detected at the extraction well. 4. Lowest and largest values of mean and variance of flux-averaged concentration at the well tend to occur at the same time. This observation is supported by our numerical findings and by analytical results illustrated in Appendix A. 5. The activation of the pumping well (or the increase of its extraction rate) induces contam- inant dilution with fresh water at the well. As a consequence, the detected contaminant concentration tends to be reduced. 6. The choice of the type of engineering control to the temporal sequence of well pumping rates could represent a key factor in quantifying the uncertainty of the contaminant concentra- tion detected at the well. This observation is supported by considering that uncertainty associated with detected BTCs at the well increases for highly variable pumping regimes (compare scenario III against scenarios I and II). We show that this controlling role of the temporal extraction schedule on uncertainty has direct consequences to risk analysis. The selection of a dynamic pumping regime has a clear influence on the temporal evolution of risk, as defined in 5.7, at the well, i.e., pumping rate fluctuations induce a temporally oscillating risk pattern. The mean flux-averaged solute concentration and the probability of exceedance of a given threshold value show a similar temporal evolution when the variance of flux-averaged concentrations at the pumping well is small. In these cases, larger mean concentration values correspond to larger probability of exceedance of a given concentra- tion threshold. Otherwise, this correspondence was not detected in our cases for increased concentration variance. Our findings suggest that risk analysis should incorporate dynamic pumping rates to reflect realis- tic operational practices and be able to provide important information to risk managers, including realistic quantifications of the uncertainty associated with risk at the pumping well. 95 Note that the results and conclusions presented in this study are confined to a hydrogeological setting whose randomly heterogeneous conductivity is characterized by given sets of parameters. Theresultscanbeimpactedbyotherfactors, including, e.g., thedistancebetweenthecontaminant sourceandthereceptor, thedimensionsofthesourcezone, massreleaseconditions, theratioofthe baseflowdischargewithrespecttothepumpingrate, chemicalreactionsandtheconceptualization of the hydrogeological heterogeneity. The influence of variable pumping rates on probabilistic risk analysis depends on the interactions between temporal fluctuations induced by the action of pumping, local scale dispersion, and spatial heterogeneity. Quantifying the effects of the overall spectrum of scenarios is complex due to the interplay of three distinct time scales defined by the pumping well operation, advection and local scale dispersion. The relative importance of the temporally variable pumping rate on risk will depend on the ratio between these characteristic time scales and should be subject of further investigation. As a future projection of the study, one can incorporate in our methodological framework adverse human health effects [4,66] and other well vulnerability criteria, including, e.g., the time taken to breach a certain quality objective (e.g., drinking standard) at the well, the total time of well failure (i.e., non-compliance with a quality objective), and the time taken to recover the well from failure. In the next chapter we explore the importance of accounting for realistic groundwater ex- traction in groundwater contamination predictions, for different representations of the natural variability of K. 96 Chapter 6 Solute concentration at a well in non-Gaussian aquifers under constant and time-varying pumping schedule The content of this chapter is published in Journal of Contaminant Hydrology [345]. 6.3 Introduction Probabilisticcharacterizationsofsubsurfacecontaminanttransportwithinwellfieldsunderpracti- cal operational conditions are cornerstones of modern groundwater management and risk analysis best practices. Modeling set-up and results rely on an appropriate assessment of multiple factors, including the spatial variability of hydrogeological variables and the adopted engineering controls (e.g., groundwater pumping and/or injection operations). In this study, we aim at exploring the feedback between transient versus uniform pumping schedules and Gaussian or non-Gaussian na- ture of the heterogeneous transmissivity (T) field in influencing contaminant concentration and its uncertainty at extraction wells. Therefore, this chapter analyzes the combined effects of (a) structural heterogeneity of properties characterizing geological media and (b) planned sequences of pumping cycles on solute concentrations recovered at pumping wells, including an appraisal of the corresponding uncertainties. 97 A sequence of predefined pumping intervals is typically scheduled by water management agen- cies to achieve a trade-off between maximization of the benefits to anthropogenic activities and minimization of the environmental footprint of the groundwater withdrawal process. Nonethe- less, most studies focusing on probabilistic analyses of subsurface contaminant transport within well fields are limited to scenarios associated with constant extraction practices. For example, a number of studies [49,88,95,:9,338,38;,436,437] employ numerical Monte Carlo simulations to quantify the uncertainty in the extension of well protection areas within randomly heterogeneous aquifers, under steady-state background groundwater flow and in the presence of a single or mul- tiple wells operating at a constant rate. Additional examples of studies that consider the way solute transport is driven by one or more wells pumping at a constant extraction/injection rate include [5;,325,37:–382,392,3;8]. Despite its importance, the impact of transient pumping on solute transport has received much lessattention[44,45,33:,438]. Liberaetal.[344]systematicallyinvestigatetheeffectsofscheduling of pumping operations (i.e., considering transient vs. constant in time pumping) on contaminant solute breakthrough curves (BTCs) detected at the pumping well in a spatially heterogeneous multi-Gaussian log-conductivity field. These authors show that a transient pumping strategy can markedly affect the temporal pattern of BTCs. The results of this study elucidate the importance of the pumping sequence on uncertainty quantification of solute transport, risk analysis and contaminated site management. Transient flow effects on the delineation of wellhead protection areas (WHPA) or capture zones have been considered in a series of works. For instance, [386] highlight the importance of considering seasonal variations in pumping operations for WHPA delineation. Reilly and Pullock [388] observe that transient flow conditions should be considered to properly characterize transport of solutes released near the boundary of the well. Festger and Walter [87] and Jacobson et al. [326] evaluate the effects on capture zones of temporal variations in the direction of the hydraulic gradient. Jacobson et al. [326] illustrate how the uncertainty in the magnitude and direction of the mean regional flow influences the extent of time-dependent 98 capture zones. Neupauer et al. [369] study chaotic advection of a solute in an aquifer where a transient flow field is induced by injection and extraction sequences (see also [385]. In addition to pumping activities, the ubiquitous heterogeneity of hydraulic conductivity, K, or transmissivity T, is known to affect solute transport [55]. High costs associated with site characterization contribute to hamper exhaustive reproductions of K-fields. This contributes to uncertainty associated with the quantitative description of transport scenarios. A common procedure adopted in stochastic hydrogeology is to assume a multivariate Gaussian distribution forthespatialfieldoflog-conductivity,Y = ln(K). Spatialvariabilityofhydrogeologicalattributes is, however, known to be more complex than described by a Gaussian model [9:,:5,;6,324,355, 445,448,452]. In this context, high-resolution data analysis performed at the Macrodispersion Experiment (MADE) site, at the Columbus Air Force Base (Mississippi, USA), indicate that highly heterogeneous aquifers could be characterized in terms of non-Gaussian Y fields [358]. Fogg et al. [94] study alluvial heterogeneity in the Livermore Valley (California, USA) and question the ability of the multi-Gaussian assumption for Y to render an adequate description of the investigated area. Rubin and Journel [3:7] point out that multi-Gaussian models can potentially fail in representing connected paths of extreme permeability values which might take place in the subsurface. This is also observed in other studies [328,38:,3:9] and the influence of such features on contaminant transport, risk assessment and groundwater remediation strategies has been highlighted in several works [59,5;,327,3;:, and references therein]. Methods to generate syntheticrandomfieldsthatreflectaspectsofhydrogeologicstructureand/orarchitecture, someof which may render such fields non-Gaussian, have been proposed (e.g., [7:] and references therein). A critical element which is emerging from a variety of studies is that the assumption of Gaus- sianity forY is not consistent with features displayed by the sample probability distribution (and main statistical moments) of increments ΔY (s) = Y (x)−Y (y) between two vector locations x andy (s =x−y, denoting separation scale or lag). A common manifestation of this phenomenon 99 is that while frequency distributions of Y often exhibit mild peaks and light tails, those of in- crements ΔY (s) are typically symmetric with peaks that grow sharper, and tails that become heavier, as s =ksk decreases [34:,359,375,397,399]. Hydrogeologic variables that have been shown to exhibit such behaviors include log-permeabilities of porous and fractured geologic me- dia [34:,375,397,399,3;7], neutron porosities in deep boreholes [::], and soil composition data and hydraulic parameter estimates [:;,;2] in a deep vadose zone. Manifestations of similar sta- tistical scaling of a variety of Earth, environmental, ecological, biological, physical, astrophysical, and financial variables are reported, among others, by [367]. As stated above, these features are clearly non compatible with a description of Y which is based on a Gaussian distribution model. Painter [375] proposes to adopt L´ evy-stable distributions to characterize permeability hetero- geneity. Strebelle [425] proposes an algorithm that utilizes multiple-point statistics inferred from training images to model conductivity fields. Linde et al. [348] tackle the same problem through the use of training images based on hydrogeological facies mapping of nearby outcrops. Haslauer et al. [;6] employ non-Gaussian copula-based K models to study solute macrodispersion. A sta- tistical framework that captures the disparate, scale-dependent distributions of Y and ΔY in a unified and consistent manner is offered by relying on the Generalized sub-Gaussian (GSG) model introduced by [398,39:]. In this context, Riva et al. [394] explore analytically lead-order effects that non-Gaussian heterogeneity described by the GSG model has on the stochastic description of flow and transport under mean uniform steady-state flow in an unbounded, two-dimensional domain. In light of the above, it is relevant to ask the following question: what is the feedback between transient pumping operations and the non-Gaussian nature of Y in influencing solute concentra- tion and its uncertainty at the extraction well? Key research and operational questions driving our study are the following: how important are non-Gaussian conductivity features of the kind revealed when considering consistency between distributions of Y and its increments in transport 9: when pumping is in operation? Does the pumping scheme overshadow the significance of non- Gaussian Y fields on statistics of solute concentration? Through a suite of computational studies, this work also enables to analyze the relative impact of the conductivity structure (i.e., Gaussian vs. non-Gaussian) on the solute concentration observed at the operating well in the presence of constant and time-varying pumping rates. 6.4 Problem Formulation We consider a fully saturated two-dimensional (2D) confined porous formation identified by a Cartesian coordinate system, with vector location indicated by x = (x,y). A uniform-in-the- mean base flow q 0 takes place along the x-direction, and the transmissivity field, T, is spatially heterogeneous. The porous formation porosity and storativity, respectively denoted as φ and S, are considered to be constant. A pumping well is operating with an extraction rateQ w (t), with t indicating time, at locationx w = (x w ,y w ). Figure 6.3 illustrates a sketch of the setting analyzed. Flow within the hydrogeological system is governed by: S ∂h(x,t) ∂t =∇· [T (x)∇h(x,t)] +Q w (t)δ(x−x w ), (6.3) where h denotes hydraulic head and δ represents the Dirac delta function. West (left) and east (right) boundaries of the porous formation (see Figure 6.3) are characterized by fixed hydraulic head values and no-flow boundary conditions are prescribed on the north and south boundaries of the domain. A hazardous non-reactive solute is instantaneously released at time t 0 within an area A 0 (see Figure 6.3) with constant concentration c 0 . The spatio-temporal evolution of the solute plume is assumed to be governed by the advection-dispersion equation: ∂C(x,t) ∂t −∇· [D∇C(x,t)−v(x,t)C(x,t)] = 0; c(x,t 0 ) =c 0 for x ∈A 0 . (6.4) 9; Figure 6.3: Sketch of the problem analyzed. where C indicates solute concentration and v is the Darcy-scale velocity (i.e., v = q/φ, q indi- cating specific discharge). Local-scale dispersion is given by the tensor D, with components D x and D y , respectively along the x and y-direction (Figure 6.3). Note that Q w in equation (6.3) is negative if extraction occurs and positive in case of injection. In our analysis, the prescribed head and flow boundary conditions are located sufficiently far away from the solute transport area to avoid boundary effects [3:5]. Inthefollowing, weinvestigatethewayconcentrationbreakthroughcurves(BTCs)areaffected by constant and transient pumping conditions. This analysis is carried out for Gaussian and/or non-Gaussian spatially random log-transmissivity fields as described in chapter 6.5. :2 Symbol Significance Units Values L sx , L sy Aquifer Size m 170, 150 Δ x , Δ y Grid Size m 1/8I, 1/8I J Mean head gradient – 0.59% φ Porosity – 0.2 α x Longitudinal dispersivity m 0.01 α y Transversal dispersivity m 0.0001 x w , y w Location of pumping well m 114.5, 74.5 Q w Constant well pumping rate m 3 /d 0.3 Q w Variable well pumping rate m 3 /d 0.8/0/0.4/0/0.8/0/0.4/ 0/0.8/0/0.4/0/0.8/0/0.4/0 V w Volume of pumped water m 3 1440 C ∗ Concentration threshold g/m 3 10 P e Péclet number – 800 Table 6.3: Main parameters employed in the study. 6.5 Methodology 6.5.3 Domain configuration and numerical implementation We solve equations (6.3) and (6.4) within a 2D system of size of Ω =L sx ×L sy , whereL sx = 170 m, L sy = 150 m (see Fig. 6.3). Values of the main parameters adopted in this study are listed in Table 6.3. The values listed in Table 6.3 were selected for the purpose of illustration, all computational results being presented in dimensionless form. In agreement with the main results of previous works [33;,363], the domain is discretized by a uniform grid formed of square elements of size Δx = Δy = 1/8 I, with I being the characteristic integral scale of the random log-transmissivity, Y = ln(T ), field. Details of the models adopted to describe the random nature of Y are provided in chapter 6.5.4. The solute is instantaneously released over an aerial source zone A 0 of size L sx ×L sy (see Figure 6.3). The barycenter of A 0 is located at distance L (measured along the x-direction) from an operating pumping well. Following the work of Libera et al. [344], which was based on an analysis of pumping strate- gies employed by groundwater management to satisfy diverse societal needs (e.g., drinking and irrigation), we adopt the pumping operation shown in Figure 6.4. :3 0 1000 2000 3000 4000 t[d] 0 0.2 0.4 0.6 0.8 1 Q w [m 3 /d] S I (a) 0 1000 2000 3000 4000 t[d] 0 0.2 0.4 0.6 0.8 1 Q w [m 3 /d] S II (b) Figure 6.4: Pumping strategies (a) constant flow rate S I and (b) variable flow rate S II . A constant in time extraction strategy, here denoted byS I , is depicted in Figure 6.4a. Figure 6.4b depicts the pattern of a withdrawal strategy that varies in time according to a predefined sequence, indicated byS II . Note that, as explained in chapter 6.3, most available studies refer to constant in time pumping scenarios that may not accurately represent realistic operations. The selected strategies (i.e., S I and S II ) are characterized by the extraction of the same volume of groundwater across the simulation time. For a given Y field, we employ the well-tested codes MODFLOW [;4] and MT3DMS [456], respectively to solve the transient groundwater flow equation (6.3) and the advection-dispersion equation (6.4). Note that the groundwater flow equation (6.3) and the numerical model (i.e., MODFLOW) adopted do not consider possible additional effects (e.g., non-Darcian flow, skin ef- fects, or storage) at the well. Solute transport is computed through the Method of Characteristics. To quantify uncertainty in the concentration C at the well, we employ a numerical Monte Carlo :4 (MC) framework. Our analysis is based on a set of 10, 000 MC simulations for each of the inves- tigated scenarios. The choice of the size of the MC sample is based on a statistical convergence analysis (details not shown). 6.5.4 Random Y Model As mentioned in chapter 6.3, a variety of works emphasize the importance of adopting non- Gaussian Y fields in solute transport studies. In our analysis, we employ the GSG model intro- duced by [398,39:] to generate multiple realizations of Y. As stated in the Introduction, this model has the key ability to embed in a unique theoretical framework the scale-dependent non- Gaussian features of the main statistics ofY and its increments taken at diverse lags, which have been documented in a variety of studies. The GSG model is here only briefly summarized for the sake of completeness, additional details being provided by [398,39:]. We write zero-mean random fluctuations, Y 0 (x) =Y (x)−hYi, as Y 0 (x) =U(x)G(x) (6.5) where denotes ensemble mean (expectation), G(x) is a single- or multi-scale Gaussian random field andU(x) is a non-negative subordinator independent ofG(x). The subordinatorU consists of statistically independent identically distributed (iid) non-negative random values at all points x. We consider the subordinator U to be log-normally distributed (other choices being possible) according toU≡ lnN(0, (2−α) 2 ), whereα< 2. Note that in this case the PDF ofY 0 in equation (6.5) coincides with the classical normal-lognormal model (see, e.g., [::] and references therein). To investigate the effects of Gaussian versus non-Gaussian Y fields on solute transport we consider two forms of Y: sub-Gaussian (see equation (6.5)), Y SG , and Gaussian, Y G . The latter is obtained from the former when α→ 2. For the purpose of comparison, we consider these two forms of the Y field to have equal mean values, i.e. hYi =hY SG i =hY G i = 0, variances :5 σ 2 Y =σ 2 Y G =σ 2 Y SG and integral scales I =I Y G =I Y SG . The Gaussian random function, G(x), in equation (6.5) constitutes a truncated fractional Brownian motion (tfBm) with truncated power variogram [74]: γ 2 G (s) =γ 2 (s;λ u )−γ 2 (s;λ l ), (6.6) where γ 2 (s;λ m ) = Aλ 2H m 2H " 1−exp − s λ m + s λ m 2H Γ 1− 2H, s λ m # , m =l,u. (6.7) Quantities λ u and λ l in equations (6.6-6.7) are the lower and upper cutoff scales of the vari- ogram model, respectively proportional to the length scales of data support and domain size; H is the Hurst coefficient and A is a constant. This choice of variogram model is consistent with documented scaling phenomena, including power-law scaling of sample structure functions (in- cluding the variogram ofY) in midranges of lags and nonlinear scaling of power-law exponent with order of sample structure function [::–;3,377,397,399,3;6,3;7]. As such, equation (6.6) allows bridging across scales by analyzing jointly data characterized by diverse support/measurement scales across windows (observation domains) of diverse size at a site [368]. We simulate three diverse collections (ensembles) of non-GaussianY SG fields (each constituted by 10, 000 realizations) distinguished by three values of the parameterα = 1.2, 1.5, 1.8, represent- ing strong to relatively mild departure from a Gaussian behavior) and one set of Gaussian Y G fields (also with 10, 000 realizations). Input parameters used to generate Y G and Y SG fields are listed in Table 6.4 (input values employed to generateY G are listed in the column labeledα→ 2). Each Gaussian and non-Gaussian collection/ensemble of log-transmissivity realizations is cou- pled with both pumping strategies adopted (S I andS II , see Figure 6.4) within the numerical MC framework. Therefore, we perform a total number of 80, 000 MC simulations. :6 Symbol Significance Units Values α=1.2 α=1.5 α=1.8 α→ 2 I Integral scale of Y m 8.00 8.00 8.00 8.00 σ 2 Y Variance of Y – 3.00 3.00 3.00 3.00 A Constant – 0.0587 0.1753 0.3185 0.3575 H Hurst coefficient – 0.33 0.33 0.33 0.33 λ u Upper cutoff scale m 34.58 22.67 17.98 17.20 λ l Lower cutoff scale m 1.00 1.00 1.00 1.00 Table 6.4: Parameters characterizing the transmissivity field. 6.6 Results and discussion Thischapterisstructuredasfollows: westartbypresentingthetemporalevolutionoftheloworder statistics (mean and variance) of the contaminant concentration, C, recovered at an observation well placed, for simplicity, at the same location of the pumping well. Then, we focus on the analysis of the peak concentration, C p , observed at the well. The latter represents an important quantity for the management and remediation of polluted areas and can be used as proxy for dilution [89]. All analyses are performed for the two pumping strategies, namely S I (constant in time) and S II (variable in time), illustrated in Figure 6.4. Diverse values of the parameter α are considered (see chapter 6.5.4) for each pumping scenario. We present all results in dimensionless form. The concentration is normalized by C ∗ , which represents a contaminant concentration threshold, as established, for instance, by environmental protection agencies (e.g., EPA). Here we set, without loss of generality to our methodological approach,C ∗ = 10 g/m 3 . The selected value of the critical concentration is in line with the EPA’ s Maximum Contaminant Level (MCL) for nitrate [42;], and is employed for the purpose of illustration in our study. Time is normalized by I/v 0 , with v 0 = q 0 φ, q 0 = T G J and T G = exp(hYi). The pumping rate Q w is normalized by Q 0 , where Q 0 =q 0 L y represents the uniform-in-the-average water flow discharge. ThelongitudinalPécletnumber,P e =Iv 0 /D x , issetto:22(seeTable6.3). Parameter values employed in this synthetic analysis allow illustrating the interplay between pumping well operations and natural heterogeneity of the porous formation. :7 6.6.3 Temporal evolution of low-order moments of solute concentration at the well Figure 6.5 depicts the temporal Monte Carlo based pattern of the dimensionless mean of the solute concentration,hCi/C ∗ , observed at the well location. 0 5 10 15 tv 0 I [−] 0 0.2 0.4 0.6 0.8 1 1.2 1.4 hCi C ∗ [−] (a) S I 0 0.2 0.4 0.6 0.8 1 Q w Q 0 [−] α = 1.2 α = 1.5 α = 1.8 α→ 2 Qw/Q0 0 5 10 15 tv 0 I [−] 0 0.2 0.4 0.6 0.8 1 1.2 1.4 hCi C ∗ [−] (b) S II 0 0.2 0.4 0.6 0.8 1 Q w Q 0 [−] α = 1.2 α = 1.5 α = 1.8 α→ 2 Qw/Q0 Figure 6.5: Temporal evolution of the normalized mean concentrationhCi observed at the pump- ing well for three values of α = 1.2, 1.5, 1.8 and for a Gaussian Y field (α→ 2). Results for pumping strategy (a) S I and (b) S II . Results depicted in Figure 6.5a, b respectively refer to the pumping scheme S I (constant withdrawal) and S II (time-dependent withdrawal). Each mean concentration BTC presented in Figures 6.5a and 6.5b refers to a given value of α. Regardless the value of α, the mean concentration BTC displays a unimodal behavior under conditions associated with S I and a multimodal pattern when the scheme S II is active. The multimodal pattern observed in Figure 6.5b descends from the observation that a time-dependent pumping rate (S II ) induces temporal oscillations of the mean contaminant BTC. The impact of variable pumping rate on concentration :8 statistics was analyzed by [344] for Gaussian Y fields. Figure 6.5b suggests that solute dilution is induced during the time periods when pumping is active. Increased pumping rates would lead to an enlargement of the catchment region (e.g., [35]). Hence, an increased volume of clean water would be captured (on average) at the well together with the solute plume. Therefore, the mean concentration decreases within these intervals because of the mixing of polluted water with clean water. We also note that the largest values ofhCi occur within the strongly non-Gaussian Y field characterized by α = 1.2 for both scenarios S I and S II . As shown in Figure 6.5, the maximum value ofhCi decreases asα increases towards 2, i.e. as theY field tends to be Gaussian. We note that the differences inhCi between the Gaussian case (α→ 2) and the GSG setting characterized by α = 1.8 are negligible. Based on the results of Figure 6.5, we conclude that the pumping regime (S I orS II ) controls the pattern of the temporal evolution of the mean contaminant BTC, regardless of the value of α. The general shape of the temporal evolution ofhCi is, in fact, very similar for all values of α considered, as observed in both Figures 6.5a and 6.5b. Figures6.6aand6.6bdepictthedimensionlessvarianceofthesoluteconcentration,Var[C]/(C ∗ ) 2 , versus dimensionless time for all values of α analyzed, respectively forS I and S II . The pumping strategy clearly controls the general pattern of the temporal evolution of Var[C], regardless the value of α, i.e., a unimodal pattern of Var[C] is identified for constant pumping (Figure 6.6a) while a multimodal behavior of Var[C] is induced by temporal variability in the pumping rate (Figure 6.6b). We note that Var[C] in Figure 6.6b decreases when pumping is active, indicating that the variability across the MC ensemble is smaller when the solute is attracted to the well by the start of pumping. In this situation, the likelihood that the solute plume is captured by the well increases and the variability of concentration values at the well decreases. We then observe thatVar[C] increases approximately by an order of magnitude under regimeS II (compare Figure 6.6b and Figure 6.6a), these results being in line with the conclusions of [344]. On these bases, one can see that the choice of the pumping extraction operation (e.g., constant in time S I , as considered in most literature works, versus transient S II (which is more realistic) has a key role :9 0 5 10 15 tv 0 I [−] 0 5 10 15 20 25 30 Var[C] C ∗2 [−] (a) S I 0 0.2 0.4 0.6 0.8 1 Q w Q 0 [−] α = 1.2 α = 1.5 α = 1.8 α→ 2 Qw/Q0 0 5 10 15 tv 0 I [−] 0 5 10 15 20 25 30 Var[C] C ∗2 [−] (b) S II 0 0.2 0.4 0.6 0.8 1 Q w Q 0 [−] α = 1.2 α = 1.5 α = 1.8 α→ 2 Qw/Q0 Figure 6.6: Temporal evolution of the normalized concentration varianceVar[C] observed at the pumping well for three values ofα = 1.2, 1, 5, 1.8 and for a GaussianY field (α→ 2). Results for pumping strategy (a) S I and (b) S II . :: on quantification of the uncertainty associated with the concentration at the well. Amongst other factors, this is also related to the observation that a constant pumping scheme always controls the same portion of the flow field at all times. Otherwise, a transient pumping schedule enables to extend the influence of the well to diverse portions of the heterogeneous system, depending on time. As such, the effect of the aquifer heterogeneous structure plays an enhanced role under pumping scenarioS II than in the presence ofS I , resulting in an overall increase ofVar[C] at the well. We then observe that the increase of the concentration variance is magnified for the lowest values of α, i.e., as the departure of the GSG fields from Gaussianity increases. We remark that the differences in Var[C] between the Gaussian case (α→ 2) and the setting characterized by α = 1.8 are negligible, similar to what we observed forhCi. In summary, the analysis of Figures 6.5 and 6.6 lead to the conclusion that the pumping scheme selection (S I or S II ) clearly influences the temporal pattern (unimodal or multimodal) of the pollutant BTCs lead-order statistics (as expressed by mean and variance). The actual magnitude of the first and second moment of C is controlled by both the pumping scheme and the structural representation (Gaussian or non-Gaussian) of Y. 6.6.4 Statistical analysis of the peak concentration Here we analyze key statistical features of the peak value of the solute concentration,C p , observed at the well. This quantity is an important environmental performance metric (EPM) for risk analysis [5:] and, as previously stated, can also be used as a proxy for dilution [89]. 6.6.5 Outliers in the peak concentration distribution We present the box plots of C p in Figure 6.7 for pumping operational setting S I and for the diverse values ofα analyzed. Figure 6.8 depicts corresponding results for scenarioS II . We recall that the thickness of the box plots corresponds to the lag between the first and third quartiles of the probability distribution. Close inspection of Figures 6.7 and 6.8 evidences a considerable :; Figure 6.7: Pumping strategy S I : Box plots of C p for three values of α = 1.2, 1, 5, 1.8 and for a Gaussian Y field (α→ 2). Figure 6.8: Sketch of the problem analyzed. ;2 number of outside values (or outliers), identified in red and corresponding to the observations that fall outside the whiskers (represented by horizontal segments connected through dashed lines to the boxplots of Figures 6.7 and 6.8) under the action of both extraction schemes. Note that the upper whisker corresponds to the largest value observed given that the length of the dashed lines is 1.5 times the interquartile range, the same criteria applies to the lower whisker. Comparing Figures 6.7 and 6.8 suggests that the range of C p values is broadest when S II is active. In this case the largest values of normalized C p are roughly three times larger than the corresponding extreme values of C p computed under scheme S I . We also note that the number of such high values generally tends to decrease as the nature of the underlying log-transmissivity field tends to Gaussian (i.e., increasing α) for both pumping scenarios. This evidence suggests that the nature of the heterogeneous structure of Y influences the distribution of C p . These observations are consistent with the fact that a non-GaussianY structure increases the likelihood oftheoccurrenceofwell-connectedzonesoflowandhighconductivity(asmanifestedthroughhigh peaks of increment PDFs at short lags, whose effects are increasingly pronounced with departure from the Gaussian behavior). This specific feature is allowed to emerge in a stronger way in the presence of transient pumping than for constant extraction as already noticed in chapter 6.6.3 for the mean and variance of C. 6.6.6 Probability density function of the maximum concentration Figure 6.9 depicts the sample PDF, p(C p ), of the peak concentration detected at the well for pumping scenario S I and all values of α investigated. For completeness, the sample cumulative distribution function (CDF) of C p , P (C p ), is also depicted. Corresponding plots for scenario S II are included in Figure 6.:. The action of the transient pumping regime S II contributes to distribute the observed values ofC p across a wider range than that documented forS I (compare Figures 6.9 and 6.:). As such, the PDF of C p for S II generally encompasses a broader range of values and is characterized by ;3 0 2 4 6 8 1 0 1 2 1 4 1 6 1 0 -4 1 0 -3 1 0 -2 1 0 -1 1 0 0 1 0 0 1 0 1 0 .0 0 .2 0 .4 0 .6 0 .8 1 .0 C p / C * P (C p ) C p / C * S I p (C p ) = 1 .2 = 1 .5 = 1 .8 → 2 Figure 6.9: Pumping strategy S I : Peak concentration PDF p(C p ) for three values of α = 1.2, 1, 5, 1.8 and for a Gaussian Y field (α→ 2). See inset for peak concentration CDF P (C p ). 0 5 1 0 1 5 2 0 2 5 3 0 3 5 4 0 1 0 -4 1 0 -3 1 0 -2 1 0 -1 1 0 0 1 0 0 1 0 1 1 0 2 0 .0 0 .2 0 .4 0 .6 0 .8 1 .0 C p / C * P (C p ) S II = 1 .2 = 1 .5 = 1 .8 → 2 p (C p ) C p / C * Figure 6.:: Pumping strategy S II : Peak concentration PDF p(C p ) for three values of α = 1.2, 1, 5, 1.8 and for a Gaussian Y field (α→ 2). See inset for peak concentration CDF P (C p ). ;4 longer positive tails than the PDF of C p resulting from S I . These observations are consistent with the results depicted in the inset plots of Figures 6.7 and 6.8. The positive tails of the PDFs are quantitatively affected by the parameter α. Non-Gaussian Y fields are characterized by an increased probability of observing higher C p at the well, when compared to Gaussian Y fields (α→ 2), thus yielding enhanced tailing for p(C p ). This behavior is consistent with our earlier observationsaccordingtowhichtheGSGnatureoftheY fieldleadstoanincreasedlikelihoodthat solute can be conveyed through connected paths of high conductivity, thus yielding an increased tailing in the PDFs ofC p (i.e., higherC p values). Similar to what observed in chapter6.6.3, results for the Gaussian Y field virtually coincide with those obtained for α = 1.8. This result further emphasizes the challenges of distinguishing between these types of fields (for relatively large values of α) solely on the basis of system responses (e.g., in this cases, concentrations detected at the well). One can note a bimodal shape for p(C p ) in both Figures 6.9 and 6.:. This feature can be attributedtotheobservationthatverylowornoconcentrationsignalsareobservedatthepumping well (i.e., the solute plume does not hit the well) across some MC realizations, a significant portion of the plume being captured in other MC realizations (e.g., [38]). Note that while the observation thatp(C p )tendstobebimodalinthepresenceofpumpingwellsisasignificantresult, thisbimodal pattern in the PDF can change in the presence of other factors, including, e.g., increased travel distance between contaminant source and operating well and change of Péclet number. While of definite interest, these analyses are outside the scope of our current contribution. 6.6.7 Average of the maximum concentration We quantify here the impact of the pumping operation schedule and of the Y field structure on the average value of the peak concentrationhC p i. We do so by computing the relative change ;5 1.2 1.5 1.8 2 α[−] 0 0.1 0.2 0.3 0.4 0.5 η[−] S I S II Figure 6.;: Relative impact of non-Gaussianity on the mean peak concentration C p measured through η, see equation 6.8, for three values of α = 1.2, 1, 5, 1.8 and pumping strategy S I (light grey) and S II (dark grey). ofhC p i obtained across the collection of Y SG MC realizations with respect to the corresponding result associated with a Gaussian (Y G ) field as: η = hC p (t;α,S i i−hC p (t;α→ 2,S i )i hC p (t;α→ 2,S i )i , i =I,II (6.8) Figure 6.; depicts η versus α for scenario S I (light grey) and S II (dark grey). These results suggest thathC p i is more sensitive to the value of Îś when the spatially heterogeneous flow field is stressed according to scheme S II than it does for S I . As shown in Figure 6.;, the magnitude of η is larger for scenario S II (dark grey) and decreases as α increases, i.e. transitioning from a GSG to a Gaussian Y field. The response of the system due to a Gaussian Y field (α→ 2) is virtually indistinguishable from that associated with values of α≥ 1.8. 6.7 Conclusions This chapter investigates the impact of the model employed to describe the random spatial het- erogeneity of the aquifer log-transmissivity field (Y) on the statistics of the solute concentration (C) at a pumping well in the presence of two distinct pumping regimes. We consider a Gaussian ;6 and a Generalized sub-Gaussian (GSG, see equation (6.5)) model to describe the randomly het- erogeneous Y field. In the following, we briefly summarize the key conclusions emerged from the analysis. Thepumpingschemeinfluencestheshapeofthetemporalevolutionofmean,hCi,andvariance, Var[C], of C whereas the choice of the Y structural representation, as quantified by the value of α in this study, controls their magnitude. Transient pumping produces a multimodal behavior (whereas constant pumping results in a unimodal pattern) ofhCi and Var[C]. The multimodal behavior ofhCi is characterized by its lowest values taking place during stress periods of water pumping, the decrease inhCi being due to contaminant dilution with fresh water. We also observe thatVar[C] tendstodecreaseduringpumping timeintervals. This behaviorisrelated totheeffect of the well operation which increases the likelihood that the plume is captured by attracting water and hence results in a decreased variability of C at the well. Values of Var[C] are roughly one order of magnitude larger under transient pumping than in the presence of constant extraction. As such, engineering control (as manifested through selected pumping schedules) plays a marked role in the uncertainty associated with C. The highest values of solute peak concentration, C p , are prone to be observed at a well op- erating according to a time-varying schedule. This feature is amplified in the presence of values of α associated with an increased departure of the GSG Y field from a Gaussian behavior. The PDF of C p , p(C p ), is characterized by a bimodal shape for all cases analyzed in this study. Our analysis shows that statistical moments (and PDFs) of C and C p obtained within a GSG Y field identified by relatively large α values, i.e. α = 1.8, and a Gaussian Y field are virtually indistinguishable. This result is consistent with the recent findings of [394]. These authors explore analytically lead-order effects that non-Gaussian heterogeneity described by the GSG model have on the stochastic description of flow and transport under uniform in the mean flow in two-dimensional unbounded randomly heterogeneous media. Their results indicate that differences between lead-order flow and transport moments associated with GSG and GaussianY ;7 fields tend to diminish asα approaches 2, becoming virtually unnoticeable forα≥ 1.8. Similar to these authors, our results indicate the existence of a threshold value forα above which the effects associated with the non-Gaussian nature of the heterogeneous conductivity structure are virtually undetectableintheconcentrationBTCsrecordedatthewell. Avalueofα = 1.8canbeconsidered as a threshold above which the impact of the Y distribution (i.e. Gaussian vs non-Gaussian) is shadowed when compared to the influence of the engineering control (i.e., groundwater pumping rate) selected. Note that while these results appear to indicate that commonly employed Gaussian models could reproduce key transport features even in the presence of non-Gaussian Y fields, they also suggest that it would be difficult to differentiate between Gaussian and non-Gaussian Y fields on the basis of such moments whenα is close but not equal to 2. Such a distinction can be validly drawn only by analyzing Y data and their increments jointly, as suggested by Riva et al. [398]. The outcomes of our work associated with the feedback between engineering factors (i.e., transient versus uniform pumping rates) and efforts aimed at the characterization of aquifer het- erogeneous structure (through Gaussian or Sub-Gaussian models) on the behavior of contaminant BTCs are of potential interest to direct technical and economical efforts towards an optimal man- agement of groundwater resources. For example, costs linked to an increase of the well pumping rate could be justified by the production of water characterized by low contaminant concentra- tions, which in turn leads to decreased water treatment costs. Our results also suggest that pumping operations can control the temporal patterns of risk and might overshadow the impact of the type of aquifer heterogeneity (as embedded in the functional format of the probability density function characterizing hydraulic conductivity) on BTCs at pumping wells. In this frame- work, there could be circumstances in which enhanced efforts should be allocated towards an improved optimal planning of the pumping regime as opposed to a detailed characterization of some features of the heterogeneous properties of an aquifer. Such allocation of resources is key to reduce the uncertainty in risk metrics and is well aligned with goal-oriented site characterization frameworks (see [5:] and references therein). As a future research outlook, it would be of interest ;8 to extend our analysis to investigate transport in realistic systems of increased level of complexity that incorporate stochastic fluctuations for water demand. Additionally, extendingthefindingsofourworktothreedimensionalaquifers’configurationsin thepresenceofpumpingwellsoperatingwithtransientratesisfocusoffutureresearch. Webelieve thatincreasingthedimensionalityofthesystem,i.e. froma 2Dtoa 3Dconfiguration,wouldenable to capture more realistic flow paths that would potentially enhance solute mixing. As shown in Dentz et al. [6:], the uncertainty of the overall solute dispersive behavior and its self-averaging properties are affected by the dimensionality of the flow field. The dilution enhancement induced by the additional degree of freedom within a three-dimensional setting would yield a decrease of the solute concentration variability across Monte Carlo realizations with an ensuing decrease of the associated variance. Varying the dimensionality of the flow field (i.e., considering a three- dimensional system) might also affect the scaling behavior of the contaminant BTCs observed at the operating well (e.g., [37:,37;]) due to increased connectivity of the permeability field [73]. In the upcoming chapter 7 we present a real-case study where we shift our attention from pumping operations to the impact of variable climatic conditions on subsurface contaminant fate and transport at contaminated sites under sustainable remediation. ;9 Chapter 7 Climate change impact on residual contaminants under sustainable remediation The content of this chapter is published in Journal of Contaminant Hydrology [343]. 7.3 Introduction Subsurface contamination is recognized as a critical issue in many communities. There are more than a thousand Superfund sites in the U.S. categorized under the U.S. Environmental Protection Agency (EPA), which present large plumes of organic solvents, heavy metals and radionuclides [433]. In addition, the EPA estimated the presence of more than 450, 000 brownfields in the U.S., contaminated by hazardous substances, pollutants or contaminants [432]. The practice of soil and groundwater remediation has been changing in recent decades, tran- sitioning from intense soil removal and active treatment solutions towards passive remediation techniques such as enhanced biodegradation or monitored natural attenuation [78]. The latter less-intensive remediation approaches, often identified as sustainable remediation techniques, have been recognized as more advantageous since they can reduce negative side effects that often ac- company intensive remediation (e.g., ecological disturbance, construction noise/traffic, intensive energy use and emissions of greenhouse gases). Additionally, sustainable remediation is coupled ;: with more innovative and attractive end-use scenarios with restricted subsurface use and longer institutional control. In most cases, a portion of contaminants are sequestered in the subsurface for a long period of time, while natural or enhanced biogeochemical processes occur to reduce contaminant concentrations. Within this context, it is critical to assess the long-term stability of residual contaminants subject to sustainable remediation practices and to ensure that the latter will not pose significant environmental and human health risk. Recently, [36:] raised the concern that climate change may pose a major risk in environmental remediation; especially with regard to the fate of residual contaminants under sustainable reme- diation. A hydrological shift has been identified as one of the key drivers influencing such risk and uncertainty. In changing climate, precipitation regimes (including amount and timing) are expected to change significantly, and extreme events, such as heavy rains and prolonged droughts, could become more frequent. Climate models also predict increasing temperatures, which would affect water budgets and reduce infiltration due to increased evapotranspiration. These climatic changes are occurring while groundwater is becoming increasingly important for drinking and irrigation purposes [82]. Despite the critical need to evaluate risks associated with climate change, there is only a limited number of studies that address the effects of climate change on contaminant transport and environmental remediation. While the impact of climate change has been investigated extensively fromtheperspectiveofwaterresources[:4,:8,35:,383],alimitednumberofstudieshaveaddressed water quality issues [439]. Moreover, most literature focuses on surface water [9;,3;2,434,435, 447], due to its visibility and accessibility [:8], while the studies on groundwater are mostly focused on agricultural effluents at the regional scale [3:,57,72,9;,342,374,378,3;;,446,447]. Current literature lacks site-scale hydrological studies that can guide sustainable remediation under changing climate conditions within risk and performance assessments as well as within regulatory frameworks. ;; This study aims at evaluating the effects of climate-driven hydrological shifts on residual contaminants in vadose zones and groundwater under sustainable remediation. We simulate groundwater flow and contaminant transport through unsaturated and saturated domains. We assume that the effect of changing precipitation and temperature can be represented by pertur- bations/shifts of natural recharge through the aquifer system. The impacts are evaluated on the basis of different decision metrics relevant to public health risk, regulatory compliance and site closure such as contaminant concentrations at monitoring wells and exports from site boundaries. We demonstrate our approach at the Department of Energy (DOE)’s Savannah River Site (SRS) F-Area Seepage Basins, South Carolina (SC), USA, where soil and groundwater were contam- inated by various metals and radioactive contaminants during the Cold War Era. For brevity, the F-Area Seepage Basins are referred to hereafter as just F-Area. Extensive subsurface char- acterization and dataset at the SRS F-Area, including hundreds of wells, geophysics data and various hydrological tests, enabled the development of a subsurface model that can be considered as a testbed for flow and transport studies [34,56,92,3:;,442]. The SRS hence provides a unique opportunity to investigate the potential consequences of climate change on residual contamination in realistic settings that display multiple representative features of other polluted sites. 7.4 Background 7.4.3 Conceptual model We consider a general conceptual hydrological model of a contaminated site characterized by residual contamination in the vadose zone and groundwater (Figure 7.3). This model extends to the SRS F-Area [34], the description of which is provided in Section 7.4.4. Initially, contaminants are discharged through a seepage basin located on the top of the model domain. The contaminant plume migrates vertically through the vadose zone, and laterally downgradient in the aquifer. At some sites, the plume reaches surface water bodies (e.g., a creek or a river) located close to the site 322 Figure 7.3: Schematic illustration of the hydrological conceptual model under investigation, rep- resenting a vertical two-dimensional cross-section driven along the middle line of the contaminant source zone. boundary, through water seepage. In order to reduce contaminant migration through the vadose zone, the contaminant source zone (i.e., seepage basin) is often capped with low-permeability material, such as clay or silt. However, residual contaminants located in the vadose zone could becomeapersistentcontaminantsourcetothegroundwaterplume[453], asshownintheschematic illustration of Figure 7.3. Contaminant concentrations are typically measured at groundwater monitoring wells to ensure the plume stability and to meet the regulatory compliance. The number of wells and frequency of sampling (e.g., quarterly or annually) are determined in agreement with a regulatory agency (e.g., the EPA). Concentrations are often compared to the Maximum Contaminant Level (MCL) recommendedbytheEPA.PredictingthetimeofexceedanceofMCL,forexample, isanimportant EnvironmentalPerformance Metric (EPM) to plan the site closure or site transfer [68]. In parallel, risk (and/or performance) assessments consider a variety of metrics and pathways, including contaminant mass flux/export at control planes (CPs), e.g., site boundaries, creeks or rivers [8,5:,66,67,356]. On-site concentrations are often used to evaluate human health risk doses through direct ingestion/drinking pathways, while the export at CPs are employed to quantify risk of the downgradient population. 323 The changes in precipitation and temperature associated with climate change are expected to affect contaminant plumes in groundwater systems and residual contaminants in the vadose zone in a complex manner. This impact can be evaluated by considering a perturbation/shift of the natural recharge through an aquifer in a long or short time frame [36;]. For example, an increase in precipitation results in higher aquifer recharge, while a decrease in precipitation, or higher temperature, hence increased evapotranspiration, leads to lower aquifer recharge. Higher recharge then (3) increases vadose-zone flow, which mobilizes sequestered contaminants, (4) raises the groundwater table and increases hydraulic gradients, resulting in enhanced plume mobility in groundwater, and (5) enhances mixing of the plume with clean water, which leads to higher dilution. On the other hand, a decrease of recharge has opposite effects, i.e, decreases the plume mobilitywhilereducingmixinganddilution. Theimpactofclimatechange-drivenalteredrecharge rates on different decision and performance metrics relevant to environmental remediation could potentially create trade-offs that should be quantitatively evaluated. 7.4.4 F-Area site description The SRS is located in South-central South Carolina, USA, approximately 100 mi (i.e., 161 km) away from the Atlantic Ocean and occupies an area of about 800 km 2 . The site was used to produce nuclear materials such as plutonium and tritium ( 3 H), for nuclear weapons during the Cold War Era. The F-Area, located in the north-central part of the SRS, included three unlined discharge basins: F-1, F-2 and F-3 [34]. The basins received approximately 7.1 billion liters of acidic, low-level radioactive waste solutions from processing irradiated uranium between 1955 and 1988 [93,32:]. The waste solution presented various radionuclides such as uranium isotopes, 90 Sr, 129 I and 99 Tc, among which 3 H, is the largest dose contributor. After the waste discharge operationwasterminatedin 1988, theF-Areabasinswereclosedandcappedwithlow-permeability material [34]. Currently, an acidic contaminant plume extends from the basins approximately 600 m downgradient to the groundwater seepage near the Fourmile Branch [34]. Enhanced natural 324 Hydrostratigraphic Unit Porosity φ [−] Permeability k [m 2 ] α [kg −1 m s 2 ] n [−] m [−] S rl [−] p [−] UUTRA 0.39 5×10 −12 4×10 −4 2 0.5 0.18 1 TCCZ 0.39 1.98×10 −14 5.1×10 −5 2 0.5 0.39 1 LUTRA 0.39 5×10 −12 5.1×10 −5 2 0.5 0.41 1 Table 7.3: Parameters adopted in the numerical model presented in this study. The table reports values for porosity (φ [−]), permeability (k [m 2 ]), water retention-curve parameters (α [kg −1 m s 2 ], n [−], m [−]), residual liquid pore saturation (S rl [−]), Mualem (1976) [364] parameter (p [−]) (for details, see [34]) attenuation is currently under way using a funnel-and-gate system which consists of groundwater flow barriers to decrease the groundwater gradient, and base injection to neutralize pH and to immobilize uranium [426]. Hydrogeology at this site has been characterized extensively in many studies [34,93,3:;,442]. There are three hydrostratigraphic units within the Upper Three Runs Aquifer: an Upper Aquifer zone (UUTRA), a Tan Clay Confining Zone (TCCZ), and a Lower Aquifer zone (LUTRA). The UUTRA and LUTRA are mainly composed by clean sand, while the TCCZ is a low-permeable mixed sand-and-clay layer. The piezometric head measurements indicate that the UUTRA and LUTRA units are hydrologically connected. The bottom of the LUTRA consists of a competent clay layer confining unit, the continuous Gordon Confining (GC) unit, which separates the deeper aquifer (Gordon Aquifer) from the upper two units. The historical monitoring data collected at the SRS have shown that the contaminant plume migrates within the UUTRA and LUTRA [3:;]. 7.5 Methodology 7.5.3 Flow and transport simulations We employ the two-dimensional (2D) flow and transport model adopted in [34], i.e., a 2D domain approximately 2600 m long and 100 m deep along the groundwater flow line, passing through the middle of the F-3 basin of the SRS. This model has been calibrated and verified using site data [34]. Figure 7.4 illustrates the 2D cross-section model domain. 325 50 m 2600 m Fourmile Branch Creek Seepage Face FSB110D Seepage F-3 Basin Recharge FSB95D No Flow UUTRA LUTRA TCCZ No Flow No Flow GC Gordon c.z. Low er UTRA a.z Tan Clay c.z. Upper UTRA a.z GC TCCZ Hy drostratigraphic Units UUTRA TCCZ LUTRA Figure 7.4: Two-dimensional model domain adopted in our study. The model includes the vadose zone and three hydrostratigraphic units (i.e., UUTRA, LUTRA and TCCZ) defined in Section 7.4.4. We assume homogeneous average hydrogeological proper- ties within each unit (see Table 7.3), whose values are compiled from available site investigation reports. Table 7.3 specifies porosity, permeability and capillary pressure/saturation data for the vadose zone [34,93,384]. Please refer to Appendix A of [34] for the parameters description. The system is considered to be advective dominated, therefore mechanical dispersion and molecular diffusion transport processes are neglected. We simulate simple 3 H decay with a half-life of 12.3 years. No-flow boundary conditions are assigned along the two vertical sides of the 2D-cross sec- tion (see Figure 7.4) according to the groundwater divides [34,92]. An impervious flow boundary condition (i.e. no-flow) is set at the bottom of the computational domain, since the confining unit at this location is highly clayey and continuous [34]. Groundwater flow and contaminant transport are simulated by means of the numerical code Amanzi which describes coupled vadose zone and groundwater flow as well as reactive transport [34,98,43:,43;]. Amanzi uses the mimetic finite 326 difference method for the Richards equation and dispersion operators, which preserves fundamen- tal mathematical and physical properties in discrete schemes [4;]. To discretize the advection operator, it uses monotone first-order and second-order MUSCL schemes with new limiters that improve accuracy on unstructured meshes [349]. Amanzi has been benchmarked against other flow and transport models as well as analytical solutions for a wide range of hydrological problems. We perform numerical simulations within the time frame 1955− 2100, i.e., from the beginning of the discharge operation at the SRS. Our model is discretized into 164160 cells and we adopt a longitudinal mesh spatial resolution of 1.25 m and a variable vertical spatial resolution, with a maximum value of 2 m. A steady-state simulation is carried out to establish the groundwater table with a given recharge value at the top of the modeled domain of 4.74× 10 −6 mm/s (i.e., 150 mm/yr) before 1955. This value represents the estimated recharge based on the rainfall records and runoff estimations [93], which is kept constant over the entire domain for all the simulation timeframe, exceptforthedischarge basin. After 1955, we perform transient simulations employing a constant 3 H discharge value at the seepage basin between 1955 and 1988. The average 3 H concentrationandmassdischargeratearerespectively 2.17×10 −9 mol/kgwandaround 1.11×10 −4 mm/s, as in [34]. These values coincide with the average 3 H concentration and discharge rate of historical data between 1955 and 1988 (see Table : of [34]). After the basin is capped in 3;::, the recharge through the basin is assumed to be three-order of magnitude lower than the value of the surrounding regions for all the simulations performed in our study. These recharge values are identified as baseline recharge conditions in our study. We model different recharge scenarios which present increased or decreased recharge values with respect to the baseline recharge conditions described above. We analyze the impact of recharge shifts on different decision (or performance) metrics: the temporal evolution of 3 H con- centration at the same two observation wells of [34], shown in Figure 7.4: (3) the source-zone UUTRA well (FSB95D) and (4) the downgradient UUTRA well (FSB110D). Both wells span the upper aquifer. Please note that we computed the contaminant concentration by taking the 327 average value of the concentration reported at a set of observation points equally spaced within a given well in the upper aquifer. In addition, we evaluate the 3 H export at the CP to the Fourmile Branch Creek (see Figure 7.4), which is the main risk pathway at this site. The CP is defined at the seepage face (indicated in Figure 7.4) where the groundwater flow reaches the surface. The concentrations are compared with the EPA’s MCL for 3 H. Although the MCL criteria is not used for compliance purposes at the SRS, it has been used to evaluate the timeframe for the site closure and transfer [68]. 7.5.4 Modeling scenarios Previous studies on the SRS have reported an overall trend towards greater rainfall in the re- gion [84]. According to the National Climate Assessment, South Carolina is expected to see precipitation increases of 10%− 20% by 2100 (see Fig. 9.7 in [76]), as well as more extreme pre- cipitation event by up to two-folds (see Fig 9.8 in [76]). Even if the amount of precipitation is not necessarily equal to the recharge, we assume that we can investigate the impact of climate change by simulating a range of different recharge values. The range would also account for the uncer- tainty associated with the climate projections. We mainly focus on increased natural recharge conditions, although, for completeness, we also investigate the impact of reduced recharge. Figure 7.5 provides a schematic illustration of the recharge scenarios simulated in our analysis. We con- sider a baseline recharge, denoted asR B , and indicated by the black line in Figure 7.5. To study the effect of climate change-induced increased/decreased recharge on contaminant transport, we definetobeaperturbationfromthebaselinerechargeandR P toidentifytheperturbedrecharge. We develop four perturbed recharge scenarios with respect to the baseline recharge conditions, corresponding to = 0, described in the Section 7.5.3. The recharge perturbation starts at a certain time, set to year 2020, and is illustrated with a specific color for each scenario in Figure 7.5. In the following, we provide a detailed description of the scenarios and of the values adopted for and we employ t ∗ to indicate year 2020. 328 1955 1988 2020 2100 R B (1 - ) Stage I Radioactive W aste Discharge Stage II Capped Basin: Residual Contaminant Baseline Recharge (R B ) R B (1 + ) R B (1 + ) Time (+ Cap Failure) Recharge Figure 7.5: Schematic illustration of the perturbed recharge scenarios analyzed in the study: (1) theconstantpositiverechargeshiftstartingin 2020(greenline), (2)theconstantnegativerecharge shift starting in 2020 (blue line), (3) the one-year extreme recharge in 2020 (magenta line) and (4) the hypothetical cap failure and constant positive recharge shift starting in 2020 (green line). The baseline recharge scenario is denoted by the horizontal black line. 329 Constant positive recharge shift The recharge conditions of the first scenario (portrayed in green, from year 2020, in Figure 7.5) are illustrated through the following equation: R P (t) =    R B , t<t ∗ (7.3) R B (1 +), t≥t ∗ (7.4) through this scenario we simulate a constant positive shift (i.e., an increase) of recharge from 2020 until the end of the simulation, set to year 2100. The perturbation assumes the following range of values: = [0.1, 0.2, 0.3, 0.4, 0.5], we therefore simulate a recharge increase of +10% to +50% compared to baseline recharge conditions. Constant negative recharge shift By the same token, the second scenario (blue line from year 2020 in Figure 7.5) is illustrated through the following equation: R P (t) =    R B , t<t ∗ (7.5) R B (1−), t≥t ∗ (7.6) this scenario is characterized by a constant negative shift of recharge (i.e., a decrease) of -10% to -50% ( = [0.1, 0.2, 0.3, 0.4, 0.5] from equation (7.6)) within the timeframe 2020− 2100. One-year extreme recharge The third scenario (magenta line from 2020 in Figure 7.5) is described by: R P (t) =    R B , if t<t ∗ (7.7) R B (1 +), if t ∗ ≤t≤t ∗ + 1 (7.8) thisscenarioconsidersasignificantincreaseofrechargewithinyear 2020,mimickinganextreme precipitation event of a factor of 1.5, 2, 5 and 10 folds compared to baseline conditions. The range ofperturbationvaluesofthisthirdscenario, accordingtoequation(7.8), arethen: = [0.5, 1, 4, 9]. 32: Cap failure and constant positive recharge shift The cap failure is also evaluated in con- junction with the natural recharge shift. Although low-permeability material is used for capping the F-Area basins, increased vegetation or other mechanisms could threaten the integrity of the source-zone capping structure [449–44;]. In the fourth scenario, we assume that the cap will fail at the same time of the recharge increase (i.e., in year 2020). Note that the assumption of cap failure in 2020 is solely for the purpose of the modeling exercise and not a prediction of actual cap failure. The perturbed recharge conditions of this scenario are then the same of the first scenario, i.e.: R P (t) =    R B , if t<t ∗ (7.9) R B (1 +), if t≥t ∗ (7.:) The cap failure is represented by increased source-zone recharge to the level of the surrounding region, therefore we hypothesize a complete failure of the containment structure. Although such a sudden failure is unlikely to happen, we assume instant failure to evaluate the most extreme case. The fourth scenario simulates the failure of the capping structure in arbitrarily assumed year 2020 under baseline recharge conditions and under the increased recharge conditions of the first scenario, characterized by = [0.1, 0.2, 0.3, 0.4, 0.5], from 2020 to 2100. In reality, such precipitation/recharge changes are expected to occur gradually or randomly rather than through the step change assumed in this analysis. However, such simple represen- tations would facilitate our fundamental understanding of various complex impacts of climatic changes on contaminant concentrations and exports. For the conditions explored in this study (Figure 7.5) we expect that the changes in the recharge rates (as described in the scenarios above) will impact the concentration breakthrough curves (BTCs) by local dilution as well as by affecting the rate at which the contaminant mass is released from the vadose zone to the groundwater system. Figure 7.6 schematically pictures 32; Concentra on Time Phase I Phase II Recharge Time Baseline recharge Increased recharge Figure 7.6: Schematic illustration of the impact of increased recharge on the concentration break- through curve (BTC) at an observation well located downgradient from the source zone. The recharge rate increase causes a slight dilution (Phase I), followed by a “rebound” effect in the BTC (Phase II) due to contaminants’ mobilization. The continuous blu curve indicates the BTC under baseline recharge whereas the dashed red curve corresponds to the concentration BTC under a change in the recharge conditions (see inset figure). the contaminant BTC under baseline conditions and perturbed recharge conditions. We expect that, increased recharge (relative to the baseline values) causes a slight dilution (Phase I in Figure 7.5) followed by a “rebound” effect in the concentration BTC (Phase II in Figure 7.5) due to the mobilization of the solute mass located in proximity of the source (Figure 7.6). 7.6 Results and Discussion 7.6.3 Spatiotemporal Dynamics of the Contaminant Plume Prior to investigating the effects of recharge on transport observables at specific wells and at the CP, we provide snapshots of the simulated contaminant plume in the coupled vadose zone- groundwater system (Figure 7.7) for the baseline case. Figure 7.7 shows the plume migration at different years, namely 1961 (Figure 7.7b), 1988 (Figure 7.7c), 2000 (Figure 7.7d) and 2033 (Figure 7.7e). 332 Figure 7.7: Illustration of the simulated tritium plume in the SRS. The initial condition is dis- played in panel (a). The plume snapshots are shown for years (b) 1961, (b) 1988, (c) 2000 and (d) 2033. 333 As observed in [34], the contaminant plume first moves downward in the vadose zone, and then spreads laterally below the groundwater table (Figure 7.7b). During the operation, most of the plume migrates within the upper aquifer, although a part of the plume penetrates the TCCZ and reachestheloweraquifer(Figure7.7c). Aftertheoperationends, thecleanwaterfrontarrivesfrom upstream and pushes the plume downgradient (Figure 7.7d). The contaminant plume displays a stratified distribution since the residual contaminants in the vadose zone constitute a persistent contamination source and the low-permeability TCCZ becomes a secondary source. As displayed in Figure 7.7e, in year 2033 the vadose zone and TCCZ continue to be the residual contaminant sources. 7.6.4 Contaminant concentration and export BTCs We first compare the temporal evolution of the 3 H concentrations [mol kgw −1 ] at the observation wells and of the 3 H export [mol y −1 ] at the CP among different recharge scenarios (Figures 7.8- 7.;) for 0−100 years forward from the assumed recharge perturbations and cap failure (i.e., within the timeframe 2020− 2100). Under baseline recharge conditions (indiated by = 0 in Figures 7.8-7.;), 3 H concentrations and export generally decrease within the timeframe 2020− 2100, given that the peak concentration and export occur closer to the contaminant discharge operation timeframe (1955− 1988), i.e., before year 2020. When the positive constant shift of recharge occurs (Figure 7.8), the concentration initially decreases for approximately 2− 10 years at the source-zone well (Figure 7.8a) and for around 5− 20 years at the downgradient well (Figure 7.8b). This slight decrease of concentration is an outcome of dilution, attributed to the presence of more water in the system which enhances the mixing of the plume with clean water. We then observe a “rebound” in the concentrations due to the additional recharge, during which the concentration increases with respect to the concentration BTC obtained for the baseline case ( = 0). This is attributed to the fact that the residual contaminants in the vadose zone are mobilized under higher recharge. As a consequence, 334 2020 2040 2060 2080 2100 Time [y] 0 1 2 3 4 5 6 C [mol kgw -1 ] 10 -12 Positive recharge shift (a) Source-zone well =0 =0.1 =0.2 =0.3 =0.4 =0.5 2020 2040 2060 2080 2100 Time [y] 0.5 1 1.5 2 2.5 3 C [mol kgw -1 ] 10 -12 (b) Downgradient well =0 =0.1 =0.2 =0.3 =0.4 =0.5 2020 2040 2060 2080 2100 Time [y] 5 10 15 Export [mol y -1 ] 10 -6 (c) Control plane =0 =0.1 =0.2 =0.3 =0.4 =0.5 Figure7.8: Temporal evolution of 3 H concentration at the: (a) source-zone well, (b) downgradient well and (c) 3 H export at the CP for the baseline scenario, indicated by = 0, and the constant positive recharge shift scenario, characterized by = 0.1, 0.2, 0.3, 0.4, 0.5. The results of different recharge perturbations are portrayed in different colors. The thin horizontal dashed black line represents the MCL for 3 H. 335 2020 2040 2060 2080 2100 Time [y] 1 2 3 4 5 6 C [mol kgw -1 ] 10 -12 Negative recharge shift (a) Source-zone well =0 =0.1 =0.2 =0.3 =0.4 =0.5 2020 2040 2060 2080 2100 Time [y] 0.5 1 1.5 2 2.5 3 C [mol kgw -1 ] 10 -12 (b) Downgradient well =0 =0.1 =0.2 =0.3 =0.4 =0.5 2020 2040 2060 2080 2100 Time [y] 5 10 15 Export [mol y -1 ] 10 -6 (c) Control plane =0 =0.1 =0.2 =0.3 =0.4 =0.5 Figure7.9: Temporal evolution of 3 H concentration at the: (a) source-zone well, (b) downgradient well and (c) 3 H export at the CP for the baseline scenario, indicated by = 0, and the constant negative recharge shift scenario, characterized by = 0.1, 0.2, 0.3, 0.4, 0.5. The results of different recharge perturbations are portrayed in different colors. The thin horizontal dashed black line represents the MCL for 3 H. 336 2020 2040 2060 2080 2100 Time [y] 5 10 15 C [mol kgw -1 ] 10 -12 One-year extreme recharge (a) Source-zone well =0 =0.5 =1 =4 =9 2020 2040 2060 2080 2100 Time [y] 1 2 3 4 C [mol kgw -1 ] 10 -12 (b) Downgradient well =0 =0.5 =1 =4 =9 2020 2040 2060 2080 2100 Time [y] 0.5 1 1.5 2 2.5 3 Export [mol y -1 ] 10 -5 (c) Control plane =0 =0.5 =1 =4 =9 Figure7.:: Temporal evolution of 3 H concentration at the: (a) source-zone well, (b) downgradient well and (c) 3 H export at the CP for the baseline scenario, indicated by = 0, and the one-year extreme recharge scenario, characterized by = 0.5, 1, 4, 9. The results of different recharge perturbations are portrayed in different colors. The thin horizontal dashed black line represents the MCL for 3 H. 337 2020 2040 2060 2080 2100 Time [y] 0 2 4 6 8 C [mol kgw -1 ] 10 -11 Cap failure (a) Source-zone well =0 =0 (c) =0.1 (c) =0.2 (c) =0.3 (c) =0.4 (c) =0.5 (c) 2020 2040 2060 2080 2100 Time [y] 10 -15 10 -13 10 -11 10 -9 C [mol kgw -1 ] 2020 2040 2060 2080 2100 Time [y] 0.5 1 1.5 2 C [mol kgw -1 ] 10 -11 (b) Downgradient well =0 =0 (c) =0.1 (c) =0.2 (c) =0.3 (c) =0.4 (c) =0.5 (c) 2020 2040 2060 2080 2100 Time [y] 10 -15 10 -13 10 -11 10 -9 C [mol kgw -1 ] 2020 2040 2060 2080 2100 Time [y] 5 10 15 Export [mol y -1 ] 10 -6 (c) Control plane =0 =0 (c) =0.1 (c) =0.2 (c) =0.3 (c) =0.4 (c) =0.5 (c) 2020 2040 2060 2080 2100 Time [y] 10 -10 10 -8 10 -6 10 -4 Export [mol y -1 ] Figure7.;: Temporal evolution of 3 H concentration at the: (a) source-zone well, (b) downgradient well and (c) 3 H export at the CP for the baseline scenario, indicated by = 0, the cap failure under the baseline scenario, indicated by = 0 (c) , and the cap failure under the constant positive recharge shift scenario, indicated by = 0.1 (c) , 0.2 (c) , 0.3 (c) , 0.4 (c) , 0.5 (c) , with (c) indicating the capping failure. The results of different recharge perturbations are portrayed in different colors. The thin horizontal dashed black line represents the MCL for 3 H. The inset shows the log-scale plot. 338 more solute mass reaches the wells and the CP. The concentration increase during the “rebound” phase is more pronounced at the source-zone well (Figure 7.8a) than at the downgradient well (Figure 7.8b) given that the latter is located further away from the source zone and consequently less influenced by the mobilization of residual 3 H. These results are in agreement with previous theoretical analysis which showed that near-source locations are more sensitive to changes in the solute mass release at the source zone [58]. By comparing the effects of different magnitudes of recharge shifts on well concentrations, we notice that a recharge increase characterized by = 0.1 does not influence the concentration BTC significantly. On the other hand, Figure 7.8a and 7.8b show that, for larger recharge, the timing of the “rebound” effect happens before and the corresponding peak concentration is higher (i.e., see concentration BTCs within the range = 0.2− 0.5). We also point out that the concentration decrease after the peak is generally more rapid as increases (see BTCs produced by = 0.4− 0.5 in Figure 7.8a-b). Indeed, under higher recharge, a bigger part of the residual 3 H is mobilized and flushed out of the aquifer system earlier in time. The export, on the other hand, slightly tends to increase but does not change significantly even under the maximum increase of recharge, indicated by = 0.5 (Figure 7.8c). The minor effects on the export are attributed to the fact that the export is an integrated measure that incorporates the upwelling groundwater from the lower aquifer which is less affected by the recharge changes. Decreased recharge (Figure 7.9) produces minor effects on well concentrations compared to increased recharge (compare Figure 7.8a-b to Figure 7.9a-b). Lower recharge causes a small increase in the concentration for around 5− 20 years at the source-zone well (Figure 7.9a) and for approximately 10− 20 years at the downgradient well (Figure 7.9b). Figure 7.9a depicts a higher peak concentration as increases. This effect takes place because a lower amount of water in the system decreases the dilution potential of the aquifer hence leading to an increase the concentration. The differences between the concentration BTCs obtained for different values at both observation wells in Figure 7.9a-b are less visible than in Figure 7.8a-b. We observe 339 that increased recharge has a larger impact on the source-zone well (Figure 7.9a) than on the downgradient well (Figure 7.9b), whereas decreased recharge causes more uniform effects between the two observation wells. This is because the mobilization of 3 H caused by increased recharge has more impact closer to the source (Figure 7.9a), whereas the decrease of dilution, emerging from decreased recharge and resulting in higher concentration (Figure 7.9a-b) has a more uniform effect on the aquifer system. Finally, minor effects are displayed on the export (Figure 7.9c), as noted also for increased recharge conditions (see Figure 7.8c). However we notice that the 3 H export tends to slightly decrease under negative recharge shifts given that less residual 3 H is mobilized from the vadose zone. We next observe the outcomes of the third scenario (i.e., the one-year extreme recharge event) in Figure 7.:. One-year of extreme recharge significantly increases the well concentrations during the“rebound”phaseoveranextendedperiodoftime, i.e., foraround 10−20years(Figure7.:a-b). Wethenobservethatashortextremeeventcouldaffectthewellconcentrationsforseveraldecades. At the source-zone well (Figure 7.:a), prior to the “rebound” phase, the concentration decreases slightly after the perturbation event (in 2020) due to dilution, particularly under the most extreme scenario, indicated by = 9. On the other hand, at the downgradient well (Figure 7.:b), the concentration decrease due to dilution is more pronounced given that the plume has traveled a longer distance. Five-to-ten years after the initial dilution, the concentration increases, during the “rebound” phase. The extreme one-year precipitation indicated by = 0.5 does not significantly affect the concentrations, especially at the downgradient well, where even the recharge scenario characterized by = 1 does not produce significant changes when compared to baseline conditions ( = 0). We also observe that higher recharge shifts produce higher peak concentrations in the “rebound” phase due to the mobilization of the solute mass at the source zone (Figure 7.:a-b) as well as higher initial dilution. The effect of recharge changes on the export (Figure 7.:c) is minor compared to the influence on the concentrations in the wells (Figure 7.:a-b), although we observe a peak of the export around 2020 when = 9. 33: Figure 7.; reports the results obtained for the scenario characterized by the failure of the capping structure under baseline conditions and constant positive recharge shifts starting in 2020. Figure 7.; shows that the concentration increase at the wells is significantly larger than the increase observed for undamaged source-zone capping conditions. The concentration “rebound” happens after 5− 8 years at the source-zone well (Figure 7.;a) and after 10− 15 years at the downgradient wells (Figure 7.;b). The export at the CP also shows a visible increase after around 2035 (Figure 7.;c). This is due to the fact that, under higher recharge and no capping structure, the bulk of residual contaminants in the vadose zone migrate to the groundwater system. We notice that the increase of concentration/export happens earlier in time at the well closer to the contaminant source zone, later at the downgradient well and even further when considering the export at the CP, located at the right boundary of the domain. Moreover, as the recharge value increases, from = 0.1 to = 0.5, the peak concentration/export, caused by the additional recharge, becomes higher and occurs earlier in time (see Figure 7.;). The inset log-scale plot of Figure 7.;a shows that the extreme concentration increase resulting from the cap failure is followed by a rapid decrease of the concentration, as the 3 H plume is flushed out of the system earlier when the capping structure fails, moving the contamination problem downgradient. This effect is amplified as increases. The same observations apply to the downgradient well and to the export (see inset log-scale plot of Figure 7.;b-c). We finally notice that the difference between the BTCs given by different values is more pronounced when we assume that the cap fails (Figure 7.;). The results of an additional recharge scenario, characterized by a range of positive recharge shifts within a timeframe shorter than a year, is presented in the Supplementary Material. These outcomes confirm that our findings also apply to a smaller time scale of hydrological shift. 33; 7.6.5 ImpactofrechargeperturbationsonkeyEnvironmentalPerformance Metrics We quantify the impact of climate change-driven recharge shifts on key EPMs (e.g, peak con- taminant concentration, early and late arrival times, time of exceedance of MCL). Given that the recharge shifts considered in our study start from year 2020 (after the discharge operation at the SRS F-Area terminated, i.e. during a contaminant concentration/export descending/tailing phase), we do not consider the peak concentration but we analyze the maximum change on the concentration and export BTCs induced by the recharge perturbations. Moreover we investigate the effect of recharge shifts on the time of exceedance of the MCL for 3 H. 7.6.5.3 Maximum difference between baseline and perturbed BTCs We identify the maximum normalized difference between the contaminant concentration or export BTC obtained under baseline and perturbed recharge conditions as: γ c = 1 C b (t dmax ) max t |C b (t)−C p (t)|×100, (7.;) where C b (t) and C p (t) in (7.;) are the contaminant concentration [mol kgw −1 ] or contaminant export [mol y −1 ] of the baseline (subscript “b”) and the perturbed (subscript “p”) recharge sce- narios respectively, whereas C b (t dmax ) is the baseline concentration or export taken at the time where the difference between baseline and perturbed BTCs is maximum. The metric γ c (7.;) is plotted in Figure 7.32 as a function of the recharge perturbation, whose values are described in Section 7.5.4. Note that, from (7.;), γ c = 0 for the baseline scenario, indicated by = 0. We recall that only a limited range of recharge perturbations were simulated in our analysis, therefore the comments below are based on a linear interpolation between γ c and the analyzed values. When the recharge is positively shifted by a constant value (Figure 7.32a), the relative change in the well concentrations, quantified by γ c (7.;), exhibit a non-linear or step-function response 342 to the perturbation . As shown in Figure 7.32a, γ c at the source-zone well drastically increases (from around 10% to almost 120%) when the recharge perturbation changes from 0.1 to 0.2, whileγ c atthedowngradientwellishighlyimpacted(itchangesfromapproximately 20%to 100%) when goes from 0.3 to 0.4. The changes are smaller during the other recharge shifts’ intervals. Our results suggest the presence of a threshold, or trigger level of recharge, above which the well concentrations are significantly affected. This trigger level is lower at the source-zone well than at the downgradient well; a recharge corresponding to = 0.1 could represent this threshold at wells closer to the source-zone, while a higher trigger level of recharge, approximately identified by = 0.3, is identified downgradient. The export, on the other hand, increases in a quasi-linear manner with the recharge perturbation (Figure 7.32a). As shown in the concentration time series (see Figure 7.8c), the effect of recharge shifts on the export is significantly smaller compared to that on the well concentrations. Indeed a recharge shift corresponding to = 0.5 results in a value ofγ c of approximately 20% when observing the CP export, while a recharge shift given by = 0.5 produces values of γ c of approximately 80%− 100% at the observation wells. Figure 7.32b shows theγ c - relationship for the negative recharge shift scenario. Decreasing therechargerategenerallyhasasmallerimpactonγ c withrespecttoincreasedrechargeconditions (compare the values of γ c in Figure 7.32a and 7.32b), according to the results in Figure 7.9. We identify non-linear responses and trigger levels of recharge perturbation at the downgradient well and at the CP. The trigger levels are very small at the downgradient well (e.g., less than = 0.1) and higher for the export (i.e., = 0.2). The maximum difference between the baseline case and the perturbed BTC (computed via equation (7.;)) is larger at the downgradient well than at the source-zone well (see Figure 7.32b). This is different than what was observed from the results of the positive shift scenario. When comparing Figures 7.32a and 7.32b we observe that, under positive shifts of recharge (Figure 7.32a), γ c evaluated at the source-zone is more sensitive to smaller values. In other terms, we could say that the source-zone well responds first to the perturbation . This is because the source-zone well is more impacted by the mobilization of 3 H 343 0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 100 120 Positive recharge shift (a) Source-zone well Downgradient well Control Plane 0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40 50 Negative recharge shift (b) Source-zone well Downgradient well Control Plane 0 0.5 1 4 9 0 100 200 300 400 500 600 700 One-year extreme recharge (c) Source-zone well Downgradient well Control Plane 0 0 (c) 0.1 (c) 0.2 (c) 0.3 (c) 0.4 (c) 0.5 (c) 0 1000 2000 3000 4000 5000 6000 7000 Cap failure (d) Source-zone well Downgradient well Control Plane Figure 7.32: Maximum difference between contaminant concentrations and export BTCs of the baselinescenarioandtheperturbedrechargescenarios(γ c )versusperturbation()fortheconstant positive recharge shift scenario (a), the constant negative recharge shift scenario (b), the one- year extreme recharge scenario (c), the cap failure and constant positive recharge shift scenario (d), where the superscript (c) indicates the capping failure. Results at the source-zone well are illustrated in red, at the downgradient well are pictured in blue and at the CP are indicated in green. 344 givenitsproximitytothesource. However, when the recharge decreases(Figure7.32b), changesin γ c are more prominent at the downgradient well for smaller perturbations, i.e., the downgradient well responds first than the source-zone well. Indeed, the changes in the concentration are more uniform between the source-zone well and the downgradient well (compare Figures 7.9a and 7.9b) under lower recharge because they happen due to the concentration (or less dilution) as opposed to the mobilization of contaminants. Moreover, these changes take place at a later time downgradient, when the baseline concentration is smaller, resulting in higher γ c downgradient. Under the extreme one-year recharge scenario (Figure7.32c), the trigger-level is larger (around = 1) when compared to the results obtained for the constant positive recharge shift scenario at both observation wells (compare Figures 7.32a and 7.32c). We observe that the γ c - relationship is quite similar for the two observation wells until = 4. However, when = 9, the recharge increase impact at the downgradient well becomes quite significant. In quantitative terms, this impact is approximately two times that at the source zone well (compare red and blue lines in Figure 7.32c). Finally, the impact of the extreme event on the contaminant export is quite limited and increases linearly with the perturbation (Figure 7.32c). However, a significant increase of the export is detected for the most extreme one-year precipitation event investigated in this study (i.e., when = 9, γ c equals 65%). The cap failure scenario (Figure 7.32d) also shows a non-linear response when evaluating the sensitivities of the wells’ concentrations and the export to the perturbation . We observe an extreme increase of γ c caused by the signification mobilization of the residual contaminants entrapped in the vadose zone when the capping structure fails, confirming the importance of the latter containment system. Successive increases of recharge do not significantly affect the value of γ c , i.e., γ c reaches a plateau under capping failure conditions. The almost constant value of γ c is generally higher at the source-zone well than at the donwgradient well and at the CP. Our analysis then suggests that the failure of the capping can be identified as a trigger situation after which major changes in both well concentrations and on the export are expected. 345 0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40 Positive recharge shift (a) Source-zone well Downgradient well 0 0.1 0.2 0.3 0.4 0.5 0 20 40 60 80 Negative recharge shift (b) Source-zone well Downgradient well 0 0.5 1 4 9 0 50 100 150 One-year extreme recharge (c) Source-zone well Downgradient well 0 0 (c) 0.1 (c) 0.2 (c) 0.3 (c) 0.4 (c) 0.5 (c) 0 100 200 300 400 Cap failure (d) Source-zone well Downgradient well Figure 7.33: Difference of time of exceedance of the MCL for 3 H between the baseline scenario and the perturbed recharge scenarios (γ t ) versus perturbation () for the constant positive recharge shift scenario (a), the constant negative recharge shift scenario (b), the one-year extreme recharge scenario(c), thecapfailureandconstantpositiverechargeshiftscenario(d), wherethesuperscript (c) indicates the capping failure. Results at the source-zone well are illustrated in red and results at the downgradient well are pictured in blue. 7.6.5.4 Difference of time above MCL between baseline and perturbed BTCs In this Section, we identify the metricγ t as the normalized difference of the time of exceedance of the MCL for 3 H between the baseline recharge scenario and the perturbed recharge scenarios as: γ t = 1 t(C b (t)>MCL) |t(C b (t)>MCL)−t(C p (t)>MCL)|×100. (7.32) with t(C b (t) > MCL) and t(C p (t) > MCL) respectively being the time [d] during which the contaminant concentration exceeds the EPA’s MCL for 3 H for the baseline (subscript “b”) and the perturbed (subscript “p”) recharge scenarios. As before, the metric γ t (7.32), expressed in percentage, is plotted as a function of for each recharge scenario in Figure 7.33 and, from (7.32), γ t = 0 for the baseline scenario (i.e., when = 0). Similarly to before, our observations are based on a linear interpolation between γ t and the range of values considered. 346 Under the positive shift of recharge scenario (Figure 7.33a), γ t shows a non-linear response to the recharge perturbation at the source-zone well (in red). In this case, γ t increases until its maximum value under a perturbation = 0.2 and = 0.3 but decreases afterwards. We then notice that the intermediate values of produce the highest influence on the timeframe of MCL exceedance, as compared to the highest values of. Indeed, positive recharge shifts initially dilute the plume in the tailing phase, and increase the 3 H mobilization, leading to an increase of the concentration during the “rebound” phase, but afterwards flush the plume out of the system, i.e., at a later stage the 3 H concentration values reach the 3 H MCL faster (see BTCs characterized by = 0.4− 0.5 in Figure 7.8a). The higher the recharge, the higher is the concentration “rebound” peak and the faster the plume is flushed out of the aquifer, therefore the maximum increase of the time above MCL is produced by intermediate value of recharge (for example given by = 0.2, 0.3). On the other hand, γ t (7.32) increases almost linearly with the permanent positive recharge perturbation at the downgradient well (in blue) but the maximum change of γ t (7.32) with respect to the baseline recharge scenario is only around 15%. Decreased recharge (Figure 7.33b) also results in a non-linear response ofγ t to. Compared to thepositiverechargeshiftscenarioinwhichthesource-zonewellrespondsfirst(Figure7.33a),when the recharge is negatively shifted (Figure 7.33b), the downgradient well responds first. A decrease ofthenaturalrechargedoesnotsignificantlyinfluencethetimewhenthe 3 Hconcentrationexceeds MCL at the source-zone well (in red) given that the minor concentration increase takes place when the baseline concentration exceeds the MCL. However we identify the value of = 0.2 as the trigger level of recharge decrease after which the time above MCL is significantly affected at the downgradient well (in blue). The remaining recharge steps produce less significant changes at the downgradient well. One-year of extreme recharge produces relatively smaller changes in the time above MCL as compared to the effect of permanent shifts of recharge when < 1 (compare Figure 7.33c with 347 Figure 7.33a). This scenario does not influence much the time of MCL exceedance at the source- zone well (in red), where a maximum γ t around 20% is identified when = 4. Nevertheless, the response at the downgradient well, also non-linear, shows a trigger level represented by = 1 above which the time of exceedance of MCL is highly impacted. This significant impact takes place because the increase of concentration downgradient, which results in values above the 3 H MCL, happens later in time, compared to the upgradient well, when the baseline concentration is already below the MCL. When the capping structure fails (Figure 7.33d), we notice thatγ t reaches the maximum value at both wells, similarly to what observed when looking at the γ c - results. The value of γ t then decreases at both wells when the recharge increases in the presence of no cap. This decrease is almost linear at the source-zone well. In the overall, we observe that the capping failure is identified as a trigger condition which causes an important increase of both EPMs (as quantified by γ c and γ t ), therefore a substantial risk increase. When comparing the results in Figures 7.32 and 7.33 we understand that different recharge’s thresholds/trigger levels are identified depending on the metric of interest for environmental reme- diation (for instance the maximum concentration/export BTC’s change or the time of exceedance of MCL’s change caused by the recharge perturbations), the location of the observation location and the measured variable (e.g., contaminant concentration or contaminant export). 7.7 Summary and Conclusions In this paper we investigate the impact of climate change-driven aquifer’s recharge shifts on residual contaminants in soil and groundwater subject to sustainable remediation. For the sake of simplicity we assume that the changes in precipitation and temperature can be translated into changes of the natural aquifer’s recharge. We establish the evaluation methodology, including the development of potential future scenarios on the basis of national climate assessments, site-specific 348 model developments, and evaluation criteria. We illustrate four scenarios characterized by a range of variable recharge values: (1) constant positive recharge shift after a certain year (2) constant negative recharge shift after a certain year, (3) one-year of extreme recharge, (4) cap failure and constant positive recharge shift. Our methodology is demonstrated by simulating the 3 H plume migration within the US DOE’s nuclear reservation Savannah River Site F-Area. We employ the unsaturated-saturated flow model Amanzi, calibrated and verified using site data. In summary, our results generally show that changes of the recharge regime (even small) can significantly affect contaminant concentrations. The most noticeable outcome is the concentration “rebound”effecttakingplace, afteraninitialslightdilution, underhigherrechargeand/orcapping failureconditions, whichisgivenbythemobilizationofthecontaminantmassfromthesourcezone and its transfer to the aquifer system. The concentration “rebound” effect is more pronounced and happens earlier as the recharge perturbation increases. Decreased recharge conditions could also cause a small concentration increase attributed to a decrease in the dilution potential of the aquifer. On the other hand, the 3 H export at the CP is only minimally influenced by recharge shifts, except for some extreme recharge scenarios. Trigger levels of recharge which highly impact the concentration BTC at the wells are identified. These threshold/trigger levels depend on the observation location and on the EPM under investigation, quantified through γ c and γ t . The latter display a non-linear response to the perturbation. For example, it is interesting to observe that the most significant influence on the time of exceedance of MCL under positive recharge shifts (i.e., first scenario) is identified under recharge values in the middle of the range considered given that higher values cause a steeper decrease of the concentration after the “rebound” phase, i.e., the concentration reaches the MCL faster. Our results suggest that close monitoring of wells concentrations should be adopted during precipitation (connected to higher aquifers’ recharge) and drought (connected to lower aquifers’ recharge) periods, however the actual risk of the downgradient population, as quantified through the export, could be under control even when the well concentrations are remarkably impacted. 349 Our findings constitute then immediate contribution to guide groundwater monitoring in the presence of increased climatic variabilities, particularly in explaining concentration anomalies. For example, without a proper understanding, the concentration increase, due to higher or lower recharge, could be mistaken, for instance attributed to additional leaks or contaminant sources. In addition, ouranalysisindicatesthatsource-zonewellsarecriticaltoearlydetectmobilizedresidual contaminants under increased recharge or cap failure conditions. It would be advantageous to have more frequent sampling or in situ monitoring in the proximity of the contaminant source zones as an early warning system [3;3]. Aside from monitoring the contaminant concentration, characterizing the hydraulic fluxes in the vicinity of the source zone can also aid in understanding the macro-dispersive behavior of the plume and corresponding risks to the environment and public health [65,322]. Our work also emphasizes the importance of properly maintaining the capping structurenotonlytosequesterresidualcontaminantsbutalsotoreducetheuncertaintyassociated with climate variability, in fact the difference among recharge scenarios is smaller when the cap is undamaged. We finally highlight that currently, simplified models adopting conservative assumptions are often used for performance and risk assessments at contaminated sites. Conservative approaches usually assume the worst case scenario, associated with higher recharge rates or higher per- meability values to increase the plume mobility. The modeling scenarios investigated in this work, however, call into question the appropriateness of such conservative approaches, in fact the non-trivial trade-offs arising from the interplay between dilution and contaminants’ mobilization require the use of more realistic and accurate flow and transport simulations, achieved through proper calibration processes, as well as probabilistic risk assessments [8,66,344,357,3;9]. Overall, the trade-offs identified in our work must be evaluated with respect to the specific time, location and performance metric under investigation. Understanding these trade-offs would enable better allocation of available resources towards reducing uncertainty in decision making [5:,67]. 34: Our work could be expanded by considering a more complex geochemistry setting, as well as surface water processes (e.g., evapotranspiration), and land model components (e.g., vegeta- tion) of the SRS F-Area. Additionally, this study would benefit from an inclusion of different conceptualization of the aquifer’s heterogeneous properties together with remediation strategies located at the site (e.g., pump and treat, funnel and gate systems). In particular, the effect of geological heterogeneity should be explored in more detail given that it can significantly augment contaminant plume mixing rates [63]. Supplementary material As explained through the paper, given the intrinsic uncertainties of climate predictions, our study considers a wide range of natural recharge values. Even so, the smallest time scale of recharge shift considered is one year, in the third scenario. Therefore we decided to simulate flow and transport for a fifth scenario characterized by a smaller time scale of recharge variability, i.e., 3 months (one trimester). In this case, the perturbed recharge, indicated as R P , is expressed as: R P (t) =    R B , if t<t ∗ (7.33) R B (1 +) n , if t ∗ ≤t≤t ∗ + 1 (7.34) In equation (7.34), always indicates the recharge perturbation and t ∗ = 2020, whereas n indicates the trimester number within year 2020, i.e., n = 1 for the first trimester of 2020, n = 2 for the second trimester of 2020, n = 3 for the third trimester of 2020 and n = 4 for the fourth trimester of 2020. Therefore the fifth scenario simulates a positive shift of recharge by a factor every three months within year 2020, starting in 2020 and ending in 2021, indeed the recharge is set back to baseline values starting in year 2021. The recharge shift by a factor in 2020 is applied with respect to the recharge value of the precedent trimester and it starts from baseline recharge in 2020. The perturbation assumes the following range of values: = [0.1, 0.2, 0.3, 0.4, 0.5]. 34; 2020 2040 2060 2080 2100 Time [y] 1 2 3 4 5 6 7 C [mol kgw -1 ] 10 -12 Positive recharge shift (a) Source-zone well =0 =0.1 =0.2 =0.3 =0.4 =0.5 2020 2040 2060 2080 2100 Time [y] 0.5 1 1.5 2 2.5 3 C [mol kgw -1 ] 10 -12 (b) Downgradient well =0 =0.1 =0.2 =0.3 =0.4 =0.5 2020 2040 2060 2080 2100 Time [y] 5 10 15 Export [mol y -1 ] 10 -6 (c) Control plane =0 =0.1 =0.2 =0.3 =0.4 =0.5 Figure7.34: Temporalevolutionof 3 Hconcentrationatthe: (a)source-zonewell,(b)downgradient well and (c) 3 H export at the CP for the baseline scenario, indicated by = 0, and the positive three months recharge shift scenario, characterized by = 0.1, 0.2, 0.3, 0.4, 0.5. The results of different recharge perturbations are portrayed in different colors. The thin horizontal dashed black line represents the MCL for 3 H. 352 The temporal evolution of the 3 H concentrations [mol kgw −1 ] at the observation wells and of the 3 H export [mol y −1 ] at the CP for the fifth scenario within the timeframe 2020-2100 is presented in Figure 7.34. Under positive three months recharge shifts (Figure 7.34), the concen- tration initially slightly decreases for approximately 2− 8 years at the source-zone well (Figure 7.34a) and for around 8− 10 years at the downgradient well (Figure 7.34b) because of dilution, similarly to what observed for the first scenario (i.e., under a constant positive recharge shift). Also in this case we identify a concentration “rebound” phase afterwards, during which the con- centration increases at both wells (Figure 7.34a-b). At the source-zone well, under the recharge shifts indicated by = 0.3− 0.5 we notice a subsequent decrease of concentration followed by a concentration “rebound” and so forth, caused by the different recharge increases during year 2020. We also observe that, as the perturbation is higher, the initial concentration decrease is generally shorter and more pronounced and the subsequent concentration peak is higher (Figure 7.34a-b), however when = 0.1− 0.2 only minor effects with respect to the baseline concentration are identified, especially at the downgradient well. When comparing Figure 7.34a and Figure 7.34b we observe that the BTC’s pattern (i.e., slight initial decrease followed by a concentration rebound) is less pronounced and happens later in time in Figure 7.34b according to the fact that the downgradient well is further away from the source. In line with the outcomes of the other recharge scenarios, the recharge perturbation only causes minor changes of the contaminant export to the creek (Figure 7.34c). The results of the fifth scenario confirm that the key findings of our study hold even when the recharge shifts are more frequent in time, which could represent a more realistic situation. In fact, we generally observe that, even when the time variability of the recharge shift is shorter (i.e., three months instead of a year or more years), the positive shift of recharge triggers a pattern in the concentration BTC characterized by a slight decrease followed by a concentration “rebound”phase. Thisconfirmsthat,nomatterwhatisthetimescaleoftherechargeperturbation, increasing aquifers’ natural recharge initially causes a slight decrease of the concentration followed 353 by a concentration rebound effect. Indeed, the higher volume of water that enters the system initially results in dilution, and later produces an increase of the concentration attributed to the mobilization of residual 3 H under additional water. The upcoming chapter presents ongoing research aiming at improving the prediction of the required residence time of reclaimed water into aquifers during MAR operations and at studying the connected uncertainties. 354 Chapter 8 Quantification of solute arrival time uncertainty in the planning of Managed Aquifer Recharge The content of this chapter constitutes work in progress that will lead to a journal paper. 8.3 Introduction Itiswellestablishedthatgroundwaterconstitutesaprimarysourceofwaterworldwideandaround 40% of California’s annual water supply [82]. Quantifying the availability of subsurface water is not trivial and often aquifers are subject to critical conditions of overdraft (i.e., more water is withdrawn from the aquifers than can be replenished by precipitation and other means). Indeed, in the state of California, USA, the annual rate of overdraft is estimated to be 1− 2 million acre-feet and there is concern that overdraft has grown in the last years because of drought, espe- cially in the Central Valley region [444]. Overdraft can lead to significant adverse environmental and economic impacts, including increased extraction costs, costs of well deepening or replace- ment, land subsidence, water quality degradation. California’s over drafted basins must achieve groundwater sustainability by 2040 or 2042 according to the 2014 Sustainable Groundwater Man- agement Act, released by the Department of Water Resources (DWR). To better illustrate this point in August 2015, DWR released a list of 21 basins subject to overdraft. DWR is currently 355 defining sustainable management and establishes different criteria to achieve the balance between extraction and recharge of aquifers in the long term [354]. To overcome challenges associated with groundwater overdraft novel solutions have emerged, for instance, effective conjunctive use of water sources (i.e., surface water and groundwater) which isfundamentaltoachievesustainabilityofwaterresources. Conjunctiveuseisdonebysubstituting surfacewaterwithgroundwater. Inordertomanageclimateextremes(i.e., longtermdroughtsand intense floods), that are projected to increase in the future due to climate change, conjunctive use ofgroundwateriscomplementedwithmanagedaquiferrecharge(MAR)thatconsistsinrecharging subsurface aquifers with surface water. MAR can be seen as an extension of conjunctive use because, instead of substituting surface water with groundwater, it recharges and stores surface water into the subsurface for subsequent recovery or environmental benefits. Between the water sources for MAR we find stormwater, rivers and lakes water, reclaimed and desalinated water and also groundwater located in other aquifers [422]. MAR presents the advantage of low energy requirements and natural purification that occurs because of subsurface processes like mechanical filtering, sorption and biodegradation [422]. The use of recycled or reclaimed water for aquifer recharge or other purposes is often adopted and is becoming increasingly popular worldwide [427]. However, there are multiple water quality issues connected to this practice. Indeed, reclaimed wastewater may introduce organic or micobiological contaminants in the subsurface and the latter maybecapturedbyproductionwellsatalaterstage. Betweenthecontaminantsofconcernwefind emerging contaminants, which are pharmaceuticals and other organic wastewater contaminants (OWCs) documented in water resources around the world [32] and not yet regulated. The state of California follows a series of regulation to plan MAR practices utilizing reclaimed water and subsequent water extraction [427]. As explained above, reclaimed water injected in the subsurface undergoes transformation processes such as dilution and natural degradation which may decrease contaminant concentrations. The state of California established a residence (or retention) time of all reclaimed water inserted into an aquifer before reaching production wells. 356 This established residence time ensures that degradation and dilution mechanisms occur, as well asawaterproductionofonlyaspecifiedfractionofthereclaimedwater(regardlessoftheresidence time). Finally it is established that the extraction wells must be located at a certain distance from the reclaimed water recharge area. At the current stage, engineeers and hydrogeologists rely on numerical models or tracer ex- periments in order to establish the residence times in reclaimed water systems. For example, it is permissible to use simple seepage velocity based on Darcy’ s law calculations or determinis- tic numerical models during the planning stages of the MAR project in order to receive initial regulatory approval. The residence time estimates are then subject to certain safety factors to offset their uncertainty. The deterministic numerical models employed are either groundwater flow models that only account for water velocity or simple transport models that also incorporate advection and dispersion. Within three months of the aquifer recharge project’ s initial operation (or beforehand, if feasible) an empirical tracer test, using either an added (artificial) or intrinsic tracer, is performed to establish/verify the arrival time estimates and the resulting “boundary zones” where drinking water well construction would be prohibited [423]. The retention time shall be the time representing the difference from when the water with the tracer is applied at the groundwater replenishment reuse project to when either: 2% of the initially introduced tracer concentration has reached the downgradient monitoring point, or 10% of the peak tracer unit value observed at the downgradient monitoring point reached the monitoring point [423]. Cur- rent regulations on the residence time of reclaimed water are then based on deterministic studies, thus neglecting uncertainty, and often consider simplified conditions, for example they do not properly account for the heterogeneity of subsurface hydrogeological properties. Additional sci- entific studies are needed to refine the current regulations, in particular stochastic models that inglobe hydrogeological heterogeneity could quantify the uncertainty connected to the predicted residence time. 357 Planning MAR operation is challenging because of different factors: the natural variability of the subsurface hydrogeological properties, the presence of external factors as climatic variability or engineered variability of production operations. It is well known that the heterogeneity of hydrogeological properties controls contaminants’ arrival times [37,58,5:,66,:6,3:4,3:6,3::]. Thehighcontrastofpermeabilityvaluesleadstothepresenceoffastflowchannels(e.g., connected structure) which control first time arrivals [3:2] and the overall dispersive behaviour of a potential solute plume stemming from a MAR operation. The goal of this chapter is to quantify the uncertainty of reclaimed water travel times under different climatic and engineered scenarios which would help better planning and designing MAR systemsthat, atthecurrentstate, considerafixedresidencetimeofthewaterinthesubsurfacebe- foreextraction. Weaimatansweringthefollowingresearchquestions: whatistheprobabilitythat the travel time of a solute is below the established threshold? What is the uncertainty connected with the determination of the residence time? How are different factors (e.g., climatic variability, water pumping) occurring in MAR operational areas affecting arrival times’ uncertainty? 8.4 Model formulation We consider a three-dimensional unconfined regional aquifer identified by a Cartesian coordinate system x=(x, y, z) and illustrated in Figure 8.3. The aquifer is characterized by spatially hetero- geneous, locally isotropic, hydraulic conductivity K and uniform porosity φ. An areal recharge attributed to precipitation is applied to the top of the domain and managed aquifer recharge by means of reclaimed water is performed through a perculating basin. The boundary conditions are such that the base flow is uniform-in-the-mean and at steady-state. Five pumping wells located 358 Figure 8.3: Conceptual model of the problem under investigation: a three-dimensional unconfined aquifer subject to a distributed areal recharge and managed aquifer recharge (MAR) through a perculating basin. Five extraction wells are operating. Emerging contaminants move along the flow direction and reach a control plane located downgradient. at x i , i = 1,.., 5, are operating at a fixed extraction rate Q w . Flow is governed by the following partial differential equation: ∂ ∂x h K(x) ∂h(x) ∂x i + ∂ ∂y h K(x) ∂h(x) ∂y i + ∂ ∂z h K(x) ∂h(x) ∂z i − 5 X i=1 Q w δ(x−x i )−R 0 = 0, (8.3) and the specific discharge field is given by Darcy’s law: q(x) = −K(x)∇h(x), (8.4) where q indicates the specific discharge (m d −1 ), h (m) is the hydraulic head, K (m d −1 ) is the hydraulic conductivity and R 0 (d −1 ) is the aquifer recharge per unit length. As explained in Section 8.3 reclaimed water contains contaminants that could create potential human health and environmental threats. We then assume that the contamination source zone coincides with the MAR basin and that non-reactive contaminants are instantaneously released at time t 0 within the MAR area. Contaminant transport is governed by the advection-dispersion equation: ∂C(x,t) ∂t −∇· [D∇C(x,t)−v(x)C(x,t)] = 0, (8.5) 359 where C indicates the concentration of mobile solute, v is the Darcy-scale velocity (i.e., v = q/φ) and D is the dispersion tensor assumed to be uniform. The following Section presents the numerical implementation followed to solve equations (8.3)-(8.4) and (8.5). 8.5 Numerical implementation For our simulations, we consider an aquifer with dimensions 10000 m× 6000 m× 100 m along x, y and z respectively. Given that equations (8.3)-(8.5) represent partial differential equations characterized by spatially variable coefficients, the flow and transport problems are solved nu- merically. For such reasons, our physical domain is discretized into 100× 60× 50 cells for a total of 3× 10 5 numerical cells. The top of the domain is subject to a uniform recharge rate R = 8.64× 10 −5 m d −1 (i.e., 3 cm yr −1 ), which correspond to typical values encountered in southern California [427]. The MAR operational basin has dimensions of 1600 m× 800 m along the longitudinal and transversal direction and its centroid is located atx=(2950 m, 2950 m, 0 m). Five extraction wells located at a fixed longitudinal distance from the source zone, transversally equi-spaced and aligned with the source zone centroid are operating at a fixed pumping rate (see Figure 8.3). We apply no-flow boundary conditions at the lateral side edges and at the bottom of the domain whereas we consider depth-dependent fixed hydraulic heads along the upstream and downstream boundaries that lead to a mean longitudinal gradientJ. The areal recharge through the top part of the domain and the head differential between upstream and downstream edge create then a “diagonal” regional flow from top-upstream to bottom-downstream. Due to lack of a detailed site characterization, the hydraulic conductivity is regarded as a Random Space Function (RSF) [3:4]. In our work, the log-conductivity, namely Y = ln(K) is characterized by its mean valuehYi, its variance σ 2 Y and correlation scales λx, λy and λz and the spatial covariance of Y is modeled through an exponential correlation model. Values of the geostatistical parameters are provided in Table 8.3. Given the uncertainty in the Y field, we 35: Symbol Significance Units Values L x , L y , L z Domain Size m 10 4 , 6× 10 3 , 100 Δ x , Δ y , Δ z Cell dimension m 100, 100, 2 J Mean longitudinal gradient - 2× 10 −3 hKi Average hydraulic conductivity m d −1 2 λ x , λ y , λ z Correlation length of Y m 800, 800, 16 σ 2 Y Variance of Y - 1, 0.25 φ Porosity - 0.25 R Areal recharge rate m d −1 8.64× 10 −5 Q w Pumping rate m 3 d −1 500, 1000 α x Longitudinal dispersivity m 0.01 α y Transverse dispersivity m 0.001 Table 8.3: Model discretization and physical parameters employed in the study. cast our problem in a statistical Monte Carlo (MC) framework. We randomly generate each Y field using a sequential Gaussian simulator [389]. For each Y field, the flow problem equations (8.3)-(8.4) are solved through the finite-difference model MODFLOW 2005 [;4] and the transport equation (8.5) through the GPU accelerated Random Walk Particle Tracking code PAR 2 [3:3]. We performed convergence analysis to set the optimal number of inserted contaminant particles as well as the number and duration of the simulations’ time step (δt). We show in Figure 8.4 that the mean temporal evolution of the fraction of particles recovered at a control plane (CP) located at the same x-coordinate of the pumping wells over 500 MC simulations is the same with 10 4 inserted particles andδt=100 days and with 10 5 particles andδt=1 day. Based on the results displayed in Figure 8.4, we decided to insert 10 4 particles with a numerical time stepδt = 100 d. The total simulation time is 4× 10 5 days, around 1000 years, which is the time taken to recover 95% of the solute particles injected. Each particle carries a fraction of the solute mass injected. Particles are moved according to advection and local-scale dispersion. The velocity field obtained through the solution of equations (8.3)-(8.4) is employed to compute advection while local-scale dispersion is simulated by adding a random displacement to the particle. For additional details on the solution of the transport equation (8.5) please refer to Rizzo et al., 423; [3:3]. The inserted 35; 0 100 200 300 400 500 600 700 800 900 1000 time [yr] 0 20 40 60 80 100 fraction of particles [%] Figure 8.4: Mean temporal evolution of the fraction of particles recovered at the control plane over 500 MC simulations for time step (δt) of 1 day and 10 5 inserted particles (blue curve) and for time step of 100 days and 10 4 inserted particles (orange circles). For this comparison we considered K fields characterized by σ 2 Y = 1, recharge R, and Q w =1000 m 3 d −1 applied to the five extraction wells. particles are recovered at a CP located atx = 6950 m. Table 8.3 reports the parameters adopted to numerically solve groundwater flow (equations 8.3-8.4) and transport (equation 8.5). We analyze the outcomes of the MC simulations by computing arrival times of particles at the CP. We consider different scenarios characterized by different heterogeneities of the natural system as well as different external factors (climatic or engineered). In particular we consider both σ 2 Y = 1 and σ 2 Y = 0.25, the areal recharge value presented in Table 8.3, indicated as R, as well as a decrease of recharge of one order of magnitude, representative of drought conditions, and two values of extraction rate for the five extraction wells: 500 m 3 d −1 and 1000 m 3 d −1 . As illustrated before (see Section 8.3), our goal is to better understand the uncertainty con- nected with the design of MAR operations. In particular we compute the time of arrival of 2% of the solute mass (t 2% ), following the tracer test regulations of the State of California [423], 5% of the solute mass (t 5% , identified as early arrival time), 50% of the solute mass (t 50% ) and 95% of the solute mass (t 95% , identified as late arrival time) at the CP of interest. Because of the randomness in Y, contaminants’ arrival time is a random function that must be characterized statistically by its probability density function (PDF) and corresponding cumulative distribution function (CDF). The CDF provides values of risk reliability and risk exceedance (through the 362 complementary CDF) with respect to a regulatory threshold value, e.g., the fixed residence time of reclaimed water in the subsurface established by a water management agency, which corre- sponds to the minimum time required to achieve natural purification of the reclaimed water. For example, risk, identified by ξ, is defined in terms of the probability that the travel time of a con- taminant from the source zone to the CP is below a threshold established by a water management agency: ξ =P [t≤t ∗ ], (8.6) where t ∗ is the fixed residence time and P denotes probability. Moreover, we show the impact of different natural variability, climatic and engineered con- ditions on the travel time uncertainty by means of the travel time probability density function (PDF). The results of our computations are presented in the following Section. 8.6 Computational results 8.6.3 Residence time reliability In Figure 8.5 we present the CDF of dimensionless arrival time of 2% of the contaminant mass (i.e., t 2% ) at the CP location for σ 2 Y = 1, recharge R = 8.64× 10 −5 m d −1 , and Q w = 1000 m 3 d −1 for the five extraction wells. All results described in the following refer to extraction wells operating at Q w = 1000 m 3 d −1 unless otherwise specified. Figure 8.5 serves as an illustrative example to show the importance of quantifying the reliability (see equation 8.6) of an established time threshold (e.g., the residence time t ∗ ), information not considered during the planning of MAR operations at the current state. For comparison, we include the computational results ob- tained for a deterministic aquifer characterized by a homogeneous K field. For our homogeneous case, we set K = 2 m d −1 (i.e., average conductivity value). By assuming homogeneous hydro- geological properties, we obtain a dimensionless t 2% of 1.68 (see Figure 8.5), The first arrivals 363 0 1 2 3 4 5 6 tU/ 0 0.2 0.4 0.6 0.8 1 CDF Figure 8.5: Cumulative distribution function (CDF) of dimensionlesst 2% forσ 2 Y = 1, rechargeR, and Q w = 1000 m 3 d −1 . CDF, displayed in Figure 8.5, shows that the probability of the dimensionless t 2% being lower than 1.68 is approximately 90% (i.e., ξ≥ 0.9). Therefore there is a relatively high probability that the contaminant will reach the extraction wells’ location after the established residence time necessary to achieve natural purification. Additional information on the statistical characterization of arrival times under different sce- narios are provided in the next Section. 8.6.4 ImpactofarealrechargeandY varianceonarrivaltimesuncertainty Figure 8.6 presents the PDF of t 5% (a), t 50% (b) and t 95% (c) for standard recharge conditions of southern California,R (see Table 8.3), for two levels ofY-field heterogeneity, namelyσ 2 Y = 1 and σ 2 Y = 0.25. The analysis is repeated and shown in Figure 8.6d-f for a lower recharge characterized by one order of magnitude less, i.e., for R = 8.64× 10 −6 . Contaminant plumes traveling in geological formations displaying larger variability in the Y- field, epitomized by the Y variance, will display increased macro-scale spreading (i.e. macrodis- persion). The latter implies an increase of the probability of finding faster first arrival times, as well at t 50% . This effect is emphasized in the presence of higher areal recharge (compare Figure 8.6a-c with Figure8.6d-f obtained for a lower recharge). Indeed, increasing recharge leads to more water in the system which changes the hydraulic gradient and results in a faster mobilization of 364 0 100 200 300 400 500 600 700 800 900 t 5% [yr] 0 0.5 1 1.5 2 2.5 PDF 10 -5 R=8.64 10 -5 m d -1 2 Y =1 2 Y =0.25 0 100 200 300 400 500 600 700 800 t 5% [yr] 0 0.5 1 1.5 2 2.5 PDF 10 -5 R=8.64 10 -6 m d -1 2 Y =1 2 Y =0.25 0 200 400 600 800 1000 1200 t 50% [yr] 0 0.5 1 1.5 2 2.5 PDF 10 -5 2 Y =1 2 Y =0.25 0 100 200 300 400 500 600 700 800 900 1000 t 50% [yr] 0 0.5 1 1.5 2 2.5 PDF 10 -5 2 Y =1 2 Y =0.25 0 200 400 600 800 1000 1200 t 95% [yr] 0 0.5 1 1.5 2 2.5 PDF 10 -5 2 Y =1 2 Y =0.25 0 200 400 600 800 1000 1200 t 95% [yr] 0 0.5 1 1.5 2 2.5 PDF 10 -5 2 Y =1 2 Y =0.25 Figure 8.6: Probability density function (PDF) oft 5% (a),t 50% (b),t 95% (c) for higherY hetero- geneity (σ 2 Y = 1, continuous curves) and lower Y heterogeneity (σ 2 Y = 0.25, dashed curves) and recharge conditions R. PDF of t 5% (d), t 50% (e), t 95% (f) for higher Y heterogeneity (σ 2 Y = 1, continuous curves) and lower Y heterogeneity (σ 2 Y = 0.25, dashed curves) and lower recharge conditions. contaminants. Similar effects are observed on t 95% , where a wider PDF is observed for σ 2 Y = 1 when opposed to the results obtained fromσ 2 Y = 0.25. The effect of heterogeneity is enhanced for lower recharge conditions (Figure 8.6d-f). By comparing 8.6a-c to 8.6d-f, we generally observe that the impact of the variance in Y (i.e., heterogeneity) on arrival times is emphasized under lower recharge conditions, meaning that different climatic inputs could emphasize or smooth the impact of natural variability on transport, as also shown in the previous chapters. 365 0 200 400 600 800 1000 1200 arrival time [yr] 0 0.5 1 1.5 2 PDF 10 -5 Q w =1000 m 3 /d t 5% t 50% t 95% 0 200 400 600 800 1000 1200 arrival time [yr] 0 0.5 1 1.5 2 PDF 10 -5 Q w =500 m 3 /d t 5% t 50% t 95% Figure 8.7: Probability density function (PDF) oft 5% ,t 50% andt 95% forQ w =500 m 3 d −1 (a) and Q w =1000 m 3 d −1 (b) and σ 2 Y = 1. 8.6.5 Impact of well pumping rate on arrival time uncertainty Figure 8.7 presents the PDF oft 5% ,t 50% andt 95% under fixed recharge conditions and for pump- ing wells operating at a pumping rate of Q w =500 m 3 d −1 (a) and Q w =1000 m 3 d −1 (b). As expected, the probability of finding faster arrival times increases with larger values of Q w . Hy- draulic gradients are enhanced with an increase of pumping rates which in turn leads to faster arrival times. This is particularly observed for t 5% and t 50% . The effect of Q w on t 95% is not so pronounced as shown in Figure 8.7. This is due to the fact that the trailing edge of the plume, associated with t 95% , is mainly controlled by the continuity of low conductive areas and therefore lower velocities. We also observe that higher pumping rates decrease the spread of the travel time PDF, in particular for t 5% and t 50% . 8.7 Summary In summary, our preliminary results show the contribution that stochastic modeling of flow and transport processes can add to the design of MAR operations, in particular in increasing the 366 reliability of the residence time estimation and in quantifying the uncertainty of the estimation. MAR operations often employ reclaimed water, which is treated to a certain standard but may contain residual contamination, for example emerging contaminants, which are often times not regulated. Giventhepotentialharmthatcontaminatedwastewatercancreate,watermanagement agencies establish a residence time of water in the aquifer before successive extraction in order to ensure that the natural degradation mechanisms that contaminants undergo while migrating into aquifers occur, and that the quality of the extracted water respects environmental standards. Potential negative consequences of withdrawing contaminated reclaimed water at a later stage are the increase of water treatment costs as well as bad effects on human health, in particular in the presence of emerging contaminants, whose maximum contaminant level has not yet been established by environmental protection agencies. Therefore, it is very important to ensure that the extracted water is not contaminated, i.e., that there is sufficient time between the injection of reclaimed water into aquifers and the extraction. The estimation of the residence time of reclaimed water underground is then a key factor. Water managers establish the residence time on the basis of simple calculations based on Darcy’s law as well as deterministic groundwater flow models and simple contaminant transport models. To account for the degree of uncertainty they employ safety factors and they perform tracer tests. Current regulation then neglect uncertainty stemming from the incomplete characterization of the aquifers’ hydrogeological properties, and consider simplified conditions. In our work we simulated the migration of reclaimed contaminated water inserted in the subsurface through an aquifer model that could be applied to an aquifer located in southern California and subject to MAR operations. We considered a regional model subject to realistic recharge and pumping rates. The contribution of our study consists in quantifying the impact of subsurface heterogeneity variability, as well as of variability of recharge conditions and pumping ratesontheuncertaintyconnectedwitharrivaltimesofreclaimedwatertoanobservationlocation. The latter constitute important information to reinforce the current estimation of the residence 367 time. Our study considers the presence of realistic external factors, whose magnitude can change and aims at providing risk estimates. The work presented in this chapter is currently ongoing and will be expanded. We are now working on computing the contaminant concentration at the five extraction wells and comparing it to initial concentration values to understand whether or not the probability that the retention time is sufficient. The analysis on the concentration field is also useful to estimate the dilution potential of the MAR system. 368 Chapter 9 Conclusions Groundwater contamination constitutes a severe issue worldwide. Challenges connected with un- derstanding and modeling groundwater flow and contaminant transport processes are connected to different factors. First of all, hydrogeological properties of aquifers (e.g., hydraulic conductivity and porosity) are heterogeneous and a detailed spatial characterization is unfeasible because of the high cost associated with collecting subsurface data and the limited resources available. More- over, aquifers are stressed by multiple external factors, both stemming from natural or engineered variability. Between those external factors we find changing climatic conditions, groundwater pumping and groundwater recharge. In the overall, the principal elements controlling subsurface contaminant transport in our society are enclosed in: natural variability of the subsurface envi- ronment, engineered variability (i.e., water management practices as pumping and recharge) and climatic variability. The novelty of this thesis lies in considering the interaction between these three different areas and its effect on contaminant transport. To the best of our knowledge, most works in the literatureregardtheproblemasasinglesystemwhereasourstudyconsidersa“multi-system”type of approach by investigating multiple potential scenarios in which different elements stemming from the three above areas are keyed together and affect the risk originating from subsurface contamination. The following paragraphs will conclude this thesis by briefly summarizing the 369 main findings of each chapter on published research or ongoing work and finally highlighting the key contributions of our work in modern day society. The second chapter presents the impact of the coupled spatial variability in porosity and hydraulic conductivity on contaminant transport and probabilistic risk analysis. Previous works in the literature focused solely on the effects of conductivity variability on flow and transport. Inglobing the spatial heterogeneity of the porosity field, positively correlated to the hydraulic conductivity field, in contaminant transport predictions affect both arrival times and peak mass fluxes at observation locations as it decreases the macrodispersion of the plume leading to later first arrival times, earlier late arrivals and increased peak mass fluxes, which correspond to higher health and environmental risks. We also observe that the relevance of incorporating the spatial variability of the porosity can depend on the metric of interest: this relevance is emphasized when investigating contaminants’ arrival times over peak mass fluxes at the observation locations. Additionally the impact of porosity variability on mean increased lifetime cancer risk caused by exposure to chlorinated solvents is only minor but we identify an important impact on cancer risk uncertainty, which is reduced when the heterogeneity in porosity is included in transport models. The third chapter explores the effect of multiple temporally variable pumping regimes on key stochastic features of contaminant transport in randomly heterogeneous aquifers. This chapter gives insights on a realistic approach for the evaluation of the risk associated with contamination of groundwater extracted at operating wells. We show that transient pumping strategies can markedly affect the temporal evolution of contaminant concentration breakthrough curves and statistical moments by inducing multiple peaks of the concentration signal at the well, according to the pumping signal. Additionally, activating or increasing water withdrawal from pumping wells induces contaminant dilution with fresh water at the well which leads to a concentration decrease and therefore lowers the associated risks. The uncertainty of the concentration observed at the pumping well is also controlled by the selected pumping schedule, and this uncertainty shows a remarkable increase for highly transient pumping schedules. Finally, risk at the well, 36: given by the probability to exceed maximum contaminant levels established by environmental protection agencies, shows a temporally oscillating behavior under transient pumping. We also show that both aquifer heterogeneity and Péclet number control the magnitude of risk whereas the pumping schedule dictates the temporal evolution of the risk. The fourth chapter investigates the impact of the modeling approach employed to describe the random spatial heterogeneity of the aquifer hydraulic conductivity on the statistics of the solute concentration at pumping wells. In particular while simulating the hydraulic conductivity realizations we consider both usually employed Gaussian geostatistical models and non-Gaussian models. The latter are connected with increased probability of finding well connected zones of high conductivity (fast-flow channels) or low conductivity zones. This study shows the interplay between the role played by the model employed to represent the natural variability of hydrogeo- logical properties and the contribution given by the engineered pumping variability on statistics of concentration breakthrough curves. We indicate that the temporal pattern of the concentration signal at an extraction well is controlled by the variability of the extraction schedule, whereas the geostatistical model employed to picture the conductivity field controls the magnitude of the concentration breakthrough curve. In particular, as the departure of the non-Gaussian models from a Gaussian behavior increases, we observe higher mean concentration values. Our results also highlight the existence of a threshold value of a geostatistical parameter above which the effects associated with the non-Gaussian nature of the heterogeneous conductivity structure are undetectable in the concentration breakthrough curve observed at the well. The fifth chapter presents a real case study investigating the impact of climate change-driven aquifer’s recharge shifts as well as failure of remediation approaches (specifically the capping sys- tem) on residual subsurface tritium contamination, at the US DOE’s nuclear reservation Savannah RiverSiteF-Area(SouthCarolina, USA),subjecttosustainableremediationtechniques. Weshow that changes of the recharge regime (no matter how small) can significantly affect contaminant concentrations. In particular, higher recharge and/or capping failure conditions cause an initial 36; slight dilution followed by a significant concentration increase motivated by the high mobilization of contaminant mass from the source zone and its transfer to the aquifer system as the amount of water in the system is augmented. The influence of decreased recharge conditions is less relevant compared to increased recharge and is manifested through a slight increase of the concentration attributed to a decrease in the dilution potential of the aquifer. On the other hand, the tritium export outside the site’s boundary is only minimally influenced by recharge shifts, except for some extreme recharge scenarios. Our analysis identified trigger levels of recharge which highly impact the concentration breakthrough curves at the observation wells. The last chapter presents preliminary findings on the importance of stochastic flow and trans- port modeling to improve the planning and design of groundwater artificial recharge operations, often employing reclaimed water. Our conceptual model could be applied to southern Califor- nia’s aquifers subject to managed aquifers recharge practices as it represents a regional aquifer where realistic groundwater recharge and groundwater extraction take place. A key element in the design of managed aquifer recharge is the calculation of the residence time of reclaimed water underground, in order to ensure that the natural degradation of potential residual contaminants in the injected water occur. In our preliminary research we quantify the uncertainty stemming from the lack of characterization of the conductivity field which reflects on the uncertainty in solute arrival times. The latter information is fundamental as it can help better estimate the res- idence time of reclaimed water underground as well as provide probability bounds on the current regulations. We explore multiple areal recharge and hydraulic conductivity representations and quantify the effect on the probability density function of different arrival times. For example, we show that higher variability of the conductivity field results in earlier arrival, effect emphasized under higher recharge conditions, which contribute in a faster mobilization of the plume. In the overall, the research line presented in this thesis has a direct contribution in helping risk and water managers making more informative decision on resources allocation to direct technical and economical efforts towards an optimal management of water resources. Available resources 372 need to be prioritized and assigned in a defensible manner to ensure that unacceptable identified risks are reduced to acceptable levels. The multiple scenarios explored in this thesis, intercon- necting different areas that influence contaminant transport in modern day society, showed cases in which, for instance, some factors could play a major role over others and therefore need proper consideration and resources allocation. For instance costs derived from increasing well pumping rates could be justified by the production of water characterized by lower contaminant concen- trations, which could result in a decrease of water treatment costs. The fact that the pumping operations can control the temporal patterns of risk and might overshadow the impact of the type of aquifer heterogeneity on concentration signals at pumping wells leads to circumstances in which enhanced efforts should be allocated towards an improved optimal planning of the pumping regime as opposed to a detailed characterization of some features of the heterogeneous properties of an aquifer. Moreover, our research help directing monitoring efforts at contaminated site, for example we showed that close and more frequent monitoring of wells concentrations should be adopted during precipitation and drought periods, however a close monitoring of contaminants’ export outside sites’ boundaries, which could constitute the risk pathway to the downgradient population, may not be necessary. Our findings also contribute in explaining contaminant con- centrations anomalies that could emerge under climatic changing conditions. The results of our research show non-trivial trade-offs arising from the interplay between dilution and contaminants’ mobilization which require the use of realistic and accurate flow and transport simulations as well as probabilistic risk assessments. Our stochastic modeling approach highlights the contribution of different factors in increasing or reducing contaminant transport predictions’ uncertainty, whose quantification is a key element. Our research then enables better allocation of usable resources towards reducing uncertainty in decision making (e.g., in risk assessment or aquifer remediation activities) which is well aligned with goal-oriented site characterization frameworks. 373 Reference List [3] W. Aimrun, M. Amin, and S. Eltaib. Effective porosity of paddy soils as an estimation of its saturated hydraulic conductivity. Geoderma, 343(5-6):3;9–425, 4226. [4] R.AndričevićandV.Cvetković. Evaluationofriskfromcontaminantsmigratingbyground- water. Water Resour. Res., 54(5):833–843, 3;;8. [5] R. Andričević and V. Cvetković. Relative dispersion for solute flux in aquifers. J. Fluid Mech. Mechanics, 583:367–396, 3;;:. [6] G. E. Archie. Introduction to petrophysics of reservoir rocks. AAPG bulletin, 56(7):;65– ;83, 3;72. [7] S. Assouline, M. Möller, S. Cohen, M. Ben-Hur, A. Grava, K. Narkis, and A. Silber. Soil-plant system response to pulsed drip irrigation and salinity. Soil Sci. Soc. Am. J., 92(7):3778–378:, 4228. [8] A. L. Atchley, R. M. Maxwell, and A. K. Navarre-Sitchler. Human health risk assessment of CO4 leakage into overlying aquifers using a stochastic, geochemical reactive transport approach. Environ. Sci. Technol., 69(33):7;76–7;84, 4235. [9] J. E. Atkins and E. F. McBride. Porosity and packing of holocene river, dune, and beach sands (3). AAPG Bulletin, 98(5):55;–577, 3;;4. [:] C. D. Atkinson, J. H. McGowen, S. Bloch, L. L. Lundell, and P. N. Trumbly. Braidplain and deltaic reservoir, prudhoe bay field, alaska. In Sandstone petroleum reservoirs, pages 9–4;. Springer, 3;;2. [;] M. Barahona-Palomo, M. Riva, X. Sanchez-Vila, E. Vazquez-Sune, and A. Guadagnini. Quantitative comparison of impeller-flowmeter and particle-size-distribution techniques for the characterization of hydraulic conductivity variability. Hydrogeol. J., 3;(5):825–834, 4233. [32] K. K. Barnes, D. W. Kolpin, E. T. Furlong, S. D. Zaugg, M. T. Meyer, and L. B. Barber. A national reconnaissance of pharmaceuticals and other organic wastewater contaminants in the United States —I) Groundwater. Sci. Total Environ., 624(4-5):3;4–422, 422:. [33] D. W. Barr. Discussion of “Goodbye, Hazen; Hello, Kozeny-Carman” by W. David Carrier III. J. Geotech. Geoenviron. Eng., 353(:):3279–327:, 4227. [34] S. A. Bea, H. Wainwright, N. Spycher, B. Faybishenko, S. S. Hubbard, and M. E. Denham. Identifying key controls on the behavior of an acidic-U ( VI ) plume in the Savannah River Site using reactive transport modeling. J. Contam. Hydrol., 373:56–76, 4235. [35] J. Bear. Hydraulics of groundwater. McGraw-Hill Inc., New York, 3;9;. 374 [36] J. Bear and M. Jacobs. On the movement of water bodies injected into aquifers. J. Hydrol., 5(3):59–79, 3;87. [37] A. Bellin, P. Salandin, and A. Rinaldo. Simulation of dispersion in heterogeneous porous formations: Statistics, first-order theories, convergence of computations. Water Resour. Res., 4:(;):4433–4449, 3;;4. [38] A. Bellin and D. Tonina. Probability density function of non-reactive solute concentration in heterogeneous porous formations. J. Contam. Hydrol., ;6(3):32;–347, 4229. [39] M. Bianchi and C. Zheng. A lithofacies approach for modeling non-Fickian solute transport in a heterogeneous alluvial aquifer. Water Resour. Res., 74(3):774–787, 4238. [3:] J. P. Bloomfield, R. J. Williams, D. C. Gooddy, J. N. Cape, and P. Guha. Impacts of climate change on the fate and behaviour of pesticides in surface and groundwater-a UK perspective. Sci. Total Environ., 58;(3-5):385–399, 4228. [3;] P. C. Carman. Fluid flow through granular beds. Trans. Inst. Chem. Eng., 37:372–388, 3;59. [42] P. C. Carman. Flow of gases through porous media. 3;78. [43] W. D. Carrier III. Goodbye, Hazen; Hello, Kozeny-Carman. J. Geotech. Geoenviron. Eng., 34;(33):3276–3278, 4225. [44] L. C. Chang, C. A. Shoemaker, and P. L. F. Liu. Optimal time-varying pumping rates for groundwater remediation: Application of a constrained optimal control algorithm. Water Resour. Res., 4:(34):5379–5395, 3;;4. [45] Y. Chen, C. Lu, and J. Luo. Solute transport in divergent radial flow with multistep pumping. Water Resour. Res., 6:(4), 4234. [46] V. Ciriello, V. Di Federico, M. Riva, F. Cadini, J. De Sanctis, E. Zio, and A. Guadagnini. Polynomial chaos expansion for global sensitivity analysis applied to a model of radionu- clide migration in a randomly heterogeneous aquifer. Stochast. Environ. Res. Risk Assess., 49(6):;67–;76, 4235. [47] V. Ciriello, A. Guadagnini, V. Di Federico, Y. Edery, and B. Berkowitz. Comparative analysis of formulations for conservative transport in porous media through sensitivity- based parameter calibration. Water Resour. Res., 6;(;):7428–7442, 4235. [48] J. F. Clark, G. B. Hudson, M. L. Davisson, G. Woodside, and R. Herndon. Geochemical imaging of flow near an artificial recharge facility, orange county, california. Groundwater, 64(4):389–396, 4226. [49] B. E. Cole and S. E. Silliman. Capture zones for passive wells in heterogeneous unconfined aquifers. Groundwater, 57(3):;4–;:, 3;;9. [4:] V. Cvetkovic, A. Shapiro, and G. Dagan. A solute flux approach to transport in heteroge- neous formations: 4. Uncertainty analysis. Water Resour. Res., 4:(7):3599–35::, 3;;4. [4;] L. B. da Veiga, K. Lipnikov, and G. Manzini. The mimetic finite difference method for elliptic problems, volume 33. Springer, 4236. [52] G. Dagan. Solute transport in heterogeneous porous formations. J. Fluid Mech., 367:373– 399, 3;:6. 375 [53] G. Dagan. Statistical theory of groundwater flow and transport: Pore to laboratory, labo- ratory to formation, and formation to regional scale. Water Resour. Res., 44(;S), 3;:8. [54] G. Dagan, V. Cvetkovic, and A. Shapiro. A solute flux approach to transport in heteroge- neous formations: 3. the general framework. Water Resour. Res., 4:(7):358;–3598, 3;;4. [55] G. Dagan and S. P. Neuman. Subsurface flow and transport: A stochastic approach. Cam- bridge University Press, 3;;9. [56] M. Dai, J. M. Kelley, and K. O. Buesseler. Sources and migration of plutonium in ground- water at the Savannah River Site. Environ. Sci. Technol., 58(39):58;2–58;;, 4224. [57] A. Darracq, F. Greffe, F. Hannerz, G. Destouni, and V. Cvetkovic. Nutrient transport scenarios in a changing Stockholm and Malaren valley region, Sweden. Water Sci. Technol., 73(5-6):53–5:, 4227. [58] F. P. J. de Barros. Evaluating the combined effects of source zone mass release rates and aquifer heterogeneity on solute discharge uncertainty. Adv. Water Resour., 339:362–372, 423:. [59] F. P. J. de Barros, A. Bellin, V. Cvetkovic, G. Dagan, and A. Fiori. Aquifer heterogeneity controls on adverse human health effects and the concept of the hazard attenuation factor. Water Resour. Res., 74(:):7;33–7;44, 4238. [5:] F. P. J. de Barros, S. Ezzedine, and Y. Rubin. Impact of hydrogeological data on measures of uncertainty, site characterization and environmental performance metrics. Adv. Water Resour., 58:73–85, 4234. [5;] F. P. J. de Barros, D. Fernàndez-Garcia, D. Bolster, and X. Sanchez-Vila. A risk-based probabilistic framework to estimate the endpoint of remediation: Concentration rebound by rate-limited mass transfer. Water Resour. Res., 6;(6):3;4;–3;64, 4235. [62] F. P. J. de Barros and A. Fiori. First-order based cumulative distribution function for solute concentration in heterogeneous aquifers: Theoretical analysis and implications for human health risk assessment. Water Resour. Res., 72(7):623:–6259, 4236. [63] F. P. J. de Barros, A. Fiori, F. Boso, and A. Bellin. A theoretical framework for modeling dilution enhancement of non-reactive solutes in heterogeneous porous media. J. Contam. Hydrol., 397:94–:5, 4237. [64] F. P. J. de Barros, A. Guadagnini, D. Fernàndez-Garcia, M. Riva, and X. Sanchez-Vila. Controlling scaling metrics for improved characterization of well-head protection regions. J. Hydrol., 6;6:329–337, 4235. [65] F.P.J.deBarrosandW.Nowak. Onthelinkbetweencontaminantsourcereleaseconditions and plume prediction uncertainty. J. Contam. Hydrol., 338(3-6):46–56, 4232. [66] F. P. J. de Barros and Y. Rubin. A risk-driven approach for subsurface site characterization. Water Resour. Res., 66(3), 422:. [67] F. P. J. de Barros, Y. Rubin, and R. M. Maxwell. The concept of comparative information yield curves and its application to risk-based site characterization. Water Resour. Res., 67(8), 422;. [68] M. Denham and C. Eddy-Dilek. Influences on effective decay rates of radionuclides in groundwater: F-Area Seepage Basins, Savannah River Site. WM4239 Conf. Phoenix, Ari- zona, USA, 4239. 376 [69] M. Dentz and J. Carrera. Effective dispersion in temporally fluctuating flow through a heterogeneous medium. Physical Review E, 8:(5):258532, 4225. [6:] M. Dentz and F. P. J. de Barros. Dispersion variance for transport in heterogeneous porous media. Water Resour. Res., 6;(8):5665–5683, 4235. [6;] M. Dentz, T. Le Borgne, A. Englert, and B. Bijeljic. Mixing, spreading and reaction in heterogeneous media: A brief review. J. Contam. Hydrol., 342:3–39, 4233. [72] G. Destouni and A. Darracq. Nutrient cycling and N4O emissions in a changing climate: the subsurface water system role. Environ. Res. Lett., 6(25722:):3–9, 422;. [73] M. Di Dato, A. Fiori, F. P. J. de Barros, and A. Bellin. Radial solute transport in highly heterogeneous aquifers: Modeling and experimental comparison. Water Resour. Res., 4239. [74] V. Di Federico and S. P. Neuman. Scaling of random fields by means of truncated power variograms and associated spectra. Water Resour. Res., 55(7):3297–32:7, 3;;9. [75] P. M. Doyen. Permeability, conductivity, and pore geometry of sandstone. J. Geophys. Res., ;5(B9):994;–9962, 3;::. [76] D. R. Easterling, K. Kunkel, J. Arnold, T. Knutson, A. LeGrande, L. R. Leung, R. Vose, D. Waliser, and M. Wehner. Precipitation change in the United States. 4239. [77] M. D. Einarson and D. M. Mackay. Peer reviewed: predicting impacts of groundwater contamination. Environ. Sci. Technol., 57(5):88A–95A, 4223. [78] D. E. Ellis and P. W. Hadley. Sustainable remediation white paper−Integrating sustainable principles, practices, and metrics into remediation projects. Remediat. J., 3;(5):7–336, 422;. [79] R. Enzenhoefer, W. Nowak, and R. Helmig. Probabilistic exposure risk assessment with advective–dispersive well vulnerability criteria. Adv. Water Resour., 58:343–354, 4234. [7:] O. Falivene, P. Arbus, A. Gardiner, G. Pickup, J. A. Muoz, and L. Cabrera. Best practice stochasticfaciesmodelingfromachannel-fillturbiditesandstoneanalog(thequarryoutcrop, eocene ainsa basin, northeast spain). AAPG bulletin, ;2(9):3225–324;, 4228. [7;] C. Fallico. Reconsideration at field scale of the relationship between hydraulic conductivity and porosity: the case of a sandy aquifer in south italy. The Scientific World Journal, 4236, 4236. [82] J. S. Famiglietti. The global groundwater crisis. Nat. Clim. Chang., 6(33):;67–;6:, 4236. [83] R. Fay and M. Mumtaz. Development of a priority list of chemical mixtures occurring at 33:: hazardous waste sites, using the hazdat database. Food and chemical toxicology, 56(33-34):3385–3387, 3;;8. [84] B. Faybishenko, H. Wainwright, M. Denham, M. Amidon, M. Millings, G. Flach, and C. Eddy-Dilek. Basic climatic and hydrological data mining and analytics of the Savan- nah River site F-Area. 4238. [85] D. Fernàndez-Garcia, T. H. Illangasekare, and H. Rajaram. Differences in the scale- dependence of dispersivity estimated from temporal and spatial moments in chemically and physically heterogeneous porous media. Adv. Water Resour., 4:(9):967–97;, 4227. 377 [86] V. Feron, C. Hendriksen, A. Speek, H. Til, and B. Spit. Lifespan oral toxicity study of vinyl chloride in rats. Food and cosmetics toxicology, 3;:539–555, 3;:3. [87] A. D. Festger and G. R. Walter. The capture efficiency map: The capture zone under time-varying flow. Groundwater, 62(8):83;–84:, 4224. [88] L. Feyen, K. J. Beven, F. De Smedt, and J. Freer. Stochastic capture zone delineation within the generalized likelihood uncertainty estimation methodology: conditioning on head observations. Water Resour. Res., 59(5):847–85:, 4223. [89] A. Fiori. The lagrangian concentration approach for determining dilution in aquifer trans- port: Theoretical analysis and comparison with field experiments. Water Resour. Res., 59(34):5327–5336, 4223. [8:] A. Fiori, S. Berglund, V. Cvetkovic, and G. Dagan. A first-order analysis of solute flux statistics in aquifers: The combined effect of pore-scale dispersion, sampling, and linear sorption kinetics. Water Resour. Res., 5:(:), 4224. [8;] A. Fiori and G. Dagan. Concentration fluctuations in aquifer transport: A rigorous first- order solution and applications. J. Contam. Hydrol., 67(3):35;–385, 4222. [92] G. Flach. Groundwater flow model of the general separations area using Porflow (U). WSRC-TR-4226-22328.WestinghouseSavannahRiverCompanyLLC,SavannahRiverSite, Aiken, SC 4;:2:. 4226. [93] G. P. Flach, S. a. Crisman, and F. J. Molz. Comparison of single-domain and dual-domain subsurface transport models. Ground Water, 64(8-9)::37–:4:, 4226. [94] G. E. Fogg, C. D. Noyes, and S. F. Carle. Geologically based model of heterogeneous hydraulic conductivity in an alluvial setting. Hydrogeol. J., 8(3):353–365, 3;;:. [95] S. Franzetti and A. Guadagnini. Probabilistic estimation of well catchments in heteroge- neous aquifers. J. Hydrol., 396(3-4):36;–393, 3;;8. [96] D. Franzmeier. Estimation of hydraulic conductivity from effective porosity data for some indiana soils. Soil Sci. Soc. Am. J., 77(8):3:23–3:25, 3;;3. [97] R. A. Freeze and J. A. Cherry. Groundwater, 826 pp, 3;9;. [98] M.Freshley,S.S.Hubbard,G.Flach,V.Freedman,D.Agarwal,B.Andre,Y.Bott,X.Chen, J. Davis, B. Faybishenko, I. Gorton, C. Murray, D. Moulton, J. Meyer, M. Rockhold, A. Shoshani, C. Steefel, H. Wainwright, and S. Waichler. Phase II Demonstration ASCEM United States Department of Energy. 4234. [99] E. Frind, J. Molson, and D. Rudolph. Well vulnerability: a quantitative approach for source water protection. Groundwater, 66(7):954–964, 4228. [9:] J. Fu and J. J. Gómez-Hernández. Uncertainty assessment and data worth in groundwater flow and mass transport modeling using a blocking markov chain monte carlo method. J. Hydrol., 586(5):54:–563, 422;. [9;] M. Futter, R. Helliwell, M. Hutchins, and J. Aherne. Modelling the effects of changing climate and nitrogen deposition on nitrate dynamics in a scottish mountain catchment. Hydrol. Res., 62(4-5):375–388, 422;. [:2] L. W. Gelhar. Stochastic subsurface hydrology from theory to applications. Water Resour. Res., 44(;S), 3;:8. 378 [:3] L. W. Gelhar. Stochastic subsurface hydrology. Prentice-Hall, 3;;5. [:4] D. Gellens and E. Roulin. Streamflow response of belgian catchments to ipcc climate change scenarios. J. Hydrol., 432(3):464–47:, 3;;:. [:5] J. J. Gómez-Hernández and X.-H. Wen. To be or not to be multi-gaussian? a reflection on stochastic hydrogeology. Adv. Water Resour., 43(3):69–83, 3;;:. [:6] H. Gotovac, V. Cvetkovic, and R. Andricevic. Flow and travel time statistics in highly heterogeneous porous media. Water Resour. Res., 67(9), 422;. [:7] A. Grau. Material balance: Quartz cement vs. internal sources of silica, East Brae Field, offshoreUnitedKingdom: Ph.D.Thesis, Colorado Schoolof Mines, GoldenColorado. 4222. [:8] T. R. Green, M. Taniguchi, H. Kooi, J. J. Gurdak, D. M. Allen, K. M. Hiscock, H. Trei- del, and A. Aureli. Beneath the surface of global change: Impacts of climate change on groundwater. J. Hydrol., 627(5):754–782, 4233. [:9] A. Guadagnini and S. Franzetti. Time-related capture zones for contaminants in randomly heterogeneous formations. Groundwater, 59(4):475–482, 3;;;. [::] A.Guadagnini, S.Neuman, T.Nan, M.Riva, andC.Winter. Scalablestatisticsofcorrelated randomvariablesandextremesappliedtodeepboreholeporosities. Hydrol. Earth Syst. Sci., 3;(4):94;, 4237. [:;] A. Guadagnini, S. Neuman, M. Schaap, and M. Riva. Anisotropic statistical scaling of vadose zone hydraulic property estimates near maricopa, arizona. Water Resour. Res., 6;(34)::685–:69;, 4235. [;2] A. Guadagnini, S. Neuman, M. Schaap, and M. Riva. Anisotropic statistical scaling of soil and sediment texture in a stratified deep vadose zone near Maricopa, Arizona. Geoderma, 436:439–449, 4236. [;3] A. Guadagnini, M. Riva, and S. Neuman. Extended power-law scaling of heavy-tailed random air-permeability fields in fractured and sedimentary rocks. Hydrol. Earth Syst. Sci., 38(;):546;, 4234. [;4] A. W. Harbaugh. MODFLOW-4227, the US Geological Survey modular ground-water model: the ground-water flow process. US Department of the Interior, US Geological Survey Reston, 4227. [;5] T. Hashimoto, J. R. Stedinger, and D. P. Loucks. Reliability, resiliency, and vulnerability criteria for water resource system performance evaluation. Water Resour. Res., 3:(3):36–42, 3;:4. [;6] C. Haslauer, P. Guthke, A. Bárdossy, and E. Sudicky. Effects of non-gaussian copula-based hydraulic conductivity fields on macrodispersion. Water Resour. Res., 6:(9), 4234. [;7] A. E. Hassan. Water flow and solute mass flux in heterogeneous porous formations with spatially random porosity. J. Hydrol., 464(3):3–47, 4223. [;8] A. E. Hassan, J. H. Cushman, and J. W. Delleur. Significance of porosity variability to transport in heterogeneous porous media. Water Resour. Res., 56(;):446;–447;, 3;;:. [;9] A. E. Hassan and M. M. Mohamed. On using particle tracking methods to simulate trans- port in single-continuum and dual continua porous media. J. Hydrol., 497(5):464–482, 4225. 379 [;:] C. V. Henri and D. Fernàndez-Garcia. Toward efficiency in heterogeneous multispecies reactive transport modeling: A particle-tracking solution for first-order network reactions. Water Resour. Res., 72(;):9428–9452, 4236. [;;] C. V. Henri, D. Fernàndez-Garcia, and F. P. J. de Barros. Probabilistic human health risk assessment of degradation-related chemical mixtures in heterogeneous aquifers: Risk statistics, hot spots, and preferential channels. Water Resour. Res., 73(8):62:8–632:, 4237. [322] C. V. Henri, D. Fernàndez-Garcia, and F. P. J. de Barros. Assessing the joint impact of dnapl source-zone behavior and degradation products on the probabilistic characterization of human health risk. Adv. Water Resour., :::346–35:, 4238. [323] B. X. Hu, M. M. Meerschaert, W. Barrash, D. W. Hyndman, C. He, X. Li, and L. Guo. Examining the influence of heterogeneous porosity fields on conservative solute transport. J. Contam. Hydrol., 32:(5):99–::, 422;. [324] L. Hu, Y. Zhao, Y. Liu, C. Scheepens, and A. Bouchard. Updating multipoint simulations using the ensemble kalman filter. Comput. Geosci., 73:9–37, 4235. [325] P. Indelman and G. Dagan. Solute transport in divergent radial flow through heterogeneous porous media. J. Fluid Mech. Mechanics, 5:6:37;–3:4, 3;;;. [326] E. Jacobson, R. Andricevic, and J. Morrice. Probabilistic capture zone delineation based on an analytic solution. Groundwater, 62(3)::7–;7, 4224. [327] A. Journal and F. Alabert. Non-gaussian data expansion in the earth sciences. Terra Nova, 3(4):345–356, 3;:;. [328] A. G. Journel and C. V. Deutsch. Entropy and spatial disorder. Mathematical Geology, 47(5):54;–577, 3;;5. [329] V. Kapoor and P. K. Kitanidis. Concentration fluctuations and dilution in aquifers. Water Resour. Res., 56(7):33:3–33;5, 3;;:. [32:] T. H. Killian, N. L. Kolb, P. Corbo, and I. W. Marine. Environmental information doc- ument, F-Area seepage basins. Report No. DPST :7-926.E.I. du Pont de Nemours & Co, Savannah River Laboratory, Aiken SC 4;:2:, 3;:8. [32;] P. K. Kitanidis. Prediction by the method of moments of transport in a heterogeneous formation. J. Hydrol., 324(3-6):675–695, 3;::. [332] P. K. Kitanidis. Introduction to geostatistics: applications in hydrogeology. Cambridge University Press, 3;;9. [333] P. K. Kitanidis and P. L. McCarty. Delivery and mixing in the subsurface: processes and design principles for in situ remediation, volume 6. Springer Science & Business Media, 4234. [334] C. Knudby and J. Carrera. On the relationship between indicators of geostatistical, flow and transport connectivity. Adv. Water Resour., 4:(6):627–643, 4227. [335] C. E. Koltermann and S. M. Gorelick. Fractional packing model for hydraulic conductivity derived from sediment mixtures. Water Resour. Res., 53(34):54:5–54;9, 3;;7. [336] J. Kozeny. Über kapillare Leitung der Wasser in Boden. Royal Academy of Science, Vienna, Proc. Class I, 358:493–528, 3;49. 37: [337] A. Kreft and A. Zuber. On the use of the dispersion model of fluid flow. The International Journal of Applied Radiation and Isotopes, 52(33):927–92:, 3;9;. [338] M. Leeuwen, A. P. Butler, C. te Stroet, and J. A. Tompkins. Stochastic determination of well capture zones conditioned on regular grids of transmissivity measurements. Water Resour. Res., 58(6):;6;–;79, 4222. [339] M. Leeuwen, C. te Stroet, A. P. Butler, and J. A. Tompkins. Stochastic determination of well capture zones. Water Resour. Res., 56(;):4437–4445, 3;;:. [33:] S. Leray, J.-R. de Dreuzy, L. Aquilina, V. Vergnaud-Ayraud, T. Labasque, O. Bour, and T. Le Borgne. Temporal evolution of age data under transient pumping conditions. J. Hydrol., 733:777–788, 4236. [33;] P. C. Leube, F. P. J. de Barros, W. Nowak, and R. Rajagopal. Towards optimal allocation of computer resources: Trade-offs between uncertainty quantification, discretization and model reduction. Environmental modelling & software, 72:;9–329, 4235. [342] R. Li and J. W. Merchant. Modeling vulnerability of groundwater to pollution under future scenarios of climate change and biofuels-related land use change: A case study in North Dakota, USA. Sci. Total Environ., 669:54–67, 4235. [343] A. Libera, F. P. J. de Barros, B. Faybishenko, C. Eddy-Dilek, M. Denham, K. Lipnikov, D. Moulton, B. Maco, and H. Wainwright. Climate change impact on residual contaminants under sustainable remediation. J. Contam. Hydrol., 32573:, 423;. [344] A.Libera, F.P.J.deBarros, andA.Guadagnini. Influenceofpumpingoperationalschedule on solute concentrations at a well in randomly heterogeneous aquifers. J. Hydrol., 768:6;2– 724, 4239. [345] A. Libera, F. P. J. de Barros, M. Riva, and A. Guadagnini. Solute concentration at a well in non-gaussian aquifers under constant and time-varying pumping schedule. J. Contam. Hydrol., 427:59–68, 4239. [346] A. Libera, C. V. Henri, and F. P. J. de Barros. Hydraulic conductivity and porosity hetero- geneity controls on environmental performance metrics: Implications in probabilistic risk analysis. Adv. Water Resour., 423;. [347] S. Lin. Longitudinal dispersion in porous media with variable porosity. J. Hydrol., 56(3- 4):35–3;, 3;99. [348] N. Linde, T. Lochbühler, M. Dogan, and R. L. Van Dam. Tomogram-based comparison of geostatistical models: Application to the macrodispersion experiment (made) site. J. Hydrol., 753:765–778, 4237. [349] K. Lipnikov, D. Svyatskiy, and Y. Vassilevski. A monotone finite volume method for advection-diffusion equations on unstructured polygonal meshes. J. Comput. Phys., 44;(33):6239–6254, 4232. [34:] H. H. Liu and F. J. Molz. Comment on “Evidence for non-Gaussian scaling behavior in heterogeneous sedimentary formations” by Scott Painter. Water Resour. Res., 55(6):;29– ;2:, 3;;9. [34;] C. Lu, P. L. Bjerg, F. Zhang, and M. M. Broholm. Sorption of chlorinated solvents and degradation products on natural clayey tills. Chemosphere, :5(33):3689–3696, 4233. 37; [352] Z. Lu and D. Zhang. On stochastic study of well capture zones in bounded, randomly heterogeneous media. Water Resour. Res., 5;(6), 4225. [353] E. Luijendijk and T. Gleeson. How well can we predict permeability in sedimentary basins? Deriving and evaluating porosity–permeability equations for noncemented sand and clay mixtures. Geofluids, 37(3-4):89–:5, 4237. [354] D. MacEwan, M. Cayar, A. Taghavi, D. Mitchell, S. Hatchett, and R. Howitt. Hydroeco- nomic modeling of sustainable groundwater management. Water Resour. Res., 75(5):45:6– 4625, 4239. [355] G. Mariethoz, P. Renard, and J. Straubhaar. The direct sampling method to perform multiple-point geostatistical simulations. Water Resour. Res., 68(33), 4232. [356] R. Maxwell and W. Kastenberg. Stochastic environmental risk analysis: an integrated methodology for predicting cancer risk from contaminated groundwater. Stochast. Environ. Res. Risk Assess., 35(3):49–69, 3;;;. [357] R. M. Maxwell, S. F. Carle, and A. F. B. Tompson. Contamination, risk, and heterogeneity: On the effectiveness of aquifer remediation. Environ. Geol., 76(:):3993–39:8, 422:. [358] M. M. Meerschaert, M. Dogan, R. L. Dam, D. W. Hyndman, and D. A. Benson. Hydraulic conductivity fields: Gaussian or not? Water Resour. Res., 6;(:):6952–6959, 4235. [359] M. M. Meerschaert, T. J. Kozubowski, F. J. Molz, and S. Lu. Fractional laplace model for hydraulic conductivity. Geophysical Research Letters, 53(:), 4226. [35:] H. Middelkoop, K. Daamen, D. Gellens, W. Grabs, J. C. Kwadijk, H. Lang, B. W. Parmet, B. Schädler, J. Schulla, and K. Wilke. Impact of climate change on hydrological regimes and water resources management in the rhine basin. Clim. change, 6;(3):327–34:, 4223. [35;] J. Molson and E. Frind. On the use of mean groundwater age, life expectancy and cap- ture probability for defining aquifer vulnerability and time-of-travel zones for source water protection. J. Contam. Hydrol., 349(3):98–:9, 4234. [362] R. H. Morin. Negative correlation between porosity and hydraulic conductivity in sand- and-gravel aquifers at cape cod, massachusetts, usa. J. Hydrol., 538(3-6):65–74, 4228. [363] M. Moslehi, R. Rajagopal, and F. P. J. de Barros. Optimal allocation of computational resources in hydrogeological models under uncertainty. Adv. Water Resour., :5:4;;–52;, 4237. [364] Y. Mualem. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res., 34(5):735–744, 3;98. [365] P. H. Nelson. Permeability-porosity data sets for sandstones. The Leading Edge, 45(33):3365–3366, 4226. [366] P.H.Nelsonetal. Permeability-porosityrelationshipsinsedimentaryrocks. The log analyst, 57(25), 3;;6. [367] S. P. Neuman, A. Guadagnini, M. Riva, and M. Siena. Recent advances in statistical and scaling analysis of earth and environmental variables. In Advances in Hydrogeology, pages 3–47. Springer, 4235. [368] S. P. Neuman, M. Riva, and A. Guadagnini. On the geostatistical characterization of hierarchical media. Water Resour. Res., 66(4), 422:. 382 [369] R. M. Neupauer, J. D. Meiss, and D. C. Mays. Chaotic advection and reaction during engineered injection and extraction in heterogeneous porous media. Water Resour. Res., 72(4):3655–3669, 4236. [36:] S. O’Connell and D. Hou. Resilience: a new consideration for environmental remediation in an era of climate change. Remediat. J., 48(3):79–89, 4237. [36;] S. O’Connell and D. Hou. Surfactant-oxidation co-application for soil and groundwater remediation. Remediat. J., 48(3):79–89, 4237. [372] J. Odong. Evaluation of empirical formulae for determination of hydraulic conductivity based on grain-size analysis. J. Am. Sci., 5(5):76–82, 4229. [373] S. Oladyshkin, F. de Barros, and W. Nowak. Global sensitivity analysis: a flexible and efficient framework with an example from stochastic hydrogeology. Adv. Water Resour., 59:32–44, 4234. [374] J. E. Olesen, T. R. Carter, C. H. Díaz-Ambrona, S. Fronzek, T. Heidmann, T. Hickler, T. Holt, M. I. Minguez, P. Morales, J. P. Palutikof, M. Quemada, M. Ruiz-Ramos, G. H. Rubæk, F. Sau, B. Smith, and M. T. Sykes. Uncertainties in projected impacts of climate change on European agriculture and terrestrial ecosystems based on scenarios from regional climate models. Clim. Change, :3(Suppl. 3):345–365, 4229. [375] S. Painter. Evidence for non-gaussian scaling behavior in heterogeneous sedimentary for- mations. Water Resour. Res., 54(7):33:5–33;7, 3;;8. [376] M. N. Panda and L. W. Lake. Estimation of single-phase permeability from parameters of particle-size distribution. AAPG bulletin, 9:(9):324:–325;, 3;;6. [377] M. Panzeri, M. Riva, A. Guadagnini, and S. Neuman. Theory and generation of conditional, scalable sub-gaussian random fields. Water Resour. Res., 74(5):3968–3983, 4238. [378] M. J. Park, J. Y. Park, H. J. Shin, M. S. Lee, G. A. Park, I. K. Jung, and S. J. Kim. Projection of future climate change impacts on nonpoint source pollution loads for a for- est dominant dam watershed by reflecting future vegetation canopy in a soil and water assessment tool model. Water Sci. Technol., 83(:):3;97–3;:8, 4232. [379] J. C. Parker and E. Park. Modeling field-scale dense nonaqueous phase liquid dissolution kinetics in heterogeneous aquifers. Water Resour. Res., 62(7), 4226. [37:] D. Pedretti, D. Fernàndez-Garcia, D. Bolster, and X. Sanchez-Vila. On the formation of breakthrough curves tailing during convergent flow tracer tests in three-dimensional hetero- geneous aquifers. Water Resour. Res., 6;(9):6379–6395, 4235. [37;] D.Pedretti, D.Fernàndez-Garcia, X.Sanchez-Vila, D.Bolster, andD.Benson. Apparentdi- rectional mass-transfer capacity coefficients in three-dimensional anisotropic heterogeneous aquifers under radial convergent transport. Water Resour. Res., 72(4):3427–3446, 4236. [382] D. Pedretti and A. Fiori. Travel time distributions under convergent radial flow in hetero- geneous formations: Insight from the analytical solution of a stratified model. Adv. Water Resour., 82:322–32;, 4235. [383] L. Pfister, J. Kwadijk, A. Musy, A. Bronstert, and L. Hoffmann. Climate change, land use change and runoff prediction in the rhine–meuse basins. River research and applications, 42(5):44;–463, 4226. 383 [384] M. Phifer, M. Millings, and G. Flach. Hydraulic property data package for the e-area and z-area soils. Cementitious Materials and Waste Zones, Washington Savannah River Company, Savannah River Site, Aiken, South Carolina, 4228. [385] A. N. Piscopo, R. M. Neupauer, and J. R. Kasprzyk. Optimal design of active spread- ing systems to remediate sorbing groundwater contaminants in situ. J. Contam. Hydrol., 3;2:4;–65, 4238. [386] T. S. Ramanarayanan, D. E. Storm, and M. D. Smolen. Seasonal pumping variation effects on wellhead protection area delineation. Water Resour. Bull., 53(5):643–652, 3;;7. [387] P. S. C. Rao, J. W. Jawitz, C. G. Enfield, R. Falta Jr, M. D. Annable, and A. L. Wood. Technology integration for contaminated site remediation: clean-up goals and performance criteria. Groundwater quality: Natural and enhanced restoration of groundwater pollution, 497:793–79:, 4223. [388] T. E. Reilly and D. W. Pollock. Sources of water to wells for transient cyclic systems. Groundwater, 56(8):;9;–;::, 3;;8. [389] N. Remy, A. Boucher, and J. Wu. Applied geostatistics with SGeMS: a user’s guide. Cam- bridge University Press, 422;. [38:] P. Renard and D. Allard. Connectivity metrics for subsurface flow and transport. Adv. Water Resour., 73:38:–3;8, 4235. [38;] M. Riva, A. Guadagnini, and F. Ballio. Time-related capture zones for radial flow in two di- mensional randomly heterogeneous media. Stochast. Environ. Res. Risk Assess., 35(5):439– 452, 3;;;. [392] M. Riva, A. Guadagnini, and M. De Simoni. Assessment of uncertainty associated with the estimation of well catchments by moment equations. Adv. Water Resour., 4;(7):898–8;3, 4228. [393] M. Riva, A. Guadagnini, D. Fernandez-Garcia, X. Sanchez-Vila, and T. Ptak. Relative importance of geostatistical and transport models in describing heavily tailed breakthrough curves at the lauswiesen site. J. Contam. Hydrol., 323(3):3–35, 422:. [394] M. Riva, A. Guadagnini, and S. P. Neuman. Theoretical analysis of non-gaussian hetero- geneity effects on subsurface flow and transport. Water Resour. Res., 75(6):4;;:–5234, 4239. [395] M. Riva, L. Guadagnini, and A. Guadagnini. Effects of uncertainty of lithofacies, conduc- tivity and porosity distributions on stochastic interpretations of a field scale tracer test. Stochast. Environ. Res. Risk Assess., 46(9):;77–;92, 4232. [396] M. Riva, L. Guadagnini, A. Guadagnini, T. Ptak, and E. Martac. Probabilistic study of well capture zones distribution at the lauswiesen field site. J. Contam. Hydrol., ::(3):;4–33:, 4228. [397] M. Riva, S. P. Neuman, and A. Guadagnini. Sub-gaussian model of processes with heavy- tailed distributions applied to air permeabilities of fractured tuff. Stochast. Environ. Res. Risk Assess., 49(3):3;7–429, 4235. [398] M. Riva, S. P. Neuman, and A. Guadagnini. New scaling model for variables and increments with heavy-tailed distributions. Water Resour. Res., 73(8):6845–6856, 4237. 384 [399] M.Riva,S.P.Neuman,A.Guadagnini,andM.Siena. Anisotropicscalingofbereasandstone log air permeability statistics. Vadose Zone Journal, 34(5), 4235. [39:] M. Riva, M. Panzeri, A. Guadagnini, and S. P. Neuman. Simulation and analysis of scalable non-gaussian statistically anisotropic random functions. J. Hydrol., 753:::–;7, 4237. [39;] M. Riva, X. Sanchez-Vila, and A. Guadagnini. Estimation of spatial covariance of log conductivity from particle size data. Water Resour. Res., 72(8):74;:–752:, 4236. [3:2] C. B. Rizzo and F. P. J. de Barros. Minimum hydraulic resistance and least resistance path in heterogeneous porous media. Water Resour. Res., 75(32)::7;8–:835, 4239. [3:3] C. B. Rizzo, A. Nakano, and F. P. J. de Barros. Par4: Parallel random walk particle tracking method for solute transport in porous media. Computer Physics Communications, 423;. [3:4] Y. Rubin. Applied stochastic hydrogeology. Oxford University Press, 4225. [3:5] Y. Rubin and G. Dagan. Stochastic analysis of boundaries effects on head spatial variability in heterogeneous aquifers 4. impervious boundary. Water Resour. Res., 47(6):929–934, 3;:;. [3:6] Y. Rubin and G. Dagan. Conditional estimation of solute travel time in heterogeneous formations: Impact of transmissivity measurements. Water Resour. Res., 4:(6):3255–3262, 3;;4. [3:7] Y. Rubin and A. G. Journel. Simulation of non-gaussian space random functions for mod- eling transport in groundwater. Water Resour. Res., 49(9):3933–3943, 3;;3. [3:8] P. Salandin and V. Fiorotto. Solute transport in highly heterogeneous aquifers. Water Resour. Res., 56(7):;6;–;83, 3;;:. [3:9] X. Sánchez-Vila, J. Carrera, and J. P. Girardi. Scale effects in transmissivity. J. Hydrol., 3:5(3-4):3–44, 3;;8. [3::] X. Sanchez-Vila and A. Guadagnini. Travel time and trajectory moments of conservative solutes in three dimensional heterogeneous porous media under mean uniform flow. Adv. Water Resour., 4:(7):64;–65;, 4227. [3:;] D. S. Sassen, S. S. Hubbard, S. A. Bea, J. Chen, N. Spycher, and M. E. Denham. Re- active facies: An approach for parameterizing field-scale reactive transport models using geophysical methods. Water Resour. Res., 6:(32):3–42, 4234. [3;2] D. Schiedek, B. Sundelin, J. W. Readman, and R. W. Macdonald. Interactions between climate change and contaminants. Mar. Pollut. Bull., 76(34):3:67–3:78, 4229. [3;3] F. Schmidt, H. M. Wainwright, B. Faybishenko, M. Denham, and C. Eddy-Dilek. In situ monitoring of groundwater contamination using the kalman filter. Environ. Sci. Technol., 74(35):963:–9647, 423:. [3;4] C. Scholz, F. Wirner, J. Götz, U. Rüde, G. E. Schröder-Turk, K. Mecke, and C. Bechinger. Permeability of porous materials determined from the Euler characteristic. Physical Review Letters, 32;(48):486726, 4234. [3;5] J.-O. Selroos and V. Cvetkovic. Mass flux statistics of kinetically sorbing solute in heteroge- neous aquifers: Analytical solution and comparison with simulations. Water Resour. Res., 52(3):85–8;, 3;;6. 385 [3;6] M. Siena, A. Guadagnini, M. Riva, B. Bijeljic, J. P. Nunes, and M. Blunt. Statistical scaling of pore-scale lagrangian velocities in natural porous media. Physical Review E, ;2(4):245235, 4236. [3;7] M. Siena, A. Guadagnini, M. Riva, and S. Neuman. Extended power-law scaling of air permeabilities measured on a block of tuff. Hydrol. Earth Syst. Sci., 38(3):4;–64, 4234. [3;8] E. R. Siirila and R. M. Maxwell. Evaluating effective reaction rates of kinetically driven solutesinlarge-scale, statisticallyanisotropicmedia: Humanhealthriskimplications. Water Resour. Res., 6:(6), 4234. [3;9] E. R. Siirila and R. M. Maxwell. A new perspective on human health risk assessment: Development of a time dependent methodology and the effect of varying exposure durations. Sci. Total Environ., 653:443–454, 4234. [3;:] S. Silliman and A. Wright. Stochastic analysis of paths of high hydraulic conductivity in porous media. Water Resour. Res., 46(33):3;23–3;32, 3;::. [3;;] A. M. S. Sjoeng, O. Kaste, and R. F. Wright. Modelling future NO5 leaching from an upland headwater catchment in SW Norway using the MAGIC model: II. Simulation of future nitrate leaching given scenarios of climate change and nitrogen deposition. Hydrol. Res., 62(4-5):439–455, 422;. [422] C. Sprenger, N. Hartog, M. Hernández, E. Vilanova, G. Grützmacher, F. Scheibler, and S. Hannappel. Inventory of managed aquifer recharge sites in europe: historical develop- ment, current situation and perspectives. Hydrogeol. J., pages 3–36, 4239. [423] State Water Resources Control Board. Regulations Related to Recycled Water. 4237. [424] F. Stauffer, S. Attinger, S. Zimmermann, and W. Kinzelbach. Uncertainty estimation of well catchments in heterogeneous aquifers. Water Resour. Res., 5:(33), 4224. [425] S. Strebelle. Conditional simulation of complex geological structures using multiple-point statistics. Mathematical Geology, 56(3):3–43, 4224. [426] T. K. Tokunaga, J. Wan, and M. E. Denham. Estimates of vadose zone drainage from a capped seepage basin, F-Area, Savannah River Site. Vadose Zone Journal, 33(5), 4234. [427] A. F. Tompson, S. F. Carle, N. D. Rosenberg, and R. M. Maxwell. Analysis of groundwater migration from artificial recharge in a large urban aquifer: A simulation perspective. Water Resour. Res., 57(32):4;:3–4;;:, 3;;;. [428] U. S. Environmental Protection Agency. Risk Assessment Guidance for superfund, Volume 3: Humanhealthmanual(PartA),Rep.EPA/762/3-:;/224. Technicalreport,EPA/762/3- :;/224, 3;:;. [429] U. S. Environmental Protection Agency. Health effects assessment summary tables, fy3;;9 update, environ. criteria and assess. off., off. of health and environ. assess., off. of res. and dev., cincinnati, ohio. 3;;9. [42:] U. S. Environmental Protection Agency. Anaerobic biodegradation rates of organic chem- icals in groundwater: A summary of field and laboratory studies, off. of solid waste, wash- ington, d. c. 3;;;. [42;] U.S.EnvironmentalProtectionAgency. Nationalprimarydrinkingwaterregulations. EPA :38-F-2;-226, May 422;. 386 [432] U. S. Environmental Protection Agency, https://www.epa.gov/brownfields. [433] U. S. Environmental Protection Agency, https://www.epa.gov/superfund. [434] A. Van Bokhoven. The impact of climate change on the water quality of the rhine river. 4228. [435] M. Van Vliet and J. Zwolsman. Impact of summer droughts on the water quality of the meuse river. J. Hydrol., 575(3):3–39, 422:. [436] M. D. Varljen and J. Shafer. Assessment of uncertainty in time-related capture zones using conditional simulation of hydraulic conductivity. Groundwater, 4;(7):959–96:, 3;;3. [437] S. Vassolo, W. Kinzelbach, and W. Schäfer. Determination of a well head protection zone by stochastic inverse modelling. J. Hydrol., 428(5-6):48:–4:2, 3;;:. [438] V. V. Vesselinov. Uncertainties in transient capture-zone estimates of groundwater supply wells. Journal of Contemporary Water Research & Education, 359(3):3–9, 4229. [439] A. Visser, J. Kroes, M. T. van Vliet, S. Blenkinsop, H. J. Fowler, and H. P. Broers. Cli- mate change impacts on the leaching of a heavy metal contamination in a small lowland catchment. J. Contam. Hydrol., 349(3):69–86, 4234. [43:] H. Wainwright, B. Faybishenko, S. Molins, J. Davis, B. Arora, Pau G., J. Johnson, G. Flach, M. Denham, C. Eddy-dilek, D. Moulton, K. Lipnikov, C. Gable, T. Miller, and M. Freshley. Effective long-term monitoring strategies by integrating reactive transport models with in situ geochemical measurements 38434. WM4238 Conf., pages 3–37, 4238. [43;] H. Wainwright, S. Molins, J. Davis, B. Arora, B. Faybishenko, H. Krishnan, S. Hubbard, G. Flach, M. Denham, and C. Eddy-Dilek. Using ASCEM modeling and visualization to inform stakeholders of contaminant plume evolution and remediation efficacy at F-Basin Savannah River, SC. WM4237 Conf., pages 3–36, 4237. [442] H. M. Wainwright, J. Chen, D. S. Sassen, and S. S. Hubbard. Bayesian hierarchical ap- proach and geophysical data sets for estimation of reactive facies over plume scales. Water Resour. Res., 72(8):6786–67:6, 4236. [443] J. E. Warren, F. F. Skiba, et al. Macroscopic dispersion. Society of Petroleum Engineers Journal, 6(25):437–452, 3;86. [444] Water Education Foundation. The 4236 Sustainable Groundwater Management Act: A Handbook to Understanding and Implementing the Law. 4237. [445] X.-H. Wen and J. J. Gómez-Hernández. Numerical modeling of macrodispersion in hetero- geneous media: a comparison of multi-gaussian and non-multi-gaussian models. J. Contam. Hydrol., 52(3):34;–378, 3;;:. [446] P. G. Whitehead, R. L. Wilby, R. W. Battarbee, M. Kernan, Wade, and A. J. A review of the potential impacts of climate change on surface water quality. Hydrol. Sci. J., 76(3):323– 345, 422;. [447] R. Wilby, P. Whitehead, A. Wade, D. Butterfield, R. Davis, and G. Watts. Integrated modelling of climate change impacts on water resources and quality in a lowland catchment: River kennet, uk. J. Hydrol., 552(3):426–442, 4228. 387 [448] M. Willmann, J. Carrera, and X. Sánchez-Vila. Transport upscaling in heterogeneous aquifers: Whatphysicalparameterscontrolmemoryfunctions? Water Resour. Res.,66(34), 422:. [449] W. Worthy, M. D. Abkowitz, and J. H. Clarke. A systematic approach to the evaluation of rcra disposal facilities under future climate-induced events. Remediat. J., pages 93–:3, 4237. [44:] W. Worthy, J. H. Clarke, and M. D. Abkowitz. Near-surface disposal performance assess- ment: Modeling monthly precipitation and temperature in various climate environments. Remediat. J., pages ;;–32:, 4235. [44;] W. Worthy, J. H. Clarke, and M. D. Abkowitz. Surfactant-Oxidation Co-Application for soil and groundwater Remediation. Remediat. J., 48(4):323–32:, 4235. [452] T. Xu and J. J. Gómez-Hernández. Inverse sequential simulation: A new approach for the characterization of hydraulic conductivities demonstrated on a non-gaussian field. Water Resour. Res., 73(6):4449–4464, 4237. [453] J. M. Zachara, P. E. Long, J. Bargar, J. A. Davis, P. Fox, J. K. Fredrickson, M. D. Freshley, A.E.Konopka, C.Liu, J.P.McKinley, M.L.Rockhold, K.H.Williams, andS.B.Yabusaki. Persistence of uranium groundwater plumes: Contrasting mechanisms at two DOE sites in the groundwater-river interaction zone. J. Contam. Hydrol., 369:67–94, 4235. [454] A. Zarlenga, F. P. J. de Barros, and A. Fiori. Uncertainty quantification of adverse human health effects from continuously released contaminant sources in groundwater systems. J. Hydrol., 763::72–:83, 4238. [455] C. Zheng and G. D. Bennett. Applied contaminant transport modeling. Wiley Inter-Science, New York, 4224. [456] C. Zheng and P. P. Wang. MT3DMS: a modular three-dimensional multispecies transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems; documentation and user’s guide. Technical report, Alabama Univ University, 3;;;. 388 Appendix A Low-order moments of the flux-averaged concentration Thisappendixprovidestheoreticaldevelopmentssupportingtheobservedcorrespondencebetween the low and high values in the mean and variance of the flux-averaged concentration (see Figures 5.6-5.7). We start by considering the migration of a solute plume originating far from the operat- ing pumping well from a source zone of volume V 0 located upstream of the well. To simplify the derivation, we quantify the mean and variance of mass flux at a control plane located far away from the well and spanning the wellhead protection area (WHPA) along a direction normal to the mean background groundwater flow (see Figure A.3). Figure A.3: Schematic representation of the problem analyzed using the Lagrangian framework. Under this setting, we assume that the pumping rate at the well is constant such that Q w (t) = Q w and develop theoretical expressions for the moments of flux-averaged concentra- tions under a uniform-in-the-mean flow condition. As such, the key purpose of our subsequent developments is to establish an analogy with the observations stemming from the numerical re- sults illustrated in chapter 5. We obtain an approximation for the mean and variance of the 389 flux-averaged concentration at the control plane by making use of the Lagrangian framework de- veloped by [54]. For further details on the Lagrangian approach, the reader is referred to Chapters ; and 32 of [3:4]. Given the purpose of the analogy, we neglect the effects of local-scale dispersion (e.g., transport is purely advective) in the following derivations and discussion. We start by defining C(t) as the flux-averaged contaminant concentration [M/L 3 ] measured at the control plane described above (denoted by CP) C(t) = Q m (t) Q w . (A.3) We respectively denoteQ m (t) [M/T ] andQ w [L 3 /T ] as the solute mass and volumetric flow rate across the CP at time t. SinceQ w is deterministic in our work, the statistics ofC(t) depend only on the statistics of Q m (t). For initial concentration c o instantaneously and uniformly injected within a volume V 0 , the mean of Q m (t) can be expressed as: hQ m (t|˜ a)i =C 0 Z V0 g 1 (t|˜ a)d˜ a. (A.4) with˜ a indicating the initial location (withinV 0 ) of a solute particle released in the aquifer,g 1 (t|˜ a) being the solute travel time probability density function (PDF) from the source to the CP. The second moment of Q m (t) can be expressed as: h[Q m (t|˜ a)] 2 i =C 2 0 Z V0 Z V0 g 2 (t,t|˜ a 0 ,˜ a 00 )d˜ a 0 d˜ a 00 , (A.5) where g 2 (t,t|˜ a 0 ,˜ a 00 ) is the two-particle travel time PDF. The variance of Q m (t) is obtained by evaluating Var[Q m ] =h[Q m ] 2 i−hQ 2 m i. For low heterogeneity, e.g. σ 2 Y < 1, the PDF g1 can be represented by a lognormal distribution [4:,54,3;5]. This result has been verified numerically [:6,3:4]. Under this assumption, g1 scales as: g 1 (t)∼ 1 t e − ln(t) 2 (A.6) By the same token, the two-particle travel time PDF g 2 scales as: g 2 (t)∼ 1 t 2 e − ln(t) 2 (A.7) We observe that √ g 2 evolves like g 1 . Thus, the temporal evolution of the mean and variance of Q m (t) (and consequently of C(t)) are similar (see (A.4) and (A.5)). This theoretical finding, albeit under the simplified conditions here considered, constitutes an additional support to the temporal coincidence between maxima and minima of the mean and the variance ofC(t) observed in our work. Note that the travel time scaling we derive is strictly valid for low values of log- conductivity/transmissivity variance. It can nevertheless serve as an approximation for moderate to high variances as seen for example in the work of [3:8] who found a good quality agreement between longitudinal velocity covariances obtained through numerical Monte Carlo simulations and first-order theory for log-conductivity variance as large as 4 under uniform-in-the-mean flow. 38: Appendix B Supplementary material In Figures SM42-SM43 we illustrate the temporal evolution of the contaminant BTCs at the pumping well when the aquifers hydraulic transmissivity is homogeneous (set to the deterministic valueof 1m 2 /d). FigureSM42referstotheconstantpumpingscheduleofScenarioIII andFigure SM43 to the dynamic schedule of Scenario III. We notice that, when transmissivity is spatially homogeneous, the shape of the BTC observed at the operating well has some differences under steady-state or transient pumping regimes. These differences are related to the asymmetry of the concentration BTC and the peak of the mean concentration. The concentration profile becomes slightly skewed in the presence of time-varying pumping operations, as observed in Figure SM43. In Figure SM43 we observe the concentration signal at the well when the latter is not active. When pumping is activated, there is a lag period of response. This observation contributes to explain the behavior of the concentration signal at the well during inactive times. With reference to the increase of concentration, it can be observed that: (a) the activation of the extraction well attracts the plume as well as clean water which acts as a dilution mechanism; (b) when the pumping rate is low or the well is not operating, the contaminant concentration observed at the well is less averaged/diluted, thus resulting in higher concentration. This is a well-known effect, i.e., capture zones induced by the action of pumping lead to increased dilution [67,77]. Figure SM3: Mean concentrationhCi observed at the pumping well for the uniform pumping strategy of Scenario I. 38; Figure SM4: Variance of the concentrationVar[C] observed at the pumping well for the uniform pumping strategy of Scenario I. Figure SM5: Mean concentrationhCi observed at the pumping well for the transient pumping strategy of Scenario I. 392 Figure SM6: Variance of the concentrationVar[C] observed at the pumping well for the transient pumping strategy of Scenario I. Figure SM7: Mean concentrationhCi observed at the pumping well for the uniform pumping strategy of Scenario II. 393 Figure SM8: Variance of the concentrationVar[C] observed at the pumping well for the uniform pumping strategy of Scenario II. Figure SM9: Mean concentrationhCi observed at the pumping well for the transient pumping strategy of Scenario II. 394 Figure SM:: Variance of the concentrationVar[C] observed at the pumping well for the transient pumping strategy of Scenario II. Figure SM;: Coefficient of variation CV of the concentration observed at the pumping well for the uniform pumping strategy of Scenario I. 395 Figure SM32: Coefficient of variation CV of the concentration observed at the pumping well for the transient pumping strategy of Scenario I. Figure SM33: Coefficient of variation CV of the concentration observed at the pumping well for the uniform pumping strategy of Scenario II. 396 Figure SM34: Coefficient of variation CV of the concentration observed at the pumping well for the transient pumping strategy of Scenario II. FigureSM35: ProbabilityofexceedanceoftheconcentrationthresholdP (C >C ∗ )fortheuniform pumping strategy of Scenario I. 397 Figure SM36: Probability of exceedance of the concentration threshold P (C >C ∗ ) for the tran- sient pumping strategy of Scenario I. FigureSM37: ProbabilityofexceedanceoftheconcentrationthresholdP (C >C ∗ )fortheuniform pumping strategy of Scenario II. 398 Figure SM38: Probability of exceedance of the concentration threshold P (C >C ∗ ) for the tran- sient pumping strategy of Scenario II. 399 Figure SM39: Transmissivity field T 1. 39: Figure SM3:: Transmissivity field T 2. 39; Figure SM3;: Transmissivity field T 3. 3:2 Figure SM42: Mean concentrationhCi observed at the pumping well for the constant pumping strategy of Scenario III and a homogeneous aquifer. Figure SM43: Mean concentrationhCi observed at the pumping well for the transient pumping strategy of Scenario III and a homogeneous aquifer. 3:3

Abstract (if available)

Abstract Our research study focuses on better understanding the eﬀect of multiple interconnected factors on subsurface contamination. Through the use of stochastic computational modeling and fundamental ﬂow and mass transport physics we evaluate contaminant transport and the uncertainty associated with risk predictions in multiple scenarios, characterized by different heterogeneity configurations of the porous media and/or by the presence of external factors stressing the system. External factors are connected to climatic conditions (e.g, precipitation, evapotranspiration, natural recharge) and engineering stresses (e.g, groundwater pumping). The research line presented in this thesis has a direct contribution in helping risk and water managers making more informative decision on resources allocation to direct technical and economical eﬀorts towards an optimal management of water resources. Available resources need to be prioritized and assigned in a defensible manner to ensure that unacceptable identified risks are reduced to acceptable levels. The multiple scenarios explored in this thesis, interconnecting different areas that inﬂuence contaminant transport in modern day society, showed cases in which some factors could play a major role over others and therefore need proper consideration and resources allocation. Our stochastic modeling approach also highlights the contribution of different factors in increasing or reducing contaminant transport predictions’ uncertainty, whose quantification is a key element. Our research then enables better allocation of usable resources towards reducing uncertainty in decision making (e.g., in risk assessment or aquifer remediation activities) which is well aligned with goal-oriented site characterization frameworks.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

$Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure$

PDF

Hydraulic fracturing and the environment: risk assessment for groundwater contamination from well casing failure

PDF

Efficient stochastic simulations of hydrogeological systems: from model complexity to data assimilation

PDF

Hybrid physics-based and data-driven computational approaches for multi-scale hydrologic systems under heterogeneity

PDF

On the transport dynamics of miscible fluids in porous media under different sources of disorder

PDF

Thermally driven water treatment with membrane distillation: membrane performance, waste heat integration, and cooling analysis

PDF

Developing high-resolution spatiotemporal methods to model and quantify water use for energy

PDF

Developing frameworks to quantify the operational and environmental performance of energy systems within the context of climate change

PDF

Application of the adaptive river management approach to Ayamama River in Turkey

PDF

Regeneration of used petroleum-based lubricants and biolubricants by a novel green and sustainable technology

PDF

Stochastic oilfield optimization under uncertainty in future development plans

PDF

Machine-learning approaches for modeling of complex materials and media

PDF

Wastewater-based epidemiology for emerging biological contaminants

PDF

Investigating the role of urban emission sources on redox-active PM compounds and the chemical analysis of the standardized diesel exhaust particles

PDF

Flow and thermal transport at porous interfaces

PDF

Integrating systems of desalination and potable reuse: reduced energy consumption for increased water supply

$Chemical and toxicological characteristics and historical trends of size-fractioned particulate matter from traffic-related emissions in Los Angeles$

PDF

Chemical and toxicological characteristics and historical trends of size-fractioned particulate matter from traffic-related emissions in Los Angeles

PDF

Identifying and mitigating the effects of urban heat islands in California

PDF

Chemical and toxicological characteristics of particulate matter in urban environments with a focus on its sources, associated health impacts and mitigation policies

PDF

Defending industrial control systems: an end-to-end approach for managing cyber-physical risk

PDF

Efficient processing of streaming data in multi-user and multi-abstraction workflows

Asset Metadata

Creator Libera, Arianna (author)

Core Title Groundwater contaminant transport predictions in a sustainable water management scenario

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Engineering (Environmental Engineering)

Publication Date 08/15/2019

Defense Date 05/08/2019

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag groundwater contaminant transport,OAI-PMH Harvest,sustainable water management

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor de Barros, Felipe (committee chair), Moghaddam, Mahta (committee member), Sanders, Kelly (committee member)

Creator Email libera@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-211979

Unique identifier UC11663206

Identifier etd-LiberaAria-7803.pdf (filename),usctheses-c89-211979 (legacy record id)

Legacy Identifier etd-LiberaAria-7803.pdf

Dmrecord 211979

Document Type Dissertation

Format application/pdf (imt)

Rights Libera, Arianna

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA