Auxiliary Methods The BEC model was described by Moore et al. (2002; 2004). Here we review only some key aspects of the model relevant for the present study. The BEC model includes several phytoplankton functional groups including diatoms, calcifiers, diazotrophs, and picoplankton (Moore et al., 2002: 2004). The calcifiers are included implicitly as a variable fraction of our small phytoplankton group, according to environmental factors (Moore et al., 2004). These functional groups are believed to be important drivers of ocean biogeochemical cycles and are included specifically to promote the study of ocean biogeochemistry - climate interactions (Le Quéré et al., 2005; Hood et al., 2006). Our diazotroph functional group is parameterized largely based on observations of Trichodesmium spp. (Moore et al., 2002; 2004). To the extent that other diazotrophs may respond differently to environmental forcings, the model may not capture their dynamics well. Very little is known at present about how non-Trichodesmium species respond to environmental forcings. Several studies have modeled N fixation rates but did not include both Fe and P as potentially limiting factors (Hood et al., 2001; 2004; Fennel et al., 2002; Coles et al., 2004). We assume fixed C/N/P ratios for each phytoplankton group and for the sinking export of particulate organic matter (POM). These ratios are set close to the Redfield values (117C/16N/1P, Anderson and Sarmiento, 1994), with the exception of the diazotrophs that have an N/P ratio of 50 mol/mol, and the same C/N ratio as the other groups. We note that the N/P ratio of Trichodesmium spp. can vary considerably in different environments (Krauk et al., 2006), a feature not currently included in the model. Sinking particles are implicitly treated and assumed to sink and remineralize instantly within the water column at the same grid location where they are generated. The remineralization profiles are determined by a variant of the mineral ballast model proposed by Armstrong et al. (2002). In this model the sinking flux of mineral ballast (biogenic silica, calcium carbonate, and mineral dust particles) is presumed to transport some POM to the deep ocean which is remineralized with the same profiles as the ballast material (see Moore et al., 2004 for details). All nutrients are simulated globally over all depths with no restoring of any kind (Moore et al., 2004). The parameter changes noted in Table S1 are similar to those described previously by Moore et al. (2006). The most significant effect is to increase iron scavenging in the deep ocean, while retaining similar rates in the upper ocean. The original iron related parameters worked well for the upper ocean but allowed the deep ocean to drift up to 0.6nM in most regions over longer simulations. Other changes make the diazotrophs somewhat more efficient at nutrient utilization. We have also modified the parameters associated with the remineralization profiles of soft POM and biogenic silica. These values gave a modestly better match to upper ocean observations of macronutrient concentrations. Auxiliary Results and Discussion There are very limited field observations of nitrogen fixation over seasonal timescales to compare our results with for model evaluation. Perhaps the best studied sites would be at the time series stations, BATS and HOT, where field estimates for N fixation are ~ 15 mmolN/m2/yr at BATS (Orcutt et al., 2001) and range between ~31-51 mmolN/m2/yr for HOT (Karl et al., 1997). We noted previously our low simulated values in the North Atlantic subtropical gyre. Our simulations are somewhat low at BATS due to extreme P limitation in this region ~ 3-6 mmolN/m2/yr. At HOT our control simulation was rather low at 12 mmolN/m2/yr (not surprising given the constraint on denitrification). For the simulations where the denitrification constraint was released we simulated values of 25, 74, and 37 mmolN/m2/yr for cases 1, 5 (LGM dust), and 6, respectively during year 200. Table S2 summarizes the basin scale N fixation and denitrification in each of our experiments as well as the N cycle impacts on primary production, export production, and atmospheric CO2 concentrations. Our simulations can be compared with two other global scale estimates of N fixation. Lee et al. (2002) estimated N fixation based on DIC drawdown in warm, nitrate-depleted regions, and suggested highest rates of areal N fixation in the North Atlantic, with peaks in each basin in the subtropical gyres. Deutsch et al. (in press) estimated the highest rates within the Pacific basin with maximum rates in the eastern tropical Pacific. Our own simulations had highest rates in the northern Indian basin, where high dust flux and substantial denitrification generated favorable conditions. At present there are very few observations in the Indian Ocean and over most of the Pacific basin. The best sampled region by shipboard studies is in the North Atlantic, where Capone et al., (2005) estimated a mean N fixation by Trichodesmium of ~ 87 mmolN/m2/yr. Our simulations for cases 1 and 6 give similar or slightly higher values for the tropical North Atlantic, but much lower values farther north (Figure 4). Preliminary work suggests that allowing the diazotrophs to access dissolved organic phosphorus improves our simulations in the North Atlantic, but has little impact elsewhere, and thus would not change our main conclusions here. Cases 1 and 6 also match the observed surface and sub-surface N* values fairly well in the Atlantic basin (Figures 5 and S4). Negative feedbacks that act to stabilize the ocean N cycle arise via different paths, which we broadly characterize into two groups: uncoupled and coupled. First, there are mechanisms whereby a perturbation to an individual N source/sink will be self-damped; for example in Case 1 (Figure 2), Pacific denitrification decreases with time following the initial pulse independent of any changes in N fixation in the basin simply because overlying surface nitrate concentrations, and thus export production, are lowered. Both coupled and uncoupled feedbacks are illustrated in the Indian basin in Case 1 whereby the rapid rise in N fixation balances the enhanced denitrification after about 100 years. As discussed earlier, the coupled feedbacks can occur on a range of timescales from rather short for surface water N/P to longer for changes in ocean N inventory. We can quantify the strength of negative uncoupled feedbacks (Ftotal) in our simulations by comparing peak imbalance (Imax in Tg N/yr) between N fixation and denitrification with the ending imbalance (Iend), where: Ftotal would be 1.0 with a perfect negative feedback that returned the N-cycle to net zero balance; an ocean with no stabilizing feedbacks would have a value of 0.0. We subtract 6.7 TgN/yr from the maximum and ending imbalances to account for the initial imbalance in the control simulation. In Case 1 the negative feedback strength was 0.80 at year 200, increasing to 0.86 at year 430. Case 6 had a similar value (0.81 at year 200). In the Case 5, more Fe-replete experiment, the feedback strength was 0.96, and the initial imbalance was larger (Figure 3). The feedback strength was modest in cases 2 and 3 (0.33 and 0.40, respectively). The feedback response was always weaker in the Pacific, and globally the gap between N fixation and denitrification was always still narrowing at the end of the simulations. Case 4 was unusual due to the artificial limit on denitrification, with a calculated feedback strength of 1.2. The rise in denitrification nearly matched the elevated N fixation, erasing the initial imbalance. The coupled feedback response (Fcoupled) can be quantified in an integrated sense over our 200 year simulations by comparing the net change in nitrogen inventory (with both denitrification and N fixation responses) to the net change due solely to the initially perturbed component, assuming that the other component of the N cycle remained fixed. For Case 1 for example: where and (Tg N) are the perturbations in the global time-integrated denitrification N loss and N fixation gain, respectively. Thus, if N fixation did not respond in any way to the increase in denitrification, the coupled feedback response Fcoupled would be 0.0, and a perfect, inventory maintaining feedback response would have value 1.0. For Case 1, the integrated feedback response is 0.25, indicating that if N fixation had not increased, the oceans would have lost an additional 25% more N over two centuries. For Case 2, the coupled feedback response is 0.16, and for Case 3 (reversing and in the above equation) the feedback response is 0.33, and for Case 4 the coupled feedback response is 0.26, and for Case 6 it was 0.29. Note that the strength of the integrated feedbacks are weaker than at the end of simulation (Ftotal above). This is due to the timescales of the negative feedbacks, which are not instantaneous but take several decades to centuries. For cases 5 and 6 the situation is more complicated as both denitrification and N fixation were perturbed. Calculating the Case 5 coupled feedback as for Case 1 above gives a value of 0.48, indicating that the oceans would have lost an additional 48% more N if N fixation had remained constant. This again highlights the stronger negative feedback response in the more iron-replete ocean (comparing Cases 1 and 5). Additional References Anderson, L. A., and J. L. Sarmiento, 1994. Redfield ratios of remineralization determined by nutrient data analysis. Global Biogeochem. Cycles, 8: 65-80. Armstrong, R. A., C. Lee, J. I. Hedges, S. Honjo, and S. G. Wakeham, 2002. A new, mechanistic model for organic carbon fluxes in the ocean based on the quantitative association of POC with ballast minerals, Deep-Sea Res., Part II, 49, 219–236. Coles, V.J., Hood, R.R., Pascual, M., Capone, D.G., 2004. Modeling the impact of Trichodesmium and nitrogen fixation in the Atlantic Ocean. J. Geophys. Res., 109, C06007, doi:10.1029/2002JC001754. Krauk, J.M., Villareal, T.A., Sohm, J.A., Montoya, J.P., Capone, D.G., 2006. Plasticity of N:P ratios in laboratory and field populations of Trichodesmium spp., Aquatic Microbial Ecology, 42, 243-253. Fennel, K., Spitz, Y.H., Letelier, R.M., Abbott, M.R., 2002. A deterministic model for N2-fixation at the HOT site in the subtropical North Pacific. Deep-Sea Research II, 49, 149-174. Hood, R.R., Bates, N.R., Capone, D.G., Olson, D.B., 2001. Modeling the effect of nitrogen fixation on carbon and nitrogen fluxes at BATS. Deep-Sea Research II, 48, 1609-1648. Hood, R.R., Laws, E.A., Armstrong, R.A., Bates, N.R., Brown, C.W., and others, 2006. Pelagic functional group modeling: Progress, challenges, and prospects. Deep-Sea Research II, 53, 459-512. Le Quéré, C., S.P. Harrison, I.C. Prentice, E.T. Buitenhuis, O. Aumont, and others, 2005. Ecosystem dynamics based on plankton functional types for global ocean biogeochemistry models, Global Change Biology, 11, 2016-2040. Van Mooy, B.A.S., Keil, R.G., Devol, A.H., 2002. Impact of suboxia on sinking particulate organic carbon: Enhanced carbon flux and preferential degradation of amino acids via denitrificiation. Geochimica et Cosmochimica Acta, Vol. 66, 457-465. Auxiliary Figures Figure S1. Selected global scale fluxes are shown for each case experiment. The first 200 years of Case 1 are plotted for easier comparison. The last 20 years of the control are also shown in each plot. Figure S2. The N* tracer (= ([NO3-] + [NH4+]) - 16 * [PO4-]) averaged over sub-euphotic zone depths (103-215m) is plotted over time for each of the case experiments at global and basin scales. Figure S3. The minimum O2 concentration within the water column at each location is plotted from the World Ocean Atlas 2001 data, from our control (year 2050) and case experiment simulations (year 200). Figure S4. Spatial patterns are displayed for the N* tracer (= ([NO3-] + [NH4+]) - 16 * [PO4-]) averaged over euphotic zone depths (< 103m) for the WOA2001 and each of our case experiments. 1