A structure generator for modelling the initial sediment distribution of an artificial hydrologic catchment

Artificially-created hydrological catchments are characterised by sediment structures from technological construction processes that can potentially be important for modelling of flow and transport and for understanding initial soil and ecosystem development. The subsurface spatial structures of such catchments have not yet been sufficiently explored and described. Our objective was to develop a structure generator programme for modelling the 3-D spatial distribution patterns of dumped sediments depending on the technical earth-moving and deposition processes. We are focussing in a first step on integrating sediment dumping, particle size, and bulk density modification processes on the catchment scale. For the model development, the artificially-constructed hydrological catchment “Chicken Creek” located in Lower Lusatia, Germany, served as an example. The structure generator describes 3-D technological sediment distributions at two scales: (i) for a 2D-vertical cross-section, texture and bulk density distributions are generated within individual spoil cones that result from mass dumping, particle segregation, and compaction and (ii) for the whole catchment, the spoil cones are horizontally arranged along trajectories of mass dumping controlled by the belt stacker-machine relative to the catchment’s clay layer topography. The generated 3-D texture and bulk density distributions are interpolated and visualised as a gridded 3-D-volume body using 3-D computer-aided design software. The generated subsurface sediment distribution for the Correspondence to: T. Maurer (maurer@tu-cottbus.de) Chicken Creek catchment was found to correspond to observed patterns already without calibration. Spatial aggregation and interpolation in the gridded volume body modified the generated distributions towards more uniform (unimodal) distributions and lower values of the standard deviations. The modelling approach is generally applicable to all situations where large masses of unconsolidated sediment are moved and dumped thereby allowing generation of basic soil structures and patterns of hydrological systems.


Introduction
Hydrologic catchments can be defined as discrete geo-bodies of solid rocks or sediments with a given 3-D geometry.Ideally, a catchment is delineated at the bottom by a layer of rocks or sediments with low permeability that restricts water percolation, forming the subsurface catchment (Brutsaert, 2005).At the surface, the catchment is characterised by topography and land use conditions.Problems regarding hydrologic modelling are often related to uncertainties in subsurface structures and hydraulic properties and with respect to the spatial extent and geometry of the bottom and lateral boundaries.Furthermore, discrepancies in lateral delineation between soil surface and bottom layer can occur.For better identification of catchment boundaries, hydrogeological data (i.e.spring discharge, soil and geological surveys, tracer tests, and isotopes) have been combined with modelling (e.g.Benischke et al., 2010), and digital elevation models (DEMs) have been used (Hammond and Han, 2006).Also, geophysical approaches are increasingly being Published by Copernicus Publications on behalf of the European Geosciences Union.
T. Maurer et al.: Initial sediment distribution of an artificial hydrologic catchment applied for the non-invasive exploration of subsurface structures (Reynolds, 2005).Electromagnetic methods, namely ground penetrating radar, are widely used for investigating near-surface hydrogeological problems and structures since the beginning of the 90s (e.g.Beres and Haeni, 1991).Moreover, electrical resistance tomography provides constant hydrological monitoring of structure-process interactions in aquifers (Kosinski and Kelly, 1981;Revil et al., 1998).The integration of such methods into spatial watershed modelling is, however, still not well established (Robinson et al., 2008).
The particular processes of catchment formation are important to understand the distribution and characteristics of spatial structures (e.g.along layers or facies of fluvial deposits).For instance, alluvial processes produce sedimentary structures while the structure of aquifers in karst regions depends on the distribution of fractures (e.g.Siemers and Dreybrodt, 1998).Structural differences between the alluvial and karstic domains of the Orleans valley aquifer were found to result in a "dynamically confined system" with reduced hydrochemical exchange between both domains (Le Borgne et al., 2005).
The specific catchment formation processes determine the relevant spatial structures, which in turn determine a catchment's hydraulic properties.Modelling their spatial distribution in dependency of the structure of the solid phase was subject of recent investigations.Gascuel-Odoux et al. (2010), for example, tested the impact of spatial heterogeneity of hydraulic conductivity on hydrologic model performance (Hill-Vi, Weiler and McDonnell, 2004), showing that including lateral variation in hydraulic conductivity seems to be as important as variation in soil depth for correctly representing the fluxes and dynamics in a hillslope aquifer.In a similar study, Gauthier et al. (2009) showed that only sufficiently complex scenarios of spatial heterogeneity can adequately reproduce catchment features and response variables.An evaluation of the applicability of pedotransfer functions for the spatial modelling of soil hydraulic properties (Stumpp et al., 2009) unanimously found that prediction quality increases when assuming complex spatial distributions of soil hydraulic properties.Sciuto and Diekkrueger (2010) quantified the impact of spatial variability of soil hydraulic properties and the effects of spatial discretisation for an experimental catchment.They found that aggregating soil hydraulic properties results in lower uncertainties than does using a coarser discretisation.These results stress the importance of verifying and conditioning spatial catchment models with measured data to achieve a sufficiently complex description of spatial heterogeneity.
A classical approach in geosciences for the description of structural heterogeneity is to utilise geostatistical methods.The spatial interpolation between measurement points is generally not suited for depicting sharply bounded structures like different sediment facies (Michael et al., 2010).More sophisticated methods like multi-point geostatistical techniques excel indeed at respecting local data and allow more realistic representations of spatial patterns, but are still dependent on the input of training images that already purport the spatial heterogeneity (Hu and Chugunova, 2008).Given the shortcomings of spatial modelling based on using "pure" geostatistics (Koltermann and Gorelick, 1996), alternative approaches have been developed in recent years: Grunwald et al. (2000) implemented the Virtual Reality Modelling Language in combination with soil data, topographical attributes, and kriging to create 3-D representations of soil landscapes; Zappa et al. (2006) conditioned 3-D geostatistical simulations using "model blocks" derived from soil profiles that represent facies associations to model a glaciofluvial aquifer; Moreton et al. (2002) used a physical model of the subsurface depositional stratigraphy of a braided river system to feed object-based digital spatial reservoir models, while Teles et al. (2001) used a multi-agent concept in combination with simple construction rules based on literature and observations to reproduce the structure of an alluvial plain.Sech et al. (2009) developed a surface-based spatial modelling approach that enables explicit representation of heterogeneity across a hierarchy of length scales.In this context, Murray et al. (2009) propose introducing the development of a hierarchical suite of structural models to achieve predictive power over the vast range of scales of earth-surface phenomena.
Catchment and reservoir modelling approaches distinguish between object-based methods that distribute predefined "geobodies" within a model domain and process-based techniques that describe the physical processes rather than the existing structure (Michael et al., 2010).For instance, Pyrcz et al. (2009) developed an object-based model that constructs stochastic pseudo-process-based fluvial models.Gross and Small (1998) formulated a process model to simulate the development of river facies, and Miller et al. (2008) used a process-modelling approach to generate geologically realistic structures of turbidite fans.Their process-based spatial model was conditioned with measured data using geostatistical techniques (Michael et al., 2010).Such structuregenerating methods produce only stochastically possible realizations representing the reservoir.Since the actual structures can never be reproduced adequately (de Marsily et al., 2005), consideration of uncertainties of the predicted spatial distributions is important.
Any process-based "structure generator" model that is able to reproduce the key structural elements requires basic knowledge of formation processes and of the arrangement of structural elements.Such information can most easily be collected for artificial catchments, which have been comprehensively planned and constructed such that size, geometry, sediment composition, and boundaries are much better defined than natural systems.An artificial catchment thus qualifies as a large scale research tool for investigating hydrological problems and ecosystem development (cf.Hooper, 2001).Up to now, only a few of these open air laboratories were implemented: Gu and Freer (1995) describe experiments on  2004) used mixtures of mineral material and peat placed over glacial mineral soil above a sealing saline-sodic shale surface for testing changes in hydraulic conductivities; and Nicolau (2002) investigated runoff generation and routing on artificial slopes derived from open cast mining reclamation.Depending on the size of the catchment construction, more-or-less "homogeneous" sediment distributions can be realised (Kendall et al., 2001).
Nevertheless, construction process-specific sediment heterogeneities may still have a considerable impact on hydrological processes.Knappe et al. (2004), for instance, found relatively high variability in seepage rates, which are probably reflecting the heterogeneity of forest-reclaimed overburden dumps in the central German lignite mining district.Tracer experiments by Hangen et al. (2005) on similar overburden spoil sediments in Lower Lusatia, Germany, indicated that flow processes were influenced by inclined dumping structures and lignite fragments (cf.Gerke, 2006).For modelling such overburden spoil piles, the spatial heterogeneity of sediment texture and bulk density was considered by taking technical dumping and mixing processes into account (Buczko and Gerke, 2005a).Experiments in 2-D-slab cells by Lebron and Robinson (2003) suggested that internal layering of different grain sizes can occur due to avalanching processes.It remains challenging to transfer such sediment structures into hydraulic and transport properties such that the predicted flow paths reflect the internal structure of the artificial system not only for cross-sections (e.g.Buczko and Gerke, 2005b) but for catchments.
The objective of this study was to develop a structure generator programme that models the 3-D spatial distribution of dumped sediments by considering the technical earth-moving processes and the composition of potential parent materials at the outcrop site.The purposes of our processbased model were to create the basis for a tool that (i) generates 3-D-realisations of spatial sediment distributions in constructed catchments, (ii) provides basic 3-D information for testing and comparing hydrological models (Holländer et al., 2009) and for studying structure-and-process interdependencies, and (iii) creates 3-D catchments for simulating the dynamic initial structural development using the functionality of an advanced geological modelling software.In this first step, the model will be focussing on the basic method for integrating sediment dumping, particle size, and density modification processes on the catchment scale.Other subsequently occurring catchment construction processes (i.e.compaction by bulldozing and filling of remaining volumes) will be included later.
First, we describe the modelling of the textural and density distribution in individual 2-D-vertical cross-sections (i.e.spoil cones) that result from the dumping of sediments.Then, these individual cones are horizontally distributed in the catchment along trajectories following observed spoil ridges.Eventually, by combination and spatial interpolation, the distribution of sediment properties in a representative part of the 3-D catchment volume is compiled and consequences for the final structure of the generated artificial catchment are discussed.

The artificial catchment "Chicken Creek"
The artificial hydrological catchment 'Chicken Creek' (i.e. in German: "Hühnerwasser") is located in the post-mining recultivation area of the open cast lignite mine "Welzow Süd", about 20 km south of the city of Cottbus in the Lower Lusatian lignite mining area (Fig. 1).The 6-ha catchment was constructed by Vattenfall Europe Mining (VEM) to serve as headwaters for the restored streambed of the former "Chicken Creek".It has the form of a NW-SE running, almost rectangular, slightly inclined plane with a lengthwise extension of about 450 m and a transversal extension of about 160 m.The terrain inclines from NW to SE with an average of 3.5 • .The size of the catchment is 5.9 ha at the soil surface and 6.1 ha at the bottom clay layer; thus the subsurface catchment is slightly larger than that at the surface.The catchment is draining in southeastern direction into a pond.The construction of the catchment was finished in October 2005.Since then, the soil-geosystem developed without human interference (Gerwin et al., 2009a).
The catchment was constructed using large-scale stacker technology.The material was excavated from the undisturbed excavation site of the open cast lignite mine and was delivered over a distance of several kilometres by conveyor belts.The catchment foundation underneath was dumped by stackers above the conveyor bridge spoil to form the base structure (Fig. 2).Both conveyor bridge spoil and stacker base spoil consist predominantly of coarse-textured, tertiary overburden material.Above the base structure, a 1-3 m clay base liner of low permeability (Kendzia et al., 2008) was applied, acting as the catchment's aquiclude .It consists of the quaternary "Flaschenton" ("bottle clay"), which is fetched as a by-product of mining operations.The clay base liner forms an elongated, bowl-like subsurface structure, thus purporting groundwater flow towards the pond, respectively the single catchment outlet (Gerwin et al., 2011).
North of the pond, a subsurface clay dam, running roughly in E-W direction, was constructed (Fig. 3).This clay dam has two functions: (i) to serve as a stabilisation barrier to prevent the sandy sediments to slide downhill upon the clay layer and (ii) to direct groundwater flow towards an artificial spring that discharges into the pond.North of the clay dam, quaternary material was applied using large-scale stacker technology.The uppermost covering layer consists of 1.5-3 m coarse-textured sandy overburden sediments of quaternary origin.These sandy-loamy sediments constitute the permeable catchment volume.The permeable sediment body is  explicitly delimited to the sides by the clay base liner in the form of an elevated brim, which ends about 0.5-2 m below the ground surface (Fig. 2c).The resulting surface and subsurface catchments are widely, but not totally congruent.The construction was carried out between May 2004 and October 2005.Two slightly different materials were dumped on the eastern and western part of the catchment (Fig. 3).The stacker dumping technology produced characteristic structures, i.e. spoil ridges consisting of single overlapping spoil cones (Fig. 4).This fragmentation in cones was caused by the stepwise sweep of the stacker arm.Between the eastern and western dumps remained a central trench (Figs. 3 and 4).Between April and July 2005, this depression was filled by bulldozers, scraping material from both sides into the trench, thus producing significantly different internal structures in the centre of the catchment.After bulldozing, the catchment surface was flattened to remove any artificial unevenness (September/October 2005).
During the construction phase, spatial data were collected about the setting of the clay base liner.The main data source was the photogrammetric analysis of aerial photographs carried out by the mine-surveying department of Vattenfall Europe Mining (Table 1).These data covered mainly the area around the central trench and the southern area around the pond.A complete photogrammetric coverage of the clay base liner was not possible due to the limited number of aerial photographs taken during construction.Borehole field campaigns were carried out in the months after construction to collect data about the exact position of the base liner in areas which were not covered by the photogrammetric surveys.In order to define the subsurface catchment, several borehole campaigns were carried out to determine the borders of the base liner and, in particular, its peripheral topography (Gerwin et al., 2009b;Maurer et al., 2009).Aerial photographs and corresponding DEMs from different points of time during construction ware also used to reconstruct the principal layout of construction elements and volumes (see Sect. 2.3).

Soil and sediment properties
Soil samples for texture analysis were taken between October 2005 and April 2006 in depths of 0-0.3 m and 0.3-1.0m along a 20 m × 20 m grid.Additional soil samples were taken in a 40 m × 40 m grid in depths of 0-0.5, 0.5-1.0,1.0-1.5 and 1.5-2.0 m.Altogether, 311 borehole samples were collected from raster sampling.Results showed that the material used for construction of the covering layer consists of coarse-textured sediments (Table 2) in which the sand fractions (i.e.2.0-0.063mm particle diameter) dominate with an average of about 83 % (coarse sand 13.1 %, medium sand 46.9 %, fine sand 23.6 %), and with contents of about 9 % silt (i.e.0.063-0.002mm) and 7 % clay (i.e.<0.002 mm) only for the "fine earth" particles <2 mm).The pleistocene, mainly glacigenic material (Table 3) has relatively high skeleton (i.e.>2 mm) contents of about 12 % in relation to total mass.Due to the aforementioned separate delivery of two materials for the western and eastern part of the catchment, the average texture of each side shows noticeable differences: the eastern part has a 5 % higher sand content , and accordingly lower silt (2 % less) and clay (3 % less) contents than the western part (Table 2).
During the construction phase in January 2005, additional samples were taken from 11 spoil cones in approximately 0.5 m height above the clay base liner (Gerwin et al., 2009b).These samples were used for X-ray diffraction analysis (Whittig and Allardice, 1986).Soil texture was determined after elimination of carbonates (H 2 O 2 treatment was omitted because organic matter was absent) using wet sieving and the pipette method (Gee and Bauder, 1986).The clay fraction is dominated by illite (41-60 %) and mixed-layer clay minerals (12-38 %; <50 % illite).Kaolinite (9-16 %) and vermiculite (3-4 %) were found in relatively small amounts.More soil properties can be found in Gerwin et al. (2009b).

Parent material of the outcrop at the excavation site
The parent material of the Chicken Creek sediment body was scraped off by bucket excavators located in the forefront of the Welzow Süd mining area and delivered to the stacker via conveyor belts.Texture data were obtained for most of the geological units from Vattenfall Europe Mining (Table 3).Aerial photographs show that the Chicken Creek sediments were delivered in two separated batches in September and October 2004, for the eastern and western parts.The exact location of the excavator was not recorded for these dates.

The structure generator model
The conceptual model considers the processes of parent material dumping and levelling of the catchment surface (Fig. 4).The generator considers two scales: 2-D overburden spoil cones and catchment scale.First, the base material composition is defined considering the geologic conditions at the excavation site and mixing due to conveyor belt trans-port (Fig. 4a).Then, a series of 2-D spoil cone cross-sections are constructed.For that purpose, we modified and extended the original generator (core) program (Buczko et al., 2001) to (i) include an oblique cone axis (Fig. 4b), (ii) consider adjacent spoil ridges and their impact on cone geometry (Fig. 4b), (iii) generate a given number of cross sections consecutively in the correct spatial alignments (Fig. 4c) by, (iv) transforming grid-coordinates (x, z) into a geographic coordinate system, and (v) to include a more flexible way of parameter input via Excel-sheets.Following the approach by Buczko et al. (2001), the programme calculates particle segregation at the cone flanks and compaction in the central zone.The structure generator programme was written in Visual Basic for Applications (VBA, Microsoft Corp., Redmond) to ensure a straightforward incorporation of input data stored in Excel-sheets.In a final step, data are aggregated and interpolated in a 3-D representation of the catchment (Fig. 4) using the GOCAD software (version 2009.3patch 1, Paradigm Ltd., George Town, KY).GOCAD stands for "Geological Objects Computer Aided Design" and enables interactive   3-D geologic modelling of the geometry and properties of complex subsurface objects.

Geometry of individual 2-D spoil cone cross sections
In the model, each spoil cone is virtually represented as a 2-D cross section consisting of gridded data points.The principal cone shapes and properties stored in a 2-D-cross section grid are basically governed by input parameters.Input parameters are directly derived from the recorded spatial settings or can have adjustable values.2-D-cross sections of spoil cones are arranged in sequences with an adjustable horizontal spacing along the digitised spoil ridges.Spoil ridges are digitised as curve objects consisting of segments and nodes in GO-CAD.Each node represents the position of a spoil cone cross section in the catchment model.The density of nodes on the curve object representing the spoil ridge can be adjusted in GOCAD, allowing variable horizontal distances between cross sections in the spatial model.Node coordinates are exported as ASCII files to serve as input data for the structure generator.The programme works off every 2-D cross section in the ASCII listing successively by considering the spatial framework (i.e.slope angle, spoil cone height, oblique central axis, spacing between adjacent ridges, spatial orientation of the current ridge) in the outermost loop (Fig. 5).
The basic unit with an initially (i.e.before the dumping process) homogeneous sediment composition is defined as a "layer".We assume that layers are generated when a surge of the parent material flows down the spoil cone flanks, undergoing segregation processes at the same time.Layer thickness ds m (Table 4) was defined ex ante based on (i) observations of dune avalanches (Andreotti and Bonneau, 2009) to lie in the magnitude of several centimetres, and on (ii) considerations to provide a spatial resolution of these smallest elements that is still observable after spatial aggregation.The structure generator slightly varies the texture for each layer (see below), simulating the variations in the composition of the successively dumped sediment.This results in spatial heterogeneity on the basic spoil cone scale, which becomes especially relevant when regarding high-resolution sections of the spatial model.From the layer thickness ds m , the horizontal and vertical lengths of a layer m, dx m and dz m are calculated: Values for slope angle α were defined ex-ante based on values for sand slopes taken from literature (e.g.Hoffmann et al., 2005), while spoil cone height H and back width R (i.e. the horizontal distance to the following spoil ridge) were derived based on spatial analysis in GOCAD (see Sects.2.3 and 3.1).Using these values (Table 4), the spoil cone length L can be calculated as: The initial impact point is defined as the point where the dumped material of a new spoil cone initially hits the surface.In most cases, depending on spoil cone position on the ridge and its height, this point is located on the flank of a neighbouring spoil ridge.The z-coordinate Zinit of the initial impact point is calculated first as: with β being the angle of the oblique central axis of the spoil cone cross section.The x-coordinate Xinit can then be defined as: The total number of layers l within each spoil cone is calculated as: The example in Fig. 4b has eight layers (l = 8).

Spatial alignment of cone cross sections in the catchment
Raster points in cross sections need to be assigned to geographic coordinates according to their correct alignment on the spoil ridge (Fig. 6).Since construction routine of the Vector rotation was used to align cross sections correctly on the spoil ridge.Depending on their position on a ridge, cross sections are not properly aligned towards the stacker's position.This is a consequence of vector rotation.For example, for the east side of Chicken Creek catchment, stacker positions can be either "southeast" or "northeast" relative to spoil ridges.Whether the stacker was located in the NE or SE of the catchment can not be reconstructed.Assuming a relative stacker position in the southeast, cross sections located on ridges with a "northward" component need to be mirrored or flipped, otherwise they would point in the opposite direction to the "correctly" aligned cross sections.Technical details about the alignment algorithms can be found in Appendix A.

Definition of internal structure
Spoil cone cross sections are generated using two loops, the outer loop for the spoil cone lengths, and the inner loop for the spoil cone heights.This means that cross sections are assembled column-wise from left to right.This generates a 2-D grid where, in principle, each grid point is tested (i) whether it is inside or outside of the spoil cones, (ii) whether it is on the left or right side of the central compaction zone and (iii) to which layer m it will ultimately be assigned.A grid point is located inside the spoil cone if the following conditions are true: with Zcorel and Xcorel being the relative coordinates used during grid construction.The back width R is required to calculate where the spoil cone is "cut off" by the preceding spoil cone/spoil ridge.
In our model, we assume that the more or less oblique impact of sediment dumped by a stacker also causes spoil cone growth along an oblique central axis with the angle β (Table 4).As a consequence, layer thicknesses differ for the right and left sides of this axis.Thus, in a next step, the programme checks if grid points are located to the left or to the right of the oblique central axis of the spoil cone (left side check -condition is true if): In case that the angle of the central axis, β is oblique (i.e.not equal to 0 • ), two separate solutions for the assignment of grid cells to a layer need to be applied.In case the grid point is located to the left of the oblique central axis (Eq. 10 is true) then the grid point is assigned to layer m i if two conditions are true: In case the grid point is located to the right of the oblique central axis, the two following conditions are checked: where dxr m is being defined as the extension of layer m in x-direction for layers to the right of the central axis according to: and with Zd being the distance H − Zinit, calculated using the formula:

Texture definition
In a first approach, we aim at reproducing spatial patterns of sediment heterogeneity as being observed in the post-mining area of Welzow-Süd.Here, according to an educated guess in the field, sediment properties on a spoil ridge change significantly after a distance of about 10 m-30 m.In our model, we therefore define a "basic sediment composition" P for a series of 2-D cross sections (i.e.spoil cones).The corresponding number of cross sections, N, is randomly chosen between 3 and 9. Any variation in the upper and lower boundary value for N , results in higher or lower mean lengths of spoil ridge sections (Fig. 4c) with the same basic sediment composition.The basic sediment composition, P , is determined by mixing two randomly chosen sediment types P 1 and P 2 (Table 3) according to a random ratio rnd.Each of nparticle size classes (here: skeleton, coarse, medium, fine sand and silt, clay), subscript i, in the basic sediment P i is defined by mixing P 1 i and P 2 i using the random factor rnd, which lies between 0 and 1: It is assumed that each spoil cone has a slightly different composition Psc that differs from the basic sediment composition P .Variations are, however, supposed to occur in a relatively narrow limit with the maximum value Cscmax, which is an adjustable parameter (e.g.Cscmax = 0.2 allows fluctuations of maximal 20 %, Table 4).Particle distributions are being varied at 2 scales: (i) For the spoil cone, every spoil cone particle size class Psc i * is defined using: After variation, the particle sum is uneven, so each of n particle fractions Psc i * has to be corrected to add up properly by using: (ii) For each layer l, the particle distribution is further varied.We assume that each spoil cone shows slight heterogeneities in the internal sediment composition.The initial sediment composition of each layer is varied in narrower limits as for spoil cone sediment variations using the pre-defined parameter Csclmax (Table 4): The random mixing of the often heterogeneous parent material and the application of the mixtures according to observed distances along a spoil ridge is designed to reproduce the specific patterns of stacker dumping.At this stage, the approach does not yet consider the geological configuration at the excavation site.

Segregation, compaction and bulk density distribution in individual cones
Following the establishment of one spoil cone cross section, the grid is filled layer-wise with the predefined sediment for each layer.In the model conception, the particle distribution of the material is altered by segregation processes due to gravitational deposition, ultimately resulting in inverse grading.More sophisticated theories on particle segregation processes were developed in recent years (e.g.Gray et al., 2006;Dasgupta and Priyanka, 2011); their implementation is, however, beyond the scope here because it would require considerably larger efforts for parameterization, additional data, and it would be computationally more demanding.For the calculation of particle segregation, we adopted a relatively simple empirical approach (Leibiger, 1964) that was based on observations and further developed by Buczko et al. (2001): Here, Pscl i (x,z) is the mass of the particle class i at the location with the coordinates (x,z) after segregation, H max , defined as the vertical distance from the reference elevation level H k , is calculated as follows: where dz m is the vertical extension of the current layer m.The reference elevation level H k is defined as the actual value of the relative height coordinate Zcorel.ζ i is an empirical dimensionless segregation factor gained from quaternary sediments in Lower Lusatia (Schlabendorf-Nord), ranging (Buczko et al., 2001) from −0.12 for clay particles to 0.12 for coarse sand particles and skeleton.2-D distribution of spoil bulk densities ρ b are calculated based on the pre-defined dumping height and the calculated particle size distribution from Eq. ( 21) as (Buczko et al., 2001): where ρ b0 is the initial estimated value (based on calibration data from spoil cone sampling) of the bulk density, w U and w R are weighing factors to account for the relative influences of the degree of non-uniformity and the random component on bulk density, U m is the average value and U (x,z) is the value of the degree of non-uniformity at position (x,z).Z * is a spatially uncorrelated random number between 0 and 1, while a and b are empirical factors accounting for the variation of the bulk density introduced by the grain size distribution and the random component.
According to data from spoil cone sampling, we assume values of ρ b0 = 1.2 g cm −3 for the whole spoil cone except for the central compaction area, where bulk density is calculated as a function of the changing dumping height and an elevated base value (here: 1.4 g cm −3 ) according to a linear regression equation (Matschak, 1969): where ρ b,c is the bulk density in the compaction zone in g cm −3 and H D is the pre-defined sediment drop height in metre (Table 4).
The width of the central compaction zone was defined as the input parameter Czone (Table 4), based on values from spoil cone field sampling (see Sect. 2.4).Grid points are checked if they are located inside the central compaction zone using conditional inquiries: If the condition is true, ρ b0 = ρ b,c , otherwise ρ b0 = 1.2.

Data aggregation and output
Datasets generated for hundreds of spoil cone cross sections with generic original resolutions of 1 cm 2 are far too voluminous to be handled properly with the currently available computing capabilities.Thus, prior to export to GOCAD, datasets need usually to be aggregated.Grid data of a spoil cone cross section are aggregated after the high resolution 2-D cross grid is established and correctly aligned and texture and bulk densities have been calculated.The degree of aggregation D A is chosen ex ante (Table 4).For example, D A = 3 aggregates values contained in a 3 × 3 grid cell cluster.Along the spoil cone "slopes", the aggregation routine checks whether the majority of cells in a cluster lie inside the spoil cone.If the condition is true, the values in the cluster are aggregated, otherwise the cluster is filled with a no-datavalue.In principle, D A can be chosen arbitrarily, but should practically be inside reasonable limits.A value of D A = 7 was found as the optimal compromise between desired high spatial resolution and available computing power (memory space) for our conditions.Aggregated data are then added to the opened ASCII save file.The save file contains for each aggregated grid point geographical coordinates (WGS 84), the height above sea level, skeleton content, coarse, medium and fine sand contents, silt content, clay content and bulk density values.

Determination of delineating surfaces and construction of a volume body
We constructed a gridded DEM of the clay base liner surface based on the available photogrammetric and borehole data (Maurer et al., 2009;Schneider et al., 2011).The subsurface catchment area of the clay base liner was determined based on a flow routing analysis using the Deterministic 8 algorithm (short: D8, O' Callaghan and Mark, 1984) in SAGA GIS (SAGA User Group Association, Goettingen, Germany).The D8 algorithm assigns runoff from one cell to one of the 8 neighbour cells, thereby allowing a clear definition of the catchment area.The subsurface catchment constitutes the lower delineating surface of the catchment's relevant sediment body (Fig. 7).
For the construction of the upper delineating surface, a DEM dated 26 November 2005 was available.This dataset is the earliest record of the initial soil surface, about 2 months after completion of the construction works on the Chicken Creek catchment.The DEM was recorded with a horizontal resolution of 1 m (i.e.1/30 arc seconds) during a routine photogrammetric survey by the VEM mine surveying department.The standard deviation in elevation values in VEM photogrammetric datasets, based on airplane altitude and camera configurations, is 0.15 m.To reduce obvious systematic errors, the DEM was improved by referencing the photogrammetric data to dGPS elevation data (Schneider et al., 2011).
From the delineating surfaces, we constructed a gridded volume body (Stratigraphic Grid or SGrid in GOCAD) in several steps as described in Schneider et al. (2011).As the surface and subsurface watershed do not exactly coincide, a 3-D block model representing the maximum extent of the water storage layer in horizontal direction was constructed based on the two delineating surfaces.A surface representing the storage layer's lateral boundary was then constructed based on the borderlines of the two surfaces.This surface was used to modify the gridded volume model in GOCAD using a procedure that "collapses" the volume of model cells that are located outside the delineating surfaces.The outline of the pond in November 2005 was digitised by visual interpretation of the model and the pond area was separated from the model (Schneider et al., 2011).

Spatial configuration of spoil ridges, volumes and masses
The spatial distribution of spoil ridges was derived by combining spectral information from aerial photographs with the spatial information contained in DEMs from the construction phase.Both aspects can be viewed simultaneously by overlying the DEMs with aerial photographs.Spoil ridges were manually digitised (Fig. 8a).Coverage of spoil ridges was incomplete due to the subsequent bulldozing of large areas.
We additionally used aerial photographs from later stages to trace spoil ridges locally, e.g. by discolouring of bare soil surfaces (probably as a consequence of different topsoil water contents) or vegetation patterns.Horizontal distances between spoil ridges, here denoted as back width R, were derived by constructing surfaces from spoil ridge curve objects and measuring vertical distances after rotation in the vertical.For the segmentation of individual ridges, horizontal distances between spoil ridges crests were calculated as vertical distances between surface objects.This way, distance values get assigned to the points/nodes of the spoil ridge curve.The property value was then written into an ASCII-file.
For the determination of the deposited and relocated sediment volumes, differential analysis of relevant DEMs was carried out.Applicable DEMs are the clay layer surface (s clay) and the initial soil surface of November 2005 (s 0511).Furthermore, the DEM recorded at 19 October 2004 (s 0410) depicts the point of time during construction when the complete mass of the sediment has been applied on both halves of the catchment, but have not yet been bulldozed in a noteworthy degree (Figs. 8 and 3).Thus, the differential volume (s 0511 -s clay) is that of the catchment's sediment body in its initial state after construction, (s 0410 -s clay) gives the spatial distribution (east and west) of the originally applied sediment masses, and (s 0410s 0511) gives the volumes abraded by subsequent bulldozing as well as the volume of the central trench that was filled from both sides by bulldozing (Fig. 8b and Table 5).The two latter difference volumes give also information about (i) the presumed average height of spoil ridges (large areas were already bulldozed to a certain degree, so true ridge heights remain somehow speculative) and (ii) to what extend these spoil ridges have been eventually "cut" by the final bulldozing process.

Spatial interpolation of structure generator data in GOCAD
After importing the ASCII data in GOCAD, the individual values were associated with the prepared 3-D volume body of the artificial catchment.This resulted in further aggregation of the data.Values were aggregated by calculating the arithmetic mean of all individual values contained in single cells.Interpolation between cell values was carried out using the GOCAD interpolation method Discrete Smooth Interpolation (DSI), which is basically a linear interpolation method that considers the geometry of spatial (geologic) boundaries (Mallet, 1992).Linear interpolation is assumed to be adequate given the high density of data and the regular distribution of values.

Calibration: spoil cone-internal structuring
Direct sampling on the catchment was not possible because invasive sampling methods would have caused a massive disturbance.We therefore collected reference data about internal bulk density distribution from spoil cones at the Wolkenberg site.The "Wolkenberg" is an artificial hill adjacent to the Chicken Creek catchment.It was dumped from November 2007 to March 2008 using the same stacker technology that was applied for the Chicken Creek catchment.The dumping of the Wolkenberg site offered the opportunity to take samples from rather freshly dumped spoil cones to gain information about internal bulk density distribution.For soil profile raster sampling, two adjacent spoil cones were chosen that had, by eye, a comparable sediment composition.Sampling was carried out by subsequently digging off three vertical cross sections through the spoil cones, each with a horizontal distance of 1 m, up to a maximal depth of 1.50 m.Each cross section profile was sampled using a specially designed cube-shaped shovel with a defined volume of 1 l.Samples were stored and transported in plastic bags.Sample points were arranged in a grid with a vertical spacing of 0.3 m and a horizontal spacing of 0.5 m (Fig. 9a).The total length of each cross section was 7 m.The coordinates (Gauss-Krueger) of each cross section's extremities were recorded with a differential GPS (dGPS, Trimble R8, Trimble Navigation Ltd., Sunnyvale).
For the determination of bulk-densities, overall 161 volume-samples were analysed.Fresh mass was determined by weighting, samples were air dried and weighed again to determine volumetric water contents.Skeleton was removed after air drying and weighed separately for all samples.Additionally, particle size distributions of selected samples were determined using wet sieving for the sand fractions and the classic pipette method for the silt and clay fractions (Gee and Bauder, 1986).The resulting point data were geo-referenced using the dGPS coordinates and visualised in 3-D using GOCAD.A central compaction zone exists within each spoil cone, albeit not clearly defined (Fig. 9b, c).Bulk densities range from 1.3 to 1.6 g cm −3 .The data were used to calibrate the parameters ρ b0 and ρ b,c .Heterogeneities in particle size distribution could not be detected within single spoil cones.This can probably be ascribed to the applied sampling method, which was too coarse to successfully capture sediment heterogeneity.

Construction of the spatial framework and uncertainties
The main requirement for the structure generator was the definition of the catchment borders.The 3-D-representation of the surface catchment has an area of 59 245 m 2 , the subsurface catchment (i.e. the approximate surface of the clay base liner) has an area of 60 898 m 2 .Slight uncertainties exist about the surface of the clay base liner due to the fact that -for the most parts -only borehole data exist.Surface and subsurface catchment area do not exactly overlap, so the exact lateral delineation of the catchment is not clearly defined.Connecting the surface and subsurface catchment borders to construct a lateral boundary surface seemed to be the most straightforward solution.The constructed 3-D-body of the entire catchment (lake area excluded) has a volume of 120 659 m 3 (cf.Fig. 7).Sediment mass balances are important for model validation, in particular when including further processes.Volume analysis (Table 5) indicates that the initially applied sediment masses (data from 19 October 2004) took up a 27 % larger volume than the sediment body after the terminal bulldozing.This finding can be ascribed to subsequent consolidation and compaction by bulldozers.If we assume initial mean bulk densities as were found at the comparable Wolkenberg site (1.45 g cm −3 ), compaction would have resulted in a calculated mean bulk density of about 1.84 g cm −3 .The result of this calculation roughly agrees with data from samples taken in May 2010 from the uppermost 80 cm of four locations in the catchment.Here, the mean bulk density value was about 1.8 g cm −3 (A.Badorreck, BTU Cottbus, personal communication, 18 July 2011).Presuming the calculated bulk density of 1.84 g cm −3 , the mass of the applied sediment was calculated to be in the magnitude of about 232 000 tons.The volume of the central trench was filled up during the bulldozing process.The volume discrepancy of about 40 000 m 3 between material removed from (s 0410) (42 000 m 3 from the west and 20 000 m 3 from the east part, totalling about 62 000 m 3 ) and the volume of the trench (about 22 000 m 3 , Table 5) must also be explained by compaction: assuming that the trench filling has the calculated final bulk density (1.84 g cm −3 ), about 51 000 m 3 of material with the original bulk density (1.45 g cm −3 ) have been inserted in the trench, leaving about 11 000 m 3 volume surplus on both sides, which most probably was also lost due to compaction processes.Information on the intensity of compaction and the associated technical processes thus needs to be incorporated in the future.
Analysis of dumped ridge heights showed that maximal heights of dumped sediment account for almost 12 m above the clay base liner.One tier of spoil ridges was applied in the eastern part and two tiers of spoil ridges were applied in the western part.However, the analysis also suggests that the uppermost tier in the western part, i.e. more than 6 m of sediment, was completely displaced by bulldozing.The average height difference for (s 0410 -s clay) was calculated to be 2.56 m; however, this value does not necessarily reflect mean spoil ridge heights, as most areas have already been bulldozed.Results of (s 0410 -s 0511) suggest that an average of about 0.54 m was cut off from the initial ridges.Also, maximum values of (s 0511 -s clay) indicate heights of 3.5 m (east) and 3.7 m (west).These maximum values -on the other hand, are possibly also a result of the displacement of material from peripheral spoil ridges towards the centre.It must also be presumed that stacker dumping is configured to produce constant ridge heights.We therefore assumed a constant spoil ridge height of 3 m.

Generated 3-D structures
The model scenario discussed in the following is just a single realization that was used to demonstrate the model.Results of the process-based structure generator approach are shown exemplarily for one realisation of structural heterogeneity of the eastern part of the Chicken Creek catchment (Fig. 10), denoted as volume of interest in the following.The modelled volume of the sediment body has a surface extension of 18 568 m 2 , almost 1.9 ha, which constitutes about a third of the area of the entire artificial catchment (see Sect. 2.1).Based on aerial image analysis, it was assumed that the whole volume was regularly filled with spoil ridges (Fig. 10a).The manually specified spatial location of the structural elements (spoil ridges) may fluctuate within narrow margins (1-2 m) because of subjective errors during digitisation and the incomplete spatial information available.We expect though that these uncertainties have only a small impact on the overall results, as relative configurations, distances and -roughly -the strike azimuths of spoil ridges can still be relatively well derived and extrapolated from the visible parts of spoil ridges on the aerial images.Data were imported as aggregated values.For the present realisation, a predefined (dimensionless) aggregation factor D A = 7 was used.That means that 7 × 7 data points were combined into one by calculating the arithmetic mean (see Sect. 2.2.6).One value thus covers an area of 49 cm 2 in a modelled spoil cone cross-section.The defined distance between cross-sections along the spoil ridges was 3 m, and Averaging caused slight distortions of the value distribution and statistical parameters of the original dataset.As can be seen for the example of the fine sand content, upscaling and interpolation result in slight increases of mean and median values, whereas standard deviation and variance are reduced (Table 6).Likewise, the 25th and 75th percentiles are slightly moving towards the arithmetic mean.Value distribution shifts from rudimental bimodal towards unimodal normal distribution after interpolation (Fig. 12).

Patterns of particle size distributions and bulk density
Spatial patterns of generated structural heterogeneity are similar to those observed in reality (Fig. 11): on both the aerial image and the spatial model depicted in Fig. 11, distinct surface pattern elements can be identified.The ellipse markings on both images can be described by three parameters (azimuth, length, width).A visual inspection suggests that the modelled pattern elements can be described by a similar range of parameters than the observed pattern elements.These similarities are a result of the sediment mixing and pattern imitating approach used for modelling: (i) the feature azimuth is governed by the course of stacker trajectories, (ii) their width is governed by the determined back width R, which was derived from stacker trajectory geometry, and (iii) their length is governed by the limitations of the random number N of consecutive spoil cone cross sections with the same base sediment composition.The trajectories are crucial for establishment of the basic sediment structure of the catchment.Here this information has already been available from aerial photographs and at the moment was not considered as a random process.In future model development steps, the dumping trajectories will be included in order to increase flexibility for the generation of realizations of catchment structures.The average particle size distribution (here for fine sand fraction) in the catchment model does not exactly resemble the distribution that was measured on samples from the 20 m × 20 m grid (here for fine sand, Tables 2 and 6, Fig. 12  panel A1).A direct comparison of particle size distributions is not indicated at the moment, as in the present version of the structure generator, the parent material composition is generated as a (stochastic) mixture of all eligible geologic units at the excavation site (Tables 3 and 6).The generated mixtures thus do not necessarily represent the actual mixtures that were dumped for catchment construction (cf.Table 2).This effect is intentional, as our main goal was to reproduce heterogeneity patterns in a first step.The comparison between observed average bulk densities in noncompacted spoil cones and average bulk densities in the volume of interest (Table 7) indicated an overestimation of bulk density in the realisation of the spatial model examined in this study.Here, mean values of 1.56 g cm −3 were calculated (respectively 1.53 g cm −3 in the dataset before interpolation), whereas the mean value from spoil cone sampling was 1.46 g cm −3 .Although we defined lower basic bulk density values for the impact zone (1.4 g cm −3 ) as proposed by Matschak (1969), calculated values seem to be still too high.Overestimation of bulk densities can be attenuated in future model realisations by further adjusting the relevant input parameters, e.g.stacker drop height (here H D = 5 m).One option to resolve the problem of poor density estimation is to condition the modelling to measurement data.An applicable method for conditioning spatial models with local hard data was proposed by Michael et al. (2010): they combined process-based geologic modelling with multiple-point geostatistical simulations using training images from objectbased modelling.The results were geologically realistic spatial models that were fully conditioned to measurements.Similar combinations of a process-based structure generator and geostatistical data conditioning approaches have been presented elsewhere (e.g.Teles et al., 2004;Reza et al., 2006).Such comparisons with a well-characterised study site are imperative for the evaluation of hydrological models (de Marsily et al., 2005).

Conclusions
We present a structure generator capable of reproducing characteristic heterogeneity patterns of dumped sediments on an artificially constructed hydrologic catchment.The approach is based on information about technical processes at two spatial scales including dumping and particle segregation within individual spoil cones and the distribution of spoil cones along dumping trajectories.The structure generator creates distributed texture (skeleton; fine, medium, coarse sand; silt; clay) and bulk density data in relatively high resolution for the 2-D-vertical cross-sections and spatially aggregated for the 3-D catchment scale.The "scalability" of the model is useful when analysing effects of sediment structures on flow processes with hydrological models.The 3-D structures and patterns can be used to investigate the effects of heterogeneity on different scale levels.
The results suggest that spoil cones with a compacted zone are forming characteristic spatial patterns typical for the applied technology.Particle segregation additionally modifies sediment properties.The effect of the patterns with respect to the distribution of soil hydraulic properties is clearly important for hydrological analyses.The stochastic mixing of parent materials and the consideration of stacker trajectories produces structural heterogeneities that are similar to observed spatial patterns.Moreover, modelled bulk density variations are satisfactorily reflecting the spatial dumping structure; still the somewhat overestimated absolute density values need additional calibration.We expect that these features play an important role in the catchment's hydraulic behaviour, which can be verified by future hydrological modelling based on such spatial models.The results indicate that a physically-based modelling of sediment structures is possible.Uncertainty in catchment modelling can be easily considered by introducing stochastic components at each step of the structure generation process.
We used the artificially created catchment "Chicken Creek" with relatively well-defined conditions here simply for developing the structure model.This allowed for a better testing of novel ideas of structure model development than alternative natural catchments.Despite the artificial catchment being a unique experimental site, our modeling approach is more general and in principle not restricted to this catchment: (i) the technological processes described in the structure generator model are ubiquitous: they are occurring everywhere where unconsolidated sediments are moved and dumped.This applies not only for most of the world's open cast mining areas, but also for landfalls and construction sites, among other larger scale earth moving operations.The presented structure modeling approach can be adapted to any of the above situations.(ii) Furthermore, the presented study can lead the way for similar studies that are trying to link a structure generator approach with the functionality of advanced geological modelling software.This is, of course, not limited to technogene structures, but can be adopted as a simplifying modeling of sediment structures that result from geomorphological processes or for other geological settings using results of more detailed process models.
T. Maurer et al.: Initial sediment distribution of an artificial hydrologic catchment opposite sides of rows 1 to n, flipping columns 1 to m according to the following matrix (also see Figure A1):

44Figure 3 .
Figure 3. Aerial photograph (October 2004) showing the construction layout of the Chicken Creek catchment.Spoil ridges were digitised based on this image.

Fig. 3 .
Fig. 3. Aerial photograph (October 2004) showing the construction layout of the Chicken Creek catchment.Spoil ridges were digitised based on this image.

TFig. 4 .
Fig. 4. Conceptual model of the catchment construction steps: (A) modelling of the geological conditions at the excavation site and material transport (not yet included).(B) Relevant parameters and processes considered for spoil cone construction.(C) 3-D-configuration of spoil cone cross sections.(D) Incorporation of the structure generator data in the gridded catchment model and inclusion of bulldozing and cone dcraping processes (not yet included).

Fig. 5 .
Fig. 5. Simplified flow chart diagram of the structure generator programme.Sketches on the right hand side illustrate the cone construction steps.

Fig. 6 .
Fig. 6.The spatial alignment of 2-D cross sections (black) was derived from the spatial orientation of stacker trajectories (blue).

Fig. 7 .
Fig. 7. From surface to volume body: the DEM of the clay base liner, which was constructed from photogrammetric and borehole data; was used together with the DEM of the initial ground surface to construct a volume body (Stratigraphic Grid, SGrid) in GOCAD.

Fig. 8 .
Fig. 8. (A) Location of spoil ridges (blue line) obtained from manual digitisation of aerial images and combined with the surface model of the Chicken Creek catchment (shape indicated by the red line) from the construction phase in October 2004, and (B) vertical distances between the surface model of the construction phase (October 2004) and the surface model from November 2005, when construction was finished.The volume difference obtained from the comparison of surface models allowed to calculate volumes of the dumped sediment masses and to derive mean spoil cone heights.

Fig. 9 .
Fig. 9. 3-D soil bulk density distribution obtained from a freshlydumped spoil cone location: (A) spoil cones during sampling with sampling scheme.(B) 3-D visualization of bulk density after linear interpolation (DSI) in GOCAD.(C) Spoil cone zones with highest bulk densities (compaction).

Fig. 10 .
Fig. 10.Structure generator results shown as colour-coded 3-D point data in GOCAD.(A) Cross-section of a single spoil cone, featuring a compacted central zone; (B) side view of several spoil ridges, showing the lateral alignment of single cones; (C) spoil ridges on the complete eastern half of "Chicken Creek".

Fig. 11 .
Fig. 11.Comparison of generated sediment heterogeneity patterns (left) and patterns on the eastern part of the artificial catchment associated with sediment heterogeneity (right).The ellipses mark zones with relatively low sand contents -and respectively high contents of silt and clay -on the left side, and zones with relatively dense vegetation cover and/or crusting on the right side.

Fig. 12 .
Fig. 12. Frequency distributions of cell values for the fine sand content, given in % on the x-axis for the simulated eastern part of the Chicken Creek catchment.(A) Original frequency distributions from the structure generator are bimodal.The white bars in A1 show the histogram for measurement data from raster grid sampling.(B) During aggregation of the original values in the raster grid cells of the volume body in GOCAD, frequency distributions are modified.(C) After linear interpolation (DSI) between grid cells, the bimodal distribution has disappeared.

Table 1 .
Data used for deriving spatial boundaries and determining sediment properties in the structure generator model.

Table 2 .
Texture data (fine earth fraction particle distributions) from raster soil sampling carried out in 2006, showing the differing values for the eastern and western part of the Chicken Creek catchment.The "eastern" and "western" part comprise only the areas where material was dumped with stacker technology.Courtesy of subproject Z1, SFB/TRR 38.The simulated textures (model) are for the current mixing approach after spatial interpolation.aFor the catchment north of the clay dam.

Table 3 .
Properties of the principal types of parent material from the excavation site.Values based on data courtesy of VEM.

Table 4 .
Input parameters used in the present model scenario and qualitative assessment of parameter uncertainty regarding modeling results.Csclmax, the value of ds m can have a higher or lower impact on internal cone structuring; parameter value is based on literature and theoretical considerations.
α35 ˚slope angle controls dipping of cone internal layers and spoil ridge -and thus of potential hydraulic interfaces; value is from literature and observations.β15 ˚obliqueness angle affecting the distribution of bulk densities; is expected to be more noticeable when considering smaller scales; value is derived from conceptual considerations.ρb0 1.2 g cm −3 a general shift of average bulk densities; value is indirectly based on experimental analysis of two spoil cones.

Table 5 .
Deposited and translocated sediment volumes during Chicken Creek construction as calculated in GOCAD.Positive values indicate that the surface during construction (when material was piled up) lies above the reference surface (i.e. the final bulldozed surface).Negative values indicate that the construction surface is below the reference surface (e.g.due to the filling of the central trench).

Table 6 .
Statistical values for the property "fine sand content" (in %) after aggregation and interpolation in GOCAD.Also see Fig.12a-c.

Table 7 .
Statistical values for particle size distributions (in %) and bulk density after linear interpolation (Discrete Smooth Interpolation, DSI).Size intervals are given in µm.The volume of interest (eastern part of the catchment) contained 230 587 cells.