by Michael More BSc, University of British Columbia, 2010 A thesis submitted in partial fulfillment of The requirements for the degree of Master of Science in The Faculty of Graduate Studies (Geography) The University of British Columbia (Vancouver) October 2012 © Michael More, 2012 Building a Global Sediment Database ii Abstract Global patterns of soil loss are poorly understood. While recent research has explored the physical processes that drive soil erosion, through the study of sediment transport, few spatial analyses of these processes have been conducted. One major reason why the spatial element of sediment dynamics has not been studied is a lack of data: datasets are often expensive and difficult to obtain. To meet this need, this thesis proposes the construction of a Global Sediment Database, which will be freely available to all. The Database will be updatable, and will contain a detailed Data Quality Report so researchers can determine the most effective use of the data. The Data Quality Report will also quantitatively summarize the error and uncertainty of each dataset. This thesis will also demonstrate how the Database can be used to conduct spatial analyses of sediment processes using Geographic Information Systems (GIS). To this end, various spatial interpolation methods will be explored and evaluated, using the Yellow River as an example. iii Table of Contents Abstract ....................................................................................................................................... ii Table of Contents ....................................................................................................................... iii List of Tables .............................................................................................................................. iv List of Figures ............................................................................................................................. v Acknowledgements .................................................................................................................... vi Dedication ................................................................................................................................. vii Chapter 1 Introduction .................................................................................................................... 1 1. Motivation ............................................................................................................................... 1 1.2.1 Sediment Transport......................................................................................................... 5 1.2.2 Sediment Yield ............................................................................................................... 7 1.2.3 Controls of Sediment Processes ..................................................................................... 8 1.3 Research Questions ............................................................................................................. 16 Chapter 2 Methods ........................................................................................................................ 18 2.1 Database Structure and Organization .................................................................................. 18 2.2 Data ..................................................................................................................................... 20 2.2.1 Data Gathering .............................................................................................................. 20 2.2.2 Data Quality .................................................................................................................. 23 2.3 Using the Global Sediment Database to Map Sediment Yield ........................................... 31 2.4 Using GIS to Create Interpolated Sediment Maps .............................................................. 34 Chapter 3 Analysis & Discussion ................................................................................................. 44 3.1 Background Information on the Yellow River .................................................................... 44 3.2 Rationale for Selecting Kriging Algorithms ....................................................................... 48 3.2.1 Range, Nugget and Sill ................................................................................................. 52 3.2.2 Considering Error Values ............................................................................................. 54 3.2.3 Semi-Variography ........................................................................................................ 57 3.2.4 Comparing Universal and Ordinary Semi-Variograms ................................................ 62 3.3 Analyzing Sediment Maps and Error Surfaces ................................................................... 66 3.4 Comparing Kriging Results of Other Rivers ....................................................................... 74 Chapter 4 Conclusions .................................................................................................................. 81 References ..................................................................................................................................... 83 iv List of Tables Table 2.1 Example entries from the Global Sediment Database (showing the Yellow River)..... 19 Table 2.2 Further detail on example database (Yellow River; same entries as table 2.1)............ 20 Table 2.3 Example of published tabular data (modified from Lu & Higgit, 1999)...................... 21 Table 2.4 Evaluating spatial precision.......................................................................................... 25 Table 2.5 Sample data stations from the Yangtze River: temporal resolution analysis................ 27 Table 2.6 Evaluation of temporal precision using years of data available.................................... 29 Table 2.7 Evaluating point density............................................................................................... 30 Table 2.8 Data Quality Report for selected datasets.................................................................... 31 Table 2.9 Factors of interpolation................................................................................................ 35 Table 3.1 Sediment relation details for Yellow River sub-basins................................................ 46 Table 3.2 Comparing error values of ordinary and universal kriging.......................................... 56 Table 3.3 Yellow River ordinary kriging error values................................................................. 59 Table 3.4 Yellow River universal kriging error values................................................................ 60 Table 3.5 Yellow River universal kriging error ranking.............................................................. 61 Table 3.6 Yangtze River ordinary kriging error values................................................................ 75 Table 3.7 Yangtze River universal kriging error values............................................................... 76 Table 3.8 Danube River ordinary kriging error values................................................................. 77 Table 3.9 Danube River universal kriging error values................................................................ 77 v List of Figures Figure 1.1 The components of sediment yield (modified from Church, 2006) .............................. 6 Figure 2.1 Map of sediment yield in Romania (Radoane & Radoane, 2005) ............................... 22 Figure 2.2 Data stations in Bolivia (Aalto et al., 2006) ................................................................ 23 Figure 2.3 Example of a regional sediment relation (upper reach of the Yellow River). ............. 33 Figure 2.4 Graph showing simple and ordinary kriging (modified from ESRI, 2011 .................. 38 Figure 2.5 Graph showing statistical model of universal kriging (modified from ESRI, 2011). . 40 Figure 2.6 Map of the Yellow River, showing the rectangular kriging surface. .......................... 41 Figure 2.7 Error surface for the Yellow River .............................................................................. 41 Figure 2.8 Yellow River interpolation surface, trimmed to watershed boundary ........................ 42 Figure 3.1 Yellow River reference map ........................................................................................ 44 Figure 3.2 Yellow River sediment relations ................................................................................. 47 Figure 3.3 Flow chart—preparing data for spatial analysis. ......................................................... 50 Figure 3.4 Range, nugget, sill, and partial sill in kriging (modified from ESRI, 2011). .............. 52 Figure 3.5 Example kriging algorithm of the Yellow River, with statistical details (see text). ... 53 Figure 3.6 Kriging error values—predicted values (y-axis) vs. measured values (x-axis)........... 55 Figure 3.7 Kriging semi-variogram options (modified from ESRI, 2011). .................................. 58 Figure 3.8 Ordinary circular kriging semi-variogram................................................................... 63 Figure 3.9 Ordinary Gaussian kriging semi-variogram ................................................................ 64 Figure 3.10 Universal circular kriging semi-variogram ............................................................... 65 Figure 3.11 Universal Gaussian kriging semi-variogram ............................................................. 66 Figure 3.12 Sediment map of the Yellow River: ordinary Gaussian kriging ............................... 68 Figure 3.13 Error surface map of the Yellow River: ordinary Gaussian kriging ......................... 69 Figure 3.14 Sediment map of the Yellow River: universal Gaussian kriging .............................. 72 Figure 3.15 Error surface map of the Yellow River: universal Gaussian kriging ........................ 72 Figure 3.16 Flow chart detailing selection of best kriging semi-variogram ................................. 78 vi Acknowledgements My supervisors, Marwan Hassan and Sally Hermansen, were tireless, patient, and supportive. The opportunity to work with both of them was a major factor in deciding to undertake a Master’s degree. For GIS matters, Brian Klinkenberg was always generous with his time. Jose Aparicio and Alejandro Cervantes-Larios were equally noble in humouring my pestering questions. The long hours of coding data in excel was supported by many students from GEOB 405/508 in 2010, including Alexis Carter, Ken Hamilton, Michael Currie, Ryan Sitzo, and Dave Hamilton. Data for the USA was generously provided by Andrew Simon and Lauren Klimetz. Data for China came from Marwan’s own research, while the Canadian data was donated by Michael Church. Olav Slaymaker provided a sympathetic ear as I searched for Russian data, while Maria Radoãne provided data from Romania to supplement her published figures. Joris de Vente kindly provided his PhD thesis, which contained data for Spain and Italy. Scott Graham was inspirational, especially in thinking about how the inner workings of GIS can be explored. Support from family and friends was crucial to my success. Both my parents were teachers, and their example inspired my own academic career. My mom shared tips from her days of thesis-writing, while my dad had good strategies for how to keep track of mountains of data. Spousal support is always essential in these endeavours. The patience, flexibility, and perspective from loved ones are not always recognized, but they are certainly appreciated whole- heartedly. vii Dedication for Buck, who kept me sane for mom, who taught me patience for da, who always knew where to find the sharpest knife in the dr 1 Chapter 1 Introduction 1.1 Motivation Soil is an essential natural resource (Walling, 2005; Montgomery, 2007): thus, the continuing depletion of the world’s soil is a major threat to humanity. Yet soil is arguably the most undervalued element to our survival: oil, gas, and minerals get much of the attention and resources (Cooke & Doornkamp, 1990; Montgomery, 2007). Although soil loss has been a recognized problem since at least the late 1800’s (examples cited in Cooke & Doornkamp, 1990), the issue remains problematic even now, well into the twenty-first century. Rainfall, runoff and wind are responsible for much of the world’s soil loss from year to year (e.g., Jansson, 1981; Church et al., 1999; Walling & Fang, 2003; Church, 2006). However, focusing on these processes is too narrow a perspective, as human activities greatly increase a landscape’s vulnerability to erosion. For example, deforestation is contributing to soil erosion, as vegetation acts to keep soil in place. Another example is intensive agriculture. As a result of logging, soil is no longer bound up in the root structures of the trees and shrubs, logging roads have likely been cut into the landscape, and the surface itself has been disturbed through the action of felling and removing trees. All these factors increase the likelihood of subsequent severe soil erosion (e.g., Cooke & Doornkamp, 1990; Montgomery, 2007). Agriculture is also a major factor in soil erosion. The westward American push across North America throughout the mid to late 1800’s, was spurred, to a significant degree, by soil loss (summarized in Montgomery, 2007). As valleys and slopes were cleared for farming, little attention was paid to soil retention, so overland runoff could often cause a decrease in arable land, as once fertile soil was washed downstream. Rather than address the issues, farmers would 2 move farther west, and repeat their mistake (Montgomery, 2007). Likewise, cash-cropping practices in the Amazon jungle jeopardize the forest ecosystem, as more and more areas are cleared, exacerbating soil erosion (Montgomery, 2007). Though soil-erosion processes are relatively well-understood there is still not a collective, global consensus regarding appropriate conservation strategies (Cooke & Doornkamp, 1990). The task of adequately measuring soil loss can be difficult, expensive, and time-consuming, and few countries possess the resources do so effectively over the long term. Despite calls from conservationists as early as the late 19 th century (i.e. George Marsh, 1864, cited in Cooke & Doornkamp, 1990), soil erosion continues to be a concern in many parts of the world. Most responses to this problem have taken one of the following four approaches (Montgomery, 2007). First, scientific research has been conducted to understand and predict soil-erosion dynamics. Second, engineers have developed technology to help mitigate or prevent soil loss. Third, attempts have been made to apply research results to land management policies. Finally, reviews of all these actions have been conducted, in order to assess the usefulness of various potential solutions (Montgomery, 2007). These efforts can be strengthened by understanding the spatial patterns of sediment yield in the landscape. Understanding how soil loss occurs is only part of the question: knowing where and how much sediment is moved is important, so researchers can maximize the effectiveness of their soil conservation strategies. The fact that few studies having systematically investigated the spatial concerns of sediment dynamics and soil loss underscores a more pressing concern: data is often expensive and difficult to acquire, and thus the availability of such information is low. 3 To reach a global awareness of where soil loss occurs, large spatially representative datasets are needed. While datasets exist for some parts of the world, there are no standard methods for collecting, processing, and documenting such data. While there has been a greater emphasis on the consideration of spatial scale in sediment yield analyses (especially Church et al, 1999; and Hassan et al, 2008), a detailed investigation into which specific tools might be best for spatial sediment analyses has yet to be published. Both Church et al. (1999) and Hassan et al. (2008) use interpolation to conduct their spatial analyses of sediment patterns, but there are dozens of ways of interpolating spatial data. Which is most suitable? The thesis proposes to fill both these gaps. First, we will build a Global Sediment Database, which will incorporate previously published datasets into an organized, easy-to-use format. Second, we will examine and evaluate methods of spatially analyzing sediment data. Both these goals will incorporate data quality analysis. For the Database, each data source will be quantitatively assessed enabling potential users to see at-a-glance how reliable a dataset might be. In investigating methods of spatial analysis, an example dataset—the Yellow River, in China—will be processed using Geographic Information Systems (GIS). To evaluate the appropriateness of the interpolation method, statistical measures of error values will be analyzed. This process of evaluating data quality, and systematically selecting the most appropriate analytical interpolating tools, can serve as a template for other researchers, a kind of step-by-step guide for effectively using the data that will be contained in the Global Sediment Database. First, however, we need to take a step back and explore the data that will be used to create the Global Sediment Database. Sediment yield will be the primary component of the Database. A three-part literature review will examine background information that is necessary to understand why spatial analysis of sediment yield is important, and how such research can be 4 conducted with the available data. In the first section we will introduce basic concepts of sediment transport and yield and why these values are important; second, we will review the main controls on the spatial and temporal patterns of global sediment yield; third, we will demonstrate that considering the spatial pattern of sediment processes is essential in understanding sediment dynamics in the landscape. We will also demonstrate that data quality is important to evaluate when investigating sediment yield, as data quality is a major hurdle in creating the Database, and has prevented progress on this topic. Investigating these topics provides the motivation for this project: to build a Global Sediment Database, and to investigate methods for spatial sediment yield analysis and extrapolation. 1.2 Literature review Few long-term data records for measuring and quantifying sediment processes exist for most parts of the world, since it is difficult to obtain such data. For example, installing and operating monitoring stations is expensive in terms of equipment, transportation to remote locations, and labour. For manually sampled data, there is a risk of drowning, while for automatically gathered data, instruments can be damaged or displaced by severe events. Difficulties can also arise from the data samplers themselves. The sample and recording of sediment (both bedload and suspended) is limited to the sample capacity or to the quality of the turbidity sensors. Accounting for the danger and difficulty of sampling extreme events is expensive: equipment must be reinforced, repaired or replaced; field researchers must be sent in groups, and with proper training; and extra sampling days should be planned, to account for lost data. 5 Once logistics have been arranged there is the issue of landscape complexity. Deciding where to place stations is not straightforward, nor is it possible to guarantee complete representativeness. A consequent problem is that the sensitivity of basin sediment dynamics to environmental change is not fully known, despite the comparatively extensive datasets regarding other hydrological and meteorological parameters, such as runoff, river floods, and precipitation (e.g., Church et al., 1989; Walling, 2005). Even though these fundamental processes are well- understood, this knowledge does not fully alleviate the comparatively unmonitored sediment processes of a basin. Further, the reliability of the few long-term sediment databases is open to question, given the persistent problem of obtaining accurate estimates of sediment loads (Walling & Webb, 1981, 1987, 1988; Leeks, 2005). Nevertheless, progress has been made. Well-established, well-documented data sources exist in the USA and in China, organized by national governments, with ongoing monitoring at hundreds of sites. Typically, such data stations report sediment concentration, sediment load and sediment yield (excluding dissolved material and bedload). To represent temporal patterns in sediment dynamics, there is a need for long term datasets. By using long-term datasets, we can create standards for data quality analysis and metadata documentation Thus, the Global Sediment Database will become a model for data collection and processing strategies. 1.2.1 Sediment Transport Rivers and the watersheds they drain are the fundamental components of the landscape. Fluvial processes dominate land denudation, and therefore control the shape of the earth’s surface. The physical processes that control the erosion, transport, and deposition of material from catchment surfaces and within river channels have been well-documented (e.g., Walling, 2005; Church, 2006). Sediment of different calibres moves in different ways through a stream 6 and forms distinctive deposits which influence channel form. Fine material (wash load), generally less than about 0.125mm, moves in suspension (figure 1.1) and is deposited overbank in floodplains, or in backwaters, to form vertically accreted sediment bodies. Medium to coarse sand moves intermittently in suspension (bed material load) or in traction and forms vertically and laterally accreted bedforms within the channel. Gravel moves in traction (bed material load, or bedload) and is deposited in laterally accreted bars within the channel. This classification of fluvial sediments differs from the customary one based upon mode of transport and measurement (e.g., Church, 2006), but is more meaningful from both morphological and sedimentological viewpoints. Figure 1.1 The components of sediment yield (modified from Church, 2006). Walling (2005) claimed that measuring all components of sediment yield is impractical. Even in non-flood conditions, it is often difficult to separate bedload, saltation, and suspended load, as these processes often intermingle. Measuring bedload is also problematic (see Church, 2006). Pit traps are perhaps the most representative way to sample bedload, but collecting the data (potentially thousands of tons of sediment) is expensive and difficult, and they are too small relative to the amount of sediment that passes through a system, resulting in full traps, which 7 renders their data useless (Leeks, 2005). With these two components being difficult to measure reliably, researchers tend to turn to the suspended load, which can be measured continuously throughout the year using manual or automatic samplers (Walling, 2005). Despite the lack of bedload data, measurements of only suspended sediment yield can still be useful, because most sediment moves in suspension. Several studies have attempted to quantify what percentage of total sediment yield is bed load vs. suspended load. While these estimations have varied, they have typically been in the range of ~10% for bedload (Walling & Fang, 2003; Leeks, 2005; Church, 2006). This is one reason why suspended sediment data are often used in research—the other reason being the logistical challenges of measuring bedload, as mentioned above (Church et al., 1989; Nearing et al., 2005; Walling, 2005). More detail on measuring suspended sediment yield can be found in Milliman & Meade, 1983; Jansson, 1988; Dedkov & Mozzherin, 1992; Walling & Fang, 2003; and de Vente et al., 2007). 1.2.2 Sediment Yield Sediment yield (SY) represents the suspended load of a river, but is typically used to represent the integrated yield from the relevant contributing watershed. However, the sediment that is eroded will not all necessarily be completely transported through to the terminus of a river network: some sediment will remain stored in the watershed, albeit in a different location. Again, the term Sediment Delivery Ratio (which compares the material that is delivered to a site with the material that is evacuated from it) is used to describe the fraction of material that reaches a particular point in a river system. Typically, less sediment leaves a system than is eroded: deposition within a watershed is thus usually greater in large basins, due to the increased storage potential. 8 The relative quantities of sediment erosion, transport, and deposition are not always well- known. Often it is assumed that there is a correlation between yield and the size of a basin (Walling, 1977; Church et al., 1989; de Vente et al., 2007), whereby larger basins present greater opportunities for storage, thus limiting delivery to the ocean. The next section will more thoroughly explore the significant controls on sediment dynamics. 1.2.3 Controls of Sediment Processes At the continent or landscape scale, climate is one of the main dominant controls of sediment dynamics, in addition to topography and human activity. Two examples are presented in the studies by Langbein and Schumm (1958) and Wilson (1973) of sediment yield in the U.S. Langbeing & Schumm (1956) related SY to mean annual precipitation (P), and showed that it reaches a peak at about 300 mm (semiarid climate). Their model states that SY declines for P < 300mm because there is too little runoff to erode much material, and declines for P > 300mm because of the impact of vegetation. Wilson (1973) showed that in the U.S. a peak in SY occurs at about P ~ 380 mm, and that there is a second peak at P 1,500-2,000 mm in seasonally wet climates (Mediterranean and tropical wet/dry climates). This peak is higher than that at P = 380 mm because vegetation becomes attenuated during the dry season, and virtually all the rain occurs in just a few months of the year. These models, however, are too simplistic, because they do not include other factors than control sediment yield. Langbein and Schumm (1958) in particular did not consider the role of landuse, geology, or landscape history on sediment yield. The Langbein & Schumm model was explored by Milliman & Meade (1983), who explored discharge as a key factor in explaining sediment yield magnitude. Walling (1977) similarly noted peaks in suspended sediment transport shortly after heavy precipitation, and that annual precipitation rates are also significant. Using a similar model to Langbein & Schumm 9 (1958), Walling (2005) argued that moderate rainfall leads to the highest sediment yields, between 300 and 500mm per year. However, not all data support this model; monsoon regions, for instance, have rainfall well in excess of 500mm per year, and have extremely high annual suspended yield. Jansson (1988) concluded that global trends of sediment yield could be categorized based on broader climatic variables, such as overall temperature, dryness, and continentality. In another case, Yair & Enzel (1987, cited in Goude, 1995) examined streams in the Negev desert in Israel, which is in an arid/semi-arid climate. Under these conditions, the authors concluded that there was a simple linear relationship between runoff and SY, and also that topography played an equally important role. However, the reliability of this study is questionable, since it was based on only two data points. The effects of topography on sediment yield are evident in the high SY values for mountainous parts of the world. Steeper terrain leads to higher rates of erosion in these areas, whereas flatter topography like plains has lower erosion, and a greater likelihood of deposition. It is the landscape topography that determines how a river network flows: concavities tend to trap sediment longer than sloped surfaces. Indeed, different topographic features of the landscape vary in their ability to trap or store sediment. In a similar way, different locations along a river system will store different amounts of sediment (Church, 2006). Flatter areas, like floodplains and alluvial fans, tend to store material, while banks and channel slopes tend to be more effective transport zones. Typically, material stored within the channel itself remain stored for less time than material in the floodplain or alluvial fans, as annual floods are effective in mobilizing such in-stream deposits, especially ones comprised mainly of medium or fine sediment. This interaction between 10 hydrology and topography is circular: a river channel tends to exist in the landscape because a slope is sufficiently steep for water to flow continuously. Flowing water helps dig a deeper channel, which then feeds eroded sediment into the river system. As the river shapes the surrounding topography, so the slopes, valleys, and ridges of the watershed guide the river’s path. It is the topography, hydrology, and sediment supply, which control sediment storage locations (Leeks, 2005); higher supply and flatter topography generally lead to more storage. In terms of erosion and transport, higher slopes are associated with both, and local topography will in particular determine the path transported material will take as it travels; for instance, from the hillslope towards the stream network (Church, 2006; Slaymaker, 2006). Landscape history is important to consider in studies of sediment yield. At the broadest scale, as demonstrated by Church et al. (1989) and Church & Slaymaker (1989), regions such as British Columbia, Canada, that are still influenced by their glacial history behave very differently from non-glaciated basins. The authors investigated the spatial pattern of specific fluvial clastic sediment yield (SSY) in BC, and found a positive correlation overall with drainage basin area (Ad) for basins up to 30,000 km 2 . When plotted, glaciated basins were close to the upper envelope of their data. These findings are in contrast with previous research of predominantly agricultural basins, where declining sediment has been reported in response to an increase in area, as larger basins provide more opportunities for sediment storage. In BC, however, glacial sediment deposits were left along valleys; these deposits are a secondary sediment source for the current fluvial regime (Church et al., 1989; Church & Slaymaker, 1989). The evacuation of these sediments from valley bottoms is a significant factor in the positive relationship found by Church et al. (1989) for basins smaller than 30,000 km 2 . On the other hand, the decline in sediment yield for larger basins was due to barriers that prevented 11 the sediment from entering the river system. This implies that BC’s landscape has not recovered from the legacy of glaciation (Church & Slaymaker, 1989). While topography and landscape history may determine pattern of sediment yield (or channel dynamics), geology affects the rate at which material may be eroded. SY is positively correlated with the erodibility of the surficial rocks. The Yellow River (Huanghe) in China is a striking example. The highly erodible loess material in its middle reach contributes to the high amount of sediment yield produced by the river (Wang et al., 2006; Hassan et al., 2008). The erodibility of surficial material is thus an important geologic control on sediment yield (Church et al., 1989; Slaymaker & Church, 1989; De Vente et al., 2007). In small basins, local variations in surficial material can produce a wide range of values (Church et al., 1989), which again suggests caution in examining basins of differing size (Milliman & Meade, 1983; Church et al., 1999; Walling & Fang, 2003). The geology of larger basins also needs to be considered. Smaller basins are likely to be more homogeneous in terms of lithology. Larger basins are more likely to be complex in terms of climate and geology. Larger basins are thus more likely to exhibit a complex response in terms of sediment yield. De Vente et al. (2007) presented summary data for hundreds of rivers from around the world. They defined geologic categories for these rivers: basins with harder rocks (e.g., classified as “Crystalline Mountain”, “Crystalline Volcanic”, etc.) had lower sediment yields than other basins in similar regions. Regions in other geologic classifications—those with moderately hard or softer rocks—did not show a clear trend. De Vente et al. (2007) cite Romanian data from Rãdoane & Rãdoane (2005) as another example: in examining the Romanian basins, the authors similarly found lower sediment yield values for basins with volcanic or crystalline surficial materials. Geologic faults also play a role in sediment 12 yield, as the presence of a fault can affect the entrenchment or confinement of a river channel. This was a factor in the analysis of British Columbia data in Church et al. (1989). Another factor that influences sediment yield values is vegetation, which is related to climate. Cooke & Doornkamp (1990) asserted that vegetation is a significant deterrent to soil erosion and thus sediment yield. Walling & Fang (2003) similarly found strong correlations between the type of vegetation cover and the magnitude of sediment transport, with forested areas often reporting the lowest sediment yield. Vegetation prevents soil erosion in four ways (Cooke & Doornkamp, 1990): it interrupts rainfall and stemflow; plants consume water, which further reduces waters flow; roots help keep soil in place, and conduct water away from soil surfaces; and the increased surface roughness provided by vegetation reduces runoff velocity and the resultant erosive potential. All these processes, however, are dependent on climate, and changes in climate will similarly result in vegetative changes. Vegetation type, however, is not the only relevant factor that controls how plants affect sediment transport. Cooke & Doornkamp (1990) noted that canopy density, plant height, degree of cover, root density, water consumption, and level of mulch all need to be considered. Further, as noted by Goudie (1995), the effect of vegetation changes with the time of year, with climate, and with fire history. Fires are particularly important, and in some regions, such as northern California (Cook & Doornkamp, 1990), sharp increases in sediment yield occurred soon after severe fires: not only did the removal of vegetation leave surfaces vulnerable, but the charred, hydrophobic organic matter that remained actually increased the likelihood of erosion. Severe weather events have a short-term, albeit important, role in sediment dynamics. Milliman & Meade (1983) reported that especially in dry, continental basins, upwards of 80% of the annual sediment yield can occur from severe events alone (in some cases, from only a single 13 event). Meade et al. (1990) likewise analyzed the effects extreme events had on two rivers in the USA. For the Eel River, in California, almost 170 million tons of sediment was transported in three days; under normal conditions, it would have taken 10 years to move that much material. For the Juniata River in Pennsylvania, Hurricane Agnes moved three or four years’ worth of sediment. Yet despite the potentially extreme rates of sediment transport that can occur as a result of extreme weather, such values are often omitted from broad-scale sediment dynamics studies, such as those conducted by Jansson, 1988; Church et al., 1989 and 1999; Walling & Fang, 2003; and Walling, 2005. In each of these studies, data points that were strongly affected by severe storms were either omitted or given “less weight” (Church et al., 1989). More frequently, values for storm days may be missing, or unreliable, because it can be dangerous, difficult, and expensive to gather such data (Walling, 1977). It has been suggested (Church et al., 1989) that considering consecutive years may help reduce the potential influence of extreme events. Church et al. (1989) recommended analyzing the temporal variability of a dataset, so that when data records are too short, they are compensated for, or simply omitted. This thesis adopts a similar perspective. In addition to the natural forces that affect sediment yield described above, anthropogenic forces are also significant. The Yellow River is again a good example, as the changes in agricultural practices there have had a dramatic effect on sediment yield (Wang et al., 2006). Hassan et al. (2008) further note that changes in soil conversation practices and dam construction in the Yellow River watershed both resulted in significantly reduced sediment yields over time. The intensive agriculture in the small basins of the Loess plateau had led to some of the highest annual sediment yield values in the world (Walling, 2005), although values have declined in the 14 past few decades (Wang et al., 2006). Walling & Fang (2003) similarly argued that agricultural practices remove vegetation and disturb the soil, both of which make the surface more susceptible to runoff erosion. There has been much study regarding the issue of soil loss and agriculture. Indeed, as noted below, traditional sediment yield trends are not always applicable since the vast majority of sediment yield data has been focused on heavily cultivated basins (Church et al., 1989; Leeks, 2005; de Vente et al., 2007). Finally, Walling (2005) estimated that a majority of the world’s sediment yield is initiated by human activity. Dam construction represents another key anthropogenic factor in sediment dynamics. Walling and Fang (2003), in their global study of sediment yield, ranked dams as the most important anthropogenic factor for sediment dynamics. Their qualitative assessment of 145 rivers demonstrated that dams and reservoirs disrupt about twice as much sediment as landuse changes and climate change, combined. Another global study by Milliman & Meade (1983) presented a historical analysis of anthropogenic effects on sediment dynamics. Up until the Industrial Age, these authors cited the glacial legacy as a prime factor in determining the fluvial regime of much of the northern hemisphere, which in turn dictated sediment processes. Starting in the mid-twentieth century, however, Milliman & Meade (1983) argued that dam and reservoir construction became a more important factor, because of sediment starvation downstream of dams, and the vast accumulation of material behind dams. Walling (2005) similarly reported that lakes, reservoirs, and dams collectively are responsible for up to 80% of the anthropogenic influence on sediment dynamics. Meade & Moody (2009), in their long-term study of the Missouri-Mississippi River system, found that dams trap roughly 100-150 million tons per year, and that the overall sediment transport of the system has declined from about 400 million tons per year before 1900 15 to only about 145 million tons per year in the last few decades (1987-2006). However, the authors assert that only half of this decline was caused by dam construction. Part of the decline in sediment yield resulted from a fundamental change in the watershed, from being a transport- limited system to a supply-limited system. This demonstrates that human influences have the potential not only to alter the sediment regime, but to alter the fundamental way a river system operates. An extensive study of anthropogenic changes to sediment is presented by Chu et al. (1999). These authors examined nine major rivers in China, and found that over about 50 years, the total sediment flux from these rivers decreased by ~5.0 x 10 10 tons. They attributed about half of this decrease to dams, and the rest to soil and water conservation, increased water consumption, and in-channel sand mining. The conclusions of the research summarized above suggest that human activities can greatly influence sediment yield (e.g., Milliman & Meade, 1983; Chu et al., 1999; Walling, 2005). How substantial are such changes compared to average global levels of sediment yield? An excellent summary of global values for annual sediment yield is found in Leeks (2005). These values range from 8.3 to 51.1 billion tons per year, with an average of about 20 billion tons per year. Similar values are reported in Walling & Fang (2003) and de Vente et al. (2007). An account of the sediment yield of Europe is reported by Owens & Batalla (2003): about 2 billion tons per year are eroded from Europe, about 700 million tons of which reach the ocean. Millliman & Farnsworth (2011) made a similar estimate – about 2.5 billion tons for Europe, and about 27 billion for the world. With so much material being removed from the landscape, rivers are indeed a highly significant factor in determining our landscapes’ forms. Further, the variability of these values highlights the importance of creating the Global Sediment Database. 16 1.3 Research Questions There were two goals of the above literature review: to document the factors that control sediment processes and to demonstrate the role of spatial scales in sediment yield analysis. It is clear from examples such as Church et al. (1989) and Hassan et al. (2008) that considering the spatial scale is important to understanding watershed-scale trends in sediment yield. Understanding sediment processes is important because of its role in landscape evolution, soil erosion, and nutrient fluxes, and because human activities are threatening to greatly disrupt these processes on a global scale (e.g., Chu et al., 1999; Rãdoane & Rãdoane, 2005; Wang et al., 2006). It is because these global sediment processes are not well-understood that the Global Sediment Database needs to be created. Regarding data quality, issues of temporal resolution and the spatial distribution have been discussed. A related topic to data quality is metadata – or ‘data about data’. Namely, how data were collected, by whom, when, and under what standards. Even in looking just at the few examples provided in the preceding section, it is clear that there is not sufficient consistency from study to study to attempt to truly understand the complexities of global sediment dynamics. To address this problem we propose to compile a Global Sediment Database of previously published data sources, and explicitly document the data quality and metadata of all inputs. Specifically, data quality will be numerically scored for each entry in the database. To aid in a continued understanding of world sediment processes, this database will be made available to the public. Ideally, this database will serve as a foundation to which more data can be added over time, so that information from all over the world can be easily and reliably compared. Furhter, different interpolative GIS (Geographic Information Systems) techniques will be employed in creating these sediment maps. This will reinforce the idea that data quality and 17 metadata documentation matter, as not all methods are appropriate for all data. In addition to Yellow River sediment maps, error maps will be produced, and quantitative statistical measures for each interpolation method will be employed to evaluate which is more suitable for the data. By building a Global Sediment Database, quantitatively evaluating data quality, documenting methodological choices, and investigating interpolation and visualization methods, this thesis will explore the possibilities in examining the global sediment system. The database will be a tool for future researchers to employ, as we try to understand how and why sediment dynamics are changing around the world. 18 Chapter 2 Methods The global sediment database compiles existing datasets. Currently, over 7,000 data points have been processed and entered into the database. The plan is to continue to update these datasets and make them available to researchers. Since this is the first database of its kind, no template or standard methodology exists to guide data processing and data quality evaluation of the datasets. So at this stage it is important to establish a clear structure for the database; as it expands over time, newly added data can be easily inserted into the database. The first part of this section will describe the template of the global sediment database, and will also present a method of evaluating and recording the quality of each dataset stored in the database. The second half of this section documents how the datasets were used in calculating sediment yield, and how those values used to interpolate sediment maps using Geographic Information Systems (GIS). 2.1 Database Structure and Organization An initial problem with building the Database is what kind of information to include. Some datasets contain exhaustive information, with decades of data for hundreds of points. Others only document a handful of measurements, with only average values, taken from just a few years of sampling. For consistency in comparing and contrasting the rivers it would be helpful for the Global Sediment Database to report the same information for each measuring station, and the same variables from river to river. By enabling Database users to compare information from different rivers, a greater understanding of fluvial process should be fostered. To help narrow the scope of the Database, this thesis focuses primarily on suspended sediment data, with some additional data from dams. Bedload data are difficult to obtain, and 19 rarely reported in the literature; therefore, the focus for this thesis is on suspended load. Most of the suspended sediment data used in this thesis to build the Global Sediment Database was obtained from published sources, or downloaded from governmental agencies. Table 2.1 Example entries from the Global Sediment Database (showing the Yellow River). Station River or tributary Longitude (DD) Latitude (DD) Area (km 2 ) Elevation (masl) SSY (T/km 2 /yr) Batan Bagou R. 100.5500 35.2500 3483 2557.956039 732.571501 Qushi'an Qushi'an R. 100.1000 35.3333 5721 1045.090216 181.250693 Laqu Mangla R. 100.7500 35.5833 1717 773.47944 446.132464 Tongren Longwu R. 102.0167 35.5167 2832 610.047472 213.268881 Longwuhekou Longwu R. 102.1000 35.8333 4959 1554.023216 311.646142 Shuangcheng Daxia R. 103.0500 35.4833 6144 1124.171848 181.408992 Fengjiatai Daxia R. 103.2667 35.6500 6851 4324.315 630.052547 Tangnailang Daxia R. 102.8500 35.1500 1509 660.160398 433.529091 Wuchengyi Zhuanglang R. 103.1833 36.8667 2001 352.984984 172.478753 Zhoujiacun Zhuanglang R. 103.4333 36.2167 3982 3006.676208 752.191413 To help maintain data consistency, the Global Sediment Database will keep all points in the same format, so that users can easily search for the data they want. The core of the Database will be the specific sediment yield (SSY) for all available years, station area, and station coordinates. In addition, each entry in the Database will contain metadata, including data collection agency, country, river/tributary name, and station name (see table 2.1, above). Quantitative metadata will also be included, specifying the temporal attributes of each point: length of record, number of years recorded, and longest continuous period of data (see table 2.2). 20 By providing extensive metadata in the Database, users can assess the quality of individual data points. However, the Database on its own does not assess the quality of entire datasets. For example, how accurate are the coordinates for the Yellow River? How many points are there in total? And is the data from this source reliable overall? To answer these questions a Data Quality Report will be created which will quantitatively evaluate the data quality of each dataset. To understand how data quality was measured, however, it is necessary to explore how data were acquired and processed. Table 2.2 Further detail on example database (Yellow River; same entries as table 2.1). Country Watershed Station No. years Length of record Longest continuous data China Yellow River Batan 20 28 16 China Yellow River Qushi'an 18 21 6 China Yellow River Laqu 12 12 12 China Yellow River Tongren 6 6 6 China Yellow River Longwuhekou 22 22 22 China Yellow River Shuangcheng 18 22 8 China Yellow River Fengjiatai 19 22 8 China Yellow River Tangnailang 18 22 8 China Yellow River Wuchengyi 20 22 8 China Yellow River Zhoujiacun 19 22 7 2.2 Data This section describes data collection and data quality evaluation. Data were gathered in several formats: published lists of sites, tabularized spreadsheets of data, and maps showing data station locations. For each dataset, values were processed and formatted using Microsoft Excel 2010. Additional rows were simply added as more datasets were entered. For this section, the Yellow and Yangtze Rivers will be used as the main examples, with other datasets referenced as needed. 2.2.1 Data Gathering Data sources consisted of two parts: location information for the data stations (coordinates or maps), and tabular or numerical data of the sediment yield values of those 21 stations. Most datasets were collected from previously published research, which meant that these data had to be digitized before they could be added to the Database. For the location information, if coordinates were given, they could be directly entered into a GIS. When maps of station location were provided (e.g., figure 2.1), these images were georeferenced into an existing spatial coordinate system, and the points were then digitized. For numerical data, most sources published data tables (as in table 2.3), which were manually entered into MS Excel. Table 2.3 Example of published tabular data (modified from Lu & Higgit, 1999). No. Tributaries Stations Basin area (km 2 ) Years 1 Jinsha-Yalong Zimenda 133,704 28 6 Shigu 232,651 28 16 Huatan (Qiaojia) 450,696 30 23 Ninnan 3,074 25 25 Qianxinaqiao 2,549 27 29 Meigu 1,607 25 30 Pingshand 485,099 31 31 Hengjiang 14,781 25 32 Zhutuo 694,725 26 In a few cases, however, researchers provided their original datasets (e.g., data from Water Resources Canada, from Church et al., 1999; or data from Hassan et al., 2008, for China). These data were already processed in GIS or spreadsheet format, and could be added to the Global Sediment Database with minimal processing. Similarly, data from the United States Geological Survey (USGS) was downloaded directly from the agency’s website and added to the Database. 22 The data that required the most processing were station locations. In most cases, data station coordinates were not provided, so data points had to instead be digitized by hand, from maps such as those in figures 2.1 and 2.2. These were georeferenced to an existing spatial reference system (using GIS), and relative points were digitized from the georeferenced images. The quality (spatial accuracy) of these maps varied significantly. The Romanian example (from Radoãne & Radoãne, 2005) shown in figure 2.1 was one of the better examples. The map is clear, well-symbolized, and shows secondary information such as rivers and towns to help georeference the points. Figure 2.2 shows a map of Bolivia from Aalto et al. (2006) that is of lower quality. It does not have many candidates for strong georeferencing control points, and even once the map was rectified, its low resolution reduced the precision of digitizing the data points. Map quality was one major consideration in Figure 2.1 Map of sediment yield in Romania (Radoane & Radoane, 2005). 23 evaluating data quality (discussed in detail below). In some cases, such as data from Africa, neither maps nor coordinates were available. Instead, only station names were listed. Reference maps were used to find nearby landmarks for georeferencing. 2.2.2 Data Quality The spatial data (maps of station points) from published sources were not the only data of inconsistent quality. The distribution of these points, and the tabular data associated with the station locations (sediment yield values, temporal information) also varied. For example, data from the Yellow River in China records yearly sediment yield values, as well as information on hydrologic network, and how far a data point is from the mouth of the river. The Romanian data from Radoãne & Radoãne (2005), however, records only a long-term average value. Given the wide range of data quality, a numerical system of evaluating each dataset was established. These rankings form the core of the Data Quality Report. Three metrics were used: spatial precision, Figure 2.2 Data stations in Bolivia (Aalto et al., 2006). 24 temporal resolution, and point density. These are not the only possibilities: other factors, such as quality of metadata, precision of sediment data, and overall completeness could also have been used. Including too many metrics would make the date evaluation system too complicated, and hard to reproduce for other datasets. Further, some these other metrics—especially metadata quality—are harder to quantify, and so are more difficult to assess consistently and objectively. More importantly, it became clear during the data gathering stage that the three factors we chose—spatial precision, temporal resolution, and point density—varied strongly from dataset to dataset, whereas other factors were more consistent. For example, since only yearly suspended sediment data was gathered, the precision of such data was similar for all datasets, whereas the length of record (temporal resolution) varied from as few as a handful of years for some points, to over a century’s worth of data for others. Each of the three chosen metrics will be described in the following section, with data from the Yellow and Yangtze Rivers serving as examples. Evaluating Spatial Precision Each metric was ranked out of four, and then totaled for a score out of 12 for each data source (the higher, the better). The following tables illustrate the criteria that determine each metric’s rank. Table 2.4 (below) summarizes how spatial precision was evaluated. In some cases, precise station locations are not available at all. A station name may be given, or a nearby tributary may be listed, but no coordinates or maps are available to locate the data point exactly. In short, no reliable spatial data is available at all. Such datasets get the lowest rank of 1. Some datasets were accompanied by maps; not all of these were of excellent quality. Lower-quality maps received a rank of 2, while higher-quality maps were ranked as 3. The Bolivian map from Aalto et al. (2006, see figure 2.2), received a rank of 2 because of its 25 limited use in mapping precise station locations. The Romanian example from Radoãne & Radoãne (2005, figure 2.1), on the other hand, gives greater detail, which facilitates plotting the data points. The Romanian map is ranked as 3. Note the presentation of secondary features in the Romanian map, such as rivers, towns, and lakes. These are easy to see and serve as guideposts for locating the data stations. Even when a river is partially obscured by a data point or other feature, it is easy to see where the flow continues on afterwards. Table 2.4 Evaluating spatial precision. Rank Description 4 Coordinates known and/or projected. e.g.,, ready-made GIS shapefile or spreadsheet, such as data from China or the USA 3 Digitized from a high-quality map. This means ≥6 easily matched control points with an overall root mean square error of <0.1 after geo-referencing. Ex. data from Romania, Spain, or Italy 2 Digitized from a low-quality map: control points uncertain, resolution low, and/or overall root mean square error of >0.1 after geo-referencing. Ex. data from Scotland, South America, or Siberia 1 No coordinates, and no map available; but some geographic location is known, such as station name, nearby town, sub-basin name, etc. Ex. Siret basin in Romania, data from Tunisia, Morocco In the Bolivian example, secondary features like rivers and towns are not easily visible. Instead, a shaded elevation maps underlies the data stations. This is less useful in cross- referencing station locations for three reasons. First, using only an elevation map to locate points, especially in a mountainous region like Bolivia, is challenging in itself, as distinguishing one valley or peak from another is not always straightforward. Second, the map itself is less clear, in terms of colour palette and feature clarity, so identifying individual features is difficult. Third, the white circles and lettering that indicate station location obscure the underlying elevation, adding an additional level of complexity to the task of locating data points. In addition to descriptive, qualitative evaluation of data station maps described above, a numerical evaluation was also incorporated. When maps are georeferenced into a GIS, the precision of the control points used is analyzed to estimate the total root-mean-square (RMS) error of the georeferencing process. In general, the more control points used, the better the fit 26 will be, and the lower the RMS error will be. To qualify as rank 3 in terms of data station precision, georeferenced maps must have at least 6 control points and an overall RMS error less than 0.1. Maps that fail to meet these criteria are ranked as 2. Precise latitude and longitude coordinate values are more useful than even high-quality maps. Such datasets were ranked as 4, the highest level. These points can easily be projected using our GIS, which eliminates the need of digitizing points by hand. Ranks 1, 2 and 3 are more likely to have errors, not just from the original data source, but from the digitizing process as well. Evaluating Temporal Resolution The data quality evaluation for temporal precision is shown in table 2.5 (below). The rationale for these criteria is motivated by Church et al. (1989). These authors concluded that a compromise must be made regarding the length of record for sediment yield data points. If the limit is set too high (say, ten years of data), many points will be eliminated from consideration from a dataset. If there are too few points, meaningful spatial analysis is impossible. If the limit is set too low (at two years, perhaps), then the year-to-year variation in sediment yield values will overwhelm any long-term trend in the data. Table 2.5 Sample data stations from the Yangtze River: temporal resolution analysis. Station Name Years with data Total years of data Longest continuous record (years) 5-year average (T/km 2 /year) Average of all years (T/km 2 /year) Dong River 1963-1967; 1969- 1974; 1976-1978 14 6 370 401 PaoJiang River 1965-1970; 1972, 1973 8 6 740 823 YingQianShui River 1964-1985 22 22 83 89 27 The compromise that Church et al. (1989) decided upon for their study of sediment yield in Canada was five years of continuous data. This is a somewhat arbitrary limit. As the authors themselves argued, some data stations were stable enough that even over three or four years of data the overall average was consistent, whereas at other stations longer records did not guarantee stable sediment yield values in the medium or long term. Church et al. (1989) came to this conclusion by examining a sample of individual stations, and assessing how the long-term average sediment yield varied depending on how mean years were selected for that station. A similar analysis was performed for this thesis. Table 2.5 shows an example of three data points from the Yangtze River that were randomly selected and evaluated for temporal resolution. For each data station in table 2.5 the years with data are shown. Note that the Dong and PaoJiang river have multiple, non-overlapping records. The total number of data years is also shown, and is the longest continuous stretch of data. Note for the YingQianShui River there are no breaks in the data record, whereas for the other two stations at least one break exists. Next, a five–year average of the sediment yield of each station was calculated. For the Dong River, which had two stretches of data at least five years long, the earliest continuous period was used. Finally, the average sediment yield for all years is shown in the right-hand column. For all three stations, the long-term average is not significantly different than the five-year average. Similar analyses were performed for all datasets used in this thesis. And since we adopted the recommendations of Church et al. (1989), we used a five-year continuous window of data as our standard. This is a reasonable compromise that does not exclude too many points but nevertheless helps limit the effects of year-to-year variation. Variation between data stations must also be considered. For the three example points presented in table 2.5, two stations have 6 years of continuous data, while the third has 22 years 28 of continuous data. How then, should the overall temporal resolution be evaluated for the Yangtze data? This thesis documents temporal data quality in two ways. First, point-by-point information on length of record, longest continuous period of data, and number of data years will be incorporated into the Global Sediment Database (e.g., the data shown in table 2.5). Potential future users of the Database can select data points that meet their own criteria for temporal data quality. Second, the data quality report for the Database will provide an overall evaluation of each dataset’s temporal data quality. As discussed above, this will be based on Church et al.’s (1989) compromise of five years of continuous data. Datasets will be ranked based on what proportion of their data span at least five continuous years. If no temporal metadata is available, or if fewer than 5% of points in a dataset have five years of continuous data, that dataset is ranked as 1, the lowest level. Likewise, datasets where between 6% and 25% of points meet the five-year criterion will be ranked as 2, and datasets where between 26% and 50% of points meet the five-year criterion will be ranked as 3. The highest level, rank 4, will be awarded to datasets where more than half of the points have at least five years of continuous data (these rankings are summarized in table 2.6). By providing both point-by-point and dataset-level temporal data quality evaluations, the Global Sediment Database will allow users to select the data that is most appropriate for their work. Evaluating Point Density Table 2.6 Evaluation of temporal precision using years of data available. Rank Temporal precision value 4 >50% of stations have at least 5 continuous years of data 3 26-50% of stations have at least 5 continuous years of data 2 6-25% of stations have at least 5 continuous years of data 1 0-5% of stations have at least 5 continuous years of data, or temporal resolution is unknown 29 The number of points in a dataset is also important for conducting reliable spatial analysis, such as interpolation (Bivand et al., 2008). Most interpolation algorithms require between 5 and 12 nearby points when predicting values (ESRI, 2011). For each point to have approximately that many neighbours, about 20-30 points would be needed. With fewer than 20 points, any points at the edge of the dataset will likely have very few nearby points, which would skew the interpolation. To a certain extent, such ‘edge effects’ are potentially problematic no matter how many points there are (ESRI, 2011): the outermost ones will always have fewer neighbours than central points. Nevertheless, the more points there are, the less significant such problems are likely to be. The sheer number of points, however, is not in itself sufficient to ensure adequate data coverage. The distribution and density of the points also matter. If there are large gaps between points, interpolations will be skewed, as the algorithms will have to rely more extensively on guesswork. This is because interpolation methods generally assume a level of spatial autocorrelation (Bivand et al., 2008). That is, nearby points share a geographic similarity, while distant points are less similar. Interpolation works by finding a balance between data value and data position. If large areas exist with no data points, it is more difficult to account for autocorrelation, and the resultant kriging interpolation will be less reliable. As a simple method of determining whether a dataset has adequate coverage in a given basin, point density was used to approximate these spatial concerns. Point density is defined as the number of points in a given dataset divided by the area encompassed by that data (shown in table 2.7). While it does not give a measure of gaps in data it at least helps suggest which datasets might be more useful in developing statistical prediction models at all. Point density 30 data quality information will be added to the Data Quality Report in the same manner as were the spatial and temporal precision metadata. Table 2.7 Evaluating point density. Rank Description 4 More than 200 points per million km 2 3 100-199 points per million km 2 2 50-99 points per million km 2 1 Fewer than 50 points per million per km 2 Now that the three methods of evaluating data quality have been discussed, an example of how a dataset is scored overall will be presented. As discussed above, data for the Yangtze River was available as GIS shapefiles, so precise decimal degree coordinates were known for latitude and longitude. This ranks as 4/4 for data station precision. As demonstrated by the three example stations in table 2.5, there is variability in terms of temporal data quality within the Yangtze dataset. All stations had at least five years of continuous data (stations with <5 years of data were removed before we obtained the data), while some had as many as 24 years of continuous data. Since more than half of the data points have at least five years of continuous data, the Yangtze dataset gets 4/4 for temporal resolution, too. In terms of point density, the Yangtze has 169 points, and its basin is 1,615,884 km 2 . This converts to 104 points per million km 2 , which is rank 3/4 for point density. Thus, the overall data quality of the Yangtze river is 11/12 (i.e., 4+4+3). This is the value that will be recorded in the entry for the Yangtze dataset in the Data Quality Report. Table 2.8 shows the data quality evaluation of the datasets mentioned in this section. Table 2.8 Data Quality Report for selected datasets. Dataset Station location precision (/4) Temporal resolution (/4) Point density (/4) Total (/12) Yellow R. 4 4 3 11 Yangtze R. 4 4 3 11 Romania 3 1 2 6 Bolivia 2 1 1 4 31 However, it should be reiterated that this data quality evaluation is only one of many possible methods, and that its criteria, while based on an analyses of the datasets used in the Global Sediment Database, is somewhat arbitrary. A numerical score out of 12 is a quick way of summarizing different data quality issues of a dataset. But if other users of the Database have different standards for data quality, or require additional metadata, they can refer to the Data Quality Report, and make their own decision of how to rank the data. For instance, users would see that the Yangtze dataset has coordinate information in decimal degrees, projected using the WGS 84 datum; 55% of its data has a record length of at least five continuous years; and that it has 335 points for a basin that is 1.6 million km 2 . If Database users need even more specific information, they can refer to the Database itself, and see which points have a record length of only six years, and which have more than two decades of data. Since detailed metadata is provided in both the Global Sediment Database and the Data Quality Report, users can pick and choose which datasets, and even which data points, are most useful for their purposes. By providing metadata in both the Global Sediment Database and the Data Quality Report, users can compare the quality and potential usefulness of multiple databases. For instance, the Yellow River has the same data quality rank as the Yangtze River (11/12), so researchers may decide that comparing the two would be worthwhile. 2.3 Using the Global Sediment Database to Map Sediment Yield By providing an extensive Database that records the same variables for each dataset and for each data point, we hope that information from the Database can be applied to research projects. As argued in the introduction section, making sediment maps is one possible use of suspended sediment data. This section will explain the methodology used to create a sediment yield map for the Yellow River using the Global Sediment Database. This involves several 32 steps: processing and checking specific sediment yield (SSY) data; establishing sediment rating curves for sub-basins within the Yellow River watershed; and finally interpolating the sediment data using Geographic Information Systems (GIS). First, as described in more detail in the introduction (see figures 1.1 and 1.2), recall that SSY is the amount of sediment, in tons per square kilometre, which is transported from a given location over a year. The Database reports these values for the Yellow River. In addition to considering the quality of the data points, other variables should be accounted for: contributing area (Ad), geology, land use, landscape history, hydrology, climate, vegetation, dams, and other anthropogenic factors can significantly alter the sediment dynamics of a basin (see, for example, Walling, 2005). Because watersheds are heterogeneous, in order to develop meaningful regional relations, it is recommended (Church et al., 1999; Hassan et al., 2008) that smaller, more homogeneous sub-basins be considered for analysis. First, sediment relations for each sub-basin within a dataset must be developed. For the example rivers (Yellow, Yangtze), there are relatively well agreed-upon sub-basins already established (see Hassan et al., 2008 for the Yellow, and Hassan et al., 2011 for the Yangtze), which will be used for the interpolation analysis of this thesis. 33 Figure 2.3 Example of a regional sediment relation (upper reach of the Yellow River). Once a watershed has been divided into approximately homogeneous units, a regional sediment relation is established by first plotting specific sediment yield versus contributing area (figure 2.3). A power-law relation is then calculated for this plot. The exponent from this equation, together with a data station’s specific sediment yield and area, is then used to back- calculate the scaled specific sediment yield ks, using the equation below. L/A is the specific sediment yield (tons/km 2 /year); Ad is the contributing area (km 2 ); and b is the exponent from the regression plot. While the exponent, b, represents the regional trend, which we assume to be constant, ks represents the local variability of the sediment yield. This variability is related to local conditions at each data station. If the rating curve were used by itself, it would mask much of the local variability, producing an artificially smooth surface. By instead back-calculating ks, we preserve an aspect of the local conditions at each site. ks values are then used in the GIS interpolation. Upper Reach y = 8758.5x-0.407 R² = 0.2255 1 10 100 1000 10000 1 10 100 1000 10000 100000 34 2.4 Using GIS to Create Interpolated Sediment Maps This section covers how the data from the Global Sediment Database are used to create sediment maps, using the Yellow River as an example. This includes linking the tabular ks values to digitized, georeferenced point features; choosing the correct interpolation method and statistical semi-variogram; and the execution of the interpolation procedure. The rationale of each step will be discussed below. The outcome of these steps is a series of interpolated sediment yield maps, with accompanying maps predicting the errors of these interpolations. The procedure is summarized in the following list: 1. Select appropriate data, based on desired data quality 2. Obtain (or create) relevant geographic shapefiles (data stations, river network, watershed, sub-basins, etc.) 3. Calculate ks values, and link to relevant point features 4. Select and execute appropriate interpolations 5. Produce the error surfaces associated with these interpolations 6. Trim the interpolation and error surfaces to match the extent of the data 7. Visualize the final results, for both interpolation maps and error maps First, we select appropriate data. The Yellow River dataset is ranked as 11/12 in the Data Quality Report. While this is almost an ideal ranking, it is worthwhile to examine its data quality more closely before proceeding. Like the Yangtze River, the Yellow River scores 4/4 for both temporal resolution and data station precision (all data points have at least five years of data, and decimal degree coordinates are available). In terms of point density, the Yellow River has 161 points for an area that is 874,997km 2 , or 184 points per million km 2 . This ranks as 3/4, hence why the Yellow River is rated as 11/12 overall for data quality. The data quality is high enough that the entire dataset can be used without reservation. Second, we create the needed GIS shapefiles. In addition to the sediment data from the Global Sediment Database, other spatial data are needed before interpolation can proceed. 35 Helpfully, the data provider for the Yellow River, Hassan et al., (2008) provided us with GIS shapefiles for station location, hydrologic network, watershed outline, and sub-basin polygons. Third, we calculated scaled sediment yield values (ks,) using Excel. Since the coordinates were known precisely, as were the unique station names, the Excel data could be added directly to the GIS map document containing the spatial files from step 2. The data are then projected into the same geographic coordinate system as the other shapefiles. Fourth, we need to decide which interpolation algorithm is most suited to the data. This is not a simple step, as there are dozens of options available, and not all of them are appropriate. Table 2.9, below, lists the basic factors that govern interpolations. Considering how these factors relate to sediment yield data represents a first pass in eliminating inappropriate interpolation methods. Table 2.9 Factors of interpolation. Interpolation is… …or Global Local Exact Approximate Stochastic Deterministic Abrupt Smooth Given sub-basins are being used in this thesis to back-calculate ks, global interpolations are not suitable. The back-calculation procedure already helps account for variation within a sub-basin, so merging basins during interpolation is inappropriate. Similarly, since ks represents the variability of sediment yield, which is itself related to local conditions at each data station, exact interpolation methods are preferable. This will ensure consistency between the interpolation surface and the original sediment yield relation. In short, the process of dividing data into sub-basins, producing a sub-basin level sediment relation, and performing back- calculations to help account for local conditions addresses the concerns of global vs. local and 36 exact vs. approximate. Since these factors are already accounted for, they can be set aside for the interpolation. In a similar way, there are several reasons why a stochastic interpolation method is preferable to a deterministic one. First, as explained above, a stochastic method (using a power- law fit to make the sediment yield relation) was already used in the back-calculation, so it is appropriate to continue in this direction as the underlying assumptions of stochastic and deterministic processes are different. Since the back-calculation already implicitly assumes a stochastic component to sediment yield, there is little gain in adding a deterministic interpolation. Further, a deterministic approach would require more knowledge about the key variables that affect sediment yield. While we do possess some secondary information (such as slope, aspect, discharge, elevation, hydrologic network, or physiographic region), we do not know how to fully model these variables, and to attempt to do so would be in contrast with the earlier back-calculation approach. Regarding the smoothness of the data, there are two things to consider. First, it is again significant that the back-calculation has been performed at the sub-basin level. It is thus expected that there would be abrupt demarcations between adjacent sub-basins. In other words, the interpolation is not intended to predict across watershed boundaries. Second, we have already assumed that the exponent, b, of the sediment relation is constant for a given sub-basin; this is how the back-calculation of ks was executed. By using an interpolation method that allows for smooth surfaces, we can check this assumption against our interpolation surface. Therefore, the ideal interpolation for mapping sediment yield will be flexible enough to allow for both smooth and abrupt surfaces. 37 Kriging interpolation meets all these criteria: it is exact and stochastic, and it can be adjusted to be either global or local, abrupt or smooth. In addition, kriging has several inherent strengths that are useful when interpolating riverine data: it considers spatial autocorrelation, it analyzes the overall statistical trend of nearby points, and it is less sensitive to changes in spatial scale. Finally, and most importantly, kriging allows for the generation of error surfaces, so the interpolation predictions can be checked against the original data. When data are spatially autocorrelated, data points are not exclusively independent. That is, there is potentially a relationship between the value at point A and the value at point B (for instance, point B might be downstream from A, and thus should have a higher sediment yield than A). Crucially, kriging does not immediately assume the directionality or magnitude of this relationship (as is the case with most other interpolation methods). Rather, it examines both the overall spatial trend (i.e. of all points) and the local statistical distribution of points to assess what kind of relationship exists. If the underlying pattern is known (for instance, that spatial autocorrelation follows an exponential pattern, or a linear pattern) specific statistical models (semi-variograms) can be selected to fine-tune the kriging process. Kriging also allows one to fine-tune how the regional statistical pattern is assessed. For example, one can specify the number of nearby points the regional statistical prediction should consider, or set a maximum distance (radius) of the regional statistical trend. This combination of scalability and statistical flexibility makes kriging ideal for a variable such as sediment yield, which may exhibit inconsistent, and even contradictory, degrees of spatial variation. Previously published sediments maps, such as for Canada (Church et al., 1999) and the Yellow River (Hassan et al., 2008) have used kriging to interpolate SSY. Both these examples argue that kriging is the preferred method because it is equally capable of handling both large 38 and small scales, even within a single interpolation. In other words, it is not as sensitive to the geometry of point distribution as other methods. Further, the user can select from a variety of different semi-variograms, each one assuming a different model of spatial autocorrelation. This customizability means kriging can be used in a wide variety of context, especially when auto- correlation values are known. Figure 2.4 Graph showing simple and ordinary kriging (modified from ESRI, 2011). There are three basic kinds of kriging: simple, ordinary, and universal. Each of these in turn has multiple semi-variogram options. Kriging type has a strong influence on the way the overall interpolation surface is generated. Semi-variograms tweak these overall trends to consider different patterns of auto-correlation. This section will examine these options, and narrow the possible interpolation methods given the kind of data under examination (suspended sediment yield). Simple and ordinary kriging are presented in figure 2.4 (above). Both these methods use the formula Z(s) = µ + ε(s), where s is data point location, Z(s) is the measured value (the back- calculated ks) at position s, µ is a constant function, and ε(s) is the difference between Z(s) and µ. 39 Significantly, in simple kriging µ must be known; in ordinary kriging, µ can be unknown. In other words, both methods assume that there will be a constant mean throughout the dataset. Simple kriging requires that this mean is known, while in ordinary kriging the mean is unknown. Note that µ does not represent the average of just the measured values, but of the entire landscape. In other words, if you had data point for every square metre of a watershed, µ would be the average of those points. Clearly, for sediment data, it is unreasonable to assume that µ is known. But is µ constant? Is it the same ubiquitously throughout the landscape? This question is harder to assess. Even if µ turns out to be variable, it is still instructive to perform ordinary kriging. Not only may this help demonstrate whether µ is indeed constant, but ordinary kriging may highlight underlying spatial trends that can be useful in conducting more advanced interpolation techniques. A second motive suggests the use of ordinary kriging. Both of the previously published sediment maps (Church et al., 1999; Hassan et al., 2008) used ordinary kriging. While this is not justification in itself, by using the same methods as these researchers, the effectiveness of ordinary kriging for sediment mapping can be assessed. More importantly, the resulting maps may suggest whether or not µ is constant for each dataset. If µ is indeed variable then universal kriging is appropriate. Figure 2.5 illustrates the statistical model. The equation describing universal kriging is a bit different from the equation for simple and ordinary kriging: Z(s) = µ(s) + ε(s). In this equation, µ–the constant mean– has become µ(s), a function which describes a changing mean. The error term is again ε(s). By examining both ordinary and universal kriging, this thesis will assess whether the example datasets exhibit constant or variable averages. Further, by exploring more than one interpolation method, this thesis will demonstrate that not all interpolation methods are equally valid. Thus, 40 the choices made when executing a GIS analysis tool are significant, and should be fully documented in any published research. An example of a universal kriging interpolation is shown in figure 2.6 (next page). Figure 2.5 Graph showing statistical model of universal kriging (modified from ESRI, 2011). Now that ks values have been interpolated, error surfaces can be created (step 5). Recall from figures 2.7 and 2.8 that the kriging formula is Z(s) = µ + ε(s), where s is data point location, Z(s) is the measured value (the back-calculated ks) at position s, µ is a constant function, and ε(s) is an error value, or the difference between Z(s) and µ. These values can be mapped, as for each cell in a kriging interpolation a prediction value and an error value are created, using this formula. The same statistical toolkit that produces the interpolation surface can be used to produce an error surface for that interpolation. Thus, for each interpolation method, there will be a corresponding error surface. These pairs of maps (interpolation/error) can be analyzed in tandem to illustrate and evaluate the spatial pattern of the sediment maps. Figure 2.7 shows the error surface from the interpolation in figure 2.6. 41 Figure 2.6 Map of the Yellow River, showing the rectangular kriging surface. Figure 2.7 Error surface for the Yellow River. The interpolation and error surfaces described above should not immediately be mapped, however. Interpolation and error surfaces are produced as rectangles, which correspond to the maximum extent of the input dataset. Thus, we need to trim these surfaces to fit the boundaries of the watershed’s boundary (step 6). Figure 2.8 shows a screen capture of this, again showing 42 the Yellow River. By extracting only those portions of the interpolation surface that fall within the watershed boundary, a more reliable map can be produced. Figure 2.8 Yellow River interpolation surface, trimmed to watershed boundary. This thesis goes farther than making sediment maps and error surfaces. As each kriging method has different statistical models available (different semi-variograms), this thesis will also explore these options. Each one describes the kind of fit—µ, or µ(s)—that is applied to the input data. Further, both error surfaces and interpolation surfaces can be statistically compared. Thus, by comparing successive pairs of maps, the most ideal interpolation method for our dataset can be established. To ensure the comparability of these maps, the visualization techniques for each interpolation (step 7) need to be consistent. The colour scheme presented for the prediction surfaces in figure 2.6 and 2.8, for instance, will remain constant for similar maps. Further, the class breaks need to remain consistent so that the same shade of brown in figure 2.6, representing 2,501–5,000 T/km2/yr, is exactly the same hue, for the same range of values, in figure 2.8. A 43 consistent visual scheme will also be used for the error maps. Finally, as discussed above, it is important to keep in mind that kriging interpolation is most accurate closest to the data points— hence why the surface in figure 2.6 was clipped to match the watershed boundary, resulting in figure 2.8. The data points and the hydrologic network were thus added to these maps to help the reader orient the prediction surface in the landscape. This section described the methodology used in creating the Global Sediment Database and the Data Quality Report. The process of gathering data, evaluating data quality, calculating sediment yield, and creating sediment yield maps was also presented. The next section will expand on the topics of data quality and metadata explicitly, through a series of prediction and error maps depicting the Yellow River. By comparing statistical measures, prediction surface, and error surfaces from different kriging variations, the most appropriate interpolation method will be determined. 44 Chapter 3 Analysis & Discussion 3.1 Background Information on the Yellow River To evaluate the different kriging interpolation methods that were used to develop the sediment yield maps of the Yellow River, it is important to discuss some background information on the Yellow River basin. Geology, physiology, and climate will be described, and regional sediment yield relations (like those presented in Chapter 2) for individual homogenous sub-basins and the overall watershed will be analyzed. Figure 3.1 Yellow River reference map. Figure 3.1 shows the overall layout of the Yellow River and its main tributaries. It flows from west to east, with a delta at the Bohai Sea. The entire basin is approximately 750,000 km 2 , and Hassan et al. (2008) defined 10 homogenous sub-basins within the basin. A desert region, found in the middle of the watershed, was not included in a sub-basin, as it is generally assumed (Wang et al., 2006; Hassan et al., 2008) that the desert is not hydrologically connected to the rest 45 of the region—it drains internally, and does not connect to the main stem. We thus exclude the desert from further analysis. The Yellow River can also be divided into three physiographic regions: the Tibetan Plateau, the Loess Plateau, and the North China Plain (figure 3.1). The description of these regions in the following three paragraphs is based on published material (e.g.,, Long & Chien, 1986; Xu & Yan, 2005; Wang et al., 2006; and Hassan et al., 2008). The uppermost regions (basins 1-3 – figure 3.1) are part of the high Tibetan plateau, where the main stem is deeply incised into the underlying bedrock. The middle regions (basins 4-8) from the Loess Plateau. Loess material is the primary sediment source for the Yellow River, and is between 100m and 200m thick in this area. Overall, Loess sediment covers ~37% of the Yellow River watershed. The lower sub-basins (basins 9 and 10) flow through the North China plain. This section of the river is confined by levees, and there are few tributaries here. Much of the Yellow River basin is semi-arid (including the dry Loess regions), with an annual precipitation of 478mm (this does not include the excluded desert region, mentioned above). Despite low rainfall, the Yellow River is known for some of the highest sediment yield values in the world, due to the high erodibility of the Loess sediments. The historic hydrological regime has been substantially altered, as more than 180 large dams are present within the basin (>10,000,000 m 3 in volume). The construction of these dams may be one reason why sediment yield has been much lower in recent years (Walling, 2006). Other potential factors include soil conservation practices and changes in land use along the river. Such changes in sediment dynamics are not uniform throughout the basin. Overall, the uppermost sub-basins (1-3) have remained more consistent, while the Loess regions (sub-basins 4-7) have varied more substantially over time. The analysis of these changes is based on the 46 work of Hassan et al. (2008), who created regional sediment relation for each sub-basin. As described in Chapter 2, Hassan et al. (2008) struck a balance between sub-basin level consistency and individual station variation. The regional sediment relation exponent, b, was assumed to be constant for each sub-basin. The back-calculation procedure used both b and each station’s sediment yield data to produce scaled ks values at each station. The statistical details on the power relations established for the 10 sub-basins of the Yellow River appears in table 3.1. These values are from Hassan et al. (2008), and the following statistical discussion is based on these values. The sign of b (positive or negative) indicates whether the basin is aggrading or degrading: a positive b means sediment yield increases with area, while a negative b means sediment yield decreases with area. When b is close to zero (e.g., sub-basins 1, 6, and 8), there is no scale effect. Table 3.1 Sediment relation details for Yellow River sub-basins. Modified from Hassan et al., 2008. Original values are in log units; those here are back-transferred. Sub-basin ks intercept (T/km 2 /yr) Confidence interval of ks b (exponent) Confidence interval of b r 2 1 464.52 79.6 –0.036 0.536 0.001 2 58,210.3 10.4 –0.677 0.292 0.59 3 3.92 4.12 2.029 0.983 0.86 4 751.6 2.28 0.251 0.025 0.49 5 10.09 3.37 0.813 0.176 0.51 6 1,570.4 10.6 –0.026 0.320 0.00 7 101,624.9 1.97 –0.290 0.087 0.35 8 1,888.0 2.24 0.002 0.104 0.04 9 223.9 2.34 0.184 0.106 0.25 10 3,614.1 2.12 –0.230 0.095 0.84 The confidence interval of each of these values is also shown. For nine sub-basins, the ks is larger than the corresponding confidence interval. This indicates a high degree of reliability. For sub-basin 3, however, the confidence interval of the ks intercept is larger than the value itself (Table 3.1). This is because there are only five measurements within the sub-basin, which are widely scattered. But the r 2 is still high (0.86), because the sub-basins spans a wide range of 47 contributing area. It warrants examining this region on the maps that will be presented later in this section. In the case of b, four sub-basins (1, 2, 6, and 8) have a confidence interval larger than b itself, although these were reported as statistically significant at α = 0.1 by Hassan et al. (2008). If the kriging interpolation for these regions shows anomalies, these statistics may need to be considered. If the kriging surface suggests a different spatial pattern of sediment values for these regions, they might need to be analyzed separately. For region 1, the lack of data points is one reason the confidence interval of b is so high. However, b itself is close to zero, so the overall effect on the back-calculated ks is minimal. Sub-basins 6 and 8 similarly have b which is close to zero, so their wide confidence intervals for b are not as significant. Figure 3.2 Yellow River sediment relations. Moving from the individual sub-basins to the watershed as a whole, we find that for the majority of the Yellow River there is no relationship between sediment yield and basin area. Out of 161 data plotted in Figure 3.2, 33 lie along the main stem of the Yellow River, and have been plotted separately (in blue). The main stem stations, however, show a steep increasing trend with drainage basin area implying that they differ from the rest of the stations located on the Main Stem y = 1E-06x1.5907 R² = 0.7433 Tributaries y = 1959.7x0.0114 R² = 0.0001 1 10 100 1000 10000 100000 1 10 100 1000 10000 100000 1000000 Yellow River Sediment Yield Relations 48 landscape. As Hassan et al. (2008) pointed out, this significant trend in the main stem data suggests that sediment yield from the land surface increases systematically as one moves downstream from the confined headwaters to the Loess Plateau and into the Northern Plains. In short, there is no overall scaling effect for the Yellow River, considered as a whole. Using the exponent (b) values in table 3.1, and after examining the overall data distribution (shown in figure 3.2), Hassan et al. (2008) proceeded with the back-calculation procedure. Having similarly checked that the values from the Global Sediment Database match those of Hassan et al. (2008) our examination of the Yellow River’s data quality is complete. The ideal kriging algorithm can now be determined. 3.2 Rationale for Selecting Kriging Algorithms Recall from Chapter 2 that kriging is a spatial interpolation method for estimating values of a variable at locations where those values are not specifically known. Kriging computes statistical relations between known data points to predict data values across a landscape. There are many possible kriging algorithms that could be used to interpolated sediment yield values for the Yellow River, but which is best? Each of the kriging options explored in this section (kriging type; range, nugget, and sill; semi-variogram) can be combined in different ways to produce literally dozens of kriging algorithms. Most of these will be unsuitable. For this thesis, the “best” option means the most accurate, and most suited to represent the spatial pattern of our data. To determine the most accurate kriging algorithm the overall error values of possible methods will be examined. Each possible algorithm is called a semi-variogram (Armstrong, 1998). From the 10 possible options outlined as possibly appropriate in Chapter 2, the four with the smallest error will be examined in greater detail, by analyzing how the shape of each model 49 fits the source data. After these two rounds of elimination, two final semi-variograms will be compared. Sediment yield prediction and error maps will be created for both, so that the underlying spatial patterns of sediment yield and error values can be assessed. In this way a single kriging semi-variogram will be identified as the most suitable for the Yellow River data. Following the analysis of the Yellow River maps, two other example datasets from the Global Sediment Database will be briefly analyzed for comparison purposes. This will help suggest a decision-making scheme so that the best kriging semi-variogram can be determined for any sediment yield dataset. Before examining kriging options in detail, let us review how the data were first processed. As the preceding section demonstrates, using ‘raw’ data from the Global Sediment Database without first examining data distribution, statistics, and quality, is not recommended. The flowchart in figure 3.3 summarizes how data are prepared for analysis. 50 Figure 3.3 Flow chart—preparing data for spatial analysis. The first steps—selecting a dataset from the Global Sediment Database, and investigating its data quality—were described in the first part of this Chapter. More detailed procedures are not outlined in figure 3.3 because different projects will have different goals, and different thresholds for data quality. Recall that the Data Quality Report that accompanies the Global 51 Sediment Database contains detailed and summary data quality information, so users can put as much detail into their metadata analyses as needed. In most cases the data quality check at the beginning of this flow chart will require several iterations to ensure adequate data quality. The remaining sections of the flow chart follow two branches, dealing with spatial data and sediment data. Evaluation, analysis, and processing of sediment data was discussed in Chapter 2 and the earlier part of this Chapter (creating regional sediment relations, back- calculating ks, and so on). Regarding spatial data, the discussion on the evaluation of spatial precision (Section 2.2.2) summarizes how to deal with the main challenge of spatial data processing: obtaining accurate data point GIS files, and determining the ideal spatial coordinate system. For this thesis, the data files and coordinate system used by Hassan et al. (2008) were adopted. Given the variability of the SSY data for the Yellow River, it is important to consider the appropriateness of the kriging algorithm. As discussed in Chapter 2, there are three possible methods of kriging: simple, universal, and ordinary. Simple kriging is not appropriate for sediment data, so we have the various semi-variograms of universal and ordinary kriging to choose from. Several parameters need to be considered. First, the basic kriging features—range, nugget, and sill—will be examined. Next, the various semi-variogram options will be considered. By discussing these factors, we will narrow the possibilities of appropriate kriging algorithm for the Yellow River to four candidates. These four will be subject to further evaluation. 52 3.2.1 Range, Nugget and Sill Figure 3.4 Range, nugget, sill, and partial sill in kriging (modified from ESRI, 2011). Figure 3.4 illustrates the first of these parameters: range, nugget, and sill (ESRI, 2011; Negreiros et al., 2010; Armstrong, 1998); it is a hypothetical graph of a kriging algorithm. A data point is located at the origin (0,0) and as distance from that point increase (x-axis), the variable y likewise increases. At a certain point, the kriging graph levels off—that is, there is no further change in y with increasing distance x. This y-value at this point is called the sill. The range is the distance away from the origin at which this leveling off occurs. Data points farther away than the range are not considered to be spatially autocorrelated (as the kriging equation would continue as almost a horizontal line). Larger ranges mean that more of the landscape is subject to autocorrelation. Larger sills mean that the magnitude of the autocorrelation is high. The nugget is the vertical offset of the kriging equation from the origin. This represents the minimum autocorrelation value. In other words, even at an infinitesimally small distance from our origin point, there is some change in our variable as a result of spatial autocorrelation. Finally, a partial sill is also sometimes reported: this is simply the sill minus the nugget. 53 Unless one has specific mathematical knowledge about how a spatial variable changes with distance, it is not possible to predict these values ahead of time. That is, one cannot ‘measure’ the sill or range of an expected interpolation directly from the landscape. Once a kriging algorithm has been applied, however, one can examine how well the interpolation ‘fits’ the data, and evaluate whether the range, nugget, and sill values match the data (see figure 3.5). Figure 3.5 Example kriging algorithm of the Yellow River, with statistical details (see text). In this example, ordinary kriging was used. The y-axis shows the variation in sediment value, while the x-axis is the distance from the data point being interpolated. This graph shows an increasing trend, with higher values occurring farther away from the data station. The red dots represent other data stations’ location (x-position), and how each of these stations compares to the target station (y-position). Further, averages for these data stations (red dots) are binned at equal intervals, with the average of each bin shown as a blue cross. Figure 3.5 thus shows the surrounding data stations in three ways: in terms of real positions and sediment value deviations 54 (red dots); spatially-averaged values, moving away from the target station (blue crosses); and an overall trend line (blue line). The range, nugget, and partial sill are highlighted in the inset table, which shows the statistical properties of this algorithm. There is a great deal of scatter in this example and the algorithm seems to underestimate the sill: the high degree of scatter in the upper-right of the graph is not accounted for by this algorithm. This would suggest that this model is less accurate farther away from the target station; specifically, interpolated values will be too as compared to positions closer to the target station. 3.2.2 Considering Error Values In addition to examining the graph of the kriging algorithm, the error values of the interpolation can be examined. When the kriging algorithm is generated, the software produces several statistical measures that describe the overall reliability and potential errors of the model. Figure 3.6 uses the data from figure 3.5, and shows a graph of predicted values (on the y-axis) vs. measured values (on the x-axis). The blue line is the regression function, while the green line is the same function, but forced through the origin (i.e. with a nugget of zero). Several statistical values are computed using these regression functions, which are discussed below. 55 Figure 3.6 Kriging error values—predicted values (y-axis) vs. measured values (x-axis). By examining these values we will assess the accuracy of this interpolation algorithm and compare it with other algorithms to decide which is most appropriate. According to the documentation accompanying this kriging tool (the Geostatistical Wizard within ESRI’s ArcGIS, ESRI 2011), the mean error and mean standardized error should both be as small as possible, the root-means-square standardized error should be as close to 1.0 as possible, and the average standard error and root-mean square error should be approximately equal. To demonstrate how such values can be used to decide between two kriging options, consider table 3.2. On the left, y-axis: Predicted values x10 5 x-axis: Measured values x10 5 56 we see the error details for the same kriging function seen in figure 3.4 and 3.5: ordinary kriging, with a stable semi-variogram. On the right, we see the same details for a different kriging function: universal kriging, again with a stable semi-variogram. Each of these will be described in turn, to asses which kriging interpolation is more accurate. Table 3.2 Comparing error values of ordinary and universal kriging. Error Statistic Ordinary kriging Universal kriging Mean error 650.4 301.3 Root-mean-square error 28,362.7 30,243.0 Mean standardized error 0.00955 0.000269 Root-mean-square standardized error 0.9379 0.9628 Average standard error 33,120.9 31,988.7 As discussed in Chapter 2, kriging algorithms produce interpolated values based on an equation, where the predicted value equals the measured value plus an error term. The mean error is simply the average difference between predicted values and measured values: lower values are preferred. The mean standardized error computes the error term at each data point, and then divides it by that data point’s value. In other words, if point x has a ks value of 1,000 T/km 2 /yr, and an error value of 50 T/km 2 /yr, then the mean standardized error would is 0.05 (50 / 1,000 = 0.05). Again, lower values are preferred here. The root-mean-square standardizes error (also known as the RMS standardized error, or R 2 ) describes how well the kriging function fits the data. A perfect R 2 or 1.0 would mean that all variation within the data is explained by the kriging function. An R 2 substantially less than 1.0 means that values are systematically under- estimated; while an R 2 substantially more than 1.0 indicates systematic over-estimation. The remaining statistics, root-mean-square error and average standard error, should be approximately equal. The RMS error indicates how closely the predicted values match the measured values, while the average standard error is computed based on error terms alone. If there is a substantial amount of measurement error, these two values are less likely to be equal. Such a situation suggests a systematic error or bias of the measured data. If the values are 57 approximately the same, then the errors produced by the kriging algorithm are assumed to be stochastic, not deterministic. This is the preferred outcome. In considering table 3.2 again, the error statistics all indicate that the universal kriging algorithm is more accurate than the simple kriging algorithm. The universal method has lower mean error, lower mean standardized error, and RMS-standardized error close to 1.0, and RMS and average standard error that are approximately equal. For the ordinary kriging, these values were 33,120.9 vs., 28,362.7–a difference of 3,758.2. For the universal kriging, the difference was only 1,745.7 (31,988.7 vs. 30,243.0). This analysis of universal and ordinary kriging error values suggests that universal kriging may be more appropriate for the Yellow River data. We will see whether this is true for other basins in section 3.4, when data from the Yangtze River and from Romania are examined. But given the comparison between universal and ordinary kriging in section 2.4, the preferred suitability of universal kriging for the Yellow River makes sense. Recall that ordinary kriging requires a constant, unknown mean, while universal kriging allows the mean to vary. We are already assuming that sub-basins within the Yellow River are different from each other by creating regional sediment relations and back-calculating ks. It follows that average sediment yield values might similarly vary in different parts of the watershed. Despite the slight advantage that universal kriging has at this point, however, ordinary kriging should not be excluded until semi-variography is considered. 3.2.3 Semi-Variography Semi-variograms allow for different mathematical functions to be used as kriging algorithms. Some of these may reflect the data better than others as different models will result in different predicted values. This is especially true when the model differs significantly from 58 the data near the origin. Further, the steeper the model is near the origin, the more influence the data points closest to the target data station will have on the output interpolation. In other words, a stronger localized bias in autocorrelation occurs with semi-variograms that are steep near the origin. Figure 3.7 shows example semi-variograms. Note how most are somewhat similar farther away from the origin, but differ more substantially (in steepness and shape) close to the origin. Figure 3.7 Kriging semi-variogram options (modified from ESRI, 2011). 59 Spherical, circular, and exponential semi-variograms differ more in range than they do in shape. The spherical semi-variogram reaches its sill value quite quickly, whereas the other two are more gradual. The Gaussian option is the only one that curves inwards at first; it then reverses its shape and reaches its sill at about the same place as the circular semi-variogram. The linear semi-variogram is the simplest, but perhaps the most divergent from its data. In order to evaluate which of these semi-variograms is best suited to the Yellow River data, each of these five semi-variograms will be compared, for both universal and ordinary kriging (10 semi-variograms in total). The error values will be compared, as in table 3.2. The two best semi-variograms from each kriging method will be analyzed in more detail: the shape of these four semi-variogram graphs will be examined, to further refine which is most suitable. After these two rounds of eliminating less-suitable algorithms, the semi-variograms with the best fits will be used to generate sediment maps and error maps. These maps will be further analyzed in the hopes that a single kriging algorithm emerges as the best choice for mapping the sediment data of the Yellow River. Table 3.3 Yellow River ordinary kriging error values. Error type* Circular Spherical Exponential Gaussian Linear Mean Error 777.22 827.39 1,086.16 301.33 853.87 RMS Error 28,669.3 28,558.0 28,598.5 30,243.0 28,479.3 Mean Std. Error 0.010687 0.011987 0.016473 0.000269 0.012404 RMS Std. Error 0.95097 0.94708 0.96046 0.98283 0.94672 Avg. Std. Error 32,754.2 32,866.6 34,731.4 31,988.7 33,121.16 Avg. Std. – RMS Error 4,084.9 4,308.6 6,132.9 1745.7 4,641.9 * RMS = root mean square; Std. = standard; Avg. = average. Bolded values are the best fit; underlined are the second-best fit. Table 3.3 shows the various error statistics from the five ordinary kriging semi- variograms. Using the rationale described above, the best semi-variogram for each metric was evaluated. Note that more decimal places were kept than are probably statistically significant; error values will be rounded once the ideal kriging algorithm has been determined. Also note the 60 use of abbreviations in table 3.3, and that the final row shows the difference between the average standard error and the root-mean-square error. As noted above, these values should be approximately equal. So the semi-variogram with the smallest difference in these error values is more accurate than the others. The Gaussian semi-variogram was the most ideal based on all five error measures (these have been bolded in table 3.3). The second-best error values have been underlined. Four of these second-place error values are from the circular semi-variogram: mean error, RMS error, mean standard error, and average standard error. The RMS standard error of the circular semi- variogram is only third-best, at 0.95097. The exponential semi-variogram, however, reports an RMS standard error of 0.96046. However, the other values from the exponential semi-variogram do not fare as well: in fact, for all other error values, the exponential semi-variogram is the worst option. Thus, the Gaussian semi-variogram, as the clear winner, and the circular semi- variogram, as the runner-up, will be examined in more detail. First, however, we repeat the elimination process with the universal kriging semi-variograms (table 3.4). Table 3.4 Yellow River universal kriging error values. Error type* Circular Spherical Exponential Gaussian Linear Mean Error 771.30 855.80 1,230.68 513.06 914.47 RMS Error 28,638.0 28,492.1 28,607.6 30,318.2 28,478.4 Mean Std. Error 0.011623 0.013827 0.020226 0.0075560 0.0153084 RMS Std. Error 0.95070 0.95945 0.96005 0.98659 0.94473 Avg. Std. Error 32,725.4 32,848.5 34,722.8 31,962.5 32,291.0 Avg. Std. – RMS Error 4,087.4 4,356.4 6,115.2 1,644.3 3,812.6 * RMS = root mean square; Std. = standard; Avg. = average. Bolded values are the best fit; underlined are the second-best fit. The Gaussian semi-variogram is again the most ideal (values, as in table 3.3, are in bold). The RMS Standard error is even closer to 1.0 than the ordinary Gaussian semi-variogram (0.98659 for the universal, 0.98283 for the ordinary), and the difference between the average standard error and the RMS error is also smaller (1,644.3 for the universal, and 1,745.7 for the 61 ordinary). The mean error and mean standard error are also lower. A second-best semi- variogram is less simple to identify, as was the case with the ordinary semi-variograms. Whereas the ordinary circular semi-variogram was second-best in nearly every error statistic, the universal circular semi-variogram does not rank as highly. The mean error and mean standard error are indeed only second to the Gaussian semi-variogram, but the RMS standard error is only better than the linear semi-variogram. Oddly, though the linear semi-variogram has the RMS standard error furthest from 1.0, at 0.94473, its RMS error and average standard error are closer to each other than all the other options, save the Gaussian semi-variogram. So which model is second best overall? Table 3.5 helps narrow the decision. We created a ranking for each universal semi-variogram error statistic, where 1 is the best value and 5 is the worst value. That means the lowest mean error would rank as 1, and the highest (worse) as 5. Likewise, the difference between the average standard error and the RMS error should be as small as possible; the best value always gets assigned a rank of 1. Table 3.5 Yellow River universal kriging error ranking. Error type* Circular Spherical Exponential Gaussian Linear Mean Error 2 3 5 1 4 Mean Std. Error 2 3 5 1 4 RMS Std. Error 4 3 2 1 5 Avg. Std. – RMS Error 3 4 5 1 2 * RMS = root mean square; Std. = standard; Avg. = average. Bolded values are the best fit; underlined are the second-best fit. Depending on which statistical measure of error we emphasize, different semi- variograms may seem suitable as second-best. RMS standard error, for instance, is a simple metric that summarizes how well the kriging function accounts for variation within the predicted results. The exponential semi-variogram has the second-best value (0.96005), which means the ‘fit’ between the data and the semi-variogram is reasonably strong. However, this semi- variogram ranks as the worst for all three of the other statistical measures. The circular or 62 spherical semi-variograms are more balanced; some errors statistics rank well, others poorly. Of these two, the circular semi-variogram performs the best, as it is the only semi-variogram to have two second-place error statistics. The circular semi-variogram, then, along with eh Gaussian semi-variogram, will be selected as the two best universal semi-variograms, and will be subjected to further analysis. 3.2.4 Comparing Universal and Ordinary Semi-Variograms As explained above, the shape of a semi-variogram is important, especially when close to the data point in question. In other words, there is a significant amount of statistical weight put on the function close to the origin. If the fit between the semi-variogram and the data is not sufficiently accurate, then the whole assumption of kriging—that interpolated values vary in some predictable way in accordance with their spatial autocorrelation—is false, and the resultant prediction map is unreliable. Figures 3.8 through 3.11 show graphs of the four semi-variograms identified in the previous section that were determined to be most ideal for creating a sediment map of the Yellow River. These are the circular and Gaussian semi-variograms for both universal and ordinary kriging. For each figure, the x-axis is the distance away from the data point, and the y-axis is the back-calculated ks. 63 Figure 3.8 Ordinary circular kriging semi-variogram. The red dots in Figure 3.8 of the ordinary circular option show neighbouring data points, while the blue crosses show the average ks value at various distance intervals from the origin. The blue line is the semi-variogram itself. Note the region circled in orange. It lies close to the origin, and highlights a cloud of data points that lie substantially below the semi-variogram’s curve. The right-hand side of the graph shows a dispersed cloud of point, with the semi- variogram passing approximately through the middle of this cloud. When one considers the average points, however, there is an odd dip in the region circled in purple. A dense cluster of data points (in red) have pulled that region’s average values (blue crosses) down, so that the semi-variogram overestimates this region. While the discrepancies highlighted by the orange and purple circles are not extreme, they highlight a more significant concern: the semi- variogram’s curve does not fit the overall shape of the data very well. The circular semi- variogram is just too simple. The shape of the average points (blue crosses) suggests a more undulating equation. 64 Figure 3.9 Ordinary Gaussian kriging semi-variogram. The ordinary Gaussian semi-variogram appears in figure 3.9. Its curve fits the shape of the data points more closely at the lower end of the graph (orange circle), even though it still under-represents the middle region (purple circle), just as the circular semi-variogram did. In fact, the two semi-variograms are relatively similar in this region. At its furthest values, the Gaussian semi-variogram fits the average data values (blue crosses) slightly better than the circular semi-variogram. Overall, the shape of the Gaussian semi-variogram more closely matches the shape of the data points and, especially, the average values. This supports the earlier statistical analysis (table 3.3), which found that the ordinary Gaussian semi-variogram was the best statistical match. 65 Figure 3.10 Universal circular kriging semi-variogram. We now consider the universal kriging examples. Figure 3.10 shows the circular semi- variogram. There are again two main regions where the semi-variogram curve does not match the data, just as was the case with the two previous examples. The circular semi-variogram does not fir the slight curve of the data near the origin (orange circle), nor does it match the average values. And in the same vein as the previous semi-variograms, the middle region is again under- estimated by the semi-variogram’s curve. Two regions of the Gaussian semi-variogram in figure 3.11 need consideration. Near the orange circle, the prediction values match the shape of the data better than the circular semi-variogram in figure 3.10. But neat the purple circle, there is a potentially significant deviation in the prediction curve, although it is slightly less pronounced than in the circular semi-variogram. 66 Figure 3.11 Universal Gaussian kriging semi-variogram. The conclusion is the same with universal kriging as it was with ordinary kriging: the Gaussian semi-variogram is the best fit, both in terms of statistical measures (table 3.4 and 3.5), and in comparing the shape of the algorithm’s curve to the interpolated data (figures 3.10 and 3.11). These two options—universal and ordinary Gaussian kriging—will be analyzed cartographically in the next section to finally determine which single method is most suitable to creating a sediment map of the Yellow River. 3.3 Analyzing Sediment Maps and Error Surfaces Each sediment map (figures 3.12 and 3.14) is accompanied by its corresponding error surface (figures 3.13 and 3.15), derived from its kriging equation (interpolated value = measured value + error), using the Gaussian semi-variogram. The symbology for these maps is consistent. Since sediment values varied from less than 100 to over 200,000 T/km 2 /yr, a semi-log classification scheme was used, so that it was easy to visually identify each order magnitude on 67 the maps. Class breaks occurred at each half log interval, plus the first quarter log interval (e.g., 100, 250, 500, 1,000). The final classification scheme has twelve classes, ranging from a minimum of 0 to a maximum of over 230,000 T/km 2 /yr. Given the maximum value did not quite reach the next quarter-log mark of 250,000 T/km 2 /yr, a lower class break of 150,000 T/km 2 /yr was used instead. Variation in chroma is generally the preferred method for visualizing quantitative variables like sediment yield, but since these maps require spanning five orders of magnitude, it would be too difficult to see the subtle changes a monochromatic colour scheme would require. Our chosen palette instead is thematically linked to sediment yield, with blue vales representing low yield (‘clear’ water) and progressively yellow and brown colours representing higher sediment yield (‘cloudy’ water). This hue progression scheme allows for map clarity across the whole range of sediment values. For the error maps, the same number of classes (twelve) as used in the sediment maps was employed, so that the level of detail would remain consistent between them. If one map had far fewer class breaks, it might suggest to the map reader that the data quality or precision was different; since both maps derive from the same kriging interpolation, visually communicating the same level of accuracy is important. Since the range of values for the error maps was more narrow than the sediment maps (about 1,000 – 7,000 T/km2/yr) different classification and colour schema were used. An equal interval scale (every 500 T/km 2 /yr) defined the class breaks, while a single colour progression, from white to red, suggested decreased reliability (since red has negative connotations). To provide reference information to the map reader, the river network of the Yellow River was overlaid on both the sediment and error surfaces. The course of the main stem is 68 particularly relevant to the later discussion of patterns in both maps, hence why the main stem is thicker and bolder than the tributaries. Data station locations were also added, as these are particularly useful in analyzing the error surface maps. Errors values increase significantly as one moves farther and farther away from a data point. As noted in Chapter 2, the desert region was also highlighted, as it is assumed that minimal sediment flows in or out of that region. Finally, the city of Beijing was added to help users orient themselves to the maps’ location. We begin by discussing the ordinary kriging maps (figures 3.12 and 3.13). Figure 3.12 Sediment map of the Yellow River: ordinary Gaussian kriging. 69 Figure 3.13 Error surface map of the Yellow River: ordinary Gaussian kriging. Two regions of the sediment map (figure 3.12) immediately stand out. A central, north- south strip shows ks values in the hundreds of thousands of tons per km 2 per year. Referring to figure 3.1 at the very beginning of this chapter, we see that this high-yield area corresponds roughly to sub-basin 7. Indeed, this agrees with the results of Hassan et al. (2008), who, in their sediment map (figure 4), also found this region to have very high sediment yield. The other region of the Yellow River with very high sediment yield corresponds to sub-basin 2, in the western part of the basin. The shape of both these high-yield region generally matches the river network. In particular, the highest values of the central region (>100,000 T/km 2 /yr) can be seen to diminish moving upstream, first northwards, to moderately high values (10,000-50,000 T/km 2 /yr); then, moving west and then south, as the predicted ks values taper to less than 5,000 T/km 2 /yr. This makes sense, as this region is the Loess Plateau, and it is generally agreed (e.g., Long & Chien, 1986; Yu & Xan, 2005; Hassan et al., 2008) that these highly erodible sediments are a major factor behind why the Yellow River has such extremely high sediment yield values over time. 70 With the exception of the high predicted sediment yield values in sub-basin 2, much of the upper reaches show very low yield values. Since the river in this region is confined, as it incises into bedrock, low values are to be expected. The lower reaches similarly have lower yield values—below 2,500 T/km2/yr, except for a small region about 500 km from the river mouth, where values spike to about 10,000 T/km 2 /yr. The error surface map is relatively uniform, which indicates that there is not systematic skewing of the data based on the semi-variogram properties. Error values are lowest near data points, ranging from 1,000 to 2,000 T/km 2 /yr for most of the basin. The upper- and lower-most regions yield slightly higher error values, in the range of 3,000 to 4,000 T/km 2 /yr. Only in regions far removed from any data points do we see the highest error values, greater than 5,000 T/km 2 /yr, at the very edge of the interpolated error surface. As noted in the previous section (table 3.3), the overall fit of this semi-variogram is good: RMS standard error of 0.98283, the best of all the algorithms presented in this thesis. The relatively low error values present in most of figure 3.13 further indicate the reliability of this data, which was ranked in the Data Quality Report as 11/12. Despite these encouraging indications from the error surface map (figure 3.13) that the sediment yield map (3.11) is reliable, there are some warning signs—certain regions that are less trustworthy than others. The southernmost part of the Yellow River (roughly equivalent to sub- basin 5) shows very low interpolated sediment yield values (<500 T/km 2 /yr), yet the error values are at least double this (between 1,000 and 2,000 T/km 2 /yr). Some areas of the high Tibetan plateau are similarly problematic; this is likely because there are only a few data station in this region. Areas outside of sub-basin 2 (where sediment yield values are more than 10,000 T/km 2 /yr and the error values are less than about 3,000 T/km 2 /yr) have very low sediment yield, 71 and yet moderate levels of error (errors as high as 4,000 T/km 2 /yr and sediment yield values mostly less than 1,000 T/km 2 /yr). These patterns suggest that the regions of figure 3.12 with the highest sediment yield values are the most reliable, as the errors in these regions are several orders of magnitude smaller. Since the error values are rarely above 4,000 T/km 2 /yr in figure 3.13, we can consider any prediction values above that threshold to be reliable. The magnitude of the error values in figure 3.13 is a reflection of the data quality of the Yellow River dataset, and of the specific kriging algorithm used. In regions of moderate sediment yield (5,000 – 25,000 T/km2/yr), reliability is still high, but in regions below 5,000 T/km 2 /yr (and certainly below 1,000 T/km 2 /yr) error values are often equal to (if not greater than) the interpolated sediment yield values. This is because error values are most strongly influenced by proximity to a data station, not data station value (Armstrong, 1998). In other words, high sediment yield value does not necessarily mean high error value. So for low- or high-yield points close to a data station, error will likely be smaller than the sediment value; this means the prediction is reliable. For low-yield regions far away from a data station, the error will be higher, because of the distance factor; if the error exceeds the predicted value, the predicted value is not reliable. High-yield regions, on the other hand, may still be reliable even outside of the immediate vicinity of the data points, as in the increased error values are less than the increased sediment value. 72 Figure 3.14 Sediment map of the Yellow River: universal Gaussian kriging. Figure 3.15 Error surface map of the Yellow River: universal Gaussian kriging. Why is this the case? First, in looking back at our interpolation methods, regional trends (in the form of sub-basin level sediment relations) and local conditions (from the back- calculation) both play a role in generating the sediment map. Thus, we would expect the interpolated kriging surface to vary quite significantly from sub-basin to sub-basin, as the 73 physiography, climate, land use, etc. of the Yellow River varies. The error values, however, are more dependent on the statistical algorithm and semi-variogram chosen. As we saw from the error statistics in table 3.3, the ordinary Gaussian semi-variogram was quite consistent. Looking at the graph of this semi-variogram (figure 3.9) note that the semi-variogram’s curve matches the data quite well close to a data point, but less well as distance from the data point increases. This explains the smooth, consistent error surface of figure 3.13. Clearly, the distance from a data point has a stronger effect on the error value than the underlying sediment yield value, or we would expect much more variation in the error surface. Figure 3.14 show the universal Gaussian prediction interpolation for the Yellow River. Its overall spatial distribution is very similar to figure 3.12, which showed the ordinary Gaussian interpolation. The central region of the Loess Plateau still has the highest sediment yield, with values tapering upstream, to the edge of the Tibetan Plateau. Another peak in sediment yield is in sub-basin 2, although in figure 3.14 this high-yield region is somewhat larger than it was in figure 3.12. The Tibetan Plateau is another region that diverges slightly from the ordinary kriging interpolation. In figure 3.12, moderate sediment yield levels occurred in the west part of this region, with very low values elsewhere. In figure 3.14, the situation is reversed, with values of 5,000-10,000 T/km 2 /yr in the east of the region, and low values everywhere else (<1,000 T/km 2 /yr). The most significant difference between these two maps lies just east of the main, high-yield central region. A small tributary (sub-basin 8) runs parallel to the main stem here, eventually finding a confluence near the southern edge of the basin. Ideally, a kriging interpolation would recognize a sharp distinction between these areas, as the main stem feeds off Loess material, and thus has much higher sediment yield values than the tributary. Indeed, on 74 the ordinary Gaussian sediment map, there is a relatively sharp break that runs roughly parallel to both the main stem and the tributary. In the universal Gaussian sediment map, however, this boundary is more blurry, and moderately high sediment yield values can be found in the tributary’s sub-basin. Figure 3.15 shows the error surface of the universal Gaussian kriging. It is nearly identical to figure 3.13, the error surface of the ordinary Gaussian kriging. This corroborates the previous discussion: most of the error values result from proximity to data points, rather than inconsistencies or error within the source data, or problems and biases from the semi-variograms. Some slight changes to the error surface in figure 3.15 can be seen in the central region of the Loess Plateau, just south of the region of very high sediment yield. The error values seem to be a bit lower (maybe by around 500 Tons/km2/year) than they were in figure 3.12. Overall, the sediment yield maps and error surfaces from these two kriging methods are very similar. There is a slight argument to be made in favour of the ordinary Gaussian semi- variogram, as its error statistics (table 3.3) were slightly better than the universal semi- variogram, and the ordinary Gaussian sediment yield map showed watershed boundaries a bit more clearly. Nevertheless, both methods show strong results, and demonstrate that high quality data from the Global Sediment Database can be used to create meaningful maps. 3.4 Comparing Kriging Results of Other Rivers While the ordinary Gaussian semi-variogram was arguably the best kriging method for the Yellow River, other datasets will likely require a different algorithm. Using the same kind of error statistics employed in section 3.2.3, data from the Yangtze River in China and from the 75 Danube River in Romania will be briefly examined (data sources Hassan et al., 2010 and Radoãne & Radoãne, 2005, respectively). The Yangtze River, being immediately adjacent to the Yellow River, is a natural choice for comparison. The physiography of the basins is similar, although the Yangtze lacks the copious Loess material found in the Yellow River Basin. There are more data stations in the Yangtze as compared to the Yellow River (389 vs. 161), and they are more evenly distributed throughout the basin. The westernmost part of the watershed (the Tibetan Plateau), however, has fewer points overall, just as in the Yellow River’s case. The Yangtze River scored just as high on the data quality report as the Yellow River, with 4/4 for data station location, 4/4 for temporal resolution, and 3/4 for point density (total: 11/12). Detailed sediment analyses of the Yangtze can be found in Hassan et al. (2010 and 2011). Table 3.6 Yangtze River ordinary kriging error values. Error type* Circular Spherical Exponential Gaussian Linear Mean Error 294.3 287.1 432.0 510.8 461.9 RMS Error 527.4 529.4 546.0 523.2 539.2 Mean Std. Error 0.0026 0.0028 0.0055 0.0061 0.0057 RMS Std. Error 1.055 1.071 1.142 1.036 1.1089 Avg. Std. Error 503.5 499.7 491.5 511.1 495.2 Avg. Std. – RMS Error 23.9 29.7 54.5 12.1 44.0 * RMS = root mean square; Std. = standard; Avg. = average. Bolded values are the best fit. Both ordinary and universal kriging semi-variograms are considered for the Yangtze River. The statistical measures of error presented here are the same ones used in the analysis of the Yellow River. As before, the values with the best fit are bolded. For the Yangtze River, the RMS error, RMS standard error, and average standard error were most ideal when using the Gaussian semi-variogram. This mirrors the ordinary kriging of the Yellow River. The mean standard error is lowest with the circular semi-variogram, while the mean error is lowest with the 76 spherical semi-variogram. A similar situation is found in table 3.7, which presents the results of the universal kriging semi-variograms. The Gaussian again is the best for the RMS error, RMS standard error, and average standard error. This time, however, it is the linear semi-variogram with the lowest mean standard error, and the exponential semi-variogram with the lowest mean error. Table 3.7 Yangtze River universal kriging error values. Error type* Circular Spherical Exponential Gaussian Linear Mean Error 485.2 484.0 460.3 830.0 458.8 RMS Error 528.7 531.4 546.5 525.6 539.9 Mean Std. Error 0.0066 0.0066 0.0063 0.0125 0.0060 RMS Std. Error 1.057 1.073 1.1419 1.042 1.1095 Avg. Std. Error 503.4 499.4 491.4 510.9 495.3 Avg. Std. – RMS Error 25.3 32 55.1 14.7 44.6 * RMS = root mean square; Std. = standard; Avg. = average. Bolded values are the best fit. The similar statistical results from the Yellow and Yangtze Rivers suggest that something about the spatial pattern of sediment yield is similar in both watersheds. While the distribution and number of data points is quite different, the Gaussian semi-variogram shows the lowest overall error values in both cases, for universal and ordinary kriging. From an interpolation point of view, the similarities appear to be stronger than the differences between these two basins. The Romanian data from the Danube is quite different from the Chinese data from the Yellow and Yangtze Rivers. Given how Romania and China are in such different parts of the world, with different climate, geology, land use, and so on, such variation in the kriging algorithms is expected. The exponential semi-variogram is the clear winner for ordinary kriging (table 3.8), while the universal semi-variograms are divided between the linear and circular options (table 3.9). 77 Table 3.8 Danube River ordinary kriging error values. Error type* Circular Spherical Exponential Gaussian Linear Mean Error 811.6 773.0 779.9 798.7 810.4 RMS Error 7,227.0 7,210.0 7,087.5 7,285.7 7,198.5 Mean Std. Error 0.0285 0.0209 0.0188 0.0279 0.0211 RMS Std. Error 1.259 1.254 1.226 1.275 1.252 Avg. Std. Error 5,500.0 5,507.2 5,552.2 5,473.5 5,510.3 Avg. Std. – RMS Error 1,727.0 1,702.8 1,535.3 1,812.2 1,688.2 Table 3.9 Danube River universal kriging error values. Error type* Circular Spherical Exponential Gaussian Linear Mean Error 823.0 921.4 866.7 812.2 840.9 RMS Error 7,326.0 7,338.2 7,503.7 7,361.5 7,313.5 Mean Std. Error 0.0171 0.0235 0.0230 0.0267 0.0225 RMS Std. Error 1.238 1.377 1.268 1.278 1.369 Avg. Std. Error 5,488.9 5,496.5 5,545.5 5,462.8 5,500.3 Avg. Std. – RMS Error 1,837.1 1,841.7 1,958.2 1,898.7 1,813.2 * RMS = root mean square; Std. = standard; Avg. = average. Bolded values are the best fit. Two conclusions are suggested by the Danube kriging examples. First, there is a different underlying spatial geometry in the Danube than there was in the Yellow or Yangtze Rivers. Second, no single semi-variogram will always be most ideal; it is important to investigate the statistical measures of error in each case to assess which kriging model is most appropriate for a given data set. To conclude our exploration of these alternate datasets, consider the flowchart in figure 3.16, which summarizes how to select the best kriging semi-variogram. It continues from the flowchart in figure 3.3, which reviewed how to prepare sediment data and spatial data for use in kriging interpolation. 78 Figure 3.16 Flow chart detailing selection of best kriging semi-variogram. 79 First, it must be noted again that the GIS details in this chapter are based on using ESRI’s ArcMap 10.0 and that it is assumed that similar procedure can be followed with other software. The directions in the semi-variogram selection flowchart (figure 3.16) are thus somewhat generalized. The two key points of these flow charts are the ‘checks’. Once the statistical measures of error for each semi-variogram have been examined (e.g., like tables 3.6 to 3.9 for the Yangtze and Danube Rivers, respectively), how well does each semi-variogram perform? Is one algorithm clearly superior to the others? If so, it is safe to proceed with checking the graph of that semi-variogram. If more than one semi-variogram seems possible, or if they all seem equally inaccurate, the source data should be double-checked, and perhaps the kriging parameters should be modified. The second ‘check’ in figure 3.16 is to investigate how well the overall shape of the semi-variogram matches the data, using graphs like those found in figures 3.8 to 3.11. Within ArcMap, these graphs appear during the kriging selection process within the Geostatistical Wizard, and the graphs—with accompanying statistics, and regression lines—can be exported for more detailed consideration. If the match between the best semi-variogram and its data is not suitable, the 2 nd -best semi-variogram can be considered (just as the best two semi-variograms were considered for the Yellow River in this Chapter). If this still does not produce a satisfactory result, especially near the origin (0,0) of the graph, then the dataset as a whole should be re-evaluated and possibly re-analyzed to ensure that the quality of the data is sufficiently high to proceed with kriging interpolation. One conclusion for this Chapter is that the maps created for the Yellow River using the ordinary Gaussian semi-variograms are the best representations we could make with the data and 80 statistical tools available to us. More importantly, even with the brief examples of the Yangtze and Danube Rivers, it is clear that ‘one size fits all’ is not an appropriate analytical methodology when considering spatial interpolation. The true goal of this Chapter was to provide a methodological framework, so that potential users of the Global Sediment Database can make well-reasoned decisions about the kinds of GIS analyses they might perform. 81 Chapter 4 Conclusions This thesis has demonstrated how the Global Sediment Database can be created, and how it can be used to explore the spatial pattern of sediment yield in the landscape. The remaining question is how to make the Database and the Data Quality Report available to other researchers? The ongoing work of updating the Database and evaluating data quality must be done in a consistent way. Further, the Database itself needs to be tested in peer-reviewed literature. Material from this thesis, including data from the Database, will hopefully be used as the basis for further work. Once our framework for collecting, processing, and documenting data has been vetted by the scientific community, we will make the Database and the Data Quality Report available to all. At that point, we hope that the data analysis framework presented in Chapters 2 and 3 can serve as a starting point for researchers who wish to use the Database. The step-by-step discussion of how to process, evaluate, analyze and interpolate sediment data can be applied to any part of the Database, or indeed to any sediment yield dataset. By linking methodological decisions to underlying data quality we have shown how a systematic analysis of spatial and attribute data can produce a result with a known level of confidence. This attention to data uncertainty is expressed in the kriging error maps that accompany our analysis, and in the statistical exploration of the various kriging methods. Still, it is clear that no one method is ideal, so comparison between datasets is going to be essential in developing these methods. The goals of this thesis were to build the Global Sediment Database and to present a method for choosing the best tool for spatially analyzing sediment yield data. Neither of these goals has an endpoint: our aim is for the Database to continue to grow over time, and for even 82 more interpolation methods to be compared and contrasted for use with sediment data. Progress on these goals has at least begun, and that is an accomplishment. A third goal has also emerged from this discussion of finding the best interpolation method. GIS has often been described as a ‘black box’ —data are thrown in, and analyses are spat out. The analysis of various kriging methods in this thesis has helped to crack open that black box, and illustrate that there is much to be gained in exploration of the nuances of your data, and the many, many analysis options available. As we progress into the 21 st century and more parts of our world become geo-tagged, linked to GPS, or uploaded to the ‘cloud’, it is important that we, as geographers, maintain our critical gaze on the new media and tools that are being created. The scientific, social, and environmental implications of spatial data will change as large datasets like the Global Sediment Database become the norm. Issues of data quality and reliably will have to be considered, and the urge to produce results quickly—by using the black box of GIS, for instance—must be resisted. Justification and documentation must be provided at every step. If patience, diligence, and attention to detail can be incorporated into the analysis of these new data frontiers, then the knowledge of these new technologies can benefit everyone. 83 References Anderson, M.G. & McDonnell, J.J., Eds. 2005. Encyclopedia of hydrological sciences. John Wiley & Sons, West Sussex, UK. Armstrong, M. 1998. Basic Linear Geostatistics. Springer, New York. Bivand, R.S., Edzer, J., Gómez-Rubio, P., & Gómez-Rubio, V. 2008. Applied Spatial Data Analysis with R. Springer Science+Business Media, LLC. Chorley, J., Schumms, A., and Sugdend, E. 1984. Geomorphology. Methuen, London, UK, cited in Church, M.; Kellerhals, R. & Day, T.J. 1989. Regional clastic sediment yield in British Columbia. Canadian Journal of Earth Science, 26, 31-45. Church, M.A. 2006. Bed material transport and the morphology of alluvial river channels. Annual Review of Earth and Planetary Sciences, 34 (325-354). Church, M., Ham, D., Hassan, M., & Slaymaker, O. 1999. Fluvial clastic sediment yield in Canada: Scaled analysis. Canadian Journal of Earth Sciences, 36(8). Church, M.; Kellerhals, R. & Day, T.J. 1989. Regional clastic sediment yield in British Columbia. Canadian Journal of Earth Science, 26, 31-45. Church, M. & Slaymaker, O. 1989. Disequilibrium of Holocene sediment yield in glaciated British Columbia. Nature, 337. Cooke, R.U. & Doornkamp. J.C. 1990. Geomorphology in environmental management. Clarendon Press, Oxford. Dedkov, A.P. & Mozzherin, V.T. 1992. Erosion and sediment yield in mountain areas of the world. In Erosion, Debris Flows, and Environment in Mountain Regions. Walling D.E., Davies, T.R. & Hasholt, B., Eds. IAHS Publication No. 209, IAHS Press: Wallingford, UK. de Vente, J., Poesen, J., Arabkhedri, M., & Verstraeten, G. 2007. The sediment delivery problem revisited. Progress in Physical Geography, 31(2), 155-178. Environmental Systems Resources, Inc (ESRI). 2011. ArcGIS 10.0 desktop help: Geostatisical Analyst. http://help.arcgis.com/en/arcgisdesktop/10.0/help/index#//003100000004000000.htm. Accessed July 21 st , 2012. Fournier, F. 1960. Climat et Erosion, cited in Anderson, M.G. & McDonnell, J.J. 2005. Encyclopedia of hydrological sciences. John Wiley & Sons, West Sussex, UK, and also cited in Goude, A. 1995. The changing earth. Blackwell Publishers, Oxford. Goudie, A. 1995. The changing earth. Blackwell Publishers, Oxford. 84 Hassan, M.A. Church, M, Xu, J.X. & Yan, Y.X. 2008. Spatial and temporal variation of sediment yield in the landscape: Example of Huanghe (Yellow River). Geophysical Research Letters, 35(6). Hassan, M.A., Church, M., Yan, Y., & Slaymaker, O. 2010. Spatial and temporal variation of in-reach suspended sediment dynamics along the mainstem of Changjiang (Yangtze River), China. Water Resources Research, 46. Hassan, M.A., Church, M., Yan, Y., Slaymaker, O., & Xu, J. 2011. Suspended sediment balance for the mainstem of the Changjiang (Yangtze River) in the period 1964-1985. Hydrological Processes, 2011. Jansson, M.B. 1988. A global survey of sediment yield. Geografiska Annaler, Series A, 70. Leeks, G. 2005. Measuring sediment loads, yields, and source tracking. In Encyclopedia of hydrological sciences, Anderson, M.G. & McDonnell, J.J., eds. John Wiley & Sons, West Sussex, UK. Long, Y., and N. Chien. 1986. Erosion and transportation of sediment in the Yellow River basin, International Journal of Sediment Research, 1, 1– 38. Lvovich, M.I., Karasik, G.Y., Btrstseva, N.L., Medvedeva, G.P., & Maleshko, A.V. 1991. Contemporary Intensity of the World Land Intracontinental Erosion. USSR Academy of Sciences, Moscow. Meade, R.H., Yuzyk, T.R., and Day, T.J. 1990 . Movement and storage of sediment in rivers of the United States and Canada. The Geology of North America: Surface Water Hydrology. Vol 1, Geological Society of America. Boulder, CO. Milliman, J.D. & Farnswoth, K.L. 2011. River discharge to the coastal ocean. Cambridge University Press, Cambridge, UK. Milliman, J.H. & Meade, R.H. 1983. World-wide delivery of river sediment to the oceans. Journal of Geology, 91(1). Montgomery, D.R. 2007. Dirt: The erosion of civilizations. University of California Press, Los Angeles. Negreiros, J., Painho, M., Aguilar, F., & Aguilar, M. 2010. Geographical Information Systems principles of ordinary kriging interpolator. Journal of Applied Sciences, 10 (11). Nearing, M., Renard, K., & Nichols, M. 2005. Erosion predication and modelling. In Encyclopedia of hydrological sciences, Anderson, M.G. & McDonnell, J.J., Eds. John Wiley & Sons, West Sussex, UK. 85 Radoãne, M. and Radoãne, N. 2005: Dams, sediment sources and reservoir silting in Romania. Geomorphology 71, 112–25. Simon, A. & Klimetz, L. 2008. Magnitude, frequency, and duration relations for suspended sediment in stable ("reference") southeastern streams. Journal of the American water resources association, 44(5). Slaymaker, O. 2006. Towards the identification of scaling relations in drainage basin sediment budgets. Geomorphology, 80. Strakhov, A.N. 1967. Principle of Lithogensis (3 vols.), cited in Goude, A. 1995. The changing earth. Blackwell Publishers, Oxford. Swarzenski, P.W. & Campbell, P.L. 2005. On the worldwide riverine transport of sediment – Associated contaminants to the ocean. In Encyclopedia of Hydrological Science, Anderson, M.G. & McDonnell, J.J., eds. John Wiley & Sons, West Sussex, UK. Vanmaercke, M., Poesen, J., Verstraeten, G., De Vente, J., & Ocakoglu, F. 2011. Sediment yield in Europe: Spatial patterns and scale dependency. Geomorphology, 130. Walling, D.E. 1977. Assessing the accuracy of suspended sediment rating curves for a small basin. Water Resources Research, 13. Walling, D.E. 1983. The sediment delivery problem. Journal of Hydrology, 65. Waling, D.E. 2005. Sediment yields and sediment budgets, In Encyclopedia of hydrological sciences, Anderson, M.G. & McDonnell, J.J., eds. John Wiley & Sons, West Sussex, UK. Walling, D.E. & Webb, B.W. 1981. The reliability of suspended sediment load data. Erosion and Sediment Transport Measurement, IAHS Publication No. 133, IAHS Press: Wallingford, UK. Walling, D.E. & Webb, B.W. 1983. Patterns of sediment yield: a global overview. In Background to paleohydrology, Gregory K.J., ed. Wiley: Chichester, UK. Wang, S.J., Hassan, M.A., & Xie, X.P. 2006. Relationship between suspended sediment load, channel geometry and land area increment in the Yellow River Delta. Catena, 65(3). Xu, J., and Y. Yan. 2005. Scale effects on specific sediment yield in the Yellow River basin and geomorphological explanation, Journal of Hydrology, 307, 219– 232. Yair, A. & Enzel, Y. 1987. The relationship between annual rainfall and sediment yield in arid and semi-arid areas: the case of Northern Negev, Israel, cited in Goude, A. 1995. The Changing Earth. Blackwell Publishers, Oxford.