Mathematical vector framework for gravity-specific land surface curvatures calculation from triangulated irregular networks

ABSTRACT Land surface curvature (LSC) is a basic attribute of topography and influences local effects of gravitational energy and surface material transport. However, the calculation of LSCs based on triangulated irregular networks (TINs) has not been fully studied, which restricts further geoscience studies based on TIN digital elevation models (DEMs). The triangular facets and vertices of a TIN are both expressions of the land surface; therefore, based on their adjacency relationship, the LSCs can be calculated. In this study, we propose a mathematical vector framework for LSC calculation based on TINs. We define LSCs from the perspectives of the curvature tensor, slope and normal contour direction vectors, and then provide the calculation operators for LSCs based on both TIN triangular facets and vertices. Next, based on a mathematically simulated surface, we find that the TIN-based method exhibits similar effects on the scale as the grid-based methods and very low error sensitivity. In addition, based on different real landform cases with various data sources, we perform experiments involving land surface concavity–convexity and hillslope unit classification by using the TIN-based method. The results show that the TIN-based method can enhance the performance of TINs in landform classification over grid-based DEM methods. The proposed mathematical vector framework for LSC calculation can improve other geoscience studies based on TINs.


Introduction
The concept of curvature, originally derived from mathematics, is typically used to quantify the degree of bending of a surface or a line (Wilson 2018). In geoscience studies, curvatures specifically refer to land surface curvatures (LSCs), which reflect characteristics of land surface structures and morphologies in different directions (Minár, Evans, and Jenčo 2020). The system of LSCs has been developed and successfully applied in geomorphology, physical geography and geology, environmental sciences, hydrologic processes, remote sensing, ecology, soil science, etc (Lv et al. 2017;Xiong et al. 2021). Krcho (1973) proposed the general theory of LSC from a gravitational perspective, and then Shary (1995) developed this theory including gravity-invariant curvatures and defined the various relationships between them. Gravity-invariant curvatures, such as maximum, minimum, Gaussian and mean curvatures, are considered in classical differential geometry; they describe the local shape of a 2-dimensional manifold in 3-dimensional space and are also called geometric curvatures. However, the land surface is in the gravitational field, and surface material movement is simultaneously influenced by both the terrain shape and gravity (Mitášová and Hofierka 1993;Evans 2012). Therefore, the gravitational processes at the Earth's surface are difficult to assess by using gravityinvariant curvatures. Therefore, gravity-specific curvatures, such as profile, tangential and plan curvatures, have been defined and are commonly used in geoscience. Gravity-specific curvatures dominate in geoscience compared to gravity-invariant curvatures (Minár, Evans, and Jenčo 2020). In recent years, a mix of gravity-specific and gravity-invariant curvatures has been used in land surface classification (Dekavalla and Argialas 2017) and the geomorphological mapping of submarine areas (Moskalik et al. 2018). The definitions of the curvatures mentioned above are further discussed in the methodology section.
Due to the continuity of terrain representations and aiming at a simple data structure, regularly sampled grid-based digital elevation models (Grid-DEMs) are generally the basis for LSCs calculations. The calculation of LSCs is related to the second-order derivatives of the land surface. Under the assumption that the surface is at least C 2 smooth, the local surface fitting method in a 3 × 3 window is generally used to calculate the first-and second-order partial derivatives to obtain the LSCs (Hengl and Evans 2009). The full quadratic, constrained quadratic, and incomplete quartic functions are the most commonly used local surface fitting methods (Evans 1980;Zevenbergen and Thorne 1987;Shary 1995). In addition, benefitting from the equidistant sampling of Grid-DEMs, numerical differentiation methods are also popular for calculating LSCs (Moore et al. 1993). According to previous studies, these LSC calculation methods based on Grid-DEMs can provide sufficiently accurate results, especially for mathematically simulated surfaces (Schmidt, Evans, and Brinkmann 2003). However, real landforms are complex and not as smooth as mathematically simulated surfaces. Hu et al. (2021a) reported that the LSCs results are often suboptimal in complex terrain feature areas due to various errors and uncertainties in Grid-DEM generation. Numerous studies have shown that these LSCs calculation methods are very sensitive to DEM error (Florinsky 1998;Jochen, Evans, and Brinkmann 2003). In addition, geoscience research involves many multiscale characteristics; thus, algorithms should be less affected by the landform scale.
The triangulated irregular network (TIN) is introduced as another DEM structure, mainly for the following two considerations (Wolf 2004): first, a Grid-DEM cannot provide land surface details and limit data redundancy simultaneously; second, the level of detail of land surface expression should match the scale of geographic objects and processes rather than the sampling rules. In contrast to Grid-DEMs, TINs can provide sufficient detail in complex terrain and simplify representations in gradual change areas across the land surface by using various sampling densities and strategies. As the measurement techniques are developed, the TINs in previous studies based on sparsely sampled surfaces can be derived from dense point clouds from LIDAR or photogrammetric techniques (Wilson 2018). Some scholars have suggested that TINs be used in place of Grid-DEMs in land surface analysis and GIScience, and related studies have rapidly expanded in recent years. Currently, TINs play important roles in land surface visualization, terrain feature extraction, visibility analysis, geomorphic evolution modeling, spatially distributed hydrological models, etc (Tucker et al. 2001;Yang, Shi, and Li 2005;Li, Wang, and Yang 2008;Refice, Giachetta, and Capolongo 2012;Wu et al. 2014;Noh and Howat 2015;Barnhart et al. 2019). Accordingly, analysis methods based on TINs should be updated and developed. This shift is more than a data structure change, and the most important difference is reflected in the progress of analytical methods and ideas. The problems of scale have always been an important issue in geoscience research, and TINs and the corresponding analysis methods may provide a feasible solution. Terrain parameter calculations can be performed based on TINs, but there are still many deficiencies in existing research. Different from Grid-DEMs, a TIN is constructed based on the topological relations among vertices, edges and triangular facets with high complexity. Thus, the local surface fitting and finite difference methods developed for Grid-DEMs are difficult to directly apply to TINs. For example, slope and aspect calculation algorithms based on triangular facets have been studied for some time, but algorithms specific to TIN vertices have not been studied as much (Krcho 1992;van Kreveld 1997;Hu et al. 2021b).
In addition, LSCs calculation methods based on TINs, especially the methods for calculating gravityspecific curvatures, are still insufficient. van Kreveld (1997) used a curve fitting method to describe the plan and profile curvature calculations but did not perform further experiments. The computer graphics community has broadly explored gravity-invariant curvatures. The most recent study (Stupariu 2021) compared different gravity-invariant curvature calculation methods based on TINs derived from terrain point cloud data but did not consider gravity-specific curvatures. With different research objectives and disciplines, the computer graphics community does not focus on gravity-specific curvatures, although they are important in geoscience. Hence, LSCs calculation methods based on TINs are still insufficient, which limits the application of TIN models in geoscience.
Mathematical vector theory and the corresponding operations have been studied in geomorphometry and GIScience. Ritter (1987) proposed a vector-based slope and aspect algorithm, and then Hodgson and Gaile (1996) used a vector operation to generate the statistics of surface orientation. Li and Hodgson (2004) abstracted the vector field data model and operations for a raster data structure. The most recent studies concerning vector operations have mainly focused on second-order derivative calculations and land surface concavity-convexity quantification (Hu et al. 2020(Hu et al. , 2021a. These studies were based on Grid-DEMs. In this study, we propose a mathematical vector framework for LSCs calculation from TINs. In our framework, the different curvatures in the LSC system can be mutually derived. Then, our framework supports curvature calculations that involve both vertices and triangular facets. Finally, replaceable weighting methods are introduced into our framework. Our research is an extension of mathematical vector theory and vector operations in geomorphometry and GIScience.

Interpretations of LSCs
First, we provide a review of the interpretations of commonly used curvatures in geoscience. Considering a point P on an assumed smooth surface S, there is a tangent plane T P S and a corresponding normal vector n in space at that point (Figure 1(a)). An arbitrary plane passing through point P and the normal vector n will intersect the surface and form a curve. There are countless such curves, but two curves perpendicular to each other yield the maximum curvature K max and minimum curvature K min at point P (Do Carmo 1976). Similarly, this plane can also pass through the aspect vector and produce a curve, which is called the slope line in the literature (Figure 1(a)) (Minár, Evans, and Jenčo 2020;Hu et al. 2021a). The curvature of the slope line at point P is the profile curvature (K n ) s . Additionally, there is another curve perpendicular to the slope line, which is called the normal contour in the literature (Minár, Evans, and Jenčo 2020); this curve can be used to define the tangential curvature (K n ) c (Figure 1(a)). The plan curvature (also named contour curvature) (K p ) c is the curvature of the contour at point P and is not defined by using the normal vector n (Figure 1(a)) but can be derived from the relationship between (K n ) c and the land surface slope (Jochen, Evans, and Brinkmann 2003). The curvatures mentioned above are the most common and widely used in geoscience, and other curvatures can be derived based on them.

Normal curvatures and curvature tensor
From the perspective of differential geometry, bending at a surface in space is a self-dependent property that does not need to be represented and defined by curves. Considering point P again, the surface may exhibit different degrees of bending in different directions in space. Differential geometry uses normal curvatures to describe the spatial bending of the surface. Normal curvatures depend on the point location and the corresponding tangent vectors (i.e. the vectors in the tangent plane T P S) at the surface. K max ,K min , (K n ) s , and (K n ) c are surface normal curvatures and correspond to the maximal, minimal, slope, and normal contour directions, respectively (Figure 1(b)). These direction vectors are tangent vectors in the tangent plane T P S of point P. The plan curvature (K p ) c is not a normal curvature; thus, it expresses terrain plan bending rather than spatial bending. The normal curvature K n in the direction of the unit-length vector T can be calculated at a smooth surface S by using equation (1) as follows (Taubin 1995): where T = sT 1 + tT 2 is a tangent unit-length vector at point P on the surface, and {T 1 ,T 2 } is an orthonormal basis of the tangent plane T P S at point P on the surface. Note that II is the second fundamental form of surface S in differential geometry and is defined as a quadratic form with basis-invariance, i.e. the choice of basis is arbitrary (Do Carmo 1976). In equation (1), II is expressed as a symmetric matrix representation under the orthonormal basis {T 1 ,T 2 }, which is called the Weingarten matrix. In computer graphics, II is called the curvature tensor, and normal curvatures can be generated from it (Taubin 1995). In this study, we use the curvature tensor to calculate the LSCs based on TINs.

Land surface curvature calculation from TINs
Generally, the vertices of a TIN are irregularly distributed elevation sampling points on the land surface, and the facets of a TIN are abstractions and simplified representations of the local land surface. The TIN vertices and facets depict the land surface shape in a complementary manner, both are significant in geoscience. Thus, LSCs calculations should include both the vertices and facets of TINs. Inspired by the work of Rusinkiewicz (2004), who estimated the gravity-invariant curvatures in computer graphics, we propose a mathematical vector framework for LSCs calculations, including gravity-invariant and gravity-specific curvatures, based on both TIN vertices and triangular facets. It should be clarified that the definition of II and relative mathematical derivations hold only for smooth surfaces, however, as Rusinkiewicz (2004) pointed out, we can approximate them in the discretized case (such as a TIN) by using finite differences. The sub-sections below will provide the mathematical derivations and implemented finite differences in detail.

Calculation of the curvature tensor based on TIN triangular facets
The curvature tensor II (hereinafter referred to as II) for a triangular facet should be calculated first. This calculation is based on the relationship between II and the surface normal vector. Multiplying II by direction vector s in the tangent plane T P S yields the directional derivatives D s n of the surface normal vector at a smooth surface S (Rusinkiewicz 2004): There are three edge vectors s 1 ,s 2 , and s 3 for a given TIN triangular facet (Figure 2(a)). Based on the finite difference principle, D s n can be discretized as the difference in TIN vertex normal vectors (i.e. n 1 ,n 2 , and n 3 in Figure 2(a)) in the edge directions. The calculation of TIN vertex normal vectors involves a weighting method for normal vectors of triangular facets, and the method used in this study is based on the approach used by Max (1999) and Hu et al. (2021b). They used the area of each triangular facet divided by the square of the length of the two edges that touch the vertex. This weighting method is recommended but optional. Equation (2) should be implemented in the local coordinate system for each facet. The local coordinate system for each facet consists of u f = UnitVector(s 1 ), and v f = UnitVector (n f × u f ) (Figure 2(a)). Therefore, a set of linear simultaneous equations can be obtained based on equation (2) at a smooth surface S, as shown in equation (3) (Rusinkiewicz 2004): Equation (3) includes six equations and the three unknowns of e, f, and g; thus, it is an overdetermined set of equations that can be solved by using the leastsquares method. As Rusinkiewicz (2004) indicated, the finite difference approximation and least-squares approach provide a high degree of accuracy in many common cases. Thus, II expressed in the (u f ,v f ) coordinate system of a triangular facet can be calculated.

Calculation of the curvature tensor at TIN vertices
Assume that each vertex has a unique orthonormal coordinate system that consists of u p = UnitVector (s 1 × n p )), and v p = UnitVector (n p × u p ) (Figure 2(b)). In one case, the triangular facet normal vector n f equals the vertex normal vector n p , which means that the triangular facet coordinate system (u f ,v f ) and vertex coordinate system (u p ,v p ) are coplanar. Following equation (4) below, II can be expressed based on the directions of u p and v p at a smooth surface S (Rusinkiewicz 2004): where II is expressed in the triangular facet coordinate system (u f ,v f ). However, in most cases, (u f ,v f ) and (u p , v p ) are not coplanar. Accordingly, we should rotate the coordinate system of (u p ,v p ) to be coplanar with the coordinate system of (u f ,v f ). This rotation involves the cross product of vectors n f and n p and is based on equation (5) (Goldman 2011): À cos θÞ þ cos θ rxryð1 À cos θÞþrz sin θ rxrzð1 À cos θÞ À ry sin θ rxryð1 À cos θÞ À rz sin θ r 2 y ð1 À cos θÞ þ cos θ ryrzð1 À cos θÞþrx sin θ rxrzð1 À cos θÞþry sin θ ryrzð1 À cos θÞþrx sin θ r 2 z ð1 À cos θÞ þ cos θ 0 @ 1 A ; (5) where r equals the cross product of n f and n p ; θ is the angle between n f and n p ; and R is the rotation matrix. After applying this rotation to (u p ,v p ), (u p ,v p ) will be coplanar with (u f ,v f ). Thus, we can calculate II in terms of the corresponding (u p ,v p ) on each triangular facet by using equation (4), and the next step is weighting, i.e. how much the II on each triangular facet should be allocated at the corresponding vertices. Rusinkiewicz (2004) referenced the work of Meyer et al. (2003), who used "Voronoi area" weighting for triangles with various shapes and sizes. This weighting method is also optional.

Derivation of common land surface curvatures based on the curvature tensor
As mentioned in section 2.1.2, the maximum, minimum, profile and tangential curvatures are normal curvatures and can be derived directly from curvature tensor II. The other LSCs are basic curvature combinations or derivations. In this section, we provide operators for these commonly used curvatures. Equation (1) involves a case of f = 0 with a special orthonormal basis {T 1 ,T 2 }. In this case, the specific vectors {T 1 ,T 2 } are called principal directions of the surface at point P. The corresponding directional curvatures are the principal curvatures, i.e. the maximum and minimum curvatures. From differential geometry and linear algebra, the principal curvatures and directions can be calculated by the eigenvalues and corresponding eigenvectors of the II (Rusinkiewicz 2004).
Gravity-specific curvature calculations based on the curvature tensor have not been discussed in previous studies. Different from the directions of principal curvatures, the directions of the profile and tangential curvatures are specific. Thus, the objective is to express the slope and normal contour directions in the coordinate systems of (u f ,v f ) or (u p ,v p ). Assuming that a normal vector n = (n x ,n y ,n z ) could be a triangular facet or vertex vector, the corresponding aspect vector is n aspect = (n x ,n y , 0). Additionally, the normal contour direction T nc = UnitVector (n aspect × n), and the slope direction T ns = UnitVector (n × T nc ). Thus, the commonly used tangential, profile, and plan curvatures in geoscience can be obtained for a smooth surface S from equation (6) as follows: where (u, v) can be a triangular facet or a vertex coordinate system; the profile curvature is (K n ) s , the tangential curvature is (K n ) c , and the plan (contour) curvature is (K p ) c . The operators for LSCs calculations are oriented to curvature tensor II without TIN facet or vertex limitations. Accordingly, the LSCs of the TIN can be calculated based on II for a triangular facet or a vertex. This algorithm uses 1-ring of facets around each vertex to calculate the curvature tensor and LSCs, according to the pseudocode presented in Figure 3.

Derivation of other land surface curvatures
Previous studies have also proposed many other LSCs for practical application, although they are not used as commonly as the profile, tangential, and plan curvatures in geoscience. The relationships among these LSCs were studied by Shary (1995) and Minár, Evans, and Jenčo (2020). Table 1 consolidates the work of Minár, Evans, and Jenčo (2020). All these LSCs can be derived from the basic minimum K min , maximum K max , profile (K n ) s , and tangential (K n ) c curvatures. Thus, the proposed mathematical vector framework can support the relevant calculations for the LSC system.

Comparison methods
The comparison in this study focused primarily on grid-and TIN-based methods and attempted to assess the advantages and disadvantages of the proposed mathematical vector framework. For grid-based methods, the ZEVENBERGEN and EVANS methods were selected because they are widely used in land surface analysis. The ZEVENBERGEN method uses an incomplete quartic surface as the fitting function in a local 3 × 3 window, and the EVANS method uses a full quadratic surface. Details of the two methods can be found in previous studies (Evans 1980;Zevenbergen and Thorne 1987;Jochen, Evans, and Brinkmann 2003). Uniform and nonuniform sampling methods were used to build TINs and generate different triangle shapes. Uniform sampling for a TIN means isosceles right triangles (i.e. dividing grids by diagonals), while the shape of the triangles for random sampling is unconstrained. The influence of the sampling type on LSCs calculation based on TINs is explored in the next section. In short, for uniform sampling, the EVANS, ZEVENBERGEN and TIN-based methods can all be used because a Grid-DEM is actually a result of uniform sampling, and a TIN-based method can be used for nonuniform sampling. In addition, a TIN simplification method called quadric error metrics (QEM) is used in our research to produce the optimized TIN. Details about the QEM method can be found in the related works (Garland and Heckbert 1997;Feciskanin 2012). The influence of the TIN triangle shape will be explored.
Because the mathematically simulated surface (see the materials section for details) can provide analytical results, accuracy assessments were based on this surface for the ZEVENBERGEN, EVANS, and TIN-based (including uniform sampling, random sampling, and QEM simplification) methods. In the ranges of x and y given above for the mathematically simulated surface, we varied the cell size from 6 m to 20 m at an interval of 2 m to generate a series of Grid-DEMs, and these uniform points were also used to construct the TINs. In random sampling, the number of sampling points in each cell size was the same as that for the i.e. with 40,401,22,801,14,641,10,201,7,396,5,776,4,489 and 3,721 random sampling points. With the create random points tool in ArcGIS, a series of TINs was built by using the Delaunay triangulation method based on random points on the mathematically simulated surface. For the QEM simplification method, the initial TINs are constructed by uniform and random sampling with a 5 m cell size number of points respectively. Then, a series of simplified and optimized TINs are generated by using the QEM method, and the same number of sampling points is kept as above. The different cell sizes and sampling densities represent different scales of the mathematically simulated surface, and this approach can help to explore the scale effects on the calculation results in our research. Moreover, error sensitivity analysis is also an important part of the comparison of different methods. We can use the simulation error for the mathematically simulated surface to conduct this comparison experiment. The parameters of the normally distributed error are the mean and standard deviation, and the standard deviation is related to the error intensity. We set the standard deviation from 0.05 to 0.5 in 0.05 intervals to generate a series of different intensities of normally distributed error and then exert it to the cell size 6 m Grid-DEM, uniform TIN, and random TIN and It describes the bending degree of a surface in the Euclidean space (Gauss 1828) and can quantify the divergence of the sediment flow rate. Gaussian K g = K max * K min It is often used in geology and cartography due to the invariance of a curve length with zero Kg, and can identify specific conical or cylindrical landforms. Unsphericity K u = (K max -K min )/2 It quantitatively describes the proximity between the surface and the sphere (Shary 1995). Casorati K c = sqrt(((K max ) 2 + (K min ) 2 )/2) It expresses the magnitude of surface bending regardless of its shape (Casorati 1900;Florinsky 2017). Total accumulation K a = (K n ) s * (K n ) c It reflects the relative accumulation area of surface material migration (Shary 1995).
It expresses an excess of mass flow energy that can be directly proportional to denudation (and inversely proportional to accumulation) (Minár, Evans, and Jenčo 2020). Horizontal excess It describes the extent to which tangential curvature is larger than minimum curvature (Shary, Sharaya, and Mitusov 2002).
It describes the extent to which profile curvature is larger than minimum curvature (Shary, Sharaya, and Mitusov 2002). Total ring It is a transformation (square) of the contour geodesic torsion that is one of the basic trios that expresses a twisting (gravity discordance) of the land surface (Minár, Evans, and Jenčo 2020). Cross sectional* Z cc = (K n ) c /cos(S) It is the rate of slope change in the direction perpendicular to the downslope direction (Wood 1996). Longitudinal* Z ss = (K n ) s /(cos 3 (S)) It is the rate of slope change along the downslope direction (Wood 1996). Flow line (streamline)* |(K p ) s | = sqrt(K r * (1 + tan 2 (S))/ tan 2 (S)) It is a directional derivative for characterizing a change of aspect in the direction of a slope line and is used to describe the swing degree of the flow path (Krcho 1992).
* S = arctan(sqrt(n x 2 + n y 2 )/n z ) is the slope; the cross-sectional and longitudinal curvatures (Z cc and Z ss ) are the second-order derivatives of elevation and are associated with curvatures as their slope dependent subforms (Minár, Evans, and Jenčo 2020).
simplification TIN with the same number of points. The simplification TINs were generated from the uniform and random sampling TINs with the number of points of the 5 m cell size. In this way, we assessed the robustness of different LSCs calculation methods based on different DEM structures. The mean absolute error (MAE) concerning LSCs between the calculated and analytical results was used for both the scale and error effect experiments mentioned above based on the mathematically simulated surface (Chai and Draxler 2014). We also drew a line chart of the MAE values for different methods and types of TINs to exhibit the trends of the scale and error effects. Note that the boundary points of Grid-DEMs or TINs have been excluded in comparison statistics.
The first case study area is the Jiuyuangou watershed ( Figure 5(a)), a typical loess hill-gully landform located in northern Shaanxi Province. The topographyof this area is fragmented because of severe gully erosion. The open NASADEM dataset (NASA 2020) with 1 arc-second spacing (approximately 30 m) in this area was used. The compound method (Zhou and Chen 2011;Chang 2019) that integrated the maximum z-tolerance and the river network ( Figure 5(a)) was used to generate a hydrologically constrained TIN based on a Grid-DEM. In this research, we set a tolerance value of 5 m in the maximum z-tolerance method.
The second case study area is the Liujiaping gully ( Figure 5(b)), which is the orographically right branch of the Jiuyuangou watershed ( Figure 5(a)). The DEM data for this area were provided by the Shaanxi Bureau of Surveying and Mapping including an original contour map with 1:10000 plotting scale, a TIN converted from the original contours ( Figure 5(b)), and a Grid-DEM 5 m cell size converted from the TIN by using linear interpolation.
The third case study area is a small watershed ( Figure 5(c)), called Qiaogou, which is located in the Liujiaping gully ( Figure 5(b)). The elevation information for this area was based on point cloud data produced by unmanned aerial vehicle (UAV) acquisition and photogrammetric techniques (Figure 5(c)). After downsampling, the original density point cloud was converted into a sparse point cloud with a minimum spatial distance of 1 m. This sparse point cloud was directly converted into a TIN by using the Delaunay triangulation method. In addition, the dense point cloud was converted to a Grid-DEM with a 1 m cell size for comparative experiments.
The basic information for the three real landforms is provided in Table 2. We matched the appropriate data sources at different scales, and all the data sources can be converted into TINs. The methods of calculating LSCs based on the Grid-DEMs and TINs in the three different areas were further investigated.

Results
The scale and error effect experiments performed to calculate the profile and tangential curvatures were based on the mathematically simulated surface. As noted in the above sections, these two gravityspecific curvatures are fundamental to the LSCs system and are widely used in geoscience. In this section, we show the advantages and limitations of LSCs calculations based on TINs under the proposed mathematical vector framework. In the terms of Minár et al. 2013, the first experiment (section 4.2) reveals differences in the method error (computational accuracy for errorless DEM) in various scales. The second one (section 4.3) shows how compared methods deal with data error (the error rate of DEM) that dominate in our example on the cell sizes of 6 m and coarser. A balanced influence of method error and data error (using finer cell size and different DEM data sources) could be further explored.

Profile and tangential curvature results based on TINs
Based on the mathematically simulated surface with 1,442,401 random sampling points (which corresponds to uniform sampling with an interval of 1 m), the frequency curves of the profile and tangential curvatures   show that the changes are relatively gradual in most areas of the surface because the corresponding absolute curvature values are less than 0.01 (Figure 6(a,b)). This range of curvature values is similar to that for the Jiuyuangou watershed (Figure 6(c,d)), but the resolution of the data source used for Jiuyuangou is approximately 30 m. As a subregion of the Jiuyuangou area, the Liujiaping gully displayed an absolute curvature value of less than 0.05 in most areas based on the curvature frequency curves (Figure 6(e,f)). The data source for the Liujiaping gully is a contour map with a 1:10000 plotting scale, which is equal to an approximately 5 m resolution in the Grid-DEM. The most special area is the Qiaogou subregion in the Liujiaping gully, where point cloud data were used. The frequency curve results indicate that the absolute curvature value was less than 0.5 in most parts of Qiaogou ( Figure 6(g,h)), which reflects a complex land surface structure. The more detailed DEM can capture more subtle forms in the framework of the nested hierarchy of landforms (Minár and Evans 2008;Evans 2012), which results in higher curvature values. Thus, the results show that the LSCs are very sensitive to the scale of the analysis, defined by the resolution of the DEM. In addition, the results actually are influenced by the accuracy issues caused by various DEM sources and different TIN conversion methods.

Assessment of method error and scale effects
By using the analytical results of the mathematically simulated surface as reference data, we assessed the scale effects on the curvature calculation results based on the fluctuations in method error. Tables 3 and 4 show that the accuracy results of the ZEVENBERGEN method were better than those of the EVANS method based on Grid-DEMs. This is a reasonable observation as in previous studies (Jochen, Evans, and Brinkmann 2003), due to the ZEVENBERGEN polynomial fitting more precise input points than the EVANS polynomial. The two grid-based methods provide more accurate profile and tangential curvature results than the results of the TIN method based on four types of TINs. For the results of the TIN-based method, TINs with uniform sampling show the best accuracy, TINs with QEM simplification take second place, the TINs with random sampling are the worst. Specific to the results of QEM simplification, TINs simplified from uniform sampling perform better than TINs simplified from random sampling. Figure 7 shows the line charts of the MAE (equal to the method error in this experiment) as the number of sampling points decreases based on Tables 3 and 4, where similar trends of the six curves imply that the TIN-and grid-based methods have similar scale dependence. As a previous study pointed out (Minár et al. 2013), in addition to being based on the Grid-DEM, the power functions can also approximate the dependence of the method error (MAE in Tables 3 and 4) on the number of sampling points used in the TIN. Nevertheless, the LSCs calculation based on TIN is still of significance because method error is in practice usually less important than data error. Moreover, many geographical processes and phenomena need to be described and quantified at different scales, and the TIN structure can support multiple-scale terrain sampling. Method error of the TIN-based method is greatly affected by the shape of TIN triangles. Obviously, the TIN optimization used by QEM does not optimally play with the curvature calculation  * ⅰ, ⅱ, ⅲ, ⅳ, Ⅴ, and VI represent the ZEVENBERGEN, EVANS, TIN-based (uniform sampling), TIN-based (random sampling), TIN-based (QEM simplification for uniform sampling), and TIN-based (QEM simplification for random sampling) methods, respectively.  presented here. More analysis concerning the influence of the shape of TIN triangles can be found in the discussion section.

Assessment of the data error effects
The cell size of Grid-DEM and uniform sampling TIN is set as 6 m and kept invariant in this experiment. Accordingly, the vertex number of the random sampling TIN and simplification TINs is the same as that of a 6 m cell size Grid-DEM. The standard deviations of normally distributed error are increasing and its relation to the cell size leads to the high dominance of data error over method error in the MAE values. Tables 5 and 6 show that the MAE values of the results of the ZEVENBERGEN and EVANS methods based on Grid-DEMs increase rapidly as the normally distributed errors increase, while the MAE values of the TIN-based method with uniform sampling, random sampling, and QEM simplification remain low. The observation that the EVANS method provides lower error sensitivity than the ZEVENBERGEN method has been reported in a previous study (Schmidt, Evans, and Brinkmann 2003). The TIN-based method displays the best capacity for error resistance on the uniform sampling TIN, a slightly worse capacity on the QEM simplification TIN, and the worst capacity on the random sampling TIN. Specific to the results of QEM simplification, TINs simplified from uniform sampling show better characteristics than TINs simplified from random sampling. Figure 8 shows the line charts of the MAE as the normally distributed error intensity increases based on Tables 5 and  6, where the trends of error sensitivity for different methods and structures are obviously displayed. The curves for the uniform sampling TIN and QEM simplification TINs based on uniform and random sampling are similar, and they have low MAE values, which shows that the optimized TIN also has advantages in LSCs calculation over the random sampling TIN on DEMs with high error. In addition, due to the fact that large cell size will also weaken the accuracy of the LSCs calculation, the error sensitivity of different methods is expected to be different when the cell size of DEM becomes finer. Generally, high accuracy is the basic requirement for any method, but typically, the data quality must be high. However, measurement error cannot be avoided, especially in regions with complex terrain, and error-free data are impossible to obtain (Schmidt, Evans, and Brinkmann 2003). Hence, the error sensitivity of a method has become increasingly important in many applications. In the following section, we demonstrate the application potential of the proposed TIN-based LSCs calculation method for real landforms.
sampling TIN are equal. Hence, more suitable weighting methods will help to improve the calculation accuracy based on irregular triangles. On the other hand, although the isosceles right triangles (i.e. dividing grids by diagonals) show advantages in the LSCs calculation, it is still not clear which shape of triangles is optimal for the weighted method used in this paper. Perhaps equilateral triangles converted from regular hexagonal grids can be yet more suitable, and further research should be conducted. In traditional grid-based land surface analysis, the analysis extent is expanded usually experimentally (i.e, by employing a larger window) to achieve generalized results (Florinsky 2009). This situation is different in TIN-based analysis because the scales of terrain expression can be adjusted. In this way, terrain analysis can focus on a 1-ring of facets around each vertex or only a triangle facet, if a representative and optimized TIN for a certain scale of the land surface is available. As Feciskanin and Minár (2021) pointed out, TIN multiple-scale simplification for terrain analysis is becoming increasingly important. The representative triangles can be long and narrow, however, the TIN optimization aims to avoid it. From this view, developing a weighting method less sensitive to triangle shapes in the curvature computation process can be helpful to land surface analysis.

Preliminary case studies for geomorphological classification based on TINs
Application in the classification of land surface concavity-convexity.
In geoscience, land surface variations can be classified into concave and convex units. These two types of units control the direction of flow, the transport of materials, and the accumulation of sediment (Li et al. 2020). The degree of land surface concavity-convexity can be directly quantified by using the mean curvature (Romstad and Etzelmüller 2012), but we focus only on concave and convex classifications in this section to demonstrate the application of the TINbased method. We simply define the areas with K mean > 0 as convex units and the areas with K mean < 0 as concave. The narrow band at approximately 0 (i.e. planar units) is not defined because its range is difficult to determine. The traditional processes are to use membership functions to avoid crisp values (Schmidt and Hewitt 2004;Qin et al. 2009). In addition, due to its low error sensitivity, the EVANS method is used for Grid-DEMs in the three areas and compared with the proposed TIN-based method. At different scales, concave and convex terrain units have different geographical meanings. At the scale of the Jiuyuangou watershed based on NASADEM data, the concave units are usually valley bottoms, and the convex units are hills and ridges. Although the EVANS method provides low error sensitivity, the result is still greatly affected by the quality of the open DEM (Figure 9(a)). In contrast, the TIN-based method provides a reasonable hill and gully distribution pattern and ensures the continuity of ridges and hydrological networks (Figure 9(b)). At the scale of the Liujiaping gully based on a 1:10000 plotting scale contour map, the concave units involve gullies of different sizes and some concave lower hill slopes, and the convex units are mainly peaks, ridges, and some convex upper hill slopes. The Grid-DEM for the Liujiaping gully is converted from a TIN by using an interpolation method.
This process increases the uncertainty of the Grid-DEM data (Gosciewski 2013), which is reflected in the fragmented valleys and ridges in the results of the EVANS method (Figure 9(c)). The TIN-based method avoids the interpolation process and performs analyses directly with the TIN. Consequently, the classification results generally match the actual terrain features of the Liujiaping gully and are better than those of the EVANS method (Figure 9(d)). At the scale of the Qiaogou small watershed based on a point cloud data source, the concave units are road cuts, agricultural terrace surfaces, and rills on hill slopes; the convex units are mainly hilltops, scarps, and agricultural terrace ridges. In theory, based on the high-resolution Grid-DEM data derived from dense point clouds, the EVANS method can provide a reliable classification result. However, as mentioned in the above sections, measurement errors cannot be entirely avoided. The classification results of the EVANS method are not visually superior to those of the TIN-based method (Figure 9(e,f)). Furthermore, this classification difference can also be due to the different degrees of generalization between the structures of Grid-DEM and TIN.
Application in the classification of hillslope units. In this case, we perform a preliminary experiment that involves the classification of hillslope units based on the TIN for the Liujiaping gully. The principle of classification references the work of Drăguţ and Blaschke (2006) and is simplified in this research. The slope gradient is first used to classify the terrain as flat (slope gradient < 2°), hillslope (2° ≤ slope gradient ≤ 45°), or steep slope (slope gradient > 45°). Then, to simplify the problem, we use the positive or negative profile and tangential curvatures to classify the hillslope shapes as nose slopes, negative contacts, head slopes, and shoulders. The details of these principles are shown in Table 7. With this classification method, the results based on the Grid-DEM (Figure 10(a)) remain fragmented as observed for the terrain concavity-convexity classification result (Figure 9(c,d)). This result does not support people's basic understanding of hillslope units. In contrast, the classification result obtained with the TIN-based method (Figure 10(b)) based on the TIN is generally reasonable. Figure 10(c) shows the classification area difference based on a TIN and Grid-DEM for various hillslope units. Note that this classification method that uses the crisp values of terrain attributes has potential limitations because the thresholds of the slope and LSCs are scale dependent. In fact, the classification of hillslope units is complex. The crisp values of terrain attributes are usually not sufficiently representative and are replaced with membership functions (Schmidt and Hewitt 2004;Qin et al. 2009). In addition,  the object-based image analysis (OBIA) technique has been used to segment terrain into the analysis units for further processing. However, multiresolution segmentation is usually uncontrolled and relies on parameter adjustment. The triangular facets that consist of vertices are objects, and the Voronoi polygons in which the sampling points are located also are objects. The shape, size, and position of each object are controllable. The proposed TIN-based LSC calculation method can directly calculate the basic terrain attributes for triangular facets and Voronoi polygon objects. This approach presents a new opportunity for the classification of hillslope units in the future.

Implications of the mathematical vector for LSCs definition and calculation in geoscience
LSCs are widely used in natural hazard susceptibility assessments, natural potential evaluation, the mapping of landforms, land surface segmentation and classification, the processing of LIDAR data, topographic visualization, and other geoscience tasks (Jasiewicz and Stepinski 2013;Garosi et al. 2018;Khosravi et al. 2018;Chen et al. 2018;Rahmati, Reza Pourghasemi, and Melesse 2016;Chen et al. 2017). However, as Minár, Evans, and Jenčo (2020) noted, the use of LSCs in most applications is driven by GIS software capabilities but not supported by theory, which may lead to confusion in concept and definition. On the one hand, the lack of understanding can be attributed to the diverse handling and interpretation of LSCs in GIS software (Schmidt, Evans, and Brinkmann 2003). On the other hand, the LSCs system is both theoretically and mathematically complex (Florinsky 2017). Although a comprehensive definition of the LSCs system has been established in geoscience (Evans 1980;Shary 1995;Shary, Sharaya, and Mitusov 2002;Minár, Evans, and Jenčo 2020), the corresponding theory still does not involve LSC calculations based on non-Grid-DEMs. As stated in the introduction section, the applications of TINs in geoscience become increasingly common, which means that the LSCs calculation method based on TINs is urgent and necessary. Our proposed framework re-explores LSCs from the perspective of mathematical vector theory, and uses the concept of curvature tensor to re-explain and recalculate the commonly used LSCs. Our research can improve the general understanding of the LSCs system and support LSCs calculations based on TINs, point clouds or other structures. Our proposed framework is not only an improvement on and supplement to the current LSCs system but also an advancement in land surface analysis based on the structure of TINs. Due to the widespread lack of analytical methods for this data model, the application potential of TINs in geoscience research has been underestimated. Based on the theory of mathematical vectors, we preliminarily explore the application of TINs in land surface analysis. We believe that this approach will provide new ideas for other geoscience studies.

Conclusion
Considering the research gap regarding LSC calculations based on TINs, in this study, we propose a mathematical vector framework to enhance LSC system theory. In this framework, LSC can be calculated based on both triangular facets and vertices, and the selection of weighting methods in the framework is flexible. We use the concept of the curvature tensor to interpret and calculate the commonly used LSC, which provides a new perspective in geoscience research. We also investigate the capacity of the TIN-based method to perform LSCs calculations and compare it with grid-based methods. Based on a mathematically simulated surface, we reach the following conclusions. First, the TIN-based method has similar scale effects to the grid-based methods of EVANS and ZEVENBERGEN. Second, the TIN-based method is less error sensitive than the grid-based methods by the EVANS and ZEVENBERGEN polynomials for the high error DEMs. Third, the shape of the TIN triangles exerts a great influence on the LSCs calculation, which means that the accuracy of LSCs calculation can be further improved for optimized TIN using weighting methods less sensitive to the shape of the TIN. Based on three real landform units and respective different data sources, we discuss the possible applications of the TIN-based method, e.g. the classification of land surface concavity-convexity and hillslope units. We find that the TIN-based method can produce structurally and visually better classification results than the grid-based method. This qualitative comparison reflects the potential of using TINs in multiscale geoscience research and the capacity of the proposed TIN-based LSC calculation methods. Our proposed mathematical vector framework for LSCs calculations from TINs is a preliminary approach to mitigate the multiple-scale problem in geoscience. In addition, this research integrates mathematical vector and geographic theories and provides an important reference for geoscience research.