Short-term dispersal of Fukushima-derived radionuclides off Japan: modeling efforts and model-data intercomparison

. The Great East Japan Earthquake and tsunami that caused a loss of power at the Fukushima nuclear power plants (FNPP) resulted in emission of radioactive isotopes into the atmosphere and the ocean. In June of 2011, an international survey measuring a variety of radionuclide isotopes, including 137 Cs, was conducted in surface and subsurface waters off Japan. This paper presents the results of numerical simulations speciﬁcally aimed at interpreting these observations and investigating the spread of Fukushima-derived radionuclides off the coast of Japan and into the greater Paciﬁc Ocean. Together, the simulations and observations al-low us to study the dominant mechanisms governing this process, and to estimate the total amount of radionuclides in discharged coolant waters and atmospheric airborne radionuclide fallout. The numerical simulations are based on two different ocean circulation models, one inferred from AVISO altimetry and NCEP/NCAR reanalysis wind stress, and the second generated numerically by the NCOM model. Our simulations determine that > 95 % of 137 Cs remaining in the water within ∼ 600 km of Fukushima, Japan in mid-June 2011 was due to the direct oceanic discharge. The estimated strength of the oceanic source is 16.2 ± 1.6 PBq, based on minimizing the model-data mismatch. We cannot make an accurate estimate for the atmospheric source strength since most of the fallout cesium had left the


Background
In March of 2011 the Great East Japan Earthquake, a massive offshore earthquake, and resulting tsunami and aftershocks that wreaked devastation upon northern Japan also caused a loss of power at the Fukushima Dai-ichi nuclear power plants (FNPP) situated on the coast at 37 • 25 N, 141 • 2 E. The subsequent destabilization caused direct emission of radioactivity through explosive release of radioactive materials into the atmosphere and leakage of coolant water discharge into the ocean (Chino et al., 2011;Butler, 2011;NSCJ, 2011).The atmospheric contamination and fallout were measured both near the site (Uematsu et al., 2012;MEXT, 2011a;Toyoshima et al., 2011) and globally (e.g.Preparatory Commission for the Comprehensive Nuclear-Test-Ban Treaty Organization, see Fig. 1 of Kristiansen et al., 2012; Table 5 of Stohl et al., 2012;National Atmospheric Deposition Program;Fig. 1 of Wetherbee et al., 2011), and while the nearshore (within 30 km of the coast) ocean discharge was monitored by both the Japanese government (MEXT, 2011b) and the power company (TEPCO, 2011) for the first few Published by Copernicus Publications on behalf of the European Geosciences Union.
months, the spread of radiation into the greater Pacific went largely unobserved (Buesseler et al., 2011).

KOK cruise
Not quite three months after the first leakage into the ocean (30 May-16 June 2012), a cruise of international collaborators from 13 different institutions aboard the R/V Ka'imikai-o-Kanaloa (KOK) occupied an offshore survey area (Fig. 1b).Extending out 600 km from the coast (34-37 • N, 142-147 • E), this cruise took samples to measure a variety of radionuclide isotopes including 137 Cs and 134 Cs in surface and subsurface waters and in biota (Buesseler et al., 2012).The KOK cruise also included a physical oceanography component focused on the dynamics of the ocean waters advecting and mixing the radioactive discharge.As part of this component of the project, 24 drifters were deployed to capture circulation features contemporaneous with the cruise and to follow the evolution of surface water parcels into the future.Hydrographic profiles (temperature, salinity) accompanied by the oxygen, nutrient and chemical isotope sampling were performed, and underway, direct velocity measurements were procured using the shipboard Acoustic Doppler Current-meter Profiler (SADCP).

Goal
Here, we use numerical modeling to investigate the dominant mechanisms governing the short-term spread of radiation within the North Pacific to interpret the physical components of the KOK observations, and to place the radioactive isotope concentration estimates (Buesseler et al., 2012) in the context of the ocean circulation.The observed data set from the KOK cruise and the numerical models used are presented in the next section.Section 3 discusses the simulated spread of radioactive isotopes, specifically 137 Cs, from both the direct ocean discharge and from the airborne fallout using observation-based 2-D velocity fields.Section 4 considers the effects of 3-D circulation on the spreading of 137 Cs.The last section concludes with a discussion of the results.

KOK cruise data
To investigate the Lagrangian particle dispersion from the Fukushima region, 24 satellite-tracked surface drifters equipped with GPS, temperature and drogue tension sensors were released in June 2012 as part of the KOK cruise.The drifters were deployed at 20 locations, with 2 sets of double releases and one triplet release at the near-shore stations.Figure 1a shows drifter tracks from the time of release through 3 August 2012.All buoys were initially attached to drogues and advected by near-surface currents at 15 m depth.In early August, after traveling southward with a weak coastal cur-rent, 8 drifters grounded on the Japanese coast.Another 2 drifters lost their drogues within a few months to a year of their initial release.This analysis uses the quality-controlled drifter data, optimally interpolated to 6 h intervals1 .
CTD casts were performed and water samples for the radionuclide isotope analysis were collected at 50 stations spanning the area extending from 50 up to 600 km off Japan.Additional underway water sampling was performed between the stations.Collected samples revealed the presence of 137 Cs and 134 Cs (Buesseler et al., 2012).As indicated by the close-to-unity ratio of 134 Cs to 137 Cs isotopes, which have quite different half-lives (roughly 2 and 30 yr, respectively), it was concluded that all 134 Cs found in KOK samples was Fukushima-related (Buesseler et al., 2012).This result is consistent with Honda et al. (2012) who also reported that the ratio of 134 Cs to 137 Cs isotopes was about 1 in many seawater samples taken one month after the accident.A map of 134 Cs concentrations in surface waters resulting from the KOK cruise is shown in Fig. 1b, and the vertical profiles of 134 Cs are shown in Fig. 1c.Note that although we use 134 Cs to identify Fukushima waters, the short-term spreading of the two Cs isotopes in the ocean would be similar, so they can be used interchangeably in the model.

Source functions
There were two distinct sources of Fukushima-related contamination to the ocean: the localized direct discharge of radioactive coolant waters into the near-shore ocean (here referred to as the "oceanic source"), and more widespread fallout of airborne radionuclides that were released to the atmosphere during FNPP explosions and then precipitated into the ocean (referred to as the "atmospheric source").The concentration of 137 Cs in the coolant waters was monitored and the resulting time series has been reported in Buesseler et al. (2011).We use this measured time series (Fig. 2 top) to represent the oceanic source in our modeling efforts.What is less well constrained is the total amount of the released 137 Cs or, equivalently, the magnitude of the direct discharge source.Available direct discharge source estimates vary from 3.5 PBq according to Tsumune et al. (2011) to 22 PBq reported by Bailly du Bois et al. (2012).Our results suggest that the intermediate source strength of about 16.2 PBq gives the best agreement between the KOK observations and our modeling efforts, so we use this value for our simulations in Sect.3.2.We will come back to the question of total source strength in Sect.3.4.
The temporal and spatial distributions of the precipitated 137 Cs over the Pacific Ocean were, for the most part, unmeasured.The source function for the airborne 137 Cs thus has to be deduced from numerical simulations, such as those of Chino et al. (2011), Morino et al. (2011), or Stohl et ).However, we use the total atmospheric source strength that gives the best agreement with the KOK measurements.The highest limit of our estimated atmospheric source strength is 11 PBq, which corresponds to about 40 % of that estimated by Stohl et al. (2012).We use this source strength of 11 PBq for our simulations in Sect.3.3.More details on the source strength estimation will be given in Sect.3.4.

Modeling the spreading of 137 Cs
Over the last two years there have been a number of numerical modeling studies aimed at assessing the aftermath of the Fukushima disaster (e.g., Tsumane et al., 2011Tsumane et al., , 2013;;Kawamura et al., 2011;Estournel et al., 2012;Miyazawa et al, 2012;Masumoto et al., 2012;Bailly du Bois et al., 2012).Forward numerical models have allowed simulation of the spread of radioactive waters into the ocean.Inverse modeling efforts have been useful in constraining the source amplitudes of both the oceanic direct discharge and atmospheric fallout sources.The main goals of this investigation differ from this prior work in two respects.First, we focus on providing the physical oceanographic interpretation of the patterns observed by the KOK survey (Fig. 1b).Second, we investigate the general utility of data-based velocity fields for realistically reproducing the dispersal of passive tracers in the ocean over a time scale of several months.We use two different estimates of the oceanic circulation to simulate the spreading of contaminated waters.The first one is observation-based and combines daily nearsurface geostrophic currents (on a Mercator 1/3 • grid) from AVISO with 6-hourly Ekman velocities (on a 2 • × 2 • grid) based on the NOAA NCEP/NCAR wind stresses2 .Wind stress at 10 m height, τ , was converted to ocean velocity at 15 m depth, u Ek and v Ek , using the Ralph and Niiler (1999) formula u Ek + iv Ek = βe −iθ /(fρ) (τ x + iτ y )/ √ |τ |, where ρ = 1.027 kg m −3 is the assumed seawater density, f is the Coriolis parameter, θ = 55 • is the rotation angle of the Ekman current, and β = 0.065 s −1/2 .The pros of this observation-based velocity are its reliability and its global coverage; the cons are the sparse temporal and spatial resolution, its two-dimensional nature, and the absence of ageostrophic components (other than Ekman).To account for lateral diffusion and the influence of the unresolved scales, a small stochastic velocity (taken, with random sign, from the normal distribution with a standard deviation of 5 cm s −1 ) was added to the sum of geostrophic and Ekman velocities.We have carefully checked that the results of our simulations  are not sensitive to the specifics of this stochastic velocity field, but we list these details for the sake of completeness.
The second velocity field that we used comes from the Navy Coastal Ocean Model (NCOM), a high-resolution numerical model (Barron et al., 2004(Barron et al., , 2006).A regional model with 3 km horizontal resolution was nested with open boundaries within the HYCOM global 1/8 • model.The model had a hybrid vertical coordinate system, with 15 z levels at the top and 35 density-defined levels underneath, for a total of 50 vertical levels.The model was forced with wind and heat fluxes from the Coupled Ocean/Atmosphere Mesoscale Prediction System for the Western Pacific (COAMPS WPAC).Tides at the boundaries were provided by the Oregon State Model (Egbert and Erofeeva, 2002).SSH, SST and available T −S profile data from the Naval Oceanographic Office were assimilated into this model using optimal interpolation.The model was run each day, assimilating the data from the previous day and providing a 48 h forecast.One-day segments in the beginning of each run were then stacked together to create a longer time series covering the time interval from Biogeosciences, 10, 4973-4990, 2013 www.biogeosciences.net/10/4973/2013/mid-March to the end of June 2011.There were discontinuities in the velocity field between the end of one day and the beginning of the next day, but these were small and for our purposes did not create any known issues.Unlike our observation-based circulation, NCOM 3-D velocities varied with depth throughout the water column, and had the advantage of better spatial and temporal resolution than provided by observations.The disadvantages of NCOM were the smaller model domain and the limited ability to match specific measured oceanic features, such as the exact position of the Kuroshio Current and the mesoscale eddies present during the spring of 2011, with the model estimate of the circulation field.
Comparison between the observation-and NCOM-based velocity fields is presented in Fig. 3. Compared to the observation-based model, NCOM generally overestimates the mean Kuroshio velocity (by up to 0.6 m s −1 in some areas), but slightly underestimates the mean currents throughout the rest of the domain.The general shape of the mean Kuroshio Current is captured relatively well in NCOM, but the exact position of meanders and mesoscale eddies is misplaced slightly in NCOM compared to observations.The variability in the eddy velocities is of the same order in the two models, with the highest variability in the general area of the Kuroshio Extension.
Spreading of 137 Cs is modeled using a Lagrangian framework by repeatedly releasing large numbers of simulated water parcels inside the source domain, over the full duration of the source time series.These water parcels are advected by the velocity fields described above and their trajectories are estimated using a fixed-step (RK4 for runs with the stochastic velocity component) or variable-step (RK4(5) for runs without the stochastic velocity component) Runge-Kutta integration scheme with bilinear interpolation in time and space between grid points.The exponential decay of 137 Cs concentration from the initial source value, with half-life of 30.16 years, is applied to estimate the concentration of 137 Cs following each water parcel3 .Our model does not take into account the small fraction of 137 Cs in the particulate phase.The Lagrangian model provides an intuitive framework that illuminates the physical mechanisms by which the contaminated waters were brought from the source region to their position at any given time.The disadvantage of this framework is its numerical intensity due to the large number of the released water parcels and the simplified way in which the calculation treats the diffusion process (diffusivity is assumed to be spatially uniform and isotropic).As the number and frequency of the released parcels increases, the resulting Lagrangian distribution of 137 Cs approaches that estimated from an Eulerian calculation with the corresponding value of diffusivity.
To validate the process, the results of our simulations were tested to ensure that they were not overly sensitive to the further increases in the number of the released parcels.Using twice the number of fluid parcels leads to only a small change in the simulated June concentrations of Cs (see Appendix B).The lack of sensitivity to such an increase suggests that the observed distribution of 137 Cs is fairly close to the limiting value (i.e., for infinitely many parcels).
3 Short-term dispersal of 137 Cs

Real and simulated drifters
To validate, and gain some confidence in, the velocity fields, we start this section by comparing the measured drifter trajectories (Fig. 4a) to those computed using the observationbased and numerically simulated velocities (Fig. 4c and  e).Since the focus here is on the short-term dispersal of 137 Cs and on comparison with the KOK cruise observations, drifters are only compared to our models until 30 June.By that time, out of 24 drifters released, 11 were entrained in the Kuroshio and moving eastward with the current, 1 was captured by the cyclonic mesoscale eddy south of the Kuroshio, and 12 were still recirculating near-shore and remained west of 144 • E. These statistics are fairly well captured by the observation-based simulated drifters (Fig. 4c), in which 14 drifters were entrained into the Kuroshio, 7 recirculated near-shore, 2 ended up south of Kuroshio (one of which was trapped by the model eddy representing the observed eddy seen in the upper panel, 4a), and 1 drifter headed north.The smaller number of simulated drifters staying in the near-shore area is not surprising given that the near-shore circulation features are not well resolved by altimetry.The overall shape and extent of the distribution of the simulated drifters (Fig. 4c) is similar to that observed (Fig. 4a) with the Kuroshio meanders and the cyclonic mesoscale eddy at the correct locations.Note, however, that the simulated drifters do not travel as far east as the real drifters.Comparison of real and NCOM-simulated drifters (Fig. 4e) shows marginal agreement: 6 drifters entrained into the Kuroshio, 1 drifter went south of Kuroshio (this one, however, shows clockwise rotation), and 17 drifters stayed near-shore.The axis of the Kuroshio Current and the positions of the meanders and the eddy are displaced slightly compared to Fig. 4a and  c.The Kuroshio in the NCOM model seems to be less efficient in entraining and transporting fluid, while the nearshore recirculation features are extremely effective in keeping fluid from moving eastward.To make the comparison between the real and simulated drifters more quantitative, the domain is divided into 1   and observation-based maps is 0.95, and between the real and NCOM-based map is 0.89, suggesting that both approaches are reasonable, but that the observation-based approach provides a slightly more realistic outcome.

Spreading of 137 Cs from the ocean source
The spreading of a passive tracer in the ocean is governed by two processes, advection (stirring) and diffusion (mixing) (Eckart, 1948).Advection tends to distort the initial distribution of a tracer by stretching and folding, which produces elongated filaments or streaks that are then stretched and folded again.As the streaks become longer, narrower and more convoluted, gradients of concentration get larger and diffusion comes into play.Concentrated streaks start to fade, diffuse away and finally disappear, leaving behind a wellmixed homogenized fluid.In an ideal situation, when the advective velocity is known at all scales, mixing is simply due to the molecular diffusion.In reality, however, velocity is rarely known at small scales so the un-or under-resolved submesoscale features have to be accounted for by increased diffusion rather than by advection.
The 137 Cs concentrations resulting from the 2-D observation-based model have dimensions of Bq m −2 and should be interpreted as depth-averaged values, c tot .Comparison between the mean vertical profile of 137 Cs and the ARGO-based estimates of mixed-layer depth, z ML indicates that the bulk (> 90 %) of the observed 137 Cs lies above the April mixed layer (see Fig. 1c).Thus, the KOK observations suggest that, to a first approximation, the surface concentrations, c surf (in Bq m −3 ), can be obtained from the depthaveraged values by dividing the latter by the mixed-layer depth, i.e., c surf = c tot /z ML .In our simulation we rely on z ML estimates from ARGO profiles (Holte and Talley, 2009;Holte et al., 2010), but climatological mixed layers estimated from historical hydrographic data yield similar results.Due to the rapid warming and restratification of surface waters in late spring, the mixed-layer depth over the western Pacific Ocean in April/June is roughly one-half/one-fifth the March value.The ARGO-based March, April and June mean z ML estimates averaged over the KOK area are indicated by gray lines in Fig. 1c.The value of c surf is inversely proportional to z ML and thus when z ML is half as deep, c surf is doubled.In situations with temporally changing mixed-layer depths, such as the western Pacific waters in late spring, it is the deepest mixed layers that should be used for the conversion, because once 137 Cs has mixed down, it cannot un-mix and reconcentrate in the surface ocean, even when the mixed layer shoals later in spring.Additionally, this particular area is subject to little upwelling so fluid parcels generally do not rise up through the water column from depth.Rather, subsurface Cs left below the bottom of the receding mixed layer will be advected by subsurface currents at depth and will flow, to   a first approximation, along subsurface isopycnal surfaces.

Biogeosciences
Since the peaks of the atmospheric and oceanic sources occurred on 18 March and 7 April, respectively, we use the March/April z ML values in modeling the spreading of 137 Cs from the atmospheric/oceanic sources.The oceanic and atmospheric sources are markedly different from each other and so are the resulting 137 Cs distributions (Figs. 5 and 6).As discussed in Sect.2.2.we use an oceanic/atmospheric source strength of 16.2/11 PBq.Notice, however, that the spreading pattern does not depend on the choice of the source strength.The oceanic source discharged large concentrations of 137 Cs into the Pacific Ocean directly off the coast of Japan.In our numerical simulations, the source domain was taken to be a 20 × 20 km 2 square centered at the FNPP location.The discharge peaked around 7 April 2011, but the coolant waters continued to leak from the reactor for several months, although with smaller 137 Cs concentrations (Buesseler et al., 2011).The oceanic source is thus prolonged in time but compact in space.
The discharged waters, stirred by the oceanic currents, form a wiggly streak of high 137 Cs concentration, extending approximately eastward off the coast of Japan and into the Pacific (Fig. 5 top left).This is in general agreement with the modeling results of Miyazawa et al. (2012) who reported that the modeled distribution of 137 Cs from the beginning of April and onward is strongly affected by the eastward-flowing Kuroshio Current.With time, the streak becomes increasingly convoluted due to the action of advection, and fuzzier as diffusion starts to mix together the high-concentration streak water with its lower-concentration surroundings (Fig. 5 top left through bottom right).By mid-June (the time of the KOK cruise), the contaminated waters reach 170 • E and 42 • N, staying mostly to the north of the Kuroshio (Fig. 5 bottom right).The homogenization process is well underway by this time, especially at the eastern end of the distribution, although the streaks have not been completely wiped away and are still quite strong in certain areas.Notice that even in June, the streak is still attached to the coast by one end because the coolant waters continued to leak out of the reactor.The amount of 137 Cs left in the KOK area at the time of the cruise is about 17.5 % of the total oceanic discharge, with 14 % in the offshore KOK area spanned by stations 1 to 17.The southern edge of the distribution is roughly aligned with the Kuroshio Current, suggesting that this strong current acts as a transport barrier, particularly over the western part of the Pacific Ocean where it is most strongly defined.Interestingly, Rypina et al. (2010) found that in the North Atlantic, the Gulf Stream also acts as a transport barrier.The similarity in these results suggests that plausibly it is typical for all strong oceanic currents to prevent cross-stream transport.
The distribution of 137 Cs also shows elevated values near the perimeters of several eddies, whereas concentrations near the eddy cores are generally small.Rypina et al. (2010) similarly observed an increase in phytoplankton concentration at the perimeter of coupled eddies (dipoles) in the Philippine Archipelago.Fluid exchange processes are different at the perimeter of an eddy compared to its core.At eddy perimeters the so-called "chaotic advection" mechanism is at play (Mancho et al., 2006;Rypina et al., 2010), which leads to the sequential stretching and folding of fluid elements and facilitates the exchange with the outer waters.Within the Lagrangian core of an eddy, on the other hand, a more "regular" state of motion exists where the behavior of fluid parcels is similar to that found in steady flows.Similar to the steady case, fluid parcels near the core follow paths around the center of the eddy without much exchange with outside fluid.This is why eddies are generally able to transport the water in their cores over long distances without "spilling" it along the way (Chelton et al., 2007;Early et al., 2011;Beron-Vera et al. 2013).The ability of oceanic eddies to trap fluid in their Lagrangian cores has been a subject of an ongoing research.
The three properties discussed above -the streakiness of the distribution, the role of the Kuroshio Current in preventing the southward spreading of 137 Cs, and the high 137 Cs concentrations at the perimeters of eddies -are responsible for several key features in the observed 137 Cs distribution.We will come back to these points in Sect.3.4.

Spreading of 137 Cs from the atmospheric source
The atmospheric source has a much wider spatial pattern than the oceanic source (Fig. 2), but the peak of the atmospheric source occurred about 3 weeks prior to the oceanic source peak, and the temporal duration of the atmospheric fallout    filamentation at these latitudes.Unlike the oceanic source that lies entirely to the north of Kuroshio, the southern edge of the initial fallout distribution of 137 Cs lies to the south of the Kuroshio (Fig. 2b).Thus, the barrier effect of the current from this area to the east where they have been vigorously stirred and mixed with the surrounding water.The northern part of the 137 Cs distribution is less distorted from its initial shape than the southern part, with a large-scale blob of high 137 Cs concentration remaining at around 152 • E and 45 • N even as late as mid-June.The near-shore coastal area in front of the FNPP has low concentrations of 137 Cs and, unlike for the oceanic source, the streaks do not reach the coast in this region.The amount of 137 Cs left in the KOK study area at the time of the cruise is about 4.5 % of the atmospheric source, with most (3.5 %) in the offshore KOK area spanned by stations 1-17.

Source amplitudes
With the numerically estimated distributions of 137 Cs from the two sources and the KOK data in hand, the source amplitudes for the atmospheric and oceanic sources can be estimated.For various amplitudes of both the atmospheric and ocean sources, A ocean and A atm , we construct the simulated 137 Cs map by "sampling" the numerical distributions of 137 Cs from the two sources at the times and positions of the KOK stations.The source amplitudes, A ocean and A atm , are then determined by minimizing the error function which represents the standard L 1 -norm and quantifies the difference between the measured (c i KOK ) and simulated (c i ocean and c i atm ) concentrations of 137 Cs, averaged over the stations.Over the range of source amplitudes 0 < A ocean < 35 PBq and 0 < A atm < 25 PBq, the error function Er has a minimum at A ocean =16.2 ± 1.6 PBq (estimates range from 9.1 PBq to 17.8 PBq with more values near the upper limit) and A atm = 0.5 ± 2.7 PBq (estimates range from 0 PBq up to 11 PBq with most values near the lower limit).The error bars correspond to one standard deviation (and value ranges indicate the smallest and largest estimates, respectively) computed using the bootstrapping procedure (Tukey, 1958) of sequentially eliminating one station and reevaluating the expression for the remaining stations.Note that theoretically a situation is possible where the error function could have two (or more) minima, e.g. one with a higher atmospheric and lower oceanic source estimate, and another with a lower atmospheric and higher oceanic source estimate.In practice, however, in every case the error function has a single minimum over the considered range of source amplitudes (see Appendix A).
Our estimates of the oceanic source strength, ranging from 9.1 PBq to 17.8 PBq are significantly higher than the estimate of 3.5-3.6PBq obtained by Tsumune et al. (2011Tsumune et al. ( , 2013)), and are roughly consistent with estimates of Bailly du Bois et al. (2012), who reported a range of 10-34 PBq for the direct discharge estimates.Our atmospheric source strength estimates, 0-11 PBq, correspond to less than half of the Stohl et al. (2012) estimates, who reported a range of 23-50 PBq for the atmospheric fallout.
Our numerical simulations predict that most (> 95 %) of 137 Cs left in the KOK area in mid-June came from the direct oceanic discharge.This number was obtained by combining the lower/upper limit of our oceanic/atmospheric estimate with the percentage of 137 Cs remaining in the KOK area (17.5 % and 4.5 % for the oceanic and atmospheric releases, respectively; see Sects.2013) who concluded, based on the radium isotope ratios, that the water sampled during the KOK cruise shows a strong coastal signature.This conclusion is further supported by the 90 Sr measurements at the KOK stations (Casacuberta et al., 2013).The distribution of 90 Sr, which is much less volatile than 137 Cs and thus came almost entirely from the oceanic source, is similar to that for 137 Cs, suggesting that both chemicals came from the same source.
It is important to note that our estimate for the atmospheric source amplitude is much less reliable than that for the oceanic source because most of the fallout atmospheric 137 Cs had moved out of the survey area by mid-June.Thus, the oceanic isotope concentrations observed during the KOK cruise, which are due almost entirely to the oceanic discharge of FNPP coolant waters, are ill-suited for constraining the total magnitude of the atmospheric input.

Physical oceanographic view on the key features of the measured 137 Cs map
Using modeled concentrations of 137 Cs in the mixed layer at times corresponding to select stations (Fig. 7), we can interpret some of the key features of the measured 137 Cs map.Since ≥ 95 % of 137 Cs remaining in the KOK area in June is derived from the direct oceanic discharge, we focus on the oceanic source in this subsection.First, the strong and direct influence of the Kuroshio Current, which acts as a transport barrier for the southward progression of 137 Cs from the oceanic source, manifests itself in the low concentrations at stations 1-3 (top row in Fig. 7) and 17-21 (2nd row in Fig. 7).
All of these stations lie to the south of the instantaneous Kuroshio core, and thus the contaminated water from the oceanic source cannot reach them.This pattern agrees with the modeling results of Tsumane et al. (2013) (their Figs. 15 and 16), who also found the smallest Cs concentrations at southernmost KOK stations.Second, the largest 137 Cs concentrations are found not at the stations closest to the FNPP but at stations 37, 46 and 41, 42 (bottom row in Fig. 7).The elevated 137 Cs values at these stations are associated with the action of the near-shore semi-permanent eddy that was present in the area from April to July.This strong and robust cyclonic eddy was entraining 31 Fig. 7. Modeled mixed-layer 137 Cs concentration from the oceanic source (in color) at times corresponding to select stations.Black dot shows the station position.The number outlined in the lower left corner indicates station number (see Figure 1d).the contaminated water that leaked from the reactor into the cyclonic motion around the eddy perimeter.This circulation pattern not only built up the concentrations near the perimeter, but also delayed the eastward advection of 137 Cs.The 137 Cs-rich waters did not penetrate into the center of the eddy, so its core remained relatively 137 Cs-free.
The third interesting observation concerns the westernmost stations 33 through 37, which sampled the streak of coolant water supplied by the continued leakage from the reactor in June (3rd row in Fig. 7).Again, this observation is in general agreement with Fig. 15 from Tsumane et al. ( 2013), which also shows the meridional elongated streak of contaminated water extending parallel to the coast about 40 km offshore (although the exact location of the streak seems to be shifted slightly to the east in their simulations compared to our modeling results and to the KOK observations).The levels of 137 Cs at these stations are consistently elevated, but are not the highest values measured because the concentration of 137 Cs in the coolant water in June was already much reduced compared to its peak value in April.In agreement with the KOK measurements, our modeling results suggest that the distribution of 137 Cs in mid-June was still filamentary and that the diffusion processes had not yet wiped away the high gradients in the 137 Cs field.
The station-by-station comparison between the measured and simulated 137 Cs concentrations at the KOK stations is shown in Fig. 8. Here, the detection limit for 137 Cs is 1.5 Bq m −3 (Buesseler et al., 2012) so we have cut off our simulated 137 Cs concentrations at the same level for a fair comparison.Note also that unlike our simulated concentrations, which correspond to Fukushima-derived 137 Cs, the measured values include both Fukushima-derived and pre-existing 137 Cs (the latter is at levels of 1-2 Bq m −3 ; Buesseler et al., 2011Buesseler et al., , 2012)).The two curves in Figure 8 are in qualitative agreement with each other, showing minima/maxima at same stations (i.e., at same geographical locations).The quantitative agreement is better for stations 1-27 than for the near-shore stations 28-51.This is plausibly because the distribution is extremely filamentary close to shore, more so than for the off-shore stations. 137Cs concentrations in the near-shore area span several orders of magnitude, with very sharp changes in concentration from one station location to the next (Fig. 1b).Our observation-based model is not able to quantitatively match all of this variability due to its coarse spatial and temporal resolution, and the oversimplified www.biogeosciences.net/10/4973/2013/Biogeosciences, 10, 4973-4990, 2013  representation of diffusion processes.Despite some quantitative mismatch, the model is still useful in explaining the qualitative features of the observed Cs distribution.To conclude, despite several limitations, the numerical modeling provided a means for interpreting the KOK observations in the context of the physical ocean circulation.

Effects of the third dimension on the short-term dispersal of 137 Cs
The modeling results in the previous section relied on the two-dimensional near-surface observation-based velocity fields and ignored the effects of vertical velocity, vertical mixing and vertical shear.In this section we use the NCOM numerical model to investigate the effects of these processes.We again restrict our attention to the spreading of 137 Cs from the oceanic source only.As in the previous section, we use the measured time series of 137 Cs concentration in the coolant waters with a source strength of 16.2 PBq and assume that the source occupies (horizontally) a 20 × 20 km 2 square centered at the FNPP location.In our 3-D analysis, however, the source concentrations are assigned not only to the water parcels at the surface, as in the previous section, but also at 10 evenly distributed subsurface layers spanning the upper 10 m of the water column.We then advect water parcels using 3-D NCOM velocities, and as before, apply the exponential decay of 137 Cs with half-life of 30.16 yr following each water parcel.To account for vertical mixing within the near-surface mixed layer4 , the depths and 137 Cs concentrations for water parcels within the mixed layer are adjusted every 6 h.In particular, the domain is divided into 0.1 • × 0.1 • degree lateral bins.In each bin the average mixed-layer con-centration is computed and then assigned to all parcels located within the mixed layer in that bin.The depths of the mixed layer parcels are also adjusted to be randomly distributed over the mixed layer in each bin.The number of parcels is kept the same throughout this procedure.
The resulting 137 Cs distributions at select dates are shown in Fig. 9.This figure suggests that from April until the end of June, most of the 137 Cs stayed within the surface mixed layer, occasionally penetrating down to greater depths at select times and locations (middle and right panels).There are two distinct processes by which 137 Cs can reach below the local mixed-layer depth.First, since the mixed-layer depth is spatially non-uniform, 137 Cs can be mixed down to the bottom of the mixed layer at a remote location, where the mixed layer is deep, and then be advected laterally by currents at depth to a location with a shallower mixed layer.Second, in regions of strong downwelling, 137 Cs can be brought below the mixed layer by the strong vertical velocities, and then again advected along the deep isopycnals by the lateral currents.
For comparison, Fig. 10 shows the mixed-layer distributions of 137 Cs from the 2-D NCOM simulation, where we advect parcels using the 2-D near-surface NCOM velocities averaged over the top 25 m of the water column, and then divide the 2-D concentrations by the mixed-layer depth at that location.The mixed-layer 137 Cs fields from the 2-D and 3-D distributions are quite similar, suggesting that the effects of the vertical velocity and vertical shear on the short-term dispersal of 137 Cs are not particularly strong.
Overall, the NCOM-based 137 Cs distributions in Figs. 9 and 10 are qualitatively similar to those from the observationbased model (Fig. 5) in their shape, extent and major features.Both models predict that initially the discharged waters formed an elongated streak extending eastward from the FNPP.With time, high concentrations of 137 Cs progress further eastward into the Pacific, staying largely to the north of Kuroshio.Eventually the streaks become more convoluted and start to diffuse away.Compared to the observation-based model (Fig. 5), the streaks in the NCOM-based 137 Cs distributions (Figs. 9 and 10) are slightly wider and more diffusive (fuzzier).This is likely due to the action of rapidly varying small-scale features that are resolved by NCOM but not by the observation-based model.As we have already learned from the simulated drifters (Fig. 4), the exact position of the Kuroshio Current is displaced in NCOM simulations (Figs. 9 and 10) compared to the observation-based model (Fig. 5).Furthermore, although most of the 137 Cs remains to the north of the Kuroshio, the transport barrier associated with the Kuroshio is leakier in NCOM.

Summary and discussion
In this research, we used numerical modeling and observations from the KOK cruise to investigate the short-term  with the mean and standard deviation of 16.2 ± 1.6 PBq for the oceanic source, and 0-11 PBq with the mean and standard deviation of only 0.5 ± 2.7 PBq for the atmospheric source.These values were determined by minimizing the data-model mismatch over the KOK survey area.It is important to stress however that the KOK observations are not well-suited for constraining the atmospheric source strength because most of the fallout 137 Cs had left the survey area by mid-June.Therefore, the reported atmospheric amplitude, especially its mean value, is possibly severely underestimated.The oceanic source amplitude, on the other hand, is well-constrained by the observed KOK 137 Cs concentrations, which are due almost entirely to the oceanic discharge.
Our oceanic source amplitude is more than four times that estimated by Tsumane et al. (2011Tsumane et al. ( , 2013)), and about three times that estimated by Kawamura et al. (2011), Estournel et al. (2012) and Miyazawa et al. (2012).However, it agrees well with the estimate of 11 to 16 PBq reported by Charette et al. ( 2013), is consistent with the 14.8 PBq estimate obtained by Masumoto et al. (2012) based on the JCOPE model results, and is within the 10-34 PBq range reported by Bailly du Bois et al. (2012).As a consistency check, combining the total 137 Cs inventory of 2 PBq reported by Buesseler et al. (2011) with our model-based percentage of 137 Cs left in the KOK area from the oceanic (17.5 %) and atmospheric (4.5 %) sources, and making use of the notion that about 95 % of the measured 137 Cs inventory comes from the oceanic source, one can obtain an estimate of the total source amplitude of about 11 PBq for the oceanic and 2.2 PBq for the atmospheric sources.These numbers agree with our estimates, which were based on minimizing the data-model mismatch.Our highest atmospheric source strength is less than 50 % of the lower limit of the Stohl et al. (2012) estimate.They reported a range of 23 to 50 PBq for the atmospheric fallout.Finally, it is important to note that the deficiencies of our circulation model associated with its two-dimensional nature and coarse spatial and temporal resolution may potentially influence our estimates of the source strengths.Based on the similarity between the three-dimensional and two-dimensional NCOM runs (Sect.4), it is concluded that the neglect of vertical velocities does not cause significant distortion of the mixed-layer distributions of radionuclides in the KOK cruise area in June 2011.The lack of spatio-temporal resolution in the near-shore area, on the other hand, may cause the radioactive waters to leave the domain too soon, suggesting that our estimate may be an upper limit of the source strength.That being said, most estimates are subject to some deficiencies and our estimates provide the advantage that they present a good perspective on the range of values resulting from different models, methods and assumptions.
Several features of the measured 137 Cs field were explained through the modeling component of our investigation.First, the absence of 137 Cs at the southernmost stations was attributed to the Kuroshio Current, which acts as a transport barrier preventing the southward progression of 137 Cs.Second, the largest concentrations of 137 Cs were shown to be associated with the semi-permanent near-shore eddy that was entraining 137 Cs-rich coastal waters from the FNPP vicinity and retaining them through stirring around the eddy perimeter.This is in general agreement with Masumoto et al. (2012) who also commented on the importance of the mesoscale current structures for the spreading of Fukushima-derived radionuclides.Finally, the intermediate 137 Cs concentrations at the westernmost meridional column of stations were explained by the fact that these stations contained the more recent coolant water that had continued to leak from the reactor after the initial discharge.The pronounced non-uniform spatial structure of the measured 137 Cs fields is consistent with our modeling results, which suggest that the June distribution of 137 Cs was still streaky and patchy and was dominated by advection processes rather than by diffusion.
The 3-D effects of the vertical velocity and vertical shear were investigated through a comparison of the 137 Cs distributions from the fully 3-D NCOM velocities with those from the 2-D NCOM velocities averaged over the upper water column.The near-surface 137 Cs distributions from these two  runs are similar, suggesting that vertical shear and vertical velocity play a secondary role in the short-term spreading of Fukushima-derived radionuclides.It also suggests that one could rely on the near-surface 2-D velocities to reproduce the major features of the short-term spreading of 137 Cs.
In this research we have focused on the short-term dispersal of Fukushima-derived radionuclides, and on the assessment and interpretation of the KOK observations.The longterm spreading of the radionuclides throughout the Pacific Ocean, as well as the comparison and interpretation of other available measurements, are subjects of future study.
Although many attempts have been made to estimate the input source functions, which are of fundamental importance for assessing the aftermath of the Fukushima disaster, different methods yield different answers.As more independent estimates become available, it might be expected that the values will converge towards a single reliable number.We hope that our own investigation of the oceanic source amplitude will help to achieve this goal.

Fig. 1 .
Fig. 1.(a) Color-coded tracks of 24 KOK and 12 Mirai cruise drifters.Initial positions are marked by black and gray dots for the KOK and Murai drifters, respectively.(b) Concentration of 134 Cs in surface water.Gray shading shows the average position of the Kuroshio Current.(c) KOK-observed water column distribution of 134 Cs.Mean vertical profile of 134 Cs is shown in red.KOK area averaged mixed-layer depths in March, April and June estimated from Argo are shown by gray lines.(d) KOK station numbers and positions.

Fig. 2 .
Fig. 2. (top) The time-series of the concentration of 137 Cs in the coolant waters.Figure adapted from Buesseler et al. (2011).(bottom) Cumulative 137 Cs that had precipitated over the North Pacific by 20 April, 2011. Figure adapted from Stohl et al. (2012).Black arrows indicate mean Kuroshio Current.

Fig. 2 .
Fig. 2. (top) The time series of the concentration of 137 Cs in the coolant waters.Figure adapted from Buesseler et al. (2011).(bottom) Cumulative 137 Cs that had precipitated over the North Pacific by 20 April 2011.Figure adapted from Stohl et al. (2012).Black arrows indicate mean Kuroshio Current.

Fig. 3 .
Fig. 3. (a) Time-mean currents for the observation-based model; (b) time-mean NCOM velocities averaged over the top 25 meters; (d) standard deviation of eddy velocities for the observation-based model; (e) standard deviation of eddy velocities for the NCOM model; (c, f) differences in means and standard deviations between the two models.(c) = (a) minus (b) and (f) = (d) minus (e).The period of data for this calculation is March 14, 2011 to July 31, 2011.

Fig. 3 .
Fig. 3. (a) Time-mean currents for the observation-based model; (b) time-mean NCOM velocities averaged over the top 25 meters; (d) standard deviation of eddy velocities for the observation-based model; (e) standard deviation of eddy velocities for the NCOM model; (c, f) differences in means and standard deviations between the two models.(c) = (a) minus (b), and (f) = (d) minus (e).The period of data for this calculation is 14 March 2011 to 31 July 2011.

Fig 4 .
Fig 4. Comparison between the real and simulated drifter tracks from deployment through June 30, 2011.Panels (a), (c) and (e) on the left show the observed, the observation-based, and the NCOM-based trajectories, respectively.Panels (b), (d) and (f) on the right show the number of drifters visiting each 1x1 o bin for the observed, the observation-based and the NCOM-based drifters, respectively.

Fig. 4 .
Fig. 4. Comparison between the real and simulated drifter tracks from deployment through 30 June 2011.(a), (c) and (e) on the left show the observed, the observation-based, and the NCOM-based trajectories, respectively.(b), (d) and (f) on the right show the number of drifters visiting each 1 × 1 • bin for the observed, the observation-based and the NCOM-based drifters, respectively.

Fig. 5 .
Fig. 5. Modeled 137 Cs concentration within the mixed layer on select dates derived from the observation-based model for the oceanic source.Colors are saturated at 500 Bq/m 3 .The source domain is taken to be 20x20 km 2 square centered at the FNPP.Source strength is 16.2 PBq.

Fig. 5 .Fig. 6 .
Fig. 5. Modeled 137 Cs concentration within the mixed layer on select dates derived from the observation-based model for the oceanic source.Colors are saturated at 500 Bq m −3 .The source domain is taken to be 20 × 20 km 2 square centered at the FNPP.Source strength is 16.2 PBq.

Fig. 6 .
Fig. 6.Simulated 137 Cs concentration within the mixed layer on select dates from the atmospheric source using the observation-based model.The spatial and temporal patterns for the atmospheric source are from Stohl et al. (2012) with a source strength of 11 PBq.Green arrows indicate the modeled Kuroshio Current position.

Fig. 7 .
Fig. 7. Modeled mixed-layer 137 Cs concentration from the oceanic source (in color) at times corresponding to select stations.Black dot shows the station position.The number outlined in the lower left corner indicates station number (see Fig. 1d).

Fig. 8 .
Fig.8.Measured and simulated (using the observation-based model) estimates of137 Cs surface concentration at the times and geographical locations corresponding to the KOK stations.Here, the detection limit for 137 Cs is 1.5 Bq m -3(Buesseler et al., 2012), therefore to simplify the comparison, the simulated 137 Cs concentrations have been cut off at the same level.

Fig. 8 .
Fig. 8. Measured and simulated (using the observation-based model) estimates of137 Cs surface concentration at the times and geographical locations corresponding to the KOK stations.Here, the detection limit for 137 Cs is 1.5 Bq m −3(Buesseler et al., 2012), therefore to simplify the comparison, the simulated 137 Cs concentrations have been cut off at the same level.

Fig. 10 .
Fig. 10.Modeled 137 Cs concentration within the mixed layer on select dates from the oceanic source using 2-D NCOM velocities averaged over the top 25 meters.Dashed lines show the modeled domain.Source characteristics as in Fig. 3.

Fig. 10 .
Fig. 10.Modeled 137 Cs concentration within the mixed layer on select dates from the oceanic source using 2-D NCOM velocities averaged over the top 25 m.Dashed lines show the modeled domain.Source characteristics as in Fig. 5.

Fig. 11 .
Fig. 11.Error function Er from Eq. (1) as a function of the oceanic and atmospheric source strengths, Aatm and Aocean, for 50 cases resulting from the bootstrapping procedure of sequentially eliminating one station at a time.The position of the minimum is indicated as a small white region in each panel.Colorbars are slightly different for different subplots.To put the mean error values in perspective, the mean measured Cs concentration averaged over all stations is 415 Bq m -3 .

Fig. A1 .
Fig. A1.Error function "Er" from Eq. (1) as a function of the oceanic and atmospheric source strengths, A atm and A ocean , for 50 cases resulting from the bootstrapping procedure of sequentially eliminating one station at a time.The position of the minimum is indicated as a small white region in each panel.Color bars are slightly different for different subplots.To put the mean error values in perspective, the mean measured Cs concentration averaged over all stations is 415 Bq m −3 .