3D GIS Modeling of Soft Geo-Objects: Taking Rainfall, Overland Flow, and Soil Erosion as an Example

Abstract

In physics, objects can be divided into rigid and soft objects according to the object deformation capacity. Similarly, geo-object can also be classified into rigid geo-objects (e.g., building, urban) and soft geo-objects (e.g., mudflow, water, soil erosion). There are three types of approaches for 3D GIS modeling, i.e., surface-based, volume-based, and hybrids in terms of geometry. These approaches are suitable for representing rigid geo-objects, but they are not suitable to simulate the intrinsic properties of the soft geo-object, i.e., dynamics and deformation. And so far there are few GIS modeling methods for simulation of soft geo-objects. GIS flow elements (FEs) and GIS soft voxels (SVs) were proposed for 3D modeling of soft geo-objects. GIS flow elements can realistically represent the dynamics and stochastics of soft geo-objects, while GIS soft voxels simulate deformation of soft geo-objects. The authors discuss the implementation and computer programming of GIS flow elements and GIS soft voxels in this study. GIS FE and SV have been successfully applied in a case study toward the simulation of the process of rainfall, overland flow, and soil erosion. A software system has been designed and developed, which has the functions of data management, model computation, and 3D simulation.

1. Introduction

Nowadays, severe shortage of water resources, ecological destruction and environmental pollution, global changes, natural disasters, etc., are the key issues of geosciences. Main research objects of these key issues such as water, polluted air, and mudflow are soft geo-objects. Modeling and 3D visualization of soft geo-objects is emerging research area. In computer graphics, several approaches have been applied for simulation of soft objects, e.g., the particle system approach and the metaball approach [1]. However, the particle system approach and the metaball approach are driven by physical force objects, which are not suitable to simulate geographic process that is driven by more complex geomodeling [2]. Mitasova et al. used densities of particles to sample rainfall excess and sediment transportation of sand and clay. Compared to traditional sampling method, their method showed several advantages [3]. For example, it can be easily extended into arbitrary dimensions, and is fairly straightforward to be implemented in a multiscale framework with data adaptive capabilities. But the method has obvious limitations. Firstly, it cannot accurately represent the dynamic change of sample points’ velocity and direction over space and time because all sample points have the same size. Furthermore, the geographic phenomena of soil separation and fusion during the process of soil erosion have not been represented. Soil erosion is a naturally occurring process in land which refers to wearing away of a field’s top soil by natural forces of water and wind. How the sediment transports, how the soil separates, and how much soil losses, all these are driven by geomodeling.

2. Methodologies

In GIS existing research studies focus more on rigid objects such as mountains, roads, and buildings, and few methods have been proposed for the modeling and 3D visualization of soft geo-objects. GIS flow elements (FEs) and GIS soft voxels (SVs) were proposed and developed for 3D modeling of soft geo-objects. GIS FEs can realistically represent the dynamics and stochastics of soft geo-objects while GIS SVs can simulate deformation of soft geo-objects.

2.1. GIS FE

A GIS FE is a basic simulation unit and spatially corresponds to a pixel in remotely sensed imagery. It is characterized by the position (i.e., x, y, z coordinates), velocity, and direction of a soft object, but volume is neglected. GIS FE is driven by geomodeling with the objective to simulate the dynamic feature in soft geo-objects. Although a GIS FE is within a pixel, it can take many appearances such as a point, a line segment, or a surface depending on the natural appearance of the soft geo-object. For example, the FE is with line segment shape for rainfall simulation. The velocity and direction of a FE can well reflect the real dynamic change. Besides the above-mentioned fundamental attributes, a FE is able to carry more properties such as color and texture that help distinguish and characterize an object. This provides flexibility for extension.

2.2. GIS SV

A GIS SV is a basic unit for deformation simulation. Similar to FE, a GIS SV is also based on a pixel from remotely sensed imagery and it also features with position, velocity, and direction, which are controlled by geomodeling. But a GIS SV carries volume information. A SV is covered by an isosurface that can well represent volume shape and surface deformation. A GIS SV has potential to carry more information by adding colors and textures or by modeling its internal structure.

2.3. Calculation of basic parameters

Here the basic parameters include direction, velocity, shape, and volume.

8-neighborhood tracing algorithm is used to calculate the direction of a FE. Eq. (1) is used to compute the velocity of the FE (V) where M₁ is geoscientific model and p₁, p₂, …, p_j are parameters affecting V. Suppose that V is directly proportional to the length of a GIS FE (L), we get (Eq. 2) where L is a proportional coefficient with a value greater than 0. (Eq. 3) is used to control the shape of the GIS SV where g represents isosurface, d is the length of a pixel, h is the average thickness of the soft geo-object, and r² = x² + y² + z² in which (x, y, z) are 3D coordinates of the critical point of the GIS SV. Eq. (4) calculates the volume Vol of the GIS SV.

$V=M_1\left(p_1,p_2,\dots \dots ,p_j\right)$E1

$V=\lambda L$E2

$\mathsf{g}\left(\mathsf{r}\right)\mathit{=}{\left(\mathsf{r}\mathit{‐}\frac{\sqrt{\mathsf{2}{\mathsf{d}}^{\mathsf{2}}\mathit{+}{\mathsf{h}}^{\mathsf{2}}}}{\mathsf{2}}\right)}^{\mathsf{2}}$E3

$Vol=d^{2}h$E4

The following section will introduce the application of using GIS FE and GIS SV theories in modeling and 3D simulation of rainfall, overland flow, and soil erosion.

3. GIS FE-based simulation of rainfall

The objective is to simulate raindrops falling from the sky to the ground surface.

3.1. Raindrop dynamics

Force objects on a raindrop include gravity, air buoyancy, air resistance, wind force, and the kinematic equation for a raindrop in vertical direction is [4]:

$m\frac{dv}{dt}=\mathit{mg}-F_R-F_B$E5

where m represents raindrop mass, v represents raindrop falling velocity, t represents time, g means gravitational acceleration, F_R means air resistance, and F_B means air buoyancy. Figure 1 shows the force objects of a raindrop.

Figure 1.

A schematic view of raindrop dynamic analysis.

3.2. Criteria for raindrop GIS FE representation

In this study raindrop GIS FE representation meets artificial rain experiments criteria [5]:

• Raindrop particle size distribution is close to natural rainfall. Natural rainfall raindrop sizes range from near zero to about 7 mm. The median particle size of an erosive rain storm is between 1 and 3 mm. Raindrop diameters normally increase with the increase in rainfall intensity.

• Raindrop impact velocity is close to the natural raindrops. Raindrop impact velocity, from droplet velocity near zero to the maximum raindrop velocity of more than 9 m/s. The landing speed of an ordinary raindrop with a diameter of 2 mm is 6–7 m/s.

• Rainfall intensity is close to the natural rainfall. Natural rainfall intensity from near zero to a few millimeters per minute. In general, low rainfall intensity is not important to soil erosion, and the frequency of high rainfall intensity is very low, so that the importance is limited. Common rainfall intensity of 0.2–2 mm/min is usually the most important rainfall intensity.

• Throughout the study area, raindrop characteristics and strength is fairly uniform.

• The rainfall is continuously simulated in the study area.

• The impact angle of most raindrops is not too much deviated from the vertical line.

3.3. Raindrop GIS FE representation

The representation includes geometry representation and dynamic representation.

Raindrop is considered to have a shape of a combination of a taper and a semisphere with white color and it is transparent. The initial position, x and z are randomly created and y position is on the top of the viewpoint. The velocity and direction are determined based on raindrop dynamics. The raindrop object has its lifespan, which ends when raindrop collides with DEM.

3.4. Programming

The development platforms are Visual Studio C++ and OpenGL.

3.4.1. Define raindrop array

Create raindrops by meeting artificial rain experiment criteria. The raindrop array stores raindrop total number, color, transparency, and coordinates information. 3D coordinates of the start point of rain line are generated with a random function. Coordinates of the end point of the rain line are calculated based on dynamic analysis.

3.4.2. Create animation

Create raindrops animation by meeting artificial rain experiment criteria as well and the pseudo code is as follows:

For i from 1 the total number of raindrops

{

Translate the start and end points of the rain line to a new position along z axis based on the initial velocity setting

1. Add an acceleration increment to the initial velocity of the raindrop

2. If a rain line touches or penetrates DEM, then

1. {

2. The life of the rain line is ended

3. Reinitialize a new rain line and set the initial velocity

4. }

}

4. GIS FE-based simulation of overland flow

The objective is to simulate the velocity and direction of overland flow.

4.1. Flow streamline dynamics

Saint Venant kinematic equation of unsteady flow of water is used as the governing equation [6]:

$\Delta M_1+\Delta M_2=P_U-P_L+W_x-T$E6

or it is written as:

$V\frac{\partial V}{\partial x}+\frac{\partial V}{\partial t}+g\frac{\partial y}{\partial x}=g\left(i_0-i_f\right)$E7

where P_U and P_L represent pressure on the upstream face and downstream face, respectively; W_x represents the gravity component in water flow direction; T is friction resistance, ΔM₁ and ΔM₂ are local momentum change and transport momentum change, respectively, x is distance in water flow direction, t is time, y is the depth of water, i₀ and i_f are bottom slope and friction slope, respectively. Flow streamline dynamic analysis is shown in Figure 2.

Figure 2.

A schematic view of flow streamline dynamic analysis.

4.2. Compute the velocity of overland flow

Overland flow velocity is an equation of discharge per unit width and slope angle [7]:

$V=K{q}^n{S}^m$E8

where K, m, n are parameters. Compared to laminar flow and turbulence flow equations, normally the value of m/n is between that of laminar and turbulent flow. By nonlinear regression analysis, we deduced the equation

$V=21.8811q^{0.4667}{S}^{0.2510}$E9

where S and q represent the slope angle and the discharge per unit width, respectively. The water discharge was computed using the equation [8]:

$q=x\left(I-f\right)cos\beta$E10

where x represents the average slope length from the slope top, β represents the average slope gradient, I represents the rainfall intensity, and f represents the infiltration rate.

4.3. Compute overland flow direction

Water and sediment discharge computation based on grid DEM is usually determined using single flow path algorithm, i.e., the method of determining the maximum gradient. For a 3 × 3 window, the center cell has eight neighbors. The water and sediment flow direction coding of each cell is based on the digital coding method in Refs. [8, 9]. For example, if the flow direction of water and sediment of a grid unit as the center of the window is due west, i.e., water and sediment in the center of the window flow into the adjacent cell 4, then the flow direction value of the center cell is 4.

In the algorithm, following values are combined for identification of landform structure.

• DEM elevation values, stored in the array ALT[i, j].

• Flow direction values, stored in the array PTR[i, j].

• Flow streamline coordinates, stored in the array FLOWLINE[i, j].

The above arrays have the same size. Each grid unit has a value to identify one of its attributes. In addition, the algorithm uses the following terms:

• Outflow point: if water and sediment flows out of a grid unit, then the unit is called an outflow point.

• Inflow point: if water and sediment flows into a grid unit, then the unit is called an inflow point.

• Flow point: if water and sediment flows from grid cell (i, j) to grid cell (x, y), then (x, y) is a flow point of (i, j).

In fact, in addition to certain points along watershed boundary that only have outflow, other points have both inflow and outflow. Algorithm:

In a 3 × 3 window, gradients from the center cell to its 8-neighborhoods are used to determinate flow directions.

• For cells at the edge of DEM or boundary of the study area, the flow direction of each cell is defined as the direction toward the boundary.

• For any other grid cell, calculate the cell’s elevation gradients to its eight neighbors. EG_{0, i} = Z₀ − Z_i (i = 1, 3, 5, or 7) represents the elevation gradient in horizontal or vertical directions, and EG_{0, i} = Z₀ − Z_i/sqrt(2) (i = 2, 4, 6, or 8) represents the elevation gradient in diagonal directions.

• Determinate the neighboring cell which has the maximum elevation gradient.

• Identification of isolated depression flow direction. Scan the study area using a 3x3 window, and (i, j) is the center cell of the window. Calculate the eight elevation gradients of (i, j). If the maximum gradient value is less than 0, then identify (i, j) as an isolated depression.

• Identification of outflow point flow direction. Scan the study area using a 3×3 window, and (i, j) is the center cell of the window. Calculate the eight elevation gradients of (i, j). If the maximum gradient value is greater than 0, then identify (i, j) as an outflow point.

4.4. Flow streamline GIS FE representation

The representation includes both geometry representation and dynamic representation.

Flow is considered to have a fine cylinder shape with Cambridge blue color and it is transparent. The initial position starts from the intersection point between raindrop and collision plane on DEM. Its velocity and direction is determined based on the analysis of dynamics. The lifetime of the flow streamline ends when it runs into a channel.

5. GIS FE- and SV-based simulation of sediment transport and soil erosion

The objective is to simulate sediment transport and the process of soil erosion by water.

5.1. Sediment particle dynamics

Dynamic analysis of a sediment particle is shown in Figure 3 [9]. The parameters are W as gravity, P_y as uplift, P_x as traction force, and T as upper-surface friction. In this simulation, only suspended load is considered. The velocity of suspended load in water flow direction is mostly equal to that of flow streamline [10]. Bed load and saltation load will be considered in future work.

Because spatially varied forces acting on flow streamline and sediment particles are too complicated to be precisely represented, we firstly use an 8-neighborhood tracing algorithm to compute the flow direction on hillslopes, then apply the modeling theory of remote sensing information model combined with experimental results to obtain an equation for flow velocity computation on hillslopes to simulate dynamic sediment laden flow in 3D space.

Figure 3.

A schematic view of sediment particle dynamic analysis.

5.2. Erosion voxels

We form GIS SVs on a pixel basis and name it erosion voxel to simulate the separation and fusion of soil mass. The erosion voxel combines the attributes not only from geographic pixels, but also from particles and metaballs in computer graphics.

(1) The main characteristics of a pixel-based erosion voxel are described as follows:

• The appearance of an erosion voxel is designed as a circumscribed ellipsoid of a cuboid. The width and length of the cuboid are equal to cell length, and the height of the cuboid represents the averaged thickness of soil loss per unit duration (which depends on specific time, location, and amount of rainfall);

• Represents averaged volume of soil loss per unit duration; therefore, it has statistical meaning;

• Represents the highest precision of geographic image data, and includes soil properties, vegetation coverage, slope, elevation (represented by the elevation of center point of top face of its inscribed cuboid) and other information;

• Once time reaches the specific step of duration, we separate a specific erosion voxel from a specific pixel, and it will run off with water (i.e., erosion voxel separation);

• When the velocity falls to zero, the erosion voxel will subside, and the fusion of erosion voxels will occur;

• When the erosion voxel runs into a channel, its lifetime will be ended.

(2) Quantitative expression of an erosion voxel

An erosion voxel is represented by a circumscribed ellipsoid:

$\frac{{\left(x-x_0\right)}^2}{{\left(\frac{\sqrt{2}}{2}d\right)}^2}+\frac{{\left(y-y_0\right)}^2}{{\left(\frac{\sqrt{2}}{2}d\right)}^2}+\frac{{\left(z-z_0+\frac{h}{2}\right)}^2}{{\left(\frac{\sqrt{2}}{2}h\right)}^2}=1$E11

where d represents cell length, h represents averaged thickness of soil loss per unit duration, and (x₀,y₀,z₀-h/2) are the coordinates of the center point of an erosion voxel. The volume of a GIS SV (Vol) is computed based on its inscribed cuboid (Figure 4).

$V=d^{2}h$E12

where d is the cell length and h is averaged thickness of soil loss per unit duration.

(3) Elevation

Real surface elevation is represented by the elevation of the center point of the upper surface of the inscribed cuboid of the erosion voxel.

(4) Transport routes

Considering that an erosion voxel only has statistical meaning, we treat its transport routes the same as sediment laden flow [2].

(5) Structure parameters of an erosion voxel

Include center-point coordinates, strength, color, transparency, timer, etc. Other attribute parameters such as soil type, vegetation coverage, and slope angle can be added as needed.

(6) Separation and fusion

Use the basic methods for rendering GIS SVs [11]. Considering that an erosion voxel only has statistical meaning, force analysis to contact surface and volume control are not performed.

Figure 4.

The profile of a GIS SV.

6. Results

Based on the aforementioned technology, including software module design, algorithm design, and pseudo code description, for the implementation of the new methodology, we develop a Modeling and 3D Simulation System for Water Erosion on Hillslopes (M3DSSWEH). The system includes three modules: module of DB management, module of model computation and verification, and module of 3D simulation. The system has been designed and developed based on the platform of Oracle, Visual C++, and OpenGL. It has a user-friendly interface based on human-machine interactive techniques and the advanced module design makes it flexible and easy for function extension.

Figure 5.

User interface of M3DSSWEH.

• User interface (UI)

The system UI consists of menus, toolbars, floating panels, and viewports (Figure 5).

• GIS simulation module

The module simulates rainfall, overland flow, and soil erosion process:

Rainfall simulation in this study references standard artificial rain experiments. Set raindrop intensity and raindrop diameter based on rainfall intensity and pixel size to enhance scientific simulation and further facilitate the simulation of raindrops splash effect.

Overland flow simulation is based on GIS FE, which is represented by a fine cylinder. The height of the cylinder is in direct proportion to the velocity of overland flow; the inclination of the cylinder represents the direction of the movement; the diameter of the cylinder is in direct proportion to the depth of overland flow; and the depth of cylinder color is related to sediment concentration. The deeper the cylinder color is, the greater sediment concentration it represents.

Sediment simulation is also based on GIS FE, which is represented by a sphere. The radius of the sphere is in direct proportion to sediment size, and the color of the sphere is related to soil type.

Figure 6.

The GIS simulation module of M3DSSWEH.

Soil erosion simulation is based on GIS SV and the rendering of the SVs are using the aforementioned template rendering algorithm (Figure 6).

• Data source

As an implementation, the data source is from a research area located between 39°43’37”–39°46’28”N, and 111°7’7”–111°9’14”E, with an area of 3.85 km². It is Wufendigou watershed in Inner Mongolia, belonging to China Loess Plateau region. The watershed is characterized by severe water and soil loss and this research will be applicable to improving the understanding of soil erosion in the China Loess Plateau region.

Theory and methodology of GIS FE and GIS SV will be further extended to better serve geoscience in the future.

Figure 7.

GIS simulation of terrain, light, shadow, rainfall, and overland flow.

Figure 8.

GIS simulation of soil erosion.

• System output (Figures 7 and 8)

• System advantages

The system shows rainfall, overland flow, and soil erosion simulation using GIS FEs and GIS SVs, and the results are satisfactory. Compared to traditional GIS models such as TIN, grid, tetrahedral, and octree, it is more convenient, vivid, and efficient to use GIS FE and GIS SV to simulate dynamic soft geo-objects. Combining GIS FE and GIS SV with the traditional GIS models, any geo-object in solid, liquid, or gas phase can be well represented.

• System prospects

7. Conclusions

This paper discusses the implementation and computer programming of GIS FE and GIS SV. Based on a case study in the China Loess Plateau region, we used GIS FE and GIS SV to simulate rainfall, overland flow, and soil erosion, respectively. Implementations show that the spatiotemporal changes of sediment-laden flow can be more intuitively and realistically simulated after GIS flow elements and GIS soft voxels technology are integrated into the system.

Both GIS FE and GIS SV are proposed based on a pixel from remotely sensed imagery that facilitates the data acquisition particularly with the significant use of remote sensing technique in geosciences and its integration with GIS. On the other hand, the pixel level data can be directly used in the calculation without further spatial operation so as to ensure the data accuracy. Moreover, both the GIS FE and GIS SV are driven by reliable geoscientific models. Therefore, GIS FE and GIS SV together can be used as basic units for simulating soft geo-objects and have the potential for practical use in other research areas such as flood modeling and simulation.

Acknowledgments

Special thanks to Professor Ainai Ma and Professor Shanjun Mao from Peking University, and Professor Hui Lin from the Chinese University of Hong Kong for their help and advice.