Rick Hooper Consortium of Universities for the Advancement

Rick Hooper Consortium of Universities for the Advancement of Hydrologic Science, Inc. (CUAHSI) Medford, MA, USA Querying a Landscape: Understanding the Hydrology at Panola Mountain 1

Modeling the Landscape Real World “Digital Environment” Water Sequences DNA quantity Meteorology Geomorphologist Remote sensing Hydrologist Aquatic Ecologist Vegetation Survey Biogeochemist and quality Conceptual Frameworks Valley Snowmelt Glaciated Processes? Groundwater Physical World Model -Mathematical Formulae Geographically Mapping Contribution? DOCPerifluvial Quality? Representations -Solution Oligotrophic? Backwater habitat Techniques Referenced Carbon source? Q, Redox Zones? Hyporheic exchange? Gradient, Roughness? Data • Theory/Process Knowledge Substrate Size, Stability? Hypothesis Testing Representation Thalweg? Well • sorted? Mineralogy? of. Chemistry? Perceptions this Benthic Community place • Intuition Measurements

Multi-disciplinary Description Rocks Lithology/Mineralogy Orogeny/Geomorphology Stratigraphy Soils Chemical-Physical Properties Vegetation Climate Land-use History/Human Modifications August, 2012 Aberdeen Short Course 3

Overview of a Southern Piedmont Landscape 4

Panola Mountain Research Watershed June, 2009 Isotope Hydrology and Biogeochemistry Workshop, Oregon State University 5

Panola Mountain Research Watershed June, 2009 Isotope Hydrology and Biogeochemistry Workshop, Oregon State University 6

Panola Mountain Research Watershed

Orientation 41 ha, 51 m relief, 3 ha bare rock 16°C (annual average); 1200 mm precip; 30% runoff (range 16% -45%) Granodiorite (albite) implaced within plagioclase -gneiss (amphibolite/biotite) Entisols weak Bt development on slopes Southern hardwood (oak/hickory/tulip poplar) with stands of loblolly pine Gullied landscape from past ag practices 1930’s release from ag; some pasture until 1960’s October, 2001 EMMA Shortcourse 8

Quantifying the Landscape Watershed Models and Hypothesis Testing 9

The Need for Simulation Models Quantitative conceptual model needed for advancing understanding of catchment hydrologic and biogeochemical processes Model attributes Objective measure of applicability Link scale of observations with scale of predictions

Problem Observational studies Controlled experiments impossible to perform Strict hypothesis testing impossible Spatial Scale Internal measurements ~ decimeters Catchment processes ~ hectares

Simulation Models and Catchment Studies Catchment’s well-defined boundary conditions ideal for simulation modeling Conservation of mass organizing principle Implicit (or explicit) hypothesis test Input Model Output If model predicts outputs “well, ” does that prove the model is correct?

Simulation Models as Hypothesis. Testing Tools Null hypothesis: Model is correct Inverts standard experimental design Matching observations necessary, but not sufficient, test “Power” of test to reject incorrect model critical

Decision Table TRUTH DECISION Not Null Reject Null α Type 1 Error (1 -α) Accept Null Type 2 Error (1 -β) β α: “Confidence Level” β: “Power of the Test” August, 2012 Aberdeen Short Course 14

Information Content of Data: Analysis of Birkenes Model Calibrated using rainfall/runoff and O 18 as conservative tracer Data information content determined rate, and dimensional parameters; could not determine routing

Implications Information content insufficient to use as hypothesis-testing tool Need an “inductive leap” to define new conceptual model Primary objective: hydrologic flowpaths through catchment

Proposition: Mixing Model of Soil/GW Solutions Streamwater chemical variation arises from mixing different proportions of “end members” End members have constant concentration Tracers mix conservatively

Overview of EMMA Concepts and Sample Applications 18

Fundamental Concept Explore geometry of “data cloud” Infer streamflow generation mechanisms (hydrograph separation) Infer mechanisms controlling chemistry Compare sites with one another Form data cloud by plotting solute vs. solute… 19

Time Series 20

Solute Mixing Space 21

Mixing Models “Conservative” mixing No chemical reactions Example: 30% solution A @ 5 mg/l and 70% solution B @ 10 mg/l: ▪ 0. 3*5 + 0. 7*10 = 1. 5 + 7 = 8. 5 mg/l for mixture Linear process=> principal components (eigenvalue) analysis can be used Inverse problem: observe mixture (streamwater) and evaluate sources 22

Hydrograph Separation Tracers—one application of mixing models Kinds Graphical (Quickflow/slowflow) Tracer-based ▪ Temporal (“event”/’pre-event”) ▪ Geographic (multiple source areas) 23

Graphical Separation Quick Flow, Assumed to be precipitation Slow Flow— Assumed to be Groundwater/interflow From Linsley, R. K. , M. A. Kohler, and J. L. H. Paulhus, 1975. Hydrology for Engineers, 2 nd Ed. , Mc. Graw-Hill, New York, 482 pp. 24

Tracer Separation (Temporal) In all storms, “old” water dominates hydrograph. “Old” water often assumed to be groundwater. From Hooper, R. P. and C. A. Shoemaker, 1985. A comparison of chemical and isotopic hydrograph separation. Water Resour. Res. 22: 1444 -1454. 25

Contradiction Both precipitation and groundwater cannot dominate. Resolution: incorrect source inference: quick ≠ precipitation; Old≠groundwater. Reject Hortonian overland flow (but doesn’t say what mechanism is operating. 26

Process Inference How precipitation becomes streamflow Role of groundwater and vadose water in streamflow Relative contribution of different areas of watershed Residence time of water within watershed 27

Why Tracers? Pros Integrated watershed response Easy (? ) to measure Cons No physics Incorrectly assume conservative behavior State of the art: combine tracers and hydrometric measurements 28

Assumptions for Tracers Conservative Sources significantly different concentrations Unmeasured sources have same concentration or don’t contribute significantly Sources maintain constant concentration 29

End-member Mixing Analysis Hydrograph separation using multiple tracers simultaneously N tracers define N+1 sources Use more tracers than necessary (overdefined system) to test consistency of tracers Typically use solutes as tracers 30

EMMA Terminology Source solutions called “end members” Concentration more extreme than stream Bound stream mixture Mathematical approach called “EMMA” Typical source solutions have been organic soil horizon solution, hillslope groundwater, valley bottom ground-water, throughfall, precipitation 31

Extension of EMMA Use same mathematics to compare sites to “reference” site Site can be categorized “Similar” sites could have same end members, or, more broadly, same reactions controlling “Dissimilar” sites must have chemical reactions occur to make them “similar” to reference site 32

Inverse Problem Measure Stream chemistry Groundwater chemisty Soilwater (vadose zone) chemistry Hypothesis: Streamwater can be explained as a (conservative) mixture of soilwater and groundwater end members 33

34

Mixing Diagram: Panola 35

Mixing Diagram 2: Panola 36

Problems Variability of Soil Solutions Temporal and Spatial Test: Variability within horizon over time << Variability across horizons << Streamwater variability Discretization into “boxes” Translation of end members to stream chemically unchanged 37

Soil Solution Caveats No intensive sampling of soil solutions How much sampling is adequate to assess variability? Sampling procedure itself may change chemistry 38

Results: Original and Projected End Members 39

Annual Hydrograph Separation 40

Event Hydrograph Separation 41

Inferring Watershed Structure End members define functional units of watershed Groundwater—summer baseflow, deep GW from floodplain Hillslope—winter baseflow, from uplands of catchment Organic—soil water contributing to stream only during storms

EMMA vs. Watershed Simulation Models Watershed compartments defined by data; not arbitrary More effective in organizing and querying data Lead to testable hypotheses Example: Repeated applications of EMMA to Panola Mt. Research Watershed. 43

From Mixing to Predictive Models Mixing models No chemical mechanisms Identify controlling soil environments Determine routing through these environments Long-term acidification models Time/space dilemma Combine MAGIC with EMMA 44

MAGIC Accept chemical mechanisms as true Calibrate model to soil end members (not stream as previously done) Assume groundwater is unchanging; model two upper boxes Determine routing from mixing model 45

Annual “Average” Routing Organic 0. 16 0. 84 Hillslope 0. 38 0. 36 Groundwater 0. 46 Determine routing among boxes by determining mixing proportions that yield average annual chemistry 46

Streamwater Chemistry 47

Soilwater Sulfate Chemistry 48

Soil Base Saturation 49

Event Chemistry for 2 Storms Winter frontal storm Wet antecedent conditions 62. 2 mm of rain in 24 hr; max. 6. 3 mm/hr Water yield 22% Summer convective storm Dry antecedent conditions 25. 7 mm; max 22. 4 mm in 15 minutes Water yield 3% 50

Summer Thunderstorm 51

Winter Frontal Storm 52

Results from Acidification Model Watershed structure (thin organic; low alkalinity hillslope) critical to evolution In most watershed models, structure is most arbitrary part True? Validity of end members/ mixing model Validity of chemical model 53

Testable Hypotheses Surface horizon near sulfate saturation Cation stress after sulfate saturation Rapid acidification of hillslope environment (further characterization) End member validity Spatial replication Temporally intensive sampling 54

Testing End Members Sulfate concentration declined from ~200 µeq/L to ~100 µeq/L between 1988 and 1991 55

Ca/SO 4 Mixing Diagram 56

Additional Field Data 2 tension lysimeter nests on hillslope, 199294 1 tension lysimeter nest, riparian zone, 199294 SE Tributary Stream, event sampling, 199295 Trench face samples (zero tension), 1996 57

58

Hypotheses 15 -cm tension lysimeters ≈ A-horizon after WY 90 Test spatial extent and assumes zero- tension≈tension samples Trench-face water ≈ Hillslope EM Zero-tension water moving downslope during events that contributes to streamwater 59

Hillslope/Riparian Mixing • Hillslope waters very different from “hillslope” EM • Deep hillslope tension VW≈Trenchface≠”Hillslope” EM • Riparian tension VW≈ “Hillslope” EM • Low ionic strength not artifact of outcrop 60

Huh? Change in organic EM reflected in stream But observed hillslope waters very different from “hillslope” EM Resolution: Streamwater variation reflects varying proportions of riparian zone water; hillslope contribution not evident during events 61

Implications Hydrologic routing not discernable from chemistry at watershed outlet MAGIC results not reliable… But…. what about all those SWAT models? Most of watershed not reflected in stream chemistry 62

Next Steps Focus on hillslope/riparian interface Can we see chemically distinct water entering riparian zone during storms? What is volume of riparian aquifer Contrast catchments with different alluvial aquifers Panola (dry); Sleepers River, VT (moderate), Mai, NZ (wet) 63

Using EMMA to Compare Sites Use reference site to define solute subspace Project test sites into solute subspace to determine “similarity” Extension of solute ratios to multiple dimensions Citation: Hooper, R. P. (2003), Diagnostic tools for mixing models of stream water chemistry, Water Resour. Res. , 39(3), 1055, doi: 10. 1029/2002 WR 001528. June, 2009 Isotope Hydrology and Biogeochemistry Workshop, Oregon State University 64

Measuring “Fit” to Subspace Bias Calculate residuals ei=(x’-x) Plot residuals vs. x Should have no pattern Relative bias: Σei/Average(xi) Root Mean Square Error (RMSE) Square Root (Average(Σei 2)) Relative RMSE = RMSE/Average(xi) October, 2001 EMMA Shortcourse 65

Hillslope/Stream Comparison Panola: Trenchface samples vs. “Upper Gage” (10 -ha basin) Maimai: Trenchface hollow, trenchface planar vs. K-weir (17 -ha basin) Sleepers River: Hollow well samples, planar hillslope well samples vs. W 9 B (13 -ha basin) October, 2001 EMMA Shortcourse 66

Panola Mountain 67

Maimai 68

Sleepers River Catchment W 9 B 69

Hypothesis Review Panola least similar; Maimai most similar Previous work at Panola used 41 -ha “Lower Gage”, not 10 -ha “Upper Gage” 70

Results: Rank of Panola UG Eigenvalues analysis +13% +2% +1% 71

Panola UG 2 -D Fit Lack-of-fit for both Mg and SO 4 72

Panola UG 3 -D Fit No evidence of lack-of-fit Panola UG has 3 -dimensional mixing space 73

Results: Panola Scalar Measures Statistics for projecting trench samples into 3 d UG mixing space 74

Residual Plots: Panola Si (and Na) show lack of fit, as expected, but Ca also shows lack of fit. 75

Results: Maimai Scalar Measures Hollow trench segment“fits”; planar does not, especially for Na, Ca, Si. 76

Results: Maimai Residuals 77

Results: Sleepers Scalar Measures 78

Results: Sleepers Residual Plots Bias, but no lack of fit Both bias and lack of fit also for Si 79

Multi-Site Comparison Panola hillslope did not fit stream chemistry (as expected) Hollow trench segments at Maimai and hollow wells at Sleepers did fit stream chemistry. Planar hillslope sections at Maimai and Sleepers did not fit stream. =>Hillslope form appears to be more important than wetness of catchment 80

Catchment Functional Units Riparian/Alluvial Aquifers Consistently connected to stream Dominated by inorganic rock/water reactions Upland Hollows First areas to wet up In driest areas, rarely connect to stream In wetter areas, greater connectivity, yet not consistently Upland Planar/Concave Rarely connected to stream Two Water Worlds 81

Advantages of EMMA Provides a useful framework for analyzing watershed chemical data sets Links internal observations with watershedscale observations Generates testable hypotheses that focus future field efforts 82

EMMA in Detail How do I do it? June, 2009 Isotope Hydrology and Biogeochemistry Workshop, Oregon State University 83

Steps in EMMA Bivariate solute-solute plots (mixing diagrams) to screen data, check outliers PCA of stream samples; determine rank Project end members, other sites into reference site subspace (“U-space”) Assess “fit” Mixing model: perform hydrograph separation 84

Quick Review Linear Algebra Matrix notation Matrix Operations Principal Components Analysis Standard multivariate statistical technique Geometric definition 85

Matrix Definition Matrices are n-dimensional arrays of numbers. Generally consider 2 -D arrays with dimensions of rows columns. A is a 2 2 matrix; B is a 3 2 matrix. 86

Matrix Algebra-Transpose If CT=C, matrices are symmetric 87

Matrix Multiplication Multiply rows by columns and sum the product of the elements Matrices must be conformable: number of columns of first = number of rows of second Examples (3 2)*(2 3) = (3 3) (2 3)*(3 2) = (2 2) (3 2)*(3 2) impossible 88

Matrix Multiplication Example BTB—multiply (2 3) by (3 2) = (2 2) matrix; matrices are conformable 89

Questions If B is a (3 x 2) matrix, what are the dimensions of BBT? In general, does AB = BA? Let C = BTB. Is C symmetric? 90

More Definitions A square matrix has equal numbers of rows and columns The identity matrix (I) is a square matrix with 1’s on the main diagonal and 0’s elsewhere. Iy = y. 91

Matrix Inverse The inverse of a matrix is that matrix, when multiplied by the original matrix yields the identity matrix: A-1 A = AA-1= I 92

Why would anyone do this? Systems of simultaneous equations: 3 x + 2 y + z = 21 2 x + y + 3 z = 29 5 x +2 y – 3 z = -3 93

Principal Components Analysis First step of factor analysis Mathematically, known as “eigenvalue analysis” Geometrically, translation and rotation of axes. Defines a new coordinate system for data. 94

A picture might help… Y Y’ X’ X Ellipse in original system: ax 2+bxy+cy 2 = 1 Ellipse in new system: a*x’ 2 + c*y’ 2 = 1 95

Geometry of PCA Translation: define new origin of coordinate system Rotation: first component along axis of major variation and rest orthogonal Variation: function of squared distance from origin PCA generally factors covariance or correlation matrix 96

Covariance Matrix Data are centered on origin Original units are maintained Influence of each variable proportional to its variance 97

Correlation Matrix Data are centered on mean and scaled by variance Each variable has equal variance and, so, equal weight in analysis Units are lost in analysis 98

Output of PCA Eigenvalues Proportion of variance explained by each new variable Helps to determine approximate “rank” (dimensionality) of data set Eigenvectors Define new coordinate system “Scores” are data expressed in new coordinate system 99

Dimensionality [Rank] of Data Sets Data set 1 Eigenvalues = (60%, 40%) Approximate rank=2 Data set 2 Eigenvalues= (99. 7%, 0. 3%) Approximate rank=1 Rank is indicated by the number of non-zero eigenvalues. 100

In 3 -Dimensions Eigenvalues = (99. 6%, 0. 25%, 0. 15%) 101

But in 5 -dimensions Eigenvalues=(59. 9%, 39. 9%, 0. 1%, 0. 006%, 0. 001%) 102

Orthogonal Projections Use eigenvectors to determine new coordinates (“scores”) Retain only “significant” eigenvectors Project data into lower dimensional subspace 103

Coordinate Systems 104

Problem Set: Geometry of Mixing Implications of mixing model Develop intuition of PCA Properties of data scaling Google Sheet Link: https: //docs. google. com/spreadsheets/d/1 i. CJn Ykp. M 5 v 9 IVVp. A 0 onl. Xpjv. Ib. VQWav. Ru 3 Dueg. XW bx 0/edit? usp=sharing 105

Part A. Determining pold using two different tracers Use δD and δ 18 O to determine proportion of old water separately Compare results Evaluate results 106

Temporal Hydrograph Separation Solve two simultaneous mass-balance equations for Qold and Qnew Qstream = Qold + Qnew Cstream. Qstream = Cold. Qold+ Cnew. Qnew To yield the proportion of old water 107

Answers to Part A % old water from Delta-O 18 0 100 %Old Water from Delta-D 0 108

Problem Set Part B-1 Assume the following 3 end members control stream chemistry. Where in Si-d. O-18 plane must stream samples lie? How many dimensions do data “fill”? How many non-zero eigenvalues? 109

Answer, Part B Stream samples must lie within triangle defined by end members Data fill 2 dimensions; there are 2 non-zero eigenvalues 110

Problem Set Part B-2 For the above end members, where must stream samples lie? How many dimensions do these data fill? How many non-zero eigenvalues are there? 111

Answers, B-2 Stream samples must lie on line. End members are collinear Data fill 1 dimension; only 1 non-zero eigevalue. 112

Problem Set, Part C-1 Perform a “graphical” PCA on these data • Where is the origin of the new coordinate system? • What is orientation of axes? • How much variance is explained by PC 1? 113

Answers Part C-1 114

Problem Set, Part C-2 Perform same analysis on this diagram. • How are relative magnitudes of eigenvalues different between 2 a and 2 b? 115

Answers Part C-2 Second component explains relatively more variance in Part 2 than in Part 3. 116

Problem Set, Part C-3 Perform same analysis. Compare with Part 1. How does standardization influence analysis? 117

Answers, Part C-3 Relative locations of points remain the same. Data no mostly lie between -2 and +2 for both solutes 2 nd component explains more variance than in Part 1. 118

References 119

EMMA Notation Measure p solutes; define mixing space (Sspace) to be p-dimensional X, matrix of reference site samples, (n observations p solutes); row xi (1 p) V, matrix of k retained eigenvectors (k p); define U-space is k-dimensional. 120

Orthogonal Projection The orthogonal projection, expressed in new coordinates (u), is just U=VXT The orthogonal projection, expressed in the original solutes, is X*=VT(VVT)-1 VXT What are the dimensions of U and X*? 121

Coordinate Systems 122

Assessing “Fit” of Subspace Bias Calculate residuals ei=(x’-x) Plot residuals vs. x Should have no pattern Relative bias: Σei/Average(xi) Root Mean Square Error (RMSE) Square Root (Average(Σei 2)) Relative RMSE = RMSE/Average(xi) 123

2 -D Example PC 1 explains 91% of variance; PC 2 explains 9%. 124

Project Data into 1 -D RB={0, 0}; RRMSE={8%, 4%} 125

Fit Site A to Reference Site RB={-1%, -0. 7%} RRMSE={3%, 2%} What’s “good”? Contrast values with reference site. 126

Site A Residuals Solute 1 Solute 2 Random scatter = “good” fit 127

Fit Site B to Reference Site RB={-42%, 41%}, RRMSE={45%, 43%} Note that projection doesn’t appear to be orthogonal because correlation matrix used for principal components. 128

Site B Residuals Solute 1 Solute 2 Residuals have definite structure and, therefore, exhibit lack of fit. 129

Step 1: Mixing Diagrams Generate solute-solute plots of streamwater for all combinations Plot potential end members on same graph Evaluate possible mixing models 130

Example: Stream K 1 L 131

Stream + Springs 132

Matrix Plots 133

Mixing Diagram Diagnostics Check for outliers: PCA sensitive to extreme values Check for curvature: Indicates chemical reactions (non-conservative) Assess variability of end members 134

Computer Tips Excel—Make plots manually Stats packages—allow categorical variables to code sites Use “vary symbols” to quickly show multiple sites on one plot Matrix plots are very useful R, Mat. Lab allow custom programs 135

Step 2. Assess Rank using PCA: Which matrix to factor? Covariance Matrix Preserves units of concentration Points are centered on mean Solutes influence proportional to variance Correlation Matrix Data are scaled to equal variance Unitless (e. g. , 6 µeq/L Mg= 56 µeq/L Na) Each solute tracer has equal weight 136

Standardization Issues Choice of matrix depends on interpretation Be careful about concentration range vs. analytical precision If solute variance more than 10 X others, it will dominate analysis Generally use correlation matrix 137

PCA of Panola Lower Gage First two components dominate. Try rank of 1, 2, and 3. 138

Test Fit Relative Bias = 0 for all dimensions Relative Root Mean Square Error Big decline in SO 4 from 1 - to 2 -D. Ca and SO 4 improve from 2 -D to 3 -D. 139

Sulfate 1 -D 2 -D 3 -D 140

Calcium 1 -D 2 -D 3 -D 141

Rank of Panola Lower Gage Decision is subjective Must be at least 2 -d; some further improvement for 3 -d. PC 1+PC 2 explain 94% of variation. Maximum a 3 -component mixing model can explain. 142

Results: Rank of Panola UG Eigenvalues analysis +13% +2% +1% 143

Panola UG 2 -D Fit Lack-of-fit for both Mg and SO 4 144

Panola UG 3 -D Fit No evidence of lack-of-fit Panola UG has 3 -dimensional mixing space 145

PCA Computer Tips Use any stats package that does PCA to extract eigenvectors R: Can extract eigenvectors. Packages? Be careful how eigenvectors are arranged (solute by factor or factor by solute) 146

R Commands >mycor<-cor(mydata) ‘or var() >coreig<-eigen(mycor) ‘calculates both eigenvalues and eigenvectors >sink(“myoutput. txt”) ‘re-directs output >coreig$vectors ‘types vectors to file >sink(“”) ‘output back to terminal 147

Step 3: Project into reference site mixing space Standardize end members using mean and standard deviation of reference site. Project end members into subspace defined by stream mixture Express in either standardized space (“Uspace”) or in original solute coordinates (X*-space) 148

Computer Tips for Projection R, SAS, Mathematica, Math. CAD and MS Excel can do matrix operations Here, we’ll use MS-Excel MMULT: matrix multiplication MINV: matrix inversion TRANSPOSE: matrix transpose Tricky to enter matrix formulas (SHFT-CTRLENTER) 149

Step 4. Assess fit For end members Determine spanning end members Assess difference between original and projected end members Residual analysis 150

Results: Panola Scalar Measures Statistics for projecting trench samples into 3 d UG mixing space 151

Residual Plots: Panola Si (and Na) show lack of fit, as expected, but Ca also shows lack of fit. 152

Step 5. Perform hydrograph separation Solve equations to yield mixing proportions, g Incorporate constraint ∑ gi = 1 Screen for negative proportions Multiply proportions by EM concentrations to yield predicted conc. Compare observed and predicted chem. Evaluate hydrograph separation 153

Hydro. Desktop Free, Windows-only package that accesses CUAHSI HIS Download at http: //www. hydrodesktop. org More info at http: //wdc. cuahsi. org Includes “Hydro. R” an R plug-in August, 2012 Aberdeen Short Course 154

Excel Spreadsheet Example EMMA analysis for Panola Lower Gage in emma. xls Contains steps to perform EMMA (including hydrograph separation), but not determining eigenvectors. Need a stats package for eigenvector extraction. 155