Inverse analysis to reconstruct hydraulic conditions of non-steady turbidity currents: Application to an ancient turbidite of the Kiyosumi Formation of the Awa Group, Boso Peninsula, central Japan

This study proposes a new method of inverse analysis from ancient turbidites to the non-steady turbidity currents with consideration of multiple grain-size classes. The forward model employed in this study is based on the shallow water equation, and the initial condition of ﬂows are assumed to be the lock-exchange type condition. To obtain a solution of the inverse problem, this study employed the genetic algorithm for ﬁnding the optimized initial conditions. The present method successfully estimated the true given initial conditions of the turbidity currents from the artiﬁcial data sets of deposits created by the calculation of the forward model. The author also applied the method to the turbidite bed in the Kiyosumi Formation. As a result of inverse analysis, the obtained solution ﬁts well to the observed data of the individual turbidite, providing estimates of the ﬂow velocity, the ﬂow thickness and the sediment concentration of the turbidity current. The ﬂow thickness and velocity when the turbidity current reached at the downstream end of the study area were reconstructed to be 334.6 m, 0.98 m/s respectively at the location of the downstream end. Result of our analysis is the ﬁrst example of reconstructing a reasonable conditions of the turbidity current from an ancient turbidite observed in the ﬁeld, and the method id expected to be applied in various regions in the future.

A remaining problem in this field of research is the quantitative understanding of the developmental process of turbidity currents at natural scales, which is essential to estimate the distribution of ancient turbidites from the limited amounts of information provided by cores and seismic profiles (Takano, 2016). Although several in-situ measurements of turbidity currents on deep sea floors have been reported (Shepard, 1965;Inman et al., 1976;Dengler et al., 1984;Xu et al., 2004;Vangriesheim et al., 2009;Arai et al., 2013;Andrieux et al., 2013), hydraulic conditions of turbidity currents such as flow velocity and sediment concentration, however, remain unclear (Kubo, 1995;Falcini et al., 2009;Lesshafft et al., 2011); conducting such in-situ measurements is quite difficult because of their highly destructive nature and infrequent occurrences (Naruse and Olariu, 2008;Talling et al., 2015). The in-situ measurements cost and it is easer to observe ancient turbidites than turbidity currents in fieldwork.
In recent years, there have been several attempts to reconstruct the hydraulic conditions of ancient turbidity currents from geologic records (Stow and Bowen, 1980;van Tassell, 1981;Bowen et al., 1984;Komar, 1985;Allen, 1991;Hiscott, 1994;Kubo, 1995;Kubo et al., 1998;Baas et al., 2000;Naruse and Olariu, 2008;Falcini et al., 2009;Lesshafft et al., 2011). Previous reconstructions of turbidity currents from their deposits were mostly based on the concepts of the capacity or competence of a flow (Dorrell et al., 2013). The competence of a flow is usually defined as the largest grain size that the flow can transport, while the capacity of a flow is defined as the amounts of sediment particles that the flow can potentially transport at a given flow conditions (Kuenen and Migliorini, 1950). Initially, the critical velocities of particle motion (i.e. the competence of a flow) inferred from turbidites were used for estimating the paleo-flow conditions of turbidity currents (Komar, 1985;Kubo, 1995;Kubo et al., 1998). However, these methods based on critical velocity estimate only the minimum values of flow velocity, because deposition can occur even at velocities much higher than the critical velocity of particle motion when the capacity of sediment transport in a turbidity current is exceeded by the volume flux of sediment (Hiscott, 1994). Instead, Hiscott (1994) emphasized the role of the capacity of a flow to estimate paleoconditions of turbidity currents, but the methodology for reconstruction has not been implemented yet. Further discussions about relations between the competence and capacity of a flow were given in (Dorrell et al., 2013). Analyses of the sedimentary structures of turbidites can be conducted as alternative methods to reconstruct flow conditions of turbidity currents (Baas et al., 2000). For example, (Dorrell et al., 2011) suggested that poorly graded division in turbidites is a product of high-concentration suspension. However, analyses of the sedimentary structures only provide a range of flow conditions (Ohata et al., 2017), not the spatio-temporal development of flows.
Here we propose a new method for inverse analysis to reconstruct the paleohydraulic conditions of turbidity currents based on ancient turbidites. We focus on the possibility of elucidating the history of the entire turbidity flow current based on the properties of the ancient turbidite, such as sedimentary structures and grainsized distribution. The methods of inverse analysis have been proposed in several recent studies (Naruse and Olariu, 2008;Falcini et al., 2009;Lesshafft et al., 2011); hydraulic parameters are estimated based on optimization of the input parameters of numerical models to fit the results of calculations with field observations of turbidites. Flows calculated with optimized parameters can be regarded as reconstructions of the actual turbidity currents that emplaced the ancient turbidites. Inverse analyses have also been applied to other gravity flows, such as pyroclastic flows (Rossano et al., 1996), tsunamis (Jaffe and Gelfenbaum, 2006;Soulsby et al., 2007;Mitra et al., 2020), and tidal channels (Masuda and Nakayama, 1988), to reconstruct hydraulic conditions from the corresponding deposits. Inverse analysis of turbidity currents is still in an early stage of development, and there is much remains room for improvment. For instance, Falcini et al. (2009) estimated the hydraulic conditions of turbidity currents from ancient turbidites of the Laga Formation in the Central Apennines, Italy. Their forward model employed the assumption of steady flow; that is that the depositional rate did not vary over time. However, ancient turbidites are characterized by graded bedding, which suggests strongly non-steady conditions. Consequently, the applicability of the model of Falcini et al. (2009) may be limited to the specific cases of actual ancient turbidites that show no grading. In contrast, Lesshafft et al. (2011) employed a direct numerical simulation (DNS) model of the vertical two-dimensional Navier-Stokes equations for inverse analysis from turbidites. However, the calculation cost of their method is so high that application of the method is not feasible with field-scale data obtained from outcrops or boreholes.
To resolve the problems associated with conducting inverse analysis of naturalscale turbidity currents, we employed a non-steady model with consideration of multiple grain-size classes as a forward model, which can describe the spatio-temporal behavior of a turbidity current that deposits a typical turbidite with graded bedding. Our model of turbidity currents is based on that of Parker et al. (1986) with modifications to account for multiple grain-size sediments. The active layer concept from Hirano's sediment continuity model Hirano (1971) was used for calculating entrainment and deposition of multiple grain-size. A combination of the shallow water equations and an active layer were used for flume experiments (e.g., Suzuki et al., 1976) and numerical model of river bed degradation (Blom, 2008;Stecca et al., 2016). Because computational cost of the forward model is significantly lower than that of two-dimensional models, this method can be applied to field-scale problems. In the forward model, the "lock-exchange" model was employed as the initial setting for numerical simulation with various conditions. Furthermore, the forward model was tested for sensitivity through examples of forward model calculations. For inverse analysis, an objective function is defined by sum of squares of deviations between the observation result and the numerical calculation results. In this inverse calculation, the initial hydraulic conditions that minimize the objective function were explored with the genetic algorithm. The optimum solution was then treated as the paleo-hydraulic conditions. For an application of this inverse analysis at fieldscale, a turbidite deposit in the Kiyosumi Formation, Boso Peninsula, Japan, was investigated.

Forward model of non-steady turbidity currents
Governing equations of non-steady turbidity currents A depth-averaged model of a non-steady one-dimensional turbidity current is described herein. This forward model is based on that of Kostic and Parker (2006), the results of which were verified by the experimental results of Garcia (1990). This model is classified as a "three-equation" model after Parker et al. (1986), based on the equations representing fluid mass conservation, streamwise momentum conversation and mass conversation of suspended sediment (Fig. 1).
Our model accounts for transport and deposition of sediment showing non-uniform grain-size distribution that is discretized to multiple grain-size classes (equations 1 to 3). The Exner equation describes mass conservation of sediment in bed for each grain-size class (equation 4). In addition, the continuity equation of sediment of each grain-size class in the active layer (Hirano, 1971) is represented in order to calculate the entrainment rate of bed sediment of each grain-size class (equation 5). As presented by Kostic and Parker (2006), these relations are expressed as follows: where x is the bed-attached streamwise coordinate, and t is time. In the above relations, H, U and C i denote the thickness of a turbidity current, the layer-averaged velocity and the layer-averaged volumetric concentration of suspended sediment of the ith grain-size class, respectively. In this study, the grain-size distribution of sediment is discretized to multiple grain-size classes at regular intervals. The parameter C T indicates the layer-averaged volumetric concentration of total suspended sediments (thus C T = C i ). The acceleration of gravity and the slope gradient are denoted g and S, respectively. The parameter u * denotes the shear velocity. Properties of sediment particles are described by the parameters R, w i and λ p , which represent the submerged specific density of sediment, the fall velocity of a sediment particle of the ith grain-size class, and the porosity of bed sediment, respectively. The parameter η i denotes the volume per unit floor area of bed sediment of the ith grain-size class, and η T indicates the total bed elevation (thus η T = η i ). The parameters related to the active layer formulation are the active layer thickness L a and F i , which indicates the volume fraction of the ith grain-size distribution in the active layer. The parameters e si , e w and r 0 are the entrainment rate of sediment of the ith grain-size class into suspension, the entrainment rate of the ambient water to the flow, and the ratio of near-bed suspended sediment concentration to the layer-averaged value. These parameters require empirical relations to close the equations, which are described in the next section. As defined in equation 1, the volume of the turbidity current increases downstream in this model because of the entrainment of ambient water to the flow (Fig. 1). Equation 2 is the momentum conservation equation, which indicates that the turbidity current is driven by the fluid pressure and density difference between the turbid fluid containing suspended sediment and the ambient fluid (Fig. 1B). In this model, it is assumed that the suspension is dilute enough to justify employing the Boussinesq approximation in Equation 2. In addition, hindered settling is ignored in this model. Equation 3 is the relation of the sediment conservation, in which both settling from suspended sediment and entrainment from bed sediment are assumed to occur concomitantly ( Fig. 1). Equation 4 is the Exner equation representing the mass conservation of bed sediment (Fig. 1). Equation 5 is the relation of the sediment mass conservation of the ith grain-size class in the active layer (Fig. 1). This equation describes the temporal change of grain-size distribution in the active layer (Hirano, 1971). The governing equations described above can be recast in dimensionless form as follows (Kostic and Parker, 2003a;Kostic and Parker, 2003b): wherex andt are a dimensionless streamwise coordinate and dimensionless time, respectively. The parameterû * denotes the dimensionless shear velocity, andŵ is the dimensionless fall velocity of the sediment. These dimensionless parameters are given by the following relations, respetively: The flow thickness H, the depth-averaged velocity U , and the depth-averaged volumetric concentration of suspended sediment C i are targets to for numerical estimation in this forward model. Corresponding dimensionless forms of these parameters areĤ,Û , andĈ i , which are defined aŝ where H 0 and C i0 denote the initial values of H and C i at the upstream boundary. U 0 is the initial velocity of the flow head described below in this section. The parameter Ri 0 in Equation 7 is the inflow bulk Richardson number defined as: To solve the equations described above, a transformed coordinate system was employed in this study based on the formulation of Kostic and Parker (2006) (Fig. 2). This model of turbidity currents involves two boundary conditions for the governing equations: the upstream influx boundary (the tail of the current), and the downstream propagating boundary (the head of the current). To address these boundary conditions, a deforming grid approach was adopted, in which the moving downstream boundary is defined as the head of the turbidity current at the fixed point x * = 1. Thus, where the parameterŝ is the dimensional position of the turbidity current head, when x * = 1. After applying the transformation of Equation 17, equations 1 to 3 of the "threeequation" model of Parker et al. (1986) may be re-expressed in the following conservative form: whereṡ is the velocity of the current head andŝ is the dimensionlessṡ.

Closure of equations
Several empirical relations are required to close the equations described above. Shear velocity u * is related to the layer-averaged flow velocity according to the following relation: Here the friction coefficient c f is assumed to be constant (0.0069). The particle settling velocity for each grain-size class with a representative diameter D is calculated from the relation of (Dietrich, 1982) for natural sand, which can be expressed as where the fit parameters b 1 , b 2 , b 3 , b 4 and b 5 are 2.891394, 0.95296, 0.056835, 0.000245, and 0.000245, respectively. The dimensionless parameter e w describes the rate of entrainment of ambient water into the turbidity current from above. The empirical formulation of Fukushima et al. (1985) for e w is given by: The parameter e s is a dimensionless coefficient for characterizing the rate of sediment entrainment into suspension by a turbidity current. This parameter is obtained from an empirical relation as follows (Wright and Parker, 2004): where the constants α 1 and α 2 are given as The parameter S f denotes a friction slope, which takes the form: The parameterṡ is the velocity of the turbidity current head and is used for calculation of the head position on the transformed coordinate system at calculation time steps. In contrast, the velocity of the turbidity current head in downstream conditions in governing equations is used for the densimetric Froude number Fr d , which is given by: The parameter r 0 is the ratio of the near-bed sediment concentration to the layeraveraged suspended sediment concentration C i for the ith grain-size class. Here we assumed that r 0 is constant (1.5) based on the experimental results of Garcia (1990) . Other input parameters are: the porosity of the bed sediment λ p = 0.4 and the thickness of the active layer L a = 0.003 m.

Implementation of the forward model
The set of equations 4, 5, 18, 19, and 20, which are based on the "three-equation" model of turbidity currents extended for treatment of multiple grain-size classes, can be calculated numerically to solve the morpho-dynamic problem. Ultimately, the five unknown parameters (Ĥ,Û ,Ĉ i , η i and F i ) in the set are obtained spatiotemporally using the numerical methods. The MacCormack scheme was used for integration of the partial differential equations 18 to 20. The predictor-corrector technique was used for the ordinary differential equations 4 and 5. The partial differential equations 18 to 20 employ the transformed coordinate system, whereas the ordinary differential equations 4 and 5 use the dimensional coordinate system to avoid the apparent advective transport. The linear interpolation technique was used to convert the parameters between these two coordinate systems for each time step.
The calculation of this numerical model is terminated when any one of the following three conditions is satisfied: (a) the sediment concentration becomes extremely thin (Ĉ T /C T 0 ≤ 0.001), (b) the flow reaches the calculation domain end on the downstream side, or (c) the time of calculation exceeds the prescribed value.

Initial and boundary conditions
In this study, the "lock-exchange" model was employed for the initial conditions of turbidity currents. A turbidity current is generated by sudden release of a suspended sediment cloud. In this model, suspended sediment is initially stirred in a confined region, and a turbidity current occurs when the lock gate at the downstream end of the confined region is opened. This setting has been adopted extensively in flume experiments of non-steady turbidity currents (Bonnecaze et al., 1993). Even in numerical simulation, this model has often adopted to simplify the initial flow conditions (Necker et al., 2005;Blanchette et al., 2005;Lesshafft et al., 2011). In this study, the initial suspended sediment cloud is set to H 0 and C i0 for flow thickness and sediment concentration of the ith grain-size class, respectively (Fig. 2). Both parameters are assumed to be uniform in a range of the initial along-slope length l 0 . At the beginning, the gate separates clear water and a sediment cloud. The sediment cloud is kept at a uniform concentration, and the gate is released at the timing of t = 0. The released turbidity current flows down along the slope, driven by the density difference between the flow and the ambient water.
The initial conditions of the slope at each grid point in the transformed coordinate system are set based on linear interpolation. The initial values of the dimensionless variablesĤ andĈ i are set to unity at each grid point. Also, the initial values of the dimensionless flow velocityÛ are set to 0 at the upstream end and 1 at the downstream end. The grain-size distributions in the bed and the active layer are assumed to be uniform at the beginning of the calculation, such that the parameters η i and F i are obtained as the initial bed thickness and 1 divided respectively by the number of grain-size classes N .
Both upstream and downstream boundary conditions are required for the numerical solution of the forward model. In this model, the flow velocity is set zero at the upstream boundary, and thusÛ | x * =0 = 0. The downstream boundary propagates with the movement of the flow head, such that sediment inflow is not allowed from the outside of the calculation domain at both the upstream and downstream boundaries. For the upstream boundary conditions of the flow heightĤ and the sediment concentrationĈ i , we employed the Neumann boundary condition, which takes the form: These equations (31) and (32) may be expressed in the following discretized form: With regard to the downstream boundary conditions, the model assumption in which the Froude number of the flow head Fr d is constant (1.2) yields the following relations: where the following relations are assumed for the mass conservation of fluid and suspended sediment: In this study, initial topography is arbitrarily set with a knickpoint. The knickpoint represents a boundary between a steep section of the channel and a gentle lobe deposit section.

Inverse analysis of turbidites Definition of the objective function
In this study, hydraulic conditions of turbidity currents were estimated from ancient turbidites based on optimization of input parameters of the forward model. The forward model of this study requires the parameters of the initial conditions (thickness H 0 , length l 0 , and sediment concentration C i0 ) as input, and produces the spatial distribution of sediment volume per unit area of each grain-size class as output. Then, the initial conditions are optimized to minimize the difference between the result of numerical simulation and outcrop data. The goal of optimization for the inverse analysis is to determine the global minimum of the objective function, J, which is defined as: where η ref i,k denotes the sediment volume per unit area of the ith grain-size class observed at the kth locality. This value is obtained by: where h k and F dep i,k denote deposit thickness and the fraction of the ith grain-size class at the kth locality, respectively. In addition, η calc i,k is the sediment volume per unit area of the ith grain-size class at the kth locality calculated using the Exner equation 4.
The objective function J is the sum of squares of deviations between the observation and the results of numerical calculations standardized according to the observed values. This function is calculated across all sampling points and grain-size classes, and it is optimized using the genetic algorithm, as described in the next subsection.

Optimization of the objective function
To obtain the optimized solution of the objective function, the genetic algorithm was employed. The genetic algorithm is a method for finding the global minimum (or maximum) of multivariate functions by mimicking the processes of biological evolution. The greatest advantage of his method is its efficacy for problems with a vast unknown solution spaces that could not be explored with other methods (Ramillien, 2001). In this method, the model parameter sets are regarded as genes of organisms, and optimizes them via the natural selection mechanism of evolution as follows: (1) An initial population of genes (parameter sets) is produced randomly in the parametric space. (2) Natural selection of genes (parameter sets) works based on the fitness value (the objective function). (3) Next generation of genes is created via crossing-over and mutation. (4) As the number of generations increases, the genes (parameter sets) in the population become optimized to the environment (observation). In this study, the Global Optimization Toolbox in MATLAB (Math-Works Inc.) was used for optimizing the objective function.

Application of the inverse model to a turbidite of the Kiyosumi Formation
This new method of inversion was applied to an individual turbidite bed in the Kiyosumi Formation on the Boso Peninsula, Japan. Here the geological settings and the characteristics of the bed are explained. Then, the results of the inversion analysis are presented.

Geological setting of the Awa Group and the Kiyosumi Formation
The Miocene to Pliocene Awa Group, which has been interpreted as deposited in a forearc basin-paleoenvironment, is distributed in the central part of the Boso Peninsula (Fig. 3). The Awa Group unconformably overlies the Paleogene Mineoka Group and is unconformably overlain by the Pleistocene Kazusa Group. The Awa group is subdivided into lower and upper parts (Ishihara and Tokuhashi, 2001). The upper part of the Awa Group is composed of the Amatsu, Kiyosumi and Anno formations, which are characteristically intercalated with a variety of volcanic tuff beds (Tokuhashi et al., 1997) (Fig. 4). The Kiyosumi Formation in the eastern to central region of the Boso Peninsula is composed of a sandstone-dominated alternating succession of turbidite sandstone and hemipelagic mudstone.
In contrast, the Kiyosumi Formation in the western region of the peninsula is composed mainly of mudstone or mudstone-dominated alternating turbidite sandstone and hemipelagic mudstone. The turbidites of the Kiyosumi Formation thin drastically westward, and the formation transitions into the Inakozawa Formation in the west (Nakajima et al., 1981;Suzuki, 1995;Ishihara and Tokuhashi, 2001). The maximum thickness of the Kiyosumi Formation is about 850 m (Tokuhashi, 1979;Tokuhashi and Ishihara, 2008). Sedimentological studies of ancient turbidites in the Kiyosumi and Anno Formation have been conducted in detail based on correlation of key tuff layers (Tokuhashi and Iwawaki, 1975;Tokuhashi, 1976a, b;Tokuhashi, 1988;Takahashi et al., 1997;Kubo et al., 1998;Ishihara and Tokuhashi, 2001;Saito and Ito, 2002;Shikazono et al., 2006). According to Tokuhashi and Iwawaki (1975) and Tokuhashi (1976a,b), many individual turbidite layers were correlated across a range of 40 km from east to west and 5 km from north to south based on these key tuff layers. The paleocurrent data recorded by the turbidite structure suggest that the channel mouth was located around Lake Kameyama (Tokuhashi, 1976a(Tokuhashi, , 1976b(Tokuhashi, , 1979. investigated the three dimensional morphology and the sedimentation processes of the turbidity currents, and suggested that the Kiyosumi Formation consists of deposits on lobes of submarine fans. The paleobathymetry of the Kiyosumi Formation is estimated to have been about 2000 m based on analysis of benthic foraminifera preserved in the hemipelagic mudstone (Hatta and Tokuhashi, 1984;Kitazato, 1987). Biostratigraphic investigations of planktonic foraminifera (Oda, 1978) and calcareous nannofossils (Kanie et al., 1991;Kameo et al., 2010), as well as magnetostratigraphic analyses, (Kimura, 1974;Niitsuma, 1976) have shown that the depositional age of the formation is early Pliocene. The absolute age has also been obtained based on fission-track dating of zircon grains (Kasuya, 1987;Okada and Bukry, 1980). Takahashi et al. (1997) summarized the biostratigraphic age determined from foraminifera, and estimated that the depositional age of the formation ranges spans about 800 thousand years from 5.1 Ma to 4.3 Ma.
This study is focused on one individual bed named the G1 turbidite by Tokuhashi (1976a). This turbidite bed is in the Kiyosumi Formation, and is intercalated between two key tuff layers (Tokuhashi, 1976a), and it can therefore be correlated over 20 km. Details of the characteristics of this G1 turbidite are described below.

Methods for field survey and analysis of deposits
The individual turbidite bed investigated in this study was identified at seven sampling localities in the research area (Fig. 3). The sampling localities are distributed in a fan-shaped region with a central angle of about 60 • and a radius of about 20 km. The channel mouth is estimated to be located around the apex of the study area.
In the field survey, route maps and geological columnar sections at 1/25 scale were logged at the correlated ca. 1 m intervals at the seven localities. Thickness, strike and dip of the G1 turbidite were measured, and samples were collected at regular intervals for grain-size analysis.
Bed thickness was measured with digital caliper (19975, Shinwa Rules Corporation, Niigata, Japan) for beds 0.01-15 cm thick. At all localities, samples for grain-size analysis were collected at 2 cm interval from the bottom to the top of the turbidite sandstones. Locs. 1 and 3, however, have space intervals of 2 cm between samples. Each sample was more than 20 g in weight, and about 4 g of sediment from each sample was used for grain-size analysis. Mud-sized particles (less than 62 µm in diameter) were removed by 250-mesh sieves (62 µm) with an ultrasonic bath before the settling-tube grain-size analysis. Later, dry-weight fractions of mud that were sieved out and recorded using an electronic scale. Then, grain-size distributions of sands in each sample were analyzed using the settling-tube method (Gibbs, 1974). This method can be used to measure grain-size distributions based on the representative settling diameter of each grain-size class without any measurements of the geometric properties of particles. The settling diameters of grains reflect their hydraulic behavior, which is significant for the purposes of this study. The settlingtube grain-size analyzer used in the study consists of a cylindrical flume 180 cm height (Faculty of Science, Kyoto University) and an electronic balance (UX820H, Shimadzu Corporation, Kyoto, Japan). The software "STube" was used for the measurements (Naruse, 2005). The compositions of grains of -1 to 5 phi were measured in increments of 0.2 phi. Finally, the fraction of mud in each sample was calculated from the summation of amounts of sieved-out mud and particles measured as 4-5 phi from the settling-tube analysis.

Examples of the forward model calculation
To test the forward model, two numerical simulations of turbidity currents transporting sediment of multiple grain-size classes were conducted. In this study, the multiple grain-size classes were set to three grain-size classes. The initial conditions employed in these numerical simulations were: initial flow thickness H 0 = 25 m, initial along-slope length l 0 = 25 m and the total sediment volumetric concentration C T 0 = 0.5%. The initial volumetric concentration of suspended sediment for the case of the multiple grain-size classes is assumed to be uniform. Other initial settings are shown in Table 4. In the case of the simulation using multiple grain-size classes, grain-size distribution of sediment is discretized to three grain-size classes as follows: very coarse sand, medium sand and very fine sand. As results of the simulations of multiple grain-size classes, the following time evolution of the turbidity currents were obtained (Fig. 5). In both cases of multiple grain-size classes, the thickness of the flow was thickest at the turbidity current head, and it thickened over time (Fig. 5A). The velocities of the current increased downstream in space, and decelerated gradually in time (Fig. 5B). The total volumetric concentrations of sediment also increased downstream in space, and decreased rapidly through time because of sedimentation (Fig. 5C). The sediment concentrations of all grain-size classes decreased over time. The very coarse sands settled out most rapidly, and very fine sands remained in the turbidity current even at the distal point. As a result, the fraction of very coarse sands in the final deposit was larger than that of very fine sands. The spatial distribution of deposit thickness showed two peaks at the upstream end of the calculation domain and at the knickpoint of the slope, with exponential decrease toward the downstream end (Fig. 6). The resultant deposit shows an upward-fining trend (Fig. 7A), as well as a fining-downstream trend (Fig. 7B).

Sensitivity tests of the forward model
Sensitivities of the forward model against the initial conditions and the model parameters were tested (Table 5). In these tests, numerical simulations were conducted, with six parameters (H 0 , l 0 , C T 0 , c f , e s and r 0 ) changed around the values that were used in the example described in the previous section. The other parameters used in these tests were not varied from those of the previous example.
The results of these sensitivity tests indicate that at the maximum reach of flows within the prescribed maximum time of calculation (2000 s), the average flow thickness, flow velocity, and sediment volumetric concentration varied greatly in response to changes in the initial conditions (Table 1), and therefore the spatial distribution of thickness (Fig. 8). Both the maximum extents and the total volumes of deposit increased as the initial flow thickness H 0 , the initial along-slope length l 0 and the initial concentration C i0 increased. The maximum extent of the deposits increased from 0.0016 m to 0.0600 m when H 0 increased from 5 m to 45 m (Table 5 cases 1, 2, and 3, Fig. 8A), and it varied from 0.0058 m to 0.0407 m when l 0 was changed from 5 m to 45 m (Table 5 cases 1, 4, and 5, Fig. 8B). In addition, variation of C T 0 from 0.1% to 1% (Table 5 cases 1, 6, and 7) resulted in change of the maximum extent of the deposits from 0.0017 m to 0.0804 m. Spatial variation of grain-size distributions of the deposits were also sensitive to the initial flow conditions. Compared with the initial parameters, the spatial distributions of both thickness and grain-sizes of deposits were less sensitive to the model parameters. When the coefficient of sediment entrainment e s was varied from half to double the value obtained from the original equation of Wright and Parker (2004) (Table 5 cases 1, 8, and 9), the maximum extents of deposits did not vary greatly (0.0246 m to 0.0340 m) (Fig. 8D). They varied from 0.0282 m to 0.0289 m when r 0 was changed from 1 to 2 (Table 5 cases 1, 10, and 11), which resulted in small changes in the maximum extent of the deposits from 0.0282 m to 0.0289 m (Fig. 8E). In addition, variation of c f from 0.001 to 0.05 (Table 5 cases 1, 12 to 15) resulted in small change in the maximum extent of the deposits from 0.0211 m to 0.0557 m (Fig. 8F).

Testing method of inverse analysis
A set of artificial turbidite data was produced by the calculations of the forward model using multiple grain-size classes. The spatial distribution of sediment volume per unit area of each grain-size class of the turbidite was calculated under the initial conditions shown in Table 4.
The topography of the objective function was examined around the initial conditions shown in Table 4, by taking sediment volumes of the artificial turbidite generated with the initial conditions taken as observations and calculating the function values relative to those generated under various other conditions. The resultant contour maps are shown in Fig. 9.
Inverse analysis of the artificial data then was carried out. Several parameters such as the number of generation and population size, are required for running the genetic algorithm, and the best settings for obtaining the smallest value of the objective function were investigated by testing various combinations of model parameters (Table 2). For both artificial data sets, 100 generations and a population size of 50 yielded lower values of the objective function. Other parameters are shown in Table 2. As results, the initial thickness H 0 , the initial length l 0 , and the initial sediment concentrations for the three size classes (very coarse sand; C 1,0 : medium sand; C 2,0 : very fine sand; C 3,0 ) were 24.68 m, 28.91 m, 0.156%, 0.149%, and 0.157%, respectively, when the value of the objective function was smallest (Table 2 Case III). In contrast, for the largest value of the objective function, these values were 25.19 m, 14.11 m, 0.226%, 0.228%, and 0.280%, respectively (Table 2 Case IV).

Characteristics of the G1 turbidite
This study is focused on an individual turbidite bed intercalated with the key tuff bed Ky 18 (informally named the G1 turbidite by Tokuhashi (1976a)). The key tuff beds Ky 1 to Ky 33 are also distributed in this study area (Nakajima et al., 1981;Kubo et al., 1998;Kameo et al., 2010). The tuff layer Ky 18 consists of a couple of a lenticular acidic tuff beds and a scoria bed intercalating a single turbidite and hemipelagic mudstone (Kubo et al., 1998). The Ky 18 tuff is intercalated in the upper part of the Kiyosumi Formation (Fig. 4), and extends across the distribution area of the formation. The lenticular acidic tuff bed is white in color and is composed of silt-sized fine volcanic glass. The bed ranges from 0 to 5 cm in thickness. The scoria bed is black, and is mainly composed of granule-sized to very course sand-sized particles. This tuff bed shows upward-finning trend, and has an average thickness of about 7 cm with a range of 4 to 10 cm (Fig. 10). Ishihara (2008), Nakajima et al. (1981) and Tokuhashi (1976a,b) correlated the G1 turbidite across 40 km based on the correlation of Ky 18. The G1 turbidite is intercalated in a succession composed of sand-dominated alternations of thick turbidite sandstone and hemipelagic mudstone. Sandstone beds range from tens of centimeters to several meters in thickness, and thins in the downstream direction. They are interpreted as the deposits of lobes of the sandy submarine fan ( Fig. 10; Tokuhashi, 1976a, b).
The G1 turbidite is composed mainly of medium sandstone with minor amounts of volcaniclastics, and it can be subdivided into four divisions based on its sedimentary structures: (1) the inversely graded division, (2) the normal graded division, (3) the parallel laminated division, and (4) the muddy division, from bottom to top. In the middle of the normal graded division, scoria and pumice are abundantly scattered, and in some cases concented layers of pumice grains are observed immediately below the parallel laminated division (e.g., Loc. 3; Fig. 11B). The bottom surface of the sandstone of the G1 turbidite (e.g., Loc. 2; Fig. 11E) shows some minor features of erosion (Loc. 2 and Loc. 5; Fig.s 11D and E). The thickness of the G1 turbidite thins downcurrent, and varies from 49.7 cm at Loc. 1 to 2.0 cm at Loc. 7 (Table 6 Fig . 11). The thickness of the deposit varies drastically between Loc. 3 (43.2 cm) and Loc. 4 (3.8 cm), which is consistent with the observations of Tokuhashi (1976aTokuhashi ( , b, 1979. The normal graded and the parallel laminated divisions correspond to the A and B divisions of the Bouma sequence, or to the characteristics of the K-type turbidites of the Awa Group, which were described by Takahashi et al. (1997). The mudstone overlying the G1 turbidite is interpreted to be composed of turbidite mudstone and hemi-pelagic mudstone (Takahashi et al., 1997). The results of the grain-size analysis show patterns of vertical and horizontal variations of grain-size distribution of the G1 turbidite (Fig. 12). In the vertical direction, the variation of grain-size distribution in the vertical direction represents the four divisions of the bed mentioned above, which are recognized at most of the sampling sites. In particular, the upward fining trends showing the coarse-tail grading are distinct in the normal graded division (Fig. 12). The term "coarse-tail grading" means a pattern of grading where coarser grains selectively decrease in proportion upward (Walker, 1975;Sylvester and Lowe, 2004). In the horizontal direction, the pattern of variation also indicates that the bed fines downstream (Fig. 12).

Inverse analysis of the G1 turbidite
In this application of the method of inverse analysis, the slope inclinations were set to around 20% for the upstream region and around 0.3% for the downstream region. The knickpoint between these two regions was located at the estimated position of the channel mouth (Fig. 3). The sampling localities were assigned estimated distances from the knickpoint to the sampling point along the downstream direction ( Table 6). The other properties of the sampling localities are shown in Table 6. The sediment volume per unit area for each grain-size class of the G1 turbidite at each location was calculated from measurements of bed thickness and grain-size distributions (Table 7). Grain-size distributions of the deposit were discretized into three grain-size classes (medium, fine, and very fine sands) for the inverse calculation (Fig. 13). The fractions of very coarse and coarse sands were very small in the G1 turbidite; therefore, these two grain-size classes were eliminated in the analysis. The results are shown in Fig. 13. The range of the initial parameters was set such that H 0 changed from 200 m to 500 m, l 0 changed from 200 m to 500 m, and C T 0 changed from 0.12% to 0.45%. Medium sand (C 1 ) changed from 0.01% to 0.05%, fine sand (C 2 ) changed from 0.01% to 0.2% and very fine sand (C 3 ) changed from 0.1% to 0.2%. The setting of the genetic algorithm is shown in Tables 8.
The hydraulic conditions of the turbidity current were then reconstructed using the inverse analysis procedure described above. The final value of the objective function was 9.73. Initial conditions obtained from the inverse analysis were as follows: flow thickness H 0 = 441.09 m, flow length l 0 = 200.00 m and sediment concentration C T = 0.3030% (medium = 0.01, fine = 0.13455, and very fine = 0.15841).

Behavior of the reconstructed turbidity current
The behavior of the turbidity current that emplaced the G1 turbidite can be estimated by calculating the evolution of hydraulic conditions over time based on the initial conditions obtained from the inverse analysis. The flow thickness, flow velocity and sediment concentration at the downstream end of the sampling area (Loc. 7, about 20 km from the channel mouth) were reconstructed as 334.55 m, 0.98 m/s, and 0.0058%, respectively, when the turbidity current reached the point (Table 9). Furthermore, the flow thickness, flow velocity and sediment concentration at the downstream end after a very long duration (60,000 seconds) were found to be 158.84 m, 0.28 m/s, and 0.0028% respectively.

Validity of the forward model
The layer-averaged model used in this study has been frequently used to simulate field-scale turbidity currents, and is known to reproduce actual morphodynamic phenomena well (e.g. Fildani et al., 2006). Validity of the forward model for calculating field-scale turbidity currents is suggested by the similarities in geometry and grain-size trends of deposits between the results of the simulations and the actual individual turbidite beds correlated over long distances. This forward model produces a turbidite showing a peak in thickness located around the proximal region of the slope ( Fig. 6D and 8). Thickness of the calculated turbidite rapidly increases downstream around a knickpoint, and then gradually decreases downstream. These features resemble those of typical turbidite beds of the Boso Peninsula as summarized by Tokuhashi (1976b). Lowe (1982) also proposed a hypothetical geometry of sandy, high-density turbidites, which shows similar geometry to the turbidites calculated by the forward model. Although turbidites predicted by this forward model show another peak in thickness on the proximal side of the slope where the initial suspended sediment cloud was located (e.g. Fig. 6D), this peak results from an artificial effect caused by employment of the "lock-exchange" initial conditions (Fig.  2). Focused on the downstream side of the knickpoint, the features of the field-scale turbidite are well reproduced by the forward model. It should be noted that the region of the steeper slope is strongly affected by the artificial settings; therefore, the results of calculations for this region were excluded from the objective function calculations.
The features of turbidites predicted by the forward model also show great similarity in both the patterns of vertical and downstreamward fining compared with those of actual ancient turbidites (Fig. 7). The graded bedding associated with the A division of the Bouma sequence is a typical feature of ancient turbidites. Middleton (1966) reported this pattern of grading, called coarse-tail grading, from experimental turbidites, and it has been proved that this grading pattern is a near-universal features of turbidites (Lowe, 1982). The results of this study also indicate that deposition of coarse-grained sands is concentrated around both the proximal side and the initial stage of turbidite formation. In contrast, finer sands tend to be deposited in the broad region from the proximal to distal ends, and deposition continues from the beginning to the end of the flow's duration. These characteristics of turbidite deposition induce the graded bedding that shows coarse-tail grading features, which is especially common in coarse-grained sandy turbidites.
The sensitivity tests of the forward model also show the adequacy of the model for the inverse analysis of turbidites. The results of the sensitivity test revealed that the forward model used in this study is sensitive to the input parameters (Fig. 8), but is less sensitive to the model parameters (Fig. 8). It would be difficult to determine the optimal solution from the result of the model calculation if the variation in the input parameters does not significantly change the results of the forward model calculation. Conversely, if the results are sensitive to the input parameters, even small differences in the initial conditions may be detected from the variation of deposits. Therefore, we infer that it is generally feasible to conduct the inverse analysis for ancient turbidites using the forward model proposed in this study. In contrast, unlike the input parameters, the selection of the empirical parameters of the model did not affect the calculation results significantly. Although, there are many possible choices of empirical formulations proposed in previous studies such as the dimensionless coefficient of the rate of sediment entrainment e s , this model is robust with arbitrary selection of model parameter.

Validation of the methodology of inverse modeling
The results of tests of the inverse analysis method using the artificial data indicated that the optimization calculation method adopted in this research can adequately reconstruct hydraulic conditions of turbidity currents from turbidite deposits. This method reconstructed the initial conditions of flow as H 0 = 24.68 m, l 0 = 28.91 m, and C T 0 = 0.462% ( Table 2). The true values are H 0 = 25 m, l 0 =25 m, and C T 0 = 0.5%, respectively. Differences between the solutions of the inverse analysis and the true values were within approximately 16% (1.28% for H 0 , 15.64% for l 0 , and 7.6% for C T 0 ). Therefore, the method proposed in this research is inferred to also suitable for estimating the paleo-hydraulic conditions of actual turbidity currents associated with turbidites in the geological record. The suitability of parameters of the genetic algorithm, for inversion using the forward model employed in this study were also explored. The results indicate that, 100 generations and a population size of 50 provided best solution in the case of the artificial data examined in this study. Furthermore, reasonable results were also obtained using a wide range of parameters for the genetic algorithm, which suggests that our method does not strongly depend on the selection of settings for the optimization method. These lines of evidence suggest that reasonable estimates of hydraulic conditions of turbidity currents can also be obtained from actual ancient turbidites using this method.

Advantages of the method proposed in this study
Compared with the methods proposed in previous studies, the present method of inverse analysis is at a great advantage for reconstructing the paleo-flow velocities of ancient deposits. Previously, researchers derived the flow velocities of turbidity currents based on the critical velocity of suspended sediments . Such a method, however, can only be used to estimate the minimum flow velocity, and pointed out that the results are quite different from the actual paleo-flow velocities of turbidity currents. In contrast, recently, several methods of inverse analysis based on optimization calculations of input values for numerical models have been proposed . Compared with these previously proposed methods, the present method has the following advantages: (1) This study employed a non-steady flow model, which can reproduce the evolution of turbidity currents over time. proposed the methodology employing a steady flow model as the forward model. As mentioned in the previous section, one of the characteristic features of turbidites is graded bedding, and this structure can be formed only by non-steady flows.
(2) The model of this study considers both sediment entrainment (resuspension) from the bed by turbidity currents and turbulent mixing of suspended particles in the flows, which are essential processes for representing suspended sediment transport. estimated the paleo-flow velocity of the turbidity current using a "box" model, giving an arbitrary value of the initial flow thickness. However, their model did not account for sediment resuspension processes. also proposed an inverse analysis method based on the nonsteady layer-averaged model of turbidity currents, but they also did not consider resuspension processes, which could create problems. (3) Calculation costs of this model are low enough to apply it with field-scale data (calculation time was 319,530 seconds). proposed an inverse analysis method based on direct numerical simulation (DNS) of the Navier-Stokes equations. However, the calculation costs of their method are extremely high, and consequently it is difficult to apply their method to the field-scale data from ancient turbidites.

Characteristics of turbidity current deposited an ancient turbidite of the Kiyosumi Formation
The present method was applied to an individual ancient turbidite of the Kiyosumi Formation, the G1 turbidite, and the hydraulic conditions of the turbidity current that emplaced the G1 turbidite were reconstructed. For application of the method proposed in this study, the G1 turbidite is appropriate for the following reasons: (1) correlation of the bed is reliable because it is intercalated with easily recognizable tuff beds (Ky 18), which can be traced 30 km from east to west and 5 km from north to south; (2) the G1 turbidite shows no significant erosional structures . In this study, the estimated flow thickness, velocity and sediment concentration were 334.6 m, 0.98 m/s, and 0.0058%, respectively (Fig. 14), when the turbidity current arrived at the sampling point at the downstream end of the research area (Loc. 7). applied the "box" model of for the G1 turbidite and estimated a paleo-flow velocity 5.56 m/s. The velocity result from the present method is lower . The validity of these reconstructed values will be discussed in detail in future works, but several in-situ observations are available for comparison with our inversion results. For example, following an earthquake under the Grand Banks in 1929, the velocity of a turbidity current was estimated from the intervals between cable breaks, and this velocity reached high as 7.7 m/s . In addition, when a storm passed over Scripps Submarine Canyon in 1976, the velocity of a turbidity current was estimated to be as high as 1.9 m/s at a depth of 44 m . During a hurricane offshore of Oahu near the Kauai Channel in 1982, turbidity currents velocities were estimated to be as high as 3.0 m/s . In Monterey submarine canyon offshore of California, turbidity currents were observed using an acoustic Doppler current profiler . The measured velocity of the currents ranged from 1.7 to 2.8 m/s. reported turbidity currents in the Congo Fan channel of the Congo Canyon, and recorded turbidity currents in the lower reaches of the canyon . showed that a turbidity current was triggered by the tsunami caused by earthquake that occurred offshore of Tohoku-Oki in 2011 and estimated flow velocity as 2.4-7.1 m/s. In addition, presented that multiple surge-like turbidity currents were examined using a new imaging method in the mouth of the Squamish River in Howe Sound, British Columbia, Canada and the flows varied in speed from 0.5-3.0 m/s. Compared with these measurements, our results indicate relatively low velocities and high flow depths. Future improvements in both the in-situ measurements and forward model development may solve these discrepancy in the estimated values.

Conclusions
We propose a new method of inverse analysis of turbidity currents based on turbidite deposits, and applied this method to an individual ancient turbidite in the Kiyosumi Formation. Despite their significance in the paleoenvironmental research and resource geology, the flow properties of turbidity currents in deep-sea environments remain uncertain because in-situ measurements are made difficult by the highly destructive nature and infrequency of these occurrences. Therefore, to better understand the behavior of actual turbidity currents, we have aimed to develop a new method of inverse analysis to reconstruct the paleo-hydraulic conditions of turbidity currents from ancient turbidites. There have been a few privious studies of inverse modeling of turbidity currents; however, several problems with these studies have been pointed out. They employed an oversimplified forward model that is not suitable for reproducing typical features of ancient turbidites, or else the calculation costs of their method were too high to apply to field-scale data.
To solve these problems, we have developed a new forward model of non-steady turbidity currents that accounts for mixed grain-size sediment, and can describe the behavior of a turbidity current that deposits a typical turbidite with graded bedding. The new model employs one-dimensional shallow water equation, which is applicable to field-scale problems. The "lock-exchange" type condition is assumed as the initial setting in this model. For inverse analysis, the objective function is defined as the sum of squares of deviations between the sediment volumes of the observations and the numerical calculations. In the present inverse calculation, the initial hydraulic conditions that minimize the objective function are explored with a genetic algorithm. Tests of our inversion method using artificial data provided reasonable results, which suggest adequacy of this optimization methodology. We then applied our method to field data from a turbidite in the Kiyosumi Formation, Boso Peninsula, Japan. The Kiyosumi Formation is sand-dominated and composed of alternating of turbidite sandstone and hemipelagic mudstone, which are interpreted to be deposits of a submarine fan lobe. In this study, an individual turbidite bed intercalated between the two key-tuff layers was correlated over 20 km, and thickness and grain-size distribution of the bed were measured at an seven sampling localities. Through the inverse analysis, the hydraulic conditions of the turbidity currents that had emplaced this turbidite were estimated. The flow thickness, velocity, and total sediment concentration were reconstructed to be 334.55 m, 0.98 m/s, and 0.0058%, respectively, when the flow reached the downstream end of the sampling area. Although verification of these results will be discussed further in the future, these reconstructed values are in agreement with the hydraulic conditions of turbidity currents monitored in previous studies. However, room for improvement remains for the forward model used in this study such as design of the initial and boundary conditions, and optimization methods including the objective function.

Availability of data and material
All source codes used in this paper are available at the second author's web site [http://turbidite.secret.jp/turbidite.secret.jp/turbinversion/].

Competing interests
The authors declare that they have no competing interest.

Funding
This work was supported by the Japan Society for the Promotion of Science, KAKENHI Grant Number JP26287127.

Authors' contributions
HN proposed the topic, conceived and designed the study. NK carried out the field survey, NK and HN coded the forward and inverse model, and NK conducted numerical analysis of data set. ST helped field survey and correlation of beds. All authors read and approved the final manuscript.

Figure legends
Figure 1 Schematic diagrams of processes considered in the forward model presented in this paper. Water entrainment from ambient water, suspended sediment transport with entrainment from and settling to bed sediment, and the active layer on the surface of beds are considered.

Figure 2
Schematic diagram of a turbidity current in transformed coordination and lock-exchange condition. Position of the head of the turbidity current is denoted asŝ and that of the tail is fixed at x = 0. Initial conditions of the suspended sediment cloud are defined by the initial thickness H 0 , the initial along-slope length l 0 and the initial sediment concentration C i of the ith grain-size class. The sediment cloud is situated on a slope, assuming that the downstream end of the cloud is confined by the virtual lock gate until t = 0.      Table 5 for details of simulation settings. Figure 9 Contour maps of values of the objective function in test cases with multiple grain-size model. Log10-transformed values of the objective function describe deviations of calculated results from the artificial data of the turbidites generated by the forward model with a given set of initial conditions. The true initial conditions are shown in Table 4. Steep minima are observed around the true initial conditions. A. The initial concentration is constant. B. The initial thickness is constant.        1, 600 * 2, 000 C 1,0 ; very coarse sand: C 2,0 ; medium sand: C 3,0 ; very fine sand.

Tables
* : The duration of the case of the single grain-size class was 1,600 seconds because the terminal condition of sediment concentration was reached.