Reconstruction of sediment-transport pathways on a modern microtidal coast by a new grain-size trend analysis method

A new method for granulometric-parameter-based reconstruction of sediment-transport pathways is proposed and is termed P-GSTA (grain-size trend analysis using principal component analysis) herein. The main advantage of this method is its applicability to depositional environments involving mixed transport processes, for instance, fluvial, tidal, and wave-influenced environments. In the P-GSTA method, a linear function with all significant granulometric parameters that are summed with different weighting factors was used to infer sediment-transport direction (sediment flux pattern); the previous grain-size trend analysis (GSTA) methods considered only three parameters (mean grain size, sorting, and skewness) with equal weighting. This study chose six parameters (namely, median grain size, coefficient of variation, skewness, kurtosis, and mud and gravel log-ratios) for calculation. First, the zero values of mud and gravel fractions are replaced, and their log-ratios are defined. Then, all values are standardized. Thereafter, principal component analysis (PCA) is conducted to determine the weighting factor of each granulometric parameter. Each principal component is then interpreted, and the function best representing a sediment flux pattern is chosen from these components. Trend vectors are calculated, solely on the basis of a map interpolated from the scores of the chosen principal component, as the two-dimensional gradient of this value. The P-GSTA method proposed in this study was applied to a modern microtidal coast (tidal sand flat along the Kushida River Delta, central Japan). Sediment-transport pathways reconstructed by this method were consistent with observed sediment-transport patterns determined by field experiments using tracer sediments and geomorphologic observation; the results of the previous GSTA method were inconsistent with the observations. The proposed method also revealed additional minor depositional processes on the sand flat, namely, the deposition of fluvial-channel lags and muddy particles. Thus, this study demonstrates that the proposed P-GSTA method is a potentially powerful tool to reconstruct sediment-transport patterns even under mixed transport processes, where the estimation of the sediment-transport function is difficult.


Introduction
Depositional systems are consequences of sediment dispersal from their provenances. Progressive and selective transport of grains differentiates lithofacies deposited in different sedimentary environments (Swift et al. 2003). In particular, in coastal environments, sediment supply and spatial changes in dominant sediment transport processes (wave, fluvial, or tidal currents) play a major role in producing a unique pattern of shoreline geomorphology, offshore bathymetry, and sediment distribution patterns, in addition to the geomorphological geometry of the depositional system and the type of the depositional environment (Tanner 1962;Liu et al. 2000). The size-selective transport process of sediments results in downcurrent variation of the sediment grain-size distribution. Therefore, spatial variation in the grain size distribution of bottom sediments can be an important clue in reconstruction of sediment-transport pathways, which is important in coastal engineering and environmental protection.
Many researchers have attempted to identify spatial variations of granulometric parameters, which have been referred to as 'grain-size trends' , associated with net sedimenttransport pathways (McCave 1978;McLaren 1981;Le Roux and Rojas 2007). This is because methods to reconstruct sediment-transport pathways from granulometric patterns are inexpensive and easily available compared to direct measurements of sediment movement. In addition, there is a possibility that ancient sediment-transport patterns can be reconstructed by applying these methods to subsurface sediments. Early studies suggested that sediments become finer in the direction of transport (e.g., Pettijohn and Ridge 1932;Self 1977), while other investigations have indicated the opposite trend (e.g., McCave 1978;Nordstrom 1981Nordstrom , 1989. In order to address this disparity, McLaren and Bowles (1985) discussed the combination of three parameters (mean grain size, sorting, and skewness) on the basis of flume experiments and statistical calculations. They concluded that sediments always become better sorted during transport, whereas the mean grain size and skewness vary with specific combination patterns (FB−: finer mean grain size, better sorting value, and greater negative skewness; or CB+: coarser mean grain size, better sorting value, and greater positive skewness). Gao and Collins (1992) applied their theory to twodimensional data of grain-size distribution in an inner bay environment. In their method, after comparing three granulometric parameters (mean, sorting, and skewness values) of adjacent points, vectors of unit length (100 m) are drawn between two points if they satisfy the conditions (FB− or CB+). These vectors are calculated from every point with respect to the eight immediately neighboring grid points. Summing the vectors at each sample point produces a single-trend vector of sediment transport pathways. This approach is called grain-size trend analysis (GSTA) and has been applied in a variety of depositional environments (e.g., Pedreros et al. 1996;Cheng et al. 2004).
It has been reported, however, that there may be serious discrepancies between the GSTA model and actual sediment-transport patterns (Asselman 1999;Masselink et al. 2008). This is probably because grainsize trends do not always follow the assumptions of McLaren and Bowles (1985), especially under mixed transport processes. Grain-size trends due to sediment transport can vary due to the depositional processes, environments, and sediment provenances. Deviation from the presumed trend of even a single parameter is sufficient to result in the GSTA not accurately reflecting the trend vectors. Indeed, Flemming (1988) disputed whether sorting value would always improve during transport. In addition, it seems problematic that the GSTA model considers only three equally weighted parameters, but completely ignores other parameters representing grain-size distribution, such as kurtosis and mud/gravel fractions. In order to address these problems, Rojas et al. (2000), Rojas (2003), and Rojas and Le Roux (2003) proposed a model that combines four grain-size parameters (mean grain size, sorting, skewness, and kurtosis) into a linear granulometric facies function with differential weighting. These weighting factors are modified incrementally and iteratively for comparison with actual sediment-transport directions measured from bedforms or visual observations, so that the error between the observed vector and the generated vector is minimized. This method can determine the optimum weights of granulometric parameters, but requires detailed observational data of existing sedimenttransport patterns, and thus vitiates the abovementioned advantages of grain-size trend analysis.
To this end, the present study employs principal component analysis (PCA) (Davis 1986;Middleton 2000) for the purpose of treating multiple grain-size parameters with appropriate weightings. PCA is a multivariate technique for explaining the correlation between explanatory variables and automatically organizing these parameters into a few linear synthetic variables with different weights. It is reasonable to regard this spatial variation of grain-size parameters as resulting from the sedimenttransport processes. Therefore, it can be expected that parameters that strongly reflect sediment transport are emphasized by PCA. Hereinafter, this new method of GSTA employing PCA is referred to as P-GSTA.
Previous GSTA methods also include problems in validation of the methods, as sediment-transport trends vary according to the time scale of the observation. Sedimenttransport pathways of daily processes such as fair-weather waves can counteract those of monthly or yearly transport processes such as storms or slope failure events. Thus, the time scale of the field observation must be specified in order to validate the method for reconstructing sedimenttransport pathways. The time scale of sediment-transport patterns estimated by any GSTA method is defined by the sampling interval (Gao and Collins 1992). Subsurface sediments in shallower sampling intervals can be entrained and transported by waves or currents, whereas deeper layers are "historical" sediments that were already fixed and are not entrained by the surface processes. Nevertheless, this problem has been given little attention.
To reconstruct the ongoing sediment-transport patterns, a sampling depth should be assigned to the thickness of the active layer of the bottom-surface sediments. The active layer is a concept proposed by Hirano (1971), indicating the uppermost layer of the bottom-surface sediments, which are frequently reworked and deposited by currents and wave processes. Thus, grain-size variation of bottomsurface sediments accompanied by sediment transport occurs within this layer (Hirano 1971). It should be noted that the thickness of the active layer varies with time under non-uniform flow conditions in natural environments; it generally becomes thicker in accordance with the time development of the probability distribution of burial and reworking as functions of depth below the bed surface (Ribberink 1987;Armanini 1995;Parker et al. 2000).
A microtidal sand flat along the Kushida River Delta, Ise Bay, central Japan, was chosen as the study area. This microtidal sand flat is a largely natural environment, lacking breakwaters or dams, and is therefore suitable for observation of sediment transport. The active-layer thickness and net sediment-transport directions over a specific time scale were then measured by tracer experiments and geomorphological observation. Finally, sampling for the GSTA method (the proposed P-GSTA and previous GSTA) was conducted according to the active-layer thickness, and the results were validated with net sedimenttransport patterns recorded by field measurements.
The proposed P-GSTA method successfully reproduced the actual sediment-transport patterns, whereas the results of the previous GSTA model were inconsistent. In addition, application of the PCA model revealed details of previously unrecognized depositional processes on the sand flat. To summarize, this method is flexible because it requires no assumption on the sediment transportation processes and the influence of sediment transportation is estimated only from the characteristics of the obtained data. Thus, the method is effective under mixed transport processes, in which estimation of the sediment-transport function is difficult.

Methods/Experimental
Granulometric parameters used for P-GSTA method On the basis of grain-size measurements of samples, each of the granulometric parameters (mean grain size μ, sorting σ, skewness S k , and kurtosis k t ) is generated by the moment method (Folk 1966;Harrington 1967), as follows: where d is the representative value of each grain-size class (in increments of 0.1φ) and p is the weight fraction.
Physical characteristics of mud and gravel are evaluated as weight percentages (mud and gravel fractions: F mud and F gravel , respectively). The coefficient of variation (CV) is also employed to represent the degree of sediment sorting. CV is the dimensionless sorting value divided by mean value, calculated as follows: The P-GSTA method employs six parameters (median grain size d 50 , CV, S k , k t , and fractions of mud and gravel F mud and F gravel ) in a sediment-transport function. Here, the median grain size d 50 is employed rather than mean grain size, because median grain size is mathematically independent of the other granulometric parameters, whereas mean grain size is used in the calculation of CV, S k , and k t . In addition, CV value is employed here instead of sorting value (standard deviation) σ. Compared to sorting value σ, CV value has the advantage of avoiding mathematical correlation with mean value. When comparing grain-size distributions whose mean values are nearly identical, sorting values are a good index of the dispersion of the grain-size distribution. However, sorting values cannot be used when comparing grainsize distributions with quite different means because sorting value (standard deviations) has an inevitable positive correlation with the mean value of the grain size distributions under consideration. In such cases, CV rather than sorting value needs to be used in general sediment evaluation.

Parametric interpolation by kriging
Kriging interpolation is employed to estimate spatial variations of granulometric parameters Webster 1980a, 1980b;Kohsaka 1999). All GSTA methods have a problem in that the estimated trend vectors are potentially affected by the sampling intervals. Asselman (1999) solved this problem using kriging to interpolate observational data sets Webster 1980a, 1980b;Kohsaka 1999). Kriging is an algorithm based on least squares, which is used to estimate the spatial variation of a real-valued function, and assumes that spatial variation can be estimated from a linear combination of measured values (Kohsaka 1999). Weighting coefficients are obtained on the basis of the spatial dependence of a variable, which can be represented by a semivariogram, i.e., the scatter diagram of covariance with respect to spatial distance. The weighting of the running average is then determined by the variogram model function, which is the function fitted to the semivariogram (Burgess and Webster 1980a). Theoretically, the covariance of spatial data increases with distance and becomes a steady value over some distance Webster 1980a, 1980b;Kohsaka 1999). This limited distance is the range within which the data indicate spatial dependence. The spatial dependencies of the data are also shown in the variances between the estimations, which provide a measure of uncertainty in the interpolation values.

Data preparation
Prior to statistical analysis, granulometric parameters should be converted into appropriate forms for analysis by conventional statistical techniques.
It is common to report grain-size fractions in percentages, but such data are not independent, because they must sum to 100%. Compositional data such as mud and gravel fractions are therefore difficult to use in statistical analyses as they are (Aitchison 1986(Aitchison , 2003Aitchison and Egozcue 2005). Aitchison (1986Aitchison ( , 2003 proposed, instead, to evaluate the relative magnitudes of constituents in comparison to other compositional variables and to use logarithms to simplify the computation of the variance and covariance. After taking logarithms, the resulting transformed values are free to vary over the entire range of real numbers. However, a zero value in the compositional data prohibits log-ratio transformation, because the transformed value is −∞, and therefore cannot be included in the statistical analysis. In order to avoid this problem, a multiplicative approach (Martín-Fernández et al. 2003) replaces a zero value with a value that is sufficiently small, as follows: where F j andF j denote the original and substitute compositional values respectively; δ j , the input value of F j ; and c, the sum-constraint constant (100%). The modification of nonzero values is multiplicative, as suggested by the simple replacement rule. Generally, δ j is the maximum rounding error, but the range δ γ /5 ≤ δ j ≤ 2δ γ , where δ γ is the maximum rounding error, has been suggested to be reasonable (Aitchison 1986;Martín-Fernández et al. 2003). This nonparametric replacement is stable, i.e., independent of both data numbers and component numbers, and retains the original information of the compositional data (Aitchison 1986;Arai and Ohta 2006). In the present study, mud/gravel fractions were treated via these two transformational processes. First, zero values of mud/gravel fractions were replaced by the maximum rounding error (0.01%). Secondly, the dimensionless magnitudes of the mud and gravel fractions were defined as the mud and gravel log-ratios r mud and r gravel , as follows: whereF mud andF gravel are the substitute values for mud and gravel fractions (Eq. 6), respectively.
The units, mean values, and standard deviation of variables strongly influence the results of PCA. To avoid this problem, the following formula was used to standardize all data: where m x is mean and s x is the standard deviation of each granulometric parameter x. This operation transforms the mean and standard deviation of each variable into 1 and 0, respectively.

Weighting (principal component analysis)
The new method for defining the linear granulometric function proposed in this study is based on principal component analysis (PCA), which is a technique for explaining the correlation between multiple variables and accounting for a large variance using only a few synthetic variables (principal components) (Davis 1986;Middleton 2000).
Let x 1 , x 2 , ⋯, x p be variables, and let z 1 , z 2 , ⋯, z m (m ≤ p) be their respective weighted synthetic variables. These synthetic variables are defined as the principal components.
Here, the matrix (L i1 ...L ip ) is the eigenvector of the ith principal component. The m synthetic variables have the following two properties: (1) the correlation of any two synthetic variables is zero (orthogonal) and (2) the following inequality holds for the variances of synthetic variables Var(z i ): Var(z 1 ) ≥ Var(z 2 ) ≥ ... Var(z m ). The amplitude of the variance of principal components can be regarded as the amount of information in each principal component. The factor loading, which represents the relative importance of each variable in the principal component, is defined as the product of the root of the eigenvalue and the eigenvector. The absolute value of the factor loading is proportional to its significance.

Calculation of trend vectors of sediment transport by P-GSTA method
The procedures to obtain trend vectors of sediment transport by GSTA with principal component analysis (PCA) weighting (P-GSTA) are described herein. As discussed above, this method employs six granulometric parameters (median grain size, CV, skewness, kurtosis, and mud and gravel fractions). First, the zero values of mud and gravel fractions are replaced (Eq. 6), and mud and gravel log-ratios are defined (Eqs. 7 and 8). Then, all values are standardized (Eq. 9). After these preliminary procedures, PCA is conducted using these six transformed parameters (Eq. 10). Each principal component is then interpreted, and the first principal component is chosen as the sediment-transport function T(x, y). The trend vectors V(x, y) are then calculated, solely on the basis of the interpolated map of the scores of the chosen principal component as the two-dimensional gradient of this value, as follows: where V(x, y) is an estimated trend vector of sediment transport, and i and j are unit vectors (magnitude: 100 m) in the x and y directions, respectively.

Application of P-GSTA method to microtidal sand flat
We applied our method to a microtidal sand flat located at the mouth of the Kushida River. The headwaters of the Kushida River are in the Daikou Mountains of Japan; the river flows eastward to the lowlands and then northward, to form a delta at Ise Bay (Fig. 1). The catchment area of the river is approximately 461 km 2 . According to the Japanese Ministry of Land, Infrastructure, and Transportation (MLIT), the annual average fluvial discharge was 11.56 m 3 /s in 2010 (http://www1.river.go.jp/). The delta is lobate or arcuate and is classified as fluvialdominated or complex fluvial-wave-tide dominated on the basis of the classification of Galloway (1975).
The tidal range monitored on the left bank of the Kushida River (36°36′ 0″ N, 136°33′ 45″ E) is approximately 2 m during a spring tide (microtidal). Wave heights have been measured at the river mouth (34°22′ 16″ N, 137°07′ 40″ E) and in the enclosed section of Ise Bay (34°55′ 12″ N, 136°44′ 25″ E) (MLIT). Measurements by MLIT in 2002 at the river mouth showed that waves occur with a period of 5 to 10 s and are 1-2 m in height, increasing to 4-5 m in height during typhoons and winter storms (http://www.mlit.go.jp/kowan/ nowphas/index.html). In the investigated area, waves mainly wash from the northeast throughout the year, but also occur from the northwest or north-northwest during winter seasons (Nakajo 2004). Tidal residual current is mainly from northwest to southeast in the southwest part The geomorphology of the Kushida River Delta is classified into four components, namely, spit, salt marsh, muddy tidal flat, and tidal sand flat (Fig. 2).

Spit
The spit extends northwest from Matsunase Harbor and branches into several components. The spit is composed of well-sorted medium-to coarse-grained sand. Branching channels from the Kushida River flow behind the spit (Fig. 2). These channels are 20-30 cm in depth and transport medium-to very coarse-grained sand, including granules and pebbles.

Salt marsh
The salt marsh lies mainly on the landward side of the spit and is composed of muddy sediment and soil containing abundant organic material (Fig. 2). Sediment coarsens toward the main fluvial channel (Kushida River).

Muddy tidal flat
Muddy tidal flats are located on the landward side of the spit, adjacent to the salt marsh (Fig. 2). The sediment consists of mud or poorly sorted sandy mud, including gravels. Small-scale washover fans are observed on the muddy flat. Burrows of crustaceans are very common, and mud-flat sediment and coarse sediment of washover fans are intermixed by their benthic activities.

Tidal sand flat
The microtidal sand flat extends in front of the spit, on the right bank of the river mouth (Fig. 2). Its maximum area, during ebb spring tide, is approximately 0.4 km 2 (Fig. 2). The offshoreward development of the intertidal sand flat is most prominent around the river mouth and the southeast part of the study area (adjacent to Matsunase Harbor) (Fig. 2). The primary geomorphological features of the sand flat consist of sand bars and braided channels, and no tidal channels are observed. The sand bars are 60 to 100 cm in height and 50 to 100 m in wavelength (Fig. 2). They are arranged almost parallel to the shoreline and normal to the direction of incoming waves (Fig. 2). Breaking waves are found on the crests of these sand bars, and wave ripples are abundant on their surface. No significant migration of these bars was observed in the period from 2002 to 2011 (Nakajo 2004; Yamashita et al. 2009). Braided channels (approximately 20-30 cm deep) are observed intermittently on the sand flat, around the river mouth. These channels are temporary braided streams from the river. They are generally produced by floods and rapidly migrate and vanish. There are very few muddy particles on the sand flat, except for soon after a fluvial flooding event (Nakajo 2004; Yamashita et al. 2009;Yamashita et al. 2011). Sediment transport on the sand flat is mainly controlled by fluvial and wave activities (Yamashita et al. 2009).

Sampling and grain-size measurements
Bottom-surface sediments on the microtidal sand flat were collected during the ebb spring tides on 1 August 2011. A total of 60 sampling sites were located at horizontal intervals of approximately 100 m, both parallel and perpendicular to the spit, with extension to the river mouth and the northern margin of the sand flat (Fig. 2). Sediments were excavated by shovel in intervals, from the top of the bedforms to 8 cm below the surface, corresponding to the active layer thickness observed on 31 July 2011. Each sample (50-100 g weight) was homogenized before grain-size analysis.
The grain-size distribution of sand-sized particles in the samples was measured by the settling tube method (Gibbs 1974;Tucker 1988) after removing mud-sized (> 4φ) and gravel-sized (< − 1φ) particles by sieving. Shell fragments were also removed before measurement. Grain-size distributions were calculated using S-Tube software (Friedman et al. 1992;Naruse 2005). The conversion method for the settling velocity of grains was proposed by Ferguson and Church (2004), using an empirical formula for natural sand. Model validation by monitoring of actual sedimenttransport patterns Actual sediment-transport patterns on the tidal sand flat were monitored from 16 May to 31 July 2011, to enable validation of the GSTA models. Therefore, the sedimenttransport patterns were examined in both normal (calm) conditions and in an episodic depositional event. The thickness of the active layer was also measured; the active layer was considered to be the mixed layer of tracer grains and natural sands. The JMA reported rainfall measured in the upper reaches of the Kushida River and the water level of the river mouth (http://www.data.jma.go.jp/obd/stats/etrn/ index.php; Fig. 3). There was moderate rainfall (< 15 mm/h) during the period from 16 May to 3 June. The water level indicates a usual spring-neap-spring cycle in this period (Fig. 3). In contrast, heavy rainfall was recorded on 8 July (> 45 mm/h) and 17-19 July (maximum: > 25 mm/h). In particular, unusually high water levels during ebb tides were monitored soon after the 17-19 July rainfall (18-20 July) (Fig. 3). The rainfall during this period was due to a typhoon (No. 6, MA-ON). Fluvial flooding events caused by the rainfalls during this period were reported through an enquiry survey.
Tracer experiments and geomorphological observations revealed net sediment-transport patterns on the tidal sand flat. Blue-colored sand was employed as tracer sediment, as it can be clearly identified from natural sand, both visually and by microscopic observations. The tracer sediments were well sorted and their mean grain size was ca. 1.5 φ. The tracer sediments were situated at five localities (points 1-5) on 16 May 2011 (spring tide) (Fig. 4). GPS and marking poles recorded the localities. At each locality, the bottom surface sediment was excavated approximately 0.3 (length) × 0.3 (breadth) × 0.2 (depth) m, and the hole was filled by the tracer sediments (approx. 25 kg). Two weeks after input of the tracer sediments (3 June 2011), bottom surface sediments were collected from four surrounding points that were located at a distance of 10 m from each tracer pit. Each sample was homogenized with sufficient stirring by hand in a transparent plastic case (8 × 15 cm basal area, and 3 cm in height). The basal sediment in this plastic case was scanned at high resolution (2400 dpi), and the number of colored sand grains in the basal area was counted to estimate the volume of surface sediments transported in each direction from the original tracer site. Sampling was also conducted 2 months after the tracer sediments were set (31 July), but no tracer grains were detected.
Mean transport direction of sediments θ was estimated from observed colored sand grains, as follows: where θ i is a direction from the tracer pit to the sampling locality at which the ith grain was detected. The mean resultant length (MRL) L MR of the vector was employed as a proxy for concentration of sedimenttransport directions, which was calculated as follows: where n is the number of detected colored grains. Furthermore, the thickness of the active layer of basal Page 7 of 18 sediment was measured after excavating sediments at the source points (Fig. 5).
The sediment-transport patterns were also estimated on the basis of the observation of geomorphological development (Figs. 4 and 5) during the period from 16 May to 31 July 2011.

Results of grain-size analysis
Spatial distributions of eight granulometric parameters (mean and median grain size, sorting, CV, skewness, kurtosis, and mud and gravel fractions) were calculated and were interpolated to a grid interval of 100 m by kriging (Figs. 6 and 7). Semivariograms were calculated from the measured data, and the variogram models were fitted to them by the least-squares method (Fig. 6). The estimation error is shown as a false-color image, and the value of the parameter is shown by an overlying contour map (Fig. 7).

Mean and median grain size
The interpolated maps of the mean and median grain size show similar characteristics. The coarser sediments (< 1.0 φ) are distributed around the river mouth (Fig. 7). The grain size becomes finer in a radial pattern from the river mouth, and grain sizes finer than 1.50φ are found at the southeast part of the sand flat (Fig. 7). The variograms of the mean and median grain size reveal their high spatial dependency (Fig. 6). Therefore, they have relatively high accuracy (Fig. 7).

Sorting
The highest sorting values (> 0.56) are found in the central to southeast part of the sand flat (Fig. 7). The sorting value tends to decrease seaward, but a lower value (< 0.55) is also evident in the river mouth (Fig. 7). The semivariogram of the sorting value suggests relatively high spatial dependence (Fig. 6). However, the y intercept of the fitted curve for semivariance of the sorting value largely deviates from the origin (nugget effect), indicating that high measurement error exists (Fig. 6).

CV value
The map of the CV value indicates that higher values (> 0.7) are distributed in the river mouth. The Fig. 6 Semivariograms and fitted variogram models of granulometric parameters. Variogram models were fitted by the least-squares method, and the model that produces the fewest residual errors was chosen from among the Gaussian, spherical, and exponential models. The spherical model was used for mean and median grain-size, sorting, and gravel fractions, the exponential model was used for the CV value, and the Gaussian model was used for skewness, kurtosis, and mud fractions CV value decreases seaward (Fig. 7). The variogram of the CV value shows high spatial dependency (Fig. 6); therefore, the accuracy of the CV value is relatively high.

Skewness
The higher skewness values (> 1.4) are mainly distributed in the river mouth and are also found in the northern part of the sand flat (Fig. 7). Skewness tends to Page 10 of 18 decrease to the southeast (Fig. 7). The lowest value (< 0.2) is distributed primarily in the southeastern margin of the sand flat (Fig. 7). The accuracy of interpolation is low because of the anomalous value of the semivariogram in the range of 100 m (Fig. 6).

Kurtosis
The higher kurtosis values are mainly distributed in the river mouth (7.0-9.0), and the value tends to decrease in a seaward direction (Fig. 7). However, high values are also found in the southeastern margin of the sand flat. The accuracy of interpolation is relatively low because the variogram indicates anomalous values in the range of 100 m (Fig. 6).

Mud fractions
The spatial distribution of the mud fractions is quite random, and therefore, the accuracy of interpolation is very low (Fig. 6). Relatively high mud fractions (> 10.0%) are distributed in patches in the northeastern and southeastern parts of the sand flat (Fig. 7).

Gravel fractions
Gravels are abundant in the sediment around the river mouth (> 9.0%) and decrease seaward, becoming sparse in the southeast part of the sand flat (< 6.0%) (Fig. 7). The variogram shows high spatial dependency of the gravel fractions (Fig. 6); therefore, the accuracy of this value is relatively high (Fig. 7).

Results of principal component analysis of granulometric parameters
The contribution values, which indicate the percentages of the original information represented by each principal component of PC1 through PC3, are 56.77, 20.32, and 11.76, respectively (Table 1). This means that more than 80% of the original information is represented by PC1 through PC3. Therefore, these three principal components are interpreted further. Kriging interpolation was conducted in order to visualize the spatial distribution patterns of scores of PC1 through PC3, using the same method as for the original granulometric parameters (Figs. 8 and 9).

PC1: sediment-transport function
PC1 is negatively correlated with the gravel log-ratio (factor loading: − 0.79) and CV value (− 0.92) and is positively correlated with the median grain size (φ) (0.93) ( Table 1). Correlations with skewness (− 0.57), kurtosis (− 0.58), and mud log-ratio (0.65) have relatively weak influences on this value (Table 1). In other words, the increase of the PC1 value indicates that the sediment becomes finer, better sorted, and less gravely (Fig. 10); this means a selective decrease in coarser materials from grain-size distribution (Fig. 10). The variogram and interpolated map of the PC1 scores indicate that the spatial variation of PC1 is high (Figs. 8 and 9). The lowest value of PC1 is found around the river mouth, and this value increases radially northeastward (Fig. 9). The scatter diagram of PC1 values against PC2 values indicates that the samples take wider values of PC1 (from − 4.5 to 4.0) than values of PC2 (from − 3.0 to 2.0); no notable clusters were formed (Fig. 11). Thus, it is interpreted that PC1 reflects the suite of size-selective sediment-transport processes resulting in fining-downcurrent trends with well-sorted grain-size distributions. The gradual decrease of coarser materials, which is represented by the increase of PC1, probably reflects the progressive lowering of bed shear stress of both fluvial and wave activities (Fig. 10). Settling velocity of finer materials is low, so that they have greater mobility as suspended loads. Therefore, the finest fraction of sediments is broadly dispersed downcurrent, producing the negative skewness of grain-size distribution. The relatively high spatial variation and the radial distribution pattern from the river mouth in PC1 strongly suggest that the PC1 value represents the sediment supply from the Kushida River and the dispersal by waves or tidal currents (Fig. 6). Therefore, the decreasing pattern of PC1 can be regarded as the sediment-transport function.

PC2: indicator of lag deposit of a fluvial flooding event
PC2 is negatively correlated to the skewness (− 0.65) and kurtosis value (− 0.66), and positively influenced by gravel log-ratio (0.55) ( Table 1). Other granulometric parameters have only slight influences on PC2 (Table 1). In other words, the increase of the PC2 value indicates the more negatively skewed and platykurtic shape of grain-size distribution and gravely sediment (Fig. 10). The higher values of PC2 (1.2-2.0) cluster around the breached section of the spit and the newly formed braided channel on the sand flat (Fig. 9). This value decreases around the river mouth and offshore area (< 4.0). The spatial dependency and precision of interpolation is low because of the anomalous values of semivariance at approximately 100 m (Fig. 8). The exclusive distribution of higher PC2 values  around the newly formed braided channel suggests a linkage to the episodic deposition by the fluvial flooding event during July 2011. The platykurtic and negatively skewed shape of grain-size distribution as well as the high gravel fraction indicates that the sediments largely consist of coarser fractions and this is in agreement with the direct observation of a high gravel fraction (Fig. 10). These coarser fractions can be interpreted as lag deposits in braided channels. Coarser sediments may have been episodically deposited during periods of high shear stress flow at times of fluvial flooding and remained through the subsequent reworking by wave or tidal currents.

PC3: indicator of muddy particles
PC3 is positively influenced by the mud log-ratio (0.74), and other granulometric parameters show limited influences (from − 0.02 to 0.33) ( Table 1). This means that the increase in the PC3 score suggests muddier sediment. The variogram and interpolated map of PC3 show that the spatial dependence of PC3 is weak (Figs. 8 and 9). The spatial distribution of PC3 is very heterogeneous. Higher PC3 values (> − 0.3) are found in the northern and southeastern parts of the sand flat, and this characteristic distribution pattern is similar to the map of the mud content ( Figs. 7 and 9). The range of variation in PC3 values (from − 2.0 to 2.5) is more limited than for PC1 (from − 4.5 to 4.0) (Fig. 11). Very poor correlation between mud logratio and other granulometric parameters of sand-sized particles indicates that muddy particles are deposited by quite different processes from those associated with sandy particles. As described above, the rapid deposition of muddy particles on the sand flat by a fluvial flood event (17-19 July 2011) is recognized, and these muddy particles preferentially remained in local depressions. The heterogeneity of the spatial distribution of PC3 might reflect a relic of mud deposition caused by this event and their local resistance to erosion.

Trend vectors obtained by P-GSTA and previous GSTA methods
Grain-size trend vectors were calculated by the P-GSTA method, in which the increase of the PC1 value is  defined as a sediment-transport function (Fig. 12). Trend vectors directed to the north or northeast are dominant along the main fluvial channel and northern parts of the sand flat (Fig. 12). Northeastward-directed vectors are dominant around the breached section of the spit. In the central part of the sand flat, the vectors are directed to the southeast. In the southeastern regions, the trend vectors are directed toward the northeast (Fig. 12). To summarize, the trend vectors obtained by the P-GSTA method show radial dispersion of coastal sediment from the river mouth with a small influence from the shoreparallel current.
The previous GSTA model ( Gao and Collins 1992;Asselman 1999) was also applied for the purposes of comparison. Trend vectors obtained by the previous GSTA method showed quite diverse and inconsistent estimates of sediment transport pathways (Fig. 12). Landward-directed vectors were dominant along the main fluvial channel, while offshore (southeastward) directed vectors were found at the central part of the sand flat. In addition, transport vectors were not defined around the southeastern region of the study area because neither FB− nor CB+ patterns were satisfied in this region.

Results of tracer experiments and topographic observations
The tracer experiments clarified the actual sedimenttransport patterns during the period 16 May to 3 June (Fig. 4). Along the main fluvial channel, the transport pattern was directed to the northeast (point 5) (Fig. 7) and the degree of concentration (MRL) was relatively high (Fig. 7). In the central part of the sand flat (points 3 and 4), southward (landward) transport was dominant (Fig. 4) and the degree of concentration (MRL) was moderate (Fig. 4). In the offshore area (point 2), eastward sediment transport was dominant (Fig. 4) and the degree of concentration was moderate (Fig. 4). At the southeastern margin of the sand flat (point 1), sediment dispersed in various directions (Fig. 4); northeastward transport is dominant, according to the number of tracer grains, but the mean direction indicates east-southeastward transport (Fig. 4).
On 31 July 2011, no tracer grains were detected from the samples, probably because they had been highly dispersed. Instead, sediment-transport patterns were estimated on the basis of the observation of the geomorphological development derived from a fluvial flooding event of 17-19 July (Fig. 3). The most prominent feature was incision of the spit and formation of braided channels (Figs. 4 and 5). The spit was breached at two points around its folding point, 30-50 m in width (Figs. 4 and 5). At the breach locations, the spit was incised by shallow, braided channels (0.3-0.5 m in depth), which were connected to branched channels behind the spit (Figs. 4 and 5). These channels flowed into the tidal sand flat more than 300 m offshore (Figs. 4 and 5). Sand bodies were deposited on the sand flat 100 m off the breach point (Figs. 4 and 5). These sand bodies extended approximately 100 m northeast and were approximately 0.5 m in height (Fig. 4). They were composed of medium-to coarse-grained sands, including shell fragments, which were similar to sediments comprising the spit and the fluvial channels, and were clearly coarser than the sediments on the tidal sand flat. This probably means that these sand bodies were formed by sediments derived from the breached part of the spit and from Fig. 12 Grain-size trend vectors. The actual sediment-transport directions during normal conditions are shown as the mean direction of tracer results (see Fig. 4a). Circles indicate the sampling points. a Grain-size trend vectors calculated by P-GSTA. b Grain-size trend vectors calculated by the previous GSTA method (Asselman 1999). Some of sampling points lack the trend vector because the GSTA method does not recognize the grain-size trend when no surrounded locations satisfy the conditions FB− or CB+ sediments in the fluvial channels. Dunes of approximately 1.0 m wavelength and 0.15 m amplitude formed in the braided channels (Fig. 5). These dunes exhibited the northeastward (offshoreward) current directions around points 3 and 5 (Figs. 4 and 5). At the eastern part of the sand flat (around points 1 and 2), muddy deposits were thickly deposited in the troughs between sand bars. The average thickness of the muddy sediment was 1-3 cm, but it occasionally reached 5 cm. No other notable geomorphological developments were observed in this area.
The active layer was clearly shown as the mixed layer of tracer grains and natural sands, and the boundary with the underlying substratum was sharp (Fig. 5). The thickness of the active layer was ca. 8 cm at point 3 (Fig. 5). The thickness was not measured at other locations due to the loss of the marking poles.
Actual sediment-transport patterns on the sand flat during a summer season As discussed above, sediment-transport patterns were measured during a summer season (May to July 2011), via tracer experiments and geomorphological observation. The results showed that fluvial discharge, mainly from the braided channels that accompanied the incision of the spit, supplied sediments to the sand flat; these drifted southward at the central part of the sand flat and were dispersed eastward or northeastward by wave activity during normal conditions.
The drastic geomorphological development associated with the 17-19 July fluvial flooding event revealed episodic sediment-transport processes around the river mouth. Fluvial sediments were supplied to the sand flat toward the northeast, mainly from braided channels incising the spit, rather than from the main fluvial channel (Fig. 4). The incision of the spit and fluvial sediment supply was also observed, in October 2009, by Yamashita et al. (2011), who reported the northeastward transport of fluvial sediment into the sand flat by a fluvial flooding event.
It is estimated that the fluvial sediments supplied to the sand flat were dispersed by wave or tidal currents during normal conditions. The eastward transport in a seaward area (points 1 and 2) roughly coincides with the direction of coastal (residual) currents (Coastal Oceanography Research Committee, The Oceanographic Society of Japan 1985; Figs. 1 and 4). The landward (southward) transport at points 3 and 4 can be interpreted as representing a depositional wave condition during summer season (Hayes 1972; Fig. 4).

Comparison with trend vectors obtained by P-GSTA and previous GSTA methods
Grain-size trend vectors calculated by the P-GSTA method showed reasonable similarity to the actual sediment-transport patterns described above (Fig. 12). A north to northeastward trend along the main fluvial channel and northern parts of the sand flat coincided with the net sediment-transport directions observed by the tracer experiment at point 5. Northeastwarddirected vectors were dominant around the breached section of the spit, and coincided with the downstream direction of the braided channels formed by the fluvial flooding event of July 2011 (Fig. 12). Thus, the trend vector of P-GSTA at point 5 coincides with the sediment transport direction estimated from geomorphological observation, whereas it contradicts that derived from the tracer data. In the central part of the sand flat, the vectors were directed to the southeast and coincided approximately with the tracer result at point 4 (Fig. 12). In addition, the vectors coincided very well with the tracer result at point 2, in which they were directed eastward. In the southeastern regions, the trend vectors were directed toward the northeast (Fig. 12), which is in slight disagreement with the tracer data, indicating eastward transport (Fig. 12).
On the other hand, sediment trend vectors calculated by the previous GSTA method did not show reasonable net sediment-transport patterns compared to the results for P-GSTA (Fig. 12). Around point 2, trend vectors showed eastward transport patterns and coincided with the tracer results (Fig. 12). Southeastward vectors found at the central part of the sand flat (around point 4) were also comparable with the tracer results. However, landward-directed vectors along the main fluvial channel clearly conflict with the tracer data at point 5 (Fig. 12). In addition, transport vectors were not defined along the braided channels and in the vicinity of point 3.

Discussion
In contrast with the results obtained using the previous GSTA method, the grain-size trend analyzed using the P-GSTA method closely matches the actual sedimenttransport estimated from the tracer experiments and geomorphological features observed soon after fluvial flooding events (Gao and Collins 1992; Asselman 1999; Fig. 12). In the GSTA method, trend vectors opposite to the fluvial downcurrent direction are produced by the lower sorting values in the river mouth (Figs. 7 and 12). The coarse sediment (lower φ value) in the river mouth results in a lower sorting value due to the inevitable correlation between mean grain size and sorting value, and therefore, sorting values cannot represent the actual degree of sediment sorting, as supposed by the previous methods. The previous GSTA method performs quite poorly, in that no trend vector is defined around the braided channels, despite direct observations of the sediment transport by a fluvial flooding event in July 2011 (Figs. 3,4,5,and 12). In contrast, sediment-transport patterns are better reconstructed via the P-GSTA method, especially around the river mouth and braided channels (Fig. 12). At point 3, trend vectors are parallel to the downcurrent direction of braided channels and thus represent sediment-transport direction by the fluvial flooding event, rather than during normal conditions indicated by tracer experiments (Fig. 12). This is probably because the sampling for P-GSTA was conducted after the flooding event while the tracer experiments were conducted before the event (Fig. 3).
The disparity between trend vectors calculated by the GSTA method and observational data shows that the assumption of McLaren and Bowles (1985) is not appropriate in this case, probably because the dominant sediment-transport processes are mixed. In contrast, the P-GSTA method can detect the unique grain-size trend in each depositional setting and is therefore effective even under mixed transport processes, in which the precise sediment-transport function is difficult to determine. Principal component analysis (PCA) can easily explore the characteristic changes of the grain-size distribution curves, which result from size-selective sediment transport in a specific depositional setting, by a linear combination of multiple grain-size parameters (Eq. 10). PCA can automatically weight multiple grainsize parameters and combine them into a linear function. In addition, the magnitude of the weighting is proportional to the spatial variation within each parameter (Middleton 2000). Therefore, the resultant linear combination of grain-size parameters, summarized by the first principal component (PC1), represents the predominant deformation pattern of the grain-size distribution curves in the sampling area. In other words, the first principal component can generally be regarded as the primary function of size-selective sediment-transport in which a specific grain-size distribution curve varies spatially. If some minor transport processes exist, it is expected that the deformation pattern of the grain-size curves affected by these processes can be identified by the subsequent principal components. In the present study, PC2 and PC3 recognized two minor depositional processes (channel-lag deposits and muddy deposits resulting from episodic flooding). Indeed, the trend vectors derived from GSTA show directions incompatible with observations in the region around the spit (Fig. 12). In this region, PC2 values are large (1.2-2.0), suggesting the influence of channel-lag deposits formed by washover processes (Fig. 9). It is notable that the trend vectors of P-GSTA show directions consistent with the observation of sediment transport (Figs. 5b and 12). On the other hand, PC3 shows smaller values (< − 1.5) around the river mouth and rhythmic bars (Fig. 9), implying the presence of fine-grained flooding deposits. Again, the trend vectors derived from GSTA are incompatible with the directions of actual sediment transport (Figs. 5 and 12), which suggests the effectiveness of the P-GSTA method.
Although P-GSTA is a very flexible method, a number of precautions must be taken in applying this method. First, careful selection of the effective parameters that represent the shape of the grain-size curve is necessary. If too many irrelevant parameters are included, the results of the PCA will be meaningless. Conversely, parameters that have a mathematically inevitable correlation will be erroneously emphasized, according to their spatial variances. Parameters should preferably be dimensionless values that are as independent from each other as possible. In the present study, median grain size and CV value are employed instead of mean grain size and sorting value, which are inevitably correlated to other granulometric parameters. In addition, the consistency of the sediment-transport function should be considered in some cases. If more than two sediment-transport functions dominate in a depositional setting, PCA alone cannot detect them. To deal with these situations, before PCA, samples should be divided into several groups on the basis of the similarities of grain-size characteristics.
Furthermore, the applicability of all GSTA methods should be verified in various depositional environments. GSTA methods preferentially target depositional environments rather than erosional environments. The bottom-surface sediments sampled for GSTA methods should be grains being actively exchanged with grains in a bedload layer, not the 'historical layer' exposed by erosion. In addition, attention should be paid to topographical undulation. Masselink et al. (2008) introduced a case of failure of the previous GSTA method (Gao and Collins 1992) on a tidal shoreface and described this failure in connection with the geomorphological variations in shoreface bars and runnels. In the case of regions such as shoreface bars or runnels where topographic undulation is highly developed, the active-layer thickness varies considerably within a depositional environment. Thus, GSTA methods are preferentially applied in depositional environments where topographical undulation is relatively minor and deposition occurs throughout the surveyed region; otherwise, the active-layer thickness should be identified at all sampling locations.

Conclusions
1. Principal component analysis (PCA) revealed the grain-size trend in a modern microtidal sand flat along the Kushida River Delta, central Japan. PCA enables multiple grain-size parameters to be organized into a simple linear function with different weights on the basis of their spatial variation. 2. The grain-size trend analysis using the results of PCA (P-GSTA) revealed that sediment becomes finer, better sorted, and less gravely through major transporting processes and reconstructed sedimenttransport pathways in the microtidal sand flat. 3. Application of the P-GSTA method revealed the presence of local depositional processes on the sand flat: deposition of fluvial-channel lags and muddy particles.