QUANTIFYING FINE-SCALE VARIABILITY IN SOIL AVAILABLE WATER WITHIN HOT, DRY DOUGLAS-FIR ECOSYSTEMS. by DOUGLAS JOHN TERPSMA B.Sc. Nat. Res. Sci., Thompson Rivers University, 2016 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE IN ENVIRONMENTAL SCIENCES in the Department of Natural Resource Science Thompson Rivers University Kamloops, British Columbia, Canada October 2018 Thesis examining committee: Thomas Pypker (PhD), Associate Professor and Thesis Supervisor, Department of Natural Resource Sciences Brian Wallace (PhD), Range Soils Ecologist and Adjunct Committee Member, BC Ministry of Forests, Lands, Natural Resource Operations and Rural Development Wendy Gardner (PhD), Associate Professor and Committee Member, Department of Natural Resource Sciences Margaret Schmidt (PhD), External Examiner, Simon Fraser University, Department of Geography ii Thesis Supervisor: Associate Professor Dr. Thomas Pypker ABSTRACT Dry, grassy sites within Douglas-fir forests are difficult to regenerate following harvesting due to harsh climatic conditions and intense inter-species competition for limited moisture resources. The objective of this study was to sample soil properties of a recently harvested opening (16 ha) to provide information about post-harvest soil conditions and their relationship with water holding capacity within these dry ecosystems. Utilizing soil properties with established soil water characteristic equations, I predicted soil water holding capacity (SWHC) across the site and at four depths. We completed a topographic survey for the site using aerial light detection and ranging (LiDAR) technology to create a high-resolution (~1m) digital elevation model (DEM). We statistically compared multiple topographic variables with water retention properties via multiple linear regression and geographically weighted regression to determine what drives soil moisture distribution on finer scales than previously studied. Coarse fragments (CF) had the highest amount of variability on the site and altered SWHC the greatest compared with other measured soil properties, with a 10% increase in CF corresponding with a 4.7 mm decrease in SWHC. Additionally, geographically weighted regression was found to outperform multiple linear regression for interpolating the measured soil properties using principle component derived topographic predictor variables. However, the models only explained roughly one half or less of the variability in all soil properties (most R 2 ≤ 0.50), thereby suggesting that local soil properties be measured to gain accurate representations of any given site prior to conducting logging or site prescription treatments. To maintain sustainable timber resources within dry Douglas-fir forests in a changing climate, an improved understanding of the soil condition prior to regeneration will become increasingly important. keywords: soil, variability, interpolation, water, topography, drought, geostatistics, LiDAR iii TABLE OF CONTENTS ABSTRACT ...................................................................................................................... ii TABLE OF CONTENTS .................................................................................................. iii ACKNOWLEDGMENTS .................................................................................................. v LIST OF FIGURES ......................................................................................................... vi LIST OF TABLES.......................................................................................................... viii LIST OF ACRONYMS ..................................................................................................... xi CHAPTER 1.0: THESIS INTRODUCTION ..................................................................... 1 Natural Disturbance History ......................................................................................... 2 Forest Management Issues ......................................................................................... 4 Soil Water Relationship ............................................................................................... 6 Research Objectives ................................................................................................... 8 LITERATURE CITED ...................................................................................................... 9 CHAPTER 2.0: FINE-SCALE SOIL VARIABILITY AND ITS IMPACT ON WATER AVAILABILITY WITHIN HOT, DRY DOUGLAS-FIR ECOSYSTEMS OF BC‟S SOUTHERN INTERIOR. ............................................................................................... 11 INTRODUCTION .......................................................................................................... 11 MATERIALS AND METHODS ...................................................................................... 14 Site Description ......................................................................................................... 14 Research Design ....................................................................................................... 17 Climate and Moisture Data ........................................................................................ 17 Measured Soil Physical and Chemical Properties ..................................................... 18 Derived Soil Variables ............................................................................................... 19 Statistical Analysis ..................................................................................................... 21 Measured Soil Properties ....................................................................................... 21 Derived Soil Properties ........................................................................................... 22 RESULTS ..................................................................................................................... 24 Climate and Moisture Data ........................................................................................ 24 Measured Soil Properties .......................................................................................... 25 Derived Soil Properties .............................................................................................. 28 Statistical Analysis ..................................................................................................... 32 iv Measured Soil Properties ....................................................................................... 32 Derived Soil Properties ........................................................................................... 36 DISCUSSION................................................................................................................ 42 Soil Forming Factors ................................................................................................. 42 Soil Variability ............................................................................................................ 43 Calculated Soil Properties ......................................................................................... 45 CONCLUSION/SUMMARY ........................................................................................... 48 LITERATURE CITED .................................................................................................... 49 CHAPTER 3.0: TESTING THE APPLICATION OF SPATIAL INTERPOLATION REGRESSION METHODS FOR PREDICTING FINE-SCALE FOREST SOIL VARIABILITY. ............................................................................................................... 53 INTRODUCTION .......................................................................................................... 53 METHODS .................................................................................................................... 58 Site Description ......................................................................................................... 58 Research Design ....................................................................................................... 58 Soil Properties ........................................................................................................... 58 Mapping and Geostatistical Analysis ......................................................................... 59 RESULTS ..................................................................................................................... 62 Soil Properties ........................................................................................................... 62 Mapping and Geostatistical Analysis ......................................................................... 62 DISCUSSION................................................................................................................ 73 Interpolated Soil Variability ........................................................................................ 73 Effect of Increased Kernel Radii on Soil Interpolations .............................................. 77 Evaluation of Methodologies ...................................................................................... 78 CONCLUSION/SUMMARY ........................................................................................... 79 LITERATURE CITED .................................................................................................... 80 CHAPTER 4.0: SIGNIFICANCE AND CONCLUSION .................................................. 87 LITERATURE CITED .................................................................................................... 91 5.0 APPENDIX I ............................................................................................................ 93 v ACKNOWLEDGMENTS I would like to sincerely thank Dr. Tom Pypker for his continuous support, countless edits, outstanding kindness and overall help with constructing my thesis. You were an absolutely fantastic supervisor and I wish you continued success in the future. Additionally, I would like to deeply thank the Ministry of Forests, Lands, Natural Resource Operations and Rural Development (MFLNRORD) for providing the funding which made this project possible, and to Dr. Brian Wallace and Dr. Chuck Bulmer (FLNRORD soil researchers) for all that they have contributed. Thank you, Brian, for your nonstop mentoring and aid with soil sampling, lab methods and data analyses. You have taught me a lot about the profession of soil ecology, and have provided me with enough experience to feel confident entering the workforce. And thank you to Chuck, for numerously taking the time to travel to Kamloops from Vernon to visit with me, and for aiding with all of the mapping aspects of this project. Your expertise in geographic information systems and mapping software was a crucial addition to my work, and for that I am very grateful. The outstanding support that all three of you continue to show me with respect to my Master‟s project and my growing career has far exceeded my expectations when I first entered the program, and has not gone unnoticed. Finally, I would like to thank Dr. Wendy Gardner for taking the time to sit on my committee and for offering her support, edits and kind conversations throughout the last two years. I could not have asked for a better committee to work with, and I am very grateful for everything the four of you have done. vi LIST OF FIGURES Figure 2.1: Aerial map of the harvested Isobel lake study site showing soil pit and soil water sensor location.……………………………………………………………. ….16 Figure 2.2: Total volumetric water content (θ) (in mm) recorded daily from the Isobel Lake study site from April of 2016 until September of 2017. Based on the data within this figure the site loses roughly 1.2 mm per day of total soil water………………....25 Figure 2.3: Scatterplot of the linear regression model comparing the measured BD values (obtained from field measurements) and predicted BD values (from Saxton and Rawls (2006) equations). BD = fine fraction bulk density……………………..…28 Figure 2.4: Soil profile of pit 7, face 'A' (2.3A): representative of a Orthic Gray Luvisol also found in pits 2, 5 and 17 at the Isobel site. Soil Profile of pit 12, face „A‟ (2.3B): representative of an Orthic Eutric Brunisol discovered in the remaining 16 soil pits.29 Figure 2.5: Averaged concentrations (%) of soil hydraulic variables within each sampled depth range. PWP = permanent wilting point, FC = field capacity and AWC = available water capacity. Error bars represent standard error of the mean. Letters represent significant statistical differences (p = 0.05) using Tukey multiple comparison of means test………………………………………………...………………31 Figure 2.6: Soil volumetric water content within the upper 30 cm soil layer in relation to the calculated values of site FC, PWP and SWHC in mm. SWHC (in mm) represents the fraction of soil water available for plant absorption above the soil wilting point. FC = field capacity, PWP = permanent wilting point, AWC = available water capacity, and SWHC = soil water holding capacity………………..……………32 Figure 2.7: Estimated recommended sample size (number of soil pits to dig) to identify the level of variability in AWC at the site (within the top 30 cm of soil), using only coarse fragments (% vol.) as the explanatory variable in a linear regression model. The „pwr.f2.test‟ was utilized within the „pwr‟ package in „R‟ for determining the sample size with a significant level of 0.05 and 0.10, f2 = R2/(1-R2) (R2 determined by testing the linear model with AWC as the dependent variable, and CF as the independent), and u = number of variables used in the linear regression model (1; coarse fragments). Power was allowed to fluctuate to determine sample size at different levels of certainty…………………………………………………….....39 vii Figure 2.8: Effect of fluctuating soil variables used in the soil water characteristics equation (Saxton and Rawls 2006) on SWHC (mm). While one variable is fluctuating, the rest were held constant at their average value from the top 30 cm of soil measured at the site. Figure (2.7A) = Fluctuating coarse fragments (% w.). Figure (2.7B) = Fluctuating soil organic matter (% v.). Figure (2.7C) = Fluctuating sand content (dec. % v.). Figure (2.7D) = Fluctuating clay content (dec. % v.). Figure (2.7E) = Fluctuating density factor. SWHC = soil water holding capacity…..……….41 Figure 3.1: Fluctuating kernel radius of topographic grid layers and their associated R2 values with respect to different soil properties. A total of 12 grid layers were used in the regression models for each individual kernel radius and were used as predictor grids for explaining soil patterns within the upper 15 cm of soil. BD = fine fraction soil bulk density…………………………………………………………..………63 Figure 3.2: Comparison of R2 values for MLR and GWR methods using the 7 PCA grids as predictor variables for interpolating various soil properties within the upper 15 cm of soil. SOM = soil organic matter, CF = coarse fragments, BD = fine fraction bulk density and SWHC = soil water holding capacity…………………………..…….70 Figure 3.3: Leave-one-out cross validation results comparing the multiple linear regression (MLR) and geographically weighted regression (GWR) models for interpolating the values of fine fraction bulk density (g cm-3). A = depth 1 (0-15 cm), B = depth 2 (15-30 cm), C = depth 3 (30-60 cm) and D = depth 4 (60-100 cm). Predictor variables=7 PCA grids. Tested soil property = fine fraction bulk density (BD)(g cm-3)………………………………………………………………………………...71 Figure 3.4: Output map of the interpolated values of fine fraction bulk density (BD) (0-15 cm) using Geographically Weighted Regression (GWR) and the 7 topographic grid layers produced through principle component analysis (PCA) as independent, predictor variables……………………………………….………………………………...72 viii LIST OF TABLES Table 2.1: Averaged soil properties measured throughout the entire study site (n=60 for variables of clay, silt, sand, pH, total C, total N, and C:N for each depth range: 20 pits X 3 faces/pit at all 4 depths. n=20 for variables of CF, BD and porosity: 1 measurement per depth within all 20 pits). Values within brackets represent one standard deviation from the mean. Total C = total carbon, Total N = total nitrogen, C:N = carbon to nitrogen ratio, CF = coarse fragments (>2 mm), and BD = fine fraction bulk density. Superscripted letters represent significant differences between values within different depth ranges (assessed using the Tukey multiple comparison of means test)………………………………………………………………………..…….27 Table 2.2: Average values (in mm of water) for the calculated soil hydraulic variables of saturation, field capacity (FC), permanent wilting point (PWP), available water capacity (AWC) and soil water holding capacity (SWHC) within the top 3 measured depths throughout the entire study site (n=59). Numbers in brackets represent standard error of the mean………………………………………...……………………..30 Table 2.3: Soil physical and chemical property correlation matrix using data collected from the Isobel lake study site. The included variables are as follows: BD= fine fraction bulk density; CF= lab determined coarse fragments; pH= soil pH; C:N= carbon to nitrogen ratio; N = total nitrogen; C = total carbon; Sand, Silt and Clay= textural proportions of soil sample; Depth= soil layer depth; and Soil Type = Brunisol or Luvisol. Bold-texted and highlighted cells represent correlations at the 0.05 significance level, and bold-texted cells represent correlations at the 0.10 significance level.......................................................................................................33 Table 2.4: Coefficient of variation (CV: relative standard error) percentage of each measured soil property explained by the different levels of sampling. The average between pits shows the average variation found between each soil face of all 20 soil pits, and within pits explains the average variation found within a single pit. CV = standard deviation divided by the mean. CF = coarse fragments, BD = fine fraction bulk density, C:N = carbon to nitrogen ratio………………………………………....….35 Table 2.5: Correlation matrix demonstrating the relationships between the measured soil variables and the predicted soil hydraulic variables. Soil Type = Brunisol or Luvisol, C = Total Carbon, N = Total Nitrogen, CF = Coarse Fragments, PWP = Permanent Wilting Point, FC = Field Capacity, AWC = Available Water Capacity, and SWHC = Soil Water Holding Capacity. Bold-texted and highlighted cells represent correlations at the 0.05 significance level, and bold-texted cells represent correlations at the 0.10 significance level……………………………………………....38 ix Table 3.1: Descriptive statistics of the 12 topographic raster grid layers produced in SagaGIS using the LiDAR-derived digital elevation model for the Isobel Lake study site.…………………………………………………….……………………………………60 Table 3.2: Average values (in % vol. of water) for the calculated soil hydraulic variables of field capacity (FC), permanent wilting point (PWP), available water capacity (AWC) and soil water holding capacity (SWHC) at 4 different depths throughout the entire study site (n=60 for each depth). Numbers in brackets represent standard error of the mean. Superscript letters represent significant statistical differences (p = 0.05) between mean values using Tukey multiple comparison of means test.....…………………………………………………………….62 Table 3.3: The number of times (out of a total of 15 tested soil variables) each topographic grid layer showed a significant correlation at the p = 0.10 level with a soil property. DEM Elev = site elevation, CI = catchment index, LS = length-slope factor, TPI = topographic position index, PO = positive openness, CA = catchment area, = CS = catchment slope, MCA = modified catchment area, TWI = topographic wetness index, MRVBF = multi-resolution of valley bottom flatness, VD = valley depth..…………………………………………...…..…….………………...………..……65 Table 3.4: Correlation matrix showing the relationship between measured soil properties and derived hydraulic properties within different depth ranges of the soil profile with topographic grid layers produced in SagaGIS. Bold text and highlighted cells represent significant relationships at the 0.10 level. DEM Elev = site elevation, CI = catchment index, LS = length-slope factor, TPI = topographic position index, PO = positive openness, CA = catchment area, = CS = catchment slope, MCA = modified catchment area, TWI = topographic wetness index, MRVBF = multiresolution of valley bottom flatness, VD = valley depth, C = total carbon, N = total nitrogen, C stock = carbon stock, CF = coarse fragments, BD = fine fraction bulk density, PWP = permanent wilting point, FC = field capacity, AWC = available water capacity, and SWHC = soil water holding capacity………………........................................................................................……65 Table 3.5: Principle component analysis loadings for 7 components identified to explain over 90% of the landscape dataset variance. The 7 principle components were derived from a combination of the 12 individual topographic grid layers of DEM Elev = site elevation, CI = catchment index, LS = length-slope factor, TPI = topographic position index, PO = positive openness, CA = catchment area, = CS = catchment slope, MCA = modified catchment area, TWI = topographic wetness index, MRVBF = multi-resolution of valley bottom flatness, and VD = valley depth..68 x Table 3.6: Resulting coefficients of determination for multiple linear regression (MLR) models using individual grid layers as predictor variables (12 total), and principle components as predictor variables (7 total), for predicting different soil properties within 4 different depths. Total C = total carbon, Total N = total nitrogen, C Stock = carbon stock, CF = coarse fragments, BD = fine fraction bulk density, AWC = available water capacity and SWHC = soil water holding capacity………………..69 Table 5.1: Description and reference for the 12 selected topographic grid layers used as covariate predictor grids for soil spatial predictions.….………………....…..93 Table 5.2: Comparison of MLR and GWR regression models using the 7 created PCA grids as predictor variables for various soil properties at all 4 measured depths. Average R2 values are created by running 10 tests with training datasets using leave-one-out cross validation. N/A values means no significant predictors were identified at the p = 0.10 level.………………………………………...…………………94 xi LIST OF ACRONYMS 95% CI = 95% Confidence Interval AIC = Akaike Information Criterion AWC = Available Water Capacity BD = Soil Fine Fraction Bulk Density C = Carbon C:N = Carbon to Nitrogen Ratio CA = Catchment Area CF = Coarse Fragments CI = Convergence Index CL = Clay Loam CMD = Hargreave‟s Climatic Moisture Deficit CS = Catchment Slope CTI = Compound Topographic Index DEM = Digital Elevation Model DF = Density Factor DSM = Digital Soil Map Eref = Hargreave‟s Reference Evaporation FC = Field Capacity GIS = Geographical Information Systems GPS = Global Positioning Systems IDF = Interior Douglas-fir L = Loam LiDAR = Aerial Light Detection and Ranging LS = Length-Slope Factor LTSP = Long Term Soil Productivity MAP = Mean Annual Precipitation MAT = Mean Annual Temperature MBD = Maximum Bulk Density MCA = Modified Catchment Area MRVBF = Multi-Resolution of Valley Bottom Flatness N = Nitrogen NDT4 = Natural Disturbance Type 4 PO = Positive Openness PTF = Pedotransfer Function PWP = Permanent Wilting Point RBD = Relative Bulk Density RH = Relative Humidity SCL = Sandy Clay Loam SL = Sandy Loam SOM = Soil Organic Matter SPI = Stream Power Index SWHC = Soil Water Holding Capacity TPI = Topographic Position Index TWI = Topographic Wetness Index VD = Valley Depth VIF = Variance Inflation Factor VPD = Vapour Pressure Deficit 1 CHAPTER 1.0: THESIS INTRODUCTION Forests of the Interior Douglas-fir (IDF) biogeoclimatic zone provide crucial natural resources within British Columbia (B.C.) for multiple industries including range, wildlife, recreation and forestry. Timber resources of these forest types have been extensively managed for over a century throughout the province, especially where access is favourable (ie. valley bottoms) and timber resources are easily extractable (Heineman et al. 2003). Within the Kamloops forest district of the Southern Interior Forest Region, the driest IDF variants have received a considerable amount of attention throughout the last few decades as healthy and sustainable forest management has been a growing concern along with trends in increasingly harsh climatic conditions (Heineman et al. 2003). The driest variants of the IDF surrounding Kamloops include the IDFxh1, IDFxh2 and the IDFdk1 (Lloyd at el. 1990). The IDFxh1 is defined as the Okanagan very dry hot IDF variant, the IDFxh2 is defined as the Thompson very dry hot IDF variant, and the IDFdk1 is defined as the Thompson dry cool IDF variant within Land Management Handbook number 23 (Lloyd et al. 1990). Within these dry climates, issues surrounding postharvest regeneration, timber growth and yield, wildlife resources, pest management, and cattle grazing must all be considered when producing land management guidelines as multiple stakeholders own value within these areas, especially in valley bottoms which are in close proximity with urban development (Vyse et al. 1998). Among the driest ecosystem variants within the IDF biogeoclimatic zone is the IDFxh2 subzone (Lloyd et al. 1990), which discontinuously surrounds the city of Kamloops, BC. Within the IDFxh2, soil types commonly found include Eutric, Melanic and Dystric Brunisols as well as Gray Luvisols (Lloyd et al. 1990). The degree of soil profile development dictates its classification within the Canadian system (Soil Classification Working Group 1998), while the sufficient redistribution of clay from the upper surface horizons to the lower, sub-horizons is a commonly used measure for distinguishing between Luvisolic and lesser developed Brunisolic soils 2 within these forest types. Soils of the Brunisolic order found within semi-arid ecosystems of BC are often restricted by climatic factors including a lack of soil water to facilitate weathering (Valentine et al. 1978). On the other hand, Luvisolic soils of these ecosystems have experienced enough chemical and physical weathering to allow for sufficient amounts of clay to leach from the soil surface into the lower soil layers. Therefore, soil weathering of Luvisolic soils within semi-arid forests is often attributed to sufficient levels of soil water over time for clay leaching to occur, whereas the Brunisolic soils in similar ecosystems likely have not received the same levels of water throughout their lifetime (Valentine et al. 1978). The redistribution of clay within Luvisolic soils is an important feature to forest managers as the clay can accumulate and, if sufficient compaction occurs overtop, root and water penetration can become restricted (Bengough and Mullins 1990). Additionally, these soil layers can become very wet in the spring, potentially leading to slope stability issues when working with heavy machinery (Valentine et al. 1978). A general lack of certainty regarding the spatial distribution of these soil types found within the IDFxh2 has been an issue for forest managers in the past when writing site prescriptions. Therefore, exploring the fine-scale variability of physical and chemical soil properties within these forest types can also help uncover patterns in the characteristics which dictate soil order. Natural Disturbance History The specific site examined within this study is characteristic of an Interior Douglas-fir, very dry, hot (IDFxh2) variant. This forest type is included in the natural disturbance type 4 (NDT 4) outlined within the B.C. Biodiversity Guidebook (British Columbia Ministry of Forests and B.C. Ministry of Environment, Lands and Parks 1995). This particular NDT type is defined by grasslands, shrublands and forested communities which experienced frequent (roughly every 4-50 years), low-intensity fires, with stand replacement crown fires occurring approximately every 150-250 years throughout their history. The IDF biogeoclimatic zones often experienced 3 periodic surface fires that consumed low-lying woody fuels, increased some shrub and herb communities while selecting against others, thinned younger stands, and maintained canopy height of live tree crowns (British Columbia Ministry of Forests and B.C. Ministry of Environment, Lands and Parks 1995). This historic fire regime maintained vegetation species composition and forest stand structure as well as regulated the accumulation of flammable woody debris. However, several decades of fire exclusion within these ecosystems has caused many Douglas-fir stands to fill in with young conifers resulting in fuel build-up which increases the chance of crown instead of surface fires. This also leads to an overall loss of understory forage for cattle grazing, and has the potential to promote insect and disease damage which is of critical concern to current forest managers (British Columbia Ministry of Forests and B.C. Ministry of Environment, Lands and Parks 1995). Common mismanagement of the dry IDF forest types within the NDT4 has significantly affected the natural biodiversity. Historic areas of mismanagement include unregulated livestock grazing during early European settlement which lead to the establishment and spread of invasive species, as well as the introduction of non-native forage plants (British Columbia Ministry of Forests and B.C. Ministry of Environment, Lands and Parks 1995). In addition, human influence has destroyed shrub and tree cover in some areas which greatly affected wildlife species that rely on them for forage, protective cover and breeding purposes (British Columbia Ministry of Forests and B.C. Ministry of Environment, Lands and Parks 1995). Areas that are occupied by closed forest canopies can be used as transitory rangeland, as timber harvesting and stand-initiating fires create sustainable amounts of forage for short periods of time. Because these sites are important for both timber extraction and livestock grazing interchangeably, both must be considered carefully when performing any kind of recommendations post-disturbance (British Columbia Ministry of Forests and B.C. Ministry of Environment, Lands and Parks 1995). Throughout recent decades, the widespread use of uniform stand-level partial cutting (as opposed to clearcutting) in these particular forest types has been the accepted practice for forest managers (Vyse et al. 1998). Although this harvesting 4 technique appears to have eliminated public fears about logging practices, an increasing amount of doubt has been expressed about the continuous use of the partial cutting within the dry-belt Douglas-fir forests. For example, Vyse and others (2006) found that aggregated tree retention during harvest operations, as opposed to evenly disbursed retention (which is the historically preferred cutting method for dry Douglas-fir forests) increased seedling performance as a result of increased light availability and reduced competition for water from large, overstory trees. Because large trees require such a high demand for moisture and light, they often suppress establishing seedlings for many years through intense competition for limited resources (Simpson 2000). However, they also noted that testing over a wider range of environmental conditions and over many years would be necessary to determine the overall outcome of regenerating seedlings as different results may arise if planting occurred during a dry year, or on coarser textured soils (Vyse et al. 2006). Forest Management Issues In general, dry, grassy sites within the dry Douglas-fir ecosystems of B.C. are considered difficult to regenerate following harvesting. Harsh climatic conditions including severe summer drought, winter temperature extremes as well as summer frosts can negatively affect a seedlings chance at survival during its crucial early years of growth and establishment (Heineman et al. 2003). These forests receive less than half of the annual precipitation during the growing season, creating conditions of intense inter-species competition for limited soil moisture resources by the occupying vegetation. The abundant presence of pinegrass (Calamagrostis rubescens) in the understory also intensifies the moisture deficit as it is an efficient competitor for soil water (Nicholson 1989). It has been noted that natural regeneration is unlikely to be successful on these sites for a variety of reasons including local climate and low levels of seed production, while the success of planted seedlings (especially of Douglas-fir) is also typically low on drier IDF sites (Heineman et al. 2003). Therefore, in order to successfully manage the timber 5 resources within these forest types and to create sustainable timber yield for future generations, forest managers must take into consideration the probable water limitations for seedling growth and come up with alternative approaches to increase planting success. Researchers have previously discovered methods for improving seedling microclimate and maintaining soil water resources later into the growing season within dry Douglas-fir variants (IDFxh2) (Black and Mitchell 1990a; Fleming et al. 1998). Treatments including forest floor scalping, soil ripping and herbicide applications have all successfully conserved soil water when compared to undisturbed control plots. Additionally, within the IDFdk, low soil water levels were found to reduce seedling growth within the controls and were not a limiting factor to growth within the treated areas (Black and Mitchell 1990a). Additional studies within similar ecosystems have demonstrated the effectiveness of site preparation techniques which aim at improving the seedling microclimate and increase the early establishment of seedlings. Removing surface organic layers, altering soil porosity/density and controlling competing vegetation have all been found to influence local radiation, heat conduction as well as vapour diffusion (reducing seedling frost damage) on clearcut sites (Fleming et al. 1998). By removing the surface organic horizons, growing season soil temperatures were improved, and spring and summer frost damage was reduced. Furthermore, controlling vegetation on the grass-dominated sites reduced evapotranspiration of the competing species and increased overall soil water availability, especially in lower elevation sites which often receive less precipitation (Fleming et al. 1998). Although these prescription techniques have been shown to increase the immediate microclimate for seedling establishment and early success, the long-term effects of soil disturbance (through both harvest operations as well as site prescription methods) on timber productivity over the long term is still an unsolved mystery as long-term studies are costly and, of course, time consuming. One approach to help mitigate the residual effects of soil disturbance within the dry Douglas-fir forests of BC‟s interior is the establishment and monitoring of 6 Long Term Soil Productivity (LTSP) sites (Powers et al. 2005). Within the dry IDF zones of the Kamloops forest district there are currently three individual sites (namely O‟Connor Lake, Black Pines and Dairy Creek; all within the IDFdk2) that have been monitored continuously for the past two decades. The site prescription techniques discussed previously were implemented in different treatment plots postharvest within each site, and compared with the control plots to determine the effects of different soil disturbances on tree productivity. The overall objectives of the LTSP project are to evaluate the effects of organic matter removal and soil compaction on natural soil processes and how they may alter site productivity, as well as to develop tools for monitoring soil progress and sustainable forest management (Reid et al. 2015). A main finding from these studies is that the prescription methods are highly site specific, and the results appear to be dependent on the existing soil properties pre-harvest (Reid et al. 2015). This means that one method applied to one site may have completely different results than another site based on the soil and climatic conditions prior to disturbance and seedling regeneration. Because of this, intensive soil property studies hold considerable merit for answering the questions of postharvest soil conditions prior to planting within the fragile forests of the dry IDF. Extensive knowledge of the soils in which seedlings will be competing for resources is becoming increasingly important to ensure the sustainable use of these forest types for multiple different industries into the future. Soil Water Relationship Physical and chemical soil properties influence plant-available soil moisture, which ultimately limits seedling survival following clearcut harvesting in the dry biogeoclimatic subzones of the southern interior of B.C. (Black and Mitchell 1990a; Fleming et al. 1996). It is well-documented that specific soil characteristics (such as texture, bulk density, carbon concentration etc.) strongly control soil water dynamics (Vereecken et al. 1989; Lin et al. 2006), however the spatial distribution, and the degree of soil variability within a site is still relatively unknown throughout these 7 ecosystems. As logging operations in British Columbia continue to progress from beetle infected, pine-dominated stands, into the dry forested regions of the IDF zone (Klenner and Sullivan 2009), management practices and site prescription methods should be altered accordingly based on dissimilar soil and site characteristics in different climatic subzones. If harvested sites within the dry IDF are to be replaced by conventional planting procedures, then competition for moisture by seedlings and understory vegetation will be of critical concern. One approach to a better understanding of root zone parameters across a landscape is to measure the variability of soil physical and chemical properties at known sampled locations, and interpolate them to predict values in unsampled locations to gather a description of the entire study site (Ettema and Wardle 2002). Knowledge of the variability in soil physical and chemical properties can ultimately be used to predict available soil moisture (Kutílek 2004), which will have a major effect on the survival rate of newly planted seedlings following clearcut disturbances in dry, Douglas-fir dominated forests (Black and Mitchell 1990b; Fleming et al. 1994). By combining intensive soil sampling on the ground with spatial mapping techniques, a description of the soils occupying the tree rooting zone can be revealed. The combination of randomly spaced soil properties throughout the landscape along with spatial statistics for data interpolation can demonstrate the degree of variability that occurs in a single-order catchment, and can serve as a baseline dataset which can drive future soil models within this particular forest type. The final dataset of soil properties can be seen as a starting point for the comparison of future timber growth and yield on the study area and in similar ecosystems, as well as future soil disturbance patterns and how soil properties may change over long periods of time. The future of timber harvesting, as well as grazing opportunities within the dry forests of the IDF is uncertain. Providing a detailed understanding of the soil properties which influence plant growth and success will prove to be an asset to forest managers applying logging and regeneration strategies within the IDFxh2. By gaining a full understanding of soil parameters post-harvest, forest managers can begin to uncover the best possible procedures for increasing a seedlings chance at initial survival, as well as continuous growth into the future. 8 Research Objectives Despite a wealth of literature explaining the relationship between soil physical properties and soil hydrological function, a lack of knowledge still remains regarding the relationship between topography and soil variability and its collective effect on the variability of water content and site productivity. The overall goal of this project is to characterize the variability in soil properties at a fine-scale (~1 m) throughout a recently harvested site by utilizing field sampling and remote sensing techniques. I hypothesize that a significant relationship exists between topography and soil variability on a fine scale. In addition, I believe that this relationship can be used to accurately map the variability in soil physical properties in relation to topographic location to help determine soil water availability throughout the landscape. My thesis is separated into two stand-alone chapters. Chapter II involves the field sampling and analysis of variability in soil properties at the site, as well as the use of measured soil properties as well as hydraulic properties in relation to topographical influence. The general goals of each chapter include:  Chapter II – To assess the degree of soil variability at a recently harvested clearcut near Isobel Lake (approximately 15 km North of Kamloops, BC) by intensively sampling the landscape, and to use the measured soil data to predict soil available water. Statistical comparisons of soil available water and measured soil properties will be evaluated to determine which soil properties have the largest influence on changing the soil water potential.  Chapter III – To apply a topographically based approach for characterizing soil conditions at the site based on LiDAR topographic survey data in conjunction with soil parameters. Statistical comparisons of the collective influence of spatial soil variability and multiple topographic features on a finescale will be analyzed to determine which topographic variables significantly alter soil hydraulic parameters. 9 LITERATURE CITED Bengough, A.G. and Mullins, C.E. 1990. Mechanical impedance to root growth: a review of experimental techniques and root growth responses. Journal of Soil Science. 41: 341-358. Black, A. and Mitchell, B. 1990a. Effects of site preparation treatments on the soil moisture regime in IDFdk, MSxk and ESSFxc clearcuts - Project 3.02. Research Memo. Vancouver, BC: Ministry of Forests, Lands and Natural Resource Operations. Black, A. and Mitchell, B. 1990b. IDF zone seedling survival and growth after site preparation treatments - Project 3.02. Research Memo. Vancouver, BC: Ministry of Forests, Lands and Natural Resource Operations. British Columbia Ministry of Forests and B.C. Ministry of Environment, Lands and Parks. 1995. Biodiversity guidebook. Forest Practices Code, Victoria, B.C. Ettema, C.H. and Wardle, D.A. 2002. Spatial soil ecology. Trends in Ecology and Evolution. 17(4): 177–83. doi:10.1016/S0169-5347(02)02496-5. Fleming, R.L. Black, A. and Eldridge, N.R. 1994. Effects of site preparation on root zone soil water regimes in high-elevation forest clearcuts. Forest Ecology and Management. 68: 173–88. Fleming, R.L., Black, T.A., and Adams, R.S. 1996. Site preparation effects on Douglas-fir and lodgepole pine water relations following planting in a pinegrass-dominated clearcut. Forest Ecology and Management. 83: 47–60. Fleming, R.L., Black, T.A., Adams, R.S. and Stathers, R.J. 1998. Silvicultural treatments, microclimatic conditions and seedling response in Southern Interior clearcuts. Can. J. Soil Sci. 78: 115 – 126. Heineman, J.L., Hope, G.D., Simard, S.W., Vyse, A., Lloyd, D.L. and Miège, D.J. 2003. The effects of site preparation and harvesting practices on planted seedling productivity and microenvironment in southern interior dry, grassy IDF forests. Res. Br., B.C. Min. For., Victoria, B.C. Tech. Rep. 009. Hope, G. 2007. Changes in soil properties, tree growth, and nutrition over a period of 10 years after stump removal and scarification on moderately coarse soils in interior British Columbia. Forest Ecology and Management. 242: 625 - 635. Klenner, W. and Sullivan, T.P. 2009. Partial and clearcut harvesting of dry Douglasfir forests: implications for small mammal communities. Forest Ecology and Management. 257: 1078–86. doi:10.1016/j.foreco.2008.11.012. Kutílek, M. 2004. Soil hydraulic properties as related to soil structure. Soil and Tillage Research. 79(2): 175–84. doi:10.1016/j.still.2004.07.006. 10 Lin, H.S., Kogelmann, W., Walker, C. and Bruns, M.A. 2006. Soil moisture patterns in a forested catchment: a hydropedological perspective. Geoderma. 131(3– 4): 345–68. doi:10.1016/j.geoderma.2005.03.013. Lloyd, D., Angove, K., Hope, G. and Thompson, C. 1990. A guide to site identification and interpretation for the Kamloops Forest Region. B.C. Min. For., Res. Br., Victoria, B.C., Land Manage. Handb. No. 23. http://www.for.gov.bc.ca/hfd/pubs/docs/lmh/lmh23.htm Nicholson, A. 1989. Water relations, survival and growth of Douglas-fir seedlings at a pinegrass dominated site in south-central British Columbia. For. Can. and B.C. Min. For., Victoria, B.C. FRDA Memo 121. Powers, R.F., Scott, D.A., Sanchez, F.G., Voldseth, R.A., Page-Dumroese, D., Elioff, J.D., and Stone, D.M. 2005. The North American Long-Term Soil Productivity experiment: findings from the first decade of research. Forest Ecology and Management. 220(1-3): 31-50. Reid, A.M., Chapman, W.K., Kranabetter, J.M., and Prescott, C.E. 2015. Response of lodgepole pine health to soil disturbance treatments in British Columbia, Canada. Can. J. For. Res. 45: 1045 – 1055. Simpson, D.G. 2000. Water use of interior douglas-fir. Canadian Journal of Forest Research-Revue Canadienne De Recherche Forestiere 30:534-547. Vereecken, H., Maes, J., Feyen, J., and Darius, P. 1989. Estimating the soil moisture retention characteristic from texture, bulk density, and carbon content. Soil Science. 148 (6): 389–403. doi:10.1097/00010694-19891200000001. Vyse, A., Hollstedt, C., and Huggard, D. 1998. Managing the Dry Douglas fir Forests of the Southern Interior: Workshop Proceedings. April 29–30, 1997. Kamloops, British Columbia, Canada. Res. Br., B.C. Min. For., Victoria, B.C. Work Pap. 34/1998 Vyse, A., Ferguson, C., Simard, S.W., Kano, T., and Puttonen, P. 2006. Growth of douglas-fir, lodgepole pine, and ponderosa pine seedlings underplanted in a partially-cut, dry Douglas-fir stand in south-central British Columbia. For. Chron. 82:723-732. 11 CHAPTER 2.0: FINE-SCALE SOIL VARIABILITY AND ITS IMPACT ON WATER AVAILABILITY WITHIN HOT, DRY DOUGLAS-FIR ECOSYSTEMS OF BC’S SOUTHERN INTERIOR. INTRODUCTION The importance of plant available soil water for the successful growth and productivity of dry-forest ecosystems has gained an increasing amount of attention as researchers predict an increase in global temperatures and summer drought conditions in the future (Kirchen et al. 2017). By inhibiting photosynthesis and transpiration, soil water deficits during the growing season ultimately affect plant growth and can cause mortality in severe conditions, leading to an overall decline in productivity. Furthermore, the commercially valuable stands within interior Douglasfir forests are affected by water availability (Simpson 2000), and therefore may have increased sensitivity to more frequent summer droughts. The continuous monitoring of soil water is not only useful for projecting plant productivity of a site, but can also improve large wildfire prediction (Krueger et al. 2017). Therefore, the availability of soil water data through time, as well as its apparent variability across the landscape, will have major management implications for both forest stand production and wildfire safety protocols. Knowledge about the heterogeneity of fine-scale (m) soil properties provides crucial insight into soil pedogenesis and vegetative production. Forest growth and tree productivity, especially during summer drought, are largely dependent on local soil water availability at a site (Mathys et al. 2014; Kirchen et al. 2017), which is ultimately driven by variation in soil physical properties and organic matter (Vereecken et al. 1989; Lin et al. 2006). However, the effects of fine-scale soil variability in relation to soil water dynamics are rarely explored within forested hillsides with relatively uniform topography. Due to the high variability in tree seedling mortality within a single harvested site (Heineman et al. 2003), it is believed 12 that soil water distribution throughout any given site within dry-forest ecosystems is likely not uniform, and therefore varies from location to location. In general, dry, grassy sites within the dry Douglas-fir ecosystems of B.C. are considered difficult to regenerate following harvesting due to harsh climatic conditions and intense inter-species competition (Heineman et al. 2003). A combination of drought conditions during the growing season, extreme winter temperatures and occasional summer frost events, all in addition to competition for limited soil water reserves with efficient grass species make it difficult for Douglas-fir seedlings to establish successfully. Furthermore, the success of newly planted seedlings often exhibits considerable short-range variation on most sites within these ecosystems (Heinemann et al. 2003; Fleming et al. 1998). Therefore, to ensure sustainable timber harvests within these forests into the future, updated soil water information needs to be measured at finer scales than previously studied to refine existing soil maps (Mathys et al. 2014). Direct measurements of soil hydraulic parameters are time-consuming (Gupta and Larson 1979; Wösten et al. 2001), and therefore costly, when compared with measurements of soil texture, bulk density and organic matter. Because of the overall difficulty in collecting direct measurements of soil hydraulic characteristics, the use of pedotransfer functions (PTF‟s) for predicting hydraulic parameters has become increasingly popular in soil mapping (Wösten et al. 2001), as well as soil hydrological studies (Gupta and Larson 1979; Vereecken et al. 1989; Wagner et al. 2001; Schaap et al. 2001). While decreasing the sampling time to maintain the practicality of soil mapping methods in forest management situations, PTF‟s are commonly used to evaluate soil available water capacity (AWC), which is defined as the difference in soil water content between field capacity (FC) and permanent wilting point (PWP) (Seneviratne et al. 2010). Above the field capacity, the soil pores cannot resist water loss from gravitational drainage, and below the wilting point, the remaining water in the soil is held by the soil matrix very tightly and is therefore inaccessible to plant roots (Seneviratne et al. 2010; Hillel 1998). Frequently cited equations outlined by Saxton and Rawls (1986; 2006) have shown to be successful 13 for producing estimates of soil water parameters using easy to measure soil properties (Haghverdi et al. 2015). It is widely recognized that soil texture (sand, silt and clay contents) has a major influence on soil water characteristics in both forest and agricultural settings, however, additional variables including soil organic matter (SOM), soil bulk density and coarse fragment content have also shown to have important effects on hydraulic simulations (Saxton and Rawls 2006). The soil water characteristics equations include these variables as input parameters, and produce output predictions of matric bulk density as well as volumetric soil water at different matric potentials of field capacity and wilting point (33 kPa and 1500 kPa, respectively) corresponding to AWC (Saxton and Willey 2006). The overarching goal of this study is to better understand forest soil variability as it affects water holding characteristics that are important for forest regeneration success and productivity within semi-arid ecosystems of BC‟s southern interior. Improper soil characterization of these landscapes can have detrimental effects on the long-term utilization of the IDFxh2 for both timber and forage resources into the future, and can negatively affect local and provincial economies. Recent droughts within BC, along with increasing average summer temperatures throughout the growing season, add to the negative impacts of water stress within these ecosystems and will be ever more difficult to control and manage for in a changing climate. Furthermore, increases in temperature exacerbate water stress of Douglasfir trees by increasing both soil and atmospheric water deficits (Restaino et al. 2016). As one of the most economically and ecologically important species within BC, the successful regeneration of Douglas-fir forests within the province will become increasingly crucial to ensure a profitable forest industry into future generations. By carefully and intensively measuring in situ soil properties, an improved characterization of soil patterns found at the site will reveal the apparent fine-scale variability and can provide knowledge on what ultimately drives soil water distribution within the IDF. Here, I attempt to quantify forest soil variability on a fine-scale across a uniform hillside, where little variation in soil-forming factors exists. Despite a wealth 14 of information explaining the influence that soil properties have on soil available water (Saxton and Rawls 2006; Wösten et al. 2001; Gupta and Larson 1979), and its relationship with tree productivity within these particular ecosystems (Fleming et al. 1998; Vyse et al. 2006), the degree of fine-scale variability that occurs within a single forested catchment has yet to be explored within these forest types. In addition, the statistical modeling of soil properties in relation to one another, as well as with the derived soil hydraulic characteristics is a key aspect of this research as it provides insight on what properties ultimately drive soil water patterns within dry, Douglas-fir ecosystems. This research aims to provide an extensive fundamental dataset to help describe forest soil variability to use in the production of future environmental models within these particular forest types. The specific objectives of this research are (i) to assess the degree of soil variability within a single-order catchment by intensively sampling the harvested landscape, and classifying the soil orders based on measured soil data, (ii) to determine the amount of plant available soil water present throughout a recently harvested site by making PTF predictions (following equations of Saxton and Rawls 1986; 2006) using measured values of soil texture, carbon, coarse fragments and levels of soil density, and (iii) to statistically compare the soil water values with measured soil properties to determine which parameters have the largest influence on available soil water. MATERIALS AND METHODS Site Description The 16-ha study area is located within British Columbia‟s southern interior plateau (50°50'0.37"N 120°25'6.35"W) (Figure 2.1). The site is located within the IDFxh2 biogeoclimatic subzone (Lloyd et al. 1990). Site elevation ranges from 1010 to 1080 m (mean of 1032 m). The site is uniformly sloping (19%) and faces southeast. Prior to harvest in 2015, the dominant vegetation on the site included climax 15 stands of interior Douglas-fir (Pseudotsuga menziesii) with an understory rich in herbs and dominated by pinegrass (Calamagrostis rubescens). Average mean annual temperature (MAT) over recent decades (1981-2010) was 5.1oC. Mean annual precipitation (MAP) during this period equaled 435 mm, receiving less than one half during the growing season of May through September (209 mm). Hargreaves reference evapouration (Eref) was approximately 626 mm, while the Hargreaves climatic moisture deficit (CMD) was 335 mm. The thirty-year average annual relative humidity (RH) was 63% (Wang et al. 2016). 16 Figure 2.1: Aerial map of the harvested Isobel lake study site showing soil pit and soil water sensor locations. 17 Research Design Twenty soil pits were excavated using a backhoe creating three 100 cm x 100 cm measurable profile faces within each pit. The pits were arranged approximately along the same elevation through the mid-slope position of a hillside that was harvested in 2015. The soil pits were separated by a minimum distance of 15 meters between pits 10 and 16, and a maximum distance of 607 meters between pits 4 and 18. The pit locations were chosen based on a larger study that involves 5 individual trials of tree growth under treatments of grazing and no-grazing by range cattle. The pits were randomly placed to surround the 5 individual trials with the intention of characterizing the soil conditions throughout the extent of the harvested area. Soil subsamples were collected from three faces of each soil pit (A, B and C) at 0-15 cm, 15-30 cm, 30-60 cm and 60-100 cm depths (one sample per profile face at each depth). Samples were initially placed into plastic bags and air dried for at least two weeks before being processed. Air dried samples were crushed with a wooden rolling pin and rubber mallet to break up the soil aggregates before being sieved through a 2 mm screen to remove all coarse fragments larger than 2 mm. The samples were transferred into plastic containers to be stored for future analysis. Soil sampling began in the fall of 2015, but an early snowfall delayed the remaining sampling until August and September of 2016. The location of the soil pits was recorded using a handheld Garmin GPS system. Climate and Moisture Data A total of 12 soil moisture stations (EC5 soil moisture sensors, Decagon Devices, Pullman, WA, USA) were established across the 16-ha site in July 2015, and have been continuously monitoring soil volumetric water (θ) at 10, 25 and 50 cm below the mineral soil surface for over two years, capturing the majority of three growing seasons. One water measurement was recorded every hour and all data 18 were monitored and stored on a data logger (EM50, Decagon Devices). Two weather stations were also established within the harvested landscape to record additional climatic variables of precipitation (ECRN-100 tipping buckets, Decagon Devices), air temperature and humidity (VP-4, Decagon Devices) under disturbed sections of clearcut forest. Soil volumetric water content recorded during the 2016 and 2017 growing season was used to demonstrate the rate at which soils of this particular site deplete their water stores as a result of transpiration by competing vegetation as well as evapouration. This information was used to assess the number of days newly established seedlings will have access to available soil water before resources are depleted from the site based on the total soil water holding capacity (SWHC). Measured Soil Physical and Chemical Properties Volume of coarse fragments (%) for each soil depth was recorded in the field using visual observation methods (LMH 25, page 82) (BC Ministry of Forests and Range and BC Ministry of Environment 2010). Soil particle size distribution (% sand, silt and clay) was measured for all 240 samples (20 soil pits x 3 sub-samples x 4 depths) during February of 2017 following the hydrometer method described by Kalra and Maynard (1991). Total soil carbon and nitrogen concentrations were determined by dry combustion (Nelson and Sommers 1996), and soil pH (1:2; CaCl) was measured using standard methods (Thomas 1996). Soil organic matter (SOM) was calculated by multiplying the organic carbon concentration by a factor of 1.72 (Howard and Howard 1990). One profile face from each pit was also measured for soil fine fraction bulk density (BD, g cm-3) at all four soil depths (total of 80 BD samples) during the summer of 2017 following the excavation method (Maynard and Curran 2007). The samples obtained for bulk density analysis were processed by first air drying them in paper bags for two weeks, and sieving the samples through a 2 mm screen to 19 remove all coarse fragments and wood debris. The mass and volume of coarse fragments and plant materials within each sample was determined and subtracted from the total mass and volume of soil fine fraction to be used for the final calculations of BD. The relative weight of rock fragments larger than 2 mm in diameter was used to provide a lab estimate for the percentage of coarse fragments within each collected sample. The lab tested values of coarse fragment concentrations (% volume) were used for further statistical analysis of the soil properties and prediction of the hydraulic parameters. Derived Soil Variables The collected soil data from each face within each pit were analyzed to determine the soil classification at the order level for each profile according to the Canadian System of Soil Classification (Soil Classification Working Group 1998). Total soil porosity was calculated from the bulk density measurements as: ( * +) [1] Where BD is the fine fraction bulk density of the soil sample, and 2.65 represents the density of the individual soil particles (density of quartz rock = 2.65 g cm-3). Maximum bulk density (MBD) was determined for each sample by applying the model outlined by Zhao and others (2008) for moderate and low plastic samples, and was calculated as: [2] Where clay and silt equal the proportion of each within a soil sample (%). Relative bulk density (RBD) was calculated as the ratio between fine fraction BD measured in 20 the field and the predicted MBD (Zhao et al. 2010). The final RBD values were used to assess post-harvest compaction on the site and were related to a threshold of 0.80 outlined by Zhao and others (2010) as the approximate limit where Douglas-fir tree growth begins to decline. Available water capacity (AWC) of each soil was determined by subtracting PWP from the FC value for each soil sample. Field measurements of texture and SOM concentration (g kg-1) were used as input variables in the soil water characteristics equations for predicting FC and PWP (Saxton and Rawls 2006). The ratio of the measured fine fraction bulk density values and the predicted bulk density values from the Saxton and Rawls (2006) models were used as the compaction factor in the model and ranged from loose (0.9) to very compacted (1.3) within the model limitations (predicted density equals 1.0 = normal compaction). Available water capacity was calculated as: [3] Where θFC is the volumetric water content at field capacity (33 kPa, % vol.) and θPWP is the volumetric water content at permanent wilting point (1500 kPa, % vol.). Using the soil water characteristic equations from Saxton and Rawls (2006), θFC was calculated as: [4] and θPWP was calculated as: [5] 21 Where is 33 kPa water (first solution) of each sample, is 1500 kPa water (first solution) of each sample, S is sand content (% w), C is clay content (% w) and OM is organic matter content (% w). Soil saturation was also calculated using the soil water characteristics equations from Saxton and Rawls (2006) as: [6] Where θFC is calculated as above, θ(S-33) is the volumetric water content at saturation minus 33 kPa water (normal density) and S equals sand concentration (% w.). Soil water holding capacity (SWHC) for each soil layer was determined to account for coarse fragment content when considering available water capacity of the soils. The SWHC (mm) for each sample was calculated following the formula used by Kirchen et al. (2017) as: ∑ ( ( )) [7] Where RVi represents the volume of coarse fragments in layer i (%), θFCi and θPWPi represent the water values at field capacity and wilting point (m3 m-3), respectively, of soil layer i, and Hi represents the height of the soil layer i (mm). Statistical Analysis Measured Soil Properties The measured soil data were analyzed using descriptive statistics, correlation and stepwise multiple regression analyses to explore the interrelationships between the measured variables. The Fligner-Killeen test (Conover et al. 1981) determined if 22 the assumption of equal variances was true between the variables and soil pit number. Analysis of variance (ANOVA) test was utilized to determine if differences in measured soil properties exist between the 20 individual soil pits, as well as to test for statistical differences between the 4 measured depths. Coefficient of variation statistics were used to evaluate the variability of individual soil properties measured between different soil pits as well as measurements taken from individual faces within a single soil pit. Derived Soil Properties ANOVA followed by a Tukey-HSD post-hoc analysis was completed to determine if/where significant differences in AWC exist between the soil pits. The derived values of Saturation, FC, PWP and SWHC were analyzed using correlation analysis to determine which measured soil variables had the largest effect on the soil water availability of the site. Descriptive statistics were utilized to explain how differing soil texture properties, relative bulk density values, coarse fragment and carbon concentrations influenced the available water content at the site, as well as to give forest managers insight into how values of these crucial variables relate to soil water patterns (particularly in the top 30 cm of the soil profile). The averaged soil pit dataset was divided into either Luvisol or Brunisol soil order to test if they have significantly different soil properties. ANOVA statistics were utilized to determine if significant differences exist between the Luvisolic (pits 2, 5, 7 and 17) and Brunisolic soil pit data. Multiple logistic regression was used to model the accuracy of the soil data for predicting whether that soil will classify as either a Brunisolic or Luvisolic soil order at different sampling scales and at different depths (Soil Classification Working Group 1998). Residual plots and Cook‟s Distance statistics were utilized to assess normality of the data and to determine if large outliers existed within the dataset. 23 For linear and logistic regression models that included multiple explanatory variables, collinearity among input variables was assessed using the variance inflation factor (VIF) and was calculated as: [8] Where R2 is the coefficient of determination of multiple linear/logistic regression of each independent variable against all other independent variables (O‟Brien 2007). The predicted probabilities from the logistic regression models were assessed using cross validation methods, and by assessing the number of correctly classified soil orders obtained from the 2 x 2 classification table (Peng et al. 2002), which is an account of the number of true-positive and true-negative event outcomes (eg. soil correctly or incorrectly classified as Luvisol) that were predicted by the logistic regression model (Krueger et al. 2017). All statistical tests were performed in R_3.3.3 (R Core Team 2017) using the JGR gui package for much of the analysis (Helbig et al. 2017). A power analysis was performed to determine the estimated sample size required for capturing a significant linear relationship between AWC values and coarse fragment concentration. A linear regression power analysis was completed in „R‟ using the „pwr.f2.test‟ within the „pwr‟ package. The estimated effect size for the model was 0.637, and was calculated as: [9] Where f2 is the predicted effect size and R2 is the coefficient of determination from the linear model for predicting AWC using coarse fragment concentration values measured at the Isobel site as the explanatory variable. This analysis was completed for the data measured within the top 30 cm of soil (n=40), and the significance level for the power analysis was set at p = 0.05 and at p = 0.10. 24 RESULTS Climate and Moisture Data During the summer of 2017 (May-September) the Kamloops area (345 m a.s.l.) received a total of 45.7 mm of precipitation, with 31.1 mm falling in May, and only 8.6 mm between the months of June and August. This marked the driest year ever recorded in the area dating back to 1895 (Foulds 2017). At the Isobel Lake site (1030 m a.s.l.), a total of 103 mm of precipitation was recorded by the ECRN-100 tipping buckets throughout the entire growing season, with 75 mm falling in May, and only 29 mm falling between June and the end of August. From July 2015 to September 2017, a total of 957.6 mms of rain fell on the Isobel Lake site. This resulted in an average of 275 mm during the summer of 2016 (May through September), and of 105 mm during the summer of 2017. This is indicative of overall lower precipitation levels recently throughout the area compared with historic averages. During the same period (July 2015 and September 2017), average RH at the site was 0.704, with the winter months (October through April) averaging around 0.825, and the summer months (May-September) averaging at 0.573 (average maximum growing season [May-September] vapour pressure deficit (VPD) for 2016 = 1.488 kPa and for 2017 = 2.095 kPa). MAT between July 2015 and September 2017 was recorded to be 7.09oC (average yearly VPD = 0.298 kPa), indicating higher average temperatures than the area has experienced throughout past decades. Within the top 30 cm, soils at the Isobel Lake site lost roughly 1.2 mm of water per day during the 2017 growing season (2 years after harvest) mainly as a result of evapotranspiration by competing vegetation. Due to very low inputs from precipitation during the 2017 growing season, total soil water fell sharply within all measured depths between mid-May and late-June demonstrating the severe effects of drought in combination with competing vegetation on these sites (Figure 2.2). 25 θ (mm) Soil Volumetric Water (mm) 120 100 80 60 40 20 0 Figure 2.2: Total volumetric water content (θ) (in mm) recorded daily (within the top 30 cm of soil) from the Isobel Lake study site from April of 2016 until September of 2017. Measured Soil Properties Soil physical and chemical properties are summarized in Table 2.1.Surface soils had similar amounts of sand and silt, with lesser amounts of clay. Particle size distribution of subsurface soils was usually dominated by sand with considerable amounts of silt, and lesser amounts of clay within all four measured depths. Sand concentration ranges from an average of 40.3% (95% confidence interval [CI] = 38.9, 41.8) at the surface to 48.4% (95% CI = 46.1, 50.8) at a depth of >60cm, while silt had an opposite relationship with depth, decreasing from an average of 38.6% (95% CI = 37.1, 40.1) at the surface to 33.0% (95% CI = 31.7, 34.2) in the lowest measured layer. Considerable amounts of clay are present at all four depths, ranging from a surface average of 21.0% (95% CI = 20.1, 22.0), to 18.6% (95% CI = 16.7, 20.5) at depths greater than 60 cm. Of the 240 collected soil samples, 89.6% 26 are categorized as a loam (L), 7.9% are sandy loam (SL), 1.7% are clay loam (CL), and 0.8% (2 samples) are categorized as sandy clay loam (SCL) according to the Canadian System of Soil Classification (Soil Classification Working Group 1998). Coarse fragment percentage appeared to increase with depth, varying from 14.1% (95% CI = 9.0, 19.3) at the surface to 36.5% (95% CI = 30.4, 40.7) at the lowest measured layer. However, the second depth (15-30 cm) had the highest concentration of coarse fragments throughout the site (40.7%, 95% CI = 31.5, 49.9). Surface soils are richer in both total carbon and nitrogen contents (means of 31.45 g kg-1 [95% CI = 26.94, 35.95] and 1.85 g kg-1 [95% CI = 1.51, 2.19], respectively) than the deeper soil samples (means of 7.07 g kg-1 [95% CI = 4.45, 9.69] for carbon and 0.41 g kg-1 [95% CI = 0.27, 0.55] for nitrogen at >60 cm). However, the carbon to nitrogen ratio (C:N ratio) remains fairly constant (between 17 and 18%) throughout the entire soil profile. Additionally, soil pH slightly increases with depth from an average of 5.46 (95% CI = 5.28, 5.64) at the surface to 6.57 (95% CI = 6.31, 6.82) at depths greater than 60 cm. Two measurements of BD were determined as outliers and were removed from the analyses. These samples exceeded a BD of 2.0 which is exceptionally higher than expected as the density of rock is 2.65 g cm -3, and contained greater than 70 % coarse fragments leaving only a small proportion of soil to measure for density. BD generally increased with depth (mean= 1.03 g cm-3 at surface [95% CI = 0.95, 1.12] up to mean= 1.41 g cm-3 at greater than 60 cm depths [95% CI = 1.34, 1.49]). The resulting total soil porosity decreased with depth opposite to bulk density (mean of 61.1% at surface to 46.6% at depths greater than 60 cm). 27 Table 2.1: Averaged soil properties measured throughout the entire study site (n=60 for variables of clay, silt, sand, pH, total C, total N, and C:N for each depth range: 20 pits X 3 faces/pit at all 4 depths. n=20 for variables of CF, BD and porosity: 1 measurement per depth within all 20 pits). Values within brackets represent one standard deviation from the mean. Total C = total carbon, Total N = total nitrogen, C:N = carbon to nitrogen ratio, CF = coarse fragments (>2 mm), and BD = fine fraction bulk density. Superscripted letters represent significant differences between values within different depth ranges (assessed using the Tukey multiple comparison of means test). Depth Range (cm below soil surface) 15-30 30-60 ab 23 (2.0) 21 (2.0)ab 34 (3.0)b 32 (3.0)c 43 (3.0) a 47 (4.0) b a 5.7 (0.45) 6.2 (0.05)b Clay (%) Silt (%) Sand (%) pH 0-15 21 (2.0)a 39 (3.0)a 40 (3.0) a 5.5 (0.04)a 60-100 18 (4.0)ac 33 (3.0)b 48 (6.0) b 6.6 (0.05)c Total C (g kg-1) 31 (12.0)a 16 (6.0)b 8.9 (4.4)c 7.1 (4.3)c Total N (g kg-1) C:N CF (%) BD (g cm-3) Porosity (%) 1.9 (0.84)a 17.8 (3.8)a 14 (11.0)a 1.03 (0.18)a 61 (10.0)a 0.92 (0.36)b 18.0 (4.7)a 41 (19.0)b 1.27 (0.28)b 52 (10.0)b 0.51 (0.20)c 17.5 (4.7)a 32 (11.0)b 1.34 (0.18)bc 49 (7.0)bc 0.41 (0.22)c 16.9 (5.4)a 37 (13.0)b 1.41 (0.15)c 47 (6.0)c Measured bulk density values were verified by comparing the predicted BD values (derived from the Saxton and Rawls (2006) equations) with the measurements obtained in the field (Figure 2.3). The relationship between the measured and predicted was significant (p-value < 0.001; R2 = 0.319). Predicted BD (g cm-3) 28 1.8 1.6 1.4 1.2 1.0 y = 0.3375x + 0.9901 R² = 0.319 0.8 0.5 1.0 1.5 2.0 Measured BD (g cm-3) Figure 2.3: Scatterplot of the linear regression model comparing the measured BD values (obtained from field measurements) and predicted BD values (from Saxton and Rawls (2006) equations). BD = fine fraction bulk density. Derived Soil Properties When examining the averaged values from all three subsamples within a pit, the soils on the site are predominantly Orthic Eutric Brunisols with an average pH greater than 5.5. However, 20% (4/20) of the sampled pits (including pits 2, 5, 7 and 17) showed characteristics of Orthic Gray Luvisols, with adequate concentrations of clay eluviated from the „A‟ horizon into the „B‟ to classify it as a true Bt soil horizon (Soil Classification Working Group 1998) (Figure 2.4). When the data from each individual profile face is considered, only roughly 12% (7/60) of the faces (specifically, 2A, 2B, 5B, 7A, 8C, 11A and 17C) classify as being Orthic Gray Luvisols, while the remaining 53 soil faces have not adequately developed past the Orthic Eutric Brunisol classification. 29 2.4A 2.4B LFH LFH Ah Ah Bm Ae Bt C Ck Figure 2.4: Soil profile of pit 7, face 'A' (2.4A): representative of a Orthic Gray Luvisol also found in pits 2, 5 and 17 at the Isobel site. Soil Profile of pit 12, face „A‟ (2.4B): representative of an Orthic Eutric Brunisol discovered in the remaining 16 soil pits. The average RBD of the site was 0.702 (n=78, 95% CI = 0.671, 0.732), indicating that overall soil density is below the threshold limit known to restrict Douglas-fir seedling growth (RBD >0.72) at this site (Zhao et al. 2010). Although 18 out of the 78 samples expressed RBD values greater than 0.80, the top 15 cm of soil did not experience RBD values above 0.767 and averaged at 0.574 (95% CI = 0.527, 0.621). The greatest values of RBD were found in the bottom two soil layers from 30-100 cm (mean value of 30-60 cm = 0.745, mean value of 60-100 cm = 0.786). The two outlier BD values that exceeded 2.0 were also eliminated from the soil hydraulic calculations, and therefore the average values of PWP, FC, AWC and SWHC were calculated with n=78; n=20 in depths 1 and 3, n=19 in depths 2 and 4. The concentration of available water per cm of soil declined with depth. On average, 30 the top 15 cm of soil exhibits greater soil water holding capabilities than the underlying 15-30 cm layer, which can hold onto water better than the 30-60 cm layer below it. Available water per cm of soil volume decreased with depth on average. The SWHC was highest in the 0-15 cm depth range and lowest in the 60-100 cm layer when considering mm of water per cm of soil. Additionally, proportions of FC, PWP and AWC all declined with increasing depth in relation to soil volume (Table 2.2) (Figure 2.5). Table 2.2: Average values (in mm of water) for the calculated soil hydraulic variables of saturation, field capacity (FC), permanent wilting point (PWP), available water capacity (AWC) and soil water holding capacity (SWHC) within the top 3 measured depths throughout the entire study site (n=59). Numbers in brackets represent one standard deviation from the mean. Volume of water (mm) in soil layer: Saturation 0-15 91.6 (10.2) Depth (cm) 15-30 78.0 (15.7) 30-60 148.3 (20.1) Total (mm) 60cm 30cm 317.9 169.6 Field Cap. (FC) 49.6 (3.3) 45.4 (3.7) 83.8 (8.1) 178.8 95.0 Wilting Point (PWP) 24.0 (1.8) 23.3 (1.6) 42.5 (5.4) 89.8 47.3 AWC - Avail. Water Cap. (FC-PWP) 25.6 (2.4) 22.1 (3.3) 41.3 (4.4) 89.0 47.7 SWHC 22.0 (3.9) 13.7 (4.9) 28.3 (6.2) 64.0 35.7 18 40 16 35 14 30 12 FC (%) PWP (%) 31 10 8 6 25 20 15 4 10 2 5 0 0 0-15 15-30 30-60 60-100 0-15 20 18 16 14 12 10 8 6 4 2 0 Depth Range (cm) 70 60 Saturation (%) AWC (%) Depth Range (cm) 15-30 30-60 60-100 50 40 30 20 10 0 0-15 15-30 30-60 60-100 Depth Range (cm) 0-15 15-30 30-60 60-100 Depth Range (cm) Figure 2.5: Averaged proportion (%) of each soil hydraulic property within each sampled depth range. PWP = permanent wilting point, FC = field capacity and AWC = available water capacity. Error bars represent standard error of the mean. Letters represent significant statistical differences (p = 0.05) using Tukey multiple comparison of means test. Adding this new soil hydraulic information with the earlier figure explaining total volumetric water content at the site, a clearer representation of soil water holding capacity in relation to site field capacity and θ is revealed (Figure 2.6). 32 Figure 2.6: Average soil volumetric water content (mm) within the upper 30 cm soil layer in relation to the calculated average values of site FC, PWP and SWHC in mm. SWHC (in mm) represents the fraction of soil water available for plant absorption above the permanent wilting point. FC = field capacity, PWP = permanent wilting point, AWC = available water capacity, and SWHC = soil water holding capacity. Statistical Analysis Measured Soil Properties Soil pit number was not a significant source of variation for any of the measured soil properties. Soil layer depth showed a strong correlation (> 0.30 or < 0.30) with all the measured soil properties (Table 2.3). The strong correlation with texture is expected as the average soil pit experienced obvious patterns in sand, silt and clay with depth. Depth showed strong negative correlation with silt (-0.61), carbon concentration (-0.77) and nitrogen concentration (-0.75), and expressed positive correlations with pH (0.66) and BD (0.48). The variables of clay and CF concentration also had notable correlation with depth (-0.30 and 0.38, respectively). In addition, total C had very strong positive correlation with total N (0.92), and both silt and clay had strong negative correlations with sand content (-0.76 and -0.67, 33 respectively). Sand was negatively correlated with C and N (-0.62 and -0.56, respectively), while silt expressed positive relationships with the two (0.68 and 0.63). Furthermore, sand was positively correlated with pH (0.53), carbon concentration was negatively correlated with pH (-0.57), and fine fraction bulk density was positively correlated with CF concentration (0.59) (Table 2.3). All the measured soil variables appear to have normal distributions. Table 2.3: Soil physical and chemical property correlation matrix using data collected from the Isobel lake study site. The included variables are as follows: BD= fine fraction bulk density; CF= lab determined coarse fragment content; pH= soil pH; C:N= carbon to nitrogen ratio; N = total nitrogen content; C = total carbon content; Sand, Silt and Clay= textural proportions of soil sample; Depth= soil layer depth, and Soil Type = Brunisol or Luvisol. Bold-texted and highlighted cells represent correlations at the 0.05 significance level, and bold-texted cells represent correlations at the 0.10 significance level. Depth Clay Silt Sand Soil Type C (g kg-1) N (g kg-1) C:N pH CF Depth Clay Silt Sand Type C N C:N pH CF BD 1.00 -0.30 1.00 -0.61 0.02 1.00 0.65 -0.67 -0.76 1.00 0.00 -0.01 0.23 -0.16 1.00 -0.77 0.16 0.68 -0.62 0.08 1.00 -0.75 0.14 0.63 -0.56 0.13 0.92 1.00 -0.03 0.11 0.06 -0.11 -0.21 0.21 -0.12 1.00 0.66 -0.24 -0.50 0.53 -0.15 -0.57 -0.49 -0.25 1.00 0.38 -0.07 -0.37 0.32 0.10 -0.38 -0.43 0.17 0.08 1.00 BD (g cm-3) 0.48 -0.10 -0.37 0.34 0.01 -0.48 -0.46 -0.02 0.22 0.59 1.00 Silt was statistically related with total C concentration (p < 0.001, R2 = 0.46) and total N concentration (p < 0.001 and R2 = 0.39) based on linear regression models, however, when running the model using both C and N as explanatory variables the coefficient of determination does not increase (R2 = 0.45) as C and N co-vary. Carbon concentration explained 32% of the variation in pH values, while sand explained 27% when individually modelled. Furthermore, CF concentration alone explained 26% of the variation in BD measurements. However, stepwise multiple regressions were used to determine the best subset of soil properties for explaining the variation in BD. The resulting best model included C in addition to CF 34 concentrations. The combination of these two variables for predicting BD explained roughly 36% of the variation when tested with linear regression analysis (p < 0.001). Finally, depth explained significant amounts of variation in all of the measured soil properties within the study (all p-values < 0.001), with C:N being the only exception as it did not exhibit much variation throughout the entire soil profile of most pits. The coefficient of variation analysis showed that all soil variables had greater variation between the different soil pits, and the individual faces within a pit had less variability than the individual pits themselves (Table 2.4). The soil variables of CF and BD only had variation measured between pits as only one measurement was taken from within each pit at each depth, and therefore could not have additional variability found within a single pit. 35 Table 2.4: Coefficient of variation (CV: relative standard error) percentage of each measured soil property explained by the different levels of sampling. The average between pits shows the average variation found between each soil face of all 20 soil pits, and within pits explains the average variation found within a single pit. CV = standard deviation divided by the mean. CF = coarse fragment concentration, BD = fine fraction bulk density, C:N = carbon to nitrogen ratio. Depth (cm) 0-15 Variable Sand Clay Silt CF BD Carbon Nitrogen C:N pH 15-30 Sand Clay Silt CF BD Carbon Nitrogen C:N pH 30-60 Sand Clay Silt CF BD Carbon Nitrogen C:N pH 60-100 Sand Clay Silt CF BD Carbon Nitrogen C:N pH Coefficient of Variation (%) Avg. Between Avg. Pits Within Pit 8.59 5.55 8.89 6.02 8.48 6.57 78.0 N/A 17.5 N/A 39.8 27.7 48.3 23.1 20.5 10.4 7.30 4.97 6.81 5.18 10.3 7.20 7.64 6.27 48.6 N/A 21.8 N/A 35.2 29.8 38.6 23.9 25.1 12.8 7.87 4.51 8.21 5.83 11.1 8.95 8.52 6.61 33.8 N/A 13.2 N/A 41.4 29.8 37.0 21.7 25.1 13.9 7.90 3.81 11.7 6.14 18.9 10.2 10.6 7.38 35.9 N/A 10.8 N/A 69.1 30.0 55.4 20.4 31.8 15.6 7.53 4.18 36 Derived Soil Properties The Luvisolic and Brunisolic soil orders found at the site did not contain many significant differences with one another when testing both the measured and derived soil variables. The only variables that showed a significant difference were silt concentration (p = 0.041) and C:N ratio (p = 0.058). None of the soil hydraulic properties showed significant variation between the two soil types. Using both forward and backward stepwise logistic regression for predicting soil type (Luvisol = 1, Brunisol = 0) from the average measured soil properties (n=80), the model that produces the lowest value of Akaike Information Criterion (AIC) and the highest predictive accuracy includes the variables of silt, C:N ratio and CF concentration (p = 0.005, p = 0.020 and p = 0.040, respectively). When analyzing the VIF of each variable, all three had VIF values close to one indicating low collinearity among the variables included in the logistic model. When assessing the cross-validation estimate of accuracy, the model is 82.5% accurate at predicting the soil type. However, upon further analysis of the 2 x 2 classification table, the resulting model predicted 62 out of 64 observations correctly as Brunisols, but only classified 4 out of 16 observations correctly as Luvisols. In other words, the model classifies most observations as Brunisolic and does not have a very strong capacity to distinguish the different soil types. A similar test was done using only the data from the top 15 cm of soil (using all 3 faces of each pit, n=60) to see if surface soil properties can predict soil type at the site as well. The logistic regression model that had the highest predictive accuracy and no collinearity among input variables included the soil properties of silt (p = 0.130) and carbon concentrations (p = 0.037). The cross-validation estimate of accuracy for this model was 86.7%, which was no lower than the model which included all of the measured soil properties. Again, however, when observing the 2 x 2 classification table, the model predicted 52 out of 53 observations correctly as Brunisols, but did not correctly classify any of the 7 observations as Luvisols, once again indicating a low capability for the model to distinguish between the two. The same test was completed looking at the top 30 cm 37 of soil (n=120), and the best model discovered through stepwise logistic regression again included the variables of silt (p = 0.016), C:N ratio (p = 0.050) and CF concentration (p = 0.107). None of the three variables showed substantial collinearity and the cross-validation estimate of accuracy equaled 88.3%. The 2 x 2 classification table showed similar results as before, with a high success rate at predicting Brunisols (all 106 classified correctly) and a very low success rate at predicting Luvisols (0 out of 14 correctly classified). Despite coarse fragments not having any influence on the predicted values of FC, PWP, Saturation or AWC, it was highly correlated with all of the derived soil hydraulic properties (Table 2.5). CF concentration showed the highest correlation with AWC (-0.61) compared to soil texture (0.56 for silt, and -0.55 for sand), C concentration (0.60) and pH (-0.34). The correlation matrix demonstrates how an increase in CF concentration would decrease the amount of water a soil can contain at saturation, as well as decrease the AWC and overall SWHC of a soil. The high correlation between CF and SWHC (-0.93) was expected as it was used as an input variable in the SWHC equation. Additionally, an increase in pH (which is negatively correlated with carbon) would decrease the AWC at the site as well as SWHC. Furthermore, an increase in organic carbon concentration would ultimately increase the SWHC of the site as it directly relates to soil organic matter. CF concentration explains 36.5% of the variation in AWC values when tested with linear regression analysis, despite not being an input variable in the Saxton and Rawls (2006) equation (p < 0.001). When testing the linear model using CF as the explanatory variable for predicting AWC within the top 30 cm of soil only, 38.9% of the variation in AWC is explained (p < 0.001). 38 Table 2.5: Correlation matrix demonstrating the relationships between the measured soil variables and the predicted soil hydraulic variables. Soil Type = Brunisol or Luvisol, C = Total Carbon, N = Total Nitrogen, CF = Coarse Fragments, PWP = Permanent Wilting Point, FC = Field Capacity, AWC = Available Water Capacity, and SWHC = Soil Water Holding Capacity. Bold-texted and highlighted cells represent correlations at the 0.05 significance level, and bold-texted cells represent correlations at the 0.10 significance level. Clay Silt Sand Type C N C:N pH CF Por Clay 1.00 Silt 0.02 1.00 Sand -0.67 -0.76 1.00 Soil Type -0.01 0.23 -0.16 1.00 Carbon 0.16 0.68 -0.62 0.08 1.00 Nitrogen 0.14 0.63 -0.56 0.13 0.92 1.00 C:N 0.11 0.06 -0.12 -0.21 0.21 -0.12 1.00 pH -0.24 -0.50 0.53 -0.15 -0.57 -0.49 -0.25 1.00 CF -0.07 -0.37 0.32 0.10 -0.38 -0.44 0.17 0.08 1.00 Porosity 0.11 0.37 -0.34 -0.01 0.48 0.46 0.02 -0.22 -0.59 1.00 Saturation 0.11 0.37 -0.34 -0.01 0.48 0.46 0.02 -0.22 -0.59 1.00 PWP 0.85 0.38 -0.84 0.04 0.66 0.59 0.20 -0.48 -0.26 0.34 FC 0.59 0.55 -0.80 0.05 0.73 0.67 0.14 -0.47 -0.52 0.78 AWC 0.20 0.56 -0.55 0.05 0.60 0.56 0.04 -0.34 -0.61 0.97 SWHC 0.13 0.56 -0.50 -0.02 0.58 0.60 -0.11 -0.24 -0.93 0.79 Sat PWP FC AWC SWHC 1.00 0.34 0.78 0.97 0.79 1.00 0.84 1.00 0.47 0.87 1.00 0.41 0.74 0.83 1.00 Results from the power analysis for determining the sample size needed to capture the apparent relationship between AWC and CF concentration within the top 30 cm of soil was between 16 (p = 0.10) and 19 (p = 0.05) soil pits (using only CF concentration as the predictor variable) depending on the desired level of confidence (Figure 2.7). Sample Size (n) 39 20 18 16 14 12 10 8 6 4 2 0 p = 0.10 p = 0.05 0.5 0.6 0.7 0.8 0.9 Power Figure 2.7: Estimated recommended sample size (number of soil pits to dig) to identify a significant relationship between AWC and CF at the site (within the top 30 cm of soil), using only coarse fragments (% vol.) as the explanatory variable in a linear regression model. The „pwr.f2.test‟ was utilized within the „pwr‟ package in „R‟ for 2 2 2 determining the sample size with a significant level of 0.05 and 0.10, f2 = R /(1-R ) (R determined by testing the linear model with AWC as the dependent variable, and CF as the independent), and u = number of variables used in the linear regression model (1; coarse fragments). Power was allowed to fluctuate to determine sample size at different levels of certainty. The Saxton and Rawls (2006) model equations for predicting FC and PWP were explored using average loam textured soils from the top 30 cm of soil (sand = 41% and clay = 22%), the average soil organic matter (SOM) content from the top 30 cm of soil (4.2%) and the average density factor (DF) value from the top 30 cm (0.90). Changes in SWHC (mm) were observed with differentiating levels of coarse fragments. Using a range of coarse fragment concentrations from 10 to 90% (in increments of 10%), the total SWHC in mm was determined for the top 300 mm of soil. The average SWHC of the site rapidly decreased from 42.1 mm (10% CF) to 4.7 mm (90% CF) with increasing levels of coarse fragments (Figure 2.8A). When other variables were held constant, changes in SOM, sand content, clay content and DF significantly altered SWHC. When holding the CF concentration constant within the soil water equations (using the average CF value from the top 30 cm [41%]), and using the same average values for texture and DF as above, SOM was allowed to fluctuate (from a value of 0 to 10% concentration) to determine its 40 effect on SWHC in the top 30 cm of soil. The average change in SWHC with changing SOM equaled approximately 0.80 mm/1% change in SOM (Figure 2.8B). Next, sand content was allowed to fluctuate while holding the other variables constant at their measured averages within the model. The average change in SWHC as a result of fluctuating sand content was approximately 1.70 mm/5% change in sand (Figure 2.8C). Clay was the next variable to manipulate, and was found to alter the SWHC by a value of roughly 1.20 mm/5% change in clay (Figure 2.8D). Finally, DF was manipulated to determine its effect on SWHC, and was found to cause a change of 1.76 mm/0.1-unit increase in DF (0.9 - 1.3) (Figure 2.8E). 41 SWHC (mm) 50 40 30 20 10 0 2.8B 50 40 30 20 10 0 Coarse Fragments (% w.) SOM (% vol.) 2.8C 2.8D 50 40 30 20 10 0 SWHC (mm) SWHC (mm) SWHC (mm) 2.8A Sand (dec. % vol.) 50 40 30 20 10 0 Clay (dec. % vol.) SWHC (mm) 2.8E 50 40 30 20 10 0 Density Factor (0.9 - 1.3) Figure 2.8: Effect of fluctuating soil variables used in the soil water characteristics equation (Saxton and Rawls 2006) on SWHC (mm). While one variable is fluctuating, the rest were held constant at their average value from the top 30 cm of soil measured at the site. Figure (2.7A) = Fluctuating coarse fragment concentration (% w.). Figure (2.7B) = Fluctuating soil organic matter (% v.). Figure (2.7C) = Fluctuating sand content (dec. % v.). Figure (2.7D) = Fluctuating clay content (dec. % v.). Figure (2.7E) = Fluctuating density factor. SWHC = soil water holding capacity. 42 DISCUSSION Soil Forming Factors Many of the measured soil properties at the Isobel Lake site followed pattern distributions that are typical for soils developed over glacial till parent materials within dry, forested ecosystems. The site was dominated by moderately coarse (loam) textured soils. The elevated concentrations of silt within the top 15 cm, suggests the site is capped with a layer of aeolian material. Because the site has developed over a till soil, the high concentration of coarse fragments was also expected, although the high variability in coarse fragments that was unexplained by topography on a small scale was a notable discovery. Total carbon ranged between an average of 3% at the surface to less than 1% at 60-100 cm depths, and total nitrogen ranged from 0.17% to 0.04% from the surface to depth. The measured values for these physical and chemical soil properties are comparable to the values reported by Hope (2006) from within very similar dry IDF ecosystems which are also located just north of Kamloops, BC. The overall classification of soil order within the study site depended on the scale of measurement. Soil order within the study site showed a difference in classification when looking at the 60-individual soil faces as opposed to the averaged data from the 20 soil pits. When observing the total dataset with 60 profile faces (including all data from faces A, B and C from all 20 pits), only 12% of the individual profile faces were classified as Luvisolic soils, while 20% of the profiles fit the criteria when observing the average dataset from all 3 faces within a single pit (one averaged value for each soil property within each depth range of a single soil pit). Hence, overall classification of soil order of some of the soils within the study site depends on the scale at which the sampling is completed. The fact that the logistic regression models did not accurately distinguish between the Luvisolic and Brunisolic soil types observed at the site shows a general lack of patterned variability 43 within these soils, and intensifies the need to measure local soil properties within these sites prior to performing site prescription methods. Despite the minor differences found in soil properties that lead to a change in classification from Brunisol to Luvisol, the soil properties within each soil order are more similar than they are different, making it difficult for regression models to accurately distinguish them. The clay accumulation in the sub-layers of Luvisolic soils has potential impacts on plant productivity and hydrology. Root penetration can become reduced at lower depths because of the finer textured soil in certain situations where excess compaction has occurred. In addition, high soil water levels during snow melt have the potential to lead to decreased slope stability and can cause temporary anaerobic conditions for tree seedling roots when the soil is above the field capacity (Valentine et al. 1978). For these reasons, the management implications for work atop forest soils will differ when working with Luvisolic as opposed to Brunisolic soils. As a result of the poor accuracy of measured soil properties for predicting soil type within these ecosystems, the importance of local soil studies when preparing site prescriptions is crucial for the accurate interpretation of the landscape, and the number of soil pits examined will dictate the confidence one may have with their interpolations. Soil Variability Excluding coarse fragment concentration, soils at the Isobel site generally displayed less variability in physical properties (i.e. sand, silt, clay, BD) compared with chemical properties (i.e. carbon, nitrogen, C:N ratio). This finding is consistent with results from similar studies which look at the fine-scale variability of physical and chemical soil properties (Amador et al. 2000; Cambardella et al. 1994). However, the large variability in coarse fragment concentration found throughout the Isobel site is a notable discovery as this soil property has a strong influence on the water holding capacity of soils and can drastically alter the growing conditions of a 44 microsite within water limited ecosystems (Cousin et al. 2003; Baetens et al. 2009). It is no surprise that the concentration of rocks within the soil profile would decrease water storage capabilities, as they take up a large portion of the overall soil volume and their surface area to volume ratio is far less than soil particles within the soil separate of sand, silt and clay. Therefore, the space occupied by coarse fragments is no longer available for water storage, and specific areas that have higher levels of coarse fragments within the soil will experience lower water holding capacities. Therefore, most of the variability in SWHC throughout the Isobel site is attributed to the variability in coarse fragment concentration, and this relationship is likely true in many other ecosystems where water is limited and the soils contain moderate to high levels of coarse fragments. Coarse fragment concentration may be the most important variable at this site because it has the highest level of variability and had the largest influence on SWHC values within the upper 30 cm of soil (Figure 2.8A). The important influence coarse fragments has on water storing capabilities of a soil is certainly well established in the literature (Poesen and Lavee 1994; Hillel 1998; Saxton and Rawls 2006), however, my results provide a quantifiable estimate of exactly how much the SWHC may change depending on the amount of rock you may have within a soil profile. Additionally, due to the high variation in rock concentration throughout the landscape, the need to measure local soil properties at multiple different points throughout a site will be important for the accurate estimation of the variability in SWHC. Because SOM was the only soil property that improved the SWHC at the Isobel site, it is important to consider this property when developing management strategies within these forest types. The retention of organic matter within the topsoil should not only improve the amount of available water at the site, but also reduce the risk of further compaction. Zhao et al. (2008) found that increased levels of soil carbon (which is directly related to SOM) significantly reduced compaction (measured as soil BD) in soil samples collected from multiple different forest types including the IDF. Additionally, because the effects of soil compaction often are more 45 severe when soils are dry (Bulmer and Simpson 2005) retaining organic matter within soils of the IDFxh should have a positive effect on improving tree survival and growth. The results from the power analysis with respect to coarse fragment concentration as the explanatory variable for determining the apparent variability in AWC throughout the study site exemplifies a challenge to forest managers in terms of collecting adequate soil samples. To gain a level of confidence of over 90%, a total of 16-19 soil pits must be dug to accurately capture the variation in AWC in relation to coarse fragment concentration. In other words, if one was to conduct a similar study at a different site, and wanted to sample enough locations to accurately capture a significant relationship (at the p = 0.05, or 0.10 level) between CF concentration and AWC, 16-19 soil pits would need to be sampled (given the effect size to be 0.637; which was calculated from the relationship between CF and AWC found at the Isobel site). However, forest managers generally do not dig this many pits and are often limited to sampling one or two areas within the site and extrapolating the data across an area. My study shows that even when sampling on relatively uniform topography, the amount of variation in coarse fragments and AWC far exceeds what can be captured in one or two pits, and sites that are comparatively heterogeneous in topography may require far more samples to gain an accurate representation. Calculated Soil Properties Relative bulk density (RBD) values measured at the Isobel Lake study site were not indicative of a highly compacted forest. Zhao et al. (2010) reported that RBD values (> 0.72) in the upper 20 cm of soils in the interior of BC were correlated with a significant decline in the fifth growing season (Zhao et al. 2010). Although the average RBD values within the upper 15 cm of soil at the Isobel site averaged at 0.57, average RBD values at lower depths between 30 and 100 cm exceeded the 46 threshold of 0.72. In addition, Zhao et al. (2010) demonstrate that the effect of high RBD values increases throughout the first seven growing seasons and therefore represents an ongoing management concern that persists further than the initial seedling establishment. Caution should be practiced when planning site preparation and planting procedures at this site to ensure further compaction does not occur because measured RBD values in the surface layers of the Isobel site are approaching the 0.72 threshold found by Zhao et al. (2010). This is especially important within IDFxh ecosystems as these forest types are generally water limited rather than energy (heat) limited (Heineman et al. 2003), and, based on the results from this study, further compaction will reduce the amount of soil water available to plants. In addition, the effects of compaction were found to have a more negative effect on tree seedling growth when the soils were under dry conditions as opposed to being moist (Bulmer and Simpson 2005; 2010). Due to the extremely dry conditions recorded at the Isobel site during the 2017 growing season, further compaction of these soils could have very negative effects on the successful regeneration of these forest types as it will further negatively affect soil water availability. The calculated values for FC and PWP matched well with the soil volumetric water content that was measured throughout the previous three growing seasons. The graph illustrating the volumetric water content in relation to FC and PWP show a general agreement as the total volume of water under drought conditions (2017 growing season) leveled off just below the average predicted value of PWP (~36-38 mm/30 cm depth) from the Saxton and Rawls (2006) equations, indicating that the soil is dry to the point where roots can no longer access the remaining soil water. This result indicates that the soil water characteristic equations created by Saxton and Rawls (2006) for predicting hydraulic properties within agricultural soils is relatively accurate within forested settings as well, assuming local physical and chemical properties can be included in the equations. The author recognizes that Douglas-fir roots can grow much deeper than 30 cm over time; however, this is the 47 depth range at which seedlings will need to access moisture from the most during their first few years of establishment. Within the topsoil of dry Douglas-fir ecosystems, the temporal distribution of precipitation inputs will influence regeneration success of newly established seedlings. The average SWHC at the study site was found to be roughly 35.74 mm/30 cm of topsoil. This means that when all of the gravitational water has drained out and the soil is at field capacity, the soils can store on average 35 mm of water within a 30-cm profile. Previous research examining the water use of interior Douglas-fir suggests that daily early summer water use by trees is between approximately 1 and 1.5 mm (Simpson 2000). My results from the volumetric water data collected from 2015 to 2017 show similar results with an average decrease in total water volume of roughly 1.2 mm/day. Because the total budget of available water at field capacity is only 35 mm, and the site is losing roughly 1.2 mm of water per day due to evapotranspiration, plant roots will only have about 29 days or one month to access soil water before it becomes too dry and all remaining water is held tightly within soil micropores. This was not a concern during the 2016 growing season as sufficient precipitation events appeared to counteract the drying effect of the soils throughout the summer, but during 2017 (when precipitation was well below average throughout the later parts of the growing season) the soil water content decreased rapidly and there were not enough precipitation events to keep the volumetric water concentration above the calculated wilting point. Soil water deficits are known to cause stress to vegetative growth by inhibiting photosynthesis and transpiration and can ultimately lead to mortality under extremely dry conditions (Mathys et al. 2014). Furthermore, climate change models have projected a general increase in temperatures, which will lead to greater decreases in available soil water (Littke et al. 2018), especially in the hottest summer months of July and August. These changes in temperature without a large increase in precipitation will increase the evapotranspiration demands of Douglas-fir seedlings and will undoubtedly hinder tree growth. Additionally, the increased drying of the IDFxh into the future has major implications for increased fire and insect 48 disturbances as the forests become more water stressed (Littke et al. 2018). If forest managers are going to continue to harvest in the dry Douglas-fir ecosystems of BC‟s southern interior, regeneration strategies will have to shift to accommodate for the predicted increases in temperature and corresponding declines in available soil water – a task that may prove to be ever more difficult as time persists. CONCLUSION/SUMMARY The results of this study provide a fundamental dataset which can be used for the development of future soil models within similar ecosystems. By intensively sampling a small forested hillside, specific relationships between soil properties and water availability were determined, and the overall variability of these relationships was explored. Because coarse fragment concentration had a highly negative effect on the soils ability to retain water, it is important for forest managers to determine the stoniness of soils prior to prescribing regeneration strategies, and to ponder the potential detriments to regeneration of Douglas-fir seedlings within rocky soil substrates. Furthermore, because organic matter content was the only soil property which improved SWHC, it is important to consider this variable when planning logging and site preparation strategies within the IDFxh. The retention of organic matter within the soil surface will not only improve water availability, but also decrease the potential for compaction. In the face of increasing average global temperatures and drought conditions, it is of the upmost importance to determine the relationships between fundamental soil properties and related hydraulic characteristics to ensure the continued success of the forest industry within BC. Because these ecosystems are becoming increasingly difficult to regenerate following harvest, information about the effects of soil variability on the variability in AWC will be crucial to forest managers tasked with 49 operating within these forest types. Future research should focus on fine-scale soil variability within other dry forest ecosystems and their relationship with AWC to help solve the ongoing issue regarding a lack of small-scale soil information, and to improve existing soil maps with greater soil water information within the province. LITERATURE CITED Amador, J.A., Wang, Y., Savin, M.C., Görres, J.H., 2000. Fine-scale spatial variability of physical and biological soil properties in Kingston, Rhode Island. Geoderma 98: 83-94. Bear, J. 1972. Dynamics of fluids in porous media. Dover Publications. ISBN 0-48665675-6 British Columbia Ministry of Forests and Range and BC Ministry of Environment. 2010. Land Management Handbook –25: Field manual for describing terrestrial ecosystems – 2nd edition. Research Branch [accessed December 8, 2016]. Bulmer, C. E., and Simpson, D. G. 2005. Soil compaction and water content as factors affecting the growth of lodgepole pine seedlings on sandy clay loam soil. Can. J. Soil Sci. 85: 667-679. Bulmer, C.E., and Simpson, D.G. 2010. Soil compaction reduced the growth of lodgepole pine and Douglas-fir seedlings in raised beds after two growing seasons. Soil Sci. Soc. Am. J. 74(6): 2162-2174. Cambardella, C.A., Moorman, T.B., Nokav, J.M., Parkin, T.B., Karlen, D.L., Turco, R.F., and Konopka, A.E. 1994. Field scale variability in soil properties in Central Iowa soils. Soil Sci. Soc. Am. J. 58: 1501-1511. Conover, W.J., Johnson, M.E., and Johnson, M.M. 1981. A comparative study of tests for homogeneity of variances, with applications to the outer continental shelf bidding data. Technometrics. 23: 351–361. doi: 10.2307/1268225. Fleming, R.L., Black, T.A., Adams, R.S., and Stathers, R.J. 1998. Silvicultural treatments, microclimatic conditions and seedling response in Southern Interior clearcuts. Can. J. Soil Sci. 78(1): 115-126. https://doi.org/10.4141/S97-042 50 Foulds, C. 2017. 2017 was driest summer ever recorded in Kamloops. Kamloops This Week (News). Aberdeen Publishing group. 1365B Dalhousie Dr. Kamloops, BC. V2C5P6. Available from: https://www.kamloopsthisweek.com/2017-driest-summer-ever-recordedkamloops/ Gupta, S.C., and Larson, W.E. 1979. Estimating soil water retention characteristics from particle size distribution, organic matter percent, and bulk density. Water Resources Research. 15(6): 1633-1635. Hagvhardi, A., Leib, B.G., Washington-Allen, R.A., and Ayers, P.D. 2015. Highresolution prediction of soil available water content within the crop root zone. Journal of Hydrology. 530: 167-179. Hebda, R. 2007. Biodiversity: geological history in British Columbia. Biodiversity BC Technical Subcommittee. Available from: http://www.biodiversitybc.org/assets/Default/BBC%20Biodiversity%20and%20 Geological%20History.pdf Heineman, J.L., Hope, G.D., Simard, S.W., Vyse, A., Lloyd, D.L., and Miège, D.J. 2003. The effects of site preparation and harvesting practices on planted seedling productivity and microenvironment in southern interior dry, grassy IDF forests. Res. Br., B.C. Min. For., Victoria, B.C. Tech. Rep. 009. Helbig, M., Urbanek, S., and Fellows, I. 2017. JGR: Java GUI for R. R package version 1.8-4. https://CRAN.R-project.org/package=JGR Hillel, D. 1998. Environmental soil physics. Academic Press, San Diego, California, USA. Hope, G. 2006. Establishment of long-term soil productivity studies on acidic soils in the interior Douglas-fir zone. LTSPS Research Note #08; British Columbia Forest Service – Southern Interior Forest Region. Howard, P.J.A., and Howard, D.M. 1990. Use of organic carbon and loss-on-ignition to estimate soil organic matter in different soil types and horizons. Biology and Fertility of Soils. 9:306-310. iMapBC. 2018. DataBC, Province of British Columbia | Bureau of Land Management, Esri, HERE, Garmin, NGA, USGS, NPS, NRCan. Available from: https://maps.gov.bc.ca/ess/hm/imap4m/ Kalra, Y.P., and Maynard, D.G. 1991. Methods manual for forest soil and plant analysis. Edmonton, Alberta: Forestry Canada. doi: NOR-X-319. Kirchen, G., Calvaruso, C., Granier, A., Redon, P-O., Van der Heijden, G., Breda, N., and Turpault, M-P. 2017. Local soil type variability controls the water budget and stand productivity in a beech forest. Forest Ecology and Management. 390: 89-103. 51 Krueger, E.S., Oschsner, T.E., Quiring, S.M., Engle, D.M., Carlson, J.D., Twidwell, D., and Fuhlendorf, S.D. 2017. Measured soil moisture is a better predictor of large growing-season wildfires than the Keetch-Byram Drought Index. Soil Science Society of America Journal. 81: 490-502. Lin, H.S., Kogelmann, W., Walker, C., and Bruns, M.A. 2006. Soil moisture patterns in a forested catchment: a hydropedological perspective. Geoderma 131 (3– 4): 345–68. doi: 10.1016/j.geoderma.2005.03.013. Lloyd, D., Angove, K., Hope, G., and Thompson, C. 1990. A guide to site identification and interpretation for the Kamloops Forest Region. B.C. Min. For., Res. Br., Victoria, B.C., Land Manage. Handb. No. 23. http://www.for.gov.bc.ca/hfd/pubs/docs/lmh/lmh23.htm Mathys, A., Coops, N.C., and Waring, R.H. 2014. Soil water availability effects the distribution of 20 tree species in western North America. Forest Ecology and Management. 313:144-152 Maynard, D.G., and Curran, M.P. 2007. Bulk density measurement in forest soils. In: M.R. Carter and E.G. Gregorich, editors. Soil Sampling and Methods Of Analysis (2nd edition). CRC Press, Boca Raton, FL. P. 863-869. Nelson, D.W., and Sommers, L.E. 1996. Total carbon, organic carbon, and organic matter. In: Sparks, D.L., editors, Methods of Soil Analysis. SSSA Book Ser. 5. SSSA, Madison. p. 961–1010. O‟Brien, R.M. 2007. A caution regarding rules of thumb for variance inflation factors. Qual. Quant. 41: 673-690. doi: 10.1007/s11135-006-9018-6 Poesen, J., and Lavee, H. 1994. Rock fragments in top soils: significance and processes. Elsevier Science: Catena. 23: 1-28. R Core Team. 2017. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Available from: https://www.R-project.org/ Restaino, C.M., Peterson, D.L., and Littell, J. 2016. Increased water deficit decreases Douglas-fir growth throughout western US forests. PNAS. 113(34): 9557-9562. Saxton, K.E., Rawls, W.J., Romberger, J.S., and Papendick, R.I. 1986. Estimating generalized soil water characteristics from texture. Soil Sci. Soc. Am. J., 50: 1031-1036. Saxton, K.E., and Willey, P.H. 2006. The SPAW model for agricultural field and pond hydrologic simulation. Watershed Models. Chapter 17. Pp. 400-435. Saxton, K.E., and Rawls, W.J. 2006. Soil water characteristics estimates by texture and organic matter for hydrologic solutions. Soil Sci. Soc. Am. J., 70: 15691578. 52 Schaap, M.J., Leij, F.J., and van Genuchten, M.T. 2001. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchal pedotransfer functions. Journal of Hydrology. 251: 163-176. Seneviratne, S.I., Corti, T., Davin, E., Hirschi, M., Jaeger, E., Lehner, I., Orlowsky, B., and Teuling, A.J. 2010. Investigating soil moisture-climate interactions in a changing climate: a review. Earth Science Reviews. 99(2010): 125-161. Simpson, D.G. 2000. Water use of interior Douglas-fir. Canadian Journal of Forest Research. 30: 534-547. Soil Classification Working Group. 1998. The Canadian system of soil classification. Agric. And Agri-Food Can. Publ. 1646 (Revised). 187 pp. Thomas, G.W. 1996. Soil pH and soil acidity. p. 475-490. In D.L. Sparks et al. (ed.) Methods of soil analysis. Part 3. Chemical methods. SSSA Book Series 5. SSSA and ASA, Madison, WI. Valentine, K.W.G., Sprout, P.N., Baker, T.E., and Lavkulich, L.M. 1978. The soil landscapes of British Columbia. Ministry of Environment: Resource Analysis Branch. Victoria, BC. Vereecken, H., Maes, J., Feyen, J., and Darius, P. 1989. Estimating the soil moisture retention characteristic from texture, bulk density, and carbon content. Soil Science. 148(6): 389–403. doi:10.1097/00010694-19891200000001. Vyse, A., Ferguson, C., Simard, S.W., Kano, T., and Puttonen, P. 2006. Growth of Douglas-fir, lodgepole pine, and ponderosa pine seedlings underplanted in a partially-cut, dry Douglas-fir stand in south-central British Columbia. The Forestry Chronicle. 82(5): 723-732. Wagner, B., Tarnawski, V.R., Hennings, V., Müller, U., Wessolek, G., and Plagge, R. 2001. Evaluation of pedo-transfer functions for unsaturated soil hydraulic conductivity using an independent data set. Geoderma. 102: 275-297. Wang, T. Hamann, A. Spittlehouse, D. Carroll, C. Locally downscaled and spatially customizable climate data for historical and future periods of North America. PLoS ONE 11(6): e0156720. https://doi.org/10.1371/journal.pone.0156720 Wösten, J.H.M., Pachepsky, Y.A., and Rawls, W.J. 2001. Pedotransfer functions: Bridging the gap between available basic soil data and missing soil hydraulic characteristic. Journal of Hydrology. 251(3/4): 123-150. Zhao, Y., Krzic, M., Bulmer, C.E., and Schmidt, M.G. 2008. Maximum bulk density of British Columbia forest soils from the Proctor test: relationships with selected physical and chemical properties. Soil Sci. Soc. Am. J. 72:442-452. Zhao, Y., Krzic, M., Bulmer, C.E., Schmidt, M.G., and Simard, S.W. 2010. Relative bulk density as a measure of compaction and its influence on tree height. Can. J. For. Res. 40:1724-1735. 53 CHAPTER 3.0: TESTING THE APPLICATION OF SPATIAL INTERPOLATION REGRESSION METHODS FOR PREDICTING FINESCALE FOREST SOIL VARIABILITY. INTRODUCTION Studies that explore fine-scale behavior of soil water have become increasingly important to forest managers as average global temperatures continue to rise and drought conditions are predicted to increase into the future. Because local soil properties influence the infiltration, run-off and redistribution of water inputs, their overall makeup dictates the soil water dynamics of a site, which has immediate effects on tree and plant growth during periods of precipitation shortage (Romano and Palladino 2002). Due to the large influence that soil properties have on soil hydraulic characteristics (Saxton and Rawls 2006), the need to accurately quantify and map the spatial variability in local soil properties within a single catchment site is crucial for determining best management practices and for accurately determining available water throughout the entire landscape. By identifying soil available water on the ground, forest managers can begin to focus silvicultural efforts to maximize water use at a site and to improve seedling performance in areas that may be experiencing significant water stress during extended periods of the growing season. General topography and landscape position affect the variability of certain soil characteristics. Distribution of water, clay particles, ions and minerals are greatly affected by a soils position on the landscape (Jenny 1994; Sinowski and Auerswald 1999; Lybrand and Rasmussen 2015). Previous research has demonstrated that landscape position, and its associated microclimate conditions, are primary controls on soil development and vegetative distribution in semi-arid ecosystems where water is limited during periods of the growing season (Koch et al. 1995; Dahlgren et al. 1997; Lybrand and Rasmussen 2015). Convergent locations typically have higher concentrations of soil organic carbon as well as fine-grained soil materials relative to 54 upslope divergent sites (Applegarth and Dahms 2001; Lybrand and Rasmussen 2015). In addition, soil profiles in convergent locations are typically more developed in relation to limited water resources being concentrated in downslope positions of the landscape (Lybrand and Rasmussen 2015). Furthermore, as you move upwards along an elevation gradient, increased precipitation and temperature fluctuations occur across short distances and can result in a wide-range of soil properties. The greater precipitation, decreased temperature and decreased evapotranspiration demands at higher elevations (Lybrand and Rasmussen 2015) indirectly affects microbial activities and soil organic carbon content (Schaetzl and Thompson 2015). Therefore, the subsequent variation in soils at different landscape positions within a small catchment zone should relate to high variations in soil hydrological processes (Biswas 2014a). The influence of topography on the spatial distribution of soil hydraulic properties of field capacity (FC), permanent wilting point (PWP), and available water capacity (AWC) has been well established in the literature over large areas on coarse scales (Mohanty and Mousli 2000; Romano and Palladino 2002; Obi et al. 2014), but has yet to be tested on small hillslope scale variability within British Columbia. A study conducted by Obi et al. (2014) on coastal plain soils in Nigeria showed that FC water content was significantly correlated (p < 0.05) with slope, compound topographic index (CTI), sand content and clay ratio. Actual AWC (measured as the amount of water in the soil between FC and PWP) was significantly correlated with aspect, CTI and sand contents, while PWP depended on slope, aspect, total curvature, stream power index (SPI) as well as with sand content. Finally, they found that hydraulic conductivity depended on SPI, sand content, silt content, and the silt-clay ratio (Obi et al. 2014). This study, however, was administered over a larger area (8412 km2) than my study site (16 ha), and the quantification of small-scale soil variability in relation to topography is still relatively unknown. Soil mapping methods as well as the accuracy of digital soil maps has been improved by the availability of global positioning systems (GPS), geographical 55 information systems (GIS), remote sensing and spatial statistics. By providing more accurate spatial information at finer scales, these innovations are becoming the wellaccepted standard by pedologists and soil mappers worldwide (Brevik et al. 2016). The concept behind spatial statistics relies on the principle of spatial autocorrelation, which measures the similarity by proximity of soil samples (ie. the farther you move away from a soil sample, the greater the variability will become in a given soil property). The use of regression based spatial statistics for creating digital soil maps (DSM) requires three main components in order to be useful and complete (Minasny and McBratney 2016). First, an input component consisting of field and lab observations of local soil properties either from legacy soil maps or from new field samples is required. Next, the inference process used for the interpolation of soil properties must be identified, and includes the building of models for comparing soil observations with environmental covariates or factors affecting soil variability. Finally, an output is produced in the form of raster maps or predictions along with uncertainty of prediction (Minasny and McBratney 2016). A simple method for interpolating soil properties using digital elevation model (DEM) data is to apply multiple linear regression (MLR) to a dataset. By using MLR alone, soil values are predicted throughout the landscape based on their relationships with values from predictor raster grids created from the DEM. While MLR applies a global method for interpolating soil properties, different geostatistical approaches have been developed to account for local variability in soil properties with respect to their spatial location, and assume that the structure of the soil data will change over space. Therefore, geostatistical methods will allow different relationships within the model to exist at different points in space (Brunsden et al. 1996; Wang et al. 2012). Geographically weighted regression (GWR) is a spatial analysis technique that was specifically designed to account for the spatial heterogeneity of soil properties by measuring the relationship between the dependent and independent variables at certain points, which differ from location to location (Fotheringham et al. 1998). Thus, the model performance when using GWR will vary across the study region and is dependent on local relationships between 56 soil data and the additional topographic variables. In order to compare the methods of MLR and GWR it is important to understand the model components of each. Global regression models such as MLR follow the assumptions outlined by the following equation: ∑ [1] where yi represents a dependent variable, the x are independent variables, α represents a parameter to be estimated, ε represents an error term, and i represents a point in space at which observations on the y and x axis are recorded. GWR builds on this model by allowing local rather than strictly global parameters to be estimated and the model is rewritten as: ∑ where [2] denotes the coordinates of the ith point in space and realization of the continuous function is the at point i. The spatial variability is then captured by allowing a continuous surface of parameter values in relation to surface measurements recorded at the site (Wang et al. 2012). Thus, GWR recognizes that spatial variations in relationships between soil and topography may exist, and provides a method to measure each point (Fotheringham et al. 1998). Previous research using GWR methods for soil spatial analyses have proven useful and accurate for interpolation of soil organic matter (SOM) (Wang et al. 2012; Mishra et al. 2010), soil organic carbon (Song et al. 2016), total nitrogen (Wang et al. 2013), and soil electrical conductivity (Terrón et al. 2012). The overall goal of this study is to quantify the fine-scale soil physicochemical and hydraulic variability of hot, dry Douglas-fir ecosystems within British Columbia using topographic indices as predictor variables in regression based spatial statistics. In doing so, the influence of topographic variables can be assessed for 57 their relationship with the measured and calculated soil parameters. I will test whether or not the variation in the topographic landscape can be used to accurately interpolate measured and derived soil properties in order to create predictions of crucial soil data throughout the entire harvested area. Similar studies for quantifying soil spatial variability have been conducted on much larger areas (Obi et al. 2014; Levi and Rasmussen 2014; Levi et al. 2015) and on similar sized single-order catchments (Holleran et al. 2015). However, measured soil properties have rarely been coupled with soil hydraulic parameters across an entire harvested landscape (Romano and Palladino 2002), and has yet to be assessed within British Columbia. Soil spatial analyses over finer scales are necessary for forest managers to better understand the hydrological and vegetation responses to different soil parameters. Additionally, the statistical modeling of the soil hydraulic parameters with topographic features via regression models is a key aspect to this research as it helps to explain the fundamental patterns of soil water dynamics within recently harvested forests of dry Douglas-fir ecosystems. By doing this, the applicability of different spatial regression methods can also be assessed for their usefulness within fine scale topographic surveys. The specific objectives of this research are to (i) determine the best scale for acute topographic analysis (1, 3 or 5 m) by testing the regression models on specific soil variables while adjusting the kernel radius to allow for different levels of cell smoothing, (ii) quantify soil variability throughout the landscape at 4 depths using multiple linear regression to determine topographic influence on soil properties within different soil layers, (iii) test the usefulness of principle component analysis (PCA) to combine covariate topographic grid layers for use as predictors within the spatial regression models, and (iv) compare different spatial regression methods of MLR with GWR to determine the practicality of each method for predicting fine-scale soil variability within these ecosystems. This research was completed in an effort to improve the understanding of fine-scale soil water dynamics within water limited ecosystems of interior Douglas-fir forests. By combining measured soil properties with these advanced geostatistical methods, I will be able to assess the suitability and accuracy of these techniques for fine-scale 58 forest soil characterization of hot, dry Douglas-fir ecosystems within British Columbia. METHODS Site Description Refer to chapter 2.0 for site description information. Research Design Refer to chapter 2.0 for research design information. The location of the soil pits was noted using a handheld Garmin GPS system and a topographic survey for the site was completed using an aerial light detection and ranging (LiDAR) survey at a 1 m ground resolution and was processed using the System for Automated Geoscientific Analysis (SagaGIS) software (version 4.0.1) (Conrad et al. 2015). Soil Properties Methods for collecting soil properties are fully explained within chapter 2.0. The soil water characteristics equations developed by Saxton and Rawls (2006) were utilized to predict soil hydraulic variables from the measured soil data. The variables of soil texture, SOM, CF concentration and DF were implemented into the model to predict the values of soil saturation, FC, PWP, and AWC (= FC – PWP) (as explained in chapter 2.0). Soil water holding capacity (SWHC) for each soil layer was determined to account for CF content when considering AWC of the soils. The 59 SWHC (mm) for each sample was calculated following the formula used by Kirchen et al. (2017) as: ∑ ( ( )) [3] where RVi represents the volume of coarse fragments in layer i (%), θFCi and θPWPi represent the water values at field capacity and wilting point (m3 m-3), respectively, of soil layer i, and Hi represents the height of the soil layer i (mm). Mapping and Geostatistical Analysis Using the LiDAR digital elevation model (DEM) created in SagaGIS, environmental covariate layers (grids) were developed from the remotely sensed 1 m pixel resolution imagery collected in 2015 (Table 3.1). The geostatistical mapping tool, SagaGIS, was utilized to produce spatial maps of the field site to explain patterns of the measured and derived soil properties. The covariate layers produced for the site were derived from both reflectance and elevation indices and included elevation, positive openness (PO), LS factor (LS), convergent index (CI), catchment area (CA), catchment slope (CS), modified catchment area (MCA), topographic wetness index (TWI), slope, topographic position index (TPI), valley depth (VD), and multiresolution index of valley bottom flatness (MRVBF) (total of 12 grid layers in total) (Table 5.1, Appendix I). The same topographic grid layers were also produced using the simple filter function in SagaGIS, and with a kernel radius of 1 and 2 (giving each 1 m2 cell a smoothing distance of one and two cells in every direction, respectively), to create grid layers at 3 m and 5 m spatial scales. The grid layers at different spatial scales were then tested via linear regression on the variables of sand, clay, total C, BD and pH to determine the best scale for soil predictions within these ecosystems. 60 Table 3.1: Descriptive statistics of the 12 topographic raster grid layers produced in SagaGIS using the LiDARderived digital elevation model for the Isobel Lake study site. CI = Convergence Index; LS = Length Slope Factor; TPI = Topographic Position Index; CA = Catchment Area; CS = Catchment Slope; MCA = Modified Catchment Area; TWI = Topographic Wetness Index; MRVBF = Multi-Resolution of Valley Bottom Flatness; VD = Valley Depth; and PO = Positive Openness. Topographic Variable 1. Elevation 2. Slope 3. CI 4. LS 5. TPI 6. CA 7. CS 8. MCA 9. TWI 10. MRVBF 11. VD 12. PO Unit (m) Radians Radians - Min 1017 0.00 -76.93 0.00 -4.71 1.0 0.01 1.0 -1.56 0.00 0.00 1.15 Max 1066 0.97 93.70 12.87 6.39 880.0 0.86 1100.0 7.16 2.79 14.49 1.59 Mean 1040 0.21 0.03 2.72 0.57 142.0 0.19 166.7 2.12 0.09 1.18 1.48 SD 11.79 0.08 3.23 1.21 1.95 375.4 0.05 454.3 0.64 0.24 2.17 0.04 A correlation matrix with linear regression analysis was created for each measured depth to demonstrate the significant relationships between the sampled and derived soil properties with the 12 individual topographic grid layers developed in SagaGIS. Multiple linear regressions (MLR) were used to test the predictive accuracy of the 12 topographic grid layers combined on the individual soil properties throughout the landscape at the 4 measured depths. The strength of the prediction at each depth is reported as the R2 value of each regression model and was assessed to determine at which depth topography is the most accurate at explaining the apparent soil variability. PCA (for combining grid layers with collinearity) was tested for its usefulness on fine-scale soil predictions by comparing spatial regression models using both PCA combined grid layers and individual grid layers as the predictors for different soil properties. The 12 grid layer values were first standardized within SagaGIS and assigned a z-score as: [4] 61 where represents a raw score for each variable, that grid layer, and represents the mean value of represents the standard deviation. To determine which layers accounted for 90% of the variability within the covariate landscape dataset, the statistical software “R_3.3.3” was utilized to run PCA (using covariance matrix) on the topographic variable data created in SagaGIS. PCA methods were then run in SagaGIS to eliminate the redundant data layers in the dataset and to create distinct environmental covariates (grids) to be used in regression models for predicting variability in soil properties (Levi and Rasmussen 2014, Holleran et al. 2015). The threshold for automated predictor selection for both regression models with and without principle components was set as p = 0.10 within SagaGIS. GWR techniques were processed in SagaGIS using the 12 topographic grid layers as predictor variables and a Gaussian weighting function with a global search distance for interpolating the soil data. A comparison of the resulting coefficient of determination (R2) from each test was used to determine the overall accuracy of these models for predicting the different soil properties. Coefficients of determination (R2) values from the multiple linear regression models were then compared with the R2 values from the GWR models to determine if GWR was successful at increasing the predictive accuracy of the PCA topographic grid layers. Both the MLR and GWR models with PCA predictor grids were cross-validated using a leave-one-out method and comparing the observed versus predicted values at each point on the landscape where the data was removed. A total of 10 runs were performed on each soil variable within each of the 4 measured depths to create the observed versus predicted comparisons, and a different line of data was randomly removed from the full dataset for each cross-validation test. 62 RESULTS Soil Properties Refer to results section of chapter 2.0 for all soil property measurements at the Isobel Lake site. The concentration of available water per cm of soil declined with depth. All of the soil hydraulic variables decreased with depth on average. The SWHC was highest in the 0-15 cm depth range and lowest in the 60-100 cm layer. Calculated SWHC showed significant differences between the surface soil layer and all the lower 3 depths (p-values all < 0.05), however, SWHC was not statistically different between the lower 3 depth ranges themselves (p-values all > 0.05). Additionally, values of FC, PWP and AWC significantly declined with depth (Table 3.2). -1 Table 3.2: Average values (in vol vol ) for the calculated soil hydraulic variables of field capacity (FC), permanent wilting point (PWP), available water capacity (AWC) and soil water holding capacity (SWHC) at 4 different depths throughout the entire study site (n=60 for each depth). Numbers in brackets represent standard error of the mean. Superscript letters represent significant statistical differences (p = 0.05) between mean values using Tukey multiple comparison of means test. Depth Range (cm below soil surface) 0-15 15-30 a -1 FC (vol vol ) -1 PWP (vol vol ) AWC (vol vol-1) SWHC (vol vol-1) 0.332 (0.003) 0.161 (0.002)a 0.171 (0.002)a 0.147 (0.003)a 30-60 b 0.301 (0.003) 0.155 (0.002)a 0.147 (0.003)b 0.0910 (0.004)b 60-100 c 0.275 (0.003) 0.138 (0.002)b 0.137 (0.002)bc 0.094 (0.003)b 0.250 (0.004)d 0.121 (0.003)c 0.130 (0.002)c 0.0850 (0.002)b Mapping and Geostatistical Analysis Testing the regression models using variable kernel radii of 0, 1 and 2 did not appear to have a positive effect on the predictability of the soil properties of sand, clay, carbon, BD or pH within the top 15 cm of soil. Increasing the smoothness of the 63 topographic data slightly increased the R2 value for some variables, while it reduced R2 values for others, but the changes were minor (Figure 3.1). Additionally, increasing the kernel radius did not improve the models at any of the lower measured depths. For this reason, I chose to continue using the original LiDAR data at 1 m scale for the remainder of the topographic regression analyses. Coefficient of Determination (R2) 0.6 0.5 0.4 Sand Clay 0.3 Carbon 0.2 BD pH 0.1 0 0 1 Kernel Radius 2 2 Figure 3.1: Fluctuating kernel radius of topographic grid layers and their associated R values with respect to different soil properties. A total of 12 grid layers were used in the regression models for each individual kernel radius and were used as predictor grids for explaining soil patterns within the upper 15 cm of soil. BD = fine fraction soil bulk density. In total, the 12 topographic grid layers showed significant correlation (at p= 0.10) with a selection of 15 soil properties 228 times throughout the 4 measured depths (Table 3.3). Different topographic grid layers showed significant correlation with the numerous soil properties throughout different depth ranges. Grids that demonstrated the most correlations with the soil data within the upper 15 cm of soil did not necessarily show significant correlation with subsurface values of the same properties. For example, the layers of Elevation, CI, MCA and TPI demonstrated high correlations with soil texture within the 0-15 cm depth, whereas only the grid layer of CI had high correlation with soil texture in the 15-30 and the 30-60 cm depth ranges (Table 3.4). Additionally, the grid layers of Slope and LS showed high 64 correlation with the soil hydraulic properties in the 0-15 cm layer, while PO, TPI and CS demonstrated high correlations with soil water variables in the 15-30 cm layer. Soil hydraulic properties within the 30-60 cm layer were highly correlated with Elevation, MCA, TWI and TPI. Hydraulic properties within the lowest measured layer of 60-100 cm also showed high correlation with Elevation and TWI; however, they also showed highly significant correlations with MRVBF and Slope. Furthermore, BD correlated highest with Elevation and TPI at the soil surface (0-15 cm) and with PO and CS in the 15-30 cm depth, while MCA, TWI and MRVBF were also highly correlated with BD within the bottom 70 cm of soil (Table 3.4). 65 Table 3.3: The number of times (out of a total of 15 tested soil variables) each topographic grid layer showed a significant correlation at the p = 0.10 level with a soil property. DEM Elev = site elevation, CI = catchment index, LS = length-slope factor, TPI = topographic position index, PO = positive openness, CA = catchment area, = CS = catchment slope, MCA = modified catchment area, TWI = topographic wetness index, MRVBF = multiresolution of valley bottom flatness, VD = valley depth. 0 - 15 cm 15 - 30 cm 30 - 60 cm 60 - 100 cm Total DEM Elev Slope 9 8 7 4 7 3 8 7 31 22 CI 4 6 7 7 24 LS 7 0 2 7 16 TPI 8 7 6 3 24 PO 5 8 4 1 18 CA 4 0 6 3 13 CS 3 8 6 3 20 MCA 5 3 7 3 18 TWI MRVBF 2 3 0 2 6 4 6 8 14 17 VD 5 0 1 5 11 Total 63 45 59 61 228 Table 3.4: Correlation matrix showing the relationship between measured and derived soil properties within different depth ranges of the soil profile with topographic grid layers produced in SagaGIS. Bold text and highlighted cells represent significant relationships at the 0.10 level. DEM Elev = site elevation, CI = catchment index, LS = length-slope factor, TPI = topographic position index, PO = positive openness, CA = catchment area, = CS = catchment slope, MCA = modified catchment area, TWI = topographic wetness index, MRVBF = multiresolution of valley bottom flatness, VD = valley depth, C = total carbon, N = total nitrogen, C stock = carbon stock, CF = coarse fragments, BD = fine fraction bulk density, PWP = permanent wilting point, FC = field capacity, AWC = available water capacity, and SWHC = soil water holding capacity. 0 - 15 cm DEM Elev 0.31 -0.22 -0.16 Slope 0.09 0.09 -0.32 CI LS TPI PO -0.11 -0.07 0.53 0.10 0.29 0.18 -0.38 -0.07 -0.32 -0.18 -0.29 -0.08 CA -0.45 0.33 0.22 CS -0.24 0.14 0.19 MCA -0.48 0.34 0.26 TWI MRVBF -0.22 -0.02 0.20 -0.17 0.04 0.32 VD -0.05 0.06 -0.01 C (g kg-1) 0.03 -0.11 0.04 -0.18 0.34 0.18 -0.12 -0.19 -0.19 -0.16 -0.02 -0.14 -1 N (g kg ) C:N -0.29 0.57 -0.07 -0.05 -0.06 -0.08 0.24 -0.15 0.17 0.31 -0.03 0.42 -0.03 -0.15 0.06 -0.44 -0.09 -0.18 -0.01 -0.27 0.08 -0.24 -0.05 -0.17 C Stock (kg m-2) pH CF (% vol.) -0.28 0.05 0.29 -0.02 -0.39 0.16 0.12 -0.04 -0.25 -0.45 0.04 0.14 0.04 0.06 0.37 0.04 0.09 0.10 0.01 -0.31 -0.14 0.02 0.05 -0.35 -0.09 -0.07 -0.25 -0.11 -0.12 0.06 -0.04 0.39 -0.08 -0.04 -0.30 0.05 BD (g cm-3 ) Saturation (% vol.) PWP (% vol.) FC (% vol.) AWC (% vol.) SWHC (% vol.) -0.46 0.46 -0.08 0.17 0.34 -0.06 0.25 -0.25 -0.32 -0.34 -0.22 -0.22 0.09 -0.09 -0.18 -0.11 0.00 -0.04 0.35 -0.36 -0.26 -0.35 0.37 0.26 -0.29 0.11 0.12 -0.36 0.19 0.22 -0.28 0.19 0.23 -0.22 -0.19 0.03 0.17 -0.17 0.03 -0.01 -0.04 0.10 0.21 -0.21 -0.05 -0.13 -0.15 0.19 0.05 -0.05 0.00 0.06 0.09 0.25 0.12 -0.12 -0.11 -0.09 -0.03 -0.04 -0.01 0.01 0.20 0.10 -0.04 0.05 0.27 -0.27 -0.13 -0.23 -0.23 -0.14 Sand Silt Clay 66 15 - 30 cm DEM Elev -0.17 0.21 -0.04 Slope 0.04 0.06 -0.11 CI LS TPI PO -0.22 -0.02 0.21 -0.10 0.40 0.04 -0.04 0.20 -0.19 -0.03 -0.21 -0.12 CA -0.17 0.09 0.10 CS 0.10 -0.16 0.07 MCA -0.16 -0.02 0.22 TWI MRVBF -0.03 0.15 -0.03 -0.22 0.08 0.08 VD 0.05 -0.06 0.02 C (g kg-1) -0.02 -0.12 0.29 -0.11 0.06 0.19 0.01 0.02 0.04 0.01 -0.02 -0.20 N (g kg-1) C:N -0.31 0.38 -0.04 -0.18 0.01 -0.02 0.38 -0.21 0.00 0.15 -0.07 0.42 0.06 -0.10 0.24 -0.33 0.02 -0.02 0.05 -0.08 0.07 -0.10 -0.11 -0.17 C Stock (kg m-2) pH CF (% vol.) -0.04 -0.28 -0.07 -0.18 -0.27 -0.09 0.06 -0.12 -0.05 0.03 -0.46 -0.20 -0.26 -0.22 0.24 -0.07 0.24 0.26 0.17 0.01 -0.12 0.19 0.20 -0.42 0.16 0.23 -0.17 0.03 0.14 0.07 0.00 0.37 0.03 -0.16 -0.15 0.01 BD (g cm-3 ) Saturation (% vol.) PWP (% vol.) FC (% vol.) AWC (% vol.) SWHC (% vol.) -0.24 0.24 -0.05 0.23 0.27 0.13 -0.27 0.27 -0.16 0.19 0.27 0.20 0.09 -0.20 0.22 0.37 -0.09 0.20 -0.22 -0.37 -0.06 -0.07 -0.17 -0.03 -0.06 0.16 -0.29 -0.33 -0.03 0.20 -0.23 -0.33 -0.14 0.16 -0.27 -0.29 0.08 -0.08 0.09 -0.01 -0.06 0.06 -0.31 0.31 0.08 0.30 0.28 0.43 0.02 -0.02 0.22 0.08 -0.01 0.12 0.19 -0.19 0.07 -0.15 -0.19 -0.12 -0.02 0.03 0.07 0.02 -0.01 -0.06 -0.05 0.05 -0.07 0.01 0.05 0.00 DEM Elev -0.19 0.17 0.12 Slope 0.03 0.05 -0.11 CI LS TPI PO -0.23 0.09 -0.13 -0.16 0.22 -0.03 0.14 0.10 0.13 -0.11 0.05 0.15 CA 0.00 -0.04 0.04 CS 0.19 -0.13 -0.17 MCA 0.17 -0.19 -0.06 TWI MRVBF 0.14 0.07 -0.15 -0.08 -0.07 -0.03 VD 0.11 -0.03 -0.14 C (g kg-1) 0.11 -0.23 -0.02 -0.25 0.21 0.20 -0.09 -0.17 -0.14 -0.05 0.06 -0.21 -1 N (g kg ) C:N -0.14 0.31 -0.17 -0.17 -0.33 -0.19 0.32 -0.19 0.14 0.15 -0.06 0.39 -0.10 0.00 0.09 -0.31 -0.14 -0.04 0.01 -0.16 0.19 -0.14 -0.12 -0.26 C Stock (kg m-2) pH CF (% vol.) -0.01 -0.36 -0.13 -0.26 -0.22 0.05 -0.25 -0.27 0.12 -0.01 -0.58 -0.12 -0.38 -0.31 0.29 0.19 -0.21 0.29 -0.06 0.15 0.24 -0.01 0.33 -0.11 -0.10 0.29 0.31 -0.06 0.16 0.34 0.19 0.32 -0.19 -0.14 0.08 -0.01 BD (g cm-3 ) Saturation (% vol.) PWP (% vol.) FC (% vol.) AWC (% vol.) SWHC (% vol.) -0.44 0.44 0.13 0.37 0.45 0.32 -0.21 0.21 -0.16 0.02 0.19 0.06 -0.10 0.10 0.10 0.17 0.18 -0.10 -0.02 -0.44 -0.13 0.03 0.44 0.13 -0.16 0.09 0.18 -0.11 0.33 0.22 0.00 0.43 0.16 -0.14 0.35 -0.11 0.40 -0.40 0.01 -0.22 -0.35 -0.34 0.23 -0.23 -0.18 -0.28 -0.26 -0.06 0.51 -0.51 -0.09 -0.38 -0.51 -0.47 0.40 -0.40 -0.07 -0.30 -0.40 -0.44 0.27 -0.27 -0.01 -0.17 -0.26 0.00 0.12 -0.12 -0.18 -0.21 -0.14 -0.07 DEM Elev -0.33 0.19 0.34 Slope 0.01 0.03 -0.04 CI LS TPI PO -0.06 0.08 -0.01 -0.13 0.10 -0.06 0.08 0.12 0.00 -0.08 -0.06 0.08 CA -0.04 -0.02 0.08 CS 0.05 0.00 -0.08 MCA -0.02 -0.04 0.07 TWI MRVBF 0.14 0.07 -0.11 -0.05 -0.12 -0.06 VD 0.11 -0.17 0.00 C (g kg-1) -0.06 -0.24 -0.38 -0.24 0.05 -0.01 -0.08 -0.04 -0.10 -0.01 0.11 0.02 N (g kg-1) C:N -0.16 0.15 -0.27 -0.05 -0.52 -0.23 0.15 -0.13 0.02 0.10 -0.17 0.29 -0.05 -0.09 0.01 -0.06 -0.04 -0.19 0.12 -0.35 0.27 -0.25 0.03 -0.07 C Stock (kg m-2) pH CF (% vol.) -0.10 -0.17 -0.05 -0.21 -0.13 -0.30 -0.41 -0.17 -0.07 -0.06 -0.39 -0.10 -0.27 -0.14 -0.05 -0.29 0.39 0.14 0.16 0.06 -0.13 0.07 0.24 -0.40 0.00 0.17 -0.08 0.06 0.10 0.12 0.08 0.14 0.39 0.03 0.39 -0.22 BD (g cm-3 ) Saturation (% vol.) PWP (% vol.) FC (% vol.) AWC (% vol.) SWHC (% vol.) -0.34 0.34 0.30 0.41 0.40 0.23 -0.53 0.53 -0.08 0.19 0.44 0.41 -0.26 -0.30 -0.07 0.03 0.26 0.30 0.07 -0.03 -0.08 -0.12 -0.05 0.07 0.09 0.05 0.00 0.06 0.25 0.22 0.07 0.03 0.17 0.29 -0.22 -0.09 0.30 -0.30 0.06 -0.09 -0.24 -0.03 -0.10 0.10 -0.08 0.00 0.07 0.30 0.46 -0.46 0.04 -0.18 -0.38 -0.15 0.63 -0.63 -0.11 -0.39 -0.58 -0.37 0.69 -0.69 -0.04 -0.36 -0.60 -0.51 -0.22 0.22 0.00 0.07 0.12 0.23 Sand Silt Clay 30 - 60 cm Sand Silt Clay 60 - 100 cm Sand Silt Clay 67 Using combinations of the 12 topographic grid layers (at 1 m spatial scale) as independent variables for predicting individual soil properties displayed varying results for different variables at the different depth ranges. The soil properties of sand, silt and clay had the highest statistical correlation with the topographic variables within the upper 15 cm of soil compared to the other 3 measured depths based on the number of times a statistical correlation was found between those variables and the 12 topographic indices (Table 3.4). Furthermore, the soil properties of C and N were not as well predicted as the combined C:N ratio in most cases, and C stock showed higher correlation with topography than total C concentration alone throughout the top 60 cm of soil. Predictions of soil pH and CF concentration both improved in the lower depths compared to the top 15 cm of soil. Soil BD was predicted by topography with the highest accuracy (based on R 2) compared to all other tested soil properties. The topographic variables explained greater than 60% of the BD variation in the 15-30 and 60-100 cm depth ranges, 42% in the top 15 cm of soil, and 31.5% in the 30-60 cm range. The topographic variables were also well correlated with the derived soil hydraulic properties of AWC and SWHC, and the prediction strength of each variable was improved in the lower depth ranges compared with the top 15 cm of soil. Hydraulic properties within the depth range of 15-30 cm showed the highest correlation with topography compared with the other 3 tested depths. The number of principle components necessary to account for 90% of the variation within the DEM were retained for further geostatistical analyses and mapping of soil properties. A total of 7 distinct principle components comprised of the 12 topographic grid layers were found to explain 90% of the landscape variance (cumulative proportion of variance explained = 0.925) (Table 3.5). 68 Table 3.5: Principle component analysis loadings for 7 components identified to explain over 90% of the landscape dataset variance. The 7 principle components were derived from a combination of the 12 individual topographic grid layers of DEM Elev = site elevation, CI = catchment index, LS = length-slope factor, TPI = topographic position index, PO = positive openness, CA = catchment area, = CS = catchment slope, MCA = modified catchment area, TWI = topographic wetness index, MRVBF = multi-resolution of valley bottom flatness, and VD = valley depth. Grid Layer Explained proportion of landscape dataset variance: PC1 PC2 PC3 PC4 PC5 PC6 PC7 0.325 0.238 0.118 0.078 0.069 0.058 0.039 Total = 0.925 Combining the multiple topographic grid layers through PCA did not have a positive effect on the predictive accuracy of many of the tested models when running MLR models (Table 3.6). The coefficients of determination (R2 values) from the MLR tests were higher for most of the soil variables at each depth compared with those created using the PCA grids rather than individual raster grids. However, in some instances, the coefficient of determination was the same or slightly increased when using the 7 PCA grids. Interestingly, the models had the hardest time explaining the variability in sand, silt and clay, and were better at explaining variability in pH, CF, BD and AWC. 69 Table 3.6: Resulting coefficients of determination for multiple linear regression (MLR) models using individual grid layers as predictor variables (12 total), and principle components as predictor variables (7 total), for predicting different soil properties within 4 different depths. Total C = total carbon concentration, Total N = total nitrogen concentration, C Stock = carbon stock, LCF = lab tested coarse fragment concentration, BD = fine fraction bulk density, AWC = available water capacity and SWHC = soil water holding capacity. 0-15 cm Grids PCA 0.299 0.298 0.237 0.197 0.386 0.238 Sand Silt Clay Depth (R2) 15-30 cm 30-60 cm Grids PCA Grids PCA 0.175 N/A 0.0380 N/A 0.144 0.095 0.0330 N/A 0.032 N/A N/A N/A Total C (g kg-1) 0.130 0.0390 0.0660 0.149 -1 Total N (g kg ) C:N 0.230 0.316 N/A 0.212 0.115 0.225 0.0450 0.0920 0.198 0.168 C Stock (kg m-2) pH CF (% vol.) 0.196 0.249 0.282 N/A 0.320 0.287 BD (g cm-3) AWC (%) SWHC (%) 0.420 0.233 0.231 0.302 0.0660 0.079 0.0460 0.0500 60-100 cm Grids PCA 0.134 N/A 0.0190 N/A 0.188 N/A 0.156 0.0990 0.108 0.185 0.261 0.241 0.245 0.143 0.167 0.378 0.403 0.0910 0.0930 0.0600 0.302 0.458 0.395 0.147 0.332 0.382 0.153 0.271 0.405 0.1000 0.284 0.158 0.674 0.673 0.482 0.376 0.336 0.192 0.667 0.555 0.463 0.604 0.454 0.272 0.315 0.311 0.286 0.338 0.284 0.298 Geographically weighted regression increased the coefficient of determination (R2) for each of the tested soil variables within the upper 30 cm of soil (Figure 3.2). Additionally, the method of applying GWR for interpolating soil properties performed better than MLR alone throughout all 4 measured depths when testing the 7 PCA grid layers as predictor variables. Cross validation of both the MLR and GWR models were sufficient and displayed R2 values close to those found when running the model using all 60 data values (Figure 3.3; Table 5.2, Appendix I). The 7 PCA grids were most useful for interpolating BD overall within the 4 measured depths, and cross-validation displayed a general agreement between the measured values from the field and the predicted values created by the MLR and GWR models (Figures 3.3 and 3.4). This result demonstrates that running PCA on the raster grid layers in combination with GWR may be the optimal method for fine-scale soil 70 predictions for this site, and especially for predictions of BD. Although the R2 value of GWR models with the individual grid layers were higher than those of the GWR models using PCA grids in many instances, the models including the PCA grids are likely more robust and do not have collinearity between the topographic variables. MLR (PCA Grids) GWR (PCA Grids) Coefficient of Determination (R2) 0.50 0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.00 Sand Clay 2 SOM CF BD SWHC Figure 3.2: Comparison of R values for MLR and GWR methods using the 7 PCA grids as predictor variables for interpolating various soil properties within the upper 15 cm of soil. SOM = soil organic matter, CF = coarse fragment concentration, BD = fine fraction bulk density and SWHC = soil water holding capacity. 71 Figure 3.3: Leave-one-out cross validation results comparing the multiple linear regression (MLR) and -3 geographically weighted regression (GWR) models for interpolating the values of fine fraction bulk density (g cm ). A = depth 1 (0-15 cm), B = depth 2 (15-30 cm), C = depth 3 (30-60 cm) and D = depth 4 (60-100 cm). Predictor -3 variables = 7 PCA grids. Tested soil property = fine fraction bulk density (BD) (g cm ). 72 Figure 3.4: Output map of the interpolated values of fine fraction bulk density (BD) (0-15 cm) using Geographically Weighted Regression (GWR) and the 7 topographic grid layers produced through principle component analysis (PCA) as independent, predictor variables. 73 DISCUSSION Interpolated Soil Variability The results from this study exemplify the usefulness of topographic raster data for the interpolation of numerous measured soil data on finer scales than previously studied within forested ecosystems. By successfully modeling the various soil properties responsible for dictating AWC as well as the SWHC throughout the site, important information about which topographic features may be influencing water patterns on a fine-scale was discovered. Although many of the relationships may not be representative of cause and effect, these results do provide insight into which topographic features should be further explored to determine their ultimate influence on the variability in AWC. Significant relationships were revealed for nearly all measured and derived soil variables when using independent topographic variables as predictors and using both MLR and GWR methods. The models produced using GWR exceeded the accuracy of MLR models alone. However, in terms of practical forest management implications within these ecosystems, the prediction accuracy of soil water resources across the study site may not be particularly useful. Due to the relatively low (~ 0.50 or less) coefficient of determination (R2) values created by the various models for many of the soil properties it becomes questionable as to whether forest managers can make ground based decisions based on models that explain roughly 50% or less of the variation in the soil data. For example, the models created for predicting sand and clay concentrations and SOM explained less than 50% of the measured variability within all four depth ranges. When considering the location of specific site prescription treatments on the ground that attempt to improve tree seedling performance, a 50/50 guess of local soil characteristics at a given site is not very reliable. Nonetheless, the models produced for BD, LCF and SWHC resulted in higher R2-values, which revealed important relationships between these soil properties and certain topographic features, and can help improve forest management within the hot, dry 74 Douglas-fir ecosystems within BC in the future. Being able to identify areas of abundant coarse fragments or increased compaction based on digital elevation models of a site could improve applications of site preparation to maximize soil water use because soil bulk density (Bulmer and Simpson 2005; Gupta and Larson 1979) and coarse fragment concentration (Poesen and Lavee 1994; Liu and She 2017; Zhang et al. 2016; Ma et al. 2010) have such a large role in dictating soil available water. In terms of overall site variability, the results from this study demonstrate how many of the measured soil properties varied little across the site. This is somewhat surprising as it is generally believed that small scale variability commonly exists within small geographic areas, and is often attributed to the sites heterogenous relationship with the five-main soil forming factors outlined by Jenny (1994). Excluding coarse fragment concentration, the soils at the Isobel site generally showed less variability in physical properties (ie. sand, silt, clay, BD) compared with chemical properties (ie. carbon, nitrogen, C:N ratio). This finding is consistent with results from similar studies which look at the fine-scale variability of physical and chemical soil properties (Amador et al. 2000). However, the large variability in coarse fragment concentration found throughout the Isobel site is a notable discovery as this soil property has a strong influence on the water holding capacity of soils and can drastically alter the growing conditions of a microsite within water limited ecosystems (Cousin et al. 2003; Baetens et al. 2009). Despite the apparent low variability in some soil properties, certain terrain attributes were highly correlated with soil properties at different depths of the profile, and were consistent with findings of other soil spatial analyses conducted within much larger areas. Romano and Palladino (2002) found that the inclusion of slope gradient and slope aspect both improved the prediction strength of soil hydraulic parameters within the topsoil of a 32-square kilometer drainage in Italy. Florinsky et al. (2002) also discovered slope gradient to hold significant correlations with soil water within the 0-30 cm depth range throughout the growing season of two consecutive years, however, this relationship was not consistent at lower depths and 75 did not exhibit statistical correlations beyond 30 cm up to 1.2 m. This result is somewhat intuitive because as slope increases, the velocity of water flow also increases, so the rainfall that is received and infiltrated per unit area on the ground is decreased, while runoff area increases and re-distributes the water to lower slope locations where it either infiltrates or evapourates. Thus, increased slope gradient results in decreased overall soil moisture (Florinsky et al. 2002). Similarly, the topographic covariate grids of slope and LS-factor (a function of slope length) in my study significantly correlated with the soil properties of PWP, FC, AWC and SWHC within the 0-15 cm depth range at the Isobel Lake site. Florinsky et al. (2002) also found that the covariate layer representing catchment area held significant relationships with soil water throughout the entire 1.2 m soil profile. However, significant correlations between soil AWC and SWHC with catchment area within the 30-60 cm depth range and with only AWC in the 60-100 cm range were present in my study. Furthermore, Florinsky et al. (2002) did not discover any R2 values for soil property predictions greater than 0.37. Similarly, many of the R2 values found in my study for soil hydraulic predictions were at or below 0.37. However, depending on the topographic variable and the depth range explored, R2 values exceeded 0.40 and sometimes even greater than 0.50 (Table 3.7). Although it should be noted that the DEM used in Florinsky‟s study was produced at a 15-m resolution using a GPS technique, whereas mine was produced at a 1-m resolution using highly detailed LiDAR data, which could explain the improved accuracy of some of my regression models. In addition, the covariate grid layer representing site elevation was also selected for the prediction of SWHC within the lower 3 measured depths (15-30, 3060 and 60-100 cm) and produced R2 values greater than those within the top 15 cm of soil. This result may be due to random noise that can occur within the upper layers of a soil, as factors such as wind and rain erosion, animal disturbance, as well as potential machine disturbance during logging operations have more effect on altering the surface soil layers compared to subsoil layers. Certain covariate grid layers that describe hydraulic tendencies on the ground, including TWI and MCA, were expected to be adequate predictors of the derived soil hydraulic properties at this site. Both TWI and MCA did not correlate well 76 with AWC or SWHC within the top two measured depth layers, however, they did exhibit better relationships within the 30-60 and 60-100 cm depth ranges (Table 3.4). Additionally, the covariate TWI was selected (via stepwise regression) for use within the prediction models for SWHC within the 0-15, 15-30 and 60-100 cm depth ranges, and MCA was selected for the prediction of SWHC within the 0-15 and 3060 cm depth ranges, which supports their usefulness when determining the topographic influence on manipulating soil water resources. In general, coefficients of determination of 0.5 or less are typical for the prediction of soil properties (Ryan et al. 2000; Herbst et al. 2006), which is consistent with most of the predictions made by the different regression models within this study. Due to the high number of sampling points within a relatively small area, it was expected that topography would be able to interpolate the observed soil properties accurately. However, one reason as to why most of the models were below 50% accurate may be a result of the apparent variability in measured soil properties far outweighing the variability in topography. The study site at Isobel Lake covers a relatively uniform hill slope with very minor peaks and gullies, meaning that most of the topographic raster data was similar across the site. Many of the measured soil properties, however, were not found to be uniform and did not necessarily follow a particular pattern with topography, making it difficult for the model to interpolate. The specific soil properties of clay, silt and sand concentrations did not vary considerably between samples (as indicated by the standard deviations of each averaged variable as well as the fact that roughly 90% of the tested soils were classified as a Loam textural class). However, the properties of total C, total N, BD and CF concentration did show considerable amounts of variability throughout the site, despite the little variability in topographic features. The ability of the different regression models to more accurately predict the properties of CF, BD and total C as compared to soil texture is logical as these soil properties showed considerable amounts of short-range variability and thus likely varied more with slight topographic differences. For the most part, models which predicted soil texture (sand, silt and 77 clay), total C and total N had higher coefficients of determination within the topsoil than in the subsoil. In contrary, coefficients for the soil properties of pH, CF, BD, AWC and SWHC were improved in the lower soil layers compared to the topsoil layer. This may have occurred because the influence of topography is greater in the subsoil, while in the topsoil other forms of influence (eg. anthropogenic influence, wind exposure, rain erosion etc.) can be more prominent and can create random spatial distributions that are not well related to topography (Herbst et al. 2006). Effect of Increased Kernel Radii on Soil Interpolations Previous studies have shown that soil data often have a higher proportion of variation at short distances than landform surfaces (Burrough 1983; Gessler et al. 1995; Oueslati et al. 2013). Because soil variability commonly outweighs the variability in landscape data on small scales, adjusting the kernel radius to smooth the topographic raster data was not expected to improve the prediction accuracy of the regression equations. With the LiDAR data projected at a 1 m spatial resolution, most of the landscape variation is captured and will contain the maximum amount of detail possible compared to using a kernel radius of 2 m or 3 m which will smooth out the difference between individual grid cells, and, therefore will lose some of the specific topographic detail. Kuo et al. (1999) discovered that increasing the grid cell size of a DEM misrepresented the curvature of the landscape, therefore resulting in less accurate soil moisture predictions during dry seasons within a watershed in central New York. Because soil properties often have higher variation over shorter ranges than topography does, the added detail within 1 m digital elevation model is important for model fitting and for understanding the relationship between the distributions of certain soil properties in relation to minor changes in topography. 78 Evaluation of Methodologies The use of PCA methods within soil spatial mapping studies is a data-driven approach to reduce the dataset and to determine important covariate landscape variables (Levi and Rasmussen 2014). This approach selects the covariate data that account for the greatest range in variability within the landscape raster data, and eliminates redundant information from the models (Holleran et al. 2015). By combining the different covariate layers and reducing redundancy within the landscape dataset, the resulting coefficient of determination (R2) value was reduced for many of the tested soil properties when compared with the models using the individual grid layers themselves. Although the resulting R2 value was reduced for many model predictions, the models using the PCA grid layers are likely more robust and, therefore, more reliable when considering the ability for landscape features to predict various fine-scale soil properties because PCA grid layers account for collinearity among predictor variables. If the landscape data is not reduced to account for collinearity, then multiple raster grid layers that share similar topographic information may be selected for predictions as they will improve the output prediction accuracy of the regression models. However, the regressions do so by adding topographic information that is repeated within other selected grid layers thereby increasing error (Levi and Rasmussen 2014). Another reason why the models including PCA grid layers did not perform as well as models using individual grids for certain soil properties could be due to the scale at which the landscape data were collected. The LiDAR data from the Isobel Lake site include a very high level of detail when measured at a 1 m spatial scale, and by combining multiple grids through PCA, some important information from the complete dataset is lost. The detail that is eliminated from the PCA procedure could help to explain the measured variability of certain soil properties on a fine-scale. The fact that GWR outperformed MLR methods for soil property interpolations is not surprising, as GWR considers the non-stationary assumption with regards to soil-topography relationships. When attempting to model the spatial distribution of soil properties within small landscapes and on fine-scales similar to the Isobel Lake 79 study site, the concept of heterogeneous spatial relationships between soil properties and topographic relief is important to consider. Because the soils measured at the Isobel Lake site showed considerable short-range variability in some properties that directly influence soil water availability (ie. coarse fragment concentration, carbon concentration and fine fraction bulk density), understanding their intrinsic relationship with certain topographic parameters at different locations throughout the landscape is fundamental for accurate representations. As a result of highly variable soil properties on such a fine scale, the practice of applying a single set of regression parameters (ie. through MLR) to all soil variables across an entire hillside (that is otherwise relatively uniform in topographic shape) likely did not capture the full extent of the relationship between the environmental covariates and soil parameters. This likely resulted in a lower coefficient of determination. CONCLUSION/SUMMARY In general, GWR outperformed the MLR models in all cases, and the use of PCA derived environmental covariate layers produced more robust and reliable models than the models which utilized the individual covariates themselves. Neither MLR nor GWR methods overestimated the soil properties. Interpolated residuals for both methods show a relative agreement with soil values that are possible for this site based on the observed data. However, in terms of practical forest management objectives, the utilization of topographic grid layers for predicting areas of moisture stress and surplus on the ground is still relatively risky when considering that most models explained roughly 50% or less of the measured soil variability at this site. Nonetheless, certain topographic features had stronger relationships than others with specific soil properties and the importance of those topographic variables can be further studied in future research projects. 80 Foresters tasked with managing the hot, dry subzones of the interior Douglasfir forests surrounding Kamloops face considerable challenges as average temperatures are projected to climb and drought conditions are expected to increase into the future. With respect to tree seedling regeneration and sustainable harvest rotations, the ability to predict and quantify available soil water resources throughout proposed harvested areas will become increasingly desirable to ensure that best management practices are applied to specific areas. By combining geostatistics with readily available topographic raster data, certain relationships between soil properties and topography were revealed for these ecosystems. As highly-detailed topographic data for the province of BC becomes more and more accessible to forest managers, the value of identifying these relationships on both coarse and finescales will save time and money when conducting soil prediction surveys into the future. Further research could include applying a combined GWR-kriging approach for fine-scale soil property predictions, as this method has been shown to be successful when tested on simulated datasets compared to normal GWR (Harris et al. 2010). Additionally, a more systematic sampling design could be applied to another site within the IDFxh2 to ensure that soil pit locations capture the entire range of possible topographic distinctiveness, whereas the location of my soil pits were determined with the intention of characterizing soil properties across the landscape, but failed to consider different extremes in topographic relief. LITERATURE CITED Amador, J.A., Wang, W., Savin, M.C., and Görres, J.H. 2000. Fine-scale spatial variability of physical and biological soil properties in Kingston, Rhode Island. Geoderma. 98(1-2): 83-94. Applegarth, M.T., and Dahms, D.E. 2001. Soil catenas of calcareous tills, Whiskey Basin, Wyoming, USA. Catena. 42: 17-38. Baetens, J.M., Verbist, K. Cornelis, W.M., Gabriels, D., and Soto, G. 2009. On the influence of coarse fragments on soil water retention. Water Resources Journal. 45(7). doi: https://doi.org/10.1029/2008WR007402 81 Biswas, A. 2014a. Landscape characteristics influence the spatial pattern of soil water storage: similarity over times and at depths. Catena 116:68–77. Doi: 10.1016/j.catena.2013.12.004. Biswas, A. 2014b. Scaling analysis of soil water storage with missing measurements using the second-generation continuous wavelet transform. European Journal of Soil Science. 65(4): 594–604. doi:10.1111/ejss.12145. Boehner, J., Koethe, R. Conrad, O., Gross, J., Ringeler, A., and Selige, T. 2002. Soil regionalisation by means of terrain analysis and process parameterisation. In: Micheli, E., Nachtergaele, F., Montanarella, L. [Ed.]: Soil Classification 2001. European Soil Bureau, Research Report No. 7, EUR 20398 EN, Luxembourg. pp.213-222. Boehner, J., and Selige, T. (2006). Spatial prediction of soil attributes using terrain analysis and climate regionalization. In: Boehner, J., McCloy, K.R., Strobl, J. SAGA - analysis and modelling applications. Goettinger Geographische Abhandlungen. 115:13-27. Brevik, E.C., Calzolari, C., Miller, B.A., Pereira, P., Kabala, C., Baumgarten, A., and Jordán, A. 2016. Soil mapping, classification, and pedologic modeling: history and future directions. Geoderma. 264: 256–74. doi: 10.1016/j.geoderma.2015.05.017. British Columbia Ministry of Forests and Range and BC Ministry of Environment. 2010. Land Management Handbook – 25: Field manual for describing terrestrial ecosystems – 2nd edition [online]. Research Branch [accessed December 8, 2016]. Available from: https://www.for.gov.bc.ca/hfd/pubs/Docs/Lmh/Lmh25/Lmh25_ed2_(2010).pdf Brunsden, C., Fotheringham, S., and Charlton, M. 1996. Geographically weighted regression: a method for exploring spatial nonstationarity. Geographical Analysis. 28(4): 281-298. Bulmer, C.E., and Simpson, D.G. 2005. Soil compaction and water content as factors affecting the growth of lodgepole pine seedlings on sandy clay loam soil. Can. J. Soil Sci. 85: 667-679. Burrough, P.A. 1983. Multiscale sources of spatial variation in soil. I. The application of fractal concepts to nested levels of soil variation. European Journal of Soil Science. 34(3): 577-597. Conrad, O. 2012. SAGA-GIS Module Library Documentation (v2.3.0). System for Automated Geoscientific Analysis. Accessed: April 24, 2018. Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., Wehberg, J., Wichmann, V., and Böhner, J. 2015. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4, Geosci. Model Dev., 8, 1991-2007, doi:10.5194/gmd-8-1991-2015. 82 Cousin, I., Nicoullaud, B., and Coutadeur, C. 2003. Influence of rock fragments on the water retention and water percolation in a calcareous soil. Catena. 853: 97–114. Dahlgren, R.A., Boettinger, J.L., Huntington, G.L., and Amundson, R.G. 1997. Soil development along an elevational transect in the western Sierra, Nevada, California. Geoderma. 78: 207-236. Florinsky, I.V., Eilers, R.G., Manning, G.R., and Fuller, L.G. 2002. Prediction of soil properties by digital terrain modelling. Environmental Modelling and Software. 17: 295-311. Fotheringham, A. S., Charlton, M. E., and Brundson, C. 1998.Geographically weighted regression: a natural evolution of the expansion method for spatial data. Environment and Planning A. 30:1905–1927. Gallant, J.C., and Dowling, T.I. 2003. A multiresolution index of valley bottom flatness for mapping depositional areas. Water Resources Research. 39(12):1347-1359 Gessler, P.E., Moore, I.D., McKenzie, N.J., and Ryan, P.J. 1995. Soil-landscape modelling and spatial prediction of soil attributes. Int. J. Geographical Information Systems. 9(4): 421-432 Guisan, A., Weiss, S.B., and Weiss, A.D. 1999. GLM versus CCA spatial modeling of plant species distribution. Plant Ecology 143: 107-122. Gupta, S. C., and Larson, W. E. 1979. Estimating soil water retention characteristics from particle size distribution, organic matter percent, and bulk density. Water Resources Research. 15(6): 1633-1635. Harden, J.W. 1982. A quantitative index of soil development from field descriptions: examples from a chronosequence in central California. Geoderma. 28 (1): 1– 28. doi:10.1016/0016-7061(82)90037-4. Harris, P., Fotheringham, A.S., Crespo, R. Charlton, M. 2010. The use of geographically weighted regression for spatial prediction: an evaluation of models using simulated data sets. Math Geoscience. 42: 657-680. Hebda, R. 2007. Biodiversity: geological history in British Columbia. Biodiversity BC Technical Subcommittee. Available from: http://www.biodiversitybc.org/assets/Default/BBC%20Biodiversity%20and%20 Geological%20History.pdf Herbst, M., Diekkrüger, B., and Vereecken, H. 2006. Geostatistical coregionalization of soil hydraulic properties in a micro-scale catchment using terrain attributes. Geoderma. 132: 206-221. 83 Holleran, M. Levi, M., and Rasmussen, C. 2015. Quantifying soil and critical zone variability in a forested catchment through digital soil mapping. SOIL Journal. (1): 47-64. iMapBC. 2018. DataBC, Province of British Columbia | Bureau of Land Management, Esri, HERE, Garmin, NGA, USGS, NPS, NRCan. Available from: https://maps.gov.bc.ca/ess/hm/imap4m/ Jenny, H. 1994. Factors of soil formation. A system of quantitative pedology. Geoderma. Vol. 68. doi:10.1016/0016-7061(95)90014-4. Kalra, Y.P., and Maynard, D.G. 1991. Methods manual for forest soil and plant analysis. Edmonton, Alberta: Forestry Canada. doi:NOR-X-319. Kirchen, G. Calvaruso, C., Granier, A., Redon, P-O., Van der Heijden, G.,Breda, N., and Turpault, M-P. 2017. Local soil type variability controls the water budget and stand productivity in a beech forest. Forest Ecology and Management. 390: 89-103. Koch, G.W., Vitousek, P.M., Steffen, W.L., and Walker, B.H. 1995. Terrestrial transects for global change research. Vegetatio. 121(1/2): 53-65. Koethe, R. and Lehmeier, F. 1996. SARA - System zur automatischen reliefanalyse. User Manual, 2. Edition [Dept. of Geography, University of Goettingen, unpublished] Kolmogorov, A.N. 1941. The local structure of turbulence in an incompressible fluid at very large Reynolds numbers. Doklady Academii Nauk SSSR 30, pp. 301– 305. Kuo, W-L., Steenhuis, T.S., McCulloch, C.E., Mohler, C.L., Weinstein, D.A., DeGloria, S.D., and Swaney, D.P. 1999. Effect of grid size on runoff and soil moisture for a variable-source-area hydrology model. Water Resources Research. 35(11): 3419-3428. Levi, M., and Rasmussen, C. 2014. Covariate selection with iterative principle component analysis for predicting physical soil properties. Geoderma. 219220: 46-57 Levi, M.R., Schaap, M.G., and Rasmussen, C. 2015. Application of spatial pedotransfer functions to understand soil modulation of vegetation response to climate. Vadose Zone Journal. doi:10.2136/vzj2014.09.0126 Lichter, J., and Costello, L. 1994. An evaluation of volume excavation and core sampling techniques for measuring soil bulk density [online]. Journal of Arboriculture. 20(3). Available from: file://mytru_vfiler1/home/T/T00052707/Downloads/p0160-0164.pdf 84 Liu, D., and She, D. 2017. Can rock fragment cover maintain soil and water for saline- sodic soil slopes under coastal reclamation? Elsevier Science: Catena. 151: 213-224. Lloyd, D., Angove, K., Hope, G., and Thompson, C. 1990. A Guide to Site Identification and Interpretation for the Kamloops Forest Region: Land Managment Handbook Number 23. Kamloops, BC: BC Ministry of Forests Research Branch. https://www.for.gov.bc.ca/hfd/pubs/docs/lmh/lmh23.htm. Lybrand, R.A., and Rasmussen, C. 2015. Quantifying climate and landscape position controls on soil development in semiarid ecosystems. Soil Science Society of America Journal. 79 (1): 104. doi:10.2136/sssaj2014.06.0242. Ma, D., Shao, M., Zhang, J., and Wang, Q. 2010. Validation of an analytical method for determining soil hydraulic properties of stony soils using experimental data. Geoderma. 159: 262-269. Minasny, B., and McBratney, A.B. 2016. Digital soil mapping: a brief history and some lessons. Geoderma. 264: 301–11. doi:10.1016/j.geoderma.2015.07.017. Mishra, U., Lal, R., Liu, D., and van Meirvenne, M. 2010. Predicting the spatial variation of the soil organic carbon pool at a regional scale. Soil Sci. Soc. Am. J. 74(3): 906-914. Mohanty, B., and Mousli, Z. 2000. Saturated hydraulic conductivity and soil water retention properties across a soil-slope transition. Water Resources Research. 36(11): 3311-3324. Obi, J., Ogban, P., Ituen, U., and Udoh, B. 2014. Development of pedotransfer functions for coastal plain soils using terrain attributes. Catena. 123(2014): 252-262. O'Callaghan, J.F., and Mark, D.M. 1984. The extraction of drainage networks from digital elevation data. Computer Vision, Graphics and Image Processing. 28:323-344 Ouesleti, I., Allamano, P., Bonifacio, E., and Claps, P. 2013. Vegetation and topographic control on spatial variability of soil organic carbon. Pedosphere. 23(1): 48-58 Poesen, J., and Lavee, H. 1994. Rock fragments in top soils: significance and processes. Elsevier Science: Catena. 23: 1-28. R Core Team. 2016. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. Available from: https://www.R-project.org/ Romano, N., and Palladino, M. 2002. Prediction of soil water retention using soil physical data and terrain attributes. Journal of Hydrology. 265(2002): 56-75. 85 Romano, N., and Chirico, G. 2004. The role of terrain analysis in using and developing pedotransfer functions. Developments in Soil Science. 30: 273294. Ryan, P.J., McKenzie, N.J., O‟Connell, D., Loughhead, A.N., Leppert, P.M., Jacquier, D., and Ashton, L. 2000. Integrating forest soils information across scales: spatial prediction of soil properties under Australian forests. Forest Ecology and Management. 138: 139-157. Saxton, K.E., and Willey, P.H. 2006. The SPAW model for agricultural field and pond hydrologic simulation. p. 401–435. In V.P. Schaetzl, R.J., and Thompson, M.L. 2015. Soils: Genesis and Geomorphology. 2nd ed. New York: Cambridge University Press. doi:978-1-107-01693-4. Sinowski, W., and Auerswald, K. 1999. Using relief parameters in a discriminant analysis to stratify geological areas with different spatial variability of soil properties. Geoderma. 89 (2000): 113–28. doi:10.1016/S00167061(98)00127-X. Song, X., Brus, D., Liu, F., Li, D., Zhao, Y., Yang, J., and Zhang, G. 2016. Mapping soil organic carbon content by geographically weighted regression: A case study in the Heihe River Basin, China. Geoderma. 261: 11-22. Terrón, J., da Silva, J.M., Moral, F., and García-Ferrer, A. 2011. Soil apparent electrical conductivity and geographically weighted regression for mapping soil. Precision Agriculture. 12: 750-761. Thomas, G.W. 1996. Soil pH and soil acidity. p. 475-490. In D.L. Sparks et al. (ed.) Methods of soil analysis. Part 3. Chemical methods. SSSA Book Series 5. SSSA and ASA, Madison, WI. Travis, M.R., Elsner, G.H., Iverson, W.D., and Johnson, C.G. 1975. VIEWIT: computation of seen areas, slope, and aspect for land-use planning. USDA F.S. Gen. Tech. Rep. PSW-11/1975, 70p. Berkeley, California, U.S.A. Wang, W. 2002. Uncertainty analysis of groundwater flow and solute transport in unsaturated- saturated soil: Maricopa case. Ph.D. diss. University of Arizona, Tucson. Wang, T. Hamann, A., Spittlehouse, D.L., and Murdock, T. 2012. ClimateWNA high-resolution spatial climate data for western North America. Journal of Applied Meteorology and Climatology. 51: 16-29 Wang, K., Zhang, C., and Li, W. 2012. Comparison of geographically weighted regression and regression kriging for estimating the spatial distribution of soil organic matter. GIScience and Remote Sensing. 49(6): 915-932. 86 Wang, K. Zhang, C., and Li, W. 2013. Predictive mapping of soil total nitrogen at a regional scale: A comparison between geographically weighted regression and cokriging. Applied Geography. 42: 73-85. Yokoyama, R., Shirasawa, M., and Pike, R.J. (2002). Visualizing topography by openness: A new application of image processing to digital elevation models. Photogrammetric Engineering and Remote Sensing. 68: pp.251-266. Zhang, Y., Zhang, M., Niu, J., Li, H., Xiao, R., Zheng, H., and Bech, J. 2016. Rock fragments and soil hydrological processes: significance and progress. Catena. 147: 153-166 Zhufeng, F. 2009. Using geostatistics, pedotransfer functions to generate 3D soil and hydraulic property distributions for deep vadose zone flow simulations. MSc. Thesis in Hydrology. University of Arizona. 87 CHAPTER 4.0: SIGNIFICANCE AND CONCLUSION The apparent variability in the measured soil properties within the 16-ha site was not as large as was originally expected. Most soil physical properties generally displayed lower fine-scale variation than chemical properties, with the exception of coarse fragment concentration – which had the highest level of variability out of all measured properties. Coarse fragments have a large effect on the amount and distribution of soil water availability. Because of the heterogeneous concentration of coarse fragments within the soil across the site, the variability in soil available water will likely follow a similar pattern. Although soil properties, including organic matter content, soil texture and bulk density also affect water availability within a soil profile, their degree of influence in comparison to coarse fragment concentration is relatively low, and is outweighed by the water loss attributed to coarse fragment volume. In terms of management objectives within IDFhx2 ecosystems, the overall variation in soil water resources did not necessarily follow strong predictive patterns with topography, and did not change significantly from one end of the study site to the other. The results from the power analyses, as well as the general lack of AWC‟s predictable variability based on topographic features demonstrates that highly intensive soil sampling at any given site may be required if forest managers want to have a true understanding of the entire water holding capabilities of the soil. My study demonstrates that even when sampling on relatively uniform topography, the amount of variation in coarse fragments and available water far exceeds what can be captured in one or two pits, and sites that are comparatively heterogeneous in topography are very likely going to require even more samples to gain an accurate representation. While many of the regression equations between soil water availability and topographic variables held relatively weak relationships, the models did provide insight into topographic covariates that may be influencing water availability on fine 88 scales. Despite the uniform topography across the Isobel Lake hillside, certain terrain attributes such as slope and LS-factor displayed significant relationships with fine-scale water distribution within the top 30 cm of soil, and shows a general agreement with what we already know about topographic influence on the variability of soil properties over much larger scales. This includes the dependence of water availability on slope gradient and elevation (Florinsky et al. 2002; Romano and Palladino 2002), and that many regression models using topography for explaining variability in soil available water typically do not produce R2 values greater than 0.50 (Florinsky et al. 2002). Although these relationships are important to quantify on finescales, due to the relatively low levels of predictive accuracy (R2 values < 0.50) of most models it becomes uncertain as to whether or not attempting this type of intensive sampling procedure over such a small area will be worth the time and effort. If, for example, scalping the forest floor prior to planting is proven to be a successful method for improving soil water availability for new Douglas-fir seedlings, then it can be suggested that the entire harvested area be scalped as it will likely be less effort than conducting an exhaustive soil sampling design similar to what was completed in this study. Future research involving the modeling of soil properties using terrain attributes could improve fine-scale sampling design by uncovering important relationships between individual and combined topographic grid layers with soil properties. Improving our understanding of these crucial relationships on progressively smaller scales, we can begin to design soil sampling methods that are cost and time effective while still providing adequate information about a site‟s overall soil variability. For example, the location of the soil pits within this study were distributed in a way that I thought would capture the apparent variability in soil properties throughout the entire site, but I failed to consider the topographic placement of each pit. When attempting to quantify the true relationship between soil variation and different topographic features it is important that local extremes in topography are sampled to provide the model with the best, and most diverse, information available at the site. Within other landscapes that exhibit more distinct topographic relief, R2 values within the range of 0.39 – 0.82 are common (Odeh et 89 al. 1994; Gessler et al. 1995; Li et al. 2017). To clarify, in a typical predictive soil mapping study, the soils located in valley bottoms and ridge peaks are very important to consider as these areas have the highest divergence in topography, and will likely display different soil characteristics as a result. In addition, sampling in an area that has considerably more topographic variability would likely result in stronger relationships between the multiple topographic grid layers and certain patterns in soil variability. If future studies can begin to uncover stronger relationships between topography and soil variation than what was discovered here, then forest management within these ecosystems could benefit as less soil measurements on the ground would be required to accurately describe the sites water availability. Combining growth data of seedlings planted throughout the site over multiple years with soil data could provide better information to managers working within dry, forested ecosystems. To conclude that fine-scale soil water availability within these sites is the absolute limiting factor to seedling growth and success, it would be necessary to quantify the variability in seedling growth across the study site and to map its distribution in relation to soil water availability. Because the 2017 growing season was unusually dry, it would also be very important to monitor seedling growth over many years following their establishment to evaluate the residual effects of low or high-water availability throughout a seedlings lifetime, and to be able to quantify the effects of drought conditions on longer term seedling survival. The effects of climate change on species distribution patterns, tree growth, and overall site productivity is an issue that has gained an immense amount of focus over the last several decades throughout Canada and the rest of the world (Mathys et al. 2014; Kirchen et al. 2017; Littke et al. 2018). The future sustainability of British Columbia‟s historically successful forest industry is uncertain, and forest managers are raising concern over the unpredictability of local climate adaptations and its influence on tree regeneration strategies. Depending on the geographic region and associated climate change situation, the results of increasing temperatures will vary, thereby making it crucial to study and understand the complex interrelationships 90 shared by forests and local climate variables. In areas similar to the Isobel Lake study site, where annual drought is common during the growing season, seedling regeneration and forest productivity is likely to be greatly affected by climate change, and could significantly alter the way we manage these dry ecosystems in the future. The relationships between Douglas-fir trees with climate have been strongly linked to water limitations (Littke et al. 2018), however, the variability in shallow soil available water supply within Douglas-fir forests is relatively uncertain, and therefore justifies the need to further study these relationships. Furthermore, recent research suggests that variation in soil properties will influence the distribution and growth of many different tree species and that current maps which represent coarse-resolution interpretations of soil parameters were insufficient and need to be refined (Coops and Waring 2001; Herbst et al. 2006; Mathys et al. 2014). Temperature and precipitation greatly affect soil available water supply within Douglas-fir forests, and future climate change models predict a decline in available soil water depending on the severity of temperature increases. Therefore, forest managers must gain a deep understanding of forest soil variability and its relationship with climate, topography and above all, water availability. The quantification and long-term monitoring of variability in soil available water on both large and small areas will greatly improve our understanding of how climate change may impact different forest types through time. The results of this study provide a detailed look at the level of spatial soil variability that occurs within small areas of the dry IDF, and how that variability relates to available water. Through the effective quantification and monitoring of available soil water throughout the past 3 growing seasons, in addition to exploring its relationship with certain topographic features, forest managers are provided with an introductory look into soil water dynamics within these sites and can begin to link this information with future climate change models. Additionally, this thesis highlights the usefulness, or lack thereof, of fine-scale interpolations of measured soil properties within forested landscapes, and demonstrates the need to sample local soils to gain accurate representations of available soil water at any given site. 91 Although the variability in physical soil properties was not particularly large at the Isobel Lake site, coarse fragments did have a large negative effect on the soil water holding capacity, and therefore provides merit for more precise understandings of forest coarse fragment variability when considering the overall influence climate change may have on tree seedling regeneration and growth. By providing a large fundamental dataset of physical and chemical soil properties in relation to topography, and how each property might influence plant available water supply, soil maps within these ecosystems can be refined and provide more information to help alleviate the uncertain productivity of hot, dry Douglas-fir forests into the future. LITERATURE CITED Coops, N.C., and Waring, R.H. 2001. Assessing forest growth across southwestern Oregon under a range of current and future global change scenarios using a process model, 3-PG. Global Change Biology. 7: 15-29. Florinsky, I.V., Eilers, R.G., Manning, G.R., and Fuller, L.G. 2002. Prediction of soil properties by digital terrain modelling. Environmental Modelling and Software. 17: 295-311. Gessler, P.E., Moore. I.D., McKenzie. N.J., and Ryan, P.J. 1995. Soil-landscape modelling and spatial prediction of soil attributes. International Journal of Geographical Information Systems. 9(4): 421-432. Herbst, M., Diekkrüger, B., and Vereecken, H. 2006. Geostatistical coregionalization of soil hydraulic properties in a micro-scale catchment using terrain attributes. Geoderma. 132: 206-221. Li, A., Tan, X., Wu, W., Liu, H., and Zhu, J. 2017. Predicting active-layer soil thickness using topographic variables at a small watershed scale. PLoS ONE 12(9): 1-17.e0183742. https://doi.org/10.1371/journal.pone.0183742 Littke, K.M., Zabowski, D., Turnblom, E., and Harrison, R.B. 2018. Estimating shallow soil available water supply for Douglas-fir forests of the Pacific Northwest: climate change impacts. Can. J. For. Res. 48: 421-430. 92 Mathys, A., Coops, N.C., and Waring, R.H. 2014. Soil water availability effects on the distribution of 20 tree species in western North America. Forest Ecology and Management. 313: 144-152. Odeh, I.O.A., McBratney, A.B., and Chittleborough, D.J. 1994. Spatial prediction of soil properties from landform attributes derived from a digital elevation model. Geoderma. 63: 197-214. Romano, N., and Palladino, M. 2002. Prediction of soil water retention using soil physical data and terrain attributes. Journal of Hydrology. 265: 56-75. TWI 9. Topo. Wetness Index LiDAR LiDAR LiDAR PO MCA 8. Modified Catchment Area LiDAR 12. Positive Openness CS 7. Catchment Slope LiDAR LiDAR CA 6. Catchment Area LiDAR VD TPI 5. Topo. Position Index LiDAR 11. Valley Depth LS 4. LS - Factor LiDAR LiDAR CI 3. Convergnece Index Source LiDAR LiDAR 10. Multiresolution of Valley Bottom Flatness MRVBF Abreviation Elev Slope Raster Grid Layer 1. Elevation 2. Slope Feature Reference Description of ground elevation at each cell Conrad (2012) Description of hill slope in radians Travis et al. (1975) Relief as a set of convergent areas (channels) and Koethe & Lehmeier (1996) divergent areas (ridges) Based on slope and specific catchment area (SCA, Boehner & Selige (2006) as substitute for slope length). Compares elevation of each cell to the mean elevation of a specified neighborhood around that Guisan et al. (1999) cell Top-down processing of cells for calculation of flow O'Callaghan & Mark (1984) accumulation and related parameters Average slope perpendicular to the nearest contour Boehner et al. (2002) line at each grid point Corresponding drainage area per unit contour width Boehner and Selige (2006) [m²·m-1] Representation of divergent and convergent flow Boehner and Selige (2006) pattern in hilly terrains Uses slope and elevation to classify valley bottoms Gallant and Dowling (2003) as flat, low areas Vertical distance to a channel network base level Conrad (2012) The dominance (positive) of a landscape position related to how wide a landscape can be viewed from Yokoyama et al. (2002) any position Table 5.1: Description and reference for the 12 selected topographic grid layers used as environmental covariates for soil spatial predictions. 93 5.0 APPENDIX I 94 Table 5.2: Comparison of multiple linear regression (MLR) and geographically weighted regression (GWR) models using the 7 created principle component analysis (PCA) grids as predictor variables for various soil properties at all 4 2 measured depths. Average R values are created by running 10 tests with training datasets using leave-one-out cross validation. N/A values means no significant predictors were identified at the p = 0.10 level. SOM = soil organic matter, CF = coarse fragment concentration, BD = fine fraction bulk density, and SWHC = soil water holding capacity. GWR (7 PCA Grids) Depth (cm) 0-15 15-30 30-60 60-100 Variable Sand Clay SOM CF BD SWHC Sand Clay SOM CF BD SWHC Sand Clay SOM CF BD SWHC Sand Clay SOM CF BD SWHC Average R2 0.410 0.330 0.156 0.431 0.459 0.327 0.135 0.162 0.204 0.313 0.440 0.306 0.166 0.098 0.204 0.462 0.409 0.414 0.251 0.272 0.244 0.288 0.684 0.381 SD 0.010 0.010 0.009 0.010 0.012 0.010 0.008 0.014 0.007 0.008 0.006 0.009 0.006 0.005 0.002 0.013 0.012 0.009 0.004 0.004 0.016 0.004 0.004 0.005 MLR (PCA Grids) Average Adj. R2 0.300 0.234 0.038 0.276 0.277 0.082 N/A N/A 0.146 0.146 0.379 0.193 N/A N/A 0.049 0.381 0.335 0.300 N/A N/A 0.089 0.159 0.604 0.28 SD 0.009 0.003 0.015 0.048 0.063 0.014 0.008 0.012 0.006 0.011 0.002 0.013 0.018 0.009 0.035 0.005 0.003 0.012