Twenty Lessons Drawn From My Subsurface Hydrology Career

This personal perspective complements a Brief Autobiography I (Figure 1) wrote for the journal Ground Water (Neuman, 2008). Whereas there I described my personal life from early childhood to my career stage at the time, here I focus on key science lessons my academic research, teaching, and professional consulting have taught me. Though many of these lessons have been shared with students and colleagues throughout my career, never have I had the opportunity to collect and summarize them for a wider audience as I do here. I thank the Editors for the opportunity to do so, Professor Emeritus Steve Burgess of the University of Washington and Professor Alberto Guadagnini of the Politecnico di Milano, Italy, for helping me to do so better, and hope readers will find some of these lessons intriguing enough to explore and develop the related ideas further.


Lessons Learned While at UC Berkeley
As told in my Brief Autobiography , the first milestone of my professional education was a 1963 Bachelor of Science degree in Geology, with a minor in Physics, from the Hebrew University in Jerusalem, Israel. While searching for a career that would enable me to combine these two areas of science, I was made aware of a need for highly trained hydrologists and water resources specialists in Israel at the time. I decided to pursue such training in the United States and was fortunate to have been offered a generous research assistantship by the late Professor Paul Witherspoon (Figure 2) at the University of California in Berkeley. I of course accepted with enthusiasm (not the least due to the aura surrounding Berkeley faculty contributions to nuclear physics and the Manhattan Project) and began, in January 1964, studies toward an MS degree in Civil and Geological Engineering with Paul. My years at Berkeley were the fulfillment of a dream. The University ranked first in the Nation overall, and its Department of Civil Engineering was recognized to be the best. Paul started developing what would soon become one of the world's leading educational and research programs in quantitative groundwater hydrology. His group initially consisted of Al Freeze, Iraj Javandel, and me. Taking engineering and mathematics classes at Berkeley was a veritable intellectual challenge. Taking classes and working with Paul was an experience of a lifetime. Our most unusual, relaxed, and productive class was one in which the four of us explored jointly and individually the mathematics behind key chapters in Muskat's (1937) text on flow through porous media, learning with and from each other.
Having spent his early research career as an experimentalist, Paul had assigned to me as an MS project the development of a laboratory heat flow analog for drawdown in a low permeability caprock in response to water withdrawal through a well from an underlying aquifer. His purpose was to verify experimentally numerical studies he had done with one of the earliest finite difference simulators of porous media flow. Paul's simulations suggested that the ratio between drawdowns in the caprock and the aquifer would yield information about the hydraulic diffusivity of the caprock, needed to determine whether aquifers in some urban areas of the United States could serve safely and economically as reservoirs for the artificial injection and storage of natural gas.
Having watched Iraj struggle with a similar project, I was determined to avoid the drudgery. I proposed to Paul that it might be equally if not more interesting and productive to address the same problem mathematically. To his everlasting credit, Paul relented on condition that I complete my "mathematical experiment," which was bound to fail, in 2 weeks or else I would start working on the "real thing." Absent a choice, and with a good measure of luck, I came up with a surprisingly (to both of us) simple analytical solution to the problem which not only confirmed Paul's numerical results but also, more importantly, enabled us to rewrite and generalize them in terms of only a few dimensionless quantities. The solution gave rise to the Neuman-Witherspoon Ratio Method for the evaluation of aquitard and aquiclude properties (Neuman & Witherspoon, 1972). This experience taught us Lesson 1: Analytical solutions, when feasible, may have important advantages over physical analogs and numerical simulations: having a dimensionless form which renders the solution general rather than parameter specific; revealing dimensionless groups of parameters and space-time coordinates that control system behavior, which may otherwise (and in fact often do) remain unidentified; obviating the need to construct cumbersome laboratory apparatus and/or numerical grids and measure/compute results across the entire apparatus/grid at all times of interest; and generally rendering parameter estimation easier, more stable, and computationally efficient. By way of analogy to Freeman Dyson (2002), proposed to me by Steve Burgess, analytical solutions represent Cartesian Science which informs experimenters about what to measure, while experimental data (representing Baconian Science) inform analysts about parameter values in their equations.
Following this initial success, it was only natural for me to continue in this line of work by developing a much more general analytic solution for the hydrodynamic response of multilayered aquifer-aquitard systems to pumping for my PhD, published in 1968. Our excitement at having successfully compared this solution with axially symmetric finite element simulations (a novel numerical method pioneered in part by our group) was soon dwarfed by our exhilaration at having successfully validated the theory by means of 2 month-long pumping tests spanning five geologic formations near Oxnard, California (Neuman & Witherspoon, 1972). Incredibly the math worked! This, and a more recent reanalysis of the same data, using a more complete analytical description of the pumping tests (Li & Neuman, 2007), as well as my career-long experience have taught me Lesson 2: Analytical solutions may provide valid and convenient tools for hydraulic characterization of complex hydrogeologic settings by means of properly designed pumping tests. There is neither need nor justification for relying solely on numerical models for this purpose. The same does not necessarily apply to the characterization and analysis of flow in very complex systems and of reactive transport.
Lesson 3: Neither quantitative nor descriptive methods of hydrogeologic analysis can be relied on without employing both in a fully integrated manner; subsurface hydrology requires advanced knowledge of geology, physics of flow/transport through geologic media, and analytical and computational mathematics and

10.1029/2020CN000131
Perspectives of Earth and Space Scientists NEUMAN statistics, geophysics, and to an increasing degree biogeochemistry. I am encouraged to witness more and more university curricula striving toward balanced, advanced, and integrated training in many if not all of these and related areas.  (Neuman, 1972(Neuman, , 1974(Neuman, , 1975. The solution offered some important advantages over the semiempirical delayed yield model of Boulton (1963), not the least of which is its ability to deal with partially penetrating wells. The work has led to Lesson 4: Short-term, transient groundwater storage in unconfined aquifers is controlled by system compressive properties, just like in artesian or confined aquifers. Modeling flow in unconfined aquifers as if storage was due solely to pore drainage/imbibition, as commonly done, is therefore inappropriate for the analysis of water-level responses to rapidly varying stresses.

Lessons Learned During My Work in Israel
While at Bet Dagan, I was invited to take part in a research project headed by Dan Zaslavsky (future Water Commissioner of Israel) and Gedeon Dagan (future winner of the 1998 Stockholm Water Prize and AGU's 2005 Robert E. Horton Medal) at the Technion (Israel Institute of Technology) in Haifa, which allowed me to develop the first finite element simulator of transient saturated-unsaturated flow in the two-dimensional vertical plane and in three dimensions under axial symmetry, UNSAT2 (Neuman, 1973a). The code was soon applied to major projects such as the redesign of Guri Dam in Venezuela by Harza engineers and later formed the basis for a popular Microsoft Window based code HYDRUS-2D (Simunek et al., 1999) developed at the U.S. Salinity Lab in Riverside, California. The work taught me Lesson 5: Modeling seepage through earth dams, levies, embankments, and other earth structures without accounting for unsaturated flow, as was (and often still remains) the practice in civil and geological engineering, is neither valid nor expedient; the same holds true for the common hydrologic practice of modeling flow in unconfined aquifers without considering the vadose zone (more on this below).

Lessons Drawn From My University of Arizon Years
Inspired by the work of Emsellem and de Marsily (1971) at the Paris School of Mines in Fontainebleau, France, I came to realize that to estimate groundwater model parameters based on observed heads and fluxes, one must regularize the problem by imposing suitable constraints on the admissible parameter space (Neuman, 1973b). With time, the idea evolved into a statistical parameter estimation method based on the incorporation of prior geostatistical information about aquifer parameters into a maximum likelihood estimation criterion proposed, developed, and implemented by my former doctoral student (later Professor and VP Research at the Polytechnic University of Catalonia, Spain) Jesus Carrera (Carrera & Neuman, 1986a, 1986b, 1986c. The enduring success of this approach suggests the following. Lesson 6: There is a fundamental difference between the estimation of parameters in groundwater and surface water models. Whereas groundwater model parameters (such as hydraulic conductivity, transmissivity, storativity) are considered to have reasonably well-defined physical meanings, the same is generally not the case with surface model parameters. The physical nature of groundwater model parameters allows, in principle and often in practice, estimating them independently of the model by means of local methods such as slug and pumping tests. Projecting such local measurements geostatistically (more on this below) onto the model domains and treating them as prior information in the manner of Carrera and Neuman, allows one to estimate with definable reliably a much larger number of parameters than would otherwise be allowable. Hence, concepts and methods of parameter estimation developed in the groundwater modeling area may not be directly applicable to the surface water area, and vice versa.
In spring 1974 I received and gladly accepted an invitation from Paul Witherspoon to take over his classes at UC Berkeley for 1 year. In the following spring I was invited to visit the University of Arizona in Tucson by the late Sid Yakowitz (my coprincipal investigator on a joint NSF research project) where the late Gene Simpson (acting head of the Hydrology and Water Resources Department) with blessing of the late Stan Davis (incoming head of the department) promptly offered me a position as Full Professor (vacated by the untimely death of Chester Kisiel). I also received, a year later, an offer as Associate Professor with tenure from the Isotope Department of the venerable Weizmann Institute of Science in Israel which, with a heavy heart, I declined. I saw, and continue to see, this personal experience as Lesson 7: Young researchers aspiring to an academic career may benefit from following a postdoctoral stint at a good university with the development of a strong research portfolio at a reputed research laboratory prior to joining academia.
I started teaching at the University of Arizona, Tucson, in August 1975 and continued to do so with joy for 40 years, till my retirement at age 76+ years in 2015. Throughout my active professional and academic careers (I continue to serve in my U of A position as Regents' Professor Emeritus), I have tried to bring as much scientific and mathematical rigor to bear on the practice, science, and teaching of subsurface hydrology as the discipline would allow. With my students and coworkers, we continued developing analytical solutions and engineering methodologies that are widely used to evaluate the hydraulic properties of multiaquifer systems (e.g,. Li & Neuman, 2007), water table aquifers (e.g., Mishra & Neuman, 2010, 2011Tartakovsky & Neuman, 2007), and fractured rock formations (for an overview see Neuman, 2005a). Our work with Guzel Tartakovsky and Phoolendra Mishra (now Department Chair of Civil and Environmental Engineering at California State University Fullerton) has taught us Lesson 8: Flow to a well in an unconfined aquifer takes place in both the saturated zone below the water table and the unsaturated zone above. It is possible, with the aid of our analytical expressions, to evaluate hydraulic properties of both zones based on standard measurements of drawdown in the saturated zone. Potential future improvements in deep unsaturated zone monitoring should allow using such data in conjunction with saturated zone measurements to characterize unsaturated zone properties more accurately.
Extensive field and theoretical work on fractured granitic rocks (Ando et al., 2003;Neuman & Depner, 1988) and unsaturated tuffs (Illman & Neuman, 2000, 2001 led us to learn Lesson 9: Except for a few dominant fractures, flow through a fractured rock mass can often be modeled mathematically as a nonuniform, anisotropic continuum or as two or more overlapping continua. This is provided one considers explicitly the spatial variability and scale dependence of these continua using appropriate geostatistical techniques. Our group has pioneered the concept of subsurface pressure tomography and used it to image spatial variations in unsaturated fractured tuff properties by propagating controlled pneumatic pressure pulses in various directions between isolated portions of vertical and inclined boreholes, on scales of several decameters at a spatial resolution of 1 m (Vesselinov et al., 2001a(Vesselinov et al., , 2001b. Doing the same in saturated media amounts to hydraulic tomography, a concept I proposed earlier (Neuman, 1987).
Hydrogeologic parameters such as permeability and porosity have been traditionally viewed as well-defined local quantities that can be assigned unique values at each point in space. Yet subsurface fluid flow takes place in a complex geologic environment whose structural, lithologic, and petrophysical characteristics vary in ways that cannot be predicted deterministically in all relevant details. These characteristics tend to exhibit discrete and continuous variations on a multiplicity of scales, causing flow parameters to do likewise. In practice, such parameters can at best be measured at selected well locations and depth intervals, where their values depend on the scale (support volume) and mode (instrumentation and procedure) of measurement.
Estimating the parameters at points where measurements are not available entails a random error. Quite often, the measurement scale is uncertain, and the data are corrupted by experimental and interpretive errors. It became clear to me and numerous colleagues that Lesson 10: Errors and uncertainties render subsurface model parameters random and corresponding flow/transport equations stochastic.
The recognition that geology is complex, and uncertain, has prompted the development of geostatistical methods to help reconstruct it based on limited data. The most common approach is to view parameter values, determined at various points within a more-or-less distinct hydrogeologic unit, as a sample from a spatially correlated random field defined over a continuum (Dagan & Neuman, 1997). This random field is characterized by a joint (multivariate) probability density function or, equivalently, its joint ensemble (infinite sample) moments (in climate modeling, ensemble denotes a finite sample). The field fluctuates randomly from point to point in the hydrogeologic unit and from one realization to another in probability space. Its spatial statistics are obtained by sampling the field in real space across the unit, and its ensemble statistics are defined in terms of samples collected in probability space across multiple random realizations. Geostatistical analysis consists of inferring such statistics (most commonly the two leading ensemble moments, mean, and variance-covariance) from a discrete set of measurements at various locations within the hydrogeologic unit.
Once the statistical properties of relevant random parameters have been inferred from data, the next step is to solve the corresponding stochastic flow equations, conditional on the inferred parameters. The most common approach is to solve such stochastic flow equations numerically by conditional Monte Carlo simulation. This entails generating numerous equally likely random realizations of the parameter fields, solving a deterministic flow equation for each realization by standard numerical methods, and averaging the results to obtain sample moments of the solution. The approach is conceptually straightforward and has the advantage of applying to a very broad range of both linear and nonlinear flow/transport problems. It is however computationally demanding. Starting in the 1990s, my students, coworkers, and I have come to realize Lesson 11: One can derive near-exact (deterministic) integro-differential equations governing the space-time evolution of conditional (ensemble) mean heads (Neuman & Orr, 1993;Tartakovsky & Neuman, 1998a, 1998b, 1998c, solute concentrations (Morales- Casique et al., 2006a;Neuman, 1993), fluid/solute fluxes, and their second statistical moments, then solve them numerically in an approximate manner (Guadagnini & Neuman, 1999a, 1999bMorales-Casique et al., 2006b).
Our conditional moment equations, valid at all points in space-time, provide unique deterministic descriptions of flow and transport in randomly heterogeneous media, conditional on measured but uncertain parameter values. The conditional mean flow and transport equations are nonlocal in that they depend not on one but on multiple points in space-time, due to spatial and temporal correlations between parameters and state variables (heads, concentrations, and fluxes). Correspondingly, the equations are non-Darcian and non-Fickian. Hence, Lesson 12: Traditional deterministic Darcian and Fickian subsurface flow and transport equations constitute localized approximations of these nonlocal (non-Darcian and non-Fickian), stochastically based mean equations. The quality of such localized approximations (and hence, of traditional deterministic models) varies depending on conditions; see, for example, Figures 12 and 18 in Morales-Casique et al. (2006b).
Our nonlocal stochastic transport theory admits Darcy's law and analogy to Fick's first and second laws at some nominal support scale ω. This scale is not a Representative Elementary Volume (REV) in the traditional sense of Bear (1972) but a working scale at which all relevant quantities are defined and measurable. As ω is typically small relative to the domain Ω of flow and transport analyses, we refer to quantities defined on this scale as "local." Accordingly, local dispersivity (which we treat as a constant but may be defined otherwise, for example set equal to zero for studies of pure advection) is associated with a unique scale ω. In the hypothetical situation where all ω-scale medium properties and advective velocities throughout Ω are known (together with all relevant sources, initial and boundary conditions), the ω-scale transport equations are deterministic. More commonly, such quantities are known at best partially and entail measurement and/or interpretive errors. Hence, the ω-scale transport equations are stochastic.
All ω-scale quantities in our nonlocal moment equations vary in space-time in a way that is conditional on available data. In particular, our theory implies Lesson 13: Solute dispersion terms in localized (and hence traditional) transport equations vary (scale) with travel times/distances of plumes propagating through the domain Ω and the information content of the model. The more one knows about the system under investigation, the smaller do the localized dispersion terms become; in the hypothetical limit where all inputs are known deterministically, these coefficients reduce to their local (as defined previously) counterparts.
The subsurface consists of porous and fractured materials exhibiting a hierarchical geologic structure (e.g., Barton & La Pointe, 1995;Ritzi et al., 2004). The same is true of other earth environments the multiscale properties of which appear to be self-affine (monofractal) or multifractal (Molz et al., 2004;Rodriguez-Iturbe & Rinaldo, 1997;Turcotte, 1997). In particular, observed multiscale behaviors of several subsurface fluid flow and transport variables in diverse geologic settings, at various sites, and across a wide range of sampling domain scales have been found to be mutually consistent with a view of log permeability as a statistically isotropic or anisotropic self-affine random field characterized by a power variogram (or semistructure function of order 2) γ 2˜s ð Þ ∝˜s 2H where˜s is a scalar measure of distance (lag) and 0 < H ≤ 1 is a Hurst scaling coefficient Neuman & Di Federico, 2003). The power variogram scales as γ 2 r˜s ð Þ¼r 2H γ 2˜s ð Þfor any r > 0. Any random field that scales in this manner is nonstationary with stationary spatial increments, the field being Gaussian and forming fractional Brownian motion (fBm) (Mandelbrot & Van Ness, 1968). I used this idea to infer (Neuman, 1990) Lesson 14: The widely observed increase in longitudinal dispersivities (estimated based on traditional Fickian models coupled with uncalibrated flow models) with the scale of observation is consistent with a view of log permeability as a self-affine random field. Conditioning this field on measurements or through flow-model calibration reduces the rate at which dispersivities increase with observation scale.
Both phenomena are consistent with our nonlocal and localized stochastic transport theories (Neuman, 2016).
Subsurface fluid flow and solute transport take place in a multiscale heterogeneous environment. Neither these phenomena nor their host environment can be observed or described with certainty at all scales and locations of relevance. The resulting ambiguity has led to alternative conceptualizations of flow and transport and multiple ways of addressing their scale and space-time dependencies. Among these are non-Fickian transport models that postulate space-nonlocal advection-dispersion equations (ADE) containing fractional derivatives (fADE, Benson et al., 2000) or time-nonlocal equations based on a continuous time random walk (CTRW, Berkowitz et al., 2006) process. As these transport equations are not linked explicitly to underlying flow equations, these and similar postulates are aptly termed "proxy models" by Fiori et al. (2015). Work by these authors and us (Neuman & Tartakovsky, 2009) suggest Lesson 15: Proxy transport models such as fADE and CTRW have lesser predictive capabilities than models that consider jointly flow and transport physics.
Physically based transport models ubiquitously consider advective porosity (ratio between Darcy flux and advective flow velocity) to be a scalar. Yet field observations by others, theory (Neuman, 2005b), and numerical experiments (Tartakovsky & Neuman, 2008) all suggest Lesson 16: Advective porosity of geologic media may vary with advective velocity direction relative to principal directions of permeability, depending on Peclet number (measure of dispersion relative to advection).
Traditional geostatistical moment analysis allows one to infer the spatial covariance structure of hierarchical, multiscale geologic materials based on numerous measurements on a given support scale across a domain or "window" of a given length scale. The resultant sample variogram often appears to fit a stationary variogram model with constant variance (sill) and integral (spatial correlation) scale. In fact, some authors (e.g., Ritzi & Allen-King, 2007), who recognize that hierarchical sedimentary architecture and associated log hydraulic conductivity fields tend to be nonstationary, nevertheless associate them with stationary "exponential-like" transition probabilities and variograms, respectively, the latter being a consequence of the former. My coworkers and I  have learned instead Lesson 17: (a) The apparent ability of stationary spatial statistics to characterize the covariance structure of nonstationary hierarchical media is an artifact stemming from the finite size of the windows within which geologic and hydrologic variables are ubiquitously sampled, and (b) the artifact is eliminated upon characterizing the covariance structure of such media with the aid of truncated power variograms (Di Federico et al., 1999;Di Federico & Neuman, 1997), which represent stationary random fields obtained upon sampling a nonstationary fractal over finite windows.
To support our opinion, Neuman et al. (2008) noted that truncated power variograms arise formally when a hierarchical medium is sampled jointly across all geologic categories and scales within a window (Neuman, 2003a(Neuman, , 2006a; cited direct evidence that geostatistical parameters (variance and integral correlation scale) inferred on the basis of traditional variograms vary systematically with support and window scales (Gelhar, 1993;Neuman, 1994); demonstrated the ability of truncated power models to capture these variations in terms of a few scaling parameters; showed that exponential and truncated power variograms are often difficult to distinguish from each other, which helps explain why hierarchical data may appear to fit the former; noted that truncated power models are unique in their ability to represent multiscale random fields having either Gaussian or heavy-tailed symmetric Levy stable probability distributions (Chen & Hsu, 2007); detailed the way in which these models allow conditioning the spatial statistics of such fields on multiscale measurements via co-kriging (Neuman, 2006b); and illustrated these capabilities on multiscale hydraulic data from an unconfined aquifer near Tübingen, Germany.
Fractional Brownian motion (fBm) has finite, nonzero structure functions of integer orders 2q, q ≥ 1, that form power laws γ 2q˜s ð Þ ∝˜s 2qH , the exponent of which scales linearly with q. Structure functions of orders q ≥ 0 associated with a multifractal field are of the form γ q˜s ð Þ ∝˜s ξ q ð Þ where ξ(q), ξ(0) = 0, is a concave function of q (d 2 ξ/dq 2 < 0) for which there is no known universally valid expression. The literature considers self-affine (monofractal) and multifractal modes of scaling to be fundamentally different, the first arising from additive and the second from multiplicative random fields (or processes). My coworkers and I were able to demonstrate theoretically (Neuman, 2010(Neuman, , 2011 and numerically (Guadagnini & Neuman, 2011;Guadagnini et al., 2012) the following: Lesson 18: Data sampled across a finite domain from one realization of an additive random field constituting fBm give rise to positive square (or absolute) increments that behave as if the field was multifractal when in fact it is monofractal. In other words, multifractality may be (and often is) an artifact of sampling a fractal random field.
Traditional geostatistical models consider data to represent samples of multivariate Gaussian functions. Yet many earth, environmental, ecological, biological, physical, social, financial, and other variables, Y, exhibit frequency distributions that are difficult to reconcile with those of their spatial or temporal increments, ΔY.
Whereas distributions of Y (or its logarithm) are at times slightly asymmetric with relatively mild peaks and tails, those of ΔY tend to be symmetric with peaks that grow sharper, and tails that become heavier, as the separation distance (lag) between pairs of Y values decreases. No statistical model appears to capture these behaviors of Y and ΔY in a unified and consistent manner. We (Riva et al., 2015) developed Lesson 19: A new model that (a) reconciles within a unique theoretical framework the probability distributions of spatial (or temporal) variables (including physical, hydrogeological, geophysical, environmental, biological, and financial data) and the way distributions of their increments change with scale. The model is unique in (b) providing a comprehensive, self-consistent, and rigorous explanation of statistical scaling and related phenomena that have puzzled analysts for decades: power law scaling of sample structure functions (statistical moments of absolute increments) in midranges of lags, extended power law scaling (linear relations between log structure functions of successive orders) at all lags, and nonlinear (and eventually anisotropic) scaling of power law exponent with order of sample structure function.
The model has generalized sub-Gaussian form, subordinating variables to truncated fractional Brownian motion (tfBm) through the action of a lognormal subordinator (leaving open the possibility of other choices). Statistics of increment are functions of, and thus scale with, lower and upper cutoffs proportional to data resolution and sampling domain scales, respectively. These statistics can thus be scaled up or down to reflect changing resolution and/or sampling domain scales. For given cutoffs the statistics are fully defined in terms of two constant parameters and a lag-dependent correlation function that depends, in a known way, on a coefficient and a Hurst scaling exponent. Decay of this correlation function with lag was shown to be the reason why sharp peaks and heavy tails of increment probability density functions flatten and lighten with lag. Our model opens the way for conditional and unconditional simulation, and interpolation, of non-Gaussian random variables in one or more space-time dimensions. It is important to distinguish this behavior from other non-Gaussian models arising, for example, from multimodal descriptions of subsurface heterogeneity.
Our approach opens the way to modeling subsurface flow and transport stochastically in scalable, randomly heterogeneous non-Gaussian environments.
Restricted access to the subsurface and its complexity renders hydrogeologic forecasts inherently uncertain.
To the limited extent that uncertainty has been considered in groundwater practice, it has been attributed largely to errors in model parameter estimates. Yet the complex and open nature of hydrologic systems renders them prone to multiple conceptualizations and mathematical descriptions. There is now a growing recognition among hydrologists that conceptual uncertainty is often more important than errors in data and model parameters. My students and I (Lu et al., 2011;Morales-Casique et al., 2010;Neuman, 2003b;Ye et al., 2004) have learned Lesson 20: It is possible, and feasible, to quantify errors in hydrologic model outputs due to the combined effects of data, parameter, conceptual and scenarios uncertainties using methods such as maximum likelihood Bayesian model averaging. The same method can also be used to assess the worth of data and optimize data collection under such uncertainties Neuman et al., 2012;Xue et al., 2014).

Postscript
Some of the above lessons have found their way into subsurface hydrology textbooks and professional practice, while other have yet to be adopted or rejected by colleagues and future generations of subsurface hydrologists. Our lessons concerning multiscale analysis and modeling should, in my view, be of interest to scientists and engineers from many disciplines. I invite interested readers, and so likely would my coworkers, to discuss these and related subjects further with us in writing or in person. Credit for work and lessons described in this brief paper is due collectively to all my many students and collaborators from the early 1960s to this day. Though I have been unable to cite or mention all their numerous contributions, I do thank each of them for their hard work, dedication to scientific and professional excellence, and lasting friendship. I regret that not all of them were able to join the group of U of A students and coworkers pictured in the following photo (Figure 3), taken at the 2016 AGU Annual Meeting in San Francisco.