Cyclostratigraphy of the Late Miocene to Pliocene sediments at IODP sites U1425 and U1430 in the Japan Sea and paleoceanographic implications

High-resolution studies of sedimentary sequences have reconstructed paleoceanographic changes in the Atlantic and southern and equatorial Pacific Oceans during the Late Miocene and Pliocene, but comparable analyses are lacking for the North Pacific Ocean. However, continuous samples of hemipelagic sequences covering this time interval were obtained at Integrated Ocean Drilling Program sites U1425 and U1430 in the Japan Sea during expedition 346. Because the paleoceanography of the Japan Sea was sensitive to glacio-eustatic sea-level changes, changes in the sediment record are manifested as cyclic lithological changes. By using a gamma ray attenuation (GRA) density, which reflects biogenic silica content, as an indicator of sea-level changes, we developed an independent orbitally tuned chronology at each site to construct high-resolution, high-precision age models for the time span of 11.8 to 1.45 Ma. First, the 405-kyr-filtered GRA profile at each site was tuned to 405-kyr-filtered orbital eccentricity. Then, using the 405-kyr-tuned age model, the 100-kyr-filtered GRA profile was tuned to short-eccentricity cycles to adjust the 405-kyr-tuned age models. We used the resulting age models to revise the time intervals of previously published lithological units and paleoceanographic stages for the Japan Sea. Our profiles based on physical properties are a good match to the oxygen isotope record, except for 6.6–3.6 Ma, a time interval with small-amplitude fluctuations in the oxygen isotope record.


Introduction
The Late Miocene to Pliocene time encompassed the expansion and stabilization of the East Antarctic ice sheet, the onset of northern hemisphere glaciation and the intensification of the Asian Monsoon system (e.g., Zachos et al. 2001;De Vleeschouwer et al. 2017, Tada et al. 2016Betzler et al. 2018;Clift 2017;Holbourn et al. 2018). Paleoceanographic changes during this period have been reconstructed at high resolution on the basis of deep-sea sediments mainly from the Atlantic and southern and equatorial Pacific oceans (e.g., Bell et al. 2014Bell et al. , 2015Westerhold et al. 2005;Holbourn et al. 2005Holbourn et al. , 2007Holbourn et al. , 2013aHolbourn et al. , 2013bHolbourn et al. , 2018Drury et al. 2016). This research has shown that changes in Miocene and Pliocene climate and oceanography were closely associated with orbital-scale changes in ice volume and the carbon cycle (Holbourn et al. 2013a(Holbourn et al. , 2013bDrury et al. 2016). However, there are few continuous, long-term sedimentary sequences of this age suitable for high-resolution analyses from the northwest Pacific Ocean because the calcium carbonate compensation depth was shallow in the area at that time (e.g., Berger 1967;Berger et al. 1982). For that reason, paleoceanographic conditions of the Pacific Ocean during this period are not well understood.
The Japan Sea has been a semi-enclosed marginal sea since at least the Early Miocene (Iijima and Tada 1990). Its connection to the Pacific Ocean through relatively shallow straits has made the Japan Sea sensitive to glacio-eustatic sea-level changes, which are manifested as cyclic lithological changes in its sedimentary records (Tada 1994). Therefore, it appears feasible to construct a high-resolution and high-precision age model by the use of cyclostratigraphy based on lithological variations, as inferred from gamma ray data, in hemipelagic sediments from the Japan Sea.
In this study, we constructed an orbitally tuned age model for the hemipelagic sediments of the Japan Sea during the period 11.8-1.45 Ma, reconstructed temporal changes in the physical properties of cores, and compared them with global and regional benthic oxygen isotope records. We utilized core data from Integrated Ocean Drilling Program (IODP) sites U1425 and U1430 in the central and south-central Japan Sea, recovered during expedition 346 (Tada et al. 2015b(Tada et al. , 2015c. We used a revised version of the stratigraphic record derived from the cores at each site ) combined with both shipboard and newly published biostratigraphic data (Kamikuri et al. 2017) to construct high-resolution, high-precision age models at the two sites.
Physical properties measured at high resolution during the cruise were utilized to correlate the sediment records between holes and between coring sites. These physical properties were also used to support orbital tuning because sedimentary cycles were clearly evident in the high-resolution record (e.g., Lourens et al. 2005;Tian et al. 2008).

Oceanographic and geological setting
The Japan Sea, with an area exceeding 10 6 km 2 and a mean depth of 1667 m, is connected to neighboring seas and the Pacific Ocean through narrow and shallow straits (Fig. 1). The Tsushima Warm Current (TWC), the only current flowing into the sea, controls its oceanographic conditions (Gamo et al. 2014). The TWC Fig. 1 A bathymetric map of the Japan Sea modified from Tada et al. (2015bTada et al. ( , 2015c. Red circles show the locations of the studied sites U1425 and U1430 drilled during IODP exp. 346 and white circles indicate the drilling sites during ODP leg. 127/128. Warm currents flowed into the Japan Sea are also shown as red arrows. FBTWC, SBTWC, and TBTWC indicate the first, second, and third branch of the Tsushima warm current (TWC), respectively Kurokawa et al. Progress in Earth and Planetary Science (2019) 6:2 Page 2 of 20 splits into three branches as it enters the Japan Sea (Hase et al. 1999). The first branch originates in the eastern channel of Tsushima Strait and flows along the west coast of Honshu Island. The second branch originates in the western channel of Tsushima Strait and flows initially along the southern Korean Peninsula, then along the shelf edge to the west of and generally parallel to the first branch. However, this branch is also strongly influenced by eddies generated by baroclinic instability (Hase et al. 1999). The third branch originates in the western channel of Tshushima Strait, flows northward along the Korean Peninsula, and then turns to the east around 36°N and forms the Subpolar Front (Kim and Yoon 1999). In the northwestern Japan Sea, cooling and freezing caused by the winter monsoon increases the subsurface water density (Talley et al. 2003), which forms the deep Japan Sea Proper Water, characterized by its high dissolved oxygen content (Gamo et al. 1986). However, oceanographic conditions have been significantly different in the past, and the sediments of the Japan Sea give evidence of the basin's strong sensitivity to regional and global climatic and oceanographic changes (Tada 1994).
The Japan Sea first opened approximately 28 Ma, and the deposition of hemipelagic sediments began around 20 Ma (Tamaki et al. 1992;Tada 1994). In their revised lithostratigraphy of Japan Sea hemipelagic sediments obtained during Ocean Drilling Program (ODP) leg 127, Tada and Iijima (1992) divided the middle Miocene to Holocene sedimentary sequence into six lithological units. Of these, unit 4 (15.5-10.5 Ma) is composed of gray, highly siliceous rocks and has two subdivisions. Subunit  is composed of claystone with occasional parallel laminations and flattened horizontal burrows, and subunit 4A (14-10.5 Ma) consists of alternating decimeter-scale layers of dark chert and light claystone, the claystone being more strongly bioturbated than the chert. Unit 3 (10.5-6 Ma) is characterized by alternating decimeter-scale layers of non-bioturbated dark diatomaceous silty clay and bioturbated light silty clay. Unit 2 (6-2.5 Ma) is composed of heavily bioturbated, homogeneous greenish diatom ooze and diatomaceous clay. Unit 1 (2.5-0 Ma) is composed of alternating centimeter-to meter-scale layers of dark and light clay to silty clay that are biosiliceous or biocalcareous to various degrees (Tada and Iijima 1992;Tada 1994). Tada (1994) reconstructed the circulation of the paleo-Japan Sea from oceanographic information (surface and bottom water conditions, calcium carbonate compensation depth, surface productivity, channel position, and sill depth) and divided the evolution of the Japan Sea into five stages corresponding to the deposition of the five uppermost lithological units of Tada and Iijima (1992). During stages 2 and 3 (15.5-6.5 Ma), which correspond to the deposition of units 4 and 3, respectively, the surface water was predominantly cool and the deep water fluctuated between suboxic and anoxic (Tada 1991(Tada , 1994. Tada (1991) interpreted cyclic meter-scale lithological changes and apparently annual parallel laminations in the Onnagawa Formation in northeastern Honshu (equivalent to subunit 4A) to argue that terrigenous input fluctuated periodically due to orbitally driven sea-level oscillations. During stage 4 (6.5-2.5 Ma), corresponding to the deposition of unit 2, cold and oxic surface water flowed into the Japan Sea from the north and formed deep water. The consequent upwelling resulted in higher productivity in the surface water and greater deposition of biogenic silica. Stage 5 (2.5-0 Ma), corresponding to the deposition of unit 1, was characterized by large and frequent oscillations of oceanographic conditions. During interglacial sea-level highstands, the TWC was strong and resulted in the production of deep water within the Japan Sea, such that the bottom water became hyperoxic. During glacial lowstands, the Japan Sea became isolated and the TWC was restricted, leading to density stratification beneath the low-salinity surface water and the development of euxinic bottom water (Tada 1994).

Methods/Experimental
Sites and materials IODP site U1425 (39°29.44′ N, 134°26.55′ E, water depth 1909 m) is on a terrace in a northeast-southwest-oriented graben in the middle of Yamato Bank, and site U1430 (37°54.16′ N, 131°32.25′ E, water depth 1072 m) is about 300 km to the southwest on the southern upper slope of the eastern South Korean Plateau (Tada et al. 2015a(Tada et al. , 2015c (Fig. 1). The location of site U1425 is now under the influence of the third branch of the TWC at the Subpolar Front, and site U1430 is slightly south of the Subpolar Front (Hase et al. 1999). Miocene sediments were recovered at these two sites during expedition 346 (Tada et al. 2015b(Tada et al. , 2015c. The sedimentary sequence at site U1425 covers the time interval from the Late Miocene to the Holocene (Fig. 2a) without any obvious hiatus and is dominated by clay, silty clay, diatomaceous ooze, and siliceous claystone (Tada et al. 2015b). The sequence has been divided into three lithological units (units I to III) on the basis of sediment composition, in particular, the abundance of biosiliceous material (Tada et al. 2015a). At site U1430, the sedimentary sequence extends from the middle Miocene to the Holocene (Fig. 2b) and has been divided into lithological units I to IV (Tada et al. 2015c). Between 7 Ma and approximately 4.5-5 Ma, sedimentation was extremely slow or absent at this site (Tada et al. 2015c;Kamikuri et al. 2017). Kurokawa et al. Progress in Earth and Planetary Science (2019)  ray attenuation (GRA), and color reflectance (Tada et al. 2015b(Tada et al. , 2015c. On the shore, core photographs were examined to identify small faults and disturbances that could not have been identified in the physical property profiles. The initial splice was revised on the basis of several new tie points within core intervals corresponding to lithological units I and III to correct miscorrelations and avoid intervals containing sponge plugs and other discontinuities (detailed in Irino et al. 2018). This study adopted the revised splices, in which the updated core composite depth below sea floor (CCSF-D) is referred to as the revised CCSF-D in this paper, and all core depths referred to hereafter are revised CCSF-D depths. The revised CCSF-D for site U1425 is CCSF-D patched rev20170309 in Irino et al. (2018), and the revised CCSF-D for site U1430 is CCSF-D patched rev20160113 in Irino et al. (2018). To check the validity of these revised splices, we correlated the revised composite NGR data with the hostile standard gamma ray (HSGR) logging data from sites U1425 and U1430. Downhole logs are unaffected by core expansion because they were acquired in situ, and their records are perfectly continuous with depth. For this reason, core-log integration helps in identifying miscorrelations that may be present (e.g., deMenocal et al. 1992;Lofi et al. 2016).
In our comparison of the revised NGR splices and the downhole HSGR logs (acquired in hole B at each site), we found strong correlations at the meter scale (Fig. 3). This confirms the validity of the revised splice as there were no significant core gaps deeper than~80 m at sites U1425 and U1430. Correlation points between NGR and HSGR at sites U1425 and U1430 are listed in Tables 1 and 2. Correlation points occur every~5 m in the record from both sites. This resolution ensured that no orbital eccentricity cycles were missing.
It should be noted that gamma ray attenuation by the steel hole casing precluded such a comparison in the upper~80 m at both sites. However, three holes were drilled in this depth range at both sites, such that the revised splices were considered reliable despite the lack of core-log integration.

Chronology
To produce an orbitally tuned age model, we used the shipboard biostratigraphic data described in Tada et al. (2015bTada et al. ( , 2015c for sites U1425 and U1430 and additional radiolarian datums from Kamikuri et al. (2017).
The biostratigraphic age model for site U1425 extended from 9.3 Ma, the oldest age constraint within the  All biostratigraphic datums used are plotted in Fig. 2 and listed in Additional file 1. The following datums were discarded. The last occurrence (LO) and first occurrence (FO) of the planktonic foraminifer Globorotalia ikebei were determined by its rare occurrence in only one sample, and its presence was attributed to reworking (Tada et al. 2015b). The LO of the diatom Thalassiosira schraderi had a relatively large age uncertainty and was not used (Tada et al. 2015b). The LO of the radiolarian Hexacontium parviakitaensis in the Japan Sea was estimated as 2.87-2.90 Ma by Kamikuri et al. (2017), who placed this datum at 86.99 m, whereas the trend from our estimated sedimentation rates suggested that the LO of H. parviakitaensis should be at 81.01 m, and we could not rule out reworking. Finally, the FO of the radiolarian Cycladophora sphaeris (=Cycladophora sakaii in Kamikuri et al. 2004Kamikuri et al. , 2007 has age uncertainties of~0.5 Myr between the eastern North Pacific, western North Pacific, and the Bering Sea (Kamikuri et al. 2004(Kamikuri et al. , 2007. The biostratigraphic age model for site U1430 extended from 11.8 to 1.45 Ma (Fig. 2b). The LO of Hexacontium parviakitaensis was excluded for consistency with site U1425. The linear sedimentation rates (LSR) for the stratigraphic interval corresponding to 5 to1 .45 Ma were too low (~1 cm/kyr) to detect meaningful orbital-scale changes in GRA, and the radiolarian datums and the occurrence of a glauconite layer suggest a hiatus during the time interval from 7 to 5 Ma (Kamikuri et al. 2017;Matsuzaki et al. 2018;Tada et al. 2015c). For this reason, we conducted an orbital tuning only for 11.8 to 7.3 Ma.

Cyclostratigraphy
Quaternary sediments in the Japan Sea contain alternating dark-and light-colored layers that were synchronously deposited on the millennial scale, providing a tool for high-resolution correlation ). An age model for the upper part of unit I (1.45-0 Ma) was established at site U1424 by Tada et al. (2018) by correlating the onboard GRA record of Tada et al. (2015bTada et al. ( , 2015c) (regarded as a proxy of sea level) with the LR04 oxygen isotope (δ 18 O) sequence of Lisiecki and Raymo (2005), which allows ages of each of the dark layers at site U1424 to be extended to other sites. The orbitally tuned age model proposed here thus applies to the sedimentary sequence older than 1.45 Ma. The Japan Sea sediments are composed of biosiliceous clay to silty clay (Tada et al. 2015b(Tada et al. , 2015c. In biosiliceous sediments, GRA reflects the ratio of detritus to diatoms and decreases as diatom content increases because diatom frustules have high internal porosity (Tada and Iijima 1983). Diatom abundance is strongly linked to sea level because the surface productivity of the Japan Sea is controlled by the influx of nutrients through Tsushima Strait, which in turn has been modulated by glacio-eustatic sea-level changes ). GRA likewise increases as the input of detritus increases. During sea-level lowstands, for instance, more terrigenous sediment should reach the study sites in the Japan Sea as the exposed area of the paleo-Japanese islands increases and the distance from the shoreline decreases, and the consequent increase in the mass accumulation rate of detritus (Tada 1991) results in higher GRA values. Thus, GRA can be used as a sea-level proxy.
To use GRA as a proxy, noise was first removed from the composite onboard GRA data (Tada et al. 2015b(Tada et al. , 2015c. GRA values below the specific gravity of seawater (1.02 g cm − 3 ) were discarded, as were data corresponding to tephra and cracked intervals. The filtered GRA profiles for sites U1425 and U1430 are presented in Additional file 2. We next performed a spectral analysis of the filtered GRA data to assess their fit to our tuned biostratigraphic age model using REDFIT software (Schulz and Mudelsee 2002). We generated Gaussian-filtered GRA time-series for both sites and assessed their correlation with temporal variations of Earth's orbital parameters, including long-eccentricity cycles (405 kyr) and short-eccentricity cycles (100 kyr) calculated by Laskar et al. (2004). For the Gaussian band-pass filtering, we used AnalySeries 2.0.8 software (Paillard et al. 1996) to sample evenly spaced data points from unevenly spaced original data. Band-pass filters were centered at 0.0025 kyr −1 (0.0005 kyr −1 bandwidth) and at 0.01 kyr −1 (0.003 kyr −1 bandwidth) to extract 405-and 100-kyr cycles, respectively. Tuning assumed no phase lag between GRA variability and orbital parameters, because the response time of an ice sheet is not known in this time interval (Holbourn et al. 2007).
To further improve our orbitally tuned age model, we tuned the GRA data from site U1425 to the oxygen isotope record, which is sensitive to sea-level fluctuation. The site U1425 record was chosen for its continuity and for its wide range of variability, which aided the tuning process. We tuned these data to the LR04 oxygen isotope record of Lisiecki and Raymo (2005)

Results: orbital tuning
A critical part of orbital tuning is the choice of the initial tie point. For site U1425, we used the dark layer no. 21-2 projection tie points of Tada et al. (2018), which are at 1.45 Ma at 54.79-m-core depth. These have an error margin of about 4 kyr, inherited from the LR04 age model used as a tuning target. This error is negligible The results of our spectral analysis show that the GRA data at both sites have a characteristic cyclicity that can be correlated with orbital cycles. At site U1425, peaks around 0.0025 and 0.01 that are significant above the 80% confidence level correspond to the long and short-eccentricity cycles, respectively, and peaks around 0.024 may be attributable to obliquity cycles (Fig. 4a). At site U1430, the power spectrum also shows a clear peak at 0.0025, a peak at 0.0057 that corresponds to the long obliquity cycle (174 kyr), and peaks around 0.024 that correspond to the 41-kyr obliquity cycle, although a short-eccentricity peak is absent (Fig. 4b). In sum, the GRA data displayed orbital signals at both sites, which we used to construct orbitally tuned age models.
We constructed an age model for each site in two steps. First, 405-kyr-filtered GRA profiles from our biostratigraphic age model were tuned to the long (405 kyr) eccentricity cycle. Because this cycle has a stable periodicity over the last 250 Myr, it is well suited for cyclostratigraphy (Hinnov 2000;Laskar et al. 2004;Hinnov and Hilgen 2012). Second, we tuned the resulting GRA profiles to the short-eccentricity cycle by correlating eccentricity minima with 100-kyr-filtered GRA maxima under the assumption that sea-level lowstands occur during eccentricity minima that increase the influx of detritus to the Japan Sea and thus raise GRA values, as explained previously.
The 405-kyr eccentricity cycle is relatively weak in benthic oxygen isotope records younger than about 5 Ma (Zachos et al. 2001;de Boer et al. 2014), although both 405-kyr and 100-kyr eccentricity cycles are clear in benthic carbon and oxygen isotope records from the Oligocene and Miocene (Pälike et al. 2006;Liebrand et al. 2011Liebrand et al. , 2017Holbourn et al. 2013aHolbourn et al. , 2013bBeddow et al. 2016). Accordingly, we did not use the 405-kyr eccentricity cycle as a tuning target after 5.0 Ma at site U1425. The plot of 405-kyr-filtered GRA records based on the biostratigraphic age model and the long-eccentricity cycle from 9.3 to 5 Ma at site U1425 shows that the two records are nearly a half-wavelength out of phase (Fig. 5a). We tuned the GRA maxima to the nearest eccentricity minima to yield the result shown in Fig. 5b. There are 10 correlated cycles in this time interval. The mismatch in amplitudes between the GRA signals and eccentricity from 6.2 to 7.5 Ma (Fig. 5b) means that amplitude modulation of GRA by the long-eccentricity cycle is not established for the 9.3-5 Ma interval.
Using an age model based on the long-eccentricitytuned age model for 9.3-5 Ma and the biostratigraphic model for 5-1.45 Ma, we conducted short-eccentricity tuning to correlate the GRA record and the short-eccentricity (100 kyr) cycle. Both GRA and eccentricity were filtered with a 100-kyr filter using Gaussian band-pass filtering (see the "Method" section), and GRA maxima were correlated with eccentricity minima (Fig. 6a and b). The amplitudes of the two records were well matched except from 6.6 to 3.6 Ma, when the overall amplitude of GRA was small regardless of eccentricity (Fig. 6b). The depths and ages of the 405-and 100-kyr correlation points for site U1425 are listed in Tables 3 and 4 and plotted in Fig. 2a, which shows that orbitally tuned correlation points are consistent with the biological datums.
We applied the same procedure to build an orbitally tuned age model for the 11.8-7.3 Ma interval at site U1430, first correlating the 405-kyr-filtered GRA maxima with long-eccentricity minima ( Fig. 5c and d), then refining the age model by correlating peaks in the 100-kyr-filtered GRA and short-eccentricity records ( Fig. 6c and d). The 405-kyr and 100-kyr correlation points are listed in Tables 3 and 4 and plotted against depth in Fig. 2b.
We checked the validity of the tuned results by comparing their power spectra to the power spectra of the biostratigraphic age models (Fig. 4a and f ). At site U1425, GRA peaks after short-eccentricity tuning became stronger at 405, 100, and 41 kyr and exceeded the 95% confidence level, although other peaks were little changed ( Fig. 4c and e). At site U1430, cycles of long eccentricity and obliquity were stronger after orbital tuning ( Fig. 4b and d). Thus, we are confident that our two-step orbital tuning improved the accuracy of the age models at both sites because the orbital signals of the GRA power spectra were stronger across all orbital peaks (Fig. 4). In addition, we tuned the GRA record to the oxygen isotope record at site U1425 to check the validity of our assumption that GRA is sensitive to sea-level changes and to check the validity of the orbital tuning results. The results, given in Additional files 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, and 13, indicate that our orbitally tuned age model is consistent with the δ 18 O-tuned age model with a lag no greater than 40 kyr, or less than half of an eccentricity cycle.
Although the GRA power spectra before orbital tuning had peaks attributable to obliquity and precession, we did not conduct tuning with respect to obliquity or precession because the obliquity and precession cycles have relatively high uncertainties owing to effects of tidal dissipation and changes in dynamical ellipticity, and because GRA amplitudes were too small during some time intervals.

Comparison of physical properties with sea-level change
After adjusting our GRA records on the basis of our short-eccentricity-tuned age model, we compared them to a previously published sea-level record. Miller et al. (2005) compiled two high-resolution benthic oxygen isotope records covering the last 9 Myr from ODP site 846 in the equatorial Pacific (Shackleton et al. 1995) and ODP site 982 in the northern North Atlantic (Hodell et al. 2001) and converted them to a record of eustatic sea level, using a calibration ratio of 0.1‰/10 m. In this record, sea-level fluctuations over glacial-interglacial cycles were roughly 40 m before~6 Ma, gradually increased to l.E = long eccentricity, s.E = short eccentricity, O = obliquity, P = precession. a U1425 biostratigraphy based age. b U1430 biostratigraphy based age. c U1425 long-eccentricity-tuned age. d U1430 long-eccentricity-tuned age. e U1425 short-eccentricity-tuned age, f U1430 short-eccentricity-tuned age Kurokawa et al. Progress in Earth and Planetary Science (2019)  about 100 m by 1.5 Ma, then remained at that amplitude until the present (Miller et al. 2005). For this comparison, we compiled the GRA dataset from both sites, then resampled it and the sea-level record of Miller et al. (2005) at 100-kyr intervals between 9 and 1.45 Ma. The result (Fig. 7) shows a fair correlation (r = 0.63) between the two datasets for this period. Between 6.6 and 3.6 Ma, the GRA amplitudes were small and the estimated sea levels were mostly restricted between 10 and − 20 m. These lower values of estimated sea-level changes may be reflected in the low amplitude of the GRA dataset, and the correlation between GRA and estimated sea-level amplitude was slightly stronger when we excluded the data from 6.6 to 3.6 Ma (r = 0.65).
DeVleeschouwer et al. (2017) compiled a high-resolution benthic oxygen isotope record for the entire Neogene and Quaternary from multiple drill sites in the Atlantic, eastern equatorial Pacific, and South China Sea. We used this record to compare the variabilities of NGR and GRA at sites U1425 and U1430 with the oxygen isotope record of benthic foraminifers since 12 Ma. The GRA and NGR data were well correlated with the benthic oxygen isotope record before 8 Ma and after 3.3 Ma (Fig. 8), including the last 1.45 Ma ; however, the correlation was unclear between 8 and 3.3 Ma, when both datasets had orbital-scale variations of small amplitude.
In the benthic oxygen isotope record of Miller et al. (1991), excursions Mi-5 to Mi-7 occurred during the time interval covered by the Japan Sea sediments examined in this study. In their orbitally tuned age model based on Fe signals in South Atlantic shallow marine sediments, Westerhold et al. (2005) recognized the Mi-5 and Mi-6 events at 11.7 and 10.7-10.4 Ma, respectively. At site U1430, we recognized GRA and NGR maxima corresponding to these two events at 11.7 and 10.6 Ma, respectively (Fig. 8). The Mi-6 event is also clearly . Although the age of the Mi-7 event is uncertain because the oxygen isotope signal has many small maxima between 9.6 and 9.0 Ma (Andersson and Jansen 2003;Westerhold et al. 2005), GRA and NGR maxima in the Japan Sea sediments are compatible with two candidate ages proposed by Westerhold et al. (2005), 9.4 and 8.7 Ma (Fig. 8). The close correspondence between physical properties and sea-level changes such as the Mi events tends to confirm the accuracy of short-eccentricity-tuned age models. In addition, our results support the possibility that the amplitude of sea-level changes can be estimated from GRA data independent of oxygen isotope values.
Correlation between sites U1425 and U1430 Figure 8 shows the temporal variations of NGR, GRA, and LSR for the last 12 Myr at sites U1425 and U1430. The sedimentary record at site U1425 extends back to 9.3 Ma without any stratigraphic gaps, whereas the sedimentary record at site U1430 from 11.8 to 0 Ma is interrupted by a probable hiatus from~7.3 to 4.5-5 Ma. Thus, comparisons of the two records are limited to the Fig. 6 Correlation between 100-kyr-filtered GRA and 100-kyr-filtered eccentricity. a, c The 100-kyr-filtered GRA (blue) based on long-eccentricitytuned age model at sites U1425 and U1430. Eccentricity is filtered at 100-kyr (red). b, d The 100-kyr-filtered GRA (blue) based on short-eccentricity (100-kyr)-tuned age model at sites U1425 and U1430. The 100-kyr-filtered GRA maxima are tuned to 100-kyr-filtered eccentricity minima  periods 9.3-7.3 and 5-1.5 Ma (Fig. 8). We therefore compiled NGR, GRA, and LSR time-series using the age model of Tada et al. (2018) for the uppermost part of both sites, biostratigraphic data for the 5.0-1.5 Ma interval at site U1430, and our short-eccentricity-tuned age models (Fig. 2). In general, temporal changes in GRA and NGR were well correlated between the two sites from 3.5 to 0 Ma. From 5.3 to 3.5 Ma, GRA and NGR values were higher at site U1430 than at site U1425 (Fig. 8), a difference that can be explained by the fact that the sediment at site U1430 is less diatomaceous or more tuffaceous than at U1425 (Tada et al. 2015b(Tada et al. , 2015c. From 9.3 to 7.3 Ma, GRA and NGR curves at the two sites were similar. The close correspondence in physical properties suggests that paleoceanographic conditions were similar at sites U1425 and U1430. During the hiatus interval at site U1430 from 7.3 to 5.3 Ma, there were no significant changes in physical properties and LSR at site U1425 (Fig. 8) or in LSRs at ODP sites 794, 795, and 797 also located in the Japan Sea (Tada 1994). The Late Miocene hiatuses have been reported from some on-land sections in northeastern Japan (Watanabe 1994), but these correspond to the diatom datums 5D, 6A, and 6B (Akiba 1986;Yanagisawa and Akiba 1998;Tada et al. 2015aTada et al. , 2015c, which range in age from 9.3 to 7.6 Ma. The hiatus at site U1430 coincides in depth with a major regional unconformity on seismic profiles (MB3 unconformity in Horozal et al. 2017) that extends over much of the eastern South Korea plateau where site U1430 was drilled. The origin of this unconformity is unclear. Watanabe (1994) and Horozal et al. (2017) attributed it to erosion by bottom currents or a brief episode of tectonic uplift. Whatever its origin, this hiatus must reflect a local, basin-scale event because there are no drastic changes in LSRs in this time interval in other Japan Sea cores. Fig. 7 Relationship between GRA and estimated sea level. GRA data are from sites U1425 and U1430 and resampled every 100 kyr. Sea level is converted from oxygen isotope values using a calibration ratio of 0.1‰/10 m (Miller et al. 2005). Estimated sea level is also evenly sampled every 100 kyr. Regression line is applied to the all data during the last 9 Myr. Red circles are the data from 6.6 to 3.6 Ma Fig. 8 Comparison of NGR, GRA, and LSR with composite benthic isotope curve from DeVleeschouwer et al. (2017) (purple). Vertical green lines indicate good correlation between physical properties and oxygen isotope, while yellow lines are indicated that correlation is unclear. Marine isotope stage number from Lisiecki and Raymo (2005) and Miller et al. (1991) is also shown in the top panel. Variations at site U1425 are shown as blue and U1430 as red. Lithological unit (Tada et al. 2015b(Tada et al. , 2015c and paleoceanographic condition stage number (Tada 1994) are also shown and that positions are revised in this study. Parallel lamination positions are shown as blue bars at site U1425 and red bars at site U1430 Revision of lithological divisions at sites U1425 and U1430 Tada et al. (2015bTada et al. ( , 2015c) defined a set of lithological units in Japan Sea sediments, of which the latest four (I to IV) and their subunits are represented in the sediments treated in this study. The revised GRA, NGR, and logging profiles and the revised age model constructed in this study lead us to propose revisions of the boundaries of these units and subunits. Tada et al. (2015bTada et al. ( , 2015c stated that unit I is characterized by decimeter-scale intervals of alternating light and dark sediment and that the boundary between subunit IA and the underlying subunit IB is distinguished by a downward decrease in the frequency of these alternations. They placed the base of subunit IB at 50.7 m at site U1425 and at 46.9 m at Site U1430; however, this change in frequency occurs deeper, at 54.8 m at site U1425 and 50.7 m at site U1430, as shown by lithological observations as well as NGR and GRA data (Fig. 8). The same discrepancy is also observed at other IODP expedition 346 sites in the Japan Sea Irino et al. 2018). The age of the sediment at these revised boundary positions is 1.45 Ma . The NGR and GRA profiles in subunit IA (Fig. 8) display large, suborbital-scale fluctuations after 1.45 Ma that correspond to the color banding in the sediments. The base of subunit IB is defined by a downward change in sediment color and increase in diatom content (Tada et al. 2015b(Tada et al. , 2015c. The increase in diatom content corresponds to a clear decrease in NGR at 90.71 m at site U1425,~1.5 m shallower than the depth specified by Tada et al. (2015bTada et al. ( , 2015c, and 61.14 m at site U1430, the same depth specified in the original definition. The age of this boundary is 2.61 Ma according to our revised age model (Fig. 8). During deposition of unit IB, NGR and GRA fluctuated strongly, as they did in unit IA, but at a lower frequency than in unit IA.
In unit II, the boundary between subunits IIA and IIB is defined by a downward increase in diatom content, which is recorded as a downward decrease in NGR or GRA at 137.34 m at site U1425 (4.14 Ma in our age model) and 74.42 m at site U1430 (4.04 Ma in our age model). NGR data do not show a clear downward decrease at 74.42 m at site U1430, but this probably reflects the proximity of an NGR peak from a glauconite layer (e.g., Tada et al. 2015c) at~78.8-79.7 m. The age of the boundary differs slightly between the two sites; however, the age of the upper hiatus at site U1430 is based on a biodatum, which has an error of about ± 0.1 Myr. Accordingly, the boundary positions are consistent between the sites. Tada et al. (2015bTada et al. ( , 2015c distinguished unit III from unit II by its lower diatom content and darker sediment color, which is reflected in higher values of NGR and GRA (Fig. 8). The boundary between units II and III (specifically, subunit IIB and subunit IIIA) is at 263.25 m at site U1425 and 84.81 m at site U1430, and its age is 7.36-7.39 Ma according to our age model (Fig. 8). Subunit IIIA is characterized by higher NGR and GRA values, with fluctuations of greater amplitude, than unit II reflecting the increased contrast between alternating light and dark layers.
In summary, our analysis shows that the revised lower boundaries of the lithological units of Tada et al. (2015bTada et al. ( , 2015c correspond to the following ages: 1.45 Ma for subunit IA, 2.61 Ma for subunit IB, 4.04-4.14 Ma for subunit IIA, and 7.36-7.39 Ma for subunit IIB.
Lithostratigraphic correlation of Japan Sea sites and refinement of paleoceanographic stages The criteria for defining lithological units were basically the same for ODP leg 127/128 and IODP expedition 346. According to Tada and Iijima (1992), unit 1 in ODP leg 127/128 sediments is characterized by decimeter-to meter-scale alternations of dark and light silty clay to clay and is divided into subunits 1A and 1B on the basis of the frequency of dark layers and the intensity of bioturbation. Unit 2 is composed of clayey diatom ooze and diatomaceous clay that is heavily bioturbated and mottled, with alternating dark and light layers and a general darkening of sediments near its base. This unit is not subdivided. Unit 3 is characterized by decimeter-to meter-scale alternating light and dark layers, the latter being faintly laminated, and by a relatively high organic carbon content. Tada (1994) compiled the paleoceanographic data of ODP leg 127/128 and paleogeographic data on the Japanese islands by Iijima and Tada (1990) to reconstruct the paleoceanographic evolution of the Japan Sea. The proposed stages of this history correspond to the lithological units of Tada and Iijima (1992). Given that correspondence, paleoceanographic stages 5, 4, and 3 of Tada (1994) are shown in Fig. 8 as coeval with lithological units I, II, and III, respectively.
At sites U1425 and U1430, the boundary between units III and II is based on the increase in diatom content and decrease in terrigenous content and is assigned an age of~7.4 Ma. Tada (1994) described a dramatic increase in biogenic silica content around 8 Ma at the ODP leg 127/128 sites, which is probably correlated with this boundary. According to Tada (1994), during paleoceanographic stage 3, bottom water conditions fluctuated between anoxic and suboxic in association with eustatic sea-level changes. This history is consistent with the sporadic occurrence of parallel laminations in lithological subunit IIIA at sites U1425 and U1430, which suggests the suppression of bioturbation due to low oxygen levels (Fig. 8). The difference in the ages of parallel laminations between the two sites may reflect the position of euxinic or anoxic water masses. The presence of parallel laminations only at site U1425 in the upper part of unit IIIA signifies euxinic or anoxic conditions in the bottom water there and oxic or suboxic conditions at site U1430, whereas the presence of parallel laminations only at site U1430 in the lower part of unit IIIA (~9.0 Ma) may indicate the development of an oxygen minimum zone within the Japan Sea, because site U1430 was at the appropriate paleo-water depth. Another line of evidence that an oxygen minimum zone developed or expanded in the Japan Sea around 8 Ma is radiolarian data indicating the presence of North Pacific type deep-water species at that paleo-water depth (Matsuzaki et al. 2018).
Paleoceanographic stage 4 (7.4-2.61 Ma) corresponds to lithological unit II, which is characterized by relatively low values of GRA and NGR at Site U1425. Tada (1994) suggested that during this stage, terrigenous input was diluted by diatoms due to higher productivity caused by enhanced production of deep water and consequent upwelling.
During paleoceanographic stage 5 (2.61-0 Ma), NGR and GRA fluctuated strongly at millennial to orbital scales, interpreted as the result of millennial-scale changes in the East Asian summer monsoon superimposed on glacio-eustatic sea-level changes . The deposition of lithological unit I is thought to have been initiated by a drop in sea level caused by the onset of northern hemisphere glaciation (Lisiecki and Raymo 2005;Tada et al. 2018).