Flood Discharge Changes by Urbanization from Farmland Based on the Unit Flood Concept at the Kurabe River, Japan ()
1. Introduction
Because of the rapid progress of urbanization from farmland since the 1960s, flood discharge has caused serious damage by increasing due to urbanization because of a shortage of river and canal capacities. This resulted from land use changing from paddy fields to the residential areas. Therefore, increasing canal capacity and establishing a small reservoir to control peak flood discharge were frequently undertaken.
However, an estimation method for flood discharge has not been established, which is the basis for this research. We recently suffered very severe flood damage from an extremely local heavy storm. Therefore, a precise local runoff analysis, not generalized, is strongly required. This research addresses these problems.
Increasing flood discharge from urbanization has been studied by many researchers from various perspectives. An outline of the results is summarized by Kadoya [1] , who discusses “the changes in flood discharge from urbanization” since 1985 based on worldwide reports. The contents of that manuscript are divided as follows: a definition of urbanization, the change in precipitation holding capacity, changes in the hydrograph, and changes to runoff characteristics. In which he searched about one hundred reports for the city runoff analysis. Among them he discussed many research on runoff analysis based on the precise aspect i.e., physical properties of the basin. For example, Hicks (1944) and Tholin (1959) researched initially the urban runoff [2] [3] . Kadoya (1972) conducted predictive study on urbanization effect on the discharge [4] . These research is mostly limited the runoff analysis inside of city not from land use changes. On the other hand, Eagleson (1962) and Kar (2015), for example, the flood changes evaluated by unit graph changes by urbanization [5] [6] .
After the Kadoya’s report, we searched about sixteen reports by urbanization [7] [8] [9] relating the land use changes from farmland to residential areas. From the aspect, we did not found the appropriate report focusing on estimating the runoff changes from farmland to residential area. We believe that the estimation of runoff changes requires fundamentally the physical basis of the basin, i.e., type and area of paddy, upland and residential lots etc., should be taken into considered for the runoff model.
However, Watanabe and Toyokuni [10] and Toyokuni [11] [12] undertook a runoff analysis based on land use characteristics and facilities by urbanization, which was clearly based on the physical properties of the basin but not from farmland.
On the other hand, many runoff models have been developed worldwide, such as the unit hydrograph method, reservoir storage function model, tank model and kinematic wave model. However, these models except for kinematic wave model for residential areas has a limitation for estimation of runoff discharge by urbanization, because of less consideration of the physical properties of the basin i.e., less of the hydraulics aspect. To estimate the flood discharge changes from urbanization, the physical properties should be included in the model as mentioned above; if not, the runoff changes for future urbanization can’t be estimated. Therefore, our model is focused on strongly the changes of those properties from farmland. That is a remarkable feature of our research.
Namely, previous research discussed the unit flood discharge [7] [8] , defined as the discharge that flows in the upstream drainage canal, i.e., the smallest drainage canals. This model was developed using the unit flood discharge concept for various land uses, such as paddy fields and upland and residential areas, and it attempts to estimate the basin discharge, including the effect of a small flood control reservoir. The discharge changes after 1976 due to ur-banization was estimated in the Kurabe River Basin in Japan. The validity of the model was confirmed by comparing observed and estimated hydrographs in the test basin.
In addition, peak discharge by urbanization closely related not only urbanization but also rainfall intensity. However, this research is limited only runoff discharges by urbanization.
2. Method
2.1. Method of Analysis
The Kurabe River Basin was located in the center of the Tedori River Alluvial Fan Area, which is 17.50 km2 in area and is not contain a mountainous area and inside of previous research area [7] [8] [9] . The drainage canal plays an irrigation role in the upstream area, which is very rare in Japan. The upstream canals have a concrete lining, and the downstream area is a natural river. The gradient is very steep, at 1/500, as is the Tedori River Alluvial Fan Area. Therefore, this area has not suffered the backwater effects, thus, the kinematic wave approach for the runoff analysis can be available. This basin has not included mountainous area and the urbanization conducted from mainly paddy field. This feature should be notice for applying the method to the other basin.
2.2. Change in Land Use
The basin was dived into 15 blocks based on the drainage canal system, as shown in Figure 1. Table 1 describes land use in the area of the 15 blocks in each of the tested years. The area was investigated using GIS (Geographical Information System) with 100 m mesh data (National Land Numerical Information download service) [13] . The residential area contains an area of river and road that has the same runoff characteristics. The areas that did not include 100 m mesh data were corrected with data coinciding with the Land Use Investigation at Hakusan District in Ishikawa Prefecture [14] . The available data for land use changes were for the years 1976, 1989, 1991, 1997 and 2009.
2.3. Tested Precipitation
The 10-minute interval precipitation data from the Kanazawa Meteorological Branch was selected as the test events. The test data were chosen from over 80 mm precipitation data during 24 hours, resulting in 21 events selected during 2009-2014. The feature of the tested precipitation, i.e., occurred year, beginning and end time of the event, the amount of the total precipitation and the peak
![]()
Figure 1. The Kurabe River basin and sub-basin including the study area.
precipitation, describe in Table 2.
2.4. Unit Flood Discharge for Various Land Uses and Discharge in the Individual Blocks
The unit flood discharge, differ from the unit hydrograph method which is very popular approach for runoff analysis, was a new concept that was developed by the authors as in previous research [7] [8] for clearly estimating basin discharge by urbanization from farmland. The unit flood discharge is expressed in mm 10 min−1 unit. Therefore, by multiplying this by the related land use area, the flood discharge can be obtained for the smallest drainage canals.
Using the unit discharge for the individual land uses, the flood discharge in the individual blocks can be calculated as follows. Which is weighted average of land use area with our unit flood discharges in a block [7] [8] :
![]()
Table 1. Various land use changes in individual blocks and land use areas in the tested years (km2).
*Increase ratio of residential (IROR: Area of residential (2009)/Area of residential (1976).
![]()
Table 2. Feature of tested precipitation events.
(1)
Here, i is the land use code (1: Paddy, 2: Upland, 3: Residential), qt is the block discharge (mm 10 min−1), t is time step, A is the block area (km2), ai is the area of land use code i (km2), Ft,i is the unit flood discharge in an individual land use code i (mm 10 min−1). Which is not unit graph and mathematical function but schematic expression of unit flood discharge in previous reports [7] [8] .
Furthermore, the relationship of block specific discharge Qt (m3・s−1・km−2) and qt (mm 10 min−1) is Qt = qt・0.6−1.
2.5. Flood Routing of Net Drainage Work
The unit flood discharge tends to decrease with flow down from the field, in the individual blocks, in the trunk canals and in the river basin. Therefore, the basin discharge should be estimated by flood routing from the unit discharge to the river mouth following the drainage network. The flood routing can be applied with a kinematic wave method [15] because of the very steep landscape.
The kinematic wave method, which is well known, is briefly described in the following section.
2.5.1. Flood Routing Method for a Drainage Network
The fundamental formula for flood routing is as follows:
When
:
, 
,
(2)
Here, x is the distance from the upstream drainage canal, Q is the discharge in the drainage canal, q is the lateral inflow per unit length for the drainage, t is the time, and b is the distance from the upstream to the end of the drainage canal.
The relationship between the cross-sectional area of the drainage canal W and discharge Q can be expressed as W = KQP. Here, K and P are the constants of the cross section.
The calculation can be performed as a difference scheme by dividing Δt as follows:
When
:
,
(3)
When
:
,
(4)
Here, j is the order of the time increment, ∆x is the routing distance of q(t) between the ∆t, and qj is the lateral inflow discharge between the time from j − 1 to j (∆t).
When the sum of the routing distance x is over the b, after obtaining b − xm = ∆xm, setting j = m, the end of the discharge of the canal Qm and the time Δtm is obtained as follows:.
In cases where
:
,
(5)
In cases where
:![]()
,![]()
(6)
Furthermore, at the crossing point of the drainage canals, the total discharge can be obtained as the sum of the discharge at the crossed canals.
2.5.2. Model of Land Consolidation and a Block Drainage System
The field lot system principally originated from land consolidation of paddy fields. Therefore, the runoff routing should start from this system. The system consisted of a long side of 30 m × 10 lots, equaling 300 m, and a short side of 100 m. Usually, both sides of the canal have a lot; thus, the short side of the system is 200 m. Therefore, we found that a 200 m × 300 m system should be considered the model of the smallest field lots, as shown in Figure 3. Furthermore, the average land use in the system was applied by the land use in the block.
In a block system, we consider flood routing to move from the field system to the basin drainage systems. Between the systems, flood routing in the second and third drainage systems was applied in a standardized system. The second drainage system receives discharge from field systems at each 200-m interval. The canal length is 500 m because it collects the drainage from three field systems.
![]()
Figure 2. Drainage system in the Kurabe River basin.
Following this system, the third drainage system collected the discharge in 300-m intervals from the secondary drainage canals. This led the discharge to the trunk canal. Based on this standardized system, the required flow time down in the block was estimated (Figure 3).
At the trunk canal system, the drainage was estimated according to the drainage network. The flood routing in the system was performed without lateral inflow but collected the block drainage discharge at a given point.
The discharge analysis involved 1 representative event that showed relatively large amounts and specific discharges from 21 events. This event is 127 mm [event No. 3-1].
2.5.3. Effect of Flood Control by a Small Reservoir in the Block
The outlet point of newly developed residential area was set as a small reservoir
![]()
Figure 3. Land consolidation drainage unit and drainage unit in the block.
to avoid increasing the discharge load for a trunk canal or river. This aimed to control the peak discharge originating from urbanization that occurred from land use changes from paddy fields to residential areas. The reservoir was established as follows: the capacity of the reservoir should hold the discharge difference between before and after residential area development. It should be stored for at least over one hour after discharge [16] . The criteria are based on the concept that the reservoir should temporally store the increased discharge from the development of the residential area.
A popular reservoir is this area is shown in Figure 4 as an example. The reservoir was set in an inlet along drainage canals that serve as a spillway i.e., lowering part of the drainage canal banks. The outlet capacity of the stored water has approximately 40 cm of squire box. The controlled discharge flowed out to the down streams of the same canal, through the outlets.
To evaluate the effect of the reservoir, a water balance calculation should be conducted by applied rainfall event, type of reservoir and initial conditions. At
![]()
Figure 4. Structure of flood control reservoir.
the latter section, to this end, the effect of the flood control reservoir was conducted by using the standardized reservoir as an example.
On the other hand, to more simply evaluate the effect of the reservoir, the average discharge was applied to estimate the reservoir effect because the above method is more complicated and depends on the type of reservoir and initial conditions.
2.5.4. Estimation of Flow Time and Discharge
The estimation of discharge from the field to the block drainage system can be performed using Equations (2)-(5), which have a lateral inflow. The flood routing of the trunk canal performed without lateral flow is described further. In the process, the lag time (peak time precipitation to peak discharge) was also estimated.
According to the canal network (Figure 2), a start time is chosen. First, the flood routing of the field drainage system occurs; second, the block routing of the drainage network occurs by collecting the field drainage; and third, the trunk discharge collection of block discharge occurs with a time difference. There is a rare case when the latter discharge overcame the former discharge. Furthermore, if the discharge at the top of the smallest canal in the field drainage system is zero, the flood routing was not conducted; therefore, a very small discharge of 0.001 m3・s−1 was applied.
In addition, the coefficient of the drainage canals K and P was determined for flood routing by actual cross sectional area and water depth,. This was performed using the roughness coefficient of a concrete lining canal of 0.014 with the gradient of 1/500 and a natural river of 0.030 with the gradient of 1/1200. The cross-sectional constant of K and P are shown in Table 3.
![]()
Table 3. Constants of cross sections K and P.
Note: T: change points of waterway cross-section (Figure 2).
3. Result
3.1. Change of Flood Discharge Unit in a Block
In order to clearly show the changes of discharge due to urbanization, we focus on the flood specific discharge of each block of different urbanization and summarized the calculation results (Table 4).
The flood discharge unit in a block for all events described at three points, which is the whole Kurabe River Basin, indicated the most urbanized block (3-2-2) and the slowly urbanized block (3-1-2). The unit discharge was estimated using Equation (1) (Figure 5).
For the Kurabe River basin, the maximum specific discharge in the 21 events changed from 6.05 m3・s−1・km−2 in 1976, to 12.3 m3・s−1・km−2 in 2009, resulting in the ratio increasing 2.02 times. Concurrently, the mean maximum specific discharge increased from 3.13 m3・s−1・km−2 (1976) to 6.11 m3・s−1・km−2 (2009). The increased ratio was 1.95 times greater, which is quite similar to the maximum.
For the most urbanized block (3-2-2), with a rate of 10.2 m3・s−1・km−2 (1976) to 18.7 m3・s−1・km−2 (2009), the resulting the ratio was 1.84 times greater. When the mean average specific discharge changed from 5.08 m3・s−1・km−2 (1976) to 9.21 m3・s−1・km−2 (2009), the resulting ratio was 1.81 times. The reason for the small ratio was that the discharge in 1976 was already urbanized.
In the slowly urbanized block (3-1-2), with a rate of 4.47 m3・s−1・km−2 (1976) to 4.83 m3・s−1・km−2 (2009), the resulting ratio was 1.08 times greater. Concurrently,
![]()
Table 4. Maximum specific discharge in the individual blocks.
Note: Area: km2, Maximum specific discharge: m3・s−1・km−2. *Increase ratio of specific discharge (IRSD): Maximum specific discharge (2009)/Maximum specific discharge (1976).
![]()
Figure 5. Maximum flood discharge changes due to urbanization.
the mean average specific discharge changed from 1.90 m3・s−1・km−2 (1976) to 2.24 m3・s−1・km−2 (2009), resulting in a ratio that was 1.17 times greater.
As mentioned above, from the perspective of maximum specific discharge in 2009, the rapid urbanization block (3-2-1) exhibited 21.7 m3・s−1・km−2, whereas the slowly urbanized block exhibited a ratio of 4.83 m3・s−1・km−2, resulting in a ratio that was 4.49 times greater. Even as an example, we determined that urbanization induced quite a large increase in the flood discharge.
3.2. Maximum Specific Discharge per Block
The maximum specific flood discharge in the individual blocks for all events described in Table 4. The results were represented as increasing ratio of specific discharge (IROS) for a unit area in the lower row. Most high IROS was 3.26 (K-3 block), and the maximum specific discharge was 21.7 m3・s−1・km−2 (3-2-1 block), which was remarkably large.
3.3. An Example of Hydrograph Changes from Urbanization in the Basin
Figure 6 describes an example of the hydrograph changes from urbanization in the basin. The hydrograph was obtained from a total precipitation of 127 mm [event No. 3-1]. Following the urbanization trend, this figure clearly indicates the change in the hydrograph, where the peak discharge changes from 88.3 m3・s−1 in 1976, 125.2 m3・s−1 in 1991and 161.7 m3・s−1 in 2009. This indicated approximately 1.83 times the increase from 1976 to 2009.
3.4. Relationship between the Unit and Basin Flood Discharge
The change of the hydrograph from the smallest canal to the river mouth is of great concern. Figure 7 describes the changes from a discharge unit hydrograph to a basin hydrograph as an example [event No. 3-1]. The figure indicates that the peak unit discharge decreased according to the downward. The upper figure describes the feature at some block and lower figure describe the flow for all basin to the river mouth. The peak discharge for block 3-1-2 changes from 9.18 m3・s−1 to 6.67 m3・s−1. The peak discharge for all basin of 172 m3・s−1 of the hydrograph unit decreased to 161 m3・s−1 of basin discharge with a time difference of 0.33 hours.
4. Consideration and Discussion
4.1. Statistical Consideration of the Discharge Unit
Based on the 21 tested events, statistical consideration was applied simply. We used the mean values and standard deviation employed as a statistical index and a normal distribution was assumed because the sample size was limited. The probability of the mean value m + σ over was 68%, and that of m ± 2σ over was 4.5%. Here, σ is the standard deviation.
Based on this concept, the mean value is m and m + σ and m ± 2σ are shown in Figure 5. The figure describes how the large flood discharge over the past maximum discharge events will occur in the future at approximately a 4.5% probability, even under the present land use conditions.
4.2. Effect of a Small Reservoir in a Block
After the 1990s, the establishment of a small reservoir in a block was required to control increased flood discharge by residential area developments. The criteria described are mentioned above. According to the criteria, we assumed the reservoir capacity and tested the effect as an example.
The criteria were initially determined in 2005. Thus, the standard of flood
![]()
Figure 6. Changes in hydrographs due to urbanization.
![]()
Figure 7. Comparison of unit and estimated discharge hydrographs.
![]()
Figure 8. Effect of flood control in the small reservoir (m3 sec−1).
discharge control data before 2005 was 1991 in our research. Therefore, the maximum discharge in the 1991 event was selected as the flood control boundary. A representative precipitation event [event No. 3-1] in this period was chosen because the event represented the maximum discharge. The relationship between the discharge of controlled versus uncontrolled discharge by the reservoir is shown in Figure 8 which is based on 2009 land use data for the whole
basin.
![]()
Figure 9. Comparison of flood control effects by reservoir and moving average (m3 sec−1).
The effect of the reservoir is described, in that uncontrolled discharge was 12.55 m3・s−1, whereas controlled discharge was 5.08 m3・s−1. The ratio was 0.40 times, which is remarkably large as shown in Figure 8.
To simply evaluate the effectiveness of the reservoir, an hourly moving average was applied for the estimated discharge before reservoir establishment. Figure 9 presents the relationship between our method after reservoir control and
that with a moving average applied. This indicates a very close relationship, with a value of 0.912 (R2 = 0.916). As indicated by this figure, the moving average of the discharge was useful for roughly and simply evaluating the effect of the small reservoir.
4.3. Comparison of the Peak Discharge from Our Method and the Rational Formula
Capacity of drainage canals as determined simply by the “Rational Formula” using hourly precipitation. This capacity does not consider principally the increase in flood discharge from urbanization. Therefore, it is of great importance to determine whether the capacity of the drainage canal at present is sufficient or not sufficient for urbanization.
The discharge was determined by the “Rational Formula” using planning precipitation average during the lag time estimated by the “Rziha formula”. The “Rational” and “Rziha formula” are presented as follows:
![]()
(7)
Here, Qm is the flood discharge, fp is the peak runoff ratio, rp is the peak precipitation in the lag time (mm・h−1), A is the basin area (km2), T is the lag time (h), L is the distance of the drainage channel (km), W is the velocity (km・h−1),
![]()
Figure 10. Comparison of peak discharge estimated by the “Rational Formula” and our methods. Note: ●: Rational formula (Tp: Rziha formula)-Proposed method (10 min), ○: Rational formula (Tp: Rziha formula)-Proposed method (60 min), ×: Rational formula (Tp: value of Proposed method)-Proposed method (60 min).
and H is the height difference of the channel (km).
The tested event was chosen as [event No. 3-1], which had the maximum discharge. The above “Rziha formula” was applied to estimate the flood discharge at 15 points in the Kurabe Basin. On the other hand, the comparative discharge of our method was estimated with the same points using 1976, 1991 and 2009 land use. The results in Figure 10 describe individual tested years.
The figures describe that the planning discharge in the “Rational Formula” was smaller 1/(2.19 - 2.66) than our method, which used a one-hour average precipitation (plot: “○”). Furthermore, the discharge was smaller 1/(3.04 - 4.11) than in our method, which used 10 min precipitation (plot: “●”). This fact indicates that the probability of flooding damage rapidly increased recently from urbanization because of the shortage of canal capacity. If we use the lag time estimated our method in the Rational Formula, the peak discharge indicate quiet similar values [1/(1.02 - 1.09)]indicated as Figure 10 (plot: “×”).
To overcome this problem, the improvement of the rainfall intensity constant in the Rational Formula runoff coefficient, return period and improvement of lag time estimation should be considered.
4.4. Comparison of Lag-Time between Previous Research and Our Method
The problem of lag time is very important for estimating peak discharge because it has a close relationship with the estimation of rainfall intensity related to discharge.
The lag time estimated by our method, in the event [event No. 3-1] was 30 min. The event produced a maximum peak discharge of 5.52 mm 10 min−1 (9.21 m3・s−1・km−2) for the whole basin area. If the intensity was converted to one hour using a relationship where the intensity is directly proportionate to the square root of the precipitation duration time, this results in approximately 33 mm・hr−1.
In contrast, Kadoya and Fukushima [17] proposed the following experimental formula based on the Kinematic Wave Concept:
![]()
(8)
Here, C is a constant, 60 - 90 are urbanized area, A is the basin area (km2), and re is the effective precipitation (mm・hr−1).
In the above formula, using A = 17.5 km2, re = 33 mm・hr−1, and C = 60 results in Tp = 33 min and C = 90, for a total of Tp = 49 min. Comparing these with the lag time as mentioned above, C = 60 - 90 is adequate. This result appears a tendency of urbanization is in progress at the basin from center to under in the Kurabe River.
4.5. Relationship between Residential Area and Specific Discharge in the Blocks
Figure 11 describes relationship between the ratio of residential area and specific discharge in the 15 blocks. The residential area ratio expressed as IROR (Residential area in 2009/Residential area in 1976 as defined in Table 1) and the specific discharge ratio expressed as IRSD (Maximum specific discharge in 2009/ Maximum specific discharge in 1976 as defined in Table 3). This figure clearly shows that IRSD increased remarkably with IROR. This fact indicate that the relationship between residential area and specific discharge has not only the relationship itself but also the increasing ratio of residential area and specific discharge.
5. Validation of Our Method by Test Basin
5.1. Outline of the Study Area
To verify the reasonability of our proposed method, a study area was established
![]()
Figure 11. Relationship of increasing ratio between residential area and specific discharge. Note: IROR: Increase ratio of residential (Table 1), IRSD: Increase ratio of specific discharge (Table 4).
on the inside of the Kurabe Basin, as in Figure 1, which had a 3.16 km2 area. Upstream and downstream of the study area was used water gauge to observe the discharge by the Shichika Water District. The drainage channel had concrete lining, thus the estimated discharge exhibited relatively high accuracy. The land use at the study area is described in Table 1 (3-1-2, 3-1-3). The total area was 3.16 km2, in which paddy fields covered 2.32 km2, upland covered 0.06 km2, and residential areas covered 0.78 km2.
5.2. Tested Rainfall Events and Estimated Hydrograph
A large precipitation event with 130 mm was observed on the 23th of August (2013) and 136.5 mm on the 16th of August (2014). These events were applied for verification of our method.
The estimation of the hydrograph was performed as follows: first, the unit of flood discharge was calculated as aforementioned. Second, the block discharge was calculated, also the same way as in previous research. Third, the discharge over the whole tested area was estimated by flood rooting according to the
channel network at the study area.
5.3. Comparison of Observed Discharge with Estimated Discharge
The comparison was performed as show in Figure 12. Both figures, the observed
and estimated, coincided well in the two cases. From the above facts, our proposed method was validated with the observation data.
6. Conclusions
This research describes the flood discharge associated with urbanization from the perspective of the unit flood discharge, which is defined as the flow rate on end-member drainage canals. The study area was selected as the Kurabe River basin, which is 17.5 km2 in area and is located in the Hokuriku Region of Japan, which has a very steep landscape. Twenty-one rainfall events of 10-minute intervals were selected, and five urbanized years of 1976, 1987, 1991, 1997 and 2009 were tested.
From 1976-2009 (35 years), the flood discharge increased by approximately 2.03 times because of the increase in residential area from 23% to 48%, and the maximum specific discharge in a block of 21.7 m3・s−1・km2 was a remarkably large amount.
Furthermore, changes in the hydrograph associated with urbanization, the effect of a small reservoir aimed for peak discharge control of the outflow from residential areas and the effect of unit flood discharge on the basin discharge were discussed. In particular, the effect of a small reservoir on the flood control was found to be very remarkable.
Furthermore, a comparison of flood discharge estimation by the “Rational formula” and our method was performed, which indicated that the former was smaller than the latter because the former had a lag time that was too long, as es-
![]()
Figure 12. Comparison of the observed and estimated hydrographs.
timated by the “Rziha formula” compared to the latter. Moreover, the lag times
obtained using both Kadoya and Fukushima’s experimental formula and our method were compared. This result indicated that the reasonability of coefficient C was 90 - 100 in Kadoya and Fukushima. Finally, the validity of the procedure was confirmed at the study area in the basin by comparing the observed and estimated flood discharge.
Based on the above results, we concluded that our proposed method is very effective for evaluating the flood discharge changes from urbanization. In particular, to estimate the flood damage from a very heavy rainfall event locally, which has begun to occur recently very frequently, requires local and individual runoff analysis rather than general.
In the future, the effectiveness of small reservoirs for flood control should be researched further because the effect depends on the characteristics of related rainfall events and facilities. This problem is very important to solve in the future because urbanization will progress but the capacity of the main canal is limited.
Acknowledgements
We express great thanks to Shitika irrigation district for valuable data on drainage and for Dr. Keiji Takase, Dr. Hiroshi Takimoto and Youichi Fujihara, professor and associate professor of Ishikawa Prefectual University for valuable comments.