Numerical Modeling of Sedimentary Basin Formation at the Termination of Lateral Faults in a Tectonic Region where Fault Propagation has Occurred

1. Introduction

1.1. Pull-apart basin forming at the termination of lateral faults

When a sedimentary basin forms at the termination of a lateral faults, it is known as a pull-apart basin. It is well known that tectonic basins, such as pull-apart basins, are generally formed at the termination of right-lateral, right-stepping and left-lateral, left-stepping fault systems (e.g., [1]). This is mainly caused by the formation of subsidence at the fault termination by the lateral motion of the faults. Subsidence is therefore likely to be found piled up at the termination of right lateral right-stepping and left lateral left-stepping fault systems. In contrast, uplift structures are formed at the termination of right-lateral left-stepping and left-lateral right-stepping fault systems, because the terminations are located in an area where uplift is piled up, due to the lateral motion of the fault (Figure 1). Such structures are found in many places globally, and their fundamental formation mechanisms have been numerically simulated by numerous researchers (e.g., [2, 3]).

Figure 1.

Schematic illustration of sedimentary basin and uplift at the terminations of lateral fault system. In this figure, right lateral motion was assumed to each vertical fault (pink line). Vertical displacements are normalized by the maximum value of absolute values of the total vertical displacement field, and they do not have a unit.

Katzman et al. [3] attempted to restore the subsurface structures of the Dead Sea, estimated from gravity anomalies (e.g., [4]), by means of Boundary Element Modeling (BEM), and indicated that it is necessary to assume long overlapping faults and a very high Poisson's ratio in order to restore the Dead Sea pull-apart basin. Rodgers [2] attempted to simulate the formation of a pull-apart basin by means of dislocation modeling (e.g., [5]), and this was probably the first study which discussed the formation of a pull-apart basin using numerical modeling.

1.2. Dislocation modeling

In general, dislocation modeling is used for the quantitative interpretation of crustal deformation caused by earthquakes and/or volcanic activity (e.g., [6, 7]). Surface or interior displacements or strains can be calculated by considering dislocations on a plane embedded in an elastic isotropic half-space (e.g., [8, 9]).

Studies on dislocation theory and its applications began with Steketee [10] and were refined by Okada [9], who derived closed analytical solutions for surface and internal deformations or strains due to shear and tensile faults, with arbitrary dip angles in a half-space. During this period, many researchers applied dislocation theory to inhomogeneous media and considered the effects of viscoelasticity and poroelasticity (e.g., [11-18]). Such basic theories and their applications to earth science have been described in many textbooks (e.g., [19-21]). In addition, the fundamental idea of dislocation theory has been applied to theory for the interpreting the gravity and potential changes due to large earthquakes (e.g., [22-24]). In particular, Okubo's formula [23] was applied to an interpretation of gravity changes due to the 2011 Tohoku-Oki earthquake, as observed by the Gravity Recovery And Climate Experiment (GRACE) [25], and Wang et al. [26] were successful in discussing co- and post-seismic deformation independently of GPS data.

As mentioned above, many useful solutions have been derived, although all solutions cannot be referred to as being “closed analytical solutions.” Closed analytical solutions (derived by [8, 9]) for a dislocation plane embedded in an elastic isotropic half-space are often employed for the quantitative interpretation of crustal movements, because they are very simple and are successful in explaining observational data. Consequently, some useful simulation software have been developed using Okada’s formula [9], (for example Coulomb (e.g., [27, 28])) and/or stress evaluation methods (for example ΔCFF (e.g., [29])), and the method has been applied to the interpretation of crustal movements and/or the evaluation of stress changes due to earthquakes (e.g., [30-32]).

In this study, we employ the analytical solutions of Okada [8] for the numerical simulation of sedimentary basin formation.

For a typical example of a pull-apart basin formed by dislocation, it would be appropriate to use a restoration model of Beppu Bay and Figure 1 for a numerical simulation. Beppu Bay is located at the junction between the Honshu and Ryukyu arcs, and it is formed at the western end of the Median Tectonic Line (MTL), which is the largest right-lateral tectonic line in southwestern Japan (c.f., [33]). Another right-lateral fault, the Kurume-Hiji Line (KHL: [34]), is located 30 km north of the MTL. The MTL and the KHL are arranged as a right-lateral right-stepping fault system, and the potential for producing a pull-apart basin exists. Kusumoto et al. [35] therefore applied the simple dislocation plane concept of Okada [8] to this fault system and showed that Beppu Bay was formed as a pull-apart basin by the right-lateral faulting of the MTL and KHL. In addition, they suggested that two major tectonic events (namely, the formation of half-graben caused by a north–south extension, and the formation of the pull-apart basin caused by east–west compression) occurred in this region.

In Figure 1 and in [35], a single motion of the fault was assumed on the dislocation planes. However, it is well known that active faults undergo multiple movements over geological time scales. In order to reflect this effect in the simulation, as the Okada's solution [8] is based on linear elasticity, we attempted to introduce the history of fault activity into the numerical model by superimposing analytical solutions for different fault parameters on a single fault. For example, if the geological evidence suggests that the active zone of a lateral fault shifted along its strike direction over time, we express its tectonic history by superimposing analytical solutions for different fault lengths (Figure 2).

Figure 2.

Illustration reflecting fault propagations in the dislocation modeling. Analytical solutions for different fault parameters are superimposed onto a single fault plane. Vertical displacements are normalized by the maximum value of absolute values of the total vertical displacement field, and they do not have a unit.

Itoh et al. [36] introduced this simulation technique while attempting to interpret topography and paleomagnetic data for the Takayama basin, central Japan. They showed that the basic structure of the Takayama basin, and changes in the declination of the thermoremanent magnetization, can be restored by considering the cumulative activity of two right-lateral faults (the Enako Fault and the Makigahora Fault) and a reverse fault (the Harayama Fault). In addition, they showed that the reverse fault can be divided into two segments that moved independently with a time lag between their active periods. Because this modeling technique can take into account the type and amount of fault motion, we can reflect the history of fault activity based on certain geological evidence, including paleomagnetic studies, in the numerical modeling, and can discuss the tectonics in detail.

This modeling technique has also been applied to the formation of pull-apart basins located in Hokkaido by Itoh et al. [37] and Tamaki et al. [38]. Tamaki et al. [38] found that a strike-slip fault motion reaching 30 km, is required to restore the distribution and volume of the Minami-Naganuma Basin located in southern central Hokkaido.

Here, we describe the disadvantages of the dislocation modeling defined using elasticity, and the solutions for the dislocation plane are given as a range of the linear elasticity. Since fracture, flow and other non-linear phenomena seen in the general solid material are not considered in the model, it is difficult to directly compare the amount of displacement between the modeled structure and the actual structure in long-time scale modeling. In general, dislocation modeling (including visco-elasticity effects) is often employed in discussions on crustal movements over a long time-scale such as on a geological time scale (e.g., [18, 39, 40]). In addition, Finite Element Modeling (FEM), Finite Difference Modeling (FDM) and Discrete Element Modeling (DEM) can simulate the formation processes of sedimentary basins or the building processes of mountains over a geological time-scale and can quantitatively discuss the mechanisms of their formation on a time axis (e.g., [41-43]).

As mentioned above, dislocation modeling defined using a range of the linear elasticity has disadvantages in the ability to directly compare the amount of displacement between the modeled structures and the actual structures. However, when we simply discuss the essential aspects of tectonics from the distribution pattern of structures caused by fault motions, dislocation modeling is a very useful tool because it provides the pattern of displacement using easy calculations.

1.3. Aims of this study

The aims of this study are to simplify the complex formation processes of sedimentary basins through numerical simulations and to show that the simplification enables us to estimate deductively which processes cause materialization. Central Hokkaido was selected as the field in which to achieve these aims, as many sedimentary basins are distributed in this area.

As will be described later in the paper, many sedimentary basins were formed from 48 to 12 Ma in central Hokkaido, and it is difficult to discuss their formation processes using only observational data because of their complex distribution both at and below the surface. In order to simplify the complex formation processes of sedimentary basins, we attempted to restore sedimentary basins using the advanced technique already mentioned, and we evaluated the fault type and the amount of movement required to form these sedimentary basins.

In the following sections, we describe the basic background and gravity anomaly in central Hokkaido, and we attempt the restoration of the sedimentary basins.

2. Basic background of central Hokkaido

2.1. Geophysical background

Hokkaido is located on the North American plate, at a junction of the Northeast Japan arc and the Kurile arc (Figure 3). Using recent GPS observations, an east-west compressive strain field has been observed in the northern part of Hokkaido, and this strain field is considered to be caused by the convergence of the Eurasia plate with the northern part of Hokkaido functioning as a part of the plate boundary (e.g., [44]).

Figure 3.

Location map of our study area. Hokkaido is located on the North American plate, at a junction of the Northeastern Japan arc and the Kurile arc. Gray dashed line indicates the old plate boundary between the Eurasian and North American Plates. Rectangular area by gray thin line indicates the study area of Itoh and Tsuru [58].

Although the present plate boundary between the North American plate and the Eurasian plate exists in the Sea of Japan, it is known that the plate boundary was located in central Hokkaido at around 13 Ma. This period of time corresponds to the stage when the uplifting of the Hidaka Mountains began (e.g., [45]). This tectonic framework is controlled by the dextral oblique collision between the Eurasian and North American Plates and the oblique subduction of the Pacific Plate beneath the Kurile Trench. It is considered that the Kurile arc migrated into the southwestward as a forearc sliver by the oblique subduction of the Pacific Plate and that the Hidaka Mountains would be formed by collision of the Kurile arc and the Northeast Japan arc (e.g., [46-48]).

Since, consequently, it is an important area for understanding characteristics and mechanism of collision zone, numerous geophysical surveys (e.g., seismic prospecting, gravity surveys, electromagnetic surveys) have been carried out around the Hidaka Mountains in order to obtain information regarding the subsurface structures and to apply such knowledge to a tectonic discussion of the Mountains and Hokkaido (e.g., [49-54]). Using these surveys, subsurface structures which indicate a collision between the Northeast Japan arc and the Kurile arc have been obtained, and have contributed to tectonic discussions. However, in contrast, geophysical studies in the sedimentary basin area in central Hokkaido are limited in number.

2.2. Geological and tectonic background

The geological characteristics in Hokkaido are that Cenozoic strata consist of island-arc-trench systems of the Northeast Japan arc and Kurile arc, and that each Cenozoic strata distributed in the Northeast Japan arc and the Kurile arc appear in the western half and eastern half area of Hokkado, respectively. This characteristic also appears in Neogene and Quaternary strata, volcanoes and their products, and topography (e.g., [41]).

Figure 4 shows the distribution of the Paleogene strata in a north-south direction in central Hokkaido. This distribution traces the old plate boundary. This N-S elongation area is included in the Ishikari-Teshio Belt that is underlain by the Cretaceous Yezo Group, and is regarded as a typical sequence in a forearc basin setting [55]. It is known that the Paleogene sedimentary strata were deposited during the early Eocene and Oligocene in almost the entire region (e.g., [56]).

Figure 4.

Distribution of the Paleogene strata. Green and blue areas indicate distribution areas of the Paleogene sedimentary layer under and on the surface, respectively. (After Kurita and Hoyanagi [56])

In the study area, sedimentary basins and sedimentary layers were formed during 48–12 Ma and have been complexly distributed. They have been divided into 5 stages according to their formation: The Ishikari stage (48–40 Ma), The Horonai stage (40–32 Ma), The Minami Naganuma stage (34–20 Ma), and The Kawabata stage (15–12 Ma). The Ishikari stage is divided into early (48–45 Ma) and late (45–40 Ma) stages. The shape of the sedimentary basin and the distribution of sedimentary layers are shown in Figure 5.

It is known that the Ishikari Group differentially subsided and was then divided into several components [57]. Based on detailed sedimentological studies, Takano and Waseda [57] also points out that the rate of subsidence accelerated during deposition of the Ishikari Group.

The Ishikari stage is the sedimentation stage of the Ishikari Group, corresponding to the Eocene and is divided into early and later stages according to the sedimentation style. Sediments in the early Ishikari stage are distributed shallowly and widely (Figure 5A). In this stage, the sedimentary basins “A”, “B” and “C” were formed (Figure 5A), and from well data their depths are estimated to be 600 m, 500 m and 1000 m, respectively. Sediments in the later Ishikari stage are distributed deeply and narrowly (Figure 5B). In this stage, the sedimentary basins named “A” and “C” were formed (Figure 5B), and from well data their depths are estimated to be 2800 m and 400 m. Itoh and Tsuru [58] identified a NNW-SSE trending deformation zone bounded by large transcurrent faults including T1 and T2 later describing from seismic reflecting data in the northern part of the Northeast Japan forearc (Figure 3) and their right lateral motions have been indicated by the clockwise rotation of Paleogene marine sediments and by paleogeographic reconstruction. Since, as already mentioned, the present study area (the western half of Hokkaido) has same Cenozoic strata distributed in the Northeast Japan arc, study area of Itoh and Tsuru [58] and our study area are geologically continuous in the Paleogene time. Consequently, it is expected that a right lateral motion of the crust was dominant.

The Horonai stage is the sedimentation stage of the Horonai Formation, the Tappu Group, the Sankebetsu Formation and the lower Magaribuchi Formation. This stage corresponds to the Eocene and the early Oligocene. In this stage, sedimentary basins “A”, “B”, “C”, “D”, “E” and “F” were formed (Figure 5C), and from well data their depths are estimated to be 3500 m, 1200 m, 1200 m, 600 m, 300 m, and 1500 m, respectively. It is expected that a right lateral motion was dominant in this stage, because right lateral motion was also dominant in the Eocene and the late Oligocene (see below).

The Minami-Naganuma stage is the sedimentation stage of the upper Magaribuchi Formation, the Minami-Naganuma Formation, the Horomui Formation, and the upper Sankebetsu Formation. This stage corresponds to the late Oligocene and the early Miocene. In this stage, sedimentary basins “B”, “D”, ”E” and “F” were formed (Figure 5D), and from well data their depths (B, E and F) are estimated to be 2000 m, 300 m and 1500 m, respectively. The maximum depth of basin “D” is unknown because of lack of the well data and/or of outcrop section of the whole Minami-Naganuma Fromation. Itoh et al. [37], Tamaki et al. [38] and Itoh and Tsuru [59] pointed out that all basins formed in the Ishikari-Teshio Belt in this stage are pull-apart basins. Tamaki et al. [38] showed that using dislocation modeling, a 30 km right-lateral strike-slip is required to restore the actual distribution and volume of the basin. Kurita and Yokoi [60] also stated that lateral faulting was dominant in forming some of the tectonic structures during the late Oligocene.

The Kawabata stage is the sedimentation stage of the Kawabata Formation, the Ukekoi Formation, the Fureoi Formation, the Kotanbetsu Formation and the Masuporo Formation. During the Neogene, Japan was affected by the opening event of the back-arc basin of the Sea of Japan. In this stage, sedimentary basins “A”, “B”, “D”, “E”, “F1” and “F2” were formed (Figure 5E), and from well data their depths are estimated to be 2000 m, 4000 m, 4000 m, 3500 m, 2000 m and 2000 m, respectively. As mentioned above, a lateral motion of the crust was dominant during the early Neogene [37, 38, 59]. Although a building of the Hidaka Mountains in around 13 Ma has been pointed out (e.g., [45]), details are unknown.

Figure 5.

Shapes of sedimentary basins in (A) early Ishikari stage (48-45 Ma), (B) late Ishikari stage (45-40 Ma), (C) Horonai stage (40–32 Ma), (D) Minami-Naganuma stage (34–20 Ma) and (E) Kawabata stage (15–12 Ma). Isopach maps of the Horonai stage and the Kawabata stage are after Association of Natural Gas Mining and Association for Offshore Petroleum Exploration [75].

3. Bouguer gravity anomaly

Numerous geological and geophysical surveys have been carried out in the Hokkaido area, and each survey has played an important role in the understanding of crustal characteristics and tectonic events in the area. In particular, seismic prospecting has proved very useful in obtaining information relating to subsurface structures. However, seismic prospecting is almost two-dimensional, and it is difficult to intuitively understand the subsurface structures as three dimensional structures, even when provided with data from more than one profile. In contrast, the characteristics of gravity anomaly maps are easy to interpret and can be used to roughly estimate three dimensional subsurface structures from the data. Figure 6 shows the Bouguer gravity anomaly map of the study area. This map is based on the gravity mesh data by Komazawa [61]. The Bouguer density of 2670 kg/m³ was employed.

Figure 6.

Bouguer gravity anomaly map. This map is based on the gravity mesh data by Komazawa [61], and the Bouguer density of 2670 kg/m³ was assumed. Contour interval is 10 mGal.

There are negative gravity anomalies in the southern and northern parts of central Hokkaido. The negative gravity anomaly in the northern part reaches -20 mGal (Figure 6 and Area I in Figure 7) and the southern negative gravity anomaly is less than -100 mGal (Figure 6 and Area III in Figure 7). The southern negative gravity anomaly located at the curved subduction zone is the lowest in the country. From seismic prospecting, it is known that this negative gravity anomaly consists of a very thick sedimentary layer (5–8 km), with a velocity of 2.5–4.8 km/s (e.g., [50]). The sedimentary layer was formed by imbrications associated with the collision process of the Northeast Japan arc and the Kurile arc (e.g., [46-48]). In contrast, there is a positive gravity anomaly in the area of the mountains, and the mountain elevations are roughly less than 2000 m. It would not be necessary to consider isostasy for the mountains, because the mountain elevations are not very high and the gravity anomaly in this area is positive.

Figure 7.

Bouguer gravity anomaly map. Gray indicates gravity low area less than 20mGal. A-A', B-B', C-C' and D-D' show gravity anomalies along each profile of four red lines shown in the Bouguer gravity anomaly map.

A flat gravity anomaly of less than 20 mGal is distributed like a belt between the northern and southern gravity anomalies (Figure 6 and Area II in Figure 7.). Figure 7 is the gravity anomaly map that the area less than 20 mGal was painted by gray. This painted area corresponds to the area where the Paleogene strata distribute under the surface (Figure 4). We show four cross section profiles (three E-W profiles, A-C, and one N-S profile, D) of the Bouguer gravity anomaly in Figure 7. From these profiles, the gravity anomalies in the region are shown to have the characteristics as follows:

1. Gravity anomalies in the west-east direction have a regional trend which tilts toward the east. This could indicate the regional gravity field in Hokkaido.

2. Gravity anomalies less than 20 mGal have a steep gradient on the east side, while those on the west side vary gently. These patterns of gravity anomalies indicate a depression structure called a “half-graben”. Since there are many confirmed lateral faults and reverse faults in this region and no normal faults, it is considered that these patterns of gravity anomalies are caused by structures formed by the activities of lateral faults and/or reverse fault.

3. Gravity anomalies in the north-south direction are relatively high and flat at the center of Hokkaido. It is possible that the high density of metamorphic belts near this region affect the observed gravity anomalies. Another cause to be considered could be the effect of subsurface structures such as a reduction of low density materials (e.g., a thin sedimentary layer) or an increase of high density material (e.g., uplift of the mantle).

In general, gravity anomalies are caused by spatial variations of subsurface structures, and indicate a deficiency or an excess of mass under the surface. In general, high gravity indicates the existence of a mass excess or of high density materials, and low gravity indicates the existence of a mass deficiency or of low density materials. These deficiencies or excesses, of mass can be evaluated quantitatively using Gauss's theorem (e.g., [62]).

$$\Delta M=\frac{1}{2\pi G}{\displaystyle \iint \Delta g\left(x,y\right)dxdy}$$E1

Here, g(x, y) is the gravity anomaly data given on xy mesh with a constant interval. G, π and ΔM are the universal gravitational constant, circular constant and deficiency or excess of mass, respectively. Equation (1) is described by an infinite integration and it is difficult to perform an infinite integration with actual field data. Consequently, we understand this as being an approximate calculation and perform a numerical integration within a finite area (S) as follows:

$$\Delta M=\frac{1}{2\pi G}{\displaystyle {\iint }_S\Delta g\left(x,y\right)dxdy}$$E2

We applied equation (2) to three areas, I, II and III, and we attempted to estimate the magnitude of mass deficiency for the formation of a gravity anomaly less than 20 mGal in each area. In the calculations, we employed the Gauss-Legendre numerical integral formula (e.g., [63]).

As a result, mass deficiencies of 4.7×10³ Gton, 8.6×10² Gton, and 1.5×10⁴ Gton were estimated in areas I, II, and III, respectively. There are large mass deficiencies in areas I and III, where the negative gravity anomalies observed are very large and a small mass deficiency in area II. In central Hokkaido, the amount of mass deficiency is different by about two digits in both the maximum and the minimum values.

The amount of mass deficiency can be transformed into the volume (V) of sediment by the following equation, under a condition assuming an appropriate density contrast (Δρ).

$$V=\frac{\Delta M}{\Delta \rho }$$E3

As an example, when a density contrast of 300 kg/m³ is assumed, volumes of sediment of 1.6×10⁴ km³, 2.9×10³ km³ and 5×10⁴ km³ are estimated in areas I, II and III, respectively.

As mentioned above, gravity anomaly indicates also spatial variations of subsurface structures including the location of tectonic lines and/or faults. It is well known that if there is a tectonic line or a fault with a large gap in the vertical direction, the spatial distribution of the gravity anomaly varies steeply around these structures. The variation rate of the spatial distribution of the gravity anomaly is called the “horizontal gradient of gravity anomaly”, and it is given by the first derivative (e.g., [64, 65]) or the second derivative (e.g., [4, 66]). In general, the first derivative of the gravity anomaly is more practical, because the calculation used is very simple and the geophysical and geological interpretations for the calculated results are straightforward.

We employed the first derivative of the gravity anomaly defined by the following equation (4), and calculated the horizontal gradient of the gravity anomaly (Figure 8):

$$\sqrt{{\left[\frac{\partial g\left(x,y\right)}{\partial x}\right]}^2+{\left[\frac{\partial g\left(x,y\right)}{\partial y}\right]}^2}$$E4

Figure 8 shows the distribution of the horizontal gradient of the Bouguer gravity anomaly more than 2 mGal/km. The contour interval is 1 mGal/km. Although there are no continuous horizontal gradient anomalies within the area where the gravity anomaly is less than 20 mGal, the continuous horizontal gradient anomalies appear around this area. This may indicate that there are not tectonic lines including faults having large vertical deformation within this gravity low area less than 20 mGal and/or that gravity anomalies due to these tectonic lines are hidden by thick sediments, although faults with large vertical deformations actually exist.

Figure 8.

Distribution of the horizontal gradient of the Bouguer gravity anomalies more than 2 mGal/km. Contour interval is 1 mGal/km.

4. Restoration of sedimentary basins

In general dislocation modeling, the dislocation plane is assumed in the modeled crust by referring to the distribution of existing active faults and/or tectonic lines, and the surface deformations are calculated by assigning displacements on the plane. If the area for modeling is small, or if the tectonics and faults assumed for modeling are clear, such a modeling procedure is useful and practical (e.g., [35, 36, 38, 67]).

However, when details of the tectonics and/or the moved faults are not so clear (as in our study), the faults and their displacements for modeling are assumed experientially from characteristic distributions of target structures, by referring to more regional rough tectonics and fault distributions. The faults and their displacements (appropriately assumed) can then be considered as an initial model and can then be corrected by trial and error, so that the calculated results fit to the actual structures or their distribution pattern.

There are numerous small faults in central Hokkaido. As mentioned above, the details of tectonics and faulting in this area are unclear. It would, therefore, be impossible to attempt to model each fault for restoring the distribution of the sedimentary basins by trial and error. Consequently, in this study, we assumed that the dislocation plane used for the modeling was not a fault plane, but a typical or average plane of a fault zone. After trial and error, we defined the nine fault zones as shown in Table 1 and Figure 9, and employed them for numerical simulations. Each fault included in these fault zones is listed in Table 1 (with literature).

Figure 9.

Fault zones defined in this study. Each fault included in these fault zones is listed in Table 1 with literature.

No.	Fault zone	Preference (recognition as tectonic zone)	Specific faults [reference]
1	Horonobe	Research Group for Active Faults [73]	Horonobe Fault [76]; Higashino Fault [77, 78]
2	Tenpoku	Ikeda et al. [74]; Itoh et al. [37]	Enbetsugawa Tectonic Line [79]
3	Chikubetsu	this study	Shosanbetsu Fault [78]; Chikubetsu Anticlinal Fault [78]
4	Onishika	this study	Onishika Fault [80]; Tenguyama Thrust Fault [81]; Okinaiyama Fault [82]; Shimokinebetsuzawagawa Fault [82]; Shisenzawa Fault [82]; Poroshiriyama-Higashi Fault [82]; Horoshin Fault [82]
5	Rumoi-Shintotsu	Association of Natural Gas Mining and; Association for Offshore Petroleum; Exploration [75]	Shozuzawa Fault [83]; Horonuka Fault [83]; Chashinai Fault [84]; Sorachi Fault [85]
6	T1	Itoh and Tsuru [58, 59]	N/A (Most of the fault trace is in offshore area and buried with sediments.)
7	T2	Itoh and Tsuru [58, 59]	N/A (Most of the fault trace is in offshore area and buried with sediments.)
8	Hidaka-North ^(a)	this study	See footnote.
9	Hidaka-South ^(b)	this study	See footnote.

Table 1.

Fault zones. Definition of each fault zone, and relationship fault zone and specific faults. (a) Hidaka-North Fault Zone is collectively defined as the eastern margin of N-S serpentinite zone along the longitudinal mountainous range in northern Hokkaido. (b) Hidaka-South Fault Zone is assigned to the western margin of the area of the highest recent uplift rate in Hokkaido (i.e. Hidaka Mountains), of which tectonic and structural context is still controversial.

In the following subsections, we give the results of numerical simulations, accompanied by simple explanations. The faults moved during each stage and their fault parameters are listed in Table 2. In calculations, the total movement of each fault plane was expressed by fault motions of 1000 times, and a high Poisson’s ratio of 0.4 was assumed because of its application to modeling in the geological time scale (e.g., [3, 36-38, 68]).

No.	Fault zone (fault type)	Early Ishikari	Late Ishikari	Horonai	Minami Naganuma	Early Kawabata	Late Kawabata
1	Horonobe (Right lateral motion)	-	-	(88, 500)	(94, 500)	(36, 61)	-
2	Tenpoku (Right lateral motion)	-	-	(28, 24)	(14, 500)	(41, 32)	-
3	Chikubetsu (Right lateral motion)	-	-	-	-	(21, 36)	-
4	Onishika (Right lateral motion)	-	(48, 150)	(61, 63)	-	-	-
5	Rumoi-Shintotsu (Right lateral motion)	(22, 47)	-	(61, 94)	(48, 78)	(44, 73)	-
6	T1 (Right lateral motion)	(14, 500)	(14, 500)	(14, 500)	-	-	-
7	T2 (Right lateral motion)	(25, 500)	(25, 500)	(61, 500)	(44, 500)	(92, 500)	-
8	Hidaka-North (Right lateral motion)	(36, 500)	(36, 500)	(29, 500)	-	-	-
	Hidaka-North (Reverse motion)	-	-	-	-	-	(2.5, 500)
9	Hidaka-South (Reverse motion)	-	-	-	-	-	(2.5, 212)

Table 2.

Fault zones moved in each stage. Values in (A, B) shown in Table indicate total slip (A: km) amount assumed on the fault plane and fault initial length (B: km). Width, depth and dip angle of fault moved as right lateral fault are assumed to be 15km, 15km and π/2, respectively. Width, depth and dip angle of fault moved as reverse fault are assumed to be 17.32km, 15km and π/3, respectively.

4.1. Early Ishikari stage (48–45 Ma)

We attempted to restore the sedimentary basins formed in the early Ishikari stage, named “A”, “B” and “C” (Figure 5A). Results are shown in Figure 10A. Right lateral movements of fault zones were required, (including the T1 fault, T2 fault, the Rumoi-ShinTotsu tectonic line and the Hidaka-North fault), in order to restore these three sedimentary basins. The amount of movement of each fault was determined by trial and errors, and results are shown in Table 2 and are as follows:

1. T1 fault zone: 14 km

2. T2 fault zone: 25 km

3. Rumoi-ShinTotsu tectonic zone: 22 km

4. Hidaka-North fault zone: 36 km

From these amounts of displacement on each fault plane, it is estimated that during this tectonic stage, a horizontal movement reaching around 100 km occurred.

Figure 10.

Distribution pattern of sedimentary basins in each stage restored by the dislocation modeling. A: early Ishikari stage. B: late Ishikari stage (45-40 Ma). C: Horonai stage (40–32 Ma). D: Minami-Naganuma stage (34–20 Ma). E: early Kawabata stage (15–13 Ma). F: late Kawabata stage (13-12Ma). Vertical displacement amounts are given in m.

From geological observations it is known that sedimentary basin “C” is the biggest basin and reaches to a depth of 1000 m. In our modeling, the depth of the modeled basin reached 1400 m. Sedimentary basins, “A” and “B” reach about 600 m and 500 m depths, respectively, and both depths of restored basins (modeled basins) corresponding to these basins are 1000m as a result of dislocation modeling. Differences between the actual basin depth and the modeled basin depth, namely “the modeled basin depth–actual basin depth,” were +400 m, +400 m and +500 m for basins “C”, “A” and “B”, respectively. The amount of restored subsidence may be a little large.

Here, we assigned the right lateral motion to each fault as mentioned above, in order to restore the spatial patterns of basin distribution. If the amount of lateral motion of each fault is reduced to adjust to the depth component of each basin, it is not possible to restore the spatial distribution patterns of the basins.

4.2. Late Ishikari stage (45–40 Ma)

In the late Ishikari stage, sedimentary basins “A” and “C” (Figure 5B) were restored (Figure 10B). The right lateral movements of fault zones (including the T1 fault, T2 fault, Onishika fault and the Hidaka-North fault) were required in order to restore these two sedimentary basins. The amount of movement of each fault was determined by trial and error (see Table 2) and is as follows.

1. T1 fault zone: 14 km

2. T2 fault zone: 25 km

3. Onishika fault zone: 48 km

4. Hidaka-North fault zone: 36 km

From these amounts of displacement on each fault plane, it is estimated that during this tectonic stage, a horizontal movement reaching about 123 km occurred.

From geological observations, sedimentary basin “C” is known to be the largest basin and reaches a depth of 2800 m. The depth of the model basin in our modeling reached 1800 m. Sedimentary basin “A” reaches a depth of about 400 m, and the modeled basin corresponding to this had a depth of 500 m. The differences between the actual basin depth and the modeled basin depth were -1000 m and +100 m in basins “C” and “A”, respectively. The restored subsidence amount of basin “C” was smaller than the actual basin depth. If the lateral motions of the Hidaka-North fault zone and Onishika fault zone were increased to adjust to the depth component of basin “C”, it was not possible to restore the spatial distribution patterns of the basins.

4.3. Horonai stage (40 Ma–32 Ma)

In the Horonai stage, we attempted to restore six sedimentary basins, “A”, “B”, “C”, “D”, “E” and “F” (Figure 5C). Results are shown in Figure 10C. The right lateral movements of fault zones (including the T1 fault, T2 fault, the Rumoi-ShinTotsu tectonic line, Onishika fault, the Tenpoku fault, Horonobe fault and Hidaka-North fault) were required in order to restore these six sedimentary basins. The amount of movement of each fault was determined by trial and error and is shown in Table 2 and as follows:

1. T1 fault zone: 14 km

2. T2 fault zone: 61 km

3. Rumoi-ShinTotsu tectonic line zone: 61 km

4. Onishika fault zone: 61 km

5. Tenpoku fault zone: 28 km

6. Horonobe fault zone: 88 km

7. Hidaka-North fault zone: 29 km

From the amount of displacement of each fault plane, it is estimated that a horizontal movement reaching about 342 km occurred during this tectonic stage.

From geological observations, it is known that depths of the sedimentary basins “A”, “B”, “C”, “D”, “E” and “F” reach to 3500 m, 1200 m, 1200 m, 600 m, 300 m and 1500 m, respectively. In our model, depths of modeled sedimentary basins “A”, “B,” “C”, “D”, “E” and “F” reached 1000 m, 600 m, 1800 m, 1100 m, 900 m and 1000 m, respectively. The differences between the actual basin depth and the modeled basin depth were -2500 m, -600 m, +600 m, +500 m, +600 m, and -500 m in basin “A”, “B”, “C”, “D”, “E” and “F”, respectively.

4.4. Minami-Naganuma stage (34–20 Ma)

We attempted to restore sedimentary basins “B”, “D”, “E” and “F” in the Minami-Naganuma stage (Figure 5D), and the results are shown in Figure 10D. The right lateral movements of fault zones, (including the T2 fault, the Rumoi-ShinTotsu tectonic line, Tenpoku fault and the Horonobe fault), were required in order to restore these four sedimentary basins.

In this stage, Tamaki et al. [38] have restored already the Minami-Naganuma basin (corresponding to basin “B”: pull-apart basin) located in south central Hokkaido. We referred to their results and determined the amount of movement of each fault by trial and error. The amount of fault movement is shown in Table 2 and as follows:

1. T2 fault zone: 44 km

2. Rumoi-ShinTotsu tectonic zone: 48 km

3. Tenpoku fault zone: 14 km

4. Horonobe fault zone: 94 km

From the amounts of displacement on each fault plane, it is estimated that horizontal movement reaching about 200 km occurred during this tectonic stage.

From geological observations it is known that the depth of sedimentary basins “B”, “E” and “F” reached 2000 m, 300 m and 1500 m, respectively. As already mentioned, the maximum depth of basin “D” is unknown. In our model, the depths of the modeled sedimentary basins “B”, “D”, “E” and “F” reached 1200 m, 1100 m, 700 m and 1200 m, respectively. The differences between the actual basin depth and the modeled basin depth were -800 m, +400 m, and -300 m for basins “B”, “E” and “F”, respectively.

4.5. Kawabata stage (15-12 Ma)

In the Kawabata stage, we attempted to restore six sedimentary basins, “A”, “B”, “D”, “E”, “F1” and “F2” (Figure 5E), and the result is shown in Figure 10E. The right lateral movements of fault zones, (including the T2 fault, Rumoi-ShinTotsu tectonic line, Onishika-chikubetsu fault, Tenpoku fault and the Horonobe fault), were required in order to restore these six sedimentary basins. The amount of movement of each fault is determined by trial and error and is shown in Table 2 and as follows:

1. T2 fault zone: 98 km

2. Rumoi-ShinTotsu tectonic line zone: 44 km

3. Onishika-chikubetsu fault zone: 21 km

4. Tenpoku fault zone: 41 km

5. Horonobe fault zone: 36 km

From the amount of displacement on each fault plane, it is estimated that a horizontal movement reaching about 240 km occurred during this tectonic stage.

As described above, these lateral motions could restore the basic distribution pattern of six sedimentary basins formed in the Kawabata stage. However, the whole distribution pattern of sedimentary basins in this stage could not be restored. After repeated trial and error, it was found that the reverse motion of two fault zones, (including the Hidaka-north fault and Hidaka-south fault), is necessary to restore the whole distribution pattern of the sedimentary basins in this stage (Figure 10F). In fact, such a reverse motion is required to successfully restore the whole basin distribution pattern. The amount of reverse motion required is shown in Table 2 and as follows:

1. Hidaka-north fault zone: 2.5 km

2. Hidaka-south fault zone: 2.5 km

Information on the depth of the sedimentary basins “A”, “B”, “D”, “E”, “F1” and “F2” is obtained from geological observations, and depths are found to reach 2000 m, 4000 m, 4000 m, 3500 m, 2000 m and 2000 m. In our model, the depths of modeled sedimentary basins “A”, “B”, “D”, “E”, “F1” and “F2” reached 1100 m, 1800 m, 1600 m, 1900 m, 1300 m and 800 m. The differences between the actual basin depth and the modeled basin depth were -900 m, -2200 m, -2400 m, -1600 m, -700 m, and -1200m in basins “A”, “B”, “D”, “E”, “F1” and “F2”, respectively. Although the whole basin distribution pattern is good, the differences between the actual basin depth and the modeled basin depth are large in all the basins. If the lateral or reverse motion of each fault zone was increased to adjust the depth component of basins, it would not be possible to restore the spatial distribution patterns of the basins.

5. Discussion

As described in the previous sections within this paper, the distribution patterns of sedimentary basins restored by dislocation modeling are very similar to the actual distribution patterns of the sedimentary basins formed in each stage, although depth differences between the actual sedimentary basins and the restored sedimentary basins occurred because of the dislocation plane based on the linear elasticity. From these results, it is suggested that almost all the sedimentary basins in central Hokkaido can be explained as pull-apart basins, caused by right-lateral fault motions during the Paleogene. The results also show that sedimentary basins formed during the Kawabata stage were formed by a combination of right-lateral fault motions and reverse fault motions located at the western margin of the Hidaka Mountains.

Although we have simplified the distributions of each fault zone (Figure 9), it would be difficult for contiguous fault zones to move independently, simultaneously, as different faults, namely a reverse fault and a lateral fault, under their arrangement as shown in Figure 9. Consequently, we suggest that the Kawabata stage should be divided into two stages, namely the “early stage” and the “late stage”, from the viewpoint of the stress field or fault motion. By considering the continuity of tectonics or the stress field, it is found that the lateral movements were made in the early Kawabata stage and that the reverse movements were made in the late Kawabata stage.

From a geological viewpoint, it has been illustrated that the building of the Hidaka Mountains was caused by reverse fault motions (e.g., [45]), and the timing of this event has been considered as being during the late Miocene or around 13 Ma (e.g., [45, 69]). This geological view supports our results and ideas, and our results also support the tectonics constructed based on geological data. Hence, we suggest that the boundary of the late Kawabata stage is around 13 Ma.

In Figures, we show the total vertical displacement field calculated from the vertical displacement field in each stage (Figures 10). The vertical displacement fields shown in Figures 11 are normalized by the maximum value of absolute values of the total vertical displacement field restored by dislocation modeling. Thus, the displacement fields shown in Figures 11 do not have a unit.

Figure 11A illustrates the normalized displacement field map that the negative displacement areas are shown in gray. This shows the distribution of the subsurface sedimentary basins restored in this study, and the distribution is seen to be similar to the distribution of actual buried sedimentary basins formed during the Paleogene (Figure 4). Figure 11B is the normalized total vertical displacement pattern restored in this study. From this figure, it is found the deepest sedimentary basin restored is located at the center of central Hokkaido. In actual depth distribution, the basin “C” is not the deepest basin. However, this sedimentary basin has a depth reaching 6000 m and is large and deep basin.

Figure 11.

Total vertical displacement pattern. The displacement field was calculated by adding up the vertical displacements field restored in each stage by the dislocation modeling. A: Gray shows negative displacement areas, and indicates the distribution of the subsurface sedimentary basins restored in this study. B: Total vertical displacement pattern normalized by the maximum value of absolute values of the total vertical displacement field.

From results shown in Figures 10-11 and the discussion above, we conclude that the sedimentary basins formed from 48 to 12 Ma in central Hokkaido can be explained by the formation of pull-apart basins, due to right lateral motions (before 13 Ma), and by reverse motions of the Hidaka-North fault zone and the Hidaka-South fault zone after 13 Ma. Namely, although the right-lateral motions were predominant from 48 Ma to 13 Ma, the reverse motions were dominant in 13 Ma. This leads us to expect a significant change in the regional tectonic stress field during this stage. Although a change of the collision direction of the Northeast Japan arc and the Kurile arc could be cited as a possibility for this, a more accurate and quantitative future investigation is required.

We here reconsider the Bouguer gravity anomaly (Figure 6 and 7) and subsurface structures in viewpoint of tectonics mentioned above. As described in the section on the Bouguer gravity anomaly, gravity low area less than 20 mGal corresponds to the area where the Paleogene strata distribute under the surface. It is known that the southern negative gravity anomaly less than -100 mGal (area III in Figure 7) consists of a very thick sedimentary layer of 8 km, and the depth of the Moho discontinuity in this area is more than 30 km (e.g., [70, 71]), and that the negative gravity anomaly which reaches -20 mGal (area I in Figure 7) consists of a sedimentary layer of several kilometers thickness, and a Moho discontinuity of around 30 km in depth (e.g., [70, 71]). Consequently, we understand that the conspicuous gravity lows in areas I and III are caused by a thick sedimentary layer and a deep Moho. In contrast, the Bouguer gravity anomaly in area II is not negative, in spite of an area of subsidence of several kilometers depth. In actual fact, using receiver function analysis [71], it is reported that the depth of the Moho discontinuity in this area ranges from a depth between 26–31 km. The Moho in area II is about 6 km shallower than in areas I and III.

We made a simple subsurface structure model consisting of sedimentary layer, crust and mantle, based on information of subsurface structures mentioned above. We then calculated the gravity anomalies along a D-D' profile with assumption density contrasts of sedimentary layer-crust and crust-mantle being -300 kg/m³ and -500 kg/m³. We employed the two dimensional Talwani's method [72] in our calculations.

Figure 12 shows the simple subsurface structure model and the estimated gravity anomaly along the D-D' profile. The assumed subsurface structure model explains the Bouguer gravity anomalies very well, in spite of the model being very simple, and leads us to assume that the conspicuous subsidence at the surface would induce the crustal thinning-mantle uplift, via isostasy. Subsidence at the surface, caused by frequent lateral motions during the Paleogene, as shown and discussed in this paper, has reached to several km. In this study, right lateral motions of the crust reaching about 1000 km were required, in order to restore the sedimentary basins distributed in central Hokkaido, as shown in the previous section. Although the correct depths of sedimentary basins were not shown in our dislocation modeling because the model was based on linear elasticity, it is expected that deeper basins would be formed in the actual crust by lateral motions. Consequently, it is considered possible that conspicuous subsidence caused by a large lateral motion would induce the crustal thinning-mantle uplift. Or, since this area experienced a tension stress field over a period of about 35 million years in spite of the local stress caused by the pull-apart basin formation, this tension stress field might induce the crustal thinning-mantle uplift viscoelastically over a very long time scale. In either case, it is necessary to reconsider the estimated mass deficiency and volume of the sediment in area II from the Bouguer gravity anomalies, by correcting the effect of mantle uplift.

Figure 12.

Two-dimensional gravity modelling along D-D' profile in Figure 7. Blue circle shows measured Bouguer gravity anomalies and red solid line shows calculated values. Bottom figure shows the density structure. Density contrasts of sedimentary layer-crust and crust-mantle were assumed to be -300 kg/m³ and -500 kg/m³, respectively.

In this study, we suggested one tectonic model to explain the distribution of sedimentary basins located in central Hokkaido and formed between 48 Ma and 12 Ma. We then discussed the characteristics of the Bouguer gravity anomalies based on the tectonic model. In future studies, we will attempt to estimate the subsurface structure more accurately, and discuss the tectonics through simulation procedures considering a more realistic behavior of the crust and mantle (e.g., FEM, FDM, DEM and others), based on the model shown in this paper.

6. Conclusion

We simply reviewed on dislocation modeling and on its applications to geological problems (including the disadvantages of the model), and make the following two points.

1. When discussing the essential aspects of tectonics, simply from a distribution pattern of structures caused by fault motions, the dislocation modeling was a very useful tool and should be used with an assumption of a high Poisson's ratio. However, it is difficult to directly compare the displacement amounts between the modeled structures and the actual structures.

2. It is useful to superimpose analytical solutions for different fault parameters on a single fault, in order to introduce the history of fault activity into the numerical model.

We then employed the suggested dislocation modeling technique for the restoration modeling of sedimentary basins formed from 48 Ma to 12 Ma (located in central Hokkaido, Japan), and evaluated the fault types (lateral faulting or reverse faulting). As a result it was found that:

1. Sedimentary basins that were formed from 48 to 12 Ma in central Hokkaido can be explained by the formation of pull-apart basins, due to right lateral motions before 13 Ma, and by reverse motions of the Hidaka-North fault zone and the Hidaka-South fault zone after 13 Ma. Although this makes us expect a significant change in the regional tectonic stress field during this stage, its source should be investigated accurately and quantitatively in the future.

2. The distribution pattern of sedimentary basins restored in this study is similar to the actual distribution of buried sedimentary basins formed during the Paleogene in this area.

Finally, we discovered the following two points relating to the Bouguer gravity anomaly in central Hokkaido:

1. A gravity low area less than 20 mGal corresponds to the area where the Paleogene strata are distributed under the surface.

2. The Bouguer gravity anomalies in the gravity low belt less than 20 mGal can be roughly explained by the subsurface structure model that the mantle around the center of the belt was lifted upwards. Although it is considered that the conspicuous subsidence caused by large lateral motion would induce the crustal thinning-mantle uplift, this possibility should be discussed more accurately and quantitatively in the future.

Acknowledgments

This study was partially supported by the First Bank of Toyama Scholarship foundation and a Grants-in-Aid for Scientific Research (No. 21671003). We are most grateful to Ana Pantar and Book Editors for their editorial advices and cooperation.