Modeling the effects of size on patch dynamics of an inert tracer

Mesoscale iron enrichment experiments have revealed that additional iron affects the phytoplankton productivity and carbon cycle. However, the role of initial size of fertilized patch in determining the patch evolution is poorly quantified due to the limited observational capability and complex of physical processes. Using a three-dimensional ocean circulation model, we simulated different sizes of inert tracer patches that were only regulated by physical circulation and diffusion. Model results showed that during the first few days since release of inert tracer, the calculated dilution rate was found to be a linear function with time, which was sensitive to the initial patch size with steeper slope for smaller size patch. After the initial phase of rapid decay, the relationship between dilution rate and time became an exponential function, which was also size dependent. Therefore, larger initial size patches can usually last longer and ultimately affect biogeochemical processes much stronger than smaller patches.


Introduction
Since 1993, Martin's iron hypothesis that iron availability limits ocean productivity (Martin, 1990), has been confirmed by 12 mesoscale iron enrichment experiments in different high nitrate-low chlorophyll (HNLC) regions, such as the equatorial Pacific, the Southern Ocean, the subarctic Pacific, etc. (Coale et al., 1996;Boyd et al., 2000;Tsuda et al., 2003;Coale et al., 2004;de Baar et al., 2005).These experiments also revealed that iron could further influence carbon cycle inside the iron patch, including air-sea CO 2 flux and vertical carbon export (Boyd et al., 2007).However, extrapolation of these results to regional and long-term impact of additional iron on biogeochemical processes is still unclear.One im-Correspondence to: F. Chai (fchai@maine.edu)portant issue is due to the limited observational capability and logistic constrains that the research vessels can only remain within the iron patch for 30 to 40 days, which is usually not long enough to observe the complete impact of iron and the entire iron patch dispersion process (Buesseler and Boyd, 2003;Boyd et al., 2007;Chai et al., 2007).For the past 12 iron enrichment experiments, while different initial sizes of the iron patch were created (e.g.64 km 2 for IronEx-1; 50 km 2 for SOIREE; 225 km 2 for Sofex), the impact of patch size on phytoplankton bloom and carbon cycling were not fully assessed.The different locations of these iron enrichment experiments make it very difficult to address how different physical processes affect movement and dilution of an iron patch.
The movement and reshaping of the fertilized patch are controlled primarily by the circulation and diffusion processes, and the chemical and biological processes inside the patch follow the movement of the patch accordingly.Inside the iron patch, the nutrients and phytoplankton bloom dynamics and carbon cycle have been addressed using both field observations and the modeling approach (Gnanadesikan et al., 2003;de Baar et al., 2005;Fujii et al., 2005;Chai et al., 2007;Boyd et al., 2007).For the physical processes controlling a surface patch, there are generally three stages of dispersion, because different control mechanisms are found on different length scale of the patch (Garrett, 1983).The first stage is a linear growth with time at a rate proportional to small-scale diffusivity; the second stage begins when the patch reaches a length scale large enough for the stirring by shear flow that advects the patch into long streaks with an exponential growth rate with time; when the length scale of the patch is longer than that of the mesoscale eddies, the third stage begins, in which surface area begins to increase linearly again.For all of the iron enrichment experiments, the iron patches usually start with a scale of ∼ 1000 m and most of the observed patches are at their second stage (Abraham et al., 2000;Ledwell et al., 1998).The third stage of the patch can be observed only if the experiment lasts long enough P. Xiu and F. Chai: Modeling patch dynamics of an inert tracer (Sundermeyer and Price, 1998).These in-situ observations mainly focused on the expanding period of a patch, however, the detailed patch dynamics of inert tracer and longterm variations are still not well understood.Plankton dynamics and carbon cycle are influenced both in the expanding and contracting periods of an iron patch, thus it is important to investigate how physical processes affect patch movement and dispersion regarding to long-term impact of iron fertilization experiments.
Owing to the observational capability and logistic difficulty of long-term measurements of a fertilized iron patch, physical-biogoechemical modeling provides us an alternative method to address some of the questions related to the impact of iron fertilization on marine ecosystems and carbon cycle (Fujii et al., 2005;Jiang and Chai, 2004;Chai et al., 2007;Jin et al., 2008).This study examines how physical processes affect surface patch dynamics and focuses on the role of the initial patch size in controlling patch dispersion.The modeled information is useful for designing and conducting large-scale iron enrichment experiments in the future.

Model description and experimental design
The model used in this study is a three-dimensional ocean circulation model, Regional Ocean Modeling System (ROMS), which represents an evolution in the family of terrain-following vertical coordinate models.ROMS solves the hydrostatic, primitive equations with horizontal curvilinear coordinates.Wang and Chao (2004) have configured the ROMS circulation model for the Pacific Ocean (45 • S to 65 • N, 99 • E to 70 • W) at 50-km resolution, with realistic geometry and topography.For this biogeochemical modeling study, we followed the approach of Wang and Chao (2004) for setting up the circulation model, but increase the horizontal resolution to 12.5 km for the entire model domain.There are 30 levels in the vertical.Near the two closed northern and southern walls, a sponge layer with a region of 5 • from the walls is applied for temperature and salinity.The treatment of the sponge layer consists of a decay term κ(T * −T ) in the temperature equation (κ(S * −S) for salinity equation), which restores the modeled variables to the observed temperature T * (salinity S * ) field at the two closed walls.The value of κ varies smoothly from 1/30 day −1 at the walls to zero at 5 degrees away from them.Liu and Chai (2009) and Chai et al. (2009) coupled the biogeochemical Carbon, Silicate, Nitrogen Ecosystem (CoSiNE) model (Chai et al., 2002) with the Pacific ROMS to study the seasonal and interannual variability of phytoplankton productivity and carbon fluxes.
Initialized with climatological temperature and salinity from the World Ocean Atlas (WOA) 2001 (Ocean Climate Laboratory National Oceanographic Data Center, 2002), the Pacific ROMS-CoSINE model has been forced with the climatological NCEP/NCAR reanalysis of air-sea fluxes (Kalnay et al., 1996) calculated using the bulk formula for several decades in order to reach quasi-equilibrium.The ROMS model is then forced with daily air-sea fluxes of heat and freshwater derived from the NCEP/NCAR reanalysis (Kalnay et al., 1996).The blended daily sea winds with resolution of 0.25 • (Zhang et al., 2006) is used to calculate the surface wind stress based on the bulk formula of Large and Pond's (1982).The heat flux is derived from the short-and long-wave radiations, sensible and latent heat fluxes that are calculated using the bulk formula with prescribed air temperature and relative humidity.The fresh-water flux is derived from the prescribed precipitation and the evaporation converted from the latent heat release.River discharges are not included in the model configuration.
To investigate patch dynamics, an inert tracer (IT) is introduced, which is controlled only by circulation and diffusion.The behavior of this IT is very similar to the SF6 tracer used during the iron fertilization experiments (Law et al., 1998;Boyd et al., 2004;Law et al., 2006).In the physical model, a nonlocal vertical-mixing scheme called K-Profile Parameterization (KPP) was used for viscosity and diffusivity (Large et al., 1994).The KPP scheme does not assume a priori that the boundary layer is well mixed and explicitly predicts an ocean boundary layer depth.The turbulent mixing is parameterized using a nonlocal bulk Richardson number within the boundary layer and through the local gradient Richardson number and a background mixing scheme below the boundary layer.Our model experiments include six different initial sizes of IT patch (Table 1) in the eastern equatorial Pacific Ocean (3 • S, 90 • W), close to the locations of the IronEx I and IronEx II experiments.For each model experiment, IT tracer was continuously added into the upper 30 m during 14 days starting from 7 January 2004.After 14 days injection of IT, first we standardized the IT concentration to a range of 0.0∼3.0µmol m −3 to mimic the iron concentration during the iron enrichment experiments (Coale et al., 1996).To define the IT patch, we use IT concentration that is higher than 0.3 µmol m −3 as a threshold, which would remove the iron limitation in the Equatorial Pacific (Gordon et al., 1997;Coale et al., 1996;Kaupp et al., 2010).By standardizing the IT concentration and mimic iron enrichment experiments, it allows us to compare results among different sizes of the IT patch.
Finally, after 14 days injection, we followed the movement of the patch to study its dilution behavior.In order to estimate the movement and reshaping of the IT patch during each time step, Gaussian ellipse fitting technique was applied to the modeled surface IT concentration.This is because ocean flow can be locally resolved into pure strain and rotation at scales greater than 1 km, which can stretch the patch into a filament.Then the distribution of IT concentration is usually well described by a Gaussian ellipse as (e.g.Abraham et al., 2000;Law et al., 2006) where C is the IT concentration (µmol m −3 ), N is the total amount of tracer within the patch.The x is the major axis of the ellipse while y is the minor axis.The σ x (t), σ y (t) are the standard deviations of the distribution in the x and y directions, respectively.After an initial adjustment, strain and dilution rates of the surface IT patch were then calculated from the changes in the dimensions of the Gaussian ellipse, by using nonlinear fitting procedure (see methods below, Abraham et al., 2000).

Results and discussion
The modeled surface patch follows a similar pattern regardless of the initial patch size, generally moves cyclonically with its center moving toward the northwest direction during the first 60 days of the simulation (Fig. 1).The lateral dispersion is largely determined by stirring associated with eddies (Garrett, 1983;Abraham et al., 2000), which changes from an initial square area to a narrow filament.For example, the initial 12.5×12.5km patch by day 38 evolves into a filament with a length of 225.2 km, a width of 126.3 km and the resultant eccentricity of 0.83; at the same time, the initial 400×400 km patch turns into a narrow filament with a length of 1865.7 km, a width of 375.4 km and a much higher eccentricity of 0.98.During the modeled experiments, smaller patches tend to have lower IT concentration and disappear much faster than larger patches (Table 1).The "disappearing day" in Table 1 is the day when IT concentration everywhere within the model domain becomes lower than 0.3 µmol m −3 , though at some days Gaussian ellipse fit could not be ap-plied mainly due to the patch being too small.12.5×12.5km patch starts with a surface area of 17 969 km 2 , reaches its maximum of 26 719 km 2 that is 48.7% larger at day 23, and disappears at day 56.The volume of 12.5×12.5km patch follows a similar manner, reaches its maximum of 1071 km 3 that is 52.6% larger than first day volume, and disappears at day 56.For surface area, the ratio of maximum to the first day value is 1.67 for large patches (mean value of 100×100, 200×200 and 400×400 patches), which is higher than 1.44 for small patches (mean value of 12.5×12.5,25×25 and 50×50 patches).While for total volume, the ratio is 1.40 for large patches, which is a little lower than 1.50 for small patches.
For patches of all sizes, both the surface area and total volume of the patch generally increases with time at first due to the high IT concentration inside the patch.Then, after reaching its peak, it starts to decrease gradually until it disappears (Fig. 2a, b).This trend is hard to obtain from in-situ measurements and most of the field experiments, which only capture the expanding period of the patch due to the limited time of observation (Timothy et al., 2006;Boyd et al., 2004).Comparison among different initial patch sizes indicates that large patches take longer time to peak than small patches.Peak time is defined as the day when the value of the surface area or total volume reaches its maximum.Surface area peak time for both the 12.5×12.5km and the 25×25 km patch is about 23 days after the start of the experiment, while for larger patches, the patch could continue expanding until more than 40 days before reaching the peak of its surface area.Disappear day is the day when the surface concentration falls below the preset limit, which can be used to describe the www.ocean-sci.net/6/413/2010/Ocean Sci., 6, 413-421, 2010 entire life history of the patch.Both the 12.5×12.5km and the 25×25 km patch can last 56 days, the 200×200 km patch persists for about 113 days, which is twice the duration of the 12.5×12.5km patch, while the 400×400 km patch lasts three times longer (194 days) than the 12.5×12.5km patch.The disappear day for surface area is the same as that for total volume, indicating the patch mostly exists in the upper water column that is consistent with iron fertilization experiments.
Previous studies showed that, in the first few days during the patch expansion, length (major axis, 2σ x ) of the fitted Gaussian ellipse is observed to increase exponentially with time in response to the strain rate (γ ), while the width (minor axis, 2σ y ) of the fitted Gaussian ellipse remains constant, where the thinning effect of strain controls the widening tendency due to diffusion (Abraham et al., 2000;Law et al., 2006).This is usually described as the second stage of the tracer dispersion suggested by Garrett (1983).In this study, the strain rate (γ ) is calculated by applying a nonlinear fitting procedure to modeled σ x (t) and time (t).During the first 32 days in this modeling study, a robust exponential function is found to fit well to the length data (Fig. 2c).Though the magnitude of ellipse length varies in response to the different initial patch sizes, the calculated strain rates are al-most the same for all patches (γ = 0.03 day −1 ).This strain rate is a little lower than the previous values of 0.07 day −1 reported for the SOIREE iron fertilization (Abraham et al., 2000) and 0.08 day −1 estimated for the sub-Arctic Pacific fertilization experiment (Law et al., 2006), indicating different dilution rates in different areas.The difference among these results might be due to both the different observational time periods and the different physical conditions related to the high spatial and temporal variability of strain in ocean currents (Sundermeyer and Price, 1998;Law et al., 2006).Law et al. (2006) suggested that different strain rate might be caused by different initial release area during the experiments, however, the strain rate is relatively stable over different initial area patches and we do not see any clear relationship between strain rate and initial area in this study.The discrepancy might be caused by the lack of biological activities in our current modeling approach, which could affect patch dynamics significantly.Ellipse width, however, decreases within this period of time, which is attributed to the stronger convergent strain rather than the small scale diffusion or mixing (Sundermeyer and Price, 1998).As a result, increasing length and decreasing width leads to an increasingly longer and narrower filament (Fig. 2c, d).After day 32, ellipse length starts to decrease and width starts to increase for about one week.It is likely caused by the unique, local physical conditions because such fluctuation has neither been observed nor modeled before.For the entire life history of the IT patch, variation of length derived from the fitted Gaussian ellipse looks much like a second-order polynomial function with time, which is in a similar manner of variation with surface area and total volume.The width of the ellipse generally decreases with time, though with a small increase around day 41 (Fig. 2d).The standard deviation, however, is much smaller for width variation (e.g.63 km for 400×400 km patch) than the length data (231 km for 400×400 km patch), which means the width remains relatively more constant than the length.For in-situ experiments, previous studies often use either linear or exponential functions to model observed variations of surface patch.From the model results, however, it is unclear whether the observed surface patch area is an exponential (Law et al., 2006), exponential and linear (Sundermeyer and Price, 1998), or polynomial function with time, since that depends on both the observational period and the local physical conditions.
To further investigate how different size patches disperse, we define and calculate the dilution rate (r) of the surface area and total volume as where A is the area or volume, t is time, and t is the time interval (Abraham et al., 2000;de Baar et al., 2005).The unit for the dilution rate r is day −1 .Overall, two stages are found for both the surface area and total volume (Fig. 3 and Fig. 4).
In the first stage, the calculated rate decreases linearly with time for about 30 days, which can be expressed as where a 1 and b 1 are the regression coefficients, r denotes the dilution rate (day −1 ) and t is the time.During the first few days, the smaller patches (12.5×12.5 km, 25×25 km and 50×50 km) of which the length scale is below the eddy scale have a similar dilution pattern, which is different from the behavior of larger patches (200×200 km and 400×400 km).Linear regression analysis shows that smaller patch tends to have steeper slope (Table 2), which results in less stable expanding relative to larger patch.
After 30 days when dilution rate reaches zero, the patch enters its second stage.Smaller patches, such as the 12.5×12.5km and the 25×25 km, start to decrease much earlier than larger patches (400×400 km).Larger patches can maintain their zero dilution rate for a very long time (e.g. over 100 days for the 400×400 km patch).For the entire second stage, an exponential function is used here to obtain the rate data: where a 2 and b 2 are the regression coefficients.In this equation, the coefficient b 2 is the one that controls the duration time in which the patch could maintain its zero dilution rate.The higher the value of b 2 , the longer the duration time.Furthermore, the value of b 2 is also shown to be size-dependent, which can be estimated well by a linear function of the initial size both for the surface area and total volume data (Fig. 5), implying that a larger initial patch has stronger ability to maintain despite the destruction due to local physical circulation and diffusion.We also conducted the same set of experiments for 2005 at the same location and the same time of the year.Different from 2004, the peak time of the surface area for the 400×400 patch is 65 days, which is longer than year 2004 (44 days).Yet, the 12.5×12.5,25×25, 50×50, 100×100 and 200×200 patches are the same as the year 2004, i.e., 23 days.On the other hand, the disappearing day for 2005 is much longer than for 2004.For example, in the case of the 12.5×12.5patch, the disappearing day is 86 days for 2005 and 56 days for 2004; for the 200×200 patch it is 155 days for 2005 and 113 days for 2004.It is difficult to generalize the development of different stages with time as described by Garrett (1983), especially using the patch measurements alone.However, the high-resolution physical models can provide detailed information about evaluation of patch dynamics, especially for different regions and/or year-to-year variations.Two similar stages of linear and exponential relationship for the dilution rate are also found for 2005, which are similar to 2004 in terms of timing in switching between these two stages.The linear relationship exists during the first 38 days, still with steeper slope for smaller initial size patch.The b 2 values from exponential regression are also found to be correlated linearly with the initial patch size (see Fig. 5 for the 2004 case), implying that the relationship between dilution rate and initial patch size is independent of the modelling time, though the magnitudes in 2005 case are slightly different from the 2004 case.The longer switching time for 2005 than 2004 implies the patches in 2004 expand much faster than 2005.As a result, patches in 2005 simulations have much longer duration in which patch maintains its zero dilution rate.The difference between these two cases is probably caused by the different physical environment, such as surface wind, current, and diffusion effect, as well as some small-scale processes like waves and turbulence that are missed by our physical model.These results again show that local physical conditions are important factors affecting patch dispersal and movement.
While there is little discussion remaining over whether iron limits production in HNLC regions, there is still significant debate on whether these results could be extrapolated to larger spatial scale and long-term effects (Boyd et al., 2007).Buesseler et al. (2008) argued that ecological impacts of iron fertilization and CO 2 mitigation are scale-dependent, and they also recommended the use of high-resolution models to assess downstream effects beyond the study area and observation period.Chai et al. (2007) suggested that for the larger patch experiments, amount of phytoplankton biomass increase stimulated by iron is much greater than the smaller patch, which allows more phytoplankton being retained within the patch.Our results show that the patch  dilution effect is more pronounced in the smaller patches than the larger ones.In addition, the larger patches tend to have more nutrients than smaller ones, which could continue to fuel the phytoplankton growth.As a result, with the same duration of iron fertilization, the larger patches could produce broader spatial extent and longer lasting phytoplankton bloom.Bowie et al. (2001) identified that patch dilution is a major pathway for the decline of added dissolved iron during their experiments.Law et al. (2006) demonstrated that dilution rate has potential impact on the onset of phytoplankton bloom termination and export.The model results indicate that physical dilution rate is not a constant but varies with time.It is largely related to the initial patch size.We speculate that initial patch size can affect patch dilution directly, which has impact on the fate of added iron in the iron enrichment experiments, could further influence the biogeochemical dynamics inside the patch.
In summary, this study shows that in order to study longterm effects of iron fertilization on phytoplankton dynamics and carbon cycle, one must consider the initial patch size as an important factor that controls the outcome of experiments.The larger size of an initial patch tends to have smaller dilution rate, which in turn increases its ability to maintain high iron concentration inside the patch.The modeled dilution rate is very sensitive to the initial size.The total volume of materials injected and local physical conditions also play a role in determining the evolution of patch dynamics and its impact on ecosystems and carbon cycling.For ocean iron fertilization experiments, besides the size of initial iron patch and total iron injected, repeatedly adding iron to the same patch of water can also maintain high iron concentration, therefore creating long lasting effects.Modeled results provide important information that can guide large-scale iron fertilization experiments.We understand that only physical processes are included in this study, however, biogeochemical processes may change the behavior of patch dilution significantly (Neufeld et al., 2002;Elliott and Chu, 2003).Detailed field measurements and long-term monitoring, including ecosystem responses to iron fertilization and the fate of fixed carbon, would improve modeling capability and allow us to assess the efficacy of sequestrating atmospheric CO 2 to the deep sea.

Fig. 1 .
Fig. 1.Temporal and spatial variations of the modeled surface IT patch (unit: µmol m −3 ), only controlled by circulation and diffusion processes.There are six different initial size patches in the experiments, and the area are 12.5×12.5,25×25, 50×50, 100×100, 200×200, and 400×400 km 2 , respectively.Red ellipse curve in each figure is the fitted Gaussian ellipse, derived from the IT concentrations.

Fig. 2 .
Fig. 2. Simulated (a) variations of total IT volume, (b) variations of surface area, (c) variations of ellipse length calculated based on the fitted Gaussian ellipse, (d) variations of the ellipse width calculated from fitted Gaussian ellipse, (e) variations of ellipse length during the first 32 days since experiment starts.Solid curves are the exponential regression results for each initial size patch.

Fig. 3 .Fig. 4 .
Fig. 3. Modeled variations of the dilution rate with time based on the modeled surface tracer area.The modeled data are separated into two stages, the first one is fitted by a linear function, and the second one is fitted by an exponential function.

Fig. 5 .
Fig. 5. Linear relationship between b 2 , derived from the exponential function fit, and the initial patch size for both the surface area and total volume.

Table 1 .
Modeled results for the six experiments conducted in this study (Maximum area is the largest area each patch can occupy, Area peak day is the day at which the patch reach its maximum area, and the Area disappear day is the day at which IT concentration becomes lower than 0.3 µmol m −3 within that region).

Table 2 .
Coefficients of the regression analysis between dilution rate and time for both surface area and total volume data.During the first 30 days, a linear equation is used, and 30 days later, an exponential equation is used (see Eqs. 3 and 4).