J Oceanogr (2014) 70:535–547
DOI 10.1007/s10872-014-0250-8

ORIGINAL ARTICLE

Oceanic dispersion simulation of perfluoroalkyl substancesin the Western North Pacific associated with the Great
East Japan Earthquake of 2011
Yasumasa Miyazawa • Nobuyoshi Yamashita •Sachi Taniyasu • Eriko Yamazaki • Xinyu Guo •Sergey M. Varlamov • Toru Miyama

Received: 15 April 2014 / Revised: 20 September 2014 / Accepted: 25 September 2014 / Published online: 10 October 2014
Ó The Oceanographic Society of Japan and Springer Japan 2014

Abstract The Great East Japan Earthquake on 11 March
2011, followed by the tsunami and fire, resulted in serious
environmental problems in and around Japan. A hugeamount of material was discharged into the ocean after thetremendous flood damage of the tsunami. A monitoring
survey of the perfluoroalkyl substances (PFAs) found evidence showing an abrupt increase in the PFA concentrationin the ocean east of Japan in 2011 after the earthquake.
To confirm the anomalous input of two typical PFAs,
perfluorooctanoate (PFOA) and perfluorooctanesulfonate(PFOS), from the Japanese coast into the ocean, associatedwith the earthquake, we conducted a series of chemical
tracer simulations using an eddy-resolving ocean reanalysis
product: JCOPE2. The simulation model involves processes representing the emission of PFAs from the land
triggered by the tsunami flood, advection of the polluted
waters, and decay of the concentration by the backgroundoceanic turbulence. Comparison of the PFOA simulation
results with the observation confirms a spike-like input of
PFOA into the Western North Pacific after the earthquake.
Advection and diffusion by the Kuroshio Extension and themesoscale eddies play a key role in the dilution of theconcentration. Optimization of unknown simulation

Y. Miyazawa (&)  X. Guo  S. M. Varlamov  T. Miyama
Japan Agency for Marine-Earth Science and Technology,
3173-25 Showa-machi, Kanazawa-ku,
Yokohama 236-0001, Japane-mail: miyazawa@jamstec.go.jp
N. Yamashita  S. Taniyasu  E. Yamazaki
National Institute of Advanced Industrial Science andTechnology, Tsukuba, Ibaraki 305-8569, Japan
X. GuoEhime University, Matsuyama, Ehime 790-8577, Japan

parameters leads to an estimation of the total amount of theanomalous PFOA emission. In contrast, the PFOS simulations are not able to explain the observed distribution,suggesting possible differences in the oceanic transport
processes between PFOS and PFOA.
Keywords Perfluoroalkyl substances  Great East Japan
Earthquake 2011  Oceanic dispersion  In-situ
observation  Numerical simulation  Ocean reanalysis

1 Introduction
Oceanic distribution of perfluoroalkyl substances (PFAs)
including perfluorooctanoate (PFOA) and perfluorooctanesulfonate (PFOS) has been investigated (Yamashita
et al. 2004) since the first discovery of the accumulation of
PFOS in a number of marine wildlife animals (Giesy and
Kannan 2001). PFAs have been manufactured for over
60 years and are widely spread throughout the environment. An international collaborative study for the global
ocean monitoring of PFAs has been carried out since 2002
by the National Institute of Advanced Industrial Scienceand Technology in Japan, the Ocean Research Institute of
the University of Tokyo, and the Leibniz Institute ofMarine Sciences in Germany (Yamashita et al. 2005).
PFAs are important not only for risk assessment in theenvironment, but also for effective tracing of the generalcirculation of the oceans on the global and regional scale(Yamashita et al. 2008).
The international PFA monitoring group recently foundthat the subsurface concentration of PFA at sampling
points east of Japan was much higher in July 2011 after theGreat East Japan Earthquake on 11 March 2011 (EQ 3.11)
than that observed in June 2010 prior to EQ 3.11

123

536

(a)

20100622 144.505E-40.466N

(b)

0
1000
2000
3000

Depth (m)Depth (m)Fig. 1 Vertical profiles of theobserved PFOA and PFOS
concentration. a Measurement
at 144.505°E, 40.466°N on 22
June 2010. b As in a except forthe zoom-up view in the upper
1,100-m layer. c Measurement
143.873°E, 37.866°N on 31 July
2011. d As in c except forzoom-up view in upper 1,100-m
layer

Y. Miyazawa et al.

4000
5000
6000

PFOA
PFOS

7000
8000

0

25

50

75

100

125

150

20100622 144.505E-40.466N
0
100
200
300
400
500
600
700
800
900
1000
1100

PFOA
PFOS
0

25

(d)

20110731 143.873E-37.866N
0
1000

3000

Depth (m)Depth (m)

2000

4000
5000
6000

PFOA
PFOS

7000
8000

0

25

50

75

100

125

PFOA/PFOS (PG/L)(Yamazaki et al. 2012). The PFOA and PFOS profiles
obtained east of Japan by a research cruise of the R/V
Hakuho Maru (KH-10-02) in June 2010 before EQ 3.11
(Fig. 1a, b) indicate the maximum concentration at thesurface together with a few local maxima at deep levels,
implying possible effects of the airborne PFAs near thesurface and re-distribution in the deep layer caused by thedense bottom water around the slope of the Japan Trench(Yamazaki et al. 2012). The maximum concentration values in the subsurface levels of the profiles obtained before
EQ 3.11 are around 20 pg/l. The profiles obtained by acruise of the same R/V (KH-11-07) after EQ 3.11 exhibit asignificant increase in the concentration in the upper 1,000m-depth layer (Fig. 1c, d), even though the surface concentrations of the profiles obtained before EQ 3.11 werehigher than those obtained after EQ 3.11.
The tremendous flood associated with the tsunami triggered by EQ 3.11 resulted in a huge amount of marine
debris floating over a wide range of the North Pacific(Showstack 2011). Various kinds of chemical compounds
including PFAs together with the debris could have beendischarged into ocean from the east coast of Japan afterEQ 3.11; however, so far no studies have focused on
numerical simulations of these phenomena. This study
aims at numerically simulating the oceanic dispersion of
PFAs discharged into the ocean triggered by EQ 3.11.

123

75

100

125

150

PFOA/PFOS (PG/L)

PFOA/PFOS (PG/L)

(c)

50

150

20110731 143.873E-37.866N
0
100
200
300
400
500
600
700
800
900
1000
1100

PFOA
PFOS
0

25

50

75

100

125

150

PFOA/PFOS (PG/L)

It is well known that a huge amount of radionuclides
originating from the Fukushima Daiichi Nuclear Power
Plant was discharged into the Western North Pacific afterEQ 3.11, and the ocean forecast and reanalysis using
numerical models have been intensively utilized for
tracking of the oceanic dispersion of the radionuclides(e.g., Masumoto et al. 2012; Miyazawa et al. 2012, 2013).
Estimation of the source information including the total
emission amounts of the radionuclides into the environment has been a central issue for the investigation of theradionuclide dispersion (e.g., Miyazawa et al. 2013). We
utilized the JCOPE2 ocean reanalysis (Miyazawa et al.
2009) for simulation of the oceanic dispersion of PFAs
associated with EQ 3.11 and estimated the unknown simulation parameters including the source information of
PFAs by adopting the same method, based on Green’s
function approach, as used for the previous study of oceanic radionuclide dispersion (Miyazawa et al. 2013).
This article is organized as follows. Section 2 describes
the transport model of PFAs, observation data, and optimization of simulation parameters using the Green’s
function approach. Section 3 shows the simulation results
and discusses the estimation of the simulation parameters.
First, we focus on the oceanic dispersion of PFOA and later
discuss the validity of the simulation of PFOS. Section 4summarizes the results and discusses related issues.

Oceanic dispersion simulation of perfluoroalkyl substances

2 Model and data
Supposing that PFAs (PFOA and PFOS) are water-soluble
tracers transported by ocean currents, we conducted aseries of the oceanic dispersion simulations for the period
from 12 March 2011 to 30 March 2012 using a transportmodel of PFAs similar to that used for the radionuclide
simulations associated with the accident at the Fukushima
Daiichi Nuclear Power Plant (Miyazawa et al. 2012, 2013).
Considering that PFOA is more water-soluble and less bioaccumulative than PFOS (Nakata et al. 2006), at first, we
focused on the simulation of PFOA, then discussed thesimulation of PFOS. The present model includes a source
flux term representing the emission of PFAs from land into
the ocean caused by the tsunami flood after EQ 3.11 andinvolved no other source/sink term such as the half-life
decay term for the radionuclides because of the chemical
stability of PFAs.
Horizontal and vertical diffusion processes are modeledby harmonic operators with horizontal and vertical diffusion coefficients estimated by the Smagorinsky (1963)
formula and the Mellor-Yamada-Nakanishi-Niino-Furuichi
mixed layer model (Nakanishi and Niino 2009; Furuichi
et al. 2012), respectively. The Smagorinsky (1963) formula
represents the horizontal diffusion coefficient,

537

Japan (Nakajima and Koarai 2011). The PFA emission flux
with a constant value is put on the top level of the source
grid during the period from 12 to 16 March 2011. We put
an equal weight for each grid because a preliminary comparison between cases with equal weight and spatiallydifferent weights depending on the corresponding flood
areas (see shades in Fig. 2) shows the better skill of theformer case (correlation of the observation data: 0.58; astatistically significant level for 95 % significance with asample size of 81 is 0.22) as compared to the latter case(correlation: 0.51). The spatially different weights for thesource grids enhance overestimation of the concentration atsome measurement points near the coast (not shown).
The total emission amount is considered an adjustableparameter. Note that the duration of the emission (5 days)
could also be an adjustable parameter, but it is not independent of the total emission amount. In this study, wechoose not the duration but the total emission amount as an
adjustable parameter, because a preliminary simulation
with a longer emission period (10 days) shows worse skill(correlation: 0.55) than the base case with the 5-day
emission period (correlation: 0.58).
Since the true states of PFA initial conditions areunclear, four different types of vertical profiles are tested as
possible initial conditions in the simulation experiments

1
Kx;y ¼ Cs DxDy rV þ ðrV ÞT ;2 2  2 !1=2
ou
1 ov ou 2
ov
T
þ
þ
þ
ð1Þ
rV þ ðrV Þ ¼ox
2 ox oyoywhere u and v denote horizontal flow; Dx (Dy) denote a gridinterval in the zonal (meridional) direction; Cs is a positive
constant, whose order could be 0.1 (Mellor 2004). The
Smagorinsky diffusion constant Cs is considered anadjustable parameter in this study.
We examine the dispersion of PFAs in the Western
North Pacific using the daily-mean ocean current data
produced by the JCOPE2 ocean reanalysis, which covers aregion between 10.5°–62°N and 108°–180°E with a horizontal resolution of 1/12° and 46 generalized sigma levels
from the surface to the maximum depth of 6,500 m (Miyazawa et al. 2009). The transport model uses the samegrid as JCOPE2. The time step of the present simulation is
5 min, and the daily-mean horizontal current data of
JCOPE2 are linearly interpolated to each time step. The
vertical velocity required for the vertical advection term is
calculated using the continuous equation from the horizontal current data.
To represent the possible emission processes of PFAs,
we assume source grids along the east coast of Japan(Fig. 2) based on the information about the tsunami floodarea provided by the Geospatial Information Authority ofFig. 2 Model grids near the flood-damaged regions associated with
EQ 3.11. Shaded area indicates the flood area (in km2) in each local
region estimated by the Geospatial Information Authority of Japan(Nakajima and Koarai 2011)123

Y. Miyazawa et al.

Depth (m)

538
0
100
200
300
400
500
600
700
800
900
1000
1100

Table 1 Parameters, cost, and correlation in the simulations with theinitial profile ‘basic’

basic
1
2
3
0 10 20 30 40 50 60 70 80 90 100 110 120

PFOA/PFOS (PG/L)

Fig. 3 Vertical profiles of the initial conditions used for thesimulations(Fig. 3). A basic case represents an initial step profile of
PFAs with a constant positive value for the upper 500 m,
which is an adjustable parameter, and zero below 500 m
depth. Other cases are taken from the observed profiles forthe upper 500 m with zero below 500 m depth. Profiles ‘1’and ‘2’ are based on the observed PFOA and PFOS subsurface profiles, respectively (Fig. 1b). To examine theeffects of the relatively high concentration at the surface(Fig. 1), we test profile ‘3’ with a surface value change
from 18 to 120 pg/l for profile ‘2’, which is similar to atypical profile found offshore around industrialized countries (Yamashita et al. 2005). The adjustable parameter forthe profiles ‘1’, ‘2’, and ‘3’ is a multiplication constant forthe original profiles with the default value of 1.
For validation and optimization of the simulations, we
utilize the data obtained by a research cruise of the R/V
Hakuho Maru (KH-11-07), which sampled vertical profiles
during a period from July to August 2011 (Yamazaki et al.
2012). Supposing possible contamination with airborne
PFAs, we exclude the observation data at the surface (0 m)from the validation at first and then examine effects of thesurface data on the validation and optimization results.
Focusing on the upper ocean dispersion processes, we also
exclude the data obtained at depth levels greater than
1,000 m from validation.
To adjust the parameters using the Green’s function
approach (Menemenlis et al. 2005; Miyazawa et al. 2013),
we minimize a cost function,

t


*f
*f
C¼ ~
y  x R1 ~
yx
ð2Þ
! 
t
y ¼ ðyo1 ; . . .yoN Þt denote Nwhere xf ¼ xf1 ; . . .xfN and !
number concentration values of the simulation and observation, respectively; R denotes the observation error
covariance matrix whose diagonal components are givenby the observation errors, which are defined as 20 % of themeasurement values. We ignore the error covariance
between the different observation data.

123

Cases

Basic

b-P-E

b-P-I

b-P-C

Total emission amount (tons)

5.8

11.5

5.8

5.8

Initial value of the upper
500-m depth (pg/l)

20

20

40

20

Diffusion constant
Cost function value for the
PFOA data

0.1
685

0.1
1,127

0.1
1,259

0.3
643

Correlation with the
PFOA data

0.58

0.55

0.61

0.59

Cost function value for the
PFOS data
Correlation with the
PFOS data

21,56543,26971,422

20,508

0.15

0.11

0.21

0.16

Table 2 As in Table 1 except with the initial profile ‘1’
Cases

1

1-P-E

1-P-I

1-P-C

Total emission amount (tons)

5.8

11.5

5.8

5.8

Ratio to the initial profile

1

1

2

1

Diffusion constant

0.1

0.1

0.1

0.3

Cost function value
for the PFOA data

1,007

1,279

991

975

Correlation with the
PFOA data

0.51

0.51

0.50

0.50

Cost function value for
the PFOS data

10,656

29,055

27,698

9,498

Correlation with the
PFOS data

0.063

0.060

0.058

0.057

3 Results
3.1 The PFOA simulations
After performing a number of preliminary simulation
experiments, we determined first guess parameter values:(1) a total emission amount for the 5-day period from 12 to
16 March 2011 as 5.76 tons, (2) an initial constant concentration value in the upper 500-m layer for the basic
profile as 20 pg/l and a default multiplication constant of 1
for profiles 1–3, and (3) the Smagorinsky diffusion constant for the horizontal diffusion parameterization (1) as 0.1(Tables 1, 2, 3, 4).
Comparison of the baseline simulation result (‘basic’ inTable 1) with the first guess parameters and the PFOA
observation data obtained from July to August 2011 suggests a similarity between them (Fig. 4a), and the correlation value is 0.58. The Green’s function approachrequires results of sensitivity experiments with perturbed

Oceanic dispersion simulation of perfluoroalkyl substances

539

Cases

2

2-P-E

2-P-I

2-P-C

Total emission amount (tons)

5.8

11.5

5.8

5.8

Ratio to the initial profile

1

1

2

1

Diffusion constant

0.1

0.1

0.1

0.3

Cost function value for the
PFOA data

852

1,184

727

801

Correlation with the
PFOA data

0.58

0.55

0.62

0.59

Cost function value for
the PFOS data

12,74032,44436,011

11,785

Correlation with the
PFOS data

0.16

0.12

0.22

0.17

200
150
100
50
0

3

3-P-E

3-P-I

3-P-C

Total emission amount (tons)

5.8

11.5

5.8

5.8

Ratio to the initial profile

1

1

2

1

Diffusion constant

0.1

0.1

0.1

0.3

Cost function value for the
PFOA data

805

1,175

827

749

Correlation with the
PFOA data

0.61

0.56

0.65

0.62

Cost function value for
the PFOS data

15,01835,502

44,941

14,099Correlation with the
PFOS data

0.20

0.14

0.27

0.21

target parameters (Menemenlis et al. 2005). We conducted
three sensitivity experiments required for the calculation ofGreen’s function by perturbing each adjustable parameterin the initial profile of the basic type (Table 1). Increasing
the total emission amount (b-P-E) and the initial value (bP-I) results in larger cost function values, while increasing
the Smagorinsky diffusion constant (b-P-C) leads to a cost
reduction. The parameter perturbation experiments with theother initial basic profiles (Tables 2, 3, 4) show that
increasing the total emission amount (1-P-E, 2-P-E, and3-P-E) and the Smagorinsky diffusion constant (1-P-C,
2-P-C, and 3-P-C) leads to larger and smaller constantfunction values, respectively. The responses of the direct
perturbation to the initial basic profiles seem to be more
complicated: reduction (increase) of the cost function
values in 1-P-I and 2-P-I (3-P-I).
As a result of the optimization through the Green’s
function approach accounting for the error covarianceamong all three parameters, the smaller emission amount,slightly larger initial concentration value, and larger Smagorinsky diffusion constant are estimated as optimal if
compared to those of the baseline experiment: ‘A-b-s’

0

50

100

150

200

250

Simulation (PG/L)(b) 250
Observation (PG/L)

Table 4 As in Table 1 except with the initial profile ‘3’
Cases

Observation (PG/L)(a) 250

Table 3 As in Table 1 except with the initial profile ‘2’

200
150
100
50
0

0

50

100

150

200

250

Simulation (PG/L)

Fig. 4 Comparison of the observed and simulated concentration of
PFOA (in pg/l). a The baseline simulation (‘basic’ in Table 1). b The
simulation with the parameters optimized for the PFOA data (‘PFOAb-s’ in Table 6)(Table 5). Use of the other initial vertical profiles results inincreased values of the expected cost function with the
optimized parameters ‘A-1-s’, ‘A-2-s’, and ‘A-3-s’, ascompared to the basic case: ‘A-b-s’ (Table 5). The evaluated emission amounts are similar among all the cases
except for ‘A-1-s’, which shows a larger value, 6.7 tons,
than those of around 5 tons estimated for the other cases.
This is considered as compensation for the low concentration in the near surface layer in the initial profile ‘1’(Fig. 3).
A simulation with the optimized parameters for the basic
case: ‘A-b-s simulation’ (Table 7) shows the actual
reduction of the cost function (685 ? 614) similar to the
expectation (618) based on the Green’s function approach
together with the improvement of the correlation(0.58 ? 0.60), which mainly comes from the corrected
overestimation owing to the decreased emission amount in
the optimized case, as indicated in Fig. 4a, b. Although thecorrelation in the optimized case (0.60) is slightly smaller
than that in the case of b-P-I (0.61), the cost function value
in the former case (614) is sufficiently smaller than that in
the latter case (1,127). Another simulation with the optimized parameters for the initial profile ‘2’ represents thecost and correlation with comparable scores to the basic
case (‘A-2-s simulation’ in Table 7) and shows the

123

540

Y. Miyazawa et al.

Table 5 Optimization results with use of the PFOA observation data
Cases
(expected cost)

InitialprofileUse of surface
data

Emissionamount ± error(tons)

Initial value ± error

Diffusion
constant ± error

A-b-s (618)

Basic

No

4.8 ± 0.5

23 ± 0.8 (pg/l)

0.37 ± 0.06

A-1-s (888)1No6.7 ± 0.5

1.5 ± 0.07 (multiplication constant)

0.49 ± 0.06

A-2-s (621)A-3-s (629)

2
3

No
No

5.0 ± 0.5
5.1 ± 0.5

1.8 ± 0.06 (multiplication constant)
1.5 ± 0.06 (multiplication constant)

0.46 ± 0.06
0.48 ± 0.06

A-b-0 (1067)

Basic

Yes

1.1 ± 0.2

27 ± 0.8 (pg/l)

0.27 ± 0.03

A-1-0 (1499)

1Yes

1.8 ± 0.2

1.8 ± 0.06 (multiplication constant)

0.32 ± 0.04

A-2-0 (1048)

2Yes

1.0 ± 0.2

2.1 ± 0.06 (multiplication constant)

0.30 ± 0.03

A-3-0 (1048)

3Yes

1.0 ± 0.2

1.8 ± 0.05 (multiplication constant)

0.31 ± 0.03

evolution of the oceanic dispersion generally similar to thebasic case (not shown).
Figure 5 shows that the PFOA injected into the Pacific
Ocean after EQ 3.11 is transported eastward by the Kuroshio Extension together with the horizontal mixing
caused by the mesoscale eddies around the KuroshioOyashio mixed water region (Miyazawa et al. 2009;
Masumoto et al. 2012). The processes are basically similar to those associated with the radionuclide dispersion
reported by previous studies (e.g., Miyazawa et al. 2012,2013). Relatively high concentrations of larger than
50 pg/l observed east of Japan during July and August
2011 are explained by the simulation (Fig. 5d, e). Figure 5d suggests that the concentration of two measurement points at 32.0°N–157.9°E and 35.9°N–165.0°E
could be affected by the transport of the discharged PFOA
along the Kuroshio Extension, although the simulation
underestimates the absolute concentration levels. One year
after the beginning of the simulation (Fig. 5f), theincreased concentration almost disappears near Japan as aresult of the eastward transport together with the horizontal and vertical mixing. Figure 5f also indicates that a
low concentration smaller than 20 pg/l emerges around
Japan, originating from the deeper layer; this will be
discussed below.
Figure 6 shows that the high concentration reaches
100-m depth over a wide range east of Japan (Fig. 6a, b),while the high concentration is confined to areas near thecoast at 200 m depth (Fig. 6c, d). The distributions of low
concentrations smaller than 20 pg/l appear at 100 and
200 m depth near the Japanese coast in July and August
2011, and they are partially transported horizontally eastward along the Kuroshio Extension (Fig. 6c, d). The low
concentration originates from the local upwelling associated with the horizontal shear of the Kuroshio in the shelfslope south of Japan (Kasai et al. 2002). The local
upwelling is also shown in the Japan Sea and near thenorthern Kuril Islands (Fig. 6c, d).

123Comparison of the simulated vertical profiles with theobserved ones obtained near the east coast of Japan(Fig. 7a–c; also see Fig. 5e) indicates that the simulation
represents concentrations larger than 50 pg/l at the upper
200-m depth as observed. The simulation represents subsurface concentrations lower than the initial value in theupper 500 m, 23.3 pg/l, at a point around the Kuroshio
Extension region (Fig. 7d), which is caused by the transport of low-concentration waters from the coastal area ofJapan (Fig. 6c). The slight increase of the near-surface
concentration shown in Fig. 7e, f resulted from the transport of the high-concentration waters along the Kuroshio
Extension (Fig. 5d). The simulated profiles at two points(Fig. 7g, h) in the northern latitudes maintain the initialcondition even 4 months after the beginning of the simulation, indicating that the simulated dispersion does notextend there (Fig. 5d).
The optimized parameter in the Green’s function
approach is represented by the summation of the contribution rates from each observation (see Eq. 12 in Miyazawa et al. 2013). Figure 8 shows the distributions of thecontribution rate for the total emission amount (Fig. 8a),the initial concentration value in the upper 500-m watercolumn (Fig. 8b), and the Smagorinsky diffusion constant(Fig. 8c). Measurement points over a wide range in theWestern North Pacific, indicating a positive contribution
rate of more than 1 % for the emission amount parameter(Fig. 8a), suggest that the accidental discharge of PFOA
associated with EQ 3.11 actually affects the PFOA distribution in the Western North Pacific after a few months.
Four points near the Japanese coast (three: 37.0°N–
142.0°E, 38.0°N–143.5°E, 38.0°N–144.0°E) and far fromthere (one: 45.0°N–165.0°E) indicating a negative rate
contribute to the reduced estimate of the optimized emission amount as compared to the first guess (‘A-b-s’ in
Table 5). The measurement points sensitive to the initialvalue adjustment with the absolute values of a contribution
rate larger than 20 % are distributed east of 155°E, far from

Oceanic dispersion simulation of perfluoroalkyl substances

541

(a)

(d)

(b)

(e)

(c)

(f)

Fig. 5 Monthly mean distributions of the PFOA concentration (in pg/l) at 10 m depth represented by the simulation with the parametersoptimized for the PFOA data (‘PFOA-b-s’ in Table 6). Closed circles
in squares indicate the concentration measured in each month. Colorbars of the circles are the same as those of the simulation. a April
2011, b May 2011, c June 2011, d July 2011, e August 2011, f March
2012

123

542

Y. Miyazawa et al.

(a)

(c)

(b)

(d)

Fig. 6 As in Fig. 5, except at 100 m (a, b) and 200 m (c, d). a, c (b, d) denote the distributions in July (August) 2011

the source region (Fig. 8b), where the initial condition
parameter works more effectively for the representation of
the observed feature than the emission amount parameter.
The increased estimation of the optimized diffusion
parameter as compared to the first guess comes mainly
from two measurement points (38.0°N–144.0°E and
36.5°N–150.0°E), while all other points indicate negative
and/or small contribution rates within 20 % (Fig. 8c).
To check the effects of the surface data on the parameter
optimization, we conduct the optimization with the addition of the surface observation data (‘A-b-0’, ‘A-1-0’,‘A-2-0’, and ‘A-3-0’ in Table 5). The results interestingly
indicate considerably reduced emission amounts andslightly increased initial values as compared to the results
without the use of the surface data (‘A-b-s’, ‘A-1-s’, ‘A-2s’, and ‘A-3-s’). Longitudinal plots of both the observed
and simulated PFOA concentration (Fig. 9) depict that thesimulation represents a simple picture of the oceanic

123

dispersion originating from the source region: relatively
higher concentrations in the subsurface layer and at thesurface are confined near the coast west of 145°E. The
subsurface observation generally represents similar features (Fig. 9a). In contrast, the surface observation data
(Fig. 9b) are not consistent with the simple picture ofoceanic dispersion, e.g., a relatively high concentration
larger than 300 pg/l is found east of 155°E far from thesource region. The correlation between the simulation andthe observation becomes low if the observation includes
the surface data (Table 7). The contribution rate distribution for the total emission amount in the optimization caseswith the addition of the surface data (not shown) indicatesalmost negative values near the coast and negligible values
smaller than 1 % far from the coast, suggesting that thedirect discharge from the source region triggered by
EQ 3.11 does not explain the observed surface concentration well. Note that use of the initial vertical profile with

Oceanic dispersion simulation of perfluoroalkyl substances

20110803 141.3333E-36.0000N

0
200
400
600
800
1000

obs
sim

0

50

100

150

(e)
depth [m]depth [m](a)20110719 157.8540E-31.9872N

0
200
400
600
800
1000

200

0

50

50

100

150

(f)
depth [m]depth [m]

obs
sim

0

200

obs
sim

0

50

100

150

(g)
depth [m]depth [m]

obs
sim

0

200

50

100

150

150

200

20110725 165.0040E-45.0073N

(b)
obs
sim

0

(h)
depth [m]depth [m]

obs
sim

0

100

50

100

150

200

PFOA [PG/L]

20110718 150.0125E-36.4890N

0
200
400
600
800
1000

50

0
200
400
600
800
1000

PFOA [PG/L]

(d)

200

PFOA [PG/L]

20110731 143.8737E-37.8668N

0
200
400
600
800
1000

150

20110723 165.0027E-34.9992N

0
200
400
600
800
1000

PFOA [PG/L]

(c)

100

PFOA [PG/L]

20110801 143.5273E-38.2350N

0
200
400
600
800
1000

(a)

obs
sim

PFOA [PG/L]

(b)

543

200

PFOA [PG/L]

20110727 165.0075E-50.9977N

0
200
400
600
800
1000

obs
sim

0

50

100

150

200

PFOA [PG/L]

Fig. 7 Vertical profiles of the PFOA concentration. Black (red)curves denote the observation (the simulation with the parametersoptimized for the PFOA data: ‘PFOA-b-s’ in Table 6). Blue lines
denote the ‘basic’ initial profiles(c)the relatively high concentration at the surface (profile ‘3’),
which implicitly represents the possible effect of the airborne PFAs, does not affect the parameter optimization
results much.
3.2 The PFOS simulations
We conduct the parameter optimization using the PFOS
observation data in the same manner as for the PFOA data
(Table 6). The evaluated total emission amounts, initialvalues, and diffusion constants are smaller than thoseevaluated for the PFOA data in all cases. The evaluatedsmaller level of the PFOS concentration as compared withPFOA is consistent with previous studies (Yamashita et al.
2008). In contrast to the case of PFOA, the baseline simulation (‘basic’ in Table 1) shows poor capability to represent the observed PFOS data (Fig. 8a). Simulations with

Fig. 8 Contribution rates (in %) averaged in 1/2° 9 1/2° grids in the
optimization case ‘A-b-s’. a The total emission amount. b The initialvalue of the upper 500-m depth. c Smagorinsky diffusion constant

123

544

Y. Miyazawa et al.

the optimized parameters result in correlation values
smaller than the 95 % significance level of 0.22 (Table 7;
also see Fig. 10b). Comparison between the simulated andobserved PFOS profiles (Fig. 11) suggests that the simple
assumptions about the accidental emission from the land
and initial value confined in the upper layer are not valid(a) Subsurface
160
140

for the representation of the observed PFOS profiles. In
particular, some of the observed profiles (Fig. 11a–c, f, h)represent two or three local subsurface maximums of the
PFOS concentration, which seem not to be easily simulatedby the relatively simple model adopted in the present study.
Longitudinal plots of the PFOS observation data in thesubsurface layer and at the surface (Fig. 12) exhibit acomplicated feature that might not be simply explained by
the present model (see also the correlation values for thedata including the surface data described in Table 7).

PFOA (PG/L)

120
100

4 Summary and discussion

80
60
40
20
0
140

145

150

155

160

165

160

165

LongitudeObs

Sim(b) Surface
600

PFOA (PG/L)

500
400
300
200
100
0
140

145

150

155

Longitude
Obs

Sim

Fig. 9 Relations between the PFOA measurement longitudes andtheir concentration values. a Subsurface data. b Surface data. Circles(triangles) denote the data of the field observation (the simulation
with the parameters optimized for the PFOA data: ‘PFOA-b-s’ in
Table 6)

Supposing an accidental emission of PFAs (PFOA andPFOS) into the ocean caused by the tremendous flood
associated with the tsunami due to EQ 3.11, we performeda series of the numerical simulations representing theoceanic dispersion using an advection-diffusion model ofthe PFA concentration. The optimization of the parameters
for the PFOA simulation—the total emission amount, initial concentration in the upper layer, and horizontal diffusion intensity—leads to improvement of the simulation
skill representing the observed subsurface profiles of
PFOA. The correlation coefficients between the observation and the simulations with the optimized parameters forthe different types of initial vertical profiles result in arange of 0.60–0.63, and they generally succeed in representing the longitudinal contrast in the PFOA distributionthat could be formed by the accidental emission.
Based on the results, we conclude that the present
model’s formulation is consistent with the real transport
process of PFOA accidentally discharged into the ocean
triggered by EQ 3.11, which could basically be describedas the advection and diffusion processes related to the
Kuroshio-Kuroshio Extension and the mesoscale eddies.
Although the vertical entrainment of the low concentrationwater in the deep layer near the coasts seems to be a little
bit exaggerated by the present setting of the initial

Table 6 Optimization results using the PFOS observation data
Cases (expected
cost)InitialprofileUse of surface
data

Emission amount ± error(tons)

Initial value ± error

Diffusion
constant ± error

S-b-s (1153)

Basic

No

0.8 ± 0.1

3.4 ± 0.2 (pg/l)

0.07 ± 0.03

S-1-s (1285)1No1.8 ± 0.1

0.19 ± 0.02 (multiplication constant)

0.24 ± 0.03

S-2-s (1154)2No

1.0 ± 0.1

0.24 ± 0.01 (multiplication constant)

0.11 ± 0.02

S-3-s (1145)3No

1.0 ± 0.1

0.21 ± 0.01 (multiplication constant)

0.11 ± 0.02

S-b-0 (1660)

Basic

Yes

0.2 ± 0.04

4.0 ± 0.2 (pg/l)

0.09 ± 0.009

S-1-0 (1856)1Yes

0.3 ± 0.04

0.27 ± 0.01 (multiplication constant)

0.10 ± 0.009

S-2-0 (1637)
S-3-0 (1617)

2
3

Yes
Yes

0.2 ± 0.04
0.1 ± 0.04

0.30 ± 0.01 (multiplication constant)
0.26 ± 0.01 (multiplication constant)

0.09 ± 0.009
0.09 ± 0.009

123Oceanic dispersion simulation of perfluoroalkyl substances

545

Table 7 Parameters, cost, and correlation in the simulations with the parameters optimized for PFOA (A-b-s and A-2-s) and PFOS (S-b-s andS-2-s) data
Cases

A-b-s simulation

A-2-s simulation

S-b-s simulation

S-2-s simulation

Total emission amount (tons)

4.8

5.0

0.8

1.0

Initial value

23 (pg/l)

1.8 (ratio)

3.4 (pg/l)

0.24 (ratio)

Diffusion constant

0.37

0.46

0.07

0.11

Cost function value for the PFOA data

614

619

–

–

Correlation with the PFOA data

0.60

0.63

–

–

Correlation with PFOA data including the surface data

0.34

0.37

–

–

Cost function value for the PFOS data
Correlation with the PFOS data

–
–

–
–

1,163
0.17

1,153
0.19

Correlation with PFOA data including the surface data

–

–

0.09

0.10

depth [m]

80
60

obs
sim

0

20110719 157.8540E-31.9872N

0
200
400
600
800
1000

obs
sim

0

10 20 30 40 50 60

PFOS [PG/L]

20

40

60

80

100

Simulation (PG/L)(b)

20110801 143.5273E-38.2350N

0
200
400
600
800
1000

obs
sim

0

100

PFOS [PG/L](f)depth [m]0depth [m](b)2020110723 165.0027E-34.9992N

0
200
400
600
800
1000

10 20 30 40 50 60

obs
sim

0

PFOS [PG/L](c)depth [m]

60
40
20

20110731 143.8737E-37.8668N

0
200
400
600
800
1000

obs
sim

0

(g)

20110725 165.0040E-45.0073N

0
200
400
600
800
1000

10 20 30 40 50 60

obs
sim

0

PFOS [PG/L]
20

40

60

80

100

Simulation (PG/L)

Fig. 10 Comparison of the observed and simulated PFOS concentration (in pg/l). a The baseline simulation (‘basic’ in Table 1). b The
simulation with the parameters optimized for the PFOS data (‘PFOSb-s’ in Table 6)(d)depth [m]020110718 150.0125E-36.4890N

0
200
400
600
800
1000

obs
sim

0

10 20 30 40 50 60

PFOS [PG/L]

condition with zero concentration below 500 m depth, theentrainment processes themselves are possible dynamical
processes caused by the local upwelling associated with thecyclonic circulation in the frontal eddies of boundary currents (Lee et al. 1981) and/or the wind-driven coastal
upwelling (Yoshida 1955).
Our simulation does not represent the PFOS distribution
as observed well; the simulations with the parameters

10 20 30 40 50 60

PFOS [PG/L]

depth [m]

80

0

(e)

10 20 30 40 50 60

40

0Observation (PG/L)

20110803 141.3333E-36.0000N

0
200
400
600
800
1000

10 20 30 40 50 60

PFOS [PG/L](h)depth [m]Observation (PG/L)

100

depth [m]

(a)

(a)

20110727 165.0075E-50.9977N

0
200
400
600
800
1000

obs
sim

0

10 20 30 40 50 60

PFOS [PG/L]

Fig. 11 As in Fig. 7 except for the PFOS observation and thesimulation with the parameters optimized for the PFOS data (‘PFOSb-s’ in Table 6)

optimized for the PFOS observation data show a relatively
low correlation, 0.17–0.19, as compared to the PFOA
cases. The observed subsurface maxima of PFOS are
poorly represented by the PFOS simulation (Fig. 11). The

123

546

Y. Miyazawa et al.

(a) Subsurface
70

PFOS (PG/L)

60
50
40
30
20
10
0
140

145

150

155

160

165

160

165

LongitudeObs

Sim(b) Surface
70

PFOS (PG/L)

60
50
40
30
20
10
0
140

145

150

155

Longitude
Obs

Sim

Fig. 12 As in Fig. 9 except for the PFOS data. Triangles denote thedata of the simulation with the parameters optimized for the PFOS
data: ‘PFOS-b-s’ in Table 6

PFOA simulation also fails to represent the subsurface
maximum around 100–200 m depth (Fig. 7a–c). Comparison of the observed PFOA and PFOS profiles suggests thatPFOS seems to penetrate deeper than PFOA. Those results
suggest that our assumptions of modeling the emission andtransport are not completely appropriate for description of
the actual oceanic transport processes of PFAs. The difference between the capabilites of PFOA and PFOS simulations might be related to the fact that PFOS is not very
soluble in water and is more bio-accumulative as compared
to PFOA (Nakata et al. 2006). More complicated processes
including the sediment dynamics and biological effects
should be included in the transport model of PFAs andcould be investigated by future studies.
Our modeling of the oceanic transport of PFA involves
the uncertainty of the initial condition. The parameter
optimization for the different types of initial vertical profiles causes a difference in the evaluation of the total
emission amounts: 4.8–6.7 (tons) for PFOA and 0.8–1.8(tons) for PFOS. We note that the initial profile types
representing higher concentrations near the surface(‘basic’, ‘2’, and ‘3’), which are similar to the typical123

profiles around industrial countries such as Japan (Yamashita et al. 2005, 2008), clearly exhibit smaller costestimates as compared to the case with an initial profile
with a subsurface maximum concentration (‘1’). The total
emission estimates excluding the case with the profile ‘1’
show narrower ranges of the estimates: 4.8–5.1 (tons) for
PFOA and 0.8–1.0 (tons) for PFOS.
Note that all simulations with different initial verticalprofiles fail to represent the observed subsurface maxima
of both the PFOA and PFOS observation data (e.g., Fig. 7,
11). These seem to be pre-existing states and could not befrom the EQ 3.11-related signals, especially for PFOS,
although the missing mechanisms also could be responsible
for the poor capability of representing the subsurface
maxima, as mentioned above. The available observation
data on PFAs are still few and not enough to even create aclimatological view of the horizontal and vertical distributions. For a more accurate description of PFAs’ oceanic
variability, we need to clarify the climatological view ofthe horizontal and vertical distribution by accumulating
more observation data and/or by performing the simulations designed for representation of the climatologicalstates.
The optimization results depend on the prescription of
the observation error of PFAs. The present study assumes
that the observation error is given by 20 % of the measurement value. In such cases, the optimization tends tounderestimate the comparatively high concentration of theobserved PFOA because the observation data for higher
concentrations and larger observation errors are lessweighted for the cost function than those with lower concentrations and smaller errors. Assuming a constant 1 pg/lof the observation error for all PFOA data, the optimization
results in a smaller emission amount: 4.4 ± 0.06 tons, alarger initial value, 46 ± 0.2 pg/l, and a smaller Smagorinsky diffusion constant, 0.18 ± 0.01. Interestingly, theestimate of the emission amount is not affected much by
the error prescription as compared to that of the initialvalue parameter.
The addition of the surface observation data to the costfunction considerably affects the parameter optimization. A
simple emission process assuming the accidental input of
PFAs into the ocean alone is not enough to represent thecomplicated surface distribution. Figures 9b and 12b suggest that the observed PFOA (PFOS) concentration distribution at the surface has O (100 pg/l) (O (10 pg/l))variability range over the Western North Pacific, which ispossibly caused by the airborne PFOA (PFOS). Surface
boundary conditions representing the supply of the airborne
PFAs into the ocean besides the condition of no flux at thesurface adopted in the present model might be required for
improvement of the capability to reproduce the surface
distribution.

Oceanic dispersion simulation of perfluoroalkyl substances
Acknowledgments This work is part of the Japan Coastal Ocean
Predictability Experiment (JCOPE) promoted by the Japan Agency
for Marine-Earth Science and Technology (JAMSTEC). We thank Dr.
Gamo and the staff of the Atmosphere Ocean Research Institute(AORI), University of Tokyo, for collection of seawater samples.
Comments from two anonymous reviewers were helpful in improving
the earlier versions of the manuscript.

References
Furuichi N, Hibiya T, Niwa Y (2012) Assessment of turbulence
closure models for resonant inertial response in the oceanic
mixed layer using a large eddy simulation model. J Oceanogr
68:285–294
Giesy JP, Kannan K (2001) Global distribution of perfluorooctane
sulfonate in wildlife. Environ Sci Technol 35:1339–1342
Kasai A, Kimura S, Nakata H, Okazaki Y (2002) Entrainment ofcoastal water into a frontal eddy of the Kuroshio and its
biological significance. J Mar Sys 37:185–198
Lee TN, Atkinson LP, Legeckis R (1981) Observation of a Gulf
Stream frontal eddy on the Georgia continent shelf, April 1977.
Deep Sea Res 28A:347–378
Masumoto Y, Miyazawa Y, Tsumune D, Kobayashi T, Estournel C,
Marsaleix P, Lanerolle L, Mehra A, Garraffo ZD (2012) Oceanic
dispersion simulation of cesium 137 from Fukushima Daiichi
Nuclear Power Plant. Elements 8:207–212
Mellor G (2004) Users guide for a three-dimensional, primitive
equation, numerical ocean model. J Phys Oceanogr 34(3):693–698
Menemenlis D, Fukumori I, Lee L (2005) Using Green’s functions tocalibrate an ocean general circulation model. Mon Wea Rev
133:1224–1240
Miyazawa Y, Zhang RC, Guo X, Tamura H, Ambe D, Lee JS, Okuno
A, Yoshinari H, Setou T, Komatsu K (2009) Water mass
variability in the western North Pacific detected in a 15 year
eddy resolving ocean reanalysis. J Oceanogr 65(6):737–756
Miyazawa Y, Masumoto Y, Varlamov SM, Miyama T (2012)
Transport simulation of the radionuclide from the shelf to open
ocean around Fukushima. Cont Shelf Res 50–51:16–29

547
Miyazawa Y, Masumoto Y, Varlamov SM, Miyama T, Takigawa M,
Honda M, Saino T (2013) Inverse estimation of source
parameters of oceanic radioactivity dispersion models associatedwith the Fukushima accident. Biogeosciences 10:2349–2363
Nakajima H, Koarai M (2011) Assessment of tsunami flood situationfrom the Great East Japan earthquake. Bull Geospat Inf Auth Jpn
59:55–66
Nakanishi M, Niino H (2009) Development of an improved turbulence closure model for the atmospheric boundary layer.
J Meteorol Soc Jpn 87:895–912
Nakata H, Kannan K, Nasu T, Cho H-S, Sinclair E, Takemura A(2006) Perfluorinated contaminants in sediments and aquaticorganisms collected from shallow water and tidal flat areas of the
Ariake Sea, Japan: environmental fate of perfluorooctane
sulfonate in aquatic ecosystems. Env Sci Tech 40:4916–4921
Showstack R (2011) Oceanographer tracks marine debris from theJapan tsunami and other incidents. EOS Trans AGU 92(37):306
Smagorinsky J (1963) General circulation experiments with theprimitive equations. I. The basic experiment. Mon Wea Rev
91:99–164
Yamashita N, Kannan K, Taniyasu S, Horii Y, Okazawa T, Petrick G,
Gamo T (2004) Analysis of perfluorinated acids at partsperquadrillion levels in seawater using liquid chromatographytandem mass spectrometry. Environ Sci Technol 38:5522–5528
Yamashita N, Kannan K, Taniyasu S, Horii Y, Petrick G, Gamo T(2005) A global survey of perfluorinated acids in oceans. Mar
Pollut Bull 51:658–668
Yamashita N, Taniyasu S, Petrick G, Wei S, Gamo T, Lam PKS,
Kannan K (2008) Perfluorinated acids as novel chemical tracersof global circulation of ocean waters. Chemosphere70:1247–1255
Yamazaki E, Yamashita N, Taniyasu S, Petrick G, Tanhua T, Kannan
K, Gamo T (2012) Accidental input of perfluoroalkyl substancesfrom land to ocean by the disaster of tsunami, EQ3.11, Japan.
Organohalog Compd 74(B1B):103
Yoshida K (1955) Coastal upwelling off the California coast. Rec
Oceanogr Wks Japan New Ser 2:8–20

123