DEVELOPMENT AND APPLICATION OF A COUPLED HYDROGEOPHYSICAL INVERSION MODEL FOR ESTIMATING SOIL AND ROOT PROPERTIES By Alexandria S. Kuhl A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Environmental Geosciences Doctor of Philosophy 20 20 ABSTRACT DEVELOPMENT AND APPLICATION OF A COUPLED HYDROGEOPHYSICAL INVERSION MODEL FOR ESTIMATING SOIL AND ROOT PROPERTIES By Alexandria S. Kuhl Vegetation and water on the landscape are directly linked , with leaf area controlling the partition ing of evapotranspiration , a process which c reate s a microclimate in the atmosphere above the plants. Therefore, widespread alterations to land use has potentially larg e implications for the water balance at regional and global scales. Unfortunately, many challenges persist that limit our ability to model with high confidence the biophysical constraints on evapotranspiration. One of the largest unknowns in this complex s ystem is the root distribution , which is highly heterogeneous and dependent on both internal factors like the species and age of the plant , as well as external factors such as the climate and soil conditions. ive approaches to improve our understanding of the interplay between roots, water, and the soil. Geop hysical tools, such as electrical resistivity , have been employed for decades to study the properties of the Earth at scale s from tens to hundreds of meters. Advancements in observing features such as oil and gas reserves, aquifer properties, and contamina n t plumes has led the way to more recent work monitor ing shallow soil properties such as water content and salinity. In this dissertation , I advance the field of hydrogeophysics with the development of a novel modeling approach that utilizes electrical re sistivity data directly to parameterize root properties of site - specific hydr o logical model s . Building on prior research using coupled hydrogeophysical inversion methods to estimate soil hydrological properties, the model presented and applied here address es the pressing need for new tools to study root properties at the field scale. Chapter 1 provides a high - level introduction of the background and motivation for this area of research. Chapter 2 establishes the feasibility of the proposed modeling framewor k at a biofuel research site . Using site - specific soil and climate forcing data, I generated synthetic hydrological and electrical resistivity datasets using fixed soil and root parameters for a plot of maize . I then tested how well the model estimated those parameters under increasing levels of uncertainty. I found that even in the most data - poor scenario, the coupled hydrogeophysical inversion estimated the synthetic parameters with a high degree of accuracy . Chapter 3 proceeds to use the now - established model approach to estimate the root properties of two contrasting biofuel treatments , an annual rotation and a perennial grass . W e again found the model returned reasonable estimates of the root distribution and evapotranspiration estimates for both crop types. In Chapter 4 I take advantage of the unique ability afforded by this modeling approach to test whether a theoretical coarse root fraction crossing the plane of the electrode array coul d produce an amplified resistivity measurement akin to what has been observed in the field. Given those estimates were within reason, subsequent estimates of coarse root fraction in a forested ecosystem were then validated against an index of above ground biomass. A statistically significant relationship was found, providing evidence in the absence of excavated root data that resistivity data can be used to non - invasively estimate the extent and relative quantity of coarse roots . Chapter 5 concludes this wo rk by exploring the statistical relationship between above ground vegetation indices and the spatial and temporal heterogeneity in the observed resistivity data , providing the groundwork for future work modeling coarse root mass in a wide array of forest e cosystem s . iv This body of work is dedicated to the memory of a fellow Brockport Golden Eagle and the countless other women like her whose lives have been stolen by the violence of th ose who seek to control them; to the hundreds of brave survivors of abuse at Michigan State by a man they should have been able to trust ; and to the generations of women and other marginalized peoples whose potential has always far exceeded the limits of a society entrenched in sexism and racism. May this new decade be the dawn of a n era of equality and freedom . v ACKNOWLEDGMENTS A s of the time of submi ssion of this dissertation, I will have spent nearly a decade as a member of the Hydrogeology Lab at Michigan State. Many impactful events, both large and small , have transpired over the course of my time here and I find myself overwhelmed with gratitude to have h ad the opportunity to go on this particular journey in this extraordinary place with the most incredible people by my side. I would like to first and foremost thank my advisor, Dr. David Hyndman , for using his limitless capacity for positivity and kindness to build up th is community that I now cherish so deeply . Dave, y our generosity of spirit , unshakeable faith in your team, and in me, are what made this wonderful experience possible and turned this dissertation into a body of work that I am finally proud to call my own . Second, I would like to extend my deepest thanks to my other committee members , Drs. Anthony Kendall, Remke Van Dam, and Bruno Basso , for their unwavering support and enthusiasm even as the days grew longer . In particular , I would l ike to thank Anthony for being both a brilliant mentor and friend to me . Anthony, y ou have been a kindred spirit, delighting in every technical or inane revelation , and bringing great joy to my life with your humor and en ergy . I also want to thank Remke fo r his expertise in geophysics and regular check - in s , despite being half a planet away. Remke , y ou were my connection in to this program , and I am so glad that you were able and willing to be a part of this from start to finish. The quality of this work benefitted greatly from your attention to detail and broad perspective. None of this would have been possible without the immense support that I received as a student here . I would like to acknowledge the funding sources that fina ncially supported m yself , this research , and many meaningful travel opportunities ; the National Science Foundation, vi National Aeronautics and Space Administration, MSU College of Natural Science , and MSU Department of Earth and Environmental Sciences. I am also immensely grateful to my coauthor and collaborator Dr. Ste phen Hamilton , who facilitated and supported the extensive data collection at the Kellogg Biological Station over the years , and willingly s hared his vast expertise in ecology. I would like to especially thank Julia Michienzi and Amanda Liddle, two impressive young women whom I had the privilege to befriend and mentor along the way . You brought your passion and abilities to this work and br eathed new life into me and my research . I am thankful as well for every member of the EES Departmen t that I had the good fortune to encounter throughout my graduate studies . You valued and encouraged my contributions to our science and community in equal measur e , and never once made me doubt that this was where I belonged . To the staff members who keep things afloa t , thank you for your diligence and care . To my extended Hydrolab family, thank you for embracing me as a person and a scientist, and being my home away from home all these years ; I will treasure these friendships for a lifetime . To my mentors , Maura Hibbitts and Dr. Stanley Radford, thank you for seeing such potential in me as a scientist ; and to my coach Jeff Richards , thank you for teaching me to dig deep and the power of believ ing in one self . Finally, I would like to thank a few very special people who lifted me up and saw me through to the end of this epic journey: Dr s . Leanne Hancock and Sherry M artin , my personal cheerleader s , your steadfast friendship has meant the world ; Brent Heerspink and Quercus Hamlin, my accountability - buddies, you ma k e a virtual dissertation push enjoyable every day ; Leah Bennett and Brian Clemente, my forever friends, you remind me who I am and keep me laughing along the way; and last but not least, my infinitely supportive parents, Diane and Diete r Kuhl , I owe it all to your love and the life you gave me here amongst the water and the trees . vii T ABLE OF CONTENTS LIST OF TABLES i CHAPTER 1: INTRODUCTION ................................ ................................ ................................ ... 1 CHAPTER 2: QUANTIFYING SOIL WATER AND ROOT DYNAMICS USING A JOINT HYDROGEOPHYSICAL INVERSION ................................ ................................ ........................ 6 Abstract ................................ ................................ ................................ ................................ ....... 6 1. Introduction ................................ ................................ ................................ ............................. 7 2. Materials and M ethods ................................ ................................ ................................ .......... 12 2.1 Reference Model ................................ ................................ ................................ ............. 14 2.2 Synthetic Experiments ................................ ................................ ................................ .... 22 2.3 Inversion Scenarios ................................ ................................ ................................ ......... 26 3. Results and Discussion ................................ ................................ ................................ ......... 29 4. Concl usions ................................ ................................ ................................ ........................... 36 Acknowledgments ................................ ................................ ................................ .................... 37 APPENDIX ................................ ................................ ................................ ................................ ... 38 CHAPTER 3: ROOT WATER UPTAKE OF BIOFUEL CROPS REVEALED BY COUPLED ELECTRICAL RESISTIVITY AND SOIL WATER CONTENT MEASUREMENTS ............. 42 Abstract ................................ ................................ ................................ ................................ ..... 42 1. Introduction ................................ ................................ ................................ ........................... 43 2. Methods ................................ ................................ ................................ ................................ 47 2.1 Study Site ................................ ................................ ................................ ........................ 47 2.2 Climate, Soil, and Vegetation Data ................................ ................................ ................. 48 2.3 Electrical Resistivity and Soil Petrophysical Data ................................ .......................... 50 2.4 Coupled Hydrogeophysical Inversion ................................ ................................ ............. 52 3. Results and Discussion ................................ ................................ ................................ ......... 62 4. Concl usions ................................ ................................ ................................ ........................... 76 Acknowledgments ................................ ................................ ................................ .................... 78 CHAPTER 4: A NEW COUPLED HYDROGEOPHYSICAL INVERSION TO ESTIMATE COARSE ROOT FRACTION FROM ERT DATA ................................ ................................ ..... 79 Abstract ................................ ................................ ................................ ................................ ..... 79 1. Introduction ................................ ................................ ................................ ........................... 80 2. Materials and Methods ................................ ................................ ................................ .......... 82 2.1 Site Description ................................ ................................ ................................ ............... 82 2.2 Modeling Methods ................................ ................................ ................................ .......... 85 2.3 Mapping Root Mass Fraction ................................ ................................ .......................... 91 viii 3. Results ................................ ................................ ................................ ................................ ... 92 3.1 Model Calibration and Validation ................................ ................................ .................. 92 3.2 Driver Sensitivity ................................ ................................ ................................ ............ 96 3.3 Mapping Root Mass ................................ ................................ ................................ ...... 100 4. Discussion ................................ ................................ ................................ ........................... 101 Acknowledgments ................................ ................................ ................................ .................. 104 CHAPTER 5: PREDICTING BIOGEOPHYSICAL DRIVERS OF WATER CONTENT AND ELECTRICAL RESISTIVITY HETEROGENEITY IN A SUCCESSIONAL FOREST ENVIRONMENT ................................ ................................ ................................ ....................... 105 Abstract ................................ ................................ ................................ ................................ ... 105 1. Introduction ................................ ................................ ................................ ......................... 106 2. Methods ................................ ................................ ................................ .............................. 108 2.1 Site Description ................................ ................................ ................................ ............. 108 2.2 Instrumentation & Data Collection ................................ ................................ ............... 109 2.3 Model Parameterization ................................ ................................ ................................ 112 2.4 Analysis o f Variance ................................ ................................ ................................ ..... 113 2.5 Modeled ER ................................ ................................ ................................ .................. 114 2.6 Scale - Dependent Water Balance Comparison ................................ .............................. 114 3. Results and Discussion ................................ ................................ ................................ ....... 115 3.1 Hydrothermal Model Validation ................................ ................................ ................... 115 3.2 Analysis of Variance ................................ ................................ ................................ ..... 116 3.3 Modeled ER ................................ ................................ ................................ .................. 119 3.4 Water Balance Comparison ................................ ................................ .......................... 122 4. Conclusions ................................ ................................ ................................ ......................... 124 Acknowledgments ................................ ................................ ................................ .................. 125 REFERENCES ................................ ................................ ................................ ........................... 126 ix LIST OF TABLES Table 2.1. Summary of optimization residuals. ................................ ................................ ............ 31 Tabl e 2.2. Uncertainty in parameter estimates. ................................ ................................ ............. 33 Table 3.1. Summary of soil properties. ................................ ................................ ......................... 52 Table 3.2. Parameter optimization results. ................................ ................................ .................... 64 Table 3.3. Soil interface optimization results. ................................ ................................ .............. 73 Table 4.1. Summary of driver scenarios. ................................ ................................ ...................... 99 Table 4.2. Summary of driver scena rios. ................................ ................................ .................... 100 Table 5.1. Analysis of regression summary table. ................................ ................................ ...... 119 x LIST OF FIGURES Figure 2.1. Coupled hydrogeophysical model schematic. ................................ ............................ 14 Figure 2.2. 1D soil profile model schematic. ................................ ................................ ................ 15 Figure 2.3. Vrugt model root distribution. ................................ ................................ .................... 17 Figure 2.4. Climate and LAI model inputs. ................................ ................................ .................. 19 Figure 2.5. Sensitivity analysis results. ................................ ................................ ......................... 25 Figure 2.6. Optimization r esults. ................................ ................................ ................................ ... 30 Figure 2.7. Residuals of optimization. ................................ ................................ .......................... 32 Figure A.2.1. Petrophysical relationship. ................................ ................................ ...................... 39 Figure A.2.2. High sensitivity parameter estimation resu lts scenario SC3. ................................ . 40 Figure A.2.3. Low sensitivity parameter estimation results scenario SC3. ................................ .. 41 Figure 3.1. Site map. ................................ ................................ ................................ ..................... 49 Figure 3.2. Hydrological model inputs. ................................ ................................ ........................ 55 Figure 3.3. Geophysical model mesh. ................................ ................................ ........................... 60 Figure 3.4. Optimization results. ................................ ................................ ................................ ... 65 Figure 3.5. Model output and validation. ................................ ................................ ...................... 67 Figure 3.6. Modeled root distribution. ................................ ................................ .......................... 69 Figure 3.7. 1D model and observations. ................................ ................................ ....................... 71 Figure 3.8. Modeled vs. obse rved ER. ................................ ................................ .......................... 72 Figure 3.9. 2D pseudo - sections of ER model and data. ................................ ................................ 74 Figu re 3.10. Sensitivity of ET to soil interface. ................................ ................................ ............ 76 Figure 4.1. Site map. ................................ ................................ ................................ ..................... 83 xi Figure 4.2. Petrophysical data. ................................ ................................ ................................ ...... 88 Figure 4.3. Modeled electrical resistivity cross - section. ................................ ............................... 89 Figure 4.4. Core profiles. ................................ ................................ ................................ .............. 93 Figure 4.5. Model validation. ................................ ................................ ................................ ........ 95 Figure 4.6. Sensitivity to model drivers. ................................ ................................ ....................... 98 Figure 4.7. Woody root mass estimate vs. tre e index. ................................ ................................ 101 Figure 5.1. Site map. ................................ ................................ ................................ ................... 108 Figure 5.2. Climate input data. ................................ ................................ ................................ .... 112 Figure 5.3. Soil water content and temperature data. ................................ ................................ . 116 Figure 5.4. Spatial variables for linear regression. ................................ ................................ ..... 121 Figure 5.5. Hydrogeophys ical model outputs. ................................ ................................ ............ 122 Figure 5.6. Spatial variability in water balance. ................................ ................................ ......... 123 Figure 5.7. Effect of resolution on water balance. ................................ ................................ ...... 124 xii KEY TO ABB R EVIATIONS AR Apparent Resistivity DBH Diameter at Breast Height EC Electrical Conductivity ER Electrical Resistivity ET Evapotranspiration KBS Kellogg Biological Station LAI Leaf Area Index TDR Time Domain Reflectometry WRM Woody Root Mass 1 CHAPTER 1: INTRODUCTION One of the greatest challenges facing a sustainable 21 st century is understanding and mitigating the impact of a changing climate and landscape on global and regional water balances. Showing no signs of slowing, i ncreases in atmospheric carbon dioxide concentrations and global surface temperatures have already been linked to shifts in patterns of precipitation and vegetation , the full consequences of which have yet to be realized (Diffenbaugh and Field 2013) . Additionally, humans continue to heavily modify both the vegetation on the landscape and the rate of rec harge of groundwater in order to meet food production , energy , and industrial demands , with large ly unknown implications for water cycling (Rost, Gerten, and Heyder 2008) . T oday more than ever , it is therefore critical that we harness the full power of ever - improving technological capabilities to forecast how these decisions and their feedbacks , both past and future , will ultimately affect the finite natural resources necessary for life on Earth . Doing so will require e nhancing our understanding of , as well a s our ability to accurately model the full suite of physical hydrological processes that determine where water is and where it is going at any given point in time or space. W ater on Earth is partitioned into disp arate pools in the atmosphere , ice, oceans, and the terrestrial landscape (including surface water and groundwater) , with a constant cycling between them via the physical processes of recharge, precipitation, evapotranspiration, and condensation. Evapot ranspiration is the sum of two of the largest pathways by which water is returned to the in the water cycle . E vaporation is the process by wh ich water is removed from the shallow soil surface directly, while transpiration is the evaporation of water from the leaf surface . To meet the demands of photosynthesis, the stomata on the leaf surface 2 o pen to the atmosphere creating a negative pressure g radient down the stem and along the length of the roots. Soil w ater in contact with those roots travels down the gradient against gravity to the leaf surface, and water that is not used for photosynthesis is transpired. Transpiration alone is estimated to account for upwards of 8 0% of terrestrial water fluxes (Jasechko et al. 2013) . However, a lack of constraints on the biophysical drivers of transpiration hinders our ability to model the spatial and temporal variability of this important process at region al and global scales . Primary factors that control the amount and timing of transpiration include the volume of water available in the soil due to climate and drainage conditions, a s well as the extent of the root area to meet the demand of the leaf area . Each of these factors are co - dependent and highly heterogeneous even with the same species , making it difficult to broadly apply assumptions of each in modeling frameworks. While remote sensing products are available for climate, leaf area , and shallow soil water content around the globe , many challenges remain to collect ing root extent data even at the field scale. Studying r oots in - situ is akin to other subsurface characterization problems in that the significant physical limitations of d irect measurements at small scales become impossible to overcome at large scales. However, physical principles that rely on the reflection and refraction of electromagnetic waves can be taken advantage of to infer the properties of th e features that those waves pass through and importantly , are scalable . Much like we use medical imaging technology to study the human body, geophysical methods can be used to , with a wide array of technique s that spans seismology, gravity, radar, and electrical resistivity (ER) . ER methods are based on l aw which establishes that the current and voltage in a circuit are proportional to the resistance . If we consider the Earth as a resistor, and in ject a 3 known current , the measured potenti al is proportional to the resistance of the material the current passed through . For studying vadose zone hydrology, ER methods are particularly useful due to the infl uence of water content and soil texture on the resistivity of the soil. law , the resistivity of a soil or rock matrix is proportional to the porosity and temperature of the subs trate, as well as the conductivity and volume of the inte rstitial pore fluid (Archie 1942) . This relationship c an be used to create spatial and temporal maps of water content in the shallow subsurface. From these maps we can create hypotheses about what hydrological process es are driving the observed patterns in space and time. However, we are limited in our ability to test these hypotheses without further data collection . O ne solution to this conundrum is the coupled hydrogeophysical inversion approach , which reverses the typ ical ER data inversion. In this meth od, the hydrology and temperature dynamics are modeled first using readily available climate, soil, and vegetation data products and used to predict a resistivity distribution . An ER forward model predicts the potential field in the given resistivity distr ibution and is compared directly to the measurements made in the field. Minimizing the difference between the model and observations through optimization can be used to update our understanding of the physical processes that lead to the predicted resistivi ty under the assumption that a well - performing process - based model is an accurate representation of reality . While this minimally invasive approach has been demonstrated to improve the estimation of soil properties in both synthetic and field studies , to date , it has not been to our knowledge applied to estimate root properties. Identifying this gap in the literature, the overarching goal of this dissertation was two - fold; first , to test the sensitivity of a coupled hydrogeophysical inversion method to estimate root properties in a synthetic study; and second, after establishing 4 the feasibility of the approach, to apply it to infer root properties and evapotranspiration rates within two distinct ecosystems. These two ecosystems we re bot h located in southwest Michigan, USA on the property of the Kellogg Biological Station and Michigan State University and consisted of 1) a heavily managed cropping system evaluating ten different experimental biofuel treatments , a nd 2) an unmanaged transect covering the full extent of forest succession . In the process of developing and applying this novel method , the following research questions were addressed and are presented in the four subsequent chapters of this dissertation. Chapter 2 addresses the question of whether the devised approach has the necessary sensitivity to the root properties of interest , including the depth and distribution of roots below a stand of maize unde r the climate and soil conditions present at our biofuel study site . After validating that t he appropriate sensitivity for parameter estimation can be derived from this approach , Chapter 3 proceeds to apply it to estimate the root properties of two potenti al biofuel treatments: an annual crop rotation of canola, maize, and soybean , and the perennial biofuel grass miscanthus. C hapter 2 is published in Vadose Zone Journal (Kuhl et al. 2018) and has been reprinted here . C hapter 3 is currently being re viewed at Vadose Zone Journal. While the intention was to subsequently apply this same approach to the forest succession site, challenges arose due to increased heterogeneity in the forest portion of the transect that could not be explained by the model as it were . Chapter 4 therefore explores the theoretical implementation of resistive coarse woody root mass in the model framework to explain the observed high resistivity anomalies that were observed in the forest but absent in the grass portion of the transect during the non - growing season . Chapter 5 uses a statistical linear multi - regression analysis to determine the significance of other spatial variables, including surface 5 topography, leaf area index, and above - ground biomass to explain the spatial and temporal variability in the observed shallow electrical resistivity throughout the growing season . Taken t ogether , this body of work demonstrates the value of this novel coupled hydrogeophysical inversion approach to make inferenc es and test hypotheses about root water uptake processes that goes beyond traditional ER imaging . In addition to establishing the utility of the approach to estimate important parameters like the root depth and crop coefficient in an agronomic setting, we also found this approach had sensitivity to the proportion of coarse woody root mass in the upper 0.5 m of soil . New t ools such as this enables the study of root extent in - situ , improves our understanding of the biophysical drivers of evapotranspiration, and provides a method to quantify the impacts of land use change on the water balance. 6 CHAPTER 2: QUANTIFYING SOIL WATER AND ROOT DYNAMICS USING A JOINT HYDROGEOPHYSICAL INVERSION A bstract Plot - to field - scale root distribution data are relatively rare, and difficult to measure with traditional methods. Nevertheless, these data are needed to accurately model root water uptake (RWU) processes within agronomic, hydrologic, and terrestrial biosphere models. New tools a re needed to effectively observe root distributions and model dynamic root growth processes. In the past decade, geophysical tools have increasingly been employed to study the vadose zone, and hydrogeophysical inversions have shown promise to noninvasively characterize water dynamics. In such an approach, the hydrology is modeled and hydrological data is inverted with the geophysical data, constraining the geophysical inversion results, decreasing uncertainty and the number of non - unique solutions. Here, we develop and test a coupled hydrogeophysical inversion approach that uses electrical resistivity data to estimate soil hydraulic, petrophysical, and root dynamic parameters. This builds on prior research that either used a coupled hydrogeophysical inversio n to estimate soil hydraulic parameters only, or a hydrological inversion to estimate root distribution or root water uptake parameters. Our results indicate that under the conditions tested, this approach accurately captures root water dynamics and soil h ydraulic parameters. This opens up opportunities to non - invasively image a variety of root distributions and soil systems, better understand the dynamics of RWU processes, and improve estimates of transpiration for systems models. 7 1. Introduction Transpi ration is the most important pathway for the exchange of water from Earth to the atmosphere, accounting for up to 80% of terrestrial evapotranspiration (Jasechko et al. 2013) . Thus, disruptions to the plant community through climate and land - use changes wi ll likely have serious implications for regional to global water balances. To predict and mitigate the effects of those changes, agronomic, hydrologic, and terrestrial biosphere models must accurately capture the exchange of water along transpiration pathw ays. Doing so requires understanding the underlying processes that drive such exchanges. The interdependent and dynamic nature of the factors controlling transpiration, and our inability to observe the processes directly, makes transpiration challenging to appropriately represent in these models. Transpiration is fundamentally controlled by root distributions and root water uptake (RWU) processes, yet due in part to a lack of dynamic root function data, these processes are often oversimplified in models (W arren et al. 2015) . Such data is rarely available since it is challenging to observe roots in the field, especially through time (Cai et al. 2017) . Direct approaches such as excavation and root windows are limited at the field scale, and are very costly an d labor intensive. These approaches are also not as feasible for deep roots associated with woody plants, such as trees (Maeght, Rewald, and Pierret 2013) . Non - destructive methods are thus needed to understand plant functions as a response to changing conditions in a range of field settings. One viable approach is to use changes in root - zone soil moisture as a proxy for the presence of RWU processes. Exist ing methods to measure or estimate changes in soil moisture have limitations. Larger scale methods that rely on measuring atmospheric fluxes, such as eddy covariance, can get measurements at the stand level, but it is difficult to isolate transpiration fro m evaporation 8 fluxes. Below ground, tools such as time domain reflectometry, neutron probes, and capacitance probes estimate soil moisture with high temporal resolution, but they are intrusive, and barring significant installation efforts, lack the necessa ry spatial resolution to capture heterogeneities in soil properties and root densities. Lysimeters can measure the drainage out of the profile, but are costly and labor intensive, and cannot capture sub - layer dynamics (Schelle et al. 2012) . Electrical resistivity (ER), a minimally invasive geophysical technique, measures the current - induced potential field underground. ER data, comprised of measured potentials at varying dipole lengths and distances, are strongly influenced by not only the po rosity of the soil, but the saturation and electrical conductance of the pore water. A petrophysical relationship equates these variables to the bulk resistivity calculated from the measured potentials (Archie 1942) . Thus, if the resistivity and pore water conductivity are known, soil moisture can be estimated via empirical relationships fit with petrophysical parameters (PP). ER surveys are particularly well suited to investigate hydrological problems because they provide a bulk measurement influenced by a volume of media surrounding the electrodes, the dimension of which can be varied with the electrode array geometry. For example, an ER survey could be designed to have increased sensitivity to the upper 0.5 m where most RWU activity is concentrated. For c ontrast, traditional discrete methods of sampling soil moisture represent conditions at a single point, and are susceptible to both over - and under - representation of volumetric soil moisture due to features such as textural layering and macropores. Further , the behavior of current flow makes it ideal to study the subsurface in multiple dimensions, unlike point measurements that require a large installation effort to capture lateral variability. Previous research has used the relationship between ER and soi l moisture to infer the presence of roots and RWU dynamics in multiple dimensions, and at high spatial and temporal 9 resolutions (e.g., Beff et al. 2013; Junliang Fan et al. 2015; Garré et al. 2011, 2013; Jayawickreme 2008; Jayawickreme, Van Dam, and Hyndma n 2010; Ma, Van Dam, and Jayawickreme 2014; Michot et al. 2003; Robinson, Slater, and Schäfer 2012; Whalley et al. 2017) . These studies applied a traditional ER inversion method to retrieve the soil moisture distribution, wherein a potential field model is optimized to fit the measured potentials from the ER survey, and translated to a static soil moisture distribution via a petrophysical relationship. However, such traditional ER data inversions may result in non - unique and unconstrained solutions and prod uce physically unrealistic soil moisture distributions (Mboh et al. 2012) . To improve the application of ER methods to hydrological problems, researchers have proposed coupling the geophysical model with a site - specific hydrological model (Hinnell et al. 2010; Minsley et al. 2011; Singha et al. 2015) . One approach to coupled hydrogeophysical inversion involves: 1) forward modeling transient water fluxes, 2) converting the final modeled soil moisture distribution into an ER distribution using a petrophysica l relationship, 3) forward modeling the potential field to compare the modeled and measured potentials, and then 4) updating coupled model parameters to minimize differences between the observed and modeled ER data (Mboh et al. 2012) . In such an approach, changes to the soil hydraulic parameters affect the soil moisture distribution via the hydrologic model; this change in simulated soil moisture alters the modeled electrical potential field. There are three categories of parameters that are typically unkno wn in near surface hydrogeophysical problems: petrophysical parameters (PP), soil hydraulic parameters (SHP), and root parameters (RP). PP describe the relationship between simulated soil moisture and subsurface ER. SHP affect the soil water dynamics via w ater retention and infiltration models, and are widely referred to as soil hydraulic parameters. RP control modeled root physiology, 10 including the growth and distribution of the roots, along with the relationship between potential transpiration, soil moist ure, and RWU. Several studies have used a coupled hydrogeophysical inversion approach to successfully estimate SHP and transient soil moisture, in both synthetic (Hinnell et al. 2010) , and field (Mboh et al. 2012; Moreno, Arnon, and Furman 2015; Thomas et al. 2017; Tran et al. 2016) scenarios. Mboh et al. (2012), used ER data from a short (several hour) inflow experiment to estimate SHP. They found the ER data - only inversion more robust than just using cumulative inflow data, and slightly worse than combin ing both data types in the objective function. Despite this, ER data without supporting hydrological data has been demonstrated to be sensitive enough to soil moisture dynamics to reasonably estimate SHP. Using a grain - size analysis and the Rosetta databas e (Schaap, Leij, and Genuchten 2001) to initialize the SHP, Moreno et al. (2015) used nine geophysical surveys throughout a year - long period to fit select SHP in a two - layer soil. The RP in this particular study were fixed (held constant) at reference valu es. Despite the interest in this area of research, few if any studies have attempted to use a coupled hydrogeophysical inversion approach to characterize RWU dynamics. RWU models, which include the root distribution and water stress functions for RWU redu ction, have to our knowledge been parameterized with root or hydrological data inversions only. Hupet et al. (2003) used neutron probe water content data to estimate model parameters, including SHP for a 1D homogeneous soil, and rooting depth, and root len gth density for a maize crop with mixed success; the RP were less well constrained in medium textured soils. Schelle et al. (2012) completed a similar study, using daily lysimeter and matric potential data to estimate SHP in two layers, as well as a root d istribution parameter. Using neutron probe soil moisture data, Vrugt, Hopmans, and Simunek (2001) calibrated multiple parameters of a new flexible 2D root 11 distribution model and some SHP, with excellent agreement between measured and modeled soil moisture in 2D. Recently, Cai et al. (2017) demonstrated the use of a similar method to parameterize several different water stress functions for 1D RWU to fit field observations of daily water content, with good agreement between modeled and observed root length d ensity distributions. Prior studies provide convincing evidence that separately: 1) RWU models can be parameterized using inverted water content data (e.g. Hupet et al., 2003), and 2) ER data serves as a proxy for water content data, and can be used to e stimate SHP (e.g. Hinnell et al., 2010). Therefore, we hypothesize that RWU models can be parameterized in a coupled h ydrogeophysical inversion using ER data, which to - date has not been directly investigated. To test this hypothesis, we develop a method th at builds upon these two established concepts by estimating SHP and RP using a coupled hydrogeophysical inversion of ER data. We also incorporate the estimation of PP in the relationship between soil moisture and resistivity. Given the challenges posed by the lack of transient root data from the field, we sought to validate this novel approach with a synthetic 1D proof - of - concept study. In addition to validating the method, we sought to identify the limitations that sparse data and parameter measurement and sensitivity errors might impose on this approach. Informing the ER with the known hydrology of the site reduces the non - unique solutions to only those that match a physical reality. Additionally, when soil water fluxes are modeled, losses from evaporatio n can be distinguished from transpiration, as well as hydraulic redistribution. This is a particular benefit over time - differential ER methods which can only provide total gains or losses in soil moisture. Furthermore, the embedded hydrological model can t hen be used to simulate soil moisture conditions prior to, between, and after ER survey events. 12 Using ER methods in this fashion to calibrate a process - based model also allows for the model to be used to forecast soil moisture fluxes under hypothetical fut ure climate conditions. As with standard ER inversions however, the petrophysical relationship must be known, which remains a challenge to any geophysical approach used to estimate hydraulic properties (Laloy et al. 2011) . In addition, the behavior of som e SHP parameters is not independent of others, making it difficult to estimate a full spatial distribution of values. Roots themselves can also affect the ER signal, however this has primarily been shown for much larger tree roots, whose resistance is dist inguishable from the surrounding material (Amato et al. 2008a) . While there are added complexities with real field settings related to lateral heterogeneity, this does not necessarily preclude the establishment of a representative hydrological model. In this paper, we establish a framework for a robust, minimally invasive, and cost and labor efficient way to calibrate the many parameters of site - specific hydrological models. We present an overview of the model components used to build the algorithm, te st the sensitivity of the parameters, and perform a series of synthetic 1D modeling experiments to validate the algorithm. To ensure the synthetic study reflects realistic conditions, we use climate data and measured soil profile characteristics from a stu dy site at the Kellogg Biological Station in southwest Michigan, USA, described in more detail below. The approach detailed here is universally applicable, and provides a path to investigate heterogeneous root and soil systems in two and three dimensions with limited a priori information. 2. Materials and Methods We developed a coupled hydrogeophysical inversion algorithm that: 1) simulates the movement of water throughout the soil profile, 2) converts snapshots of the transient soil moisture distribution to soil resistivity using a petrophysical relationship, 3) simulates the 13 potential field and ER data using a forward resistivity model, and 4) optimizes the parameters of the models by minimizing the difference between the modeled and measured ER data. Th is approach can estimate the SHP, RP, and PP that are often challenging to directly measure, either in - situ or with lab bench experiments. Our algorithm simulates the hydrogeophysical processes using four publicly available codes (Fig ure 2 .1 ): the System A pproach to Land Use Sustainability ( SALUS ) model (Basso et al. 2006) for potential evaporation and transpiration; HYDRUS for root growth (Hartmann et al. 2018) , RWU (Feddes et al. 2001) , hydraulic redistribution, v ariably - saturated hydrology (Richards 1931) , snow hydrology and heat transport (Chung and Horton 1987) ; FWD2_5D (Pidlisecky and Knight 2008) for the electrical potential forward calculations, and the global optimization shuffled complex evolution algorithm , SCE - UA (Duan, Sorooshian, and Gupta 1992) for parameter estimation. Each of these four models and algorithms are described in more detail below. To demonstrate and validate this inversion algorithm, our experiment involves three steps in which we: A) for ward run the algorithm with a set of reference parameters to generate a reference hydrological model and synthetic 'measured' ER data, B) test the objective function sensitivity to SHP, PP, RP, and C) test the inversion algorithm under six scenarios that t est the influence of variations in data density and parameter uncertainty. We evaluate each of the six inversion results relative to the synthetic reference parameters, soil moisture, root distribution, and RWU data. 14 Figure 2.1. Coupled hydrogeophysical model schematic. Schematic of the 4 - step coupled hydrogeophysical inversion algorithm that estimates petrophysical parameters (PP), soil hydraulic parameters (SHP), and root parameters (RP) to minimize the root mean square error (RMSE) b etween measured and modeled electrical potentials. Boxes contain the components of the algorithm, while labeled arrows describe the flow of output from one component to the other. Model codes used include the System Approach to Land Use Sustainability crop model (SALUS), a 1D hydrological model (HYDRUS), a electrical resistivity forward model (FWD2_5D), and a Shuffled Complex Evolution optimization model developed at the University of Arizona (SCE - UA). 2.1 Reference Model We developed a realistic plot - sc ale 1D vertical model of maize over a three - layer soil, assuming uniform soil and root properties laterally (Fig ure 2 . 2). We based this model on a test plot at the Kellogg Biological Station Great Lakes Bioenergy Research Center site in southwest Michigan, USA (described previously by Zenone et al. (2013) ). Actual data from this site for climate conditions, sediment grain size distributions, soil moisture, soil temperature, and petrophysical relationships were used for this synthetic model. This plot was al so instrumented with electrodes for ER surveys to be used in future studies. 15 Figure 2 .2. 1D s oil p rofile m odel s chematic. Schematic diagram of the model profile with soil layers and plant roots overlain by; a local weather station (with air temperature, wind speed, solar radiation, and precipitation), three buried temperature sensors; and 30 electrodes. Hydrological modeling unsaturated flow. Some HYDRUS inputs, including potential evaporation and potential transpiration, were estimated using the SALUS crop model (Basso et al., 2006), and imported to HYDRUS. SALUS i s a daily water and nutrient balance dynamic vegetation model that uses daily climate data to calculate potential evapotranspiration using the Priestly - Taylor equation. It also models leaf - area index to differentiate evaporation from transpiration (Ritchie 1998) . Daily values of potential evapotranspiration were disaggregated proportional to hourly modeled sun Calculator. HYDRUS was selected over SALUS for modeling RWU along w ith water and 16 energy fluxes because while SALUS computes daily temperature and water balance within the root zone, HYDRUS has a finer vertical discretization and incorporates a more sophisticated hydrology algorithm. Daily soil moisture output from SALUS w as also deemed insufficient due to the sensitivity of ER to diurnal changes in soil moisture (Robinson et al., 2012). For the purposes of the study we assumed no error in the climate data inputs, soil temperature, leaf - area index, potential evapotranspira tion calculation, or soil layer boundaries. We assumed that hydraulically - significant soil layering could be accurately identified from in - situ textural classification, and that for this plot - scale study these layer boundaries are horizontal. However, the model could be adapted to estimate layer boundaries by allowing layer depths to vary. Weather stations are widely available at high spatial resolutions making it easier to model temperature dynamics and potential evapotranspiration. Crop models that focus on modeling yield such as SALUS are also well suited to modeling the leaf - area index. We chose to describe the root distribution in HYDRUS with the Vrugt model (Vrugt et al., 2001) because of its flexibility. The parameters of the reference Vrugt model we re estimated to fit the normalized root distribution output from SALUS at the end of the growing season using an unconstrained nonlinear optimization (Fig. 3). The Vrugt model: (1) where is depth; and are fitting parameters that alter the rate of exponential decay, and set the depth of maximum root density, respectively. A smooth exponential root distribut ion can be obtained by setting to zero. The calibrated root distribution ( and , Fig ure 2. 3) is similar to the root density distributions observed for maize (Tardieu et al., 1988, Jackson et al., 1996). Since the decay of the expone ntial root distribution is largely dependent on the ratio of 17 (Eq uation 1 ), we only estimated in the optimization. An analysis of different and values found that multiple combinations could closely replicate the observed root dis tribution. Assuming the value would not be known in reality, and it is the ratio of that primarily influences the shape of the curve, was fixed at 1.5 m for the tested scenarios. Although we chose the Vrugt root distribution model, any of the four root distribution models available in the HYDRUS root growth module could be parameterized with this coupled inversion algorithm. In addition, the models themselves could be within the search space of the optimization, allowing for even gr eater flexibility. Figure 2 .3. Vrugt m odel r oot d istribution. The Vrugt equation (Equation 1) solved for depths up to 1.0 m with the perturbed initial parameters (orange dash) versus those parameterized to fit the output from the SALUS model (yellow star). Note the dependent variable is on the x - axis. The Feddes model (Feddes et al., 2001) was used in HYDRUS to simulate the reduction of potential transpiration during periods of water stress, which occur when the soil moisture is outside of a prescribed range. Water stress in the Feddes model is determined by four pressure 18 head thresholds; h1, h2, h3, and h 4 . In the ideal pressure head range, between h2 and h 3, water is extracted at the potential transpiration rate. When the pressure head is beyond the ideal range, (above h2 or below h3 ), the potential transpiration is reduced by a reduction coefficient proportional to the increase or decrease in pressure head. Beyond the pressure head limits for water extraction, h1 and h4 , the reduction coefficient is 0, resulting in zero RWU. The parameter h2 can be dependent on the texture o f the soil, while h3 can have an upper and lower threshold (denoted with a subscript H or L, respectively). Reference RP values for maize, h1, h2 (soil - layer dependent) , h3 H , h3 L , and h4 (Table S2) were taken from (Wesseling et al. 1991) . Grain size ana lysis from samples at the site provided percent sand, silt, and clay, and bulk density at 0.1 m intervals. From this analysis, three distinct layers, and the boundary locations between them were identified at 0.4 and 0.8 m. The soils here are described as well - drained alfisols with a silt loam layer over loamy sand, underlain by coarse sand (Fig ure 2. 2). The soil is classified as a Typic Hapludalf. The van Genuchten - Mualem model (van Genuchten 1980) was implemented in HYDRUS to model soil water retention an d infiltration: ( 2 ) ( 3 ) ( 4 ) ( 5 ) where: and are the residual and saturated water content respectively, is a tortuosity factor (held at 0.5 for this experiment), and are fitting parameters , is the effective water content, and and are the saturated and variably saturated hyd raulic conductivities, respectively, for each soil layer. These parameters, unique for each soil texture, comprise the SHP estimated in the inversion. Reference Van Genuchten - Mualem parameters 19 were extracted from the Ro setta database (Schaap et al., 2001) corresponding to the percent sand, silt, and clay in each layer, and input to the HYDRUS model. The 1D HYDRUS model was run hourly for 12 months starting in Nov 2009, and spatially discretized at 0.02 m intervals to a depth of 1 m, then increasing geometrically by a factor of 1.25 to a depth of 6 m, for a total of 73 vertical nodes. The lower and upper boundaries of the model were controlled by free drainage and atmospheric conditions, respectively. Hourly precipitation and air temperature inputs to HYDRUS were obtained from an adjacent Michigan EnviroWeather Network station at the Long - term Ecological Research site (Fig ure 2. 4). Each HYDRUS simulation, for the reference, sensitivity, and test scenarios, was initialized with the same soil moisture distribution which reflects realistic field conditions for the model start date. Figure 2.4. Climate and LAI model inputs. Measured maximum daily air temperature (orange) and daily precipitation (blue), with SALUS - modeled le af - area index (purple) throughout the model year. The ER survey dates comprising the synthetic datasets for the Data Density scenarios (DD1, DD2, DD3), are marked in the bottom - right. All Site Characterization (SC) scenarios as well as the sensitivity anal ysis share the same survey dates with DD1. 20 Once the hydrological modeling was complete, static geophysical models were created at the desired sampling frequency (survey dates are shown on Fig ure 2. 4). The geophysical modeling consists of three steps (Fig ure 2. 1): A) converting simulated soil moisture to ER at a standard temperature (25 °C), B) correcting for distributed subsurface temperatures on the survey date (Fig ure 2. 4), and C) forward - simulating the potential field for a synthetic ER survey. Ste ps A - C are repeated for each survey date comprising the ER dataset. Laboratory resistivity experiments were run to calculate an empirical relationship between soil moisture and ER measurements. Sample material was placed in a 22.2 x 4 x 3.2 cm soil box wit h metal plates for current transmission at the far ends. The resistance was measured using two electrodes near the center of the box while gradually increasing the water content in the sample (see details of this lab approach in Jayawickreme et al. (2010)) . This was repeated for nine soil samples, collected from each of the three soil layers at three adjacent plots. The power - law petrophysical relationship: ( 6 ) where: is the resistivity at depth , is the water content at dep th and and are soil - specific empirical PP, was fit to the observations for each soil layer, (Fig ure A.2. 1). For this study, we made the simplifying assumption that there was no additional grain surface conductivity term in the petrophysical relationship for any soil layer. Given the frequency of precipitation events that flush out salts accumulated through evapotranspiration, we also assumed that the pore - water conductivity did not undergo changes with time (Jayawickreme et al. 2010). However, this coupled inversion could be modified to incorporate solute transport within the hydrological model, allowing the petrophysical relationship to be transient. 21 Subsurface resistivity is temperature dependent, thus it is important to accou nt for the below - ground temperature gradient. The site was also instrumented with high temporal resolution (every 2 hours) temperature sensors (Thermochron iButton DS1922L) at three depths (0.26, 0.66, 1.20 m). We generated realistic soil temperatures for the synthetic study using a heat transport model for soil temperatures at all nodes. Heat transport parameters for the Chung and Horton soil temperature model (Chung and Horton, 1987) in HYDRUS were calibrated using soil and air temperature data from Febru ary to May 2010. The parameters were iteratively estimated for layers one, two, and three, consecutively, with the HYDRUS model (non - growing season only) in inversion mode. ER modeled at the standard 25°C was corrected to ER at the modeled temperature for each node, using the linear model from (Hayley et al. 2007) , . ( 7 ) Once the temperature - corrected ER distribution is calculated, a synthetic ER survey can be modeled. The electrode array used to generate the synthe tic data mirrors the installation at the Kellogg Biological Station site, and is comprised of 30 surface electrodes, spaced 0.3 m apart. This array length was chosen to focus on the upper 2 m of the profile. Due to the assumptions about lateral homogeneity in the model, the 13 unique electrode geometries in dipole - dipole configuration were modeled with the FW2_5D code (Pidlisecky and Knight, 2008) creating a 1D profile of effective measurement depths. The electrical potential field forward modeling code, FW 2_5D, allows the input ER distribution to vary in the x and z directions, assuming homogeneity in the y direction. Since our hydrological model only varied in the z direction, we extrapolated our calculated 1D ER laterally. The FW2_5D model was discretized such that errors for our electrode array were 0.7% on average from the analytical solution for a homogenous Earth, yet still efficient. The x dimension grid spacing was 0.15 m (half the distance between 22 electrodes) over the span of the center of the model , increasing by a factor of two in each direction from center for a total model length of 24 m. The vertical discretization was the same as that described for the HYDRUS model above. To complete the coupled hydrogeophysical inversion, we use the synthetic ER data from each survey to estimate the SHP, PP, and RP through optimization. The robust SCE - UA global optimization algorithm minimizes an objective function : ( 8 ) between the modeled ( ) s ynthetic ER data by evolving the coupled model parameters. As the soil moisture profile output from the hydrological model is updated with each optimization iteration, the soil ER distribution also changes according to the petrophysical model. This in turn alters the modeled potential field and therefore the modeled ER data. This optimization continues until the resulting modeled ER data matches the measured ER data within the optimization specified closure criteria. SCE - UA blends several traditional optimi zation approaches to efficiently find the global minimum. It starts by randomly sampling a large population from the parameter space. This population is then divided amongst a number of complexes, which each evolve independently towards a minimum. After a given number of evolutions, the complexes shuffle, divide, and begin evolving again. This is repeated until convergence is reached. The size of the population and the number of complexes can be scaled with the complexity of the problem. 2.2 Synthetic Exper iments To test the ability of the coupled hydrogeophysical inversion algorithm to estimate SHP, PP, and RP from ER data, we created synthetic ER datasets using a single forward run of the reference model. The reference inputs and parameters used to make the synthetic data were then 23 in subsequent analyses. To increase the difficulty of the optimization, randomly distributed noise (±0.5%) was added to all synthetic ER data. This relatively low amount of error was chosen so that the impact of the tested data density and parameter error would be clearer. It also reflects the low amount of error (median 0.12%) that was observed between reciprocal measurements for this array in preliminary data analysis from the field site. We next tested the sensitivity of a bi - weekly (every two weeks) ER survey dataset to each parameter. Finally, we conducted a series of inversion experiments to retrieve the reference SHP, PP, and RP via iterative optimization starting from some perturbed initial estimate of those parameters. 2. 2.1 Parameter Sensitivity T o determine which parameters to estimate, which to fix, and how to perturb parameters for each inversion scenario, we conducted a sensitivity analysis using forward runs of the PP. Recognizing some co - dependency exists between parameters, particularly with (Huisman et al. 2010; Moreno, Arnon, and Furman 2015) , we expected the global optimization algorithm to overcome this in minimizing the objective function. For this analysis , the algorithm was run sequentially 20 times for each parameter, with values distributed evenly across a predetermined parameter space, while the rest of the parameters were held fixed at their reference values. The ER data selected to test the sensitivit y was a bi - weekly dataset extracted from the reference model run at a 14 - day interval from late April 2010 through Sept 2010 for a total of 12 surveys, comprising 13 measurements each (156 measurements). We quantified the sensitivity of each tested value b y computing the root mean square error (RMSE) between the reference synthetic ER data and the perturbed ER data modeled at the same bi - weekly frequency. 24 The upper and lower limits for each parameter space were chosen based on either: 1) physical limits f or parameters with narrow realistic ranges such as with , , , and (Schaap et al., 2001); or 2) fixed ranges from the reference for parameters with greater variability. For example, upper and lower bounds for were set plus or minus two orders of magnitude around the reference value, given the large uncertainty of that parameter. The petrophysical parameter limits were set ±10% the value used in the reference run, based on variability observed in the field data (Fig ure A.2. 1). Limits on the Feddes RP were set to ±30% of the reference value to see if sensitivity to those parameters existed at levels beyond what has been previously tested (perturbation of 1% in Hupet et al., 2003). Bounds for to vary from 0.5 to 1.5 m. 2. 2.2 Sensitivity Analysis Results The results from the sensitivity analysis show that as parameters deviated from their reference values, RMSE between the reference ER data and the perturbed ER data increased nearly linearly (Fig ure 2. 5). SHP from the first soil layer were significantly more sensitive than those from the deeper layers, with the exception of which had only slightly increased sensitivity in the first layer. As in prior research (Mboh et al., 2012), we found very little sensitivity to the Feddes RP relative to SHP, with the exception of h4 . This was not surprising, as h 4 controls the point at which RWU ceases in very dry conditions. However, we observed that the RP had roughly the same sensitivity as the PP (note that had a much larger parameter space). 25 Figure 2.5. Sensitivity analysis results. Sensitivity analysis results for a) root parameters (RP), b) petrophysical parameters (PP), and c) soil hydraulic parameters (SHP) in all three layers. Th e two sensitivity thresholds at 0.03 and 0.1 are shown as horizontal grey and black dashed lines, respectively. Note the y - scale for root mean square error (RMSE) varies by parameter. The range of tested values corresponds to the upper and lower limits use d in the parameter estimation, which varies by parameter. Most RP had no sensitivity in the tested range and are not shown. A threshold RMSE value (Eq uation 8) of 0.1 was chosen to isolate the highest sensitivity parameters from the rest, which included all the SHP and PP from the first soil layer. The eight high sensitivity parameters that fell into this category included; SHP: , , , , , PP: , , (numerical subscripts denote soil layers 1, 2, 3, where applicable) as well as the RP parameter, . A secondary threshold of 0.03 was chosen to include the remaining parameters that had lower sensitivity. This category also included eight parameters; SHP: , , , , and the RP parameter, h4. Th e remaining RP parameters (all Feddes except h4 , not shown), 26 the SHP: and , and the PP: , , were determined to have essentially no sensitivity as perturbing them had no discernible effect on the simulated ER values. 2. 3 Inversion S cenarios To test the ability of the inversion to estimate SHP, PP, and RP under different conditions, we evaluated a suite of six synthetic data inversion scenarios. These experiments were designed to test the robustness of the parameter estimation to decreasing da ta quality and quantity. We compared model output and parameter estimates for inversions using three data density (DD) types, plus three levels of site characterization (SC) with the same data density. To test the dependence of the optimization on data de nsity, we created three scenarios mimicking three frequencies of field surveys (Fig ure 2. 4). The first scenario, DD1, tests a relatively long - term (six month) ER dataset, with moderate survey frequency (bi - weekly). The dataset for DD1 was used for the para meter sensitivity analysis discussed above and all three SC scenarios (described below); it has 12 surveys with 13 measurements each, for a total of 156 measurements from April to September. A second scenario, DD2, tests the effectiveness of a lower freque ncy (monthly) ER dataset over the same six month period as the dataset in DD1. The data for this scenario was extracted from the reference run at a four week interval, for a total of six surveys with 13 measurements each (78 total measurements). This datas et is half the size of the other two DD scenarios. A third scenario, DD3, tests a short - term (one month) high - frequency (every three days) ER dataset, to simulate the effect of using data exclusively from the peak growing season, but spanning less diver se seasonal weather conditions . This dataset covers a large precipitation - infiltration event, along with the subsequent transpiration - dominated drying period. The data was 27 extracted from the reference run at a three day interval for the period from July to Aug 2010, for a total of 12 surveys, again comprised of 13 measurements each (156 measurements total). In all three DD scenarios, the eight highest sensitivity parameters are estimated, while the remaining low and no sensitivity parameters are fixed at th eir reference values. To test the robustness of the model to various levels of site characterization, we again estimated the high sensitivity parameters, but introduced varying levels of error into the other parameters. These three scenarios include: SC1 , where the low and no sensitivity parameters are assumed well known, and are fixed with minimal error (1.5% - 40%) relative to reference values; SC2 , where these parameters are moderately well known, and fixed with higher error (2.8% - 81%); and SC3 , where the low sensitivity parameters are assumed completely unknown and included in the parameter estimation along with the high sensitivity parameters. In this scenario the no sensitivity parameters are fixed at their reference values. All three SC scenarios use th e same synthetic dataset with bi - weekly ER surveys as that in DD1. In each of the six scenarios, the same starting values were used for the high sensitivity parameters being estimated. With the exception of and the PP, and , which were star ted +/ - 5%, all other parameters were started +/ - 30% away from the reference value. To determine how much fixed error to add to the low sensitivity SHP for SC1 and SC2, we considered the uncertainty in estimating each parameter of the retention function with standard field methods. This was analyzed by Baroni et al. (2010) , who provided a table of standard deviations for each parameter of a three - layer soil. We used the standard deviation from soil layer 3 in their analysis, which was most similar to the sandy soils at our site. For the PP in our model, we calculated the standard deviation across the three replicate experiments that were conducted at the study site for parameters and in each soil layer. 28 Both the low and no sensitivity parameters were fixed at a value half a standard deviation (for SC1) or a full standard deviation (for SC2) from the reference value. The large variation in uncertainty across parameters meant that some were fixed closer to the reference values than others. For example, in SC2, parameters with less uncertainty such as , were fixed 3% from the reference values, while high uncertainty parameters like were fixed 72% away. The fixed values used for these scenarios are shown in Table S2. SC3 tested the feasibility of estimating both the high and low sensitivity parameters (a total of 16). This last scenario was designed to replicate the most likely field scenario, where without conducting a relatively intensive lab bench exper iment, no a priori information would be available for the SHP. The goal was to determine if the inversion algorithm was robust to local minima from the large number of unknowns. If so, the most critical parameters affecting root water dynamics and water fl uxes would still be accurately estimated. For the SCE - UA algorithm, we selected a number of complexes almost equal to the number of unknowns. For all scenarios except SC3, seven complexes were used. Optimization continued until five consecutive shuffles (~ 200 iterations each) did not improve the objective function by 0.1%. The convergence criteria was generally met after ~4000 iterations. The number of complexes in the SCE - UA algorithm was doubled to 14 for the SC3 scenario due to the larger number of free parameters. In that case, the stopping criteria of 10,000 iterations was reached before convergence, however three consecutive shuffles had not improved the objective function at that point. Final estimated parameter values and model outputs were from the iteration with the lowest objective function value. Our study is particularly interested in root water dynamics and therefore we chose to quantify the success of our approach with four outputs from HYDRUS that capture: the root 29 fraction with depth, the cu mulative RWU, the soil - water retention function (Eq uation 2 ), and transient soil moisture. Each inversion scenario is validated by calculating the RMSE between the synthetic model and the hydrologic model with estimated parameters. 3. Results and Discussi on The capability of the hydrogeophysical inversion algorithm is graphically illustrated in Fig ure 2. 6. For this figure, scenario SC3 is chosen to best demonstrate the robustness of the algorithm when all 16 parameters were estimated. The impact of the ini tial conditions on the root water and soil moisture dynamics can be observed in comparison to the reference model in Fig ure 2. 6. The hydrological model with the starting values has a deeper root distribution, more RWU, and a steeper water retention curve. Regardless of the significant error in starting values, the inversion algorithm estimates the root dynamics and the soil moisture curve well, even prior to the start of the ER surveys. 30 Figure 2.6. Optimization results. HYDRUS output of: a) root fracti on, b) RWU from the growing season, c) van Genuchten - Mualem model retention curve, and d) soil moisture at 0.2 m with ER survey dates for SC3 in green. The estimated values from scenario SC3 are plotted versus initial values and reference values. The resul ts from the other 5 scenarios were very similar to SC3 and are not shown here. Note the HYDRUS outputs were not used in the parameter estimation process, rather they are only used to validate the approach. Regardless of the scenario, there was excellent ag reement between the reference and optimized model output (Fig ure 2. 7, Table 2. 1). Even in cases where fixed error in other SHP and RP parameters was large, the inversion found a unique solution for the root distribution within 10% of the reference distribu tion. We also observed the cumulative RWU was estimated 31 within 0.01 m for all scenarios with the exception of SC2, which was off by 0.015 m (Fig ure 2. 7). The transient soil moisture distribution was also closely matched, despite some non - uniqueness in the SHP. We observed no residuals greater than 0.05 cm 3 /cm 3 water content, even in the worst performing scenario, SC2 (Fig ure 2. 7). The soil - water retention function, which is calculated with the estimated SHP (except which is not used to calculate so il - water retention), was also estimated within 0.02 cm 3 /cm 3 of the reference model for each scenario. In Scenario SC3, the SHP perturbed by 30% to start were estimated within 5% of their reference value with the exception of . In layers 2 and 3, did not have a strong sensitivity, improving only to 28%, while in layer 1 it was improved to 15% from the reference value. in all layers however, was estimated very well (within 2% of the reference). The optimization was most sensitive to the soil moisture distribution, and although some parameters were not e stimated exactly, the soil water retention function (Fig ure 2. 7) and hydraulic conductivity, were well matched. Table 2.1 . Summary of optimization residuals. RMSE (multiplied by a factor of 100 for display purposes) between the reference and op timized estimates of the root distribution, daily soil moisture in Layer 1 (at 0.2 m), and cumulative RWU. Superscripts are used to rank each scenario within that metric (1 being the best fit). Data Density Scenarios DD1 - 2 week data 0.0314 1 0.642 1 0.424 2 DD2 - 3 day data 0.0554 4 1.16 3 0.226 1 DD3 - 4 week data 0.0336 2 1.44 5 0.505 3 Site Characterization Scenarios SC1 - low error in low sens. param. 0.177 6 1.34 4 0.714 5 SC2 - high error in low sens. param. 0.116 5 1.60 6 1.60 6 SC3 - estimate low sens. param. 0.0489 3 0.844 2 0.544 4 32 Figure 2.7. Residuals of optimization. Residuals of root fraction, RWU, van Genuchten - Mualem model retention curve, and soil moisture at 0.2 m for a - d) Data Density, and e - h) Site Characterization scenarios respectively. The low residuals for all scenarios show there was agreement with the reference values. SC1 and SC2 generally had higher residu als. Note the y - axis range differences. 33 We calculated uncertainty in each estimated parameter by calculating the ± 1 standard deviation in the parameters from the iterations with 80% or more improvement in the objective function (Vrugt et al., 2003). Ta ble 2 .2 shows the parameter uncertainty for the eight high sensitivity parameters from Scenario SC3. Most were estimated within ± 1 standard deviation of the reference value. The low sensitivity parameters estimated in SC3, were all estimated within ± 1 st andard deviation of the reference value except and . The uncertainty in the 16 estimated parameters from Scenario SC3 is shown graphically in supplemental Fig ure s A.2. 2 and A.2. 3. Each scenario converged to generally the same result, despite varying the number of estimated parameters and imposed error. Table 2.2 . Uncertainty in parameter estimates. Uncertainty in the eight high sensitivity parameters estimated in Scenario SC3: the petrophysical parameters (PP) from layer 1, and the root p arameter (RP), , and the soil hydraulic parameters (SHP) from layer 1, , , , , We calculated the standard deviation in parameters from the iterations with greater than 80% improvement in objective function. Most high sensitivity parameters were estimated within ±1 of the reference value. Reference SC3 Estimate Mean PP 16.21 16.29 16.26 0.34 - 1.01 - 0.99 - 0.98 0.020 RP 8.14 8.66 8.34 1.29 SHP 1.32 1.37 1.38 0.06 2.70 2.29 2.25 0.37 - 5.79 - 5.87 - 5.80 0.37 0.39 0.40 0.40 0.021 0.066 0.071 0.070 0.004 34 In the DD tests, the model was robust to the quantity and timing of the data used in the objective function. The high temporal frequency data (DD2) only covered one large infiltration and subsequent drying event, but was still able to estimate the SHP quit e well, and provided the best estimate of the RWU. This is likely due to the fact that RWU was dominant during the period covered by the data. However, a month - long dataset too early in the growing season may not capture later RWU since there would be no s ensitivity to the mature root distribution. We found that even with half as much data (DD2), results were still good, indicating that the approach would likely be successful even if data collection opportunities were limited by travel, time, or equipment constraints. While six datasets is very sparse compared to the daily or hourly water content data typically used in hydrological inversions, like Moreno et al. (2015), we found it sufficient for this purpose. The first scenario (DD1), with moderately spac ed data over the entire growing season, otherwise yielded the best fit of the root distribution and soil moisture, indicating that a longer term dataset at frequent intervals is preferable to a shorter or sparser one. In the SC tests, we introduced parame ter error to the inversion. As expected, even low error was detrimental to the estimation of the high sensitivity parameters, particularly for and , with SC1 being amongst the worst across all three metrics (Table 2. 1). Higher fixed error in SC2 led to the highest RMSE in two metrics and the second highest in the third and limited reduction in the objective function. Despite this, the root distribution, RWU, and soil moisture match the reference model quite well, although with generally higher res iduals than the DD scenarios (Fig ure 2. 7). The soil moisture dynamics of the first layer were not highly dependent on the dynamics of the layers below, and are thus still able to be well estimated by the inversion. While increasing the number of parameter s in the SCE algorithm for SC3 slowed the convergence time by a factor of four, we observed that the algorithm did a much better job 35 estimating the high sensitivity parameters than in the scenarios with fewer free parameters but with fixed parameter error (SC1 and SC2). As observed in Hupet et al. (2003), despite the low sensitivity of the SHP in layers two and three, fixing these properties at incorrect values considerably impeded the algorithm from reducing the objective function (Table 2. 1). In contrast, allowing those parameters to be estimated in concert with the high sensitivity parameters resulted in a better fit of the root and soil moisture dynamics (Fig ure 2. 6 and Table 2. 1). The optimization algorithm was robust to the large number of unknowns in this scenario. We note that SC3 also improved the soil moisture estimates in layer 1 from DD2 and DD3, and the root distribution estimate from DD2 (Table 2. 1). This indicates that the timing of the data played a role, as we observed that every other week data with estimated low sensitivity parameters (SC3), performed better than sparser data with perfect low sensitivity parameters (DD2 and DD3). This is likely because this dataset captures the greatest variety of soil moisture regimes. This scenario was al so able to estimate parameters in all three types of soil present at the Kellogg Biological Station field site, indicating that this approach is not as limited by texture or depth as was expected. Contrary to Hupet et al. (2003), we did not observe a limi tation in parameter estimation capacity for the medium - fine textured first layer soil. We imposed a small amount of error in the ER data, but assumed model features such as layer boundaries and basic climate inputs were considered easy to obtain at the fie ld scale, and thus were held fixed at the reference values for this study. Hinnell et al. (2010) concluded that results with low residuals as presented here are likely only obtainable with a physically representative hydrological model, and future work sho uld explore those limitations. Our results show that even with large parameter bounds and 36 conservative starting values for the SHP, the soil moisture and root dynamics could be estimated accurately. 4. Conclusions In this study, we develop and validate th e use of a novel hydrogeophysical inversion algorithm to estimate SHP, PP, and RP simultaneously for a multi - layered soil. Our results indicate that this is a promising approach. Through the accurate estimation of parameters that control root and soil mois ture dynamics within the coupled hydrogeophysical model, the synthetic transient root and soil moisture distributions across a variety of data density and site characterization scenarios were retrieved. This suggests that transient soil moisture processes depend on a unique root distribution, which is critical since it is very difficult to independently measure the effective root distribution in a field setting. While prior research has studied the use of direct soil moisture data to inversely estimate root parameters, we found that data from a simple ER electrode array collected during part of the growing season was capable of the same, even when errors were present in the petrophysical relationship and data was more limited. Relative to approaches that use water content data for the hydrological inversion, the ER data was temporally sparse, yet this did not prove to be a limitation. The high spatial coverage achievable with this non - invasive approach appears to be as useful for capturing root water dynamics . The approach was very successful in the most realistic scenario (SC3), with poor a priori site characterization of the three soil layers and many unknown parameter values. This methodology provides a minimally invasive, and cost effective approach to bet ter understand root water dynamics in a variety of settings. This is a promising result for the study of perennial and/or deep root systems that are particularly difficult to characterize with traditional methods. 37 While a realistic field setting can be exp ected to contain varied topography, soil heterogeneities, and larger root systems, ER data is well - suited to capture such variability. The capabilities of modern computing and robust parameter estimation algorithms make it possible to model increasingly co mplex systems. Future work will expand to 2D and 3D to take full successful validation of the method that identified a unique transient root distribution, we plan to repeat the inversion with ER data from the Kellogg Biological Station Great Lakes Bioenergy Research Center field site for a variety of biofuel crops. Acknowledgments Th e text of this chapter is a reprint of the material as it appears in Vadose Zone Journal in a manuscript first published April 12 2018 (doi:10.2136/vzj2017.08.0154) coauthored by David Hyndman, Anthony Kendall, and Remke Van Dam . This research was funded b y National Science Foundation grants EAR - 0911642 and 1039180, and the USDA NIFA Water CAP grant 2015 - 68007 - 23133. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily refle ct the views of the NSF or USDA. 38 APPENDIX 39 Figure A.2.1. Petrophysical relationship. Petrophysical relationships measured with soil - box experiments, using triplicate samples from three adjacent plots at the study site. The fitted power law relationships (dashed lines) were used in the model. 40 Figure A.2.2. High sensitivity p arameter estimation results scenario SC3. High sensitivity parameter estimation results for scenario SC3.The distribution of estimates from the optimization iterations with greater than 80% improvement in the objective function are plotted along with the reference synthetic value (gold), the best estimate (blue), the starting estimate (red), and + - 1 standard deviation envelope (dashed black). All parameters, except , were improved from their initial estimates. Most parameters were estimated within + - 1 standard deviation of the reference value. The SHP and , and the PP k 1 were close to the + - 1 standard deviation threshold . 41 Figure A.2.3. Low sensitivity p ar ameter estimation results scenario SC3. Low sensitivity parameter estimation results for scenario SC3 . The distribution of estimates from the optimization iterations with greater than 80% improvement in the objective function are plotted along with the re ference synthetic value (gold), the best estimate (blue), the starting estimate (red), and + - 1 standard deviation envelope (dashed black). All parameters were improved from their initial estimates. All parameters except and were estimated within + - 1 standard deviation of the reference value . 42 CHAPTER 3: ROOT WATER UPTAKE OF BIOFUEL CROPS REVEALED BY COUPLED ELECTRICAL RESISTIVITY AND SOIL WATER CONTENT MEASUREMENTS Abstract Biofuel crops, including annuals such as maize, soybean, and canola as well as high - biomass perennial grasses such as miscanthus, are candidates for sustainable alternative energy sources. However, large - scale conversion of croplands to perennial biofuel crops could have substantial impacts on regional water, nutrient, and carbon cycles due to the longer growing seasons and more extensive rooting syste ms compared to most annual crops. Yet due to the limited tools available to non - destructively study the spatiotemporal patterns of root water uptake in - situ at field scales, these differences in crop water use are not well known. Geophysical imaging tools such as electrical resistivity (ER) reveal changes in water content in the soil profile. In this study, we demonstrate the use of a novel coupled hydrogeophysical approach with both time domain reflectometry soil water content and ER measurements to compar e root water uptake and soil properties of an annual crop rotation with the perennial grass miscanthus, grown in close proximity across three growing seasons (2009 - 2011) in the humid temperate climate of southwest Michigan, USA. Our results indicate that E R data are particularly sensitive to a shallow transition from a higher - water content sandy loam to a low - water content sand, which strongly affects the amount of plant - available water. We estimated maximum depths of root water uptake between 0.72 and 1.6 m. The vertical distribution of root water uptake was notably deeper in 2009 relative to 2010 and 2011, likely due to the drought conditions that first year, which we hypothesize encouraged deeper root activity. 43 1. Introduction The amount of land used for cellulosic biofuel cropland is expected to grow in the coming decades and with it, great potential to reduce greenhouse gas emissions, limit nitrogen pollution, increase biodiversity, and meet renewable energy targets (Robertson et al. 2017) . Cellulosic biofuels include high - biomass, productive perennial grasses, woody plants, or non - grain crop residues such as maize stover that can be harvested for bioenergy production as an alternative to petroleum fuel sources. Leading candida tes of perennial grasses include miscanthus ( Miscanthus giganteus ) and switchgrass ( Panicum virgatum ), which have high water use efficiency and lower nitrogen demand relative to maize (Stenjem et al. 2019) . Socio - economic models have projected that over three million hectares of cropland in the United States will need to be converted to energy crop production to meet state and federal standards by 2030 (Oliver and Khanna 2017) . The environmental impacts of such large - scale land use change potentially incl ude altered regional water balances (Abraha et al. 2015; Hickman et al. 2010) although that depends on many factors, including the soil and climatic conditions, crop varieties and cultivars, and agronomic management. There is thus a need to anticipate how an emerging large - scale shift in land use to perennial biofuel crops may affect cropland water balances, along with resulting subsurface and surface water resources. In the face of the projected growth in biofuel crop production and its potential consequen ces for water resources, new tools are needed to measure and simulate how crop water use will respond to agricultural management, soil properties, and climate variability. Modeling these effects is currently limited by gaps in field data for model validati on and the knowledge of the biophysical processes that drive water and nutrient cycling in these nascent agricultural systems (Uhlenbrook 2007) . Long - term field measurements of water balance for various 44 candidate biofuel cropping systems are limited to a h andful of study sites, and do not cover the range of environments where these crops are likely to be grown (Hamilton et al. 2015) . Studies that have quantified the likely effects on the water balance of land use conversion to biofuels in the Midwestern Uni ted States have included watershed to regional - scale modeling scenarios as well as field - scale measurements of evapotranspiration (ET). Modeling results have suggested that replacement of annual with perennial crops, such as the candidate grasses for biofu el production, will likely increase ET, and would thus reduce groundwater recharge and overland runoff (Georgescu, Lobell, and Field 2011; Schilling et al. 2008; Vanloocke, Bernacchi, and Twine 2010) . Side by side comparison of perennial grasses and maize in Illinois found ET to be higher in the perennial grasses (Hickman et al. 2010) . However, Hamilton et al. (2015) and Abraha, Chen, Hamilton, & Robertson, (2020) found that this was not necessarily the case at a southern Michigan site, where ET from the pe rennial grasses was comparable to maize in both wet and dry years based on measurements of soil water uptake and eddy covariance, respectively. They proposed that this discrepancy could be attributed to development of water - limited conditions during most y ears in the well - drained sandy loam soils of the Michigan site. At the plant scale, root function research has mainly been limited to annual plants, and thus much less is known about the spatiotemporal patterns of root water uptake by perennial grass crop s (Ryan et al. 2016) . Understanding the likely implications for water balances of conversion from annual crops or native grasslands to biofuel cropping systems requires knowledge of how root water uptake and transpiration compare between the preexisting ve getation and the new cropping system. This matter is complicated by the variability in root distribution behavior across plant species; a study by Mann, Barney, Kyser, & DiTomaso, (2013) found that switchgrass and miscanthus had notable differences in root growth in response to 45 drought and that the rooting depth of miscanthus increased significantly under irrigated conditions. Since growing biofuel crops on marginal soils is a desirable strategy to avoid competition with food crops on prime farmland, studyi ng differences in root adaptations to dry conditions is of particular importance for accurately modeling ET. Existing methods to characterize spatiotemporal variation in root biomass and distribution are impractical beyond the scale of individual plants be cause measuring roots in - situ is labor - intensive and destructive. Coring or excavation of the root zone can provide an indication of the biomass distribution with depth, and root cores can measure fine root growth (e.g., Sprunger, Culman, Robertson, & Snap p, 2017) . However, traditional extraction and sieving of roots may not capture all fine roots nor can it discriminate between active and inactive roots. Root windows and photography have been used to monitor root growth of several prospective biofuel crops (Mann et al. 2013) , but this approach provides only a one - or two - dimensional sampling of a restricted area. Electrical resistivity (ER) measurements provide a novel and minimally - invasive hydrogeophysical method to study root - water interactions . Resistivity of the soil - plant - water continuum varies temporally with fluctuations in soil water content as well as temperature and fluid conductance. Recent research has taken advantage of t his sensitivity to identify and quantify spatial zones of root water uptake from the decrease in soil water content and concomitant increase in ER over the growing season (Bass, Cardenas, and Befus 2017; Garré et al. 2013; Jayawickreme, Van Dam, and Hyndma n 2008, 2010; Robinson, Slater, and Schäfer 2012; Vanella et al. 2018) the ER data, i.e., to transform the measured apparent resistivity values into laterally - and depth - variable values o f resistivity which can then be related to soil water content. More recently, 46 coupled hydrogeophysical inversion methods using both forward hydrologic and geophysical models have been applied to study shallow subsurface water dynamics including hydrologica l model parameterization (Kuhl et al. 2018; Moreno, Arnon, and Furman 2015; Tran et al. 2016) . Using ER to build a process - based model of the soil water dynamics enables the user to test hypotheses regarding the drivers of the observed changes, as well as the ability to test future and control scenarios. This approach also avoids non - unique and unconstrained solutions inherent to smoothing processes in traditional ER data inversions (Hinnell et al. 2010) . The coupled hydrogeophysical inversion approach al so has advantages over calibrating a hydrologic model only with measurements made by soil water content probes because it provides a 3 - dimensional indication of changes in soil water content. While soil water content measurements accurately represent in - si tu conditions, their coverage is often limited to profiles at single points in a field. Previous research has demonstrated the use of point measurements to estimate root water uptake (e.g., Hamilton et al., 2015) and calibrate hydrological models (e.g., Sc helle, Iden, Fank, & Durner, 2012) , but such approaches are typically not able to characterize within - field heterogeneities. Here, we use for the first time a coupled hydrogeophysical inversion approach to estimate ET and root water dynamics from ER and s oil water content measurements across several non - irrigated biofuel crops in southwestern Michigan, USA, over a three - year period from 2009 - 2011, working in the same experiment where Hamilton et al. (2015) estimated ET from soil water profiles. We develop hydrogeophysical models for two contrasting cropping systems, one of an annual crop rotation and another of the perennial miscanthus grass, to estimate parameters controlling root water uptake and soil petrophysics with a global optimization algorithm. We examine soil water heterogeneity across the field site and consider its implications for 47 cumulative ET for each crop throughout the model period. Although our results are specific to the observed cropping systems, this approach should be broadly applicable for any plant - soil - atmosphere system with sufficient data to build a representative hydrological model. 2. Methods 2.1 Study Site Our study site is the Great Lakes Bioenergy Research Center, an experimental facility on a glacial outwash plain in southwes tern Michigan, USA ( Figure 3. 1 temperate and humid, with mean annual precipitation of 900 mm and mean annual temperature of 9°C. The soils in the region are classified as Typic Hapludalfs, with a bimodal distribution of silt and sand in the upper 0.5 m due to loess deposition that holds more than the underlying coarse sand parent material (Luehmann et al. 2016) . The t hick glaciofluvial deposits on - site are rich in carbonates (calcite, dolomite) below the chemical weathering front at approximately 2 m depth. At this interface, CO 2 driven carbonate dissolution releases Ca 2+ , Mg 2+ , and HCO 3 - ions resulting in high porewater electrical conductivities ranging between 400 - 700 µS/cm and groundwater electrical conductivity of 550 µS/cm (J in et al. 2008) . The water table is at a depth of approximately 17 m. Established in 2008, the Biofuel Cropping System Experiment at the research center maintains five replicate blocks of ten different biofuel crop treatments on individual 28 x 40 m plots ( Figure 3. 1 b). For this study, we investigated two of the treatments in Block 1, comparing root water dynamics of an annual crop rotation (soybean, maize, and canola) (plot G3) and the perennial grass miscanthus (plot G6) ( Figure 3. 1 b) over the 2009, 2010, and 2011 growing seasons. The perennial crop was established via transplanted rhizomes on May 23 and June 4 2008, whereas the annual crops were seeded in the spring of each year on May 22 2009, April 30 48 2010, and May 4 2011. The annual crops were grown us ing conventional agronomic management for the region (Sanford et al. 2016) . All crops were grown without tillage and aboveground biomass was harvested annually. 2.2 Climate, Soil, and Vegetation Data An on - site weather station (station ID: kbs, 42.4081°N, 85.3736°W) , which contributes data to the Michigan Automated Weather Network, measured hourly precipitation, air temperature, ground surface temperature, wind speed, humidity, and solar radiation, allowing calculation of reference grass potential ET using the FAO Penman Monteith method (Allen et al. 1998) . Total growing season (May 1 to Oct 1) precipitation was highly variable across the three model years (358, 568, 462 mm, respectively), with 2009 experiencing a prolonged drought from the third week of June to the second week of August (< 10 mm in 2009 vs . 181 and 175 mm in 2010 and 2011). Soil temperature and water content profiles were collected for each plot. Temperature was measured hourly using data - logging temperature sensors (Thermochron iButton DS1922L) installed vertically at three depths (0.24, 0.64, and 1.25 m) along a 5 - cm diameter buried PVC pipe. Soil water content was recorded hourly using time - domain reflectometry (TDR) probes inserted horizontally at six depths from 0.2 to 1.25 m. For a full site description and details of the TDR installa tion and calibration methods, see Hamilton et al. (2015) . Raw hourly TDR data were filtered using a 24 - hour moving median and moving standard deviation to reduce data noise. Outliers exceeding one standard deviation from the moving median were removed from both the temperature and TDR datasets. 49 Figure 3.1. Site map. a) Great Lakes Basin locator map; the star indicates the Great Lakes Bioenergy Research Center, b) aerial photo of the ten 28x40 m experimental biofuel plots instrumented with electrical resistivity (ER) arrays (white) and time domain reflectometry (TDR) and temperature sensors (yellow) for this study; plot labels G3 (green) and G6 (blue) indicate the annual crop rotation and miscanthus plots respectively, c) view along one of the 0.3 m d eep, 12 m long ER trenches in which 40 graphite electrodes (photographed at left of trench prior to installation) were installed. Soil samples were taken from each plot to characterize texture and density of the soil profiles. Samples were collected with a bucket auger and aggregated into seven sections, from 0 - 0.1 m and then six 0.2 m intervals to a depth of 1.3 m. Percent sand, silt, and clay were measured within each sample using a Malvern Mastersizer laser. Bulk density was estimated from the dry weight of each section. Estimates of field capacity and wilting point for each soil layer were visually inferred from the 2009 TDR data, which experienced significant drying throughout July and early August. Soil physical properties are shown in Table 3. 1 . Leaf area index (LAI) values were measured for each crop in all four cardinal directions at weekly intervals from the period when aboveground growth first appeared until plant senescence in the fall during the 2010 season. For model inputs (described below), LAI data was smoothed by taking the mean of the four directional measurements, then calculating a five - week moving median. A linear interpolation was then used to downscale from weekly to hourly LAI 50 for each crop in 2010. Much less frequent LAI measurements were available for 2009 and 2011, and therefore the same interpolated LAI curve from 2010 was used to approximate 2009 and 2011 hourly LAI values, adjusting the curve forward or backward in time to account for shifts in planting day for th e annual crops and emergence (determined from plant cameras recording daily) in the perennial crop. Depth profiles of root biomass data were available only for the perennial miscanthus, while only total root biomass was measured for the annual crop rotati on. Root biomass for miscanthus was sampled across four depth zones (0 - 0.10, 0.10 - 0.25, 0.25 - 0.50 and 0.50 - 1.0 m), using a 0.06 m diameter core taken at the center, adjacent to, and in between the plants after senescence at the end of each growing season i n each of the five replicate miscanthus plots. Root mass density with depth was calculated by dividing the total root mass from the three locations (center, adjacent and interstitial) by the core length . 2. 3 Electrical Resistivity and Soil Petrophysical Data Prior to the 2009 growing season, graphite rod electrodes were permanently installed within each experimental plot. Within each plot, 40 rods, each 0.08 m long, were placed in a linear array with 0.3 meter spacing in trenches (see Figure 3. 1 b, c), perpendi cular to the crop rows, with the top of the electrodes at ~ 0.3 m depth. Electrodes were buried below the plow zone to avoid disturbance by farming activities. The rods were wired to a takeout box outside each experimental plot. The relative locations and elevations of each electrode were surveyed using a total station prior to trench backfilling. ER data were collected from the permanent arrays with an AGI SuperSting R8 approximately monthly during the growing seasons from May 4 to Sept. 22, 2009, May 3 to Sept. 20, 2010 and May 2 to Sept. 21, 2011. Each survey included 439 measurements of 51 resistance, collected in a dipole - dipole configuration, where two curren t electrodes, and were placed left of two potential electrodes, and at a range of electrode spacings to create a 2D pseudosection, where ( 1 ) and is the measured potential drop across electrodes in Volts, and is the applied current. Sensitivity to greater depths is achieved through increasingly spaced current and potential electrode pairs. The effective depth, , is approximated as 0.2 times the distance between the and electrodes, which in this study ranged from 0.6 to 9.6 m, with current pair spacing ranging from 0.3 to 1.2 m. Measurement errors caused by electrode connection failures during a survey were later removed by a two - step filtering process of values three standard deviations f rom the moving lateral median for each unique effective depth, and then repeated with one standard deviation. Across all surveys, the median data removal with this method was 20 measurements per survey, or around 5%. For the 1D parameter estimation por tion of the study, filtered data were collapsed from 2D to 1D by extracting a median value for each unique effective depth. Soil - box resistivity tests were conducted on soil samples extracted at 0.08 m increments from two 3.6 m long x 0.038 m diameter core s taken with a Geoprobe MT540 at an adjacent plot with comparable soil. Soil box resistivity was measured by drying the soil at 105°C for 24 hours, sieving at 1 mm, and adding distilled water in increments equivalent to volumetric water contents of 0.05, 0 .10, 0.15, 0.20, and 0.30 ( Table 3. 1 ). Soil box resistivity measurements within each soil texture class were averaged to estimate the relationship between ER and water content at the wilting point and field capacity for each soil texture. Only o ne sample was available from the loamy sand soil class, therefore confidence intervals are not reported for that layer. 52 Table 3.1 . Summary of soil properties. Summary of soil properties, including soil texture, % sand, silt and clay, bulk density, field capacity, wilting point, and the corresponding soil box resistivity at 25°C. Note the deeper transition from loamy sand to sand at the annual rotation plot and the increase in resistivity from sandy loam to loamy sand, as well as the distinctly lower resistivity in the carbonate - rich sand versus the carbonate - leached sand. Layer Depth [m] Annual Rotation Layer Depth [m] Miscanthu s Soil Texture % Sand % Silt % C lay Bulk Densit y Field Capaci ty Wiltin g Point Resistivit y at Field Capacity Resistivit y at Wilting Point 0 - 0.12 0 - 0.12 sandy loam 78 15 7 1.13 0.27 0.14 97 ± 27 255 ± 67 0.12 0.45 0.12 0.45 sandy loam 59 27 14 1.35 0.27 0.14 97 ± 27 255 ± 67 0.45 0.75 - loamy sand 83 11 6 1.56 0.18 0.11 621 1025 0.75 0.95 0.45 1.3 sand 90 6 4 1.65 0.12 0.08 1267 ± 177 2186 ± 185 0.95 2.0 1.3 2.0 sand 95 3 2 1.65 0.12 0.08 1267 ± 177 2186 ± 185 2.0 - end 2.0 - end carbonate sand 95 3 2 1.65 0.12 0.08 840 ± 130 1660 ± 130 2. 4 Coupled Hydrogeophysical Inversion To study the root - water dynamics of each crop type, we used a novel coupled hydrogeophysical inversion method with iterative optimization enhanced from that presented in (Kuhl et al. 2018) . For coupled hydrogeophysical inversion, rather than inverting the static ER data to obtain static soil water content distributions, a transient hydrological model is coupled with a forward geophysical model, and model parameters are optimized to match observed geophysical and hydrological data. This inversion method offe rs numerous advantages, chiefly that assumptions often required for traditional inversion such as surface - placed electrodes, topographic invariance, or material lateral homogeneity are relaxed. Indeed, the complexity of the geologic material is limited onl y by the complexity of the forward hydrologic and geophysical models. 53 The specific steps for the coupled inversion applied here are: 1) A transient hydrological model is run across all observation times, including spin - up periods needed for model equilibra tion; 2) Static 1D soil water content profiles are extracted from the transient hydrological model at times corresponding with the date and time of the ER surveys; 3) A layer - specific petrophysical relationship is applied to convert soil water content to f irst reference - temperature resistivity, then corrected for the in - situ temperature; 4) The 2D electrical potential field is then calculated in a forward ER model to compare to the measured resistances from the ER surve y and lastly; 5) The hydrologic an d petrophysical model parameters are then optimized to fit observed TDR and ER data using a global optimization algorithm. In this study we introduce a new two - step optimization routine that allows for the model to invert lateral soil type transitions usin g the ER data. For this procedure, we first optimize the hydrologic, root, and petrophysical model parameters using lateral averages of both modeled and observed ER data. For this first step, we use a single 1D hydrologic model with a fixed soil texture in terface depth inferred from soil texture data. Second, using those optimized parameters, we then run multiple 1D hydrologic models, each applied to an interval along the ER transect. Soil layer transitions are optimized for each 1D hydrologic model to acco unt for the along - transect lateral variability in soil properties. 2.4.1 Forward Hydrological Model We built a HYDRUS - 1D soil hydrological model for each crop to simulate the 1D hydrological fluxes of each plot for a three - year mod el period from Nov 1 2008 to Nov 1 2011. HYDRUS - flow with many options for simulating additional subsurface transport processes. For this study we ran HYDRUS with root growth, root water uptake (transpiration), evaporation, heat and CO 2 54 transport, and snow hydrology, as well as the Unsatchem module to model soil water ion concentrations. The model was vertically discretized into a 0.04 m grid down to 17 m depth, with an hourly time step. Upper and lower boundary conditions were specified for each of the water, heat, and CO 2 transport modules. Hourly atmospheric fluxes from the weather station, including precipitation and potential ET (partitioning method described below) with a surface runoff threshold of 2 mm, and free drainage, were used as the upper and lower hydrologic boundary conditions, respectively. Canopy storage was calculated as 0.15% of LAI (Dai et al. 2003; Dickinson et al. 1991) and was used to reduce incoming precipitation until maximum storage capacity was exceeded during individual rain events (defined as consecutive hours with measured precipitation). Hourly soil surface temperatures were used as the upper boundary condition of the heat transport model and a heat flux boundary condition was set for the bottom thermal boundary. Atmospheric concentrations of CO 2 (0.00033 m 3 /m 3 ) were used as the upper boundary condition in the CO 2 transport model, while a zero gradient conditio n was used for the bottom boundary. Potential reference grass ET, from the weather station was adjusted to a crop specific potential ET, using an estimated crop coefficient [ - ] that was temporally scaled with LAI (e q. 4), and partitioned into potential evaporation, , and potential transpiration, , using the interpolated hourly LAI curves for each crop and a shape canopy factor, (eqs. 5 - 7). During non - growing periods, was set equal to 0.4 (FAO 2002). ( 2 ) 55 ( 3 ) ( 4 ) ( 5 ) Figure 3.2. Hydrological model inputs. Select model inputs for the 2009 - 2011 model period: a) measured precipitation and b) air temperature from the weather station, c) interpolated leaf area index (LAI) data for the annual rotation (green) and miscanthus (blue) plots, and d) potential evapotr anspiration (PET) for both crops calculated from the reference grass PET (grey line) using crop coefficient values (Equations 4 - 7). The 2009 - 2011 growing seasons in the annual rotation were planted in canola, maize, and soybean, respectively. Model structu ral parameters were then specified as follows: for the soil hydraulic model, we used the Van Genuchten - Mualem (VGM) equation. Parameters of the VGM model for each layer were derived from the Rosetta database (Schaap, Leij, and Genuchten 2001) using the gra in size analysis and bulk density measurements for each plot ( Table 3. 1 ). Initial soil layer transition depths (distinct from the HYDRUS computational layers) were set as half the distance between the depths of the observed changes in grain size unless the TDR data suggested a layer boundary 56 should be adjusted (transition depths shown in Table 3. 1 ). Heat transport parameters were calibrated with soil temperature data from previous work at the study site provided in Kuhl et al. (2018) . Parameters controllin g CO 2 transport and production were also taken from the database in HYDRUS but were not modified. We used a separate HYDRUS 1 - D executable including a new root growth module (Hartmann et al. 2018) that allows increased flexibility of the modeled shape and depth of the root distribution through time. We selected the Vrugt (Vrugt, Hopmans, and Simunek 2001) equation to control the shape of the root density distribution, : ( 6 ) where is the root density at depth , is the maximum rooting depth, and and are fitting parameters that control the shape of exponential decay in the root distribution; and parameters of the equation for sinusoidal root growth (Hao, Zhang, and Kravchenko 2005) : ( 7 ) where is the potential increase in rooting depth at time , and are the time to planting and time to maturity of the plant respectively, and [m] is the potential rooting depth from time The parameter can be altered to change how quickly the plant reaches its maximum rooting depth, effectively controlling the rate of root growth. Root parameters for all crops were initialized with , and . We assumed the maximum ro ot density was near the surface and therefore all parameters except were later optimized. 57 To model major ion chemistry, ionic concentrations of the seven major ions, Ca 2+ , Mg 2+ , Na + , K + , HCO 3 - , SO 4 2 - , and Cl - of the incoming precipitation was specifi ed with the summer seasonal average values (12E - 3, 3.0E - 3, 0.91E - 3, 0.41E - 3, 13E - 3, 3.4E - 3 meq/L, respectively) recorded by the National Atmospheric Deposition Program/National Trends Network station MI26 located at the study site. Calcite (CaHCO 3 + ) and do lomite (CaMg(CO 3 ) 2 ) below 2 m were specified as a precipitate in the un - leached carbonate zone below 2 m in concentrations of 500 meq/kg each. The coefficient of molecular diffusion in free water and the longitudinal dispersion rates were set at 0.08 and 0 cm - 1 , respectively, based on values from Example 2 in the Unsatchem Manual . The electrical conductivity of the porewater, [Sm], was calculated by summing the seven major ions, output at each node McNeal, Oster, & Hatcher, (1970) for mixed salt solutions: ( 8 ) where is a linear coefficient and is a linear offset for each ion, , scaled by the total cation concentration ( ). Separately, we used the PHREEQC software to determine that ion pairing was not an important influence on solution resistivity in the relatively dilute soil porewater solutions above 2 m. An additional three year ramp up period beginning on Nov 1 2005 (repeating climate and LAI conditions from 2009 - 2011, data were unavailable prior to 2009) was used for the Unsatchem module runs to stabilize the pore water conductivity which was very sensitive to the initial ion composition specified as that of the precipitation. Because the root growth and Unsatchem modules could not run in co ncert, the final root distribution from the root growth - enabled model was used with the Unsatchem module; we assumed that difference of 58 root water uptake between the standard static HYDRUS root profile and the dynamic root growth module had a negligible ef fect on the ion content of the remaining soil water . 2.4.2 Subsurface Electrical Resistivity Subsurface ER is the combined resistivity of the soil - water matrix and is here calculated in a two - step process. First, ER is calculated at a 25°C reference tem perature, then adjusted to in - situ temperature. The ER of the soil water matrix, , at each node, , and each survey time, (Waxman and Smits 1968): ( 9 ) where is the static porosity, is grain surface conductivity, is dynamic modeled water content, is modeled conductivity of the pore water, and and are dimensionless cementation and saturation exponents. Initial estimates of and were made by assuming and fitting eq. 9 to the soil box - derived petrophysical curve given the porosity and modeled pore water content for each soil layer. These parameter estimates were later updated via optimization with the ER data. Resistivity at in - situ soil temperatures was then calculated using modeled hourly soil temperature, , at each depth, , and time . We applied a temperature correction to proposed by Hayley, Bentley, Gharibi, & Nightingale, (2007) , valid for temperature ranges between 0 - 25°C: ( 10 ) 59 2.4.3 Forward Geophysica l Model The forward geophysical model was built using the Boundless Electrical Resistivity Tomography (BERT) model for Python (Rücker, Günther, and Spitzer 2006) . The 2D model domain was 20 m wide by 20 m deep with a fine - resolution triangular mesh (element area = 0.01 m) in the upper 4 m, a coarser mesh (element area = 2 m) to the water table at 17 m, and a very coarse mesh (element area = 5 m) between 17 and 20 m. Additionally, the mesh was refined around the location of each buried electrode. The model domain was extended 4 m to the west and east of the 12 - m long ER array, and 18 m deeper than the maximum effective depth to avoid boundary edge effects. The surveyed electrode locations, surface topography, and current/potential elect rode pairs for each plot were input to BERT to replicate the dipole - dipole ER surveys ( Figure 3. 3 ). For every ER survey date for each plot, the simulated 1D vertical temperature - corrected resistivity distribution was imported to BERT assuming lateral homog eneity in both soil and root distribution to create a pseudo - 2D resistivity distribution. 60 Figure 3.3. Geophysical model mesh. The upper 2 m, and inner 12 m of the BERT model domain and triangular mesh with surveyed surface topography from west to eas t for the annual crop rotation and miscanthus plots. Note the mesh refinement around the 40 electrodes buried approximately 0.30 m below the surface. Pockets of fine - resolution mesh near the center are due to the automatic mesh generation process. 2.4.4 Op timization Finally, the coupled hydrogeophysical inversion was accomplished via automated parameter optimization, adjusting the parameters of hydrologic, and petrophysical models to match both TDR and ER data. For each experimental plot, our two - step optim ization procedure first optimizes hydrologic, root growth, and petrophysical parameters for a single 1 - D hydrologic model, which is then used to compute laterally - homogeneous resistivity for the 2 - D ER model. Then, as a second optimization step, multiple 1 - D hydrologic and petrophysical models are built, and the soil layer transition depth is optimized for each model. These models each apply to a 1.2 - m wide interval along the ER transect. Optimization of the first step uses a dual - component objective funct ion, , a weighted root - mean square of both resistivity and soil water content residuals of the form: 61 ( 11 ) where is the index of ea ch observation data point in space and time, is the total number of measurements, is the water content, is the apparent resistivity, and is the standard deviation of the dataset. In the first step, was minimized by updating the 16 parameters of interest using a Shuffled Complex Evolution global optimization algorithm (Duan, Sorooshian, and Gupta 1992) . These parameters included: 1) root growth and distribution and for each of the three growing seasons; 2) the petrophysical parameters for the loam, sandy loam, and loamy sand layers, and 3) for the loam layer only, assuming surface conductivity is negligible in the coarser grained textures. For Equation 11, we first averaged the ER data for each unique effective depth to smooth the along - transect variability in apparent resistivity. This was done for both the modeled and observed resistivity. Soil water content values were taken from the 1D hydrologic model. Th e root mean squared errors of the two datasets were weighted by their standard deviations. For soil water content, we used the mean - differenced water content (e.g. ) to minimize the effect of inter - layer variability in soil properties th at our simplified layer assumptions did not capture. Parameter estimates were determined by the mean value of the lowest 10 th percentile of values from the global optimization. 95% confidence intervals were calculated as: ( 12 ) where , is the mean, and , is the standard deviation of each parameter distribution from iterations within the lowest 10 th percentile of values from the global optimization. For our second optimization step, the hydrologic, root growth, and petrophysical parameters were fixed from the first optimization, and only the depth of the transition between 62 the sandy loam and loamy sand soil layers along the transect was estimated. Thi s transition has a strong effect on the apparent resistivity, allowing for lateral sensitivity to the soil properties. Each of the two experimental plots was divided into 10 1.2 - m wide zones, and the depth of the transition between these two layers in each zone was optimized with the full (non - averaged) 2D ER dataset in a global optimization by minimizing the resistivity component of only the resistivity component of . The TDR data was sensitive to any lateral heterogeneity and therefore was not utilized in the objective function for this portion of the study. Unique root growth parameters for each zone were not estimated with the 2D data. 3. Results and Discussion A summary of the results of the parameter estimation via global optimization is presented in Table 3. 2 between water content and electrical resistivity. For both plots, distinct distributions for the best and worst iterations are clearly observed for surfa ce electrical conductivity ( ) and the parameter of the sandy loam soil layer. For the remaining parameters, the worst iterations are distributed across the full parameter search space, indicating no single parameter estimate is significantly detrime ntal to the objective function. Among the best iterations, a normal distribution appears for most parameters, indicating convergence towards that value as the best estimate to explain the data. The TDR component of the objective function ranged between 0.5 and 1.1 while the ER component ranged between 0.2 and 2.2. There was an overall weak positive correlation between the two components (R = 0.29 and 0.35) for the annual rotation and miscanthus plots, respectively) due to the strong influence of the petroph ysics on the ER but not the TDR measurements. Amongst the root growth and water uptake parameters, the crop coefficient , which controls the magnitude of the potential ET, was the most sensitive, with a narrower distribution 63 and confidence interval th an the other parameters. The days to maturity parameter , which controls how quickly the roots reach , was the least sensitive for all years in both plots, with wide confidence intervals and estimates that span the entire search space. Notably, 200 9, which was the first year after establishment for the miscanthus crop and the year that experienced the most drought - like conditions, shows significant differences in the distribution of the maximum rooting depth, , the shape parameter , as well as in , compared to 2010 and 2011. Note that increasing and decreasing together as seen in miscanthus in 2009 has an additive effect, making the ratio of smaller, shifting root fraction downward relative to the other years (eq. 2). I n the annual rotation plot by contrast, the co - variation in and results in almost no change from year to year in the ratio of . Crop LAI for 2009 and 2011 were based on the most complete year of measurements (2010) and therefore the estim ated maximum LAI in 2009 cannot be validated with measurements. Given the importance of LAI in the calculation of potential ET (eqs. 4 - 7), t he low miscanthus parameter estimate for the 2009 growing season relative to the other years could be the model adapting to an overestimation of LAI (and thus potential T) for that period. It also could be reflective of the juvenile nature of the miscanthus crop in its first year post - establishment. This is an encouraging result as it may demonstrate a mechanism to overcome errors or gaps in LAI data. This large effect of LAI on actual ET is most notable in 2009 where canola grew quickly and was harvested early in the season, as well as in the 2011 growing season where the planting of soybean occurred much later tha n the emergence of the established miscanthus. 64 Table 3.2 . Parameter optimization results. Optimization results with initial and final estimates, with 95% confidence intervals, for each plot and each model year . Soil horizons are abbreviated as follows: sandy loam (SL), loamy sand (LS), sand (S). Note that in the miscanthus plot, the parameter of the loamy sand (LS) layer was not optimized due the absence of a substantial loamy sand layer in that plot. Root Parameters Petrophysical Parameters SL LS S SL Initial 1.2 1.2 1.2 1.0 1.0 1.0 1.0 1.0 1.0 100 100 100 - 0.80 1.5 2.6 0.01 Annua l Crop 1.6 ± 0.034 1.0 ± 0.03 1 1.1 ± 0.02 8 5.4 ± 0.34 5.3 ± 0.34 4.6 ± 0.35 2.0 ± 0.08 0 1.7 ± 0.11 1.7 ± 0.09 0 88 ± 6.4 123 ± 6.0 113 ± 5.9 - 0.22 ± 0.064 0.05 9 ± 0.12 2.5 ± 0.03 4 0.012 ± 4E - 4 Misca n - thus 0.95 ± 0.017 1.1 ± 0.02 1 1.3 ± 0.02 5 4.7 ± 0.26 6.4 ± 0.30 6.7 ± 0.35 2.2 ± 0.06 5 1.3 ± 0.07 1 1.2 ± 0.07 4 134 ± 5.1 130 ± 5.3 112 ± 5.7 - 0.063 ± 0.051 - 2.4 ± 0.03 0 0.012 ± 3E - 4 65 Figure 3.4. Optimization results. Parameter distributions for the 90 th percentile (grey) and best 10 th percentile (blue) iterations of the global optimization for the miscanthus plot. Black dashes represent the 95% confidence intervals, while the mean values are reported at the top of each distribution. Greater parameter sen sitivity is indicated by a narrower distribution between dashed lines. The petrophysical parameters and have the greatest influence on the objective function , with a clear bimodal distribution in the worst performing iterations distinct from a mo re normal distribution for the highest performing iterations. The other parameters show more random distributions spanning the entire parameter search space for the worst performing iterations. Modeled outputs of cumulative ET, root depth through time, and water content through time for both cropping systems are shown in Figure 3. 5 . Modeled cumulative ET for each growing season (May 1 Oct 1) in the crop rotation treatment was 304, 395, and 324 mm for canola, maize, and soybean, respectively, and in the perennial miscanthus it was 351, 430, and 410 mm in 2009, 2010, and 2011 ( Figure 3. 5 a). Higher ET in the miscanthus plot was driven by higher LAI and crop coefficient values and a prolonged season after the annuals were harvested, 66 supporting some previous research findings that miscanthus ET rates will exceed those of annual crops, if soil water is available. Cumulative ET for the miscanthus was in good agreement with that calculated by soil water content loss for the same plots by Hamilton et al. (2015), where they found between 2010 and 2013 average miscanthus ET was 458 ±31. mm. Higher ET in 2010 versus 2011 due to greater water availability was also observed by Hamilton et al. (2015). However, o ur model output for maize during the growing season with greatest ET (2010: 395 mm) was 15% lower than the TDR - based estimate from Hamilton et al. (2015) of 467 mm for a maize - only treatment plot. 67 Figure 3.5. Model output and validation. Selected model outputs from the 2009 - 2011 growing seasons (May 1 to Oct 1). a) Cumulative growing season evapotranspiration (ET) for the annual rotation (green) and the miscanthus (blue) plots, b) maximum depth of root fraction (eq. 2) exceeding 1% throug h time, c) average daily modeled volumetric water content (VWC, solid line) and time domain reflectometry measurements (black markers) at select depths of 0.2 m and 1.25 m through time, and d) as a 1:1 plot of modeled (Mod.) vs. observed (Obs.) at all six depths for the annual rotation plot and for miscanthus in e - f). Note the ER survey sampling dates are indicated with the vertical gold lines on c and e. The 2009 - 2011 growing seasons in the annual rotation were planted in canola, maize, and soybean, respec tively. Model parameterization of root depth through time ( Figure 3. 5 b) yielded reasonable estimates for both crops, with root fraction decreasing to <1% of the total root density at depths of 1.32, 1.20, and 1.24 m for canola, maize, and soybean res pectively, and 1.64, 0.80 and 0.72 for miscanthus for each model year. Estimates for the annual rotation fall within the ranges of temperate agricultural crops summarized in Fan, Mcconkey, Wang, & Janzen (2016) . Note that the root depth shown in Figure 3. 5 b is shallower than the parameterized values shown in 68 Table 3. 2 , which represent the depth where root fraction becomes 0. The exponential decay of root density (eq. 2) produces a long tail of nominal root fraction values which we chose to truncate at 1%. We hypothesize that the deeper root depth for both crops in 2009 is likely a response to the drought conditions, which is supported by the limited measurements of root biomass collected at the end of the growing season using soil cores at the center an d adjacent to the plant, which show a deeper distribution of root mass in 2009 than in either 2010 or 2011 ( Figure 3. 6 b). The model generally overestimates the root fraction (normalized by the maximum model and measurement values) with depth relative to th e data ( Figure 3. 5 b), however this could be explained by the decay of deeper fine roots between the growing season and the sampling period after senescence in December. No root samples were collected from below 1 m, so we cannot compare maximum rooting dep th, however given the low root density in the 0.5 - 1.0 m core length ( Figure 3. 6 a), we infer that relatively few roots grew to depths below 1 m at this location. 69 Figure 3.6. Modeled root distribution. Depth profiles of a) average log observed root mass per volume with 95% confidence interval error bars (n = 5) from four depth cores (shaded region indicates length of each core) sampled at the end of each growing season for miscanthus; b) observed root m ass per volume normalized by the maximum measurement for each of the three sample dates (open circles) and modeled normalized root fraction distribution (dashed lines) at the end of the growing season for each year. Note that the 2009 sample was not collec ted until the following spring, and the downward shift in root distribution could be attributed to root decomposition in the upper 0.1 m after the end of the growing season. 2009 was also the first year after establishment from rhizomes, which could be exp ected to have less root mass than subsequent years. Using ROSETTA parameters to model soil water retention, the absolute value of the modeled water content and the rate of drainage show good agreement with the observations in the sandy loam (0.2 m depth) a nd sandy (1.25 m depth) soils for both crops ( Figure 5 c - e). The largest non - winter discrepancy occurs in the Fall of 2010 in the miscanthus ( Figure 5 e), where shallow water content does not increase with precipitation events as seen in the annual rotation nor in the model of either crop. This may be due to substantial crop residue at the surface that limits infiltration. Aside from this, agreement in the other four depths is also good in both field 70 plots, as observed in the 1:1 plots of all observation dept hs from May 1 to Oct 1 of each model year ( Figure 3. 5 d, f). In addition to model performance through time, we show for the 2010 growing season, the 1D hydrogeophysical model outputs and observations, where available, of root distribution, water content, t emperature, and porewater electrical conductivity (inverse of resistivity), all of which are used to calculate the temperature - corrected resistivity with depth via equations 9 and 10 for each ER survey date ( Figure 3. 7 ). Overall, the model reasonably simul ates the dominant trend of decreasing water content and increasing temperature throughout the growing season and the resultant resistivity changes. Modeled pore water ion content, and hence, varies with depth and time but over a range of only about 5% (Fi gure 3. 7d), resulting in smaller percent effects on the bulk soil resistivity (eq. 9). Temporal changes in soil water content and temperature, which tend to be correlated, are the dominant factors driving observed resistivity changes, with the clearest inc rease in resistivity occurring between the wetter early season surveys (blues and green) and the drier growing season surveys (yellow and pink) ( Figure 3. 7 e). Amongst the early surveys, the most resistive in both the model and the observations is the earli est survey on May 3 2010, due to the much colder temperatures ( Figure 3. 7 c), despite a slightly higher pore water conductivity ( Figure 3. 7 d), which has the opposite effect. 71 Figure 3.7. 1D model and observations. 1D outputs of the 2010 growing seaso n for the miscanthus plot, showing the model results (solid line) and observations (circular markers) at monthly intervals from May to September, for a) normalized root density distribution, b) volumetric water content (VWC), c) temperature, d) pore water electrical conductivity (EC), e) log resistivity, and f) apparent resistivity in the upper 2 m of the model domain. The depth of the apparent resistivity curve is the effective depth (a function of electrode spacing) and is only an approximation. The stron g contrast in volumetric water content and surface electrical conductivity due to the presence of clays between the sandy loam and loamy sand soil interface at ~ 0.5 m drives the shape of the modeled resistivity and apparent resistivity curves. To summariz e how well the model predicts the changes in resistivity due to hydrogeophysical processes across the two plots and three growing seasons, we calculated the change in apparent resistivity for each unique geometry and compared the model to the observations ( Figure 3. 8 ). Generally, the modeled magnitude and direction of change agree approximately with the observations for both plots, particularly later in the growing season of each year during extended drying periods, which the model predicts quite well ( Figu re 3. 8 ). We also note that the change in resistivity is generally similar across the two crops, although there 72 are notable exceptions. In 2009, there is a stark difference between the June 1 and June 29, 2009 surveys, where the annual rotation plot data show an increase in resistivity while miscanthus shows a decrease ( Figure 3. 8 a). This discrepancy is not seen in the model, which predicts a decrease in resistivity for both plots due to a large precipitation event on June 19 2009. Similarly, the chan ge between the June 28 and July 26 2010 surveys is undermodeled, with data for both plots showing an increase in resistivity while the model predicts a small decrease ( Figure 3. 8 b). This is likely due to a large precipitation event occurring on July 22, 20 10 that caused a greater pulse of infiltrating water in the hydrological model than is observed in the data. Figure 3.8. Modeled vs . observed ER. Summary of modeled (Mod) versus observed (Obs) changes in the 1D apparent resistivity between consecutive surveys in each model year for the annual rotation (open circles) and the miscanthus (stars) plots. Positive delta values indicate an increase in resistivity with time due to drier conditions during that period, while negative values indicate a decrease i n resistivity with time due to wetter conditions. The amount of change scales with the absolute value of apparent resistivity, which increases with depth. Using the estimated 1D hydrogeophysical models for each crop to estimate the spatial heterogeneity in the depth of the terminus of the sandy loam soil layer resulted in good agreement in the 2D ER data ( Figure 3. 9 ). The trend in depths of this layer interface roughly follows the surface topography ( Figure 3. 3 ) in both plots, indicating a relatively flat laying layer 73 that is deepest at the highest elevation to the west and becomes shallower with decreasing surface elevation to the east, with a slight rise in both plots west of the center. Depth estimates ranged from 0.50 m in the first zone to 0.38 in the tenth zone with a mean of 0.48 m for the annual rotation plot and from 0.52 m in the fourth zone to 0.38 in the eighth zone with a mean of 0.45 for the miscanthus plot ( Table 3 .3 ). Table 3.3 . Soil interface optimization results. Estimated transition depth to coarser soil texture with lower water holding capacity in each of the 10 zones from 0 to 12 m (west to east; Figure 3.3). West to East distance (m) 0 - 1.2 1.2 - 2.4 2.4 - 3.6 3.6 - 4.8 4.8 - 6.0 6.0 - 7.2 7.2 - 8.4 8.4 - 9.6 9.6 - 0.8 10.8 - 12 Av g. Annual Crop Rotatio n 0.50 ± 0.005 0.49 ± 0.00 6 0.50 ± 0.00 5 0.56 ± 0.00 4 0.57 ± 0.00 4 0.48 ± 0.00 4 0.50 ± 0.00 6 0.47 ± 0.00 5 0.39 ± 0.00 4 0.38 ± 0.00 5 0.48 Miscan - thus 0. 47 ± 0.008 0. 44 ± 0.00 4 0. 46 ± 0.00 5 0. 52 ± 0.00 6 0. 47 ± 0.00 4 0. 46 ± 0.00 5 0. 42 ± 0.00 6 0.38 ± 0.00 5 0. 40 ± 0.00 4 0. 48 ± 0.01 0. 45 For both crops, the modeled ER was generally overestimated closer to the surface where resistivity was low, and underestimated at measurements with the larger electrode spacings, where resistivity was high ( Figure 3. 9 e, f), suggesting that our simplified assumption of three sets of petrophysical parameters for the three primary soil textures may not be sufficient to capture the verti cal heterogeneity. It should also be noted that the optimized petrophysical curve for the sandy loam soil was considerably lower than the soil box resistivity model; modeled wilting point resistivity values did not exceed 100 ohm - m versus the 255 ohm - m obs erved in the soil box experiments ( Table 3. 1 ). 74 Figure 3.9. 2D pseudo - sections of ER model and data. 2D pseudo - sections (shown here without topographic correction) of a) observed and b) modeled apparent resistivity from the annual rotation plot and c) observed and d) modeled apparent resistivity from the miscanthus plot on Aug 24, 2009 after spatial opti mization of the transition depth at the bottom of the sandy loam soil layer. White gaps are data outliers that were filtered and neither modeled nor included in the objective function calculation. Modeled versus observed ER for all 17 survey days for the e ) annual rotation, and f) miscanthus crops. Note that overall ER is higher for the miscanthus plot than the annual rotation, due to the absence of an intermediate loamy sand soil and a shallower transition to sand in the miscanthus plot. Conversely, the modeled resistivity ranges for the sand layer were roughly three times greater than the soil box resistivity. Above the carbonate zone that begins at 2 m depth, the only ionic input to the pore water conductivity model came from dilute precipitation, while other 75 sources of ionic input (particularly nitrate leaching in the fertilized maize and mineral weathering of aluminosilicates) that could modestly increase pore water conductivity, were not included. Error in the modeled pore water conductivity linearly propagates to the resistivity estimates and could significantly alter the estimated petrophysical parameters, although the estimated parameter can accommodate for some of this error. Additionally, while a greater number of zones may have resulted in bet ter agreement between the model and observed ER, for the purposes of this study we limited the number of zones to 10 for efficiency. To analyze the sensitivity of seasonal ET to realistic variations in depth of the interface between high and low field cap acity soils, we compared the cumulative ET in each zone to the 1D model with a fixed interface depth at the 10 - zone average of 0.48 and 0.45 m for the two cropping systems, respectively. There is a clear correlation between drought conditions and ET discre pancy with the largest difference in 2009 and the smallest in the beginning of 2010, when there was consistent rainfall throughout the growing season ( Figure 3. 10 ). When high ET demand draws down soil water content (i.e., during dry periods), the presence of a deeper sandy loam soil with more plant available water becomes more important. The greatest modeled changes in ET ( - 0.98% in 2009 (canola) and - 0.94% for miscanthus in 2011) relative to the ET at the average depth were during dry conditions where the soil interface is shallowest. Where the soil interface is deeper than the baseline, modeled ET is greater. During the 2009 growing season, the cumulative seasonal baseline ET for the canola crop was approximately 304 mm, thus this soil heterogeneity equate s to about 3 mm of difference in available soil water across the transect. The effect is less pronounced for miscanthus due to there being less variability in interface depth ( Figure 3. 10 ). However, under the drier conditions in 2011, the miscanthus crop 76 s hows a much larger decrease in ET than the soybean crop, likely due to the roots being more shallow ( Figure 3. 5 ), which limits access to deeper water reserves. Figure 3.10. Sensitivity of ET to soil interface. Percent (%) difference in cumulative wat er year (Nov 1 to Oct 31, demarcated with vertical lines) ET (mm) relative to the baseline due to the estimated variability in the depth of the sandy loam interface predicted from the ER data for the a) annual rotation and b) miscanthus crop treatments. Po sitive % change means more ET relative to the baseline, while a negative % change indicates a decrease. Line colors are cooler when the interface is deeper and warmer when it is shallower. While the sensitivity of ET is generally small (<1%), it is most pr onounced where the difference in interface depth is largest and ET demand is highest, as expected. Note that only four and three lines (respectively) are visible due to multiple zones having the same or similar estimated interface depths. 4. Conclusion s Here, we demonstrate the use of a novel coupled hydrogeophysical inversion approach to estimate root growth and root distribution properties, along with heterogeneous soil properties below non - irrigated biofuel crops in a field setting. With mo nthly ER surveys over three growing seasons and continuous hourly TDR data, we estimated the root parameters of multiple hydrological models to describe the root water dynamics of a canola - maize - soybean annual crop rotation and a perennial miscanthus grass crop. Using HYDRUS - 1D, we built for the first time to our knowledge, a forward process - based model that simulates changes in root distribution, water content, pore water conductivity, and temperature to predict changes in ER through time, using a petrophy sical relationship with estimated parameters. We found good agreement between 77 modeled and observed apparent resistivity as well as the change in resistivity through time across all 17 ER surveys in both crop treatments. While the simulated water content re asonably matched the TDR data throughout the model period at all six observation depths, the magnitudes of pulses of infiltrating water were generally overestimated by the model, resulting in greater discrepancies between the ER model and observations near those events. Reasonable estimates of ET and root depth for all crops were achieved, with the method predicting deeper roots in 2009 during drought conditions. Simulated miscanthus ET agreed well with prior findings by Hamilton et al. (2015) although sim ulated maize ET in 2010 was about 15% less than those prior data - driven estimates. The perennial cropping system had higher ET than the annual crop rotation, driven by a longer growing season, higher LAI, and a greater crop coefficient. Deeper root distrib ution did not directly lead to greater ET across crops, as maximum root depths were found to be deeper for the maize and soybean than miscanthus in 2010 and 2011, respectively. Although in - situ root distribution data for validation was only available for m iscanthus, the modeled root distribution estimates generally agreed with the data, with root mass concentrated in the upper 0.5 m. This also agrees with the shallow root distribution for rainfed miscanthus found by Mann et al. (2013) . Using the calibrated hydrological models and the 2D ER data to estimate the field scale heterogeneity in soils, we resolved the lateral variation in thickness of the upper sandy loam soil layer from 0.38 to 0.57 m. We found the relatively high water avai lability in this top soil layer to have a moderate influence on the cumulative growing season ET, particularly in 2009, when the driest conditions were observed. ER data are very sensitive to changes in soil water content and large transects or even 3D re gions can be surveyed with minimal soil disturbance. This study demonstrates that coupling ER data with TDR soil water content data and a process - based hydrological model allows 78 interpolation temporally between ER surveys and extrapolation beyond the point scale of soil water content sensors. This approach holds great promise as a tool to better understand root water dynamics of crops and other vegetation, improving our ability to make inferences about the implications of land use change on regional water b alances and informing cropping decisions. As we look towards a clean energy future, understanding and minimizing the environmental impacts of large - scale conversion of land to biofuel crops will remain a critical research question. Acknowledgments This cha pter was coauthored by Anthony Kendall, Remke Van Dam, Ste phen Hamilton, and David Hyndman. A version of the text has been submitted to Vadose Zone Journal for review. This material is based upon work supported in part by the Great Lakes Bioenergy Research Center, U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research under Award Numbers DE - SC0018409 and DE - FC02 - 07ER64494; the National Science Foundation Long - term Ecological Research Program (DEB 1832042) at the Kellog g Biological Station, and by Michigan State University AgBioResearch; and by the US 2015 - 68007 - 23133 . This work was made possible by the efforts and support of many. In particular we would like to thank Kaya Diker for his extensive work on the electrical resistivity array installation and data collection; Kevin Kahmark and Mir Zaman Hussain for assistance with collecting the soil cores; Randall Schaetzel for the grain size analysis and complimentary use of his soil lab; and Amanda Liddle for her effort preparing soil samples and conducting the soil box resistivity tests. 79 CHAPTER 4: A NEW COUPLED HYDROGEOPHYSICAL INVERSION TO ESTIMATE COARSE ROOT FRACTION FROM ERT DATA Abstract Root mass c ontributes significantly to global carbon storage and respiration, yet quantifying it remains a challenge. Electrical resistivity methods have promise for quantifying the spatial distribution of root mass in field settings due to the resistive nature of ce llulosic materials that can contrast with less resistive soil. Here, we show the use of a fully coupled hydrogeophysical model that incorporates thermal, geochemical, and hydrological data across a successional forest and grassland environment to estimate coarse woody root mass fraction in the shallow subsurface. We found a statistically significant correlation between our estimated spatially variable coarse root mass fraction and above ground biomass. Our findings suggest that this minimally invasive appro ach has the potential to monitor soil moisture dynamics and identify the spatial heterogeneity of coarse root mass and soil hydraulic properties at the field scale. 80 1. Introduction In temperate, deciduous forests, roots often comprise 20 to 30% of tre e biomass (Mokany, Raison, and Prokushkin 2006; Sinacore et al. 2017) , which represents a large portion of global carbon stocks. Tree root structure is made up of both fine (<1 mm diameter) and coarse roots, each with distinct roles in carbon and water cycling, however, coarse roots account for approximately 90% of total ro ot biomass (Fahey et al. 1988) . Despite their significance for pressing environmental concerns related to the impacts of climate change and deforestation, studying tree roots in - situ and at the field scale has many challenges. Their inaccessibility is a co nsiderable barrier to understanding feedbacks between environmental conditions and biogeochemical processes that are critical to tree health, along with root growth, distribution, and function. Geophysical methods are widely used to study subsurface prope rties that are difficult and costly to directly measure at field scales. Amongst a range of geophysical tools, electrical resistivity (ER) is particularly useful to map shallow hydrogeological features with strongly contrasting electrical properties (Loke et al. 2013) . Cellulose is an inherently resistive material, varying with moisture content and species - specific cell structure (al Hagrey 2007) . Measurements of living tree trunks in temperate regions of North America found resistivity values ranged from 1 0 2 - 10 5 m and was correlated with both the diameter and species (Gora and Yanoviak 2015) . In studies where root mass has been destructively sampled and compared to ER tomograms, strong positive correlations between subsurface resistivity and woody root mas s (WRM) were found (Amato et al. 2008a; Paglis 2013) . While there is considerable literature on ER monitoring of soil moisture changes due to root water uptake (Beff et al. 2013; Junliang Fan et al. 2015; Garré et al. 2012; al Hagrey 2007; Jayawickreme, Va n Dam, and Hyndman 2010; 81 Robinson, Slater, and Schäfer 2012) , there has been limited consideration of the influence of the cellulosic root material itself on ER measurements (Ehosioke et al. 2020) . We hypothesize that the presence of WRM in the upper 0.5 of the soil in a forest environment (Fahey et al. 1988) will result in high resistivity anomalies. We test this by accounting for the influence of physical, geochemical, and thermal differences on ER signals across a successional forest environment using a coupled hydrogeophysical model (e.g. Ferré et al. 2009; Hinnell et al. 2010; Moreno, Arnon, and Furman 2015; Tran et al. 2016) . Our novel coupled model has shown promise for directly quantifying and parameterizing hydro - bio - geo - thermal models of the subsu rface (Kuhl et al. 2018) (Kuhl et al. 2018) . In this approach, a forward hydrological and geophysical model predicts ER from a modeled VWC distribution. Benefits of this method includes the direct incorporation of hydrological, temperature, and chemistry d ata in the inversion, as well as the removal of smoothing errors and unrealistic parameter estimates often introduced through traditional ER inversions (Singha et al. 2015) . In this study, we used in - situ soil temperature and water content data as well as grain size, volumetric water content, and ER data measured from two deep (3.6 m) soil cores to build and validate 1D hydrological and petrophysical models of adjacent mature forest and unmanaged grassland instrumented with an ER array. A sensitivity analy sis was conducted to both establish a hypothetical relationship between WRM and ER, and quantify the influence of biogeophysical drivers (i.e., clay fraction, LAI, canopy storage, solar radiation, ET) across two contrasting vegetation zones on ER. Then we mapped the heterogeneous WRM fraction along the mature forest portion of the study site by inverting 2D ER data from three spring season datasets spanning three years (April 28 2017, April 18 2018, and May 16 2019) using a coupled hydrogeophysical inversio n method (Kuhl et al., 2018). The spatial distribution of WRM was 82 then evaluated relative to above ground biomass inferred from a tree diameter survey along the transect. Data collected during wet spring conditions prior to the onset of transpiration were chosen for this analysis to: a) maximize the contrast between the theorized resistive WRM and the surrounding wet soil, and b) avoid confusing resistive WRM with the increases in soil resistivity expected as ET processes cause soil drying. This research va lidates established relationships belowground biomass and aboveground biomass, and demonstrates an avenue to quantify WRM at the field scale in a host of forest settings. 2. Materials and Methods 2.1 Site Description The study was conducted in a successional forest ecosystem located at the Kellogg Biological Station (KBS) in southwest Michigan, USA. This post - glacial region is mostly flat - lying agricultural land overlying thick glacial drift sediments. Michigan has a h umid temperate climate (hot - and warm - humid continental climate zones Dfa and Dfb), with an average annual temperature of 9°C and 900 mm of precipitation at the study site. The soils at the site are classified as Typic Hapludalfs that fine upwards, with a loam A horizon, sandy clay loam Bt horizon that transitions to a sandy loam, and then loamy sand BC horizon, over a coarse sand parent material at around 1.3 m depth. Dissolution of carbonate minerals in the deeper (around 2 m depth), glacial deposits domi nate the porewater and groundwater geochemistry in the region, resulting in high fluid electrical conductivities around 550 µS/cm. An instrumented transect was established spanning from the northeast to the southwest (Figure 4.1b) in succession from a mat ure hardwood forest of oak, maple, and hickory with little understory, to a more juvenile stand of hardwood trees and dense underbrush, followed by a stand of autumn olive ( Elaeagnus umbellata ) shrubs, and terminating in an open unmanaged 83 grassland. The st udy presented here focuses on data collected at two 20 m - long sections of the transect, one centered in the mature forest (MF) (Figure 4. 1 c), and the other 120 m to the southwest in the open grassland (OG) (Figure 4. 1 d). Figure 4.1. Site map. a) Location map of the study site, the Kellogg Biological Station is indicated with a star, b) aerial view of the transect showing the agricultural nature of the region, the adjacent pond, and the relative locations of the mature forest (orange) and open grassland (yellow) portions of the transect c) photograph looking northeast in the mature forest, and d) photograph looking northeast in the open grassland. 2.1.1 Site Instrumentation and Data Collection Soil water content and soil temperature profiles wer e monitored every two hours throughout the study period from March 2017 through October 2019 at each end of the transect. Thermochron iButton DS1922L temperature sensors at 0.05, 0.25, 0.75, 1.25 and 2 m depth were installed in a 0.05 m diameter PVC pipe b uried vertically. Time domain reflectometry (TDR) water content sensors at 0.1, 0.2, 0.35, 0.5, 0.65, 0.90, and 1.25 m depth were installed in the sidewall of subsequently - backfilled soil pits dug adjacent to the transect, and data were recorded using an C ampbell CS615 datalogger. Two 0.038 m - diameter soil cores were collected on Oct 24 2019 in the MF and OG ~2 m 84 from the transect using a Geoprobe MT540 to 3.6 m depth. In - situ resistivity (at 25°C) of the soil in the cores was measured with a multimeter pr obe inserted in adjacent holes drilled into the sidewall of the plastic core liner at 0.08 m increments. The cores were subsequently cut open down their length and the soil was extracted in 0.08 m increments corresponding to the span of each in - situ resist ivity measurement, weighed and dried in a 105°C oven for 24 hours to measure gravimetric water content, then ground and sieved through 2 mm and then 1 mm sieves to measure gravel and coarse sand fractions. Volumetric water content was calculated from the g ravimetric water content using bulk density measurements obtained by dividing the mass of dry soil recovered from the volume of the core. Due to compaction from the core collection process a total of 36 soil samples were recovered with an average dry sampl e weight of 150 g. Percent sand, silt, and clay were then calculated from 10 grams of each dried sample using a Malvern Mastersizer laser. Percent organic carbon was also measured for each sample via loss on ignition at 430°C. Each 0.08 cm section of dried sample was tested for the presence of carbonate minerals using effervescence in 10% HCl acid as a proxy for the unleached carbonate interface. Electrodes for repeat ER surveys were buried in a trench 0.08 cm below the surface along the two sections of the study transect. Prior to burial, the location of each electrode was surveyed. Eight electrodes spaced 0.50 m apart and 14 electrodes spaced 1.5 m apart centered in each section were configured in Wenner array for a total of 61 measurements, with minimum a nd maximum AB spacings of 1.5 and 18 m respectively. The electrodes, made of cylindrical graphite rods 0.08 m long by 0.02 m diameter, were connected with CAT - 5e cable to a SuperSting R8/IP meter for each survey. The individual wires were attached to each rod through a small hole drilled at one end, then covered with a liquid polymer coating to protect from corrosion. Surveys were conducted on Apr 28 2017, Apr 18 2018, and May 16 2019, targeting 85 the period prior to the onset of the growing season. Air tempe rature at the time of each survey was 11.6, 0.4, and 20.2°C, respectively. A tree survey was conducted in Apr 2019 to measure the diameter at breast height and spatial location of each tree with a height over 1 m located within 10 m of the electrode array. 2.2 Modeling Methods To determine whether the influence of woody root mass in a forested environment causes elevated ER, we took the following steps: 1) calibrate hydrological and petrophysical models to in - situ conditions measured in the two cores and va lidate with TDR data, 2) test the sensitivity of the ER model to biogeophysical drivers that differ between the MF and the OG, including hypothetical woody root mass fractions, and 3) validate spatial estimates of WRM with an index of above ground biomass. To accomplish this, we used a coupled hydrogeophysical model (described in detail in Kuhl et al., 2020 (in review) ) with HYDRUS - 1D , and the Boundless Electrical Resistivity Tomography (BERT) model (Rücker, Günther, and Spitzer 2006) that solves for the potential field under an induced current. In this approach, the spatial distribution of water content ( ) and soil temperature ( ) is mo deled (using HYDRUS), and then translated to resistivity ( ) using the Waxman and Smits (1968) modification of (1) where is the porosity, is the pore water conductivity, and is the grain surface conductivity. and are the cementation and saturation exponents, respectively and typically 86 range between 1 and 2. The modeled resistivity at 25°C was temperature corrected using the equat ion: (2) where is the resistivity at 25°C and is the resistivity at discrete soil temperatures ( ) ( Hayley, Bentley, Gharibi, & Nightingale, 2007) . The resistance under an induced current is th en modeled (using BERT) and compared to the field - scale measurements of ER. 2.2.1 Hydrological Model Parameters Hourly precipitation, air temperature, ground surface temperature, wind speed, solar radiation, and humidity data and a calculated reference po tential evapotranspiration (via the FAO Penman - Monteith method) used to drive the hydrological model were retrieved from a climate station (maintained by KBS, and hosted by the Michigan Automated Weather Network, station ID: kbs) located <1 km from the sit e. The major ion concentrations of the incoming precipitation were taken from the summer seasonal average values recorded by the National Atmospheric Deposition Program National Trends Network station MI26 located at KBS. Using the 0.08 m resolution soil d ata from the cores, we aggregated by like grain size distribution into 20 - layer hydrological models in HYDRUS - 1D for both the MF and OG ecosystems. The models were run with an hourly timestep beginning on Nov 1 2010 to provide a sufficient ramp up before t he analysis period to stabilize simulated deep soil conditions. We used the Van Genuchten - Mualem (VGM) equation with hysteresis in the water retention curve, setting initial VGM model parameters of each layer using the percent sand, silt and clay measured by the grain size analysis of the soil cores and the Rosetta database (Schaap, Leij, and Genuchten 2001) . The Rosetta - the drying alpha parameter was set to 50% of the wetting alpha val ue. Along with the heat and 87 CO 2 transport, and root water uptake features of HYDRUS, we used the Unsatchem module to simulate the production and transport of the major ions, primarily Ca 2+ released by CO 2 - driven carbonate dissolution of the unleached carbo nate zone. Upper boundary conditions of the HYDRUS model included downward fluxes of infiltrating precipitation with a 0.02 m maximum head before triggering runoff, and upward fluxes of evapotranspiration. Potential evapotranspiration and canopy intercept ion and storage were estimated from leaf area index (LAI) using the method described in Kuhl et al. (in review). On - site leaf area index (LAI) data were not available for this site, therefore we compiled 500 m resolution four - day MODIS LAI data from a nine - county area surrounding the study site to estimate dynamic LAI for the forest and grassland pixels withing the site. The 30 m National Land Cover Dataset (NLCD) was overlain on the MODIS grid to identify the fraction of each land cover within each MODIS c ell. 123 cells with greater than 95% deciduous forest, and 9 cells with greater than 95% grassland were identified and used to calculate temporal LAI curves for the study period. A twenty - day moving median filter was used to smooth the values of the 123 an d 9 cells for each land cover, and a smoothing spline was used to downscale to hourly data for input to the hydrological model. Median summer (June, Jul, Aug) LAI values for the three model years were 5.03 +/ - 0.24 and 1.88 +/ - 0.89 [ - ] for the MF and OG, respectively. Confidence intervals were calculated as the variance in the MODIS cells. A constant negative flux boundary condition of 0.004, and 0.0047 cm/hr was set as the bottom hydrological boundary for the MF and OG models respectively. The flux rate w as calculated to maintain steady groundwater levels over the eight year model period, which is in agreement with estimates of annual groundwater recharge rates for southwest Michigan (Holtschlag, 1997) . The initial water table was set to 9.2 m based on the surveyed depth to the 88 adjacent pond surface. 2.2.2 Petrophysical Model Parameters The relationship between water content and resistivity (Equation 1) was initialized using measurements from each core. The measurements collected below the carbonate interf ace were distinct from those above (Figure 4. 2 ), therefore two sets of and parameters were optimized, one for each zone. We used the modeled pore water conductivity on the date of the core collection and assumed a grain surface conductivity, , of zero. Outliers more than one standard deviation from the moving median of the log resistivity (sorted by water content) were removed. Figure 4.2. Petrophysical data. Petrophysical relationship between the in - situ soil core resistivity and the volumetric water content (VWC) for the mature forest (MF) and open grassland (OG) cores, along with two fitted relationships (solid lines); one for samples below the carbonate in terface (points denoted with a gold center, and gold line), and one for samples above the interface (orange line). Measurements from the carbonate - rich zone generally plotted at lower resistivities than equally dry samples from above the interface due to t he influence of ions released through carbonate dissolution processes. 89 2.2.3 Modeling Woody Root Mass To test the effect of hypothetical WRM on ER measurements, we used the forward modeled resistivity from the hydrological and petrophysical models and replaced a fraction of randomly - - resistivity (>1000 e 4. 3 ). For simplicity, we assumed that the coarse WRM cells were uniformly distributed in the shallow zone between 0.05 and 0.65 m. We created six scenarios to test the effects of WRM: 1) zero WRM, 2) 20% WRM at a resistivity of 5000 m, and four with th e following assumptions to produce end - member scenarios: low resistivity (1000 m) roots at 10 and 40% of the shallow (< 0.65 m) soil volume, and high resistivity (10000 m) roots at 10 and 40% of the shallow soil volume. Figure 4.3. Modeled electrical r esistivity cross - section. Cross - section of the modeled electrical resistivity (ohm - m) along the mature forest (MF) portion of the transect with 25% of the shallow (0.05 - 0.65 m) soil volume randomly occupied by high - resistivity (5000 ohm - m, yellow color) m between vertices. The increase in soil matrix resistivity from ~150 m (dark blue color) to 2000 m (bright green color) corresponds to the transition from the finer soil t o the coarser sand parent material. 90 2 .2.4 Model Calibration and Validation To calibrate the hydrological and petrophysical parameters of the MF and OG, we ran the hydrological models of each ecosystem beginning on Nov 1 2010, and exported the 1D water con tent at the date and time when the soil cores were collected (Oct 24 2019). We then adjusted the residual and saturated water content parameters in each layer to adjust the modeled water content to fit the measured water content in each core. The modeled w ater content was validated with TDR measurements from the date and time of the core collection. After validation of the hydrological model, we re - ran the model with updated hydrological parameters and compared the simulated resistivity (without temperature correction) to the measured in - situ resistivity (taken at 25°C in the lab) at 0.08 m intervals. Three layers in the Bt horizon had higher clay contents and much lower measured in - situ resistivity values, therefore we updated the term (Equation 1) fro m an initial value of zero to a value that scaled with the clay content, , using a power law relationship, where (3) the and maximum parameters were then adjusted to improve the agreement between the simula ted and measured core resistivity values. We validated the model for the three spring survey days by comparing the modeled water content, temperature, and ER on each survey day to the field observations, from both the OG and MF. 2.2.5 Driver Sensitivity T o test the sensitivity of the modeled to each hydrobiogeophysical driver that may contribute to the observed differences between the MF and OG ER data, we used the OG model as our control model and altered inputs one at a time to test the effect of the MF model drivers. 91 The tested drivers included: topography, which differed slightly at the two ends of the transect; LAI, which controls evapotranspiration; canopy interception and storage; soil texture, and therefore hydrological properties, which control wat er content; soil surface temperature; petrophysical parameters, which control the relationship between water content and ER, and; pore water conductivity. For each run, a single input was varied while all other inputs were kept constant at OG values betwee n the two models. Lastly, we tested the effect of all drivers combined and the effect of hypothetical WRM fraction scenarios and compared against the MF ER data from the three spring surveys. 2.3 Mapping Root Mass Fraction After establishing the reasonabl e range of effects due to the presence of WRM, we then estimated the WRM fraction in each of ten 2 m lateral intervals in the MF section of the transect by minimizing the difference between the measured and modeled ER using a shuffled complex evolution (Du an, Sorooshian, and Gupta 1992) global optimization algorithm, where the objective function, ( 4 ) is the sum of the absolute values of residuals between the log measured ER ( ) and log modeled ER ( from three spring ER surveys (183 total measurements). We then validated the 10 spatial estimates of WRM with the assumption that below ground biomass correlates with above ground biomass. To quantify the proximity and size of the trees that were surveyed along the transect we used a calculated tree influence index ( ) from Yanai, Park, and Hamburg (2006 ), where for every 0.2 m increments along the line, : ( 5 ) 92 and is the tree number from 1 to 82, is the diameter at breast height, and and are the along - line and off - line position of each tree Assuming that tree roots do not extend outward infinitely, we truncated the summation to trees w ithin increasing ranges from three to 10 m of each point, , to test the influence of the maximum root range on the comparison. To compare our 2 - m resolution WRM estimate to an average tree influence index, we calculated the mean of the 100 0.2 m resoluti on values within each 2 - m zone for each of the ranges between three and 10 m and used a linear regression model to test for significance in the correlation between estimated WRM and . 3. Results 3 .1 Model Calibration and Validation Analysis of the soil cores show these soils fine upwards, with an organic - rich surface horizon, clay - rich Bt horizon, and sandy parent material (Figure 4. 4 b). The OG core parent material had 1 to 2% more clay than the MF core, which explained some of the d ifferences in water content. The transition from sandy loam to loamy sand seen in the percent clay distribution is clearly observed in the TDR data and the modeled water content, with a decrease in field capacity of around 10% occurring between the 35 and 50 cm depth sensors in the MF, and between the 50 and 65 cm depth sensors in the OG (Figure 4. 4 c). The Rosetta - parameterized hydrological models fit the core and TDR data well, with minor adjustments to the residual and saturated water content, which were lowered in some of the layers . Ov erall, the percent clay correlated well with water content and resistivity. Small peaks in water content observed around 2.5 and 3.5 m were explained by a slight decrease in percent sand in those interbedded layers and were replicated by the model. 93 Figure 4.4. Core profiles. Depth profiles of the core data (markers) and calibrated modeled outputs (solid lines) of the open grassland (OG - green) and mature forest (MF - blue) on the date of core collection (Oct 24 2019); a) enlarged photos from the soil cores representing examples of each soil horizon; b) percent clay content, c) volumetric water content (VWC), d) inverse of soil electrical resistivity (1/ER), e) pore water electrical conductivity (EC). Note that the box plots in panel e represent lysimeter chemistry data collected at nearby study sites at several times. Soil texture, which controls water content strongly influenced the increase in depth of the ER data, which was much more conductive (less resistive) near the surface in the wetter soil, as is expected based on the petrophysical relationship (Figure 4. 2 ). The petrophysical surface resistivity term, , was adjusted to be nearly half in the OG than the MF (0.017 versus 0.004) despite clay content being slig htly higher in the MF and water content being nearly identical in this layer. Additionally, the discrepancy in water content between the two cores that occurs between 1.5 and 2.5 m is also observed in the in - situ resistivity data. Further, while water cont ent is either consistent, or continues to decrease from 1.5 to 3.5 m, the EC actually increases, which we attribute to the presence of the high EC pore fluid below the carbonate interface coinciding with this depth (Figure 4. 4 e). Although pore water EC data was not collected from the study site during the survey period, a statistical summary of data from lysimeters across the KBS property is presented to show the general agreement of the modeled EC with the bimodal distribution 94 observed in the data due to the presence of carbonate minerals in the unleached deeper sediments (Figure 4. 4 e). The hydrological and petrophysical models calibrated from the core data were validated using the measured water content, temperature, and ER d ata on the three spring survey days (Figure 4. 5 ). Due to the timing of the measurements prior to the onset of the growing season, water content is similar across all three years and comparable between the MF and OG. The model again captures the trend of de creasing water content with depth that is driven by the soil texture (Figure 4. 5 a - c). The model also does an excellent job of modeling soil temperatures across the two models and three years (Figure 4. 5 d - f). While shade provided by the MF canopy lowers sur face temperatures in the summer by as much as 6°C relative to the OG, during the three spring surveys prior to leaf - out, shallow soil temperatures in the MF were only roughly 1°C warmer than in the OG, likely due to increased insulation from tree cover and /or leaf litter. A gap in surface temperature data for the OG site in 2019 required the use of air temperature as a proxy which led to an over - prediction of soil temperature by several degrees (Figure 4. 5 f). However, the temperature is quite variable from year to year, being much colder during the early (Apr 18) 2018 survey, and warmer in the 2017 and 2019 surveys which occur a bit later in the year, and show significant warming near the surface. 95 Figure 4.5. Model validation. Model validation of three survey days spanning three spring seasons, comparing modeled and measured a - c) soil water content, d - f) soil temperature, and g - i) apparent resistivity with depth for the open grassland (OG) and mature forest (MF) ecosystems. Temperatures were considerably colder on the 2018 survey day resulting in higher overall resistivity in both OG and MF models, although soil water content was consistently high in all three. 96 Comparison of the modeled and measured ER shows the calibrated hydrogeophysical model of the O G fits the data well across all three surveys (Figure 4. 5 g - i). The increase in ER with depth is due to the decrease in water content with decreasing clay percentage, despite the increase in pore water EC, which has the opposite effect on resistivity (Equat ion 1). Elevated resistivity in the 2018 survey relative to the other years can be attributed to the colder temperatures. The ER model of the MF however does not perform as well, with considerable mismatch in the modeled ER at the two shallowest depths. Th e model does capture the overall increase in ER between the MF and OG which can be attributed to the much lower MF grain surface conductivity term which controls the resistivity of the clay - rich loam horizon near the surface. 3 . 2 Driver Sensitivity To quantify the influence of each of the hydrobiogeophysical drivers of the MF and OG models, including a hypothetical WRM fraction, we compared the effect on ER of altering each driver one at a time (Figure 4.6) . Starting with the OG model as the control, we ran scenarios that included using the MF topography, LAI, soil texture, soil surface temperature, surface conductivity, and the pore water EC as well as a combined scenario. A summary of the residuals between the model and the MF ER data are presented in Table 4. 1 for each tested driver. The topography had almost no effect while the temperature, pore water EC and LAI each had a small positive effect, reducing residuals by ~15%. Adding a WRM fraction of 20% to the model increased the resistivity and reduced residuals by ~20%. The differences in soil texture had a moderate effect, mostly at depth where MF soils are drier, and therefore increased ER overall and reduced residuals by ~40%. The petrophysics had the largest effect, increasing the shallow 97 ER consid erably due to the large decrease in the grain surface conductivity term, and reducing residuals by almost 60%. Combining all drivers, with a WRM of zero, reduced the residuals by 70% and the addition of WRM fraction resulted in a range of effects that all increased the ER and generally decreased the residuals further (with the exception of the extremely high scenario), with the best fit occurring in the scenario with 20% WRM fraction at a root resistivity of 5000 ohm - m. The fraction of WRM had a larger eff ect than the resistivity of the roots, with a factor of four increase in WRM having a much larger effect than a factor of 10 increase in the resistivity of the root fraction. 98 Figure 4.6. Sensitivity to model drivers. Sensitivity of the ER model to eac h of the hydrobiogeophysical drivers for the a/d) 2017, b/e) 2018, and c/f) 2019 model years. The first row shows single - driver scenarios, while the second shows combined drivers with WRM scenarios. Depth averaged ER data with confidence intervals for the open grassland (OG, green) and mature forest (MF, blue) are shown along with the modeled ER for the seven different driver scenarios. 99 Table 4.1 . Summary of driver scenarios. Summary of root mean square error (RMSE) between depth - averaged mature forest (MF) apparent resistivity (AR) data and model driver scenarios for each model year. Part A compares the control open grassland (OG ) model with the addition of individual forest drivers, Part B compares the combined MF model with the addition of individual woody root mass (WRM) scenarios. A Control Topo - graphy Soil Surface Temp . LAI Pore water EC 20% roots @ 5k ohm - m Soil Texture Grain surface cond . 2017 147.7 126.1 109.9 97.6 97.28 84.9 55.5 31.9 2018 169.6 195.2 173.6 158.8 165.4 153.8 113.1 97.1 2019 151.1 158.0 126.3 121.7 132.4 123.7 102.8 74.0 Avg. 156.2 159.8 136.6 126.0 131.7 120.8 90.4 67.7 B 0 % WRM 10% WRM @ 1k ohm - m 10% WRM @ 10k ohm - m 20% WRM @ 5k ohm - m 40% WRM @ 1k ohm - m 40% WRM @ 10k ohm - m 2017 35.8 34.2 34.5 40.8 58.53 87.0 2018 71.8 60.8 55.9 42.2 36.06 47.3 2019 39.4 30.2 27.1 21.1 32.2 57.7 Avg. 49.0 41.7 39.1 34.7 41.9 64.0 100 3 .3 Mapping Root Mass 37 trees were located within 10 m of the MF transect, with an average diameter at breast height of 0.2 ± 0.06 m, and an average distance from the line of 2.98 ± 0.66 m (Figure 4. 7 a). The estimated WRM fraction in 2 m increments along the MF transect result ed in a range of estimates between 0.08 and 0.52 (mean 0.34 ± 0.08) (Figure 4. 7 c). A positive correlation between WRM and tree influence index was found for all of the tested root extent limits, however this relationship was only statistically significant (p = 0.01) at root extents greater than 7 m (Table 4. 2 ). Table 4.2. Summary of driver scenarios. Pearson correlation coefficient (PCC) and p - values between estimated along - line woody root mass (WRM) and the calculated tree influence index for increasing root extent ranges from 3 to 10 m. p < 0.01 are denoted in bold. 3 m 4 m 5 m 6 m 7 m 8 m 9 m 10 m PCC 0.19 0.33 0.24 0.32 0.46 0.78 0.83 0.84 p - value 0.59 0.36 0.51 0.37 0.19 0.008 0.003 0.002 101 Figure 4.7. Woody root mass estimate vs. tree index. a) Plan view of the trees (green circles) and electrodes (pink triangles) in the mature forest transect section. A three - meter buffer indicating the minimum root extent for calculating tree influence is indicated with the gold dash. The size of th e circular markers is scaled linearly with the diameter at breast height of each tree. Note the closer spacing of the electrodes centered on the 24 m mark, and the exaggerated scale on the x - axis relative to the y; b) 2 m average tree influence index from the 10 - m maximum root extent range scenario and c) estimated woody root mass (WRM) fraction for each 2 m section of the transect; d) estimated WRM fraction versus tree influence index for the 10 - m root extent range scenario. 4. Discussion The results of the sensitivity analysis demonstrate several important points regarding the influence of the various hydrobiogeophysical drivers on ER. Using a coupled hydrogeophysical model calibrated using in - situ data from cores taken from two different ecosystems, we found that neither differences in ET processes, soil texture, water content, , pore water EC, nor temperature alone or in combination could account for the consistently elevated and distinct heterogeneity of ER observed in the MF environment relative to the OG during non - growing season periods. Despite similar clay content observed in the shallow soil of the two ecosystems, 102 the grain surface conductivity term of the petrophysical relationship (Equation 1) relating water content and resistivit y was found to be much higher in the OG, possibly due to combination of differences in land use history, soil compaction, bioturbation rates, or microbial processes that have altered the soil structure and/or mineralization. This difference alone elevated the MF ER relative to the OG; however, neither it, nor combining all the MF drivers together resulted in as good of a match with the MF ER data as was found between the model and OG ER data. Further, the grain surface conductivity could also not reasonably explain the increased heterogeneity in the MF ER unless this property itself is more heterogeneous in the MF than the OG. Other explanations for this distinct heterogeneity in the MF may be the cumulative effect of year to year differences in water conten t and sunlight controlled by the canopy structure, which is absent in the OG, and was not incorporated into the hydrological model. Other factors such as the natural lateral variability in microtopography, depth of the loam or carbonate interfaces, or petr ophysics we would expect to persist in the OG as well, however this signal was not observed in that portion of the transect. However, incorporating reasonable hypothetical estimates of WRM into the forward model with all the other MF drivers caused increas es in ER that improved the model fit and added heterogeneity that was well aligned with the observations. The optimized WRM fraction results support our hypothesis that WRM would increase shallow ER, vary considerably spatially, and be positively correlate d with above - ground biomass. The strong positive relationship between the tree influence index and our estimate of WRM at high root extent limits implies that the coarse woody roots extend a considerable distance from the tree trunk which is supported by t he literature, where destructive sampling found stronger correlations between tree influence index and coarse root biomass at extents between 4 and 6 m (Yanai, Park, and Hamburg 2006) . Without destructive sampling at our site 103 we are not able to calibrate o ur estimated WRM fraction to the in - situ woody below - ground biomass, nor confirm the resistivity of the roots which likely varies considerably with root diameter and age as is the case above - ground (Gora and Yanoviak 2015) . Therefore, we cannot infer from our WRM estimates the true fraction of the subsoil that is comprised of WRM. We also did not account for the depth - dependent exponential decay in root biomass that is often observed in destructive sampling and likely could affect the WRM estimate, however there is a path via the forward modeling to parameterize this feature in the future. Even with the simplified approach, our results confirm prior research relating a tree influence index with WRM (or a signature of WRM that elevates ER), as well as a metho d to survey WRM at large scales, with minimal disturbance to the ecosystem, to either enable long - term monitoring, or to guide for targeted sampling of root biomass. We chose to conduct this study when conditions were as similar between the two ecosystems as possible to remove uncertainty in the potential drivers and isolate the effect of WRM. Repeating the analysis during the height of the growing season when water content and temperature conditions are most dynamic and distinct between the two would be h elpful to quantify the sensitivity of the coupled model approach to uncertainty in dynamic model inputs such as evapotranspiration, temperature, and precipitation, as well as static inputs such as the petrophysics and soil texture. We also anticipate that the estimate of WRM could be used to inform the distribution of finer roots, which may exert a large control on root water uptake and could deepen growing - season MF heterogeneities. Additionally, while we focused this study on a narrow stretch of the MF tr ansect to limit the influence of other factors in our interpretation of the data, optimizing WRM along the remainder of the successional forest and comparing to the tree influence index could yield important insights into the relative WRM fraction across t ree 104 stands of varying age. This paper presents a coupled hydrogeophysical inversion that accounts for the multiple influences on resistivity in two adjacent but distinct ecosystems across three years to test the sensitivity of ER measurements to hypothetic al WRM in a temperate deciduous forest environment. While we were not able to validate the estimates of WRM without destructive sampling at the study site, we infer from quantified above - ground biomass that the presence of WRM can explain the magnitude and heterogeneity in the ER data observed in the MF portion of the transect. The results of this research highlight the power of the coupled approach to better understand the critical zone and a path forward to study important, yet poorly constrained, subsurf ace features, such as WRM. Acknowledgments This material is based upon work primarily supported by the US Department of 2015 - 68007 - 23133, with additional support from the Great Lakes Bioenergy Research Center, U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research under Award Numbers DE - SC0018409 and DE - FC02 - 07ER64494; the National Science Foundation Long - term Ecological Research Program (DEB 1832042) at the Kellogg Biological Station, and by Michigan State University AgBioResearch. We would like to thank the following undergraduate researchers for their substantial contributions towards data collection for this study; Julia Michie nzi for her dedication to site maintenance, trips to the study site for ER data collection, and her efforts along with those of Joe Szkutnik in the tree survey; as well as Amanda Liddle for her leadership on the project collecting and analyzing the soil co res that proved to be pivotal in bringing this story together. 105 CHAPTER 5: PREDICTING BIOGEOPHYSICAL DRIVERS OF WATER CONTENT AND ELECTRICAL RESISTIVITY HETEROGENEITY IN A SUCCESSIONAL FOREST ENVIRONMENT Abstract Shifts in landcover are expected to persist througho ut the next century as humans continue to adapt to changing energy, industry, and agricultural demands. To better predict how these changes will influence the water balance at regional and global scales, it is important that we develop new tools to study t he relationship between soil water content dynamics and above ground vegetation. In this paper we present the use of electrical resistivity data to monitor heterogeneity in soil water content in a successional forest environment throughout a growing season period, and propose the use of a coupled hydrogeophysical inversion method to test hypotheses and parameterize inputs to explain the observed spatial and temporal trends. Using a linear regression analysis we found significant correlations between indices of above ground biomass, including leaf area index and tree diameter, and patterns of electrical resistivity change. Incorporation of spatial variation in evapotranspiration and canopy interception into the model however was not sufficient to explain all of the observed heterogeneity, indicating the presence of some additional process or structural features controlling water content or resistivity. 106 1. Introduction Land cover in the Midwestern United States has undergone radical shifts over time beg inning with the near complete deforestation of the region in the late 19 th century. Since that time, much of that land transitioned to agricultural production along with urban and industrial development, and many marginally - productive lands have since prog ressed through stages of forest succession. These radical shifts in land cover have played a role in widespread changes in nutrient and water cycling, as well as regional surface temperatures. To understand the magnitude and scope of those changes and pred ict future impacts in the face of ongoing climate and land use change, landscape hydrology models must be able to accurately account for the water exchanged between the Earth and atmosphere. A key component of this process is the root distribution, which e nables a plant to access available soil water to meet transpirative demand. Root distributions are highly variable across the landscape, with differences linked to factors including plant species, plant age, surface slope, water table elevation, and soil h ydraulic properties. Despite its importance to landscape hydrology modeling, the state of root phenology knowledge remains limited by inherent challenges to studying roots (Sinacore et al. 2017) . Current methods to study patterns of root growth and root d istribution involve either growing roots in controlled environments with rhizotrons (root windows) or carefully excavating roots by hand from the field. Advancements with dye tracer studies coupled with rhizotrons have helped understand how roots preferent ially uptake water. While such methods have improved knowledge of fine root phenology for annual agricultural crops, few studies have been conducted with larger perennial trees and grasses. Therefore, we need new tools that allow root water uptake processe s to be characterized and incorporated into models s at the field scale. 107 Electrical resistivity (ER) data can quantify changes in soil moisture in response to root processes including hydraulic redistribution (Robinson, Slater, and Schäfer 2012) and evap otranspiration (Dick et al. 2018; Jayawickreme, Van Dam, and Hyndman 2008; Whalley et al. 2017) . ER also has the benefit of sampling large areas with limited disturbance of the soil and root structure, and the data does not require extensive time and physi cal effort to collect, unlike root excavation. New forward modeling approaches that couple hydrological and geophysical models have shown promise for estimating the physical properties of the hydrological system, including recent work that estimated hydrau lic and thermal conductivity (Moreno, Arnon, and Furman 2015; Tran et al. 2016) in field settings. Coupled hydrogeophysical modeling methods have also been developed to estimate the distribution of root water uptake and coarse root distribution (Chapters 2 and 3). Here, we utilize a novel forward - coupled hydrogeophysical model method to develop and build a framework to test process - based hypotheses that can explain observed heterogeneity in shallow ER measurements in a successional forest transect in south west Michigan, USA. We use a linear regression modeling approach to test the relationship of lateral ER data variability with factors that may influence growing - season water content dynamics such as leaf - area index (LAI), diameter at breast height (DBH) of adjacent trees, and surface topography. Using the results of the statistical modeling, we propose several hypotheses regarding the process and extent to which each factor may influence the observed resistivity. Then we compare modeled ER heterogeneity fro m the 2017 growing season to the data and quantify how well the modeled variability correlates with the observed. Finally, we test the influence of the 2D hydrological model resolution on the site scale water balance partitioning for a five - year period fro m 2012 - 2017. 108 2. Methods 2.1 Site Description The study site for this research is located at the Kellogg Biological Station in southwest Michigan, USA where the terrain formed during the advance and retreat of the Saginaw Lobe of the Laurentide Ice Sheet. The site has several vegetation zones, in succession from: 1) well - established closed - canopy (CC) hardwood (oak and maple) forest with little understory on the northeastern end, 2) less mature open - canopy (OC) forest with a heavy understory, 3) a stand of shrubs (SH) primarily comprised of the invasive Autumn Olive species, and 4) open grassland (GR) dominated by native grasses (Figure 5. 1 c). Historical aerial imagery indicates that the CC forest has been present since at least 1938, while the remainder of the transect, including the OC forest, appears to have been abandoned from agricultural production after 1970 (Figure 5. 1 d - g). Figure 5.1. Site map. Site map a) inset map of Michigan and the Great Lakes Basin, b) inset map of the Kellogg Biological Station research complex with the transect marked in red and soil map units in orange (source: USDA NRCS), c) photo from the grassland portion of the trans ect looking northeast towards the shrub and forest, d - g) aerial imagery of the transect showing the progression of field abandonment and forest succession between 1938 and 2012. Note the sandier Oshtemo Sandy Loam (OsD) land cover is primarily forest. 109 The region is covered in tens of meters of glaciofluvial deposits that are rich in carbonates on top of bedrock. Borehole logs of the regional glacial geology by Kehew et al. (2013) as well as nearby well logs, indicate the presence of sand and gravel lenses t o more than 50 m depth. Lysimeter studies and soil pits spanning the Kellogg Biological Station property indicate the depth of the carbonate - leached weathering front ranges between 1.5 - 2.5 m, due to preferential flow and micro - topography. CO 2 - driven carbon ate dissolution at this interface contributes to high porewater conductivities between 50 - 250 µS/cm above the carbonate layer and 400 - 700 µS/cm below (Jin et al. 2008) . Groundwater conductivity is consistently near 550 µS/cm. The soils at the site are clas sified as an Oshtemo Sandy Loam (OsD) to the north in the CC forest, transitioning to a Kalamazoo Loam (KaB) along the remainder of the line to the south. Theses soils are Typic Hapludalfs that fine upwards, composed of a loam A horizon, sandy clay loam an d sandy loam Bt horizons that transition to a loamy sand BC horizon, and terminate in a gravelly sand parent material at approximately 1.3 m depth. A survey of the surface elevation of a kettle pond located 200 m east of the transect (Figure 5. 1 b) during Apr 2019 indicated the groundwater elevation was approximately 9.2 m below the lowest point along the transect. Data from a pressure transducer in a well ~1000 m east of the transect shows groundwater elevations are stable from year to year, with growing season declines of approximately 0.5 m. 2.2 Instrumentation & Data Collection The study site was instrumented with 112 electrodes (cylindrical graphite rods 0.08 m long by 0.02 m diameter) spaced 1.5 m apart for a total array length of 166.5 m, w hich spans the transition from the CC forest to GR. 48 additional electrodes were interspersed in groups of eight, at six high - resolution locations 24 m apart along the transect. All electrodes were buried in 110 a trench just below the surface and connected w ith CAT - 5e cable to a central take - out location for repeat surveys during the multi - year study. The precise location of each electrode was surveyed with a total station at the time of installation. Maximum elevation gain along the transect is 1.5 m. Primar y ER data collection consisted of 1823 measurements in a Wenner configuration, with a minimum AB spacing of 4.5 m and a maximum of 126 m. These data were collected at roughly 2 - week intervals throughout the 2017 and 2018 growing season (May to November) wi th an AGI SuperSting R8 Resistivity/IP meter. 60 high - resolution ER measurements (in Wenner configuration) were also collected during each site visit using the more closely spaced electrodes, with AB spacings of 0.75, 1.5 and 2.25 m. Adjacent to each of th e six high - resolution ER data locations, five temperature sensors (Thermochron iButton 9255) were installed vertically in PVC housing at depths of 0.05, 0.25, 0.75, 1.25 and 2 m for a total of 30 sensors. Temperature was measured hourly with 0.05°C resolut ion throughout the study period. Time domain reflectometry (TDR) sensors were also installed adjacent to the line, at three of the six high - resolution locations (24, 96, and 144 m from north to south or CC to GR) and at depths of 0.1, 0.2, 0.35, 0.5, 0.65, 0.9, and 1.25 m. The TDR sensors were installed horizontally into the sidewall of 1 x 1 x 1.25 m soil pits and calibrated with gravimetric soil moisture measurements collected from the pits during installation. The hourly data was filtered using a 24 - hour moving median and standard deviation to reduce data noise. Outliers exceeding + - one standard deviation from the moving median were removed from each dataset. Soil samples were collected with a push probe at 0.2 m intervals to a depth of 1 m at the six hi gh resolution locations along the line to determine the percent sand, silt, and clay via grain size analysis. Two 3.6 m deep soil cores collected at each end of the transect were analyzed at 111 0.08 m increments for percent, sand, silt, clay, gravel fraction, gravimetric water content, and in - situ ER. A tree survey along the transect was conducted on Apr 19 2019. Diameter at breast height [m] and the spatial location of each tree was measured for 82 trees found within 10 m perpendicular to the transect. Other data, such as climate information (Figure 5. 2), was obtained from a climate station located on the Kellogg Biological Station property hosted by the Michigan Automated Weather Network (station ID: kbs). In the absence of daily field LAI measurements for t he model period, LAI data was obtained from remote sensing products. Data from the MODIS satellite with 4 - day 500 - m resolution cells comprised of at least 95% deciduous forest or grassland per the 30 m National Landcover Dataset (NLCD) were averaged within each landcover and downscaled to produce hourly LAI estimates for the forested (CC, OC, SH) and non - forested (GR) portions of the transect. To benchmark along - line variability in canopy cover, LAI measurements were collected every 2 m with a handheld LAI - 2200C Plant Canopy Analyzer (LiCor Biosciences) on Sept 1 2017. These LAI measurements were inferred to represent peak LAI along the line and used to rescale the temporal remotely sensed LAI product at 2 m resolution. Measured LAI in non - grass zones ranged from 1.9 to 6.63 (median 4.90 [ - ]). A maximum gap fraction at 2 m resolution along the line was estimated by scaling the LAI data from 1 to 0, for a maximum gap fraction of 100% in the GR and a minimum gap fraction of 0% in the CC (median = 21%). 112 Figure 5.2. Climate input data. Climate data from Nov 1 2012 Nov 1 2017, showing average weekly reference evapotranspiration (PET ref ) (green line), average weekly temperature (red line), and weekly precipitation totals (blue bars). 2.3 Model Parameterization We used a coupled hydrogeophysical model of the transect built using HYDRUS - 1D (J. and the Boundless Electrical Resistivity Tomography (Rücker, Günther, and Spitzer 2006) models. For details of the hydrogeophysical functions and model inputs, as well as calibration of the hydrological and petrophysical parameters for this study site, see details from a prior study of two portions of the same transect presented in Chapter 3. For this study we discretized the 168 - m transect into eighty - four, 2 - m wide zones. We used historical and present - day aerial imagery along with ground observations to classify the landcover of each zone as CC, OC, SH, or GR. The tree survey was used to calculate a tree influence index, presented in Yanai, Park, and Hamburg (2006) : (1) where is the diameter at breast height, and and are the along - line and off - line coordinates of each tree, and is 0.2 m increments along the entire transec t, and is the 113 maximum extent of trees included for the summation at each point, , assuming lateral roots do not extend infinitely from each tree. Based on the results from the excavation of both fine and coarse roots in -- Yanai, Park, and Hamburg (200 6) , and their correlation with , we calculated a fine - root and a coarse - root extent by limiting the extent to four and 10 m, respectively. We then calculated the average within each two - m wide zone. Using the linear relationship estab lished previously in Chapter 3 (Equation 2), we used the average to predict the fraction of woody root mass, in each zone, (2) The median of the non - zero woody root mass fraction estimates with this method was 17% with a maximum of 49%. One anomalous zone between 40 - 42 m was located <1 m from a coppiced tree with the largest diameter of all the trees surveyed, resulting in a WRM v alue of 137% which was manually reduced to 75%. For the purposes of this study, we limited the number of soil hydrological layers to five, grouping the high - resolution (0.08 cm) core analysis by soil texture classification. Observations from both the NRC S web soil survey, deep soil cores, and shallow push - probe cores indicate the depth of the loam layer increases along the line with the transition from the Oshtemo to Kalamazoo loams (Figure 5. 1 b). We validated the dynamic hydro - thermal component of the mo del by comparing the hourly time - series TDR and soil temperature data from the CC, SH, and GR zones, to the hourly modeled outputs for the 2017 season. 2. 4 Analysis of Variance To explore the relationship between heterogeneity in different biogeophysic al drivers and the observed ER data, we performed an analysis of variance between the lateral shallow ER and spatial variables driven by site characteristics, such as the tree influence index (with a 10 and 4 m 114 limit), the LAI, and the surface topography f or each ER dataset in the 2017 growing season. We also used the log difference of three pairs of datasets that represented early and late drying periods, from Apr 28 to June 10 and July 21 to Sept 15, and a rewetting period, Sept 15 to Oct 13. 2. 5 Modeled ER After testing the strength of the statistical relationships between the spatial variance in the shallow ER and the measured site characteristics, we tested how well the coupled hydrogeophysical model reproduced the observed heterogeneity. We enforced t he site - scale variability in LAI, gap fraction, woody root mass, and assumed the sandy loam/loamy sand interface was flat - lying across the study site at a minimum depth of 1 m from the surface. We compared the shallow modeled ER to the observed ER on each survey data as well as the change in shallow ER across the three different drying and wetting focus periods. 2. 6 Scale - Dependent Water Balance Comparison Finally, using the established zonal hydrogeophysical model with the relevant drivers of site heterogeneity in place, we tested the effect of model resolution and successional forest - type heterogeneity on the site - scale water balance. We used the 2 - m zone resolution as the baseline model and then decreased t he resolution to 12 and then 42 m, as well as a scenario where the transect is broken into two zones divided at the boundary between the SH and GR. The difference between annual growing season (Apr 1 to Nov 1) precipitation and model output of evaporation, transpiration, soil water storage in the upper five - m, and canopy storage was calculated for each zone, and integrated to evaluate the average water balance across each of the four spatial scales for a five - year model period from 2012 to 2017. 115 3. Result s and Discussion 3.1 Hydrothermal Model Validation To validate the temporal dynamics of the hydrological model for the 2017 growing season, we compared the simulated and observed (TDR) water contents at each of the seven depths, in each of the three zones. The model accurately captured spikes in water content due to precipitation events, and the rate of soil moisture drawdown due to drainage in the non - growing season, and evapotranspiration during the growing season, particularly in the upper 0.65 m (Figure 5. 3 a). The TDR data and model demonstrate the increasing depth of the transition from the sandy loam to much lower field - capacity loamy sand soil along the transect. In the CC forest zone, this boundary occurs between the 0.2 and 0.35 m depths, while in the OC and GR, it is between 0.5 and 0.65 m depth. This distinction agrees with the mapped boundary of the Oshtemo Sandy Loam and Kalamazoo Loam soil units which occurs near the edge of the CC, and likely has played a role in the land use history of the si te. The lower plant available water of the sandier soil could reasonably explain the inferred early abandonment of agriculture in the CC forest portion of the transect (Figure 5. 1 ). 116 Figure 5.3. Soil water content and temperature data. Observed (black marker) and modeled a) soil water content at 20 cm depth (blue line) with daily precipitation (grey bars) and b) soil temperature at 25 cm depth (red line) for the 2017 model year. Electrical resistivity survey days are indicated with vertical gold lines. The three focus periods on change in resistivity are indicated with green arrows. Similarly, we validated the thermal model by comparing the modeled and observed temperature along the transect at each depth from Apr 1 to Nov 1 2017 (Figure 5. 3 b). The mode l performed very well, in part due to the use of the shallowest temperature sensor as the upper boundary condition 3.2 Analysis of Variance The drivers of lateral heterogeneity in ER measurements was explored by analyzing the variance. To isolate later al heterogeneity, we extracted the shallowest set of ER measurements corresponding to the shortest geometry (AB = 4.5 m) and smoothed the data with a 10 - point moving median. We then ran a multi - variate regression analysis with the tree influence for fine r oots ( ) and coarse roots ( ), LAI, and surface topography data at 2 m intervals to 117 determine the statistical relationship with the ER data. Results of the analysis indicate there is a significant relationship between the observed shallow ER and the above - ground tree and canopy structure as well as the topography of the site, prior to, during, and after the growing season. During the pre - growing season (Apr 28 2017), the best predictor of ER variability in the shallow ER data (Figure 5. 4d) was th e increased proximity and size of trees ( , as shown in Table 1. Prior to leaf - out, and in the absence of appreciable evapotranspiration rates, this appears to be due to the presence of resistive coarse woody roots in the shallow subsurface, a relati onship that has been demonstrated elsewhere (Amato et al. 2008b) . Early in the growing season (Jun 10), , LAI, and topography are positively correlated with increased lateral resistivity. Likely explanations for this trend are the increase in root wa ter uptake and canopy interception associated with the presence of greater leaf area, which decreases water content more quickly compared to areas with lower leaf area. A weaker positive correlation with surface topography was found that conflicts with the relationship observed in subsequent datasets. Later in the growing season, on July 21, as evapotranspiration dries the shallow soil, LAI is again the most significant, however during the driest period, on Sept 15, the topography is the only significant v ariable. We hypothesize this could be due to the relationship between surface topography and the thickness of the less - resistive loam layer. If the parent sand material is flat lying, increases in elevation would be linked to a thicker loam layer, which wo uld reduce the resistivity observed near the surface. Alternatively, there could be a prevailing trend in the depth of this sand layer that persists across the transect that is only partially captured by the surface topography. The last survey (Oct 13) was collected after a rainfall event that restored water content to field capacity, and evapotranspiration rates had reduced from their peak due to decreased LAI after a week of ongoing leaf senescence. In this instance, both the LAI and 118 topography were posit ively correlated with resistivity. Given the low LAI and reduction in canopy storage, this positive relationship can perhaps be attributed to interception from fresh leaf litter and woody stem area, rather than the remaining leaves in the canopy. The stro ngest relationship was observed between the change in log resistivity between the drier Sept 15dataset and the July 21 dataset two months earlier in the growing season. Across this time interval, all four spatial variables were significant predictors of wh ere greater increases in resistivity was observed. Notably, LAI and were both negatively correlated with these higher increases in resistivity, which has interesting implications. While we generally would expect resistivity to increase more where the re is greater LAI, and thus more fine roots with more root water uptake, this does not appear to be the case. One explanation is that when climate conditions become drought - like and the shallow soils have uniformly reached the wilting point, areas that hav e more vegetation are able to rely on the ecosystem function of hydraulic redistribution to survive, as has been observed elsewhere in ER studies (Robinson, Slater, and Schäfer 2012) ly become as dry earlier in the season preferentially showing signs of more drying later in the season or possibly an indication that the soils in the vicinity of the trees have different hydrological properties due to a secondary process related to the pr esence of the roots themselves. Patterns in the change of early season log resistivity followed a more predictable trend, with significant positive correlations with LAI and topography. Likewise, a significant positive trend was observed with the late sea son rewetting change in resistivity (10/13 minus 9/15), where areas showing the least decrease were correlated with higher LAI. As was the case for the static Oct 13 dataset, this could be related to interception and storage by fresh leaf litter if not tre e canopy, given the progressed state of senescence on that survey date. 119 Table 5.1 . Analysis of regression summary table. p - values meeting the p = 0.01 threshold are indicated in bold. Negative correlations are indicated with parentheses Dataset LAI Topo R 2 Fstat pvalue 4/28/17 (0.03) 4.5e - 5 0.04 0.36 0.405 8.69 1.99e - 05 6/10/17 (0.12) 2.2e - 4 7.57e - 5 0.009 0.486 12.1 5.72e - 07 7/21/17 0.888 (0.457) 0.006 (0.267) 0.21 3.47 0.014 9/15/17 (0.245) 0.0508 0.9 (0.009) 0.29 5.25 0.0013 10/13/17 0.156 0.093 0.0012 0.003 0.42 9.15 1.19e - 5 Early season drying log(6/10/17) - log(4/28/17) 0.20 (0.35) 0.0002 0.0005 0.42 9.09 1.27e - 05 Late season drying log(9/15/17) - log(7/21/17) (0.0049) 0.0007 (1.05e - 8) (4.62e - 6) 0.67 25.9 2.26e - 12 Rewetting log(10/13/17) - log(9/15/17) (0.265) 0.353 8.62e - 6 (0.034) 0.42 9.22 1.1e - 5 3.3 Modeled ER Results of the modeled ER output are shown in Figure 5. 5. We found the model captured some of the observed spatial heterogeneity and most of the temporal heterogeneity. Notably the large increase in resistivity during the driest period on Sept 15 and the associated decrease in resistivity after the Oct rewetti ng event were particularly well modeled (Figure 5. 5c). Incorporating the observed LAI, gap fraction, root mass and a fixed depth to the sandy soil layer alone was not sufficient to account for all of the spatial heterogeneity (Figure 5. 5a). The most persis tent and notable discrepancy between the observed and modeled resistivity occurs in the first 20 m of the transect, where despite a much higher LAI and lower gap fraction, the root mass 120 fraction was low, and the model did not predict an appreciable increas e in resistivity. Further, a distinct increase in resistivity around 150 m in the GR portion of the transect during the drying phases was not predicted by the model. 121 Figure 5.4. Spatial variables for linear regression. Plot of different heterogeneous data analyzed for correlation with the resistivity data: a) surface topography, b) leaf area index (LAI) and gap fraction derived from LAI, c) tree influence index with a range of 4 and 10 m, d) apparent resistivity ( m) from four surveys throughout 2017, e) change in apparent resistivity ( m) for three time periods, early drying (Apr 28 to June 10), late drying (Jul 21 to Sept 15) and rewetting (Sept 15 to Oct 13). 122 Figure 5.5. Hydrogeophysical model outputs. Output from hydrogeophysica l forward model of a) modeled (solid line) and observed (close circle) apparent resistivity (AR) along the transect for five survey days between Apr 28 and Oct 13 2017; b) modeled versus observed apparent resistivity (AR) for same; and c) modeled change (d elta) AR versus observed change in AR across three focus periods, early drying (Apr 28 to June 10), late drying (Jul 21 to Sept 15) and rewetting (Sept 15 to Oct 13). Positive change indicates an increase in resistivity from one time period to another. Dat es are displayed in mm/dd/yyyy format. 3.4 Water Balance Comparison We compared the different components of the water balance in each landscape zone using the hydrogeophysical model over a five - year period from 2012 to 2017. Differences in evapotranspiration along the transect was primarily driven by differences in LAI for the different 123 vegetation zones, which controlled the partitioning of potential evaporation and transpirati on. Within vegetation zones where LAI was the same, a thicker sandy loam layer resulted in an increase in soil water storage and greater evapotranspiration due to the higher water availability. Much lower LAI values in the grass zones caused transpiration to be the lowest and evaporation the highest of all zones. The increase in soil water storage due to the deep loamy sand layer in these zones meant that recharge was not much higher than in the forest zones with higher evapotranspiration, however the shall ow root system of the grasses and the relatively low transpiration demand did not tap into this available water. Gap fraction limited canopy storage considerably which had a direct negative affect on transpiration as it reduced shallow water content (Figur e 5. 6). Figure 5.6. Spatial variability in water balance. Modeled 2012 growing season (Apr 1 to Nov 1) water balance along the transect at 2 - m resolution. Annual transpiration (blue), canopy interception (red), and evaporation (yellow) and precipitati on (blue line). Note the absence of canopy interception and storage in the grass portion from 140 - 160 m. We also compared the effect of model resolution on the estimated water balance. A summary of the spatially integrated annual growing season water balan ce for each resolution scenario is presented in Figure 5. 7 . Reducing the resolution of the spatial parameters did not have a large effect on the modeled water balance, with only 1 - 2 cm of water (~2 - 4%) difference 124 across the four scenarios in any year. The discrepancy was larger in years with more precipitation which leads to more variability in along - line canopy storage and transpiration. Figure 5. 7 . Effect of resolution on water balance. Annual growing season (Apr 1 to Nov 1) average water balance integrated along the line for four different 2D hydrological model resolution scenarios. Stacked bars represent total transpiration, canopy storage, and evaporation, for the four scenarios with decreasing resolution from left to right. Annual growi ng season precipitation is indicated with a blue horizontal bar. 4 . Conclusions In this paper, we use ER data in conjunction with a coupled hydrogeophysical model to infer drivers of spatial and temporal water content heterogeneity in a successional forest setting. An analysis of variance demonstrated a significant relationship between above - ground canopy structure and ER both prior to and during the growing season. We used a forward hydrogeophysical model to test how much of the observed patterns could be explained by 125 incorporating the physically - based and spatially heterogenous vegetation driven processes, and found only small spatial effects due to differences in evapotranspiration and canopy interception with the proportion of coarse root mass and soil t exture accounting for a much larger effect. We then used the hydrological model to compare the along - line heterogeneity in the water balance driven by differences in soil and vegetation properties as well as an analysis of the importance of model resolutio n to appropriately capture site - scale dynamics. Acknowledgments This material is based upon work primarily supported by the US Department of 2015 - 68007 - 23133, with additional support from the Great Lakes Bioenergy Research Center, U.S. Department of Energy, Office of Science, Office of Biological and Environmental Research under Award Numbers DE - SC0018409 and DE - FC02 - 07ER64494; the National Science Foundation Long - term Ecological Research Prog ram (DEB 1832042) at the Kellogg Biological Station, and by Michigan State University AgBioResearch. 126 REFERENCES 127 R EFERENCES GCB Bioenergy : 1 13. http://doi.wiley.com/10.1111/gcbb.12239. - Term Evapotranspiration Rates for Rainfed Corn versus Perennial Bioenergy Crops in a Hydrological Processes 34(3): 810 22. Allen, R. G., L. S. Pereira Irrigation and Drainage Paper No. 56. FAO, Rome, Italy . - Elec Tree physiology 28(Pollen 2007): 1441 48. Petroleum Technology : 54 62. rmination of Soil Hydraulic Parameters and Its Hydrol. Earth Syst. Sci 14(1): 251 70. www.hydrol - earth - syst - sci.net/14/251/2010/. Bass, Benjamin, M. Bayani Cardenas, and Kevi Vadose Zone Journal 16(2): 0. https://dl.sciencesocieties.org/publications/vzj/abstracts/0/0/vzj2016.11.0108. Basso, Bruno, J Italian Journal of Agronomy 4: 677 88. - Dimensional Monitoring of Soil Water Cont ent in a Maize Field Hydrology and Earth System Sciences 17: 595 609. Root Water Uptake Models Considering Dynamic Root Di stributions and Water Uptake Vadose Zone Journal (May). Chung, Sang - Water Resources Research 23(12): 2175 86. Pridmore, and Sacha J. Mooney. 2020. 128 Geoderma 365(September 2019). Bulletin of the American Meteorological Soci ety 84(8): 1013 23. Electrical Resistivity Surveys to Assess Heterogeneity in Soil Moisture Dynamics under Journal of Hydrology 559: 684 97. Dickinson, Robert E., Ann Henderson - Sellers, Cynthia Rosenzweig, and Piers J. Sellers. 1991. Agricultural and Forest Meteorology 54(2 4): 373 88. Diffenbaugh, N.S., Science 341(August): 486 93. Optimization for Conceptual Rainfall - Runoff Mode Water Resources Research 28(4): 1015 31. Vadose Zone Journal (December): 1 29. utrient Flux Following Whole - Forest Science 34(3): 744 68. Field Crops Research 189: 68 74. http://dx.doi.org/10.1016/j.fcr.2016.02.013. - Zone Soil Water in a Journal of Hydrology 523: 47 5 88. http://linkinghub.elsevier.com/retrieve/pii/S0022169415000815. Bulletin of the American Meteorological Society 82(12): 2797 2809. http://journals.ametsoc .org/doi/abs/10.1175/1520 - 0477%282001%29082%3C2797%3AMRWUIH%3E2.3.CO%3B2. Eos Trans. AGU 90(23): 200. - Dimensional Electrical Resistivity Tomography to Monitor Root Vadose Zone Journal 10: 412. 129 Vadose Zone Journal 12(2): 1 12. http://www.scopus.com/inward/record.url?eid=2 - s2.0 - 84878260843&partnerID=40&md5=bfb19e0c7d d27803d78a4fef0cf6c59d. Vadose Zone Journal 11(4). https://www.soils.org/publication s/vzj/abstracts/11/4/vzj2011.0186 (November 6, 2013). - from Equation for Predicting the Hydraulic Soil Science Society of America Journal 44: 892 98. Georgescu, M., D. B. Lobell, and Proceedings of the National Academy of Sciences 108(11): 4307 12. http://www.pnas.org/cgi/doi/10.1073/pnas.1008779108. Gora, Evan M., and Stephen P. Yanoviak. 20 Canadian Journal of Forest Research 45(3): 236 45. - Zone, Trunk, and Moisture J ournal of Experimental Botany 58(4): 839 54. Environmental Research Letters 10(6). Hao, Xinmei, Renduo Zhang, and Alexandr Soil Science 170(3): 167 74. Vadose Zone Journal 17:170040. Temperature Dependence of Electrical Resistivity: Implications for near Surface Geophysical Research Letters 34(18): 1 5. http:/ /www.agu.org/pubs/crossref/2007/2007GL031124.shtml (March 2, 2012). Comparison of Canopy Evapotranspiration for Maize and Two Perennial Grasses Identified as Potential B GCB Bioenergy 2: 157 68. http://doi.wiley.com/10.1111/j.1757 - 1707.2010.01050.x. Water Resources R esearch 46(1): 1 14. http://www.agu.org/pubs/crossref/2010/2008WR007060.shtml (May 4, 2012). - Water - Recharge Rates in the 130 United States Geolgoical Survey Water - Supply Paper 2437 . and Multi - Journal of Hydrology 380(1 2): 62 73. http://linkinghub.elsevier.com/retrieve/pii/S0022169409006702 (Nov ember 9, 2012). Water Resources Research 39(11). Natur e 496: 347 50. http://www.ncbi.nlm.nih.gov/pubmed/23552893. - Use and Climate on Regional Hydrology and Groundwater Recharge - Jayawickreme, Dushmantha H., Remke L. Van Dam, and David W. Hyndman. 2008. 35 Geophysical Research Letters Subsurface Imaging of Vegetation, Climate, and Root - Zone Moisture Interactions . http://www.agu.org/pubs/crossref/2008/2008GL034690.shtml (May 24, 2012). nd - of Plants on Soil Moisture with Time - Geophysics 75(4): WA43 50. Pleistocene Glacial Drift (Michigan, USA): Mass Balances Based on Soil Water Geochimica et Cosmochimica Acta 72(4): 1027 42. Origin of Tunnel Valleys of the Saginaw Lobe of the Laurentide Ice Sheet; Michigan, Boreas 42(2): 442 62. Kuhl, Alexandria S., Anthony D. Kendall, Remke L. Van Dam, and David W. Hyndman. 2018. Vadose Zone Journal 17(170154): 1 13. ty in a Loamy Soil: Identification of the Appropriate Pedo - Vadose Zone Journal 10(3): 1023 33. - Current Geoelectrical Imaging Journal of Applied Geophysics 95: 135 56. ht tp://linkinghub.elsevier.com/retrieve/pii/S0926985113000499 (November 3, 2013). - Storied Soils on the Outwash Plains of Annals of th e American Association of Geographers 106(3): 551 71. 131 Variability in a Temperate Deciduous Forest: Insights from Electrical Resistivity and Environmental Earth Sciences 72(5): 1367 81. http://link.springer.com/10.1007/s12665 - 014 - 3362 - y (June 1, 2015). Maeght, Jean - - and Why Frontiers in plant science 4(August): 299. http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=3741475&tool=pmcentrez&ren dertype=abstract. Mann, J. System Dynamics of Miscanthus × Giganteus and Panicum Virgatum in Response to Bioenergy Research 6(2): 678 87. Mboh, C. M. et al. Inflow Measurements for Topsoil Hydraulic Properties Properties under Constant Head Near Surface Geophysics 10: 413 26. McNeal, B.L., J.D. Oster, and J.T. Hatcher. 1970 Solution Composition Data as an Aid to In - Soil Science 110(6): 10. Irrigated Corn C Water Resources Research 39(5): 1 20. Minsley, Burke J, Bethany L Burton, Scott Ikard, and Michael H Powers. 2011. Journal of Environmental and Engineering Geophysics 16(4): 145 64. Global Change Biology 12(1): 84 96. Moreno, Ziv, Ali Arnon, and Ale - Geophysical Monitoring of Orchard Root Zone Dynamics in Semi - Irrigation Science 33: 303 18. http://dx.doi.org/10.1007/s00271 - 015 - 0467 - 3. rgy GCB Bioenergy 9: 1476 88. International Journal of Geophysics 2013(1 ). - D Electrical Computers & Geosciences 34(12): 1645 54. http://linkinghub.elsevier.com/retrieve/pii/S0098300408001519 (January 9, 2013). Journal of Applied Geophysics 1(5): 318 33. 132 Understanding Options for Agricultural Production , Springer, 41 54. Robertson, G Philip et Science 356(1349): 1 9. Variability in Hydraulic Redistribution with in an Oak Pine Forest from Resistivity Journal of Hydrology 430 431: 69 79. http://linkinghub.elsevier.com/retrieve/pii/S002216941200100X (September 10, 2012). Cycle Advances in Geosciences 18: 43 50. - Dimensional Modelling and Inversion of Dc Resistivity Data Incorporating Topography - Geophysical Journal Int ernational 166(2): 495 505. Annals of Botany 118(4): 555 59. rnative Cellulosic Bioenergy Agriculture, Ecosystems and Environment 216: 344 55. http://dx.doi.org/10.1016/j.agee.2015.10.018. Com puter Program for Estimating Soil Hydraulic Parameters with Hierarchical Journal of Hydrology 251: 163 76. of Soil Hydraulic and Root Dis Vadose Zone Journal . Water Resources Re search 44(7): 1 12. http://doi.wiley.com/10.1029/2007WR006644 (July 9, 2013). - 1D Software Package for Simulating the One - Dimensional Movement of Water, Heat, and Multiple Solutes in Variably - Saturated Module for HYDRUS (2D/3D) - Simulating Two - Dimensional Movement of and Reactions 133 Vadose Zone Journal 15(7): 25. Architecture Differ among Species within PLoS ONE 12(10): 1 22. Singha, K., F. D. Day - Subsurface Processes with Time - Hydrological Processes 29(6): 1549 76. http://doi.wiley.com/10.1002/hyp.10280 (June 2, 2015). Sprunger , Christine D., Steve W. Culman, G. Philip Robertson, and Sieglinde S. Snapp. 2017. Renewable Agriculture and Food Systems : 1 13. ality of Water Percolating below the Root Zone Agricultural Water Management 221(September 2018): 109 19. https://doi.org/10.1016/j.agwat.2019.04.008. ide Transport and Production Water Resources Research 29(2): 487 97. - Layer Vadose Zone Journal 16(2): 0. https://dl.sciencesocieties.org/publications/vzj/abstracts/16/2/vzj2016.08.0072. Using Coupled Hydrological - Thermal - Hydrology and Earth System Sci ences 20: 3477 91. Hydrological Processes 21: 3647 50. mography to Identify Soil - Journal of Hydrology 556: 310 24. https://doi.org/10.1016/j.jhydrol.2017.11.025. Giganteus Pr GCB Bioenergy (2): 180 91. - Dimensional Root Soil Science Society of America Journal 65(4): 1027 37. Warren, Jeffrey Biosphere Models -- The New phytologist 205(1): 59 78. http://www.ncbi.nlm.nih.gov/pubmed/25263989. uctivities in Oil - Society of Petroleum Engineers Journal 8(2). 134 Internal Note . nges in Soil Water for Phenotyping Root Plant and Soil 415(1 2): 407 22. Canadian Journal of Forest Research 36(2): 450 59. - aside Gra ssland to Annual and Perennial Cellulosic Agricultural and Forest Meteorology 182 183: 1 12. http://dx.doi.org/10.1016/j.agrformet.2013.07.015.