OJMS  Vol.7 No.1 , January 2017
Recent Advances in Japanese Fisheries Science in the Kuroshio-Oyashio Region through Development of the FRA-ROMS Ocean Forecast System: Overview of the Reproducibility of Reanalysis Products
To address various fisheries science problems around Japan, the Japan Fisheries Research and Education Agency (FRA) has developed an ocean forecast system by combining an ocean circulation model based on the Regional Ocean Modeling System (ROMS) with three-dimensional variational analysis schemes. This system, which is called FRA-ROMS, is a basic and essential tool for the systematic conduct of fisheries science. The main aim of FRA-ROMS is to realistically simulate mesoscale variations over the Kuroshio-Oyashio region. Here, in situ oceanographic and satellite data were assimilated into FRA-ROMS using a weekly time window. We first examined the reproducibility through comparison with several oceanographic datasets with an Eulerian reference frame. FRA-ROMS was able to reproduce representative features of mesoscale variations such as the position of the Kuroshio path, variability of the Kuroshio Extension, and southward intrusions of the Oyashio. Second, using a Lagrangian reference frame, we estimated position errors between ocean drifters and particles passively transported by simulated currents, because particle tracking is an essential technique used in applications of reanalysis products to fisheries science. Finally, we summarize recent and ongoing fisheries studies that use FRA-ROMS and mention several new developments and enhancements that will be implemented in the near future.

1. Introduction

In the last decade, many ocean forecast systems that include some sort of data assimilation method have been developed to nowcast and forecast ocean states with a particular focus on mesoscale variability; these include HYCOM/NCODA [1] , Mercator [2] , FOAM [3] , MOVE/MRI.COM [4] , JCOPE2 [5] , and so on. Progress in this area is attributable primarily to the Global Ocean Data Assimilation Experiment project, which, crucially, supports near-real-time collection, processing, and management of oceanographic data [6] [7] . In addition, several state-of-the-art community ocean models, which can be coupled with data assimilation schemes, have been developed as open-source code (e.g., ROMS [8] and ECCO [9] ) and have contributed to making ocean nowcasting and forecasting relatively easy. Consequently, physical oceanographic studies have made rapid progress in the areas of ocean dynamics and predictability.

Utilization of outputs from ocean forecast systems is not restricted to physical oceanographic studies. Data assimilation models initialize ocean state variables for forecasts by the synthesis of spatio-temporally limited oceanographic data and ocean circulation model. These initialized state variables, which can be considered to have been optimally estimated by combining observation and simulation, are generally referred to as “analysis” or “reanalysis” data. Such assimilated products are needed for many scientific and practical applications where realistic oceanographic conditions are required: in marine navigation systems used for transportation safety (e.g., [10] ); for acoustic propagation forecasting (e.g., [11] ); and in investigations of the transport of marine debris (e.g., [12] [13] ), oil spills (e.g., [14] [15] ), and radioactive materials (e.g., [16] [17] ).

In recent years, there has also been an intensive demand for assimilated model products in fisheries science. These products have been applied from both Eulerian and Lagrangian viewpoints to fisheries science with various aims. Assimilated data have been analyzed instead of observed data to clarify the linkage between fish abundance and oceanographic conditions [18] , and they have been used for stock management applications [19] . Survival, growth, migration, and recruitment of important commercial fish have been examined by using particle-tracking models and individual-based models (IBMs) with assimilated oceanographic conditions (e.g., [20] [21] [22] [23] ). Several studies have developed high-resolution nested ocean models with data assimilation products (e.g., [24] ) in an effort to understand and predict the sudden occurrence of harmful marine organisms such as jellyfish, toxic algae, and viruses. Potential fishing areas have also been identified by examining statistical relationships between fish habitats and simulated oceanographic conditions (e.g., Habitat Suitability Index (HSI) [25] [26] ). Additionally, the spatio-temporal variability of nutrients and plankton, which is associated with food availability for fish, has been simulated by an ocean forecast system coupled with a lower-trophic-level ecosystem model (e.g., [27] [28] ).

The Japanese Islands are located in the Kuroshio-Oyashio western boundary current region, where mesoscale variations are prominent (Figure 1(a)). In addition, the islands are characterized by a narrow continental shelf (typically less than 50 km wide),

Figure 1. (a) Schematic view of the Kuroshio-Oyashio system. The main currents are denoted by arrows, and geographical names are enclosed in boxes. 1st and 2nd indicate the first and second Oyashio intrusions, respectively. Contour lines with an interval of 10 cm indicate 20-year mean sea surface height estimated from the AVISO maps of absolute dynamic topography during 1993-2012 (i.e., MDT plus 20-year mean SLA). (b) Model domains and bathymetry (contour). The entire illustrated area is the 1/2-degree model domain, and the region within the dashed rectangle is the 1/10-degree model domain.

with the result that mesoscale variations intermittently affect coastal waters through interactions with submesoscale variations [29] [30] . Moreover, mesoscale variability in the Kuroshio-Oyashio region is a key issue for Japanese fisheries and fisheries science, because many pelagic fish, which are important targets of Japanese fisheries, extensively use the Kuroshio-Oyashio region as spawning and nursery grounds. Among these are the Japanese sardine (Sardinops melanostictus) [31] , Japanese anchovy (Engraulis japonicus) [32] , Pacific saury (Cololabis sairai) [33] , Japanese chub mackerel (Scomber japonicus) [34] , Japanese common squid (Todarodes pacificus) [35] , and various other species [36] [37] . For example, Japanese sardine adults migrate southward from the Oyashio to the Kuroshio region for spawning. During winter and spring, eggs, larvae, and juveniles are transported eastward by the Kuroshio, the Kuroshio Extension, or related mesoscale flows, and then juveniles migrate northward to the Oyashio nursery region via the Kuroshio-Oyashio transition zone [38] [39] . Japanese fisheries capture sardine larvae and juveniles in the Kuroshio region during the passive transport period, and they capture adults mainly in the Oyashio region during the nursery and migration period. Hence, an exhaustive dataset across the Kuroshio-Oyashio region is required to study relationships between oceanographic conditions and individual life stages of even this single pelagic fish.

The mission of our organization, the Japan Fisheries Research and Education Agency (FRA), which is one of the National Research and Development Agencies of Japan, is to conduct a wide range of research and development activities related to fisheries and to address many fisheries problems through basic research and practical applications. To respond to the demands of multiple research projects on fisheries oceanography, management, and economics, FRA oceanographers have had to prepare an operational ocean forecast system to create a basic and versatile oceanographic dataset. Such an ocean forecast system had been also expected to extensively enable Japanese prefectural fisheries institutes as well as FRA to systematically conduct many kinds of fishery researches and establish many practical applications, which contribute significantly to sustainable fishery and safety of fishery products around Japan. Although the Japanese Oceanographic Society has several operational ocean forecast systems, including MOVE (Multivariate Ocean Variational Estimation) [4] , JCOPE2 [5] , and DREAMS [40] , none of these forecast systems is specialized to fisheries science. Therefore, the FRA has developed an ocean forecast system based on the Regional Ocean Modeling System (ROMS) [41] and specialized to fisheries science called FRA-ROMS. Here, some readers may question why and how FRA-ROMS is specialized for fisheries science. These questions are addressed in Section 4, where ongoing community efforts of the FRA are introduced.

The configurations of FRA-ROMS are described in Section 2, and the reproducibility of its reanalysis products is examined in terms of both Eulerian and Lagrangian reference frames in Section 3. In Section 4 we summarize recent and ongoing fisheries studies that have used FRA-ROMS, and discuss prospects for work in the near future.

2. Configurations of FRA-ROMS

FRA-ROMS consists of an ocean circulation model based on ROMS and three-dimen- sional variational (3D-Var) objective analysis schemes, which are not included in the original ROMS code. The configurations of FRA-ROMS are briefly explained below.

2.1. Ocean Circulation Model

FRA-ROMS consists of 1/2- and 1/10-degree ocean models connected by one-way nesting, in which volume transport across the lateral boundaries is adjusted according to the method proposed by [42] . The 1/2-degree parent model covers almost the entire North Pacific region (Figure 1(b)). This model was designed to simulate basin-scale variations associated with El Niño―Southern Oscillation and the Pacific Decadal Oscillation and impose them at the lateral boundaries of the 1/10-degree child model as external forcings. Thus, daily mean outputs from the 1/2-degree model are imposed at the lateral boundaries of the 1/10-degree model. The domain of the 1/10-degree model is limited to the Northwestern Pacific (Figure 1(b)) to simulate mesoscale variations over the Kuroshio-Oyashio region. Vertical positions in both models are defined by the specific S-coordinate function proposed by [43] , and the vertical resolution of both models is 48 levels. Since a sea-ice model is not coupled with the current version of FRA-ROMS, simulated temperature and salinity (TS) values over the Okhotsk Sea and the Bering Sea in the two models are restored to the monthly mean climatology obtained from World Ocean Atlas 2001 (WOA2001) [44] [45] on a timescale of 30 days. Numerical schemes and parameterizations that have been adapted for the two models are summarized by [29] .

The two models were integrated for 20 years from the WOA2001-based initial conditions by using climatological monthly mean heat and momentum fluxes at the sea surface [29] . Then, both models were driven by 6-hourly mean atmospheric and radiation elements (i.e., wind speed, air temperature, sea level pressure, specific humidity, precipitation, and net shortwave and downward long wave radiations) from the Japanese 25- year Reanalysis (JRA-25) dataset [46] , from which heat and momentum fluxes were sequentially estimated by the COARE (Coupled-Ocean Atmosphere Response Experiment) bulk formula [47] during each model run.

2.2. Data Assimilation Based on 3D-Var Analysis

To perform reanalysis experiments for the period from 1993 to 2012, 3D-Var analysis schemes were constructed by a technique similar to that used by MOVE, an ocean forecast system developed by the Meteorological Research Institute [4] . 3D-Var analysis schemes were originally developed for a z-coordinate ocean model, but they can be stably applied to a sigma-coordinate model such as the Princeton Ocean Model [5] . Therefore, this scheme can be expected to work robustly when coupled with S-coor- dinate ocean models such as ROMS. The 3D-Var analysis schemes applied here should also create reanalysis data more precise than two-dimensional schemes (e.g., Optimal Interpolation). Meanwhile, the 3D-Var schemes may be less perfect than four-dimen- sional (4D-Var) schemes and ensemble Kalman Filter, because background errors in the 3D-Var schemes are implicitly assumed to be stationary. In addition, heat and salt budget in reanalysis data based on the 3D-Var schemes are not completely conserved, as explained below. Nevertheless, we selected the 3D-Var analysis schemes, which had an advantage that computational cost was much smaller and were more speedy/prac- tical for fisheries science, where additional heavy computations such as PTMs are needed after data assimilation.

The following cost function () was defined and minimized by a nonlinear optimization algorithm without calculating the inverse of the covariance matrix [4] [48] [49] [50] [51] .


where, , and are column vectors whose elements are the state variables, the first guess, and the observed values, respectively. is the horizontal background correlation matrix, is the observed error covariance matrix, and is an operator that picks up values of temperature and/or salinity at the observation position. h denotes the dynamic topography estimated from, is the altimetry-derived sea- level anomaly (SLA), and is the standard SLA observation error. A vertically coupled TS empirical orthogonal function (EOF) modal decomposition was adapted for. The 10 dominant modes of the coupled TS EOF were estimated prior to the objective analysis from six sub-regions of the 1/2-degree model and nine sub-regions of the 1/10-degree model [4] .

The fourth term on the right-hand side of Equation (1) is a constraint that was proposed by [51] to avoid extremely low temperature estimates. We confirmed that without this term extremely low temperatures appeared frequently in the reanalysis data, particularly in the Kuroshio-Oyashio transition zone (not shown). The low temperatures appear because without this term the cost function basically assumes a Gaussian distribution of TS errors, even when the temperature is around 0˚C and close to the freezing point of seawater (i.e., about −1.7˚C).

Optimal TS values estimated with Equation (1) are referred to as analysis data in this paper. The analysis data were assimilated into the ocean model by incremental analysis updating [52] on a 7-day cycle as follows. The ocean model was freely integrated from an initial value for 3.5 days, and the simulated TS values at the end of the 3.5 days were used as first guess values in Equation (1). After optimal TS values (i.e., analysis data) were obtained, the difference in TS between the first guess and the analysis data was estimated. This difference was divided by the number of internal-mode time steps, each corresponding to 7 days, and the result was added to the TS governing equations as an increment. The ocean model was again integrated from the same initial value for 7 days with the increment as forcing (hereafter, forced run). These procedures were repeated sequentially on a 7-day cycle. Assimilated data were sea surface temperature (SST), altimetry-derived SLA, and in situ TS profiles (Table 1). In this study, reanalysis experiments were conducted for 1993-2012, and daily mean outputs during the forced run (i.e., reanalysis data) were analyzed.

It should be borne in mind that FRA-ROMS has two different versions, an operational version and a long-term reanalysis experimental version. The former has been used operationally since May 2012 for nowcasts and two-month forecasts at intervals of

Table 1. List of assimilation data.

7 days by assimilating available near-real-time data. The output is updated every week and posted on the FRA-ROMS website (http://fm.dc.affrc.go.jp/fra-roms/index.html). Meanwhile, the latter is used for reanalysis experiments aimed at accurate reproduction of past oceanographic conditions by assimilation of high-quality delayed-mode data, which have been massively accumulated. The experimental version can be used more effectively for examination of linkages between oceanographic and fisheries/fish conditions in past time, whereas the operational version is more effective for examining linkages in near real time. The reanalysis data presented here were produced by the experimental version and are expected to be more accurate than the products of the operational version.

3. Reproducibility of the Reanalysis Product

We examined the reproducibility of the daily reanalysis data for the 20 years from 1993 to 2012 by comparing the simulated results with previous oceanographic knowledge and observed data over the Kuroshio-Oyashio region in both Eulerian and Lagrangian reference frames.

3.1. Analyses from an Eulerian Viewpoint

We examined the main flows in the Kuroshio-Oyashio region by mapping stream functions of the simulated 20-year mean volume transport from the sea surface to a depth of 1000 m or the sea bottom (Figure 2). The transport through the Tokara Strait, which is a good indicator of the Kuroshio transport through the East China Sea, was estimated as 26 Sv (=106 m3・s−1), a value comparable to observed transports of 21.9 Sv [53] , 23.4 Sv [54] , 25 Sv [55] , and 28 Sv [56] . At the eastern mouth of the strait, this transport joins with the northward transport of 18 Sv associated with the Ryukyu Current. The

Figure 2. Stream function of the simulated 20-year mean transport from the sea surface to a depth of 1000 m or the sea bottom. Major transports referred to in the main text are emphasized by colored arrows with numbers showing the associated volume transport (Sv = 106 m3・s−1): green arrow, through the Tokara Strait; small red arrow, Ryukyu Current; large red arrow, Kuroshio; orange arrow, Kuroshio recirculation; gray arrow, Oyashio.

maximum velocity of the simulated northward flow was greater than 20 cm・s?1 and appeared at depths of 500 to 1000 m along the slope of the Ryukyu trench east of Ryukyu Islands (not shown), similar to observed results reported by [57] . In addition, the Kuroshio recirculation (31 Sv) joins the eastward transport of the Kuroshio off the southwestern coast of Japan. As a result, the maximum and net eastward transports of the Kuroshio were estimated as 76 Sv and 45 Sv, respectively. Observations by the ASUKA project [58] (Figure 1(a), ASUKA line) show an annual mean net Kuroshio transport off the southwestern coast of Japan of 42 Sv, composed of an eastward transport of 57 Sv and recirculation of 15 Sv. The reanalysis data thus tended to exhibit strong bias toward the recirculation transport. A similar but smaller bias is also present in the z-coordinate MOVE forecast system [4] . The simulated Oyashio (gray arrow) flows along the Kuril and Japan Islands. The southward transport of the Oyashio of 4 Sv was integrated from the Hokkaido coast to its maximum value along the OICE observation line [59] (Figure 1(a)) and represents the net Oyashio transport. This value of 4 Sv is comparable to but somewhat larger than the 10-year mean (1993-2002) of the net Oyashio transport of about 2.76 Sv that was estimated by combining along-track altimetry data with CTD-derived geostrophic transport at the reference level of 1000 dbar along the OICE line [59] . The discrepancy between the two values is likely attributable to the zero-velocity assumption at 1000 dbar used in estimating the observation-based geostrophic transport.

We next compared the 20-year mean current velocity at the sea surface based on the FRA-ROMS reanalysis data (Figure 3(a)) with the 20-year mean absolute geostrophic velocity at the sea surface (Figure 3(b)), which was estimated by adding the altimetry-

Figure 3. Twenty-year means of current velocity at the sea surface in (a) reanalysis data and (b) altimetry-derived absolute geostrophic velocity data.

derived SLA data to the CNES-CLS09 mean dynamic topography (MDT) [60] .

The gridded 1/4-degree SLA data, which have a weekly resolution, and the MDT data were collected from the AVISO (Archiving, Validation and Interpolation of Satellite Oceanographic data group) [61] . Overall, the simulated velocity and absolute geostrophic velocity were very consistent, particularly along the paths of the Kuroshio and Kuroshio Extension, but there were some discrepancies in the western boundary current region of the subarctic gyre. The East Kamchatka Current, which flows along the eastern coast of the Kamchatka Peninsula, and the Oyashio, where it flows along the southeastern coast of Kuril, Four, and Hokkaido Islands, were more distinct in the simulation, probably owing to the coarse resolution of the altimetry data near the coasts.

Comparison of the standard deviations between the simulated current velocity and altimetry-derived geostrophic velocity at the sea surface (Figure 4) showed that the spatial patterns of current variability were very consistent. The simulated velocity used

Figure 4. Standard deviations of current velocity at the sea surface in (a) reanalysis data and (b) altimetry-derived geostrophic velocity data. The standard deviations, which were separately estimated from the east and north velocity components and then superposed, are depicted by the color scale. The vectors show the 20-year mean velocity at the sea surface (as in Figure 3).

in this comparison was based on the daily mean and subsampled at weekly intervals to adjust the time of the simulated data to that of the altimetry-derived data. The current variability along the Yellow Sea coast of China was smaller in the reanalysis data. However, this underestimation is probably only apparent, because tidal components were insufficiently removed from the AVISO altimetry data [62] . In addition, in the subarctic region north of 40˚N, the current variability was greater in the reanalysis data. This discrepancy suggests that variations with a timescale of less than a week (e.g., in wind- induced currents within the Ekman layer [63] , which are associated with ageostrophic components), contributed more significantly to the current velocity at the sea surface in the subarctic region.

Daily SST based on reanalysis data was compared with a daily gridded SST dataset with a horizontal resolution of 0.25˚ [64] , which was created by NOAA by using optimal interpolation to merge satellite, ship, and buoy observations (product name: AVHRR_OI). The comparison was performed for the area from 20˚N to 50˚N and 120˚E to 160˚E, but excluding the Okhotsk Sea. Root-mean-square differences (RMSDs), mean differences, and correlation coefficients between the two SST datasets were estimated at monthly intervals (Figure 5(a)). Correlation coefficients (not shown) were greater than

Figure 5. (a) Monthly time series of the RMSD (thick line) and mean difference of SST (thin line) between the reanalysis data and AVHRR_OI. The difference was estimated by subtracting AVHRR_OI values from the reanalysis data. The target area was limited to the area from 20˚N to 50˚N and 120˚E to 160˚E, excluding the Okhotsk Sea. Scatter plots of the two SST datasets are depicted as a frequency distribution (color scale (unit: number of data points for each class)) for (b) January and (c) July. Spatial maps of the RMSD of the SSTs during 20 years (contours; contour interval, 1˚C) are shown for (d) January and (e) July; RMSDs exceeding 0.5˚C are emphasized by color.

0.991 and generally high throughout the analysis period. The RMSDs fluctuated in the range 0.63 - 1.10 (Figure 5(a)), and the 20-year monthly mean RMSDs exhibited two peaks, one in winter and another in summer (not shown). Small negative biases also appeared in winter (Figure 5(a)), and the bias was remarkable for SSTs less than 20˚C (Figure 5(b)). Relatively large RMSDs exceeding 2˚C were distributed mainly in the Kuroshio-Oyashio transition zone (Figure 5(d)). These large RMSDs were partly because in the reanalysis data the SST fronts tended to be sharper than they were in the AVHRR_OI data, which were spatially smoothed by the optimal interpolation. In contrast, positive biases were apparent in summer (Figure 5(a)), reflecting the fact that summertime SSTs in the reanalysis data tended to be higher than those in the AVHRR_OI data (Figure 5(c)). Relatively large RMSDs, greater than 2˚C, were related to this positive bias and showed a scattered distribution along the coasts of the Yellow Sea and along the southern coasts of the Kuril and Four Islands northeast of Hokkaido (Figure 5(e)). These regions are characterized by strong tidal currents and mixing, which were not parameterized or simulated by FRA-ROMS.

The positions of the axes of the Kuroshio and Kuroshio Extension were specified from the reanalysis data by combining the methods of [65] and [66] . Used alone, the former method could not successfully specify the axis position of the Kuroshio Extension, where mesoscale eddies interact intermittently with the mainstream axis. Therefore, the Kuroshio axis position within the 123˚E - 140˚E interval was first determined by the former method. Then, following the method of Qiu and Chen [66] , sea surface heights within the 136˚E - 140˚E interval were averaged along the pre-determined axis position. The position of the sea surface height isoline with the mean along-stream value was extracted and used as the axis position from 136˚E to 160˚E. Finally, the Kuroshio axis positions in the 136˚E - 140˚E interval estimated by the two methods were averaged.

The 20-year mean axis position of the Kuroshio and its standard deviation were compared between the reanalysis data and observation-based QBOC (Quick Bulletin of Ocean Conditions) data digitized by the Marine Information Research Center (MIRC) (Figure 6(a)). Although some regional differences could be detected, the mean and standard deviation of the reanalysis data seemed to be very consistent with those of the observation-based data. The axis-position difference between the two datasets in the interval from 126.5˚E to 140.5˚E was defined following [67] [68] as the ratio of Q to L, where Q denotes the area between the two axes and L is the axis length in the MIRC observation data. The difference fluctuated temporally but was relatively stable overall throughout the period analyzed (Figure 6(b)). The mean difference was 20.2 km, which is almost the same as the width of two grid cells in the 1/10-degree model, the minimum was 7.5, the maximum was 50.9 km, and the standard deviation was 5.5 km.

Year-to-year variations in the simulated axis position of the Kuroshio Extension from 1993 to 2012 are illustrated in Figure 7. In both the simulation and observation [66] , the Kuroshio Extension exhibited stable and unstable modes. We compared the path length and the meridional position of the axes in the simulation with the observation-based results of [69] (Figure 8). The simulated and observed values were similar, suggesting that the bimodal variation of the Kuroshio Extension was successfully simulated in the reanalysis data. The path of the Kuroshio Extension seemed to show a

Figure 6. (a) 20-year mean (thick curves) and meridional standard deviation (thin vertical bars) of the Kuroshio axis position based on the reanalysis data (red) and observation-based MIRC data (blue). (b) Time series of the difference in the Kuroshio axis position in the 126.5˚E - 140.5˚E interval between the reanalysis and MIRC (QBOC) data.

higher frequency of variation in the reanalysis data compared with the observation; this difference can be attributed to the difference in time resolution between the weekly data utilized by [69] and the daily reanalysis data.

We also examined the reproducibility in the reanalysis data of the water-mass distribution in the Kuroshio-Oyashio transition zone, where the Kuroshio and Oyashio waters are distributed and intermixed in a complicated way [70] . The Kuroshio-derived waters are mainly supplied to the transition zone by the Kuroshio Extension and the Tsugaru Warm Current (Figure 1(a)). Oceanographic conditions in the transition zone are conventionally described by using spatial patterns of isotherms of an index temperature at depths of 100 or 200 m [70] . We used this method to examine the reproducibility of the southernmost latitudes of the first and second Oyashio intrusions and of the easternmost longitude of the Tsugaru Warm Current (Figure 1(a)). The southernmost latitudes of the Oyashio intrusions can be used to represent the Oyashio water distribution (Figure 1(a)) and are closely related to the formation of fishing grounds for pelagic fishes (e.g., [71] [72] [73] ). The easternmost longitude of the Tsugaru Warm Current

Figure 7. Year-to-year variability of the simulated Kuroshio Extension axis from 1993 to 2012.

is useful not only for representing the water-mass distribution but also for assessing the capability of a mesoscale model, because from summer to autumn the Tsugaru Warm Current forms a clockwise gyre with a maximum diameter of 100 km [74] , which is about the minimum magnitude of variation that mesoscale models with a grid size of about 10 km are generally able to reproduce.

The southernmost latitudes and easternmost longitude were determined by using isotherms of 5˚C and greater than 6˚C, respectively, at a depth of 100 m, as explained in the caption of Figure 9. The southernmost latitude of the first Oyashio intrusion was very consistent between the reanalysis and observed values (Figure 9(a)). Both seasonal and interannual variations seemed to be reproduced well in the reanalysis data. The southernmost latitude of the second intrusion, however, was not precisely reproduced in the reanalysis data (Figure 9(b)), mainly because the observation data, which were obtained by an in situ ship survey, were frequently limited to coastal waters near Honshu and Hokkaido Islands; thus, the estimation accuracy of the observed temperature map used to determine the southernmost latitude of the second Oyashio intrusion was probably temporally inhomogeneous and dependent on the observed temperature distribution.

The easternmost longitude of the Tsugaru Warm Current showed good agreement between the reanalysis and observation data (Figure 9(c)), though there were some slight discrepancies. In the reanalysis data, the easternmost longitude exhibited an

Figure 8. (a) Path length (141˚E - 153˚E) and (b) mean meridional position (141˚E - 158˚E) of the Kuroshio Extension (KE). Blue lines indicate the estimates of [69] , and red lines show the reanalysis data.

eastward bias from winter to spring, when the gyre shrank, and a westward bias from summer to autumn, when the gyre expanded to a maximum, indicating that FRA- ROMS tended to underestimate the seasonal variability of the Tsugaru gyre. Because of the 1/10 degree resolution of the model, it would be difficult to improve the reproducibility of the expansion of the Tsugaru Warm Current in the reanalysis data. Dynamical downscaling should resolve this problem, as discussed in Section 4.

3.2. Analyses from a Lagrangian Viewpoint

Particle-tracking methods, including use of IBMs and super-IBMs, are among the most powerful techniques used in applications of reanalysis data to fisheries science (e.g., [75] - [81] ). However, particle-tracking errors in a Lagrangian reference frame have hardly been evaluated, mainly because no well-established error evaluation method has been shared among individual modelers developing ocean forecast systems. Hence, here we propose a method, which we call the best performance particle (BPP) method, in which particle-tracking errors of mesoscale models are estimated by comparing the positions of numerically tracked particles with those of actual surface drifters.

Figure 9. Comparisons of the southernmost latitudes of the (a) first and (b) second Oyashio intrusion and (c) the easternmost longitude of the Tsugaru Warm Current water (TWCW) between reanalysis (red lines) and observation data (blue lines). These indices were determined from about 10-day (reanalysis data) and monthly (observation data interpolated by Gaussian filters) mean temperature maps at a depth of 100 m. The Oyashio domain was initially specified on a 100-m temperature map as a cold water area surrounded by 5˚C isotherms and distributed continuously off the eastern coast of Hokkaido. The boundary of the Oyashio domain typically exhibits tongue-shaped southward intrusions. The intrusions closest and second closest to the Japanese coast are defined as the first and second Oyashio intrusion, respectively. The easternmost longitude of the Tsugaru Warm Current was determined as the easternmost position of isotherms greater than 6˚C at a depth of 100 m that pass through the Tsugaru Strait (Figure 1) off southern Hokkaido into the North Pacific.

For validation data, we used Surface Velocity Program (SVP) drifter trajectory data that had passed a sophisticated quality check performed by NOAA’S Atlantic Oceanographic and Meteorological Laboratory. Only data for drifters with a drogue were utilized. SVP drifters are designed to measure mixed layer currents in the upper ocean; thus, the center position of the drogue is set at 15 m below the sea surface, and its vertical length is about 6.44 m. Daily mean positions of the drifters were estimated from original 6-hourly data. Data from drifters with a tracking period (=N) longer than 60 days and without any missing data were selected. The data from each drifter were divided into subsets, each consisting of 60 days of daily observations. Thus, the number of the subsets for a single drifter was N-59. A total of 57,608 subsets were available for comparison with particle-tracking results in the Kuroshio-Oyashio region (Figure 10(a)).

A particle-tracking experiment was performed with FRA-ROMS for 60 days, and the results were compared with the daily drifter positions in each data subset. The 60-day period was chosen because for the first few months after spawning, the swimming capacity of eggs, larvae, and juveniles of pelagic fish is generally negligible (e.g., [21] [23] [78] [82] ). A particle-tracking model (LTRANS) specialized for ROMS output [83] was modified and applied in this study. Particles were transported by the simulated current velocity at a fixed depth without any diffusion. A group of particles was initially centered around the actual drifter position on the beginning day of each data subset. The particles were evenly spaced at intervals of 0.0267˚ (~3 km) of latitude and longitude within ±0.4˚ (~40 km) of latitude and longitude of the drifter position and at depths of 5, 10, 15, and 20 m. A total of 3844 particles (961 at each depth were tracked) for each data subset. RMSEs were estimated from difference in horizontal position between the drifter and the individual particles at daily intervals for 60 days. Finally, the single particle with the minimum RMSE value, that is, the BPP, was selected and its statistical

Figure 10. (a) Number of 60-day drifter data subsets in each grid cell (resolution, 1/2 degree ´ 1/2 degree). The data subsets were counted in the grid cell in which the initial drifter position of each data subset was located. (b) Drifter trajectory lengths integrated over 60 days. Lengths greater than 2000 km are shown by red squares, and those less than 1000 km are shown by blue squares. (c) RMSEs estimated from the position difference between the drifter and the BPP. Red squares indicate RMSEs greater than 200 km, and blue squares indicate those less than 100 km. (d) Ratios of RMSE (c) to integrated trajectory length (b). Red squares indicate ratios greater than 0.2, and blue squares indicate those less than 0.1. In (b) to (d), all values are averaged over the grid cell where the initial drifter position of each data subset was located.

properties were analyzed. Multiple particles were tracked, and the BPP among them was selected both to reduce the effect of potential bias in the reanalysis data and because some conditions that would affect drifter motion were not appropriately considered in the particle-tracking experiment, as explained in greater detail in Appendix 1.

If, under idealized conditions, the simulated current velocity includes a constant unidirectional error of 2 cm/s (3 cm/s), the position error and RMSE after particle tracking for 60 days would be 104 and 60 km (156 and 90 km), respectively. Therefore, we considered RMSEs of less than 100 km to indicate relatively accurate reproducibility of the reanalysis data. In fact, as a criterion for evaluating tracking errors, the RMSE may be rather strict, because it was often relatively large when the trajectories of the particle and drifter after tracking for 60 days were similar but their daily positions were different (not shown).

The RMSE of the BPP was weakly correlated with the trajectory length of the drifter integrated over 60 days (r = 0.45; Figure 11), the mean RMSE and mean integrated trajectory length were 149 km and 1430 km, respectively, and their ratio was about 0.1. We examined the spatial distribution of integrated trajectory lengths (Figure 10(b)), RMSEs (Figure 10(c)), and their ratios (Figure 10(d)) by mapping their values averaged over each grid cell (resolution, 1/2 degree ´ 1/2 degree). All three values were mapped in the grid cell where the initial drifter position of each data subset was located. Around the Kuroshio and Kuroshio Extension, the integrated trajectory length exceeded 2000 km, whereas east of 140˚E in the 20˚N - 30˚N latitude range and east of 150˚E in the 40˚N - 50˚N latitude range, where the current velocities were small, the integrated trajectory lengths were less than 1000 km (Figure 10(b)). The distributions of larger and smaller RMSEs (Figure 10(c)) roughly corresponded to areas where the integrated trajectories were longer and shorter, respectively (Figure 10(b)). Hence,

Figure 11. Scatter plots of the RMSE between the drifter and BPP positions versus the trajectory distance of the drifters integrated over 60 days. Scatter plots were displayed as frequency (unit of percentage) in each class (i.e., data-point number in each class divided by total data number).

regardless of the initial particle position, the ratios were mostly less than 0.1 across the Kuroshio-Oyashio region, although ratios exceeding 0.2 were scattered in the region (Figure 10(d)). This result suggests that the particle-tracking error of the BPP can be roughly estimated by multiplying the integrated trajectory length of the drifter at the particle’s initial position by 0.1.

4. Ongoing Studies and Near-Future Plans

The FRA has employed FRA-ROMS reanalysis and prediction products in multiple projects, mainly ones sponsored by the Fisheries Agency. Reanalysis and prediction data are routinely utilized for forecasting fishing and oceanographic conditions in the Kuroshio-Oyashio region (http://abchan.fra.go.jp/index2.html). Predictability of an operational version of FRA-ROMS will be demonstrated in near future work. A variety of IBMs based on the reanalysis data produced by the Fisheries Agency’s research enterprises are also being developed to clarify processes leading to year-to-year variations in recruitment and stocks, and to improve stock assessment accuracy, of important commercial fish such as Japanese sardine (Sardinops melanostictus) [84] , Japanese chub mackerel (Scomber japonicus) [84] , Japanese common squid (Todarodes pacificus) [85] , Japanese walleye pollock (Theragra chalcogramma) [84] , and chum salmon (Oncorhynchus keta) [86] . Moreover, sometimes in summer, the giant jellyfish (Nemopilema nomurai) seriously affects the Japanese fishery, so the process and pathway by which this species is transported from the East China Sea to around Japan are also being investigated and predicted by particle-tracking experiments conducted with FRA- ROMS.

Lower-trophic-level ecosystem models lacking the ROMS default options have been tuned and coupled with FRA-ROMS. For example, a simple NPZD (nutrient, phytoplankton, zooplankton, detritus) ecosystem model [87] has been adapted for examining the nutrient supply process to the shelf-slope region facing the Kuroshio south of Japan [88] . A more complicated ecosystem model, eNEMURO [89] , has also been coupled with FRA-ROMS to project future food conditions for pelagic fish and to evaluate impacts of global warming on Japanese fisheries (Okunishi and Setou, personal comm.).

Many downscaling models that can be linked to the reanalysis and prediction data of FRA-ROMS have been developed for simulating coastal waters around Japan. A resolution of 1/50 degree has been selected for studying coastal-offshore interactions in the shelf-slope region [29] [30] [75] . The simulated outputs have been applied, for example, in scientific projects aimed at promoting recovery from the unparalleled damages to fisheries caused by the 11 March 2011 tsunami. For use in small, semi-enclosed bays such as Mikawa Bay and the Yatsushiro Sea in Japan, super-high-resolution ocean models with a grid size of less than 300 m have been developed [75] [90] [91] [92] to simulate the transport of jellyfish, including vertical movements, and the occurrence of red tides by using FRA-ROMS in combination with non-passive particle-tracking models. Hasegawa et al. (personal comm.) have also developed a 1/90-degree wave-ocean-sediment coupled model that uses FRA-ROMS reanalysis data to improve understanding of processes affecting the transport of the radioactive materials that were accidentally leaked from the Fukushima Nuclear Power station after the 2011 mega-earthquake and attached to sediments on the sea bottom.

To maintain and enhance the reproducibility of FRA-ROMS, an oceanographic data acquisition system has been developed by oceanographers at the FRA. TS profiles obtained by regular monitoring by Japanese prefectural fisheries institutes have been since 2007 distributed to the Global Temperature?Salinity Profile Program (GTSPP) [93] . Recently, new types of observation instruments, such as glider and underway CTDs, have been adapted to measure fishery and oceanographic conditions. These instruments have made it possible to obtain high-resolution measurements able to resolve mesoscale and submesoscale variability. The FRA has developed a system to transfer these high-resolution measurements to FRA-ROMS in the TESAC format [94] . Such efforts cannot only improve reproducibility of the reanalysis data (i.e., initial state variables for forecasts) but also enhance predictability of the operational version of FRA- ROMS.

FRA-ROMS is one of the few operational ocean forecast systems in the world that is specialized and optimized for fisheries science. In particular, FRA-ROMS embeds a versatile subsystem with multiple functions so that model outputs can be efficiently applied to fisheries science (Okunishi, personal comm.). Using this subsystem, even a scientist without any special knowledge of physical oceanography and numerical techniques can easily extract, analyze, and illustrate necessary elements from the model outputs. Moreover, a function to conduct particle-tracking experiments is being developed. This subsystem is available to not only all FRA scientists but also registered scientists affiliated with research institutes of the Japanese prefectural fisheries. In the near future, several more useful functions will be developed and the user-friendliness of the subsystem will be enhanced.

Plans exist to upgrade FRA-ROMS. One modification will extend the model domains. The eastern and southern boundaries of the 1/10-degree model will be extended to the North American continent and to 5˚S - 20˚S, respectively. The 1/2-degree model domain will also be optimally expanded. Further, the 1/2- and 1/10-degree models will be connected by applying a two-way nesting technique. In addition, several effects (e.g., tidal mixing, sea ice, and freshwater discharge from land) not considered in the current version, will be incorporated into the model. Moreover, a high-resolution (1/50 degree) coastal model is being developed. This model, which is an extended version of the model used by Kuroda et al. [29] [30] [76] , is designed to cover the entire Japanese coastline, not only in the Kuroshio-Oyashio region but also along the Japan Sea. The use of the FRA-ROMS platform and its associated models, combined with a more user- friendly subsystem, is expected to lead to dramatic progress in Japanese fisheries science in the near future.


We express our deepest gratitude to Dr. M. Kamachi, Dr. N. Usui, and Dr. Y. Fujii of the Meteorological Research Institute (MRI) for technical support and for warm- hearted encouragement. FRA-ROMS was developed and has been implemented since 2008 on the vector and cluster supercomputing systems operated by the Agriculture, Forestry and Fisheries Research Information Technology Center (AFFRIT). FRA-ROMS was developed by the physical oceanographers of the FRA with partial support from the Fisheries Agency, but the contents of this study do not necessarily reflect the views of the Fisheries Agency. This work was also financially supported by Japan Society for the Promotion of Science KAKENHI Grant Number 23740365, a Kofu-kin research project by the FRA, and by the “The Study of Kuroshio Ecosystem Dynamics for Sustainable Fisheries (SKED)” project of the Ministry of Education, Culture, Sports, Science and Technology.

Appendix 1: Why Only the BPP Was Selected from Multiple Tracked Particles to Compare with the Drifter?

Here, we explain why multiple particles were tracked and only the BPP was selected. A group of particles was initially centered about the position of an actual drifter within a range of ±0.4˚ latitude and longitude. As explained in Section 3.1, the mean difference in the position of the Kuroshio axis between the reanalysis and the observation-based data was 20.2 km, and the standard deviation (=σ) was 5.5 km, signifying that the position of the Kuroshio jet was simulated within an error range of ±3σ (3.7 to 36.7 km). By setting the initial positions of the particles within ±0.4˚ (~40 km) of the drifter position and at several depths (5, 10, 15, and 20 m) around the depth of the drifter’s drogue (about 11 - 18 m), the effects of this initial error on the particle-tracking results should be rendered negligible.

Moreover, the mesoscale variations included in the reanalysis data were not by themselves enough for the particle transport to realistically reproduce the drifter movement. Some elements affecting drifter movements in the actual ocean were not considered by our particle-tracking experiment: anisotropic lateral diffusion due to submesoscale variations associated with sub-grid scale variations [29] [30] [76] ; tidal residual currents [95] ; Stokes drift due to wind waves [96] ; surface-wave breaking [97] ; effects of local wind stress that was not resolved in the JRA25 data [98] ; and leeway and slip [99] . These facts suggest that a tracking experiment using a single particle with its initial position fixed at the observed drifter position could not appropriately evaluate particle- tracking errors. Consequently, we tracked multiple particles and then picked the BPP, the one with the minimum RMSE value, and analyzed the trajectory of the BPP as the best candidate among all of the particles released for evaluation of Lagrangian trajectory errors.

Appendix 2: Abbreviations

There are many abbreviations in this article. Frequently-used abbreviations are summarized as follows.

AVISO: Archiving, Validation and Interpolation of Satellite Oceanographic data group

BPP: The Best Performance Particle

FRA: The Japan Fisheries Research and Education Agency

IBM: Individual-Based Model

ROMS: Regional Ocean Modeling System

SLA: Sea Level Anomaly

SST: Sea Surface Temperature

TS: Temperature and Salinity

3D-Var: Three-Dimensional Variational

Cite this paper
Kuroda, H. , Setou, T. , Kakehi, S. , Ito, S. , Taneda, T. , Azumaya, T. , Inagake, D. , Hiroe, Y. , Morinaga, K. , Okazaki, M. , Yokota, T. , Okunishi, T. , Aoki, K. , Shimizu, Y. , Hasegawa, D. , Watanabe, T. , (2017) Recent Advances in Japanese Fisheries Science in the Kuroshio-Oyashio Region through Development of the FRA-ROMS Ocean Forecast System: Overview of the Reproducibility of Reanalysis Products. Open Journal of Marine Science, 7, 62-90. doi: 10.4236/ojms.2017.71006.
[1]   Metzger, E.J., Smedstad, O.M., Thoppil, P.G., Hurlburt, H.E., Cummings, J.A., Wallcraft, A.J., Zamudio, L., Franklin, D.S., Posey, P.G., Phelps, M.W., Hogan, P.J., Bub, F.L. and DeHaan, C.J. (2014) US Navy Operational Global Ocean and Arctic Ice Prediction Systems. Oceanography, 27, 32-43.

[2]   Lellouche, J.-M., Le Galloudec, O., Drévillon, M., Régnier, C., Greiner, E., Garric, G., Ferry, N., Desportes, C., Testut, C.-E., Bricaud, C., Bourdallé-Badie, R., Tranchant, B., Benkiran, M., Drillet, Y., Daudin, A. and De Nicola, C. (2014) Evaluation of Global Monitoring and Forecasting Systems at Mercator Ocean. Ocean Science, 9, 57-81.

[3]   Storkey, D., Blockley, E.W., Furner, R., Guiavarc’h, C., Lea, D., Martin, M.J., Barciela, R.M., Hines, A., Hyder, P. and Siddorn, J.R. (2010) Forecasting the Ocean State Using NEMO: The New FOAM System. Journal of Operational Oceanography, 3, 3-15.

[4]   Usui, N., Ishizaki, S., Fujii, Y., Tsujino, H., Yasuda, T. and Kamachi, M. (2006) Meteorological Research Institute Multivariate Ocean Variational Estimation (MOVE) System: Some Early Results. Advances in Space Research, 37, 806-822.

[5]   Miyazawa, Y., Zhang, R., Guo, X., Tamura, H., Ambe, D., Lee, J.-S., Okuno, A., Yoshinari, H., Setou, T. and Komatsu, K. (2009) Water Mass Variability in the Western North Pacific Detected in a 15-Year Eddy Resolving Ocean Reanalysis. Journal of Oceanography, 65, 737-756.

[6]   Bell, M.J., Lefebvre, M., Le Traon, P.-Y., Smith, N. and Wilmer-Becker, K. (2009) GODAE: The Global Ocean Data Assimilation Experiment. Oceanography, 22, 14-21.

[7]   Bell, M.J., Schiller, A., Le Traon, P.-Y., Smith, N.R., Dombrowsky, E. and Wilmer-Becker, K. (2015) An Introduction to GODAE Ocean View. Journal of Operational Oceanography, 8, s2-s11.

[8]   Moore, A.M., Arango, H.G., Broquet, G., Powell, B.S., Weaver, A.T. and Zavala-Garay, J. (2011) The Regional Ocean Modeling System (ROMS) 4-Dimensional Variational Data Assimilation Systems: Part I—System Overview and Formulation. Progress in Oceanography, 91, 34-49.

[9]   Forget, G., Campin, J.-M., Heimbach, P., Hill, C.N., Ponte, R.M. and Wunsch, C. (2015) ECCO Version 4: An Integrated Framework for Non-Linear Inverse Modeling and Global Ocean State Estimation. Geoscientific Model Development, 8, 3071-3104.

[10]   Graff, J. (2008) The Role of Operational Ocean Forecasting in E-Navigation. International Journal on Marine Navigation and Safety of Sea Transport, 2, 259-262.

[11]   Lam, F.-P.A., Haley Jr., P.J., Janmaat, J., Lermusiaux, P.F.J., Leslie, W.G., Schouten, M.W., te Raa, L.A. and Rixen, M. (2009) At-Sea Real-Time Coupled Four-Dimensional Oceanographic and Acoustic Forecasts During Battlespace Preparation 2007. Journal of Marine Systems, 78, S306-S320.

[12]   Potemra, J.T. (2012) Numerical Modeling with Application to Tracking Marine Debris. Marine Pollution Bulletin, 65, 42-50.

[13]   Lebrenton, L.C.-M. and Borrero, J.C. (2013) Modeling the Transport and Accumulation Floating Debris Generated by the 11 March 2011 Tohoku Tsunami. Marine Pollution Bulletin, 66, 53-58.

[14]   Copeland, G.J.M. and Thiam-Yew, W. (2004) Current Data Assimilation Modelling for Oil Spill Contingency Planning. Environmental Modelling and Software, 21, 142-155.

[15]   Liu, Y., Macfadyen, A., Ji, Z.-G. and Weisberg, R.H. (2013) Monitoring and Modeling the Deepwater Horizon Oil Spill: A Record-Breaking Enterprise. The American Geophysical Union as Part of the Geophysical Monograph Series, Vol. 195.

[16]   Miyazawa, Y., Masumoto, Y., Varlamov, S.M., Miyama, T., Takigawa, M., Honda, M. and Saino, T. (2013) Inverse Estimation of Source Parameters of Oceanic Radioactivity Dispersion Models Associated with the Fukushima Accident. Biogeosciences, 10, 2349-2363.

[17]   Tsumune, D., Tsubono, T., Aoyama, M., Uematsu, M., Misumi, K., Maeda, Y., Yoshida, Y. and Hayami, H. (2013) One-year, Regional-Scale Simulation of 137Cs Radioactivity in the Ocean Following the Fukushima Dai-ichi Nuclear Power Plant Accident. Biogeosciences, 10, 5601-5617.

[18]   Igarashi, H., Ichii, T., Sakai, M., Ishikawa, Y., Toyoda, T., Masuda, S., Sugiura, N., Mahapatra, K. and Awaji, T. (2015) Possible Link between Interannual Variation of Neon Flying Squid (Ommastrephes bartramii) Abundance in the North Pacific and the Climate Phase Shift in 1998/1999. Progress in Oceanography, in Press.

[19]   Igarashi, H., Awaji, T., Toyoda, T., Masuda, S., Sugiura, N., Sasaki, Y., Hiyoshi, Y., Sakai, M., Ichii, T. and Ishikawa, Y. (2010) An Optimal Synthesis of Observations and Models by Data Assimilation: Applications to Climate Analysis and Fishery Assessment. IEEE Earthzine, 14 September 2010.

[20]   Chang, Y.-L., Sheng, J., Ohashi, K., Béguer-Pon, M. and Miyazawa, Y. (2015) Impacts of Interannual Ocean Circulation Variability on Japanese Eel Larval Migration in the Western North Pacific Ocean. PLoS ONE, 10, e0144423.

[21]   Kasai, A., Komatsu, K., Sassa, C. and Konishi, Y. (2008) Transport and Survival Processes of Eggs and Larvae of Jack Mackerel Trachurus japonicus in the East China Sea. Fisheries Science, 74, 8-18.

[22]   Itoh, S., Saruwatari, T., Nishikawa, H., Yasuda, I., Komatsu, K., Tsuda, A., Setou, T. and Shimizu, M. (2011) Environmental Variability and Growth Histories of Larval Japanese Sardine (Sardinops melanostictus) and Japanese Anchovy (Engraulis japonicus) near the Frontal Area of the Kuroshio. Fisheries Oceanography, 20, 114-124.

[23]   Nishikawa, H., Yasuda, I., Komatsu, K., Sasaki, H., Sasai, Y., Setou T. and Shimizu, M. (2013) Winter Mixed Layer Depth and Spring Bloom along the Kuroshio Front: Implications for the Japanese Sardine Stock. Marine Ecology Progress Series, 487, 217-229.

[24]   Weisberg, R.H., Barth, A., Alvera-Azcárate, A. and Zheng, L. (2009) A Coordinated Coastal Ocean Observing and Modeling System for the West Florida Continental Shelf. Harmful Algae, 8, 585-597.

[25]   Alabia, I.D., Saitoh, S., Mugo, R., Igarashi, H., Ishikawa, Y., Usui, N. Kamachi, M., Awaji, T. and Seito, M. (2015) Identifying Pelagic Habitat Hotspots of Neon Flying Squid in the Temperate Waters of the Central North Pacific. PLoS ONE, 10, e0142885.

[26]   Zhang, X., Saitoh, S., Hirawake, T., Nakada, S., Koyamada, K., Awaji, T., Ishikawa, Y. and Igarashi, H. (2013) An Attempt of Dissemination of Potential Fishing Zones Prediction Map of Japanese Common Squid in the Coastal Water, Southwestern Hokkaido, Japan. Proceedings of the Asia-Pacific Advanced Network, 36, 132-141.

[27]   Komatsu, K., Matsukawa, Y., Nakata, K., Ichikawa, T. and Sasaki, K. (2007) Effect of Advective Processes on Planktonic Distributions in the Kuroshio Region Using a 3-D Lower Trophic Model and a Data Assimilative OGCM. Ecological Modelling, 202, 105-119.

[28]   Toyoda, T., Awaji, T., Masuda, S., Sugiura, N., Igarashi, H., Sasaki, Y., Hiyoshi, Y., Ishikawa, Y., Saitoh, S., Yoon, S., In, T. and Kishi, M.J. (2013) Improved State Estimations of Lower Trophic Ecosystems in the Global Ocean Based on a Green’s Function Approach. Progress in Oceanography, 119, 90-107.

[29]   Kuroda, H., Setou, T., Aoki, K., Takahashi, D., Shimizu, M. and Watanabe, T. (2013) A Numerical Study of the Kuroshio-Induced Circulation in Tosa Bay, off the Southern Coast of Japan. Continental Shelf Research, 53, 50-62.

[30]   Kuroda, H., Hirota, Y., Setou, T., Aoki, K., Takahashi, D. and Watanabe, T. (2014) Properties of Winter Mixed Layer Variability on the Shelf-Slope Region Facing the Kuroshio— Study of Tosa Bay, Southern Japan. Ocean Dynamics, 64, 47-60.

[31]   Watanabe, Y. and Wada, T. (1998) Stock Fluctuations and Ecological Changes of the Japanese Sardine. Kouseisha-Kouseikaku, Tokyo. (In Japanese)

[32]   Aoki, I. and Miyashita, K. (2000) Dispersal of Larvae and Juveniles of Japanese Anchovy Engraulis japonicus in the Kuroshio Extension and Kuroshio-Oyashio Transition Regions, Western North Pacific Ocean. Fisheries Research, 49,155-164.

[33]   Tian, Y., Akamine, T. and Suda, M. (2003) Variations in the Abundance of Pacific Saury (Cololabis saira) from the Northwestern Pacific in Relation to Oceanic-Climate Changes. Fisheries Research, 60, 439-454.

[34]   Kamimura, Y., Takahashi, M., Yamashita, N., Watanabe, C., and Kawabata, A. (2015) Larval and Juvenile Growth of Chub Mackerel Scomber japonicus in Relation to Recruitment in the Western North Pacific. Fisheries Science, 81, 505-513.

[35]   Kidokoro, H., Goto, T., Nagasawa, T., Nishida, H., Akamine, T. and Sakurai, Y. (2010) Impact of a Climate Regime Shift on the Migration of Japanese Common Squid (Todarodes pacificus) in the Sea of Japan. ICES Journal of Marine Science, 67, 1314-1322.

[36]   Sakurai, Y. (2007) An Overview of the Oyashio Ecosystem. Deep Sea Research Part II: Topical Studies in Oceanography, 54, 2526-2542.

[37]   Yatsu, A., Chiba, S., Yamanaka, Y., Ito, S., Shimizu, Y., Kaeriyama, M. and Watanabe, Y. (2013) Climate Forcing and the Kuroshio/Oyashio Ecosystem. ICES Journal of Marine Science, 70, 922-933.

[38]   Okunishi, T., Yamanaka, Y. and Ito, S. (2009) A Simulation Model for Japanese Sardine (Sardinops melanostictus) Migrations in the Western North Pacific. Ecological Modelling, 220, 462-479.

[39]   Okunishi, T., Ito, S., Ambe, D., Takasuka, A., Kameda, T., Tadokoro, K., Setou, T., Komatsu, K., Kawabata, A., Kubota, H., Ichikawa, T., Sugisaki, H., Hashioka, T., Yamanaka, Y., Yoshie, N. and Watanabe, T. (2012) A Modeling Approach to Evaluate Growth and Movement for Recruitment Success of Japanese Sardine (Sardinops melanostictus) in the Western Pacific. Fisheries Oceanography, 21, 44-57.

[40]   Hirose, N., Takayama, K., Moon, J.-H., Watanabe, T. and Nishida, Y. (2013) Regional Data Assimilation System Extended to the East Asian Marginal Seas. Umi to Sora, 89, 43-51.

[41]   Haidvogel, D.B., Arango, H., Budgell, W.P., Cornuelle, B.D., Curchitser, E., Di Lorenzo, E., Fennel, K., Geyer, W.R., Hermann, A.J., Lanerolle, L., Levin, J., McWilliams, J.C., Miller, A.J., Moore, A.M., Powell, T.M., Shchepetkin, A.F., Sherwood, C.R., Signell, R.P., Warner, J.C. and Wilkin, J. (2008) Ocean Forecasting in Terrain-Following Coordinates: Formulation and Skill Assessment of the Regional Ocean Modeling System. Journal of Computational Physics, 227, 3595-3624.

[42]   Guo, X., Hukuda, H., Miyazawa, Y. and Yamagata, T. (2003) A Triply Nested Ocean Model for Simulating the Kuroshio―Roles of Horizontal Resolution on JEBAR. Journal of Physical Oceanography, 33, 146-169.

[43]   Song, Y. and Haidvogel, D. (1994) A Semi-Implicit Ocean Circulation Model Using a Generalized Topography-Following Coordinate System. Journal of Computational Physics, 115, 228-244.

[44]   Stephens, C., Antonov, J.I., Boyer, T.P., Conkright, M.E., Locarnini, R.A., O’Brien, T.D. and Garcia, H.E. (2002) World Ocean Atlas 2001, Vol. 1, Temperature [CD-ROM], NOAA Atlas NESDIS, Vol. 49, Edited by S. Levitus, U.S. Government Printing Office, Washington DC.

[45]   Boyer, T.P., Stephens, C., Antonov, J.I., Conkright, M.E., Locarnini, R.A., O’Brien, T.D. and Garcia, H.E. (2002) World Ocean Atlas 2001, Vol. 2, Salinity [CD-ROM], NOAA Atlas NESDIS, Vol. 50, Edited by S. Levitus, U.S. Government Printing Office, Washington DC.

[46]   Onogi, K., Tsutsui, J., Koide, H., Sakamoto, M., Kobayashi, S., Hatsushika, H., Matsumoto, T., Yamazaki, N., Kamahori, H., Takahashi, K., Kadokura, S., Wada, K., Kato, K., Oyama, R., Ose, T., Mannoji, N. and Taira, R. (2007) The JRA-25 Reanalysis. Journal of the Meteorological Society of Japan, 85, 369-432.

[47]   Fairall, C.W., Bradley, E.F., Rogers, D.P., Edson, J.B. and Young, G.S. (1996) Bulk Parameterization of Air-Sea Fluxes for Tropical Ocean-Global Atmosphere Coupled-Ocean Atmosphere Response Experiment, Journal of Geophysical Research, 101, 3747-3764.

[48]   Fujii, Y. and Kamachi, M. (2003) A Reconstruction of Observed Profiles in the Sea East of Japan Using Vertical Coupled Temperature-Salinity EOF Modes. Journal of Oceanography, 59, 173-186.

[49]   Fujii, Y. and Kamachi, M. (2003) A Nonlinear Preconditioned Quasi-Newton Method without Inversion of a First-Guess Covariance Matrix in Variational Analysis. Tellus, 55A, 450-454.

[50]   Fujii, Y. and Kamachi, M. (2003) Three-Dimensional Analysis of Temperature and Salinity in the Equatorial Pacific Using a Variational Method with Vertical Coupled Temperature-Salinity Empirical Orthogonal Function Modes. Journal of Geophysical Research, 108, No. C9.

[51]   Usui, N., Ishizaki, S., Fujii, Y. and Kamachi, M. (2011) Improving Strategies with Constraints Regarding Non-Gaussian Statistics in a Three-Dimensional Variational Assimilation Method. Journal of Oceanography, 67, 253-262.

[52]   Bloom, S.C., Takacs, L.L., Da Silva A.M. and Ledvina, D. (1996) Data Assimilation Using Incremental Analysis Updates. Monthly Weather Review, 124, 1256-1271.

[53]   Bingham, F.M. and Talley, L.D. (1991) Estimates of Kuroshio Transport Using an Inverse Technique. Deep Sea Research, 38, S21-S43.

[54]   Feng, M., Mitsudera, H. and Yoshikawa, Y. (2000) Structure and Variability of the Kuroshio Current in Tokara Strait. Journal of Physical Oceanography, 30, 2257-2276.

[55]   Rikiishi, K., Hashimoto, Y., Matsuda, H. and Michigami, M. (2003) Monitoring the Kuroshio in the Tokara Strait and Izu Island Region by Using Submarine Cables. Proceedings of 3rd International Workshop on Scientific Use of Submarine Cables and Related Technology, Tokyo, 25-27 June 2003, 95.

[56]   Nakano, T., Kaneko, I. and Takatsuki, Y. (1994) The Kuroshio Structure and Transport Estimated by the Inverse Method. Journal of Physical Oceanography, 24, 609-618.

[57]   Ichikawa, H., Nakamura, H., Nishina, A. and Higashi, M. (2004) Variability of Northeastward Current Southeast of Northern Ryukyu Islands. Journal of Oceanography, 60, 351-363.

[58]   Imawaki, S., Uchida, H., Ichikawa, H., Fukasawa, M., Umatani, S. and ASUKA Group (2001) Satellite Altimeter Monitoring the Kuroshio Transport South of Japan. Geophysical Research Letters, 28, 17-20.

[59]   Ito, S., Uehara, K., Miyao, T., Miyake, H., Yasuda, I., Watanabe, T. and Shimizu, Y. (2004) Characteristics of SSH Anomaly Based on TOPEX/POSEIDON Altimetry and in Situ Measured Velocity and Transport of Oyashio on OICE. Journal of Oceanography, 60, 425-437.

[60]   Rio, M.H., Guinehut, S. and Larnicol, G. (2011) New CNES-CLS09 Global Mean Dynamic Topography Computed from the Combination of GRACE Data, Altimetry, and in Situ Measurements. Journal of Geophysical Research, 116.

[61]   SSALTO/DUACS User Handbook (2013) (M)SLA and (M)ADT Near-Real Time and Delayed Time Products. CLS-DOS-NT-06-034 SALP-MU-P-EA-21065-CLS Issue: 3rev 6.

[62]   Morimoto, A. (2009) Evaluation of Tidal Error in Altimetry Data in the Asian Marginal Seas. Journal of Oceanography, 65, 477-485.

[63]   Rio, M.-H. and Hernandez, F. (2003) High-Frequency Response of Wind-Driven Currents Measured by Drifting Buoys and Altimetry over the World Ocean. Journal of Geophysical Research, 108.

[64]   National Climatic Data Center (2007) GHRSST Level 4 AVHRR_OI Global Blended Sea Surface Temperature Analysis. Ver. 1.0. PO.DAAC, CA, USA.

[65]   Ambe, D., Imawaki, S., Uchida, H. and Ichikawa, K. (2004) Estimating the Kuroshio Axis South of Japan Using Combination of Satellite Altimetry and Drifting Buoys. Journal of Oceanography, 60, 375-382.

[66]   Qiu, B. and Chen, S. (2005) Variability of the Kuroshio Extension Jet, Recirculation Gyre and Mesoscale Eddies on Decadal Time Scales. Journal of Physical Oceanography, 35, 2090-2103.

[67]   Kamachi, M., Kuragano, T., Sugimoto, S., Yoshita, K., Sakurai, T., Nakano, T., Usui, N. and Uboldi, F. (2004) Short-Range Prediction Experiments with Operational Data Assimilation System for the Kuroshio South of Japan. Journal of Oceanography, 60, 269-282.

[68]   Kamachi, M., Kuragano, T., Ichikawa, H., Nakamura, H., Nishina, A., Isobe, A., Ambe, D., Arai, M., Gohda, N., Sugimoto, S., Yoshita, K., Sakurai, T. and Uboldi, F. (2004) Operational Data Assimilation System for the Kuroshio South of Japan: Reanalysis and Validation. Journal of Oceanography, 60, 303-312.

[69]   Qiu, B. and Chen, S. (2010) Eddy-Mean Flow Interaction in the Decadally Modulating Kuroshio Extension System. Deep Sea Research Part II, 57, 1098-1110.

[70]   Kawai, H. (1972) Hydrography of the Kuroshio Extension. Kuroshio, Its Physical Aspects. Stommell, H. and Yoshida, K. Eds., University of Tokyo Press, Tokyo, 235-352.

[71]   Hotta, H. and Fukushima, S. (1970) Fishing Ground of Saury and Its Oceanographic Conditions in the Pacific Ocean off Northeastern Japan. Bulletin of the Japanese Society for the Science of Fish, 36, 207-215. (In Japanese)

[72]   Kubo, Y. (1954) A Ecology Study of Cololabis saira (Brevoort) in the Pacific Ocean-II. Studies on the Genital Gland. Bulletin of Fisheries Experimental Station Ibaraki, 87-97. (In Japanese)

[73]   Uda, M. (1936) Fishing Centre of “Sanma” Cololabis saira (Brevoort), Correlated with the Head of the Oyashiwo Cold Current. Bulletin of the Japanese Society for the Science of Fish, 5, 236-238. (In Japanese with English Abstract)

[74]   Conlon, D.M. (1982) On the Outflow of the Tsugaru Warm Current. La mer, 20, 60-64.

[75]   Aoki, K., Onitsuka, G., Shimizu, M., Kuroda, H., Matsuo, H., Kitadai, Y. (2015) Interregional Differences in Mortality of Aquacultured Yellowtail Seriola quinqueradiata in Related to a Chattonella Bloom in the Yatsushiro Sea, Japan, in 2010. Fisheries Science, 81, 525-532.

[76]   Kuroda, H., Takahashi, D., Mitsudera, H., Azumaya, T. and Setou, T. (2014) A Preliminary Study to Understand the Transport Process for the Eggs and Larvae of Japanese Pacific Walleye Pollock Theragra chalcogramma Using Particle-tracking Experiments Based on a High-Resolution Ocean Model. Fisheries Science, 80, 127-138.

[77]   Moon, J.-H., Pang, I.-C., Yang, J.-Y. and Yoon, W.D. (2010) Behavior of the Giant Jellyfish Nemopilema nomurai in the East China Sea and East/Japan Sea during the Summer of 2005: A Numerical Model Approach Using a Particle-Tracking Experiment. Journal of Marine Systems, 80, 101-114.

[78]   Oozeki, Y., Okunishi, T., Takasuka, A. and Ambe, D. (2015) Variability in Transport Processes of Pacific Saury Cololabis saira Larvae Leading to Their Broad Dispersal: Implications for Their Ecological Role in the Western North Pacific. Progress in Oceanography, 138, 448-458.

[79]   Petrik, C.M., Duffy-Anderson, J.T., Mueter, F., Hedstrom, K. and Curchitser, E. (2015) Biophysical Transport Model Suggests Climate Variability Determines Distribution of Walleye Pollock Early Life Stages in the Eastern Bering Sea through Effects on Spawning. Progress in Oceanography, 138, 459-474.

[80]   Rose K.A., Fiechter, J., Curchitser, E.N., Hedstrom, K., Bernal, M., Creekmor, S., Haynie, A., Ito, S., Lluch-Cota, S., Megrey, B.A., Edwards, C.A., Checkley, D., Koslow, T., McClatchie, S., Werner, F., MacCall, A. and Agostini, V. (2015) Demonstration of a Fully-Coupled End-to-End Model for Small Pelagic Fish Using Sardine and Anchovy in the California Current. Progress in Oceanography, 138, 348-380.

[81]   Wei, H., Deng, L. and Wang, Y. (2015) Giant Jellyfish Nemopilema nomurai Gathering in the Yellow Sea—A Numerical Study. Journal of Marine Systems, 144, 107-116.

[82]   Itoh, S., Yasuda, I., Nishikawa, H., Sasaki, H. and Sasai, Y. (2009) Transport and Environmental Temperature Variability of Eggs and Larvae of the Japanese Anchovy (Engraulis japonicus) and Japanese Sardine (Sardinops melanostictus) in the Western North Pacific Estimated via Numerical Particle-Tracking Experiments. Fisheries Oceanography, 18, 118-133.

[83]   North, E.W., Schlag, Z., Hood, R.R., Li, M., Zhong, L., Gross, T. and Kennedy, V.S. (2008) Vertical Swimming Behavior Influences the Dispersal of Simulated Oyster Larvae in a Coupled Particle-Tracking and Hydrodynamic Model of Chesapeake Bay. Marine Ecology Progress Series, 359, 99-115.

[84]   Fisheries Agency and Fisheries Research Agency (2016) Annual Report of Analyses to Understand Causes of Variations in Fishery Stockaround Japan. (In Japanese)

[85]   Kuroda, H., Yamashita, N., Kaga, T., Azumaya, T. and Kuroda, H. (2012): A Study of Interannual Variation in Passive Transport Process of Larvae of Japanese Common Squid, Todarodes pacificus, Reproduced by Numerical Particle-Tracking Experiments Based on FRA-ROMS. Annual Report of the Research Meeting on Japanese Common Squid, 18-19. (In Japanese)

[86]   Fisheries Agency (2015) Section 2: Examination of Factors Affecting Migration and Survival of Salmon Juveniles. Annual Report of Project to Recover the Stock of Pacific Salmon, 40-55. (In Japanese)
http://www.maff.go.jp/j/budget/yosan_kansi/sikkou /tokutei_keihi/seika_h26/ippan/pdf/ippan46_05.pdf

[87]   Oschlies, A. (2001) Model-Derived Estimates of New Production: New Results Point towards Lower Values. Deep Sea Research Part II, 48, 2173-2197.

[88]   Kuroda, H., Takasuka, A., Hirota, Y., Aoki, K. and Setou, T. (2015) Characteristics of Winter Mixed Layer Variability and Nutrient Supply Processes onto the Shelf-slope Region Facing the Kuroshio. Bulletin on Coastal Oceanography, 53, 3-9. (In Japanese with English Abstract)

[89]   Yoshie, N., Guo, X., Fujii, N. and Komorita, T. (2011) Ecosystem and Nutrient Dynamics in the Seto Inland Sea, Japan. Interdisciplinary Studies on Environmental Chemistry, Modeling and Analysis of Marine Environmental Problems, 5, 39-49.

[90]   Aoki, K., Shimizu, M., Kuroda, H., Toyokawa, M. and Yamada, S. (2012) Numerical Study on the Transport Process of Jellyfish Aurelia auritasensu lato in Mikawa Bay, Japan. Bulletin of the Japanese Society of Fisheries Oceanography, 1, 9-17.

[91]   Aoki, K., Onitsuka, G., Shimizu, M., Kuroda, H., Matsuyama, Y., Kimoto, K., Matsuo, H., Kitadai, Y., Sakurada, K., Nishi., H. and Tahara, Y. (2012) Factors Controlling the Spatio-Temporal Distribution of the 2009 Chattonella antiqua Bloom in the Yatsushiro Sea. Estuarine, Coastal and Shelf Science, 114, 148-155.

[92]   Aoki, K., Onitsuka, G., Shimizu, M., Kuroda, H., Matsuo, H., Kitadai, Y., Sakurada, K., Ando, H., Nishi, H. and Tahara, Y. (2014) Variability of Factors Driving Spatial and Temporal Dispersion in River Plume and Chattonella antiqua Bloom in the Yatsushiro Sea, Japan. Marine Pollution Bulletin, 81, 131-139.

[93]   Uehara, K., Ito, S., Miyazawa, Y., Tsunoda, T., Kakehi, S., Akiyama, S., Kusaka, A., Setou, T., Komatsu, K. and Kameda, T. (2006) Framework of a Real-Time Data Distribution System. Monthly Kaiyo, 38, 515-518. (In Japanese)

[94]   Shimizu, Y., Okunishi, T., Kakehi, S., Hasegawa, D., Wagawa, T., Igeta, Y., Honda, N., Setou, T., Kuroda, H. and Ito, S. (2016) Historical Document of Operating Ocean Gliders in Fisheries Research Agency and Prospect of Their Usage for the Future. Journal of Fisheries Technology (Submitted). (In Japanese with English Abstract)

[95]   Yanagi, T. and Inoue, K. (1994) Tide and Tidal Residual Current in the Yellow/East China Sea. La mer, 32, 153-165.

[96]   Curcic, M., Chen, S.S. and Ozgokmen, T.M. (2016) Hurricane-induced Ocean Waves and Stokes Drift and Their Impacts on Surface Transport and Dispersion in the Gulf of Mexico. Geophysical Research Letter, 43, 2773-2781.

[97]   Carniel, S., Warner, J.C., Chiggiato, J. and Sclavo, M. (2009) Investigating the Impact of Surface Wave Breaking on Modeling the Trajectories of Drifters in the Northern Adriatic Sea During a Wind-Storm Event. Ocean Modelling, 30, 225-239.

[98]   Edwards, K.P., Werner, F.E. and Blanton, B.O. (2006) Comparison of Observed and Modeled Drifter Trajectories in Coastal Regions: An Improvement through Adjustments for Observed Drifter Slip and Errors in Wind Fields. Journal of Atmospheric and Oceanic Technology, 23, 1614-1620.

[99]   Niiler, P.P. and Paduan, J.D. (1995) Wind-Driven Motions in the Northeast Pacific as Measured by Lagrangian Drifters. Journal of Physical Oceanography, 25, 2819-2830.