ABSTRACT

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.

Keywords:

ROMS, 3D-Var Objective Analysis, Reanalysis Data, Fisheries Science, Particle-Tracking Experiment

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),

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] .

(1)

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

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 (=10⁶ m³・s⁻¹), 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

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-

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

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

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

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

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

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.

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

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,

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.

Acknowledgements

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.

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. and 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. http://dx.doi.org/10.4236/ojms.2017.71006

References

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