Cs off Fukushima Dai-ichi , Japan – model based estimates of dilution and fate

Abstract. In the aftermath of an earthquake and tsunami on 11 March 2011 radioactive 137Cs was discharged from a damaged nuclear power plant to the sea off Fukushima Dai-ichi, Japan. Here we explore its dilution and fate with a state-of-the-art global ocean general circulation model, which is eddy-resolving in the region of interest. We find apparent consistency between our simulated circulation, estimates of 137Cs discharged ranging from 0.94 p Bq (Japanese Government, 2011) to 3.5 ± 0.7 p Bq (Tsumune et al., 2012), and measurements by Japanese authorities and the power plant operator. In contrast, our simulations are apparently inconsistent with the high 27 ± 15 p Bq discharge estimate of Bailly du Bois et al. (2012). Expressed in terms of a diffusivity we diagnose, from our simulations, an initial dilution on the shelf of 60 to 100 m2 s−1. The cross-shelf diffusivity is at 500 ± 300 m2 s−1 significantly higher and variable in time as indicated by its uncertainty. Expressed as an effective residence time of surface water on the shelf, the latter estimate transfers to 43 ± 16 days. As regards the fate of 137Cs, our simulations suggest that activities up to 4 mBq l−1 prevail in the Kuroshio-Oyashio Interfrontal Zone one year after the accident. This allows for low but detectable 0.1 to 0.3 m Bq l−1 entering the North Pacific Intermediate Water before the 137Cs signal is flushed away. The latter estimates concern the direct release to the sea only.


Introduction
The Tōhoku earthquake and tsunami on 11 March 2011 triggered a series of accidents in four of the six light water, boiling reactors of the nuclear power plant Fukushima Dai-ichi on the shoreline of Honshu, Japan.In the process, extensive amounts of radioactive materials were emitted into the atmosphere where their dispersion was controlled by winds and precipitation (Stohl et al., 2011;Morino et al., 2011).Further, highly contaminated water, probably originating from desperate attempts to prevent nuclear meltdowns by injecting water into the reactor containments, was discharged or leaked into the North Pacific Ocean.This unprecedented release of radioactivity to the ocean, especially the associated fraction of 137 Cs, triggered concerns about public health and marine biota (e.g.Buesseler et al., 2011Buesseler et al., , 2012) ) which, in turn, induced the plant operator TEPCO and the Japanese regulatory authority MEXT to monitor pelagic 137 Cs activities.Despite this effort, carried out under harsh conditions, huge uncertainties remain.The total amount of 137 Cs discharged directly to the sea is uncertain within several tenfolds (Bailly du Bois et al., 2012;Japanese Government, 2011;Tsumune et al., 2012).Further, while the observations suggest that the contaminated water is initially advected southwards, subsequently turn eastwards offshore where oceanic stirring and mixing processes dilute it over a huge area (Buesseler et al., 2012), a modeling study by Qiao et al. (2011) suggests the opposite, i.e. a northeast transport confined to a narrow band.
This study adds to the discussion based on integrations of a state-of-the-art global ocean general circulation model that is eddy-resolving in the region of interest.The aims are to (1) explore which of the 137 Cs discharge estimates, combined with our simulated circulation, is consistent with 137 Cs observations, and (2) to come up with a quantitative measure of dilution or effective residence time of 137 Cs on the shelf.Such a measure could be of interest to studies to come with a focus on, e.g.marine biota off Fukushima that do not want to afford the computational expenses of ocean circulation models.Finally, (3) to explore the fate of 137 Cs released after one year.
In the following section we describe our main tool, the ocean general circulation model, and the numerical tracer release experiments.In Sect.3.1 time averaged quantities of the circulation model are compared with observations.Section 3.2 compares simulated surface currents and 137 Cs activities with actual conditions.Section 4 presents a model projection of how the 137 Cs is diluted onto basin-scale, followed by quantitative estimates of diffusion in Sect. 5. Results are discussed and summarized in Sect.6 and Sect.7, respectively.

Circulation model
We use the MOM4p0d (GFDL Modular Ocean Model v.4, Griffies et al., 2005) z-coordinate, free surface ocean general circulation model.The configuration covers the entire global ocean with an enhanced meridional and zonal resolution around Japan. Figure 1a and b show the zonally and meridionally varying resolution.The vertical grid, with a total of 59 levels is shown in Fig. 1c.The bottom topography, shown in Fig. 1d, is interpolated from the ETOPO5 dataset, a 5 min gridded elevation data set from the National Geophysical Data Center (http://www.ngdc.noaa.gov/mgg/fliers/93mgg01.html).We use partial cells.
The atmospheric forcing consists of (6-hourly) wind stress, heat, and freshwater flux fields derived from the ERA-40 reanalyses from the European Centre for Medium-Range Weather Forecasts (ECMWF) (Uppala et al., 2005, and many more to follow).In addition to the heat fluxes from the ECMWF, a flux correction restores sea surface temperatures (SSTs) with a time scale of 30 days to monthly mean SSTs derived from a blend of satellite products (C. Rathbone, personal communication, 2006).Sea surface salinity is restored to the World Ocean Atlas 2005 (Antonov et al., 2006) annual mean climatology with a timescale of 90 days.The vertical mixing of momentum and scalars is parameterized with the KPP approach of Large et al. (1994).The relevant parameters are (1) a critical bulk Richardson number of 0.3 and (2) a vertical background diffusivity and viscosity of 10 −5 m 2 s −1 .We account for double-diffusive and nonlocal fluxes.The integration started from rest with initial temperatures and salinities interpolated from the World Ocean Atlas 2005 annual mean (Locarnini et al., 2006;Antonov et al., 2006) onto the model grid.
After a spinup of 5 yr, covering the period 1993 to 1998, the clock assigning the model forcing is reset to 1993 and the model is integrated for another 5 yr.In this second period the model tracers are released as described in the next section (Sect.2.2).The period 1993 to 1998 is an arbitrary choice.Ideally, we would drive the circulation with actual, realistic fluxes.But even then, due to uncertainties in the initial conditions and the highly non-linear dynamics of ocean eddies, it would be impossible to make an exact hindcast without setting up complex data assimilation machinery.In order to account for this deficiency in our approach we explore an ensemble of tracer release experiments.The experiments differ with respect to initial conditions and the forcing.
In all other respects, not described here, the model configuration is identical to the configuration without data assimilation described by Oke et al. (2005) and used, e.g. by Dietze et al. (2009).It is furthermore, except for the grid, identical to the configuration used in Liu et al. (2010).

Tracer release
In order to simulate the accidental release of 137 Cs, we embedded (online) an artificial tracer into the MOM4p0d circulation model.The tracer is released into a surface grid-box comprising 10 × 10 km next to the location of the nuclear power plant Fukushim Dai-ichi (Fig. 2).The tracer is conservative, i.e. it does not decay but behaves like a dye, subject to mixing and advection only.On timescales much shorter than the ≈30 yr half-life of 137 Cs the behavior of our artificial tracer mimics that of 137 Cs.
The total, direct release of 137 Cs from land to sea is uncertain.The Japanese Government lists 3 major releases totaling four to five peta Bq in a report to the International Atomic Energy Agency (IAEA) (Table 1).The associated fraction of discharged 137 Cs activity is reported as 0.94 peta Bq in attachment VI-1.Most of this estimated discharge is based on visual records of the flying distance of contaminated water flushing out of a duct.
In contrast, the "Institut de Radioprotection et de sûreté nucléaire" (IRSN) states that "... 2.3×10 15 Bq of 137 Cs could have been released into the sea" (IRSN, 2011), but their underlying database and reasoning is unclear.Even higher estimates, 3.5 ± 0.7 p Bq and 27 ± 15 p Bq, are published in peer reviewed journals by Tsumune et al. (2012) and Bailly du Bois et al. (2012), respectively.They are based on an adjustment of discharge such that simulations with an eddy resolving circulation model compare well with observations (Tsumune et al., 2012), and on the temporal evolution of 137 Cs content in seawater within a 50 km area around the damaged power plant (Bailly du Bois et al., 2012).While the amount differs by almost a factor of 30, all estimates agree that most of the discharge was confined to a period of the order of a week: i.e. within the period 1 to 6 April according to Japanese Government (2011)  In our simulations we release a total of 2.3 × 10 15 Bq of 137 Cs.We do not account for air-sea fluxes of 137 Cs.Choosing the IRSN (2011) estimate is a pragmatic decision triggered by the fact that it is the median and almost exactly the mean of all lower (compared to the Bailly du Bois et al., 2012) estimates.It is not a statement concerning the creditability of the IRSN (2011) estimate.Please note that a more comprehensive discussion on the amount of 137 Cs discharged to the sea and the effect of air-sea fluxes is given in Sect.6.
As for the timing of the simulated discharge our choice is also pragmatic.As mentioned above all estimates agree in that most of the discharge was confined to a period of the order of a week sometime between 26 March and 8 April.All of our simulations feature an instant, meaning from one time step to another, discharge on 1 April.An additional experiment with all 137 Cs released instantly to the ocean one week later on 8 April 1993 yields, when compared to the 1 April 1993 discharge simulation, a measure of that uncertainty that is related to uncertainties concerning the exact time(s) of 137 Cs discharge.
In order to account for our inability to simulate an exact representation of the eddy-field we integrated an ensemble of 3 Model evaluation

Time averaged surface currents, eddy kinetic energy and surface mixed layer depth
The aim is to demonstrate that the circulation model shows reasonable realism as is to be expected from todays' eddyresolving general ocean circulation models.To that end a comparison with observed temperature and salinity does not provide much insight because (1) the thermocline does not deviate much from the initial, climatological, conditions in the short (5 to 10 yr) simulation period presented here, and (2) we restore both temperature and salinity at the surface (Sect.2.1) thereby forcing a close match to observations.More meaningful is a comparison with surface geostrophic currents as estimated from space.We use analyses based on pioneering work from Ducet and LeTraon (2001) and what have become known as the "Archiving, Validation, and Interpretation of Satellite Oceanographic data (Aviso, CNES)".Aviso maps of absolute geostrophic velocities and geostrophic velocity anomalies are available on a 1/3 × 1/3 • Mercator grid with weekly resolution.According to Niiler et al. (2003) these surface geostrophic currents show a good agreement (correlated at 0.8), with estimates based on drifters in the Kuroshio Extension region (25-42 • N, 142-180 • E) in the 1990's.In Fig. 3 the Aviso data are compared with modeled surface currents.The Ekman current is not subtracted from the modeled currents.This is a minor issue because the Ekman currents are typically much smaller than the eddy currents, and the Kuroshio is essentially in geostrophic balance (Uchida et al., 1998).Hence, we conclude that most of the misfit between model, and observations offshore is due to deficient model dynamics.However, this might not apply near the coast where the assumption of geostrophy might be rendered invalid by bottom friction and degraded accuracy of Aviso data.
That said, the agreement between model and observation is remarkable in the region of interest off Fukushima.Namely the main features described by Niiler et al. (2003) are in good agreement: the Kuroshio Extension jet is deflected northwards after crossing the Izu Ridge (Fig. 1) and at a bottom trough at 35 • N, 150 • E. The Tsugaru current flowing eastwards between Honshu and Hokkaido, divides into two branches.One branch feeds a near-coastal current flowing southwards to about 36 • N. The other branch merges with the northeast oriented Oyashio frontal current.
The recirculation cells of the Kuroshio south of 36 • N are, on the other hand, not well simulated by the model.But since (as will be shown later) the released 137 Cs does not cross the Kuroshio at the surface in our simulations, this misfit should be irrelevant.
Figure 4 shows a comparison between simulated and observed (as derived from Aviso geostrophic velocity anomalies) eddy kinetic energy (EKE), which characterizes the mesoscale activity.It is defined by where u and v are deviations of the zonal and meridional surface velocities from their time mean (one year in this case), respectively.The overbar denotes a temporal average over one year, u and v are weekly snapshots (for both model and observations).The simulated EKE is similar to the observed one in the region of interest (east of 140 • E) although, overall, biased high.Typically the simulated EKE, integrated over a box bounded by 40-159 • E and 34.5-42 • N is overestimated by 20 to 70 % in the simulations.Such overestimation has been identified already in a similar study (Liu et al., 2010).Most straightforward is to argue that this misfit between the model and observations is due to inaccuracies in the data.In fact it is known that EKE estimates based on observations from space are biased low (e.g.Fratantoni, 2001).A comprehensive analysis is beyond the scope of this paper.Nevertheless, we remind that, historically, friction has been kept at the minimal level demanded by numerical stability in eddy-resolving models.In our model we apply a state-dependent horizontal Smagorinsky viscosity scheme (Griffies and Hallberg, 2000).A main advantage of this scheme is that friction, on average, can be reduced relative to traditional schemes without violating numerical stability.Hence it can not be ruled out that part of the misfit between observed and simulated EKE is caused by missing friction in our model.Further, we conclude that the misfit between observed and simulated EKE is greater than the observed interannual EKE variability.Hence we can not decide, based on EKE, which of the nominal simulation years is most representative for the actual conditions in 2011.In addition to horizontal surface currents, vertical mixing affects 137 Cs activities at the surface.On the monthly to seasonal timescales considered in this study, vertical mixing is dominated by surface mixed layer dynamics.A direct comparison of surface mixed layer depths is complicated by its extremely high variability both in space and time.This variability is, predominantly, a result of the strong eddy activity (compare Fig. 4): the formation and intensification of anticyclonic eddies is accompanied by downwelling fed by the convergence of horizontal surface currents, which deepen the surface mixed layer.For cyclonic eddies, on the other hand, the reverse holds: formation and intensification comes along with a divergence of horizontal surface currents which pulls up the isopycnals and results in shallower surface mixed layers.Figure 5 compares a regional average of modeled surface mixed layers with data obtained from profiling Argo floats; given the high variability in both model and observations, the model is consistent with the observations.

Surface currents at the time of the accident and 137 Cs activities on the shelf
We aim to simulate the dilution and cross-shelf transport of 137 Cs discharged off Fukushima Dai-ichi.This task is complicated by the temporal variability of surface currents on spatial scales including that of the size of the shelf (≈40 km).
It is straightforward to assume that this variability is considerable for the task at hand, given that the meso-scale activity in the region ranks among the highest world-wide.Hence, the consistency between simulated and observed temporal or spacial averages demonstrated in Sect.3.1 does not necessitate a realistic simulation at a given time.Or, in other words, a simulated state of the eddying circulation similar to actual conditions is a precondition for successfully simulating the accidental 137 Cs discharge in April 2011.This, fortunately, turned out not to be as hopeless as the combination of chaotic eddy dynamics and the fact that we do not assimilate realtime data suggests.
The reason why we are able to come up with a rather realistic simulated state of the circulation (despite vivid eddy dynamics) is that, apparently, only two contrasting states of the circulation prevail in the region.Both, visual inspection of sea surface height data (Aviso, CNES) and our simulations imply that there is a state where the average movement of surface water on the shelf is southbound.This state seems to be controlled by the southbound near-coastal branch of the Tsugaru current, and prevailed also at the time of the accident and in the nominal model years 1993 and 1996 (upper row of panels in Fig. 6).
The other state of circulation is showcased by observations on 12 March 1997 (lower left panel in Fig. 6) where the average movement of surface water on the shelf is northbound, apparently effected by an anti-cyclonic eddy shed off by the Kuroshio current.Similar patterns, i.e. an average northward surface velocity on the shelf effected by an instability or a meander originating from the Kuroshio is simulated in the nominal model years 1994, 1995 and 1997.We conclude that the simulations labeled with the nominal model years 1993 and 1996 are well suited to simulate the fate of 137 Cs discharged in April 2011 since they mimic the actual conditions, namely the average southbound movement of surface water on the shelf.
A comparison of simulated and observed 137 Cs concentrations (Fig. 7) at sampling stations indicated in Fig. 2 confirms that, indeed, the nominal model years 1993 and 1996 mimic actual conditions relatively well while the other nominal years do less so.The measurements shown in Fig. 7 were Three simulations (the nominal years 1993, 1996 and the additional 1993 simulation where the discharge is delayed by one week) show reasonable agreement with the observations at all stations in Fig. 7 except at those nearest to the actual release site (station 1, 2).This misfit, however, does not come as a surprise since our first "wet" model grid-cell is situated further offshore in deeper waters than the actual release site (Fig. 2).We speculate that this biases the simulated activities low.Note also that the results are relatively insensitive towards the exact time of the 137 Cs discharge because differences between the standard and the 1993 simulation with delayed discharge are small (black and grey line in Fig. 7).
All remaining simulations (nominal years 1994, 1995 and 1997) are biased substantially low in Fig. 7.This is consistent with the sea surface height analysis in Fig. 6, which indicated that those nominal years are representative of a differing (although also realistic) state of the circulation.We speculate that all nominal years, taken together, encompass the range of the variable cross-shelf transport.The latter assumption will be used in Sect. 5 to derive a quantitative estimate of the cross-shelf transport variability.
A comparison with MEXT measurements (described in Buesseler et al., 2011) further offshore (although still on the shelf) is less instructive concerning the realism of respective simulations.Most of the measurements are below the detection limit of 10 Bq l −1 and those exceeding the detection limit feature high scatter, which is generally not reproduced in neither of the simulations (Fig. 8).We conclude, nevertheless, that the 1993, 1996 and the additional 1993 simulation where the discharge is delayed by one week, are consistent with the MEXT measurements, i.e. that in general, they simulate the right order of magnitude of 137 Cs surface activities.

The fate of 137 Cs released directly from the land to the ocean
The vertical integral, from the seafloor to the surface, of the 137 Cs activities yields a two-dimensional, horizontal distribution.A measure of the average horizontal position of the 137 Cs patch is then given by the center of mass of the vertically integrated activities.This metric renders the comparison of various simulated 137 Cs patches in one figure possible.Figure 9 shows that the "realistic" simulations (nominal years 1993 and 1996) feature a center of mass that is initially southbound in accordance with what is expected from the state of the circulation as depicted in Fig. 6 and discussed in Sect.3.2.Note that Fig. 9 carries also information about the simulated temporal variability of the cross-shelf transport: it takes the center of mass of the 1996 release more than 10 weeks to cross the shelf break (200 m isobath) south of 36 • N 40 .In contrast, the center of mass of the 1997 release leaves the shelf way quicker, after less than a month, north of 37 • N 20 .
The further development of the "realistically" simulated (i.e.nominal years 1993 and 1996) 137 Cs patches is shown in Fig. 10.Three processes acting on 137 Cs leaving the shelf are apparent: (1) the Kuroshio accelerating 137 Cs eastwards into the basin; (2) mesoscale activity that disperses and mixes the 137 Cs meridionally into the Kuroshio-Oyashio Interfrontal Zone where cold low salinity Oyashio water mixes with warmer high salinity Kuroshio water (Fig. 11); and (3) the near-coastal southward branch of the Tsugaru current that, together with the Oyashio current, flushes the east coasts of Honshu and Hokkaido north of 36 • N with uncontaminated waters originating from the Sea of Japan and the subpolar gyre, respectively.
Figure 10 also shows that even our "realistic" simulations differ considerably from one another.This holds especially for the appearance of enhanced (i.e. more than 10 mBq l −1 ) activities offshore east of 145 • E. But even though details among the "realistic" simulations differ, common to all of them is that, ultimately, enhanced surface 137 Cs activities appear in the Kuroshio-Oyashio Interfrontal Zone, a formation site of North Pacific Intermediate Water (NPIW) (Yasuda et al., 1995).Further, the simulations indicate that 137 Cs activities up to 4 mBq l −1 prevail for more than a year in the region (Fig. 12).This allows for enough time for the 137 Cs signal to enter the NPIW before it is flushed away.Hence, 137 Cs discharged directly from the land to the sea off Fukushima may be a useful water mass tracer if its penetration into the NPIW can be measured in the years to come.Figure 13 provides a model-based estimate of the 137 Cs concentrations in the NPIW, defined by its salinity minimum between 300 and 800 m, one year after the accident.Note that according to Park et al. (2008), they are above (although close to) the minimum detectable activity (0.1 and 0.33 mBq l −1 for deep and surface water, respectively).

Quantitative estimates of dilution
For studies to come, e.g. with a focus on the effect of radiation on marine pelagic biota, a quantitative measure of mixing on and across the shelf break off Fukushima, might be of interest because such a measure could replace the integration of expensive eddy-resolving ocean circulation models.The concept of "Fickian diffusion" provides such a measure.It is used in numerous analyses of the mixing and stirring of a purposely (or accidentally) released substance, which can be measured in very low concentrations due to its radiative, fluorometric or electrophilic nature (e.g.Houghton et al., 2006;Sundermeyer and Price, 1998;Ledwell et al., 1998;Torgersen et al., 1997;Okubo, 1971).The caveat to bear in mind, however, is that oceanic estimates of "Fickian diffusion" are dependent on the spatial scale under consideration (e.g.Okubo, 1971;Ledwell et al., 1998 "Obs. 7.4.2011" and"Obs. 12.3.1997")shows geostrophic velocities as estimated from space (Aviso, CNES).The subsequent panels titled "Mod...." refer to simulated currents.The upper and lower row of panels feature typical, contrasting states of the circulation, respectively (see text).

Horizontal diffusivity on the shelf
In the following, we will diagnose effective (model) diffusivities in our 137 Cs release simulations.The 137 Cs tracer is instantaneously released into one grid box.Intuitively, it is clear that the temporal evolution of both the maximum concentration and the effective area occupied by the tracer carries the information to diagnose effective diffusivities.We find that the maximum concentration in the patch is not a good measure of diffusion in our model, because the steep gradients caused by the instantaneous point source give rise to numerical dispersion, until the gradients are smoothed out to some degree.As for the approach using a measure of the effective area, we follow Orre et al. (2008).We omit their derivation for the sake of brevity and merely document the calculations.The 137 Cs tracer C(x, z, t) is vertically integrated over the whole water column: (x, t) = C(x, z, t)dz. (2) x(x, y) and z denote the horizontal and vertical position, respectively.t is time.The vertical integral (x, t) is normalized with the total (i.e.integrated over the entire model domain) tracer content: A measure of the effective area occupied by the tracer is defined by the second (physical space) moment of the vertically integrated and normalized tracer field: where x com (t) is the (horizontal) position of the center of mass of (x, t).An estimate of the effective diffusivity K is then given by Figure 14 shows the temporal evolution of diffusivities diagnosed from the five 137 Cs release experiments (named 1993 to 1997).We distinguish three phases: (1) initial, relatively high (200 to 300 m 2 s −1 ) diffusivities; (2) a saddle phase (60 to 100 m 2 s −1 ); and (3) the final increase which, in the 1993, 1994 and 1997 releases, exceeds values of 2000 m 2 s −1 .We speculate that the high initial values, which are relatively similar among the ensemble members, are a consequence of numerical dispersion.With time (order of days), concentration gradients smooth out and the associated numerical dispersion is suppressed as we enter the saddle phase.The length of the phase is determined by the average movement of the tracer patch.As soon as a significant fraction of the tracer crosses the shelf-break and is gripped by the coastal southward branch of the Tsugaru current, the Kuroshio, or the vivid eddies offshore, the diffusivity increases again.This explains why the saddle phase of the 1997 release, which reaches the shelf-break quickly (Fig. 9), is so much shorter than those of the other releases.We conclude that the minimum values from 60 to 100 m 2 s −1 are representative for the diffusivity on the shelf.

Cross-shelf diffusivity
To obtain an estimate parameterizing, the cross-shelf exchange as a Fickian diffusivity, the previous analysis is not helpful because it is unclear at what point in time (if at all) the above diagnosed diffusivity is representative of the crossshelf exchange.A rough estimate can be obtained by applyequation ).The station numbers are assigned to locations in Fig. 2. The black, green, magenta, brown, orange lines correspond to simulated nominal years 1993 to 1997, respectively.The grey lines correspond also to simulated year 1993 but with 137 Cs depositions being applied one week later.Note, that at station [4], day 14 a reading of 186 Bq l −1 is not shown.
which is solved by  with the initial concentration C 0 and The timescale τ describes the rate of decrease of the total 137 Cs content on the shelf.It can be interpreted as an effective residence time of 137 Cs on the shelf as determined by the circulation.Figure 15 shows the decrease of simulated 137 Cs content on the shelf for the respective ensemble members: the content stays rather constant for several weeks, followed by a rapid decrease.The timing of the rapid decrease is correlated with the center of mass approaching the shelf break (Fig. 9).Least square fitting of Eq. ( 8) applied to the first two months of rapid decrease shows reasonable skill (Fig. 15).The residence timescale τ , averaged over all ensemble members is 43 ± 16 days (Table 2).Assuming a shelf width, x, of 40 km, this timescale corresponds to a horizontal diffusivity of 500 ± 300 m 2 s −1 .We interpret the uncertainty of ±300 m 2 s −1 in the latter estimate as a measure of temporal variability.

Discussion
The rather realistic circulation is a real asset of our model (compare Sect. 3), especially because our simulations do not, in contrast to other studies (e.g.Qiao et al., 2011), feature the long-time and almost endemic model problem of a Kuroshio separating from the coast too far north (Fig. 3).The agreement between observed and simulated 137 Cs activities (Figs. 7 and 8), however, might be misleading, either because we released an unrealistic amount of 137 Cs or because our model misses a major mechanism removing 137 Cs activity from the surface ocean.To this end we discuss in the following so far unaccounted effects of (1) 137 Cs present in seawater already prior to the accident, (2) air-sea deposition, (3) uncertainties of the direct land-sea discharge and (4) transfer of 137 Cs out of the water column into the sediment by 137 Cs adsorbed to settling biotic or abiotic particles.
1. Nakanishi et al. (2010) present 137 Cs readings up to 4 mBq l −1 in the waters of interest prior to the accident.This offset has no significant effect on the agreement between observed and simulated 137 Cs activities in Figs.7  and 8.In contrast, it is a significant contribution to the simulated activities in Figs. 10, 12 and 13.
2. Stohl et al. (2011) and Morino et al. (2011) estimate that the accumulated atmospheric 137 Cs deposition to the ocean peak at 50 to 200 k Bq m −2 .Diluted over a surface mixed layer depth of at least 20 m (Fig. 5), this may, locally and at times, raise surface activities by up to 10 Bq l −1 before the signal is diluted.While this seems to add considerable uncertainty to the model evaluation in Figs.7 and 8, we expect, however, that the average effect of unaccounted air-sea fluxes is minor because oceanic 131 I to 137 Cs ratios suggest that most of the oceanic 137 Cs is the result of a direct discharge from the land to the ocean rather than resulting from atmospheric deposition (Tsumune et al., 2012).2012) assigns activities upstream, which rank among the highest in the whole region (their Fig. 4a).This biases their early inventories high.Second, their late (i.e.mid April till July) inventories are, on the other hand, biased low because of the detection limit inherent to the respective 137 Cs measurements (10 to 50 Bq l −1 , Buesseler et al., 2011).The detection limit may well transfer to a substantial inventory underestimation of 1 to 7 p Bq (in the 60 × 120 km box assuming a surface mixed layer depth of 20 m) but is not accounted for in the analysis.We conclude that the combined effect of "early" inventories biased high and "late" inventories biased low is an overestimation of the slope of the logarithmic fit that, by extrapolation backwards in time, biases the discharge estimate of Bailly du Bois et al. (2012) high.Third, their extrapolation backwards in time implicitly assumes a constant transport of 137 Cs out of the box.This does not hold for the initial period after the discharge because, before leaving the box, the 137 Cs must first reach the box boundary.
4. Concerning the transfer of 137 Cs to the sediment, lessons learned in the Baltic Sea after the Chernobyl ac-cident on 26 April 1986 give some guidance.As summarized in a short review in the associated paper in Ocean Science Discussions (Dietze and Kriest, 2011) (but discarded here for the sake of brevity), initially (i.e. on time scales of months to years), the redistribution of 137 Cs should be dominated by physical processes such as mixing and advection.
In summary, our study seems to be consistent with available data and information (an exception being the Bailly du Bois et al. ( 2012) discharge estimate), but huge uncertainties remain.Further, it cannot be used to produce any advise on hazards of radiative contaminants in general as it treats 137 Cs activities only.To this end it is noteworthy that elements like Pu are harder to detect at harmful concentrations.

Conclusions
We set out to quantify the mixing off Fukushima, Japan between surface waters on the shelf and those offshore by simulating 137 Cs discharged in the aftermath of an accident at the nuclear power plant Dai-ichi on 11 March 2011.Our main tool was a global circulation model with an enhanced, eddyresolving horizontal resolution in the region of interest.A comparison with sea surface height as observed from space and Argo float data showed that the model reasonably simulates actual conditions.This is supported by the model's capability of reproducing observed 137 Cs activities.
A caveat is the high uncertainty of the total 137 Cs discharged directly from the land to the ocean.Obviously, if modeled 137 Cs activities agree well with observation even though the simulated release is inconsistent with actual conditions, one ended up with getting the right answer for the wrong reason.Or in other words, as both the amount discharged and the circulation are uncertain to some degree, it is problematic to infer one from the other.What remains is to merely state that (1) we released 2.3 p Bq in our simulatrions following an estimate of IRSN (2011).This is, incidentally, almost exactly the average of the 0.94 p Bq estimate by the Japanese Government (2011) and the 3.5 ± 0.7 p Bq estimate by Tsumune et al. (2012).(2) Our simulations are apparently inconsistent with the high 27 ± 15 p Bq estimate of Bailly du Bois et al. (2012).(3) We think, as outlined in Sect.6, that the latter estimate is biased high because of methodological issues.
That said, we diagnose an effective horizontal diffusivity of 60 to 100 m 2 s −1 on the shelf via a measure of the spread of the simulated 137 Cs patch with time.This estimate is consistent with observations of sea surface temperature variability on the Newfoundland shelf (Zhai and Greatbatch, 2006) and a model study by Capet et al. (2008) on the Argentinian shelf.Further, we diagnose a residence time, as effected by the simulated circulation, of 43 ± 26 days of 137 Cs on the shelf, based on the exponential decline of the simulated 137 Cs content with time.A Fickian down-gradient parameterization of the associated cross-shelf transport would necessitate a diffusivity of 500 ± 300 m 2 s −1 if applied on spatial scales similar to the actual shelf-width (40 km).The respective uncertainties are associated with the (simulated) temporal variability of the cross-shelf transport.
With respect to the fate of 137 Cs discharged, our simulations project enhanced surface activities offshore in the Kuroshio-Oyashio Interfrontal Zone which is known to be a formation site of North Pacific Intermediate Water (NPIW).Specifically, we find 10 to 50 mBq l −1 six months after the release and an order of magnitude less after one year.This allows for enough time for a 137 Cs signal entering the NPIW.It is, however, close to the minimum detectable activity (MDA) in deep sea water (0.1 mBq l −1 , Park et al., 2008).Note that the latter estimate refers only to the effect of the direct release from the land to the ocean.Air-sea deposition of 137 Cs and background concentrations prior to the accident are not accounted for.

Fig. 1 .
Fig. 1.Model grid.Panels (a) and (b) are horizontal and meridional resolution, respectively.Panel (c) shows the vertical resolution z as a function of depth.Panel (d) shows the model bathymetry in the region of interest (white areas are considered as land by the model while the grey patches outline the actual land distribution).

Fig. 3 .
Fig. 3. Surface currents, averaged over the period 1994 to 1996.The color shading denotes absolute velocity.The arrows indicate the direction of the currents.The magenta line is the 200 m isobath.Panel (a) shows simulated velocities.The actual model resolution is three times higher (in both meridional and zonal direction) than indicated by the arrows.Panel (b) shows geostrophic velocities as estimated from space (Aviso, CNES).

Fig. 4 .
Fig. 4. Eddy kinetic energy as defined by Eq. (1).The left row of panels refers to simulated nominal years indicated in the respective upper left corners.The right row refers to estimates based on geostrophic velocity anomalies as observed from space.Due to lagged data access, the 2011 observations refer to an average from September 2010 to September 2011.

Fig. 5 .
Fig.5.Surface mixed layer depth (defined as the depth where density σ 0 exceeds surface values by 0.125) averaged over the region bounded by 143 • E to 160 • E and 33 • N to 42 • N. The red crosses are calculated from a total of 1796 Argo profiles in the period 1998 to 2004.The vertical red lines are associated standard deviations calculated from all observations found in the region at a given month.The black line is calculated from simulated, daily snapshots of model year 1993.The vertical grey lines denote the associated standard deviation which represents the simulated spatial variability at a given day.The gap in November is caused by files corrupted during integration.The dashed green, magenta, brown, orange lines correspond to simulated nominal years 1994 to 1997, respectively.

Fig. 6 .
Fig.6.Surface velocities at dates indicated by respective panel titles.The color shading denotes absolute velocity.The arrows indicate the direction of the currents.The magenta line is the 200 m isobath, the red cross marks the nuclear power plant.The first column of panels(titled  "Obs.7.4.2011"and "Obs.12.3.1997")shows geostrophic velocities as estimated from space (Aviso, CNES).The subsequent panels titled "Mod...." refer to simulated currents.The upper and lower row of panels feature typical, contrasting states of the circulation, respectively (see text).

Fig. 7 .
Fig. 7. Observed (red and blue o) by TEPCO and simulated (colored lines) time series of surface 137 Cs activities at different stations as indicated by panel titles.The station numbers are assigned to locations in Fig. 2. The black, green, magenta, brown, orange lines correspond to simulated nominal years 1993 to 1997, respectively.The grey lines correspond also to simulated year 1993 but with 137 Cs depositions being applied one week later.The blue o in the upper right panel denote measurements at stations 1 & 2, close to the discharge point.

(Fig. 8 .
Fig. 8. Observed (red o and vertical bars) by MEXT and simulated (colored lines) time series of surface 137 Cs activities at different stations as indicated by panel titles.The red vertical bars denote samples hosting 137 Cs activities below the detection limit (≈ 10 Bq l −1).The station numbers are assigned to locations in Fig.2.The black, green, magenta, brown, orange lines correspond to simulated nominal years 1993 to 1997, respectively.The grey lines correspond also to simulated year 1993 but with 137 Cs depositions being applied one week later.Note, that at station [4], day 14 a reading of 186 Bq l −1 is not shown.

Fig. 9 .
Fig. 9. Movement of the center of mass of simulated 137 Cs patches.The line colors denote respective nominal years indicated in the legend.The dots on the respective tracks are separated by one fortnight.The cyan-to-blue color shading refers to the bathymetry interpolated onto the model-grid in meters.

Fig. 10 .Fig. 11 .
Fig. 10.Simulated surface 137 Cs activities.Shaded in color is that fraction that is effected by the direct land-sea discharge.Background activities present prior to the accident and effects of air-sea deposition are not accounted for.The upper and lower row of panels refer to simulations labeled by nominal years 1996 and 1993, respectively (also indicated in the upper right corner of each panel).The time information centered in each column of panels denote days elapsed since the 137 Cs discharge on 1 April of respective years.The 200 m isobath is indicated by white contours.

Fig. 12 .
Fig. 12. Simulated (nominal year 1996 simulation) surface 137 Cs concentration in units mBq l −1 on 1 April, one year after the 137 Cs release off Fukushima Dai-ichi (indicated by the black x).Only activities, which according to Park et al. (2008), exceed the minimum detectable activity in surface waters, are shown.

Fig. 14 .
Fig. 14.Simulated diffusivity diagnosed (with Eq. 5) from the effective area covered by the released 137 Cs patch.Day zero corresponds to the time of release in the respective nominal years indicated by the figure legend.

Fig. 15 .
Fig. 15.Decline of the simulated 137 Cs inventories on the shelf (dashed lines).The time axis corresponds to days elapsed since the 137 Cs was released in the nominal years indicated by the figure legend.The solid lines are least square fits to Eq. (8).

Table 1 .
Government (2011)active materials discharged to the sea at Fukushima Dai-ichi nuclear power station according to JapaneseGovernment (2011).
3. As pointed out already in Sect.2.2, the amount of137Cs discharged into the ocean is uncertain.Our simulations apply the IRSN (2011) estimate of 2.3 p Bq and (relatively insensitive to the timing of the discharge) are in reasonable agreement with observed activities as shown in Figs.7 and 8.We realize that this assessment is subjective.The same would, however, apply to more quantitative estimates of model-data misfit since they must contain some weighting of observations, and this weighting, inescapably, is also subjective to some extent.To summarize, we conclude based on Figs. 7 and 8 and the fact that simulated activities scale linearly with the total discharge that our simulations are consistent with observations if the total discharge is within the envelope of the Japanese Government (2011) (0.94 p Bq) and Tsumune et al. (2012) (3.5 ± 0.7 p Bq).Our simulations are, however, apparently inconsistent with the high 27 ± 15 p Bq estimate of Bailly du Bois et al. (2012) which we think is biased high: Bailly du Bois et al. (2012) extrapolate 137 Cs measurement in space and time to calculate weekly inventories contained by a 60 × 120 km box.Because the data is reasonably described by a logarithmic fit, they conclude that an extrapolation backwards in time by one week to the anticipated time of discharge is admissible.