Assimilating water column and satellite data for marine export production estimation

Recent advances in satellite retrieval methodology now allow for estimation of particular organic carbon (POC) concentration in ocean surface waters directly from satellite-based optical data. Because of the good coverage, these data reveal small-scale spatial and temporal concentration gradients and document the evolution of surface water POC as well as the underlying driving biogeochemical processes throughout the seasons. Water column nutrient data also reveal biogeochemical activity. However, because of the scarcity of data, the deduction of temporal changes of particle production and export is not possible in most parts of the ocean. Here we present first results from a new study combining both data streams, thereby exploiting the high spatio-temporal resolution of surface POC concentrations from satellite optical sensors with water column nutrient data having sparser coverage but providing information throughout the entire water column. We use a mediumresolution global model with steady-state 3-D circulation that has been optimized by fitting to a large number of hydrographic parameters and tracers, including CFCs and natural radiocarbon. Production and export of POC is allowed to vary monthly, and the magnitudes of the monthly export fluxes are determined by fitting the model to satellite POC data as well as water column nutrient data using the adjoint method. Two cases have been investigated: (1) the production rate of POC is set to be proportional to export production (EP) and the seasonal changes are assumed sinusoidal (meridionally varying amplitude and phase), and (2) the POC production rate is linked to primary production rates (literature). Both cases were run with the same initial state and model settings, and show total cost function decreases of 12 and 95 %, respectively. The POC misfit term alone decreased by 75 and 99.8 %. The integrated annual global POC exports of the two cases are 9.9 and 12.3 Gt C yr −1, respectively. Overall, the remaining POC and phosphate misfits of both solutions are considered too large, and the difference fields still xhibit significant systematic geographical patterns. This indicates that the present model runs are too simplistic and do not fully explain the data. Further, more refined model setups are needed.


Introduction
The ocean is one of the major carbon reservoirs on Earth, containing around 40 000 Gt C, about 50 times more than in the atmosphere (Sarmiento and Sundquist, 1992).The ocean is believed to be the ultimate sink for about 90 % of human fossil fuel emissions (Archer et al., 1998).In the surface ocean, phytoplankton fixes dissolved carbon to form particulate biomass by photosynthesis (primary production PP).A fraction of the particulate material sinks into the deep ocean due to gravity, thereby sustaining a downward flux of particulate nutrients and carbon (export flux).This process, which depletes the ocean surface of nutrients and dissolved inorganic carbon (DIC) relative to the deep water, is referred to as the biological pump (Volk and Hoffert, 1985).Drawdown of surface carbon concentrations by the biological pump leads to an increased flux of CO 2 from the atmosphere, and the overall oceanic CO 2 uptake thus depends on the strength of the biological pump.Quantification of export flux and the strength of the biological pump therefore is an important objective.
Over the past decades, extensive work has been carried out to quantify the downward carbon export in the ocean.One observational approach is the use of sediment traps (Honjo et Published by Copernicus Publications on behalf of the European Geosciences Union. X. Yao and R. Schlitzer: Assimilating water column and satellite data al., 2008;Gardner, 2000;Kahler and Bauerfeind, 2001).This is a direct way of measuring downward export fluxes by capturing and preserving the sinking material.However, the disadvantages are sparse coverage in space and time.Moreover, the catchment efficiency of sediment traps, especially in the shallow waters, are debated (Gust et al., 1994), casting doubt on absolute flux values derived from shallow traps.Additionally, there is another alternative instrument of measuring carbon export directly, it is carried on ARGO floats, called "Carbo-ARGO" floats (Bishop and Wood, 2009).Besides these direct measurement approaches, radioisotope 234 Th is also widely utilized to quantify particle export out of the surface layer (Buesseler et al., 2009;Rutgers v. d. Loeff et al., 2011).
In addition to observational approaches, modeling is a powerful way of quantifying carbon export.Different model approaches are used to estimate carbon export, such as ecosystem models (Aumont et al., 2003;Oschlies and Kähler, 2004) or coupled physical-biogeochemical models (Palmer and Totterdell, 2001).Performance of these forward models in terms of how well measured distributions are simulated highly depends on clever choices of model parameters, such as primary production and functions that relate export flux with PP.The inverse model of Schlitzer (2000Schlitzer ( , 2002) ) avoids these parameter choices and treats the geographically varying annual export fluxes as independent parameters determined by the model by exploiting historical water column data.However, this approach only yields the annual average fluxes and does not reveal seasonal changes.
Recent advances in remote sensing now allow for direct quantification of surface carbon parameters such as the concentrations of particulate organic carbon (POC) (Stramski et al., 2008) and particulate inorganic carbon (PIC) (Balch et al., 2005;Gordon et al., 2001).The satellite-derived fields are provided globally at high spatial and temporal resolution, and coverage is excellent, especially when compared to shipbased water column data.Both of these quantities are sensitive to biological production and export; therefore analysis of this new data product likely helps in revealing flux variations on small space and timescales.The disadvantage is that satellite sensors only monitor the ocean surface.A combination of satellite and water column data can build up a full 4-D view of the global ocean.
The present study uses the data-driven inverse model developed by Schlitzer (1993Schlitzer ( , 2000Schlitzer ( , 2002Schlitzer ( , 2007) ) to link water column and satellite data for estimating carbon export and its seasonality.Based on the older model version, we set up a new model combining water column phosphate and satellite POC data in the adjoint model in order to get better carbon export estimations on a monthly basis.The model is calculated in phosphor units, and using a constant Redfield ratio (Redfield et al., 1963;Anderson and Sarmiento, 1994) relating between phosphor and carbon.
This paper is structured as follows: Sect. 2 introduces the data we used in our study, Sect. 3 describes the details of the method of the adjoint model approach, Sect. 4 describes the two experiments we performed, Sect. 5 shows the analysis of the results of two experiments, and Sect.6 contains the summary and conclusion.

Data
There are two kinds of data used in this study: ship-based measurements of dissolved inorganic nutrients (C Dd ) in the water column, and satellite-derived particulate (C Pd ) data for the ocean surface water.The dissolved nutrient data cover the whole global ocean and extend from the surface to the ocean bottom.The satellite-derived particulate data provide excellent spatial and temporal coverage (much better the water column data), but data are only available for the ocean surface.Combining these two diverse data types is the novel feature of the present study.

Water column data
The water column nutrient data are taken from the World Ocean Atlas 2009 (WOA09) (Garcia et al., 2010).WOA09 is based on a compilation of all presently publicly available cruise data.These original data have been objectively analyzed, and climatological fields on a 1 • × 1 • horizontal grid and at 32 standard depth levels have been produced for annual, seasonal, and monthly coverage for the world ocean.The WOA09 nutrient and oxygen fields reveal important features that reflect the action of physical as well as biogeochemical processes.The main aim of this study is to utilize the data and infer and quantify the strengths of the underlying biogeochemical processes.The emphasis here will be on the export flux of particulate carbon from the surface ocean into the deep.
As examples of nutrient distributions Fig. 1 shows surface ocean phosphate fields for June and December.Clearly visible are the very low phosphate concentrations in subtropical regions, whereas quite high concentrations are found in the northern North Atlantic, the North and equatorial Pacific, and especially in a circumpolar belt around Antarctica.We can also see that phosphate concentrations are generally lower during the productive season (June in the Northern Hemisphere, December in the Southern Hemisphere), reflecting the nutrient drawdown due to the production of particulate material.This is also the case for the Antarctic high-nutrient belt, where austral summer (December) concentrations are much smaller than during austral winter (June).
Due to the combination of about 80 yr of measurements included in WOA09 this climatological dataset has quite good global coverage horizontally and with depth.Data coverage for any specific month (e.g., January 2000), however, is very sparse, and fields as shown in Fig. 1 cannot be constructed because of lack of data.This situation is unlikely to improve in the future because of the high logistical and financial costs of ship observations and the local nature of the obtained information.Many of the station data used for WOA09 have only been sampled in the upper water column, and, as a consequence, coverage within the upper 500 m is much better than below.Because of the especially poor data coverage in the deep, WOA09 monthly values are only provided for the upper 500 m, whereas only annual average data are provided below.Due to the relative stability of the deeper ocean, the temporal variation of phosphate concentrations is weak, such that usage of annual phosphate values is justified.We take phosphate standard errors provided by WOA09 or 0.03 µmol kg −1 (whichever is greater) as the data error σ C Dd in the data assimilation procedure described below.The 0.03 µmol kg −1 value corresponds to the typical expected error of historical phosphate data.
In this study, we use phosphate (and phosphorous) as a proxy for carbon, and utilize the fact that abundance ratios of carbon and phosphorous in marine particulate organic material are nearly constant (Redfield et al., 1963).The C / P Redfield ratio is used to infer carbon export from phosphorous export.The reason for this approach is that there are about a thousand times more phosphate data compared to carbon.

Satellite POC data
Development of new retrieval methods exploiting satellitederived optical data of the ocean surface is an active field of research, and recently, new products have been released  (Stramski et al., 2008).
that for the first time provide estimates of particulate organic carbon (POC) concentrations for the global ocean with unprecedented spatial and temporal resolution.In this study we use the product of Stramski et al. (2008) downloaded from http://oceandata.sci.gsfc.nasa.gov/cgi/l3.This POC product uses data from the SeaWiFS satellite sensor on the SeaStar satellite.SeaWiFS covers the period from 18 September 1997 until 11 December 2010 and provides 10 yr of continuous observations.The release of highresolution satellite-derived surface ocean POC data is a major step forward, and now allows for direct inference of carbon export compared to the more indirect older approaches that were based on chlorophyll.
Satellite POC data have much better temporal and spatial resolution than water column data.As mentioned before, the existing water column data allow for construction of climatological fields, but are still too sparse to derive instantaneous distributions on the basin or global scale.In contrast, satellites monitor the entire ocean and revisit individual locations typically within one or two weeks.Therefore, satellite data are revealing small-scale features (9 km resolution) and documenting their temporal evolution in great detail.In this study we use satellite data to fill up gaps in the water column data.
Figure 2 shows surface ocean POC distributions for June and December (Stramski et al., 2008).POC concentrations are generally low in the subtropical gyres, while they are quite high in the biologically productive regions (equatorial and coastal upwelling regions, subpolar and polar regions during productive summer periods).The lack of coverage in polar regions is due to light limitation during wintertime.
While providing high spatial and temporal resolution and excellent coverage, satellite data have the disadvantage of not being able to "see" below the ocean surface, and thus no subsurface POC data exist.Considering advantages and disadvantages of satellite POC and water column data, it becomes clear that they are mostly complementary, and that the combination of the two will provide information not available when using any of these sets individually.In this study we assimilate both data streams and estimate the climatological monthly export of carbon in the world.
As the amount of phosphate water column data is much larger than for carbon, the present model is formulated in phosphorous units.Therefore, before using the satellite POC data in the model, it is necessary to convert carbon to phosphorous, using a fixed uniform Redfield ratio (C / P = 106, Redfield, 1963).Instead of POC the model uses the soconverted particulate organic phosphorous (POP) denoted as C Pd .Note that the conversion does not change the structure of the features as, for instance, shown in Fig. 2. The error of the satellite POC values is estimated at ± 30 % (M.Stramska, personal communication, 2010).As data error in the model we use 30 % or 13 mg m −3 (Stramski et al., 2008), whichever is larger.The associated POP errors σ C Pd are obtained by applying the Redfield ratio above.

Model
The model used in our study is an extension of the adjoint model of Schlitzer (2007).The general model strategy is described in detail in Schlitzer (1993Schlitzer ( , 1995Schlitzer ( , 2000Schlitzer ( , 2002)).In the original model, all the tracers are in steady state except chlorofluorocarbons (CFC).Following the treatment of CFC, we set up the model with phosphate and POP as time-dependant model tracers with monthly resolution.Phosphate and POP budgets contain the effects of circulation as well as particle production, degradation, and export.

Model grid
The model is global and has a nonuniform grid with horizontal resolution ranging between 1 • × 1 • and 4 • × 5 • (Fig. 3).Finer resolution is realized near coastal regions, while coarser resolution prevails in the open ocean.The model has 2421 columns and three boundary columns adjacent to the Mediterranean Sea, the Red Sea, and the Persian Gulf.There are 26 vertical layers, with thickness progressively increasing from 60 m at the surface to 500 m at 5000 m depth.
Figure 4 shows the main processes in the model.The model has a steady-state 3-D circulation field that was obtained by fitting the model to a large set of tracer data including CFC and natural radiocarbon (C 14 ) (Schlitzer, 2007).The

Model budgets equations
The model simulates two tracers, dissolved inorganic phosphate C D and POP C P .C D and C P are coupled in the surface layer (see Fig. 4).During the productive period, phosphate C D is consumed by phytoplankton, leading to the build-up of POP C P , while the remineralization of C P leads to a release of phosphate and thus to an increase of C D .Phosphate C D budgets are formulated for all boxes of the model throughout the entire water column, while budgets of suspended POP C P only exist in the surface layer.The sinking part of C P is remineralized while settling through the water column.
The concentrations of C D and C P vary with time.The C D concentrations are transported by circulation and are affected by biological production (sink) and particle remineralization (source).When POP C P is formed during phytoplankton photosynthesis near the ocean surface, dissolved phosphate C D is consumed.The C D budget for ocean surface boxes is as follows: TheC P budget for ocean surface boxes has the form The phosphate C D budgets in the deeper ocean have the same form as Eq. ( 1) except that the last term V dC P dt vanishes because of negligible subsurface C P concentrations.
In Eq. ( 1), A is the area of an interface and has either a positive or negative sign, depending on whether a positive flow u, v, or w (eastward, northward, or upward) enters or leaves the box.V is the box volume, and ∇C Di , ∇C Dj , ∇C Dk are C D surface normal gradients on the respective interfaces.The subscripts i, j , and k indicate zonal, meridional, and vertical directions.In the present model, the horizontal flows are steady state and taken from Schlitzer (2007).The vertical velocity w is calculated by solving the continuity equation from bottom to top in each column.K h and K v are horizontal and vertical mixing coefficients.
The vertical mixing coefficient K v is affected by changes of the mixed layer depth (MLD) taken from Monterey and Levitus (1997).A deepening of MLD in the winter increases K v within the MLD, leading to an intense vertical mixing of tracers within the MLD.In terms of phosphate this brings up nutrients to the upper ocean from deeper layers, thus fueling the next spring/summer bloom to come.In the model the enhanced vertical mixing within the MLD is implemented by applying a large factor on K v for all vertical interfaces within the MLD.The upscaling factor is reduced in a smooth way to 1 for the first interface below the MLD.
q D is the sink/source due to loss/gain by sinking particles.As in Schlitzer (2007) the vertical particle flux j p (z) is assumed to follow so-called Martin curves (Martin et al., 1987).Production of particulate material occurs in the top two model layers (euphotic zone).The euphotic zone z EZ is 133 m.Sinking particles are remineralized below the euphotic zone, and the particle flux decreases according to (3) In Eq. ( 3), the exponent b determines the shape of the particle flux profile, and thus controls the depth of remineralization.Here we use b values provided by Schlitzer (2007).
The parameter a represents the particle flux at the depth of the euphotic zone, commonly referred to as export production.
The C P budget in Eq. ( 2) simply links temporal concentration changes with the net source/sink term q p .We consider two different source/sink setups.(A) the source part of q p is set proportional to the particle export flux at the base of the euphotic zone, and the sink term describes exponential decay of the particulate material with a prescribed lifetime τ .(B) the source part of q p is set proportional to the net primary production (NPP) obtained from satellite data, and the sink term again describes exponential decay of particulate material, but the lifetime is temperature dependent and model adjustable.These two setups are described in detail below.
In the first experiment (Exp A) the C P budget has the following form: The source term of C P is proportional to export production a, scaled with a factor α. The sink term represents exponential degradation of C P with a prescribed lifetime τ .The export production a follows Eq. ( 5) below.
where α 0 • p 2 e (x, y) is the export flux of the steady-state solution of Schlitzer (2007) at longitude x and latitude y, and p 2 e (x, y) is the square of the spatially varying export flux parameter, which can be adjusted in the adjoint model.Formulating the flux as square of the adjustable parameter guarantees positive export flux.s(y, t) = 1 + κ(y) • sin( (y) + 2π t) is a prescribed season factor, which controls the seasonal variations and time of bloom; t is time; and κ(y) is the relative amplitude of the seasonal export variations ranging between 0 (no seasonal variations) and 1 (maximal seasonal variations).κ(y) and (y) are simplistic, prescribed formulations ensuring weak seasonality at low latitudes and strong seasonality at high latitudes.Blooms occur in spring at low latitudes and in late summer at high latitudes.
In the second experiment (Exp B) the C P budgets has the following form: Here the source of particulate material is set proportional to satellite-derived new primary production N P (Behrenfeld and Falkowski, 1997) scaled with the square of a modeladjustable parameter p α (using the square guarantees a positive source term).
In the original satellite-derived N P dataset there are no data in the polar regions, typically for some months during wintertime; however, a complete N P field is required by the model.As a fill-in strategy we take 10 % of the maximum N P value over the given box surface area.This seems to be a reasonable strategy as wintertime productivity values are expected to be much smaller than during bloom periods.In Eq. ( 6), the sink term is the C P remineralization influenced by water temperature.The remineralization rate is proportional to the square of another model-adjustable parameter p γ times a temperature-dependent factor Q. The temperature factor Q follows Eppley (1972), and the monthly temperature data from WOA09 temperature (Locarnini et al., 2010) are used.A in Eq. ( 6) is the box surface area.

Model parameters
There are two groups of parameters in the model: independent parameters p * , which can be adjusted by the model in the course of the optimization, and dependent parameters p, which depend on the independent parameters p * , and can be calculated by performing a model simulation.
The sets of independent parameters are different for the two experiments described above: There are three types of independent parameters: export parameters p e , remineralization parameters p γ , and productivity parameters p α .Note that all independent parameters appear in squared form in the model equations to ensure proper sign of the terms.For each type, there is one independent parameter for every model column; for example, each type of independent parameter describes an entire 2-D global field.Exp B has three times as many independent parameters as Exp A, and should therefore have more flexibility in fitting the data.
Exp A only contains the export-related independent parameters p e .These parameters not only appear in the source/sink term q D of the C D budgets but also in the source term of the C P budgets.Therefore, in Exp A, varying the independent parameters p e impacts both the C D and C P budgets and tightly couples the two.Strong export (large p e ) leads to large decrease of dissolved phosphate (drawdown) and at the same time to large build-up of particulate matter.Exp A therefore has an inherent anticorrelation between the two tracers.
Runs of Exp A are started with initial p e values from Schlitzer (2007).During the run p e values are adjusted to satisfy the water column phosphate and ocean surface POP data as closely as possible.
Exp B also includes the export-related independent parameters p e but has two additional sets of independent parameters, p γ and p α , affecting the POP source and sink terms in Eq. ( 6).p 2 γ is initialized with the globally uniform value of 0.048 day −1 (Schartau and Oschlies, 2003).p 2 α is initialized with a globally uniform value of 0.75.
The set of dependent parameters is composed of the model-simulated concentrations of the two tracers phosphate and POP: The size of dependent parameter vector is determined by the number of C D and C P boxes times the number of model From the observed misfits the adjoint model determines modifications of the independent parameters that will lead to a better next simulation.

Cost function
Once C D and C P are calculated during the model forward run, the C D and C P values of the last 12 months are stored and compared with observations.The model-data misfits as well as other undesired features are accumulated in the cost function F of the model.The individual terms of cost function defined in this study are listed in Table 2.The total cost function F is the sum of all the terms listed (see Eq. 9): There are two terms accumulating squared model-data phosphate and POC misfits and a third term penalizing spatial roughness of the export-related independent parameters p e (via second derivatives in east-west and north-south directions).The phosphate misfits in F phosphate are calculated for three depth ranges separately and then summed up.The reason for the subdivision is the increased flexibility by allowing for individual weights for the individual terms.The weights used for the different experiments are identical (Table 2).
Following standard procedures, model-data misfits for phosphate and POC are calculated as error-normalized differences between model simulated and observed values and then squared.The data errors used are specified in Table 2.It is assumed that data errors are uncorrelated.Given that the climatological phosphate data are the result of very many independent observations over the last 80 yr, this assumption seems justified.The algorithms for estimating satellitederived POC values may in principle introduce correlated errors; however, the magnitude of such effects (if they exist) are yet unknown.

Adjoint equations
The objective of the adjoint model is to find the minimum of cost function F under the additional condition that all model equations are satisfied exactly.For Exp A the model equations are Eqs.( 1) and (4), while for Exp B we have Eqs.( 1) and ( 6).The first step in the derivation of the adjoint equations is to rewrite the model equations in homogeneous form E i = 0, i = 1, . . ., n e .Based on the cost function F , we then formulate the Lagrangian L of the model, Here n e is the number of model equations, and λ j are the Lagrangian multipliers of the problem.
Seeking the minimum of cost function F is equivalent to finding a model solution that satisfies Eqs. ( 11), ( 12) and (13).∂L ∂λ j = E j = 0 (11) Equation ( 11) represents the model equations, and is automatically fulfilled.One proceeds by calculating the Lagrangian multipliers λ j from the Eq. ( 12).Then the calculated Lagrangian multipliers are inserted into Eq.( 13) to calculate the gradient of the cost function F with respect to the independent parameters p * .This gradient is then passed to a descent algorithm that returns with a new and improved set of independent parameters.Repeating this process iteratively will ultimately lead to a vanishing gradient vector, and thus also satisfying Eq. ( 13).The derivation and implementation of the adjoint equations is described in detail in Schlitzer (2007).
The general sequence of computational steps is shown in Fig. 5.The advantage of the adjoint part is the adjustment of independent model parameters that is guaranteed to lead to smaller cost function values in subsequent simulations.This is vastly superior to manual parameter optimizations that in most cases do not lead to improved simulations, or, in cases with thousands of parameters, is simply not technically possible.As explained in Sect. 4 below, the adjoint model consists in a backward time-stepping procedure with the same number of steps as the forward simulation.The computational cost of an adjoint run is similar to the cost of the forward run.A full iteration as shown in Fig. 5 thus requires about twice the effort of a forward simulation.Hundreds or thousands of such iterations need to be performed to reach the minimum in large-scale problems.Our present model can be run on a typical PC, and computing time for finding optimal solutions is about a day.
The adjoint code of the present study was obtained manually, and various efficiency measures were implemented.

Model experiments
As shown in Fig. 5 the model calculations can be divided into two parts.The first part is a normal forward run (simulation) that, for a given set of independent parameters, calculates the dependent parameters (tracers) C D and C P by forward time-stepping of the tracer budget equations.This step uses traditional procedures as included in most existing physical or biogeochemical models.The second part involves running the adjoint of the forward model including comparisons of the simulated tracer values with data and the accumulation of misfits in the cost function.Ultimately the gradient of the cost function with respect to the independent parameters is obtained by time-stepping the adjoint model backwards.Once this gradient is known, a new, improved set of independent parameter values is produced using a suitable descent algorithm (here: quasi-Newton conjugate gradient algorithm).In essence the parameter improvements are determined (driven) by the data and data misfits.The adjoint model guarantees that the next simulation with improved model parameters will have more realistic tracer fields C D and C P and a lower value of the cost function.Many such iterations consisting of simulation and an adjoint model run have to be repeated until a termination criteria is satisfied and an optimal solution is reached.
In this study, we have run the two experiments Exp A and Exp B described in Sects.3.2 and 3.3 above.In each case, we solve the budgets equations for phosphate C D and POP C P using an implicit time-stepping scheme that allows for large time steps of one month.We integrate over a 10 yr period, sufficient for upper ocean phosphate and POP fields to reach equilibrium.Initial values for phosphate C D are from Schlitzer (2007), and POP concentrations C P are initialized with zero.
The simulation produces 120 snapshot fields of C D and C P concentrations, which are kept in memory for usage during the backward adjoint run.Only the simulated phosphate and POP concentrations during the last year are compared with data on a month-by-month basis.Model-data misfits are accumulated in the cost function as discussed above and listed in Table 2.Note that, in addition to the data misfit terms, a smoothness term for export-related independent parameters is also included in the cost function.This term is never dominant in Exp A or Exp B (see below).
The adjoint equations Eqs. ( 12) and ( 13) have special structure allowing for very efficient calculation of Lagrangian multipliers using a backward-in-time stepping of the adjoint equations Eq. ( 12) and finally of the gradient of the cost function with respect to the independent parameters.The computational cost for Exp B is higher than for Exp A due to the three times larger number of independent parameters.
Before starting the model runs, the correctness of the cost function gradient was checked by applying the gradient check described in the Appendix for a number of independent parameters.Failures of the test in all cases were caused by coding errors.The production runs were only run after all errors were fixed and all gradient checks succeeded.

Cost function values
Initial and final (optimal) values of the total cost function as well as for the individual terms for the two experiments Exp A and Exp B are shown in Table 3.Note that both experiments were run with identical weight factors for the individual terms, thus making the values directly comparable.The initial cost function value for Exp B (2.9 × 10 9 ) is much larger than for Exp A mostly because of the very large POC misfit term caused by poor initial values of the p α -and p γindependent parameters that determine the POC source and sink strengths.In all cases at initial stage the POC and upper ocean phosphate misfit terms are the dominant ones, while the deep phosphate misfit and spatial smoothness terms are about one order of magnitude smaller.
In all runs we find decreases of the total cost function values as a result of the iterative optimization procedure.In Exp A the reduction amounts to about 12 %, mostly due to the decrease of the POC misfit term (reduced from 4 to 1×10 7 ).The phosphate misfit terms for the surface and deep ocean remain almost unchanged, while the term for twilight zone phosphate decreases by about 12 %.The overall cost function reduction for Exp B is much larger than for Exp A and amounts to about 95 %.This large reduction is mostly due to reductions of the POC misfit term (from about 267 to 0.3×10 7 ) and the surface phosphate misfit term (from 22.4 to 10.8×10 7 ).The other terms remain almost unchanged.Compared to Exp A, Exp B has reached a much better state with the total cost function and the surface misfit term, amounting to only half the Exp A values and the POC misfit only being one-third.Obviously the simulated phosphate and POC fields of Exp B are closer to the observations than the ones for Exp A (see below).
In Exp A, the export-production-related parameters p e are the only adjustable parameters, and by design of Exp A, any adjustment of p e directly influences the simulated surface POC and phosphate concentrations in a prescribed anticorrelated way.In regions where observations show simultaneous increase or decrease of POC and phosphate concentrations, Exp A has no way of fitting both data types because a better fit of one type inevitably is associated with a poorer fit of the other.Such a situation is shown in Fig. 6a for a location in the central North Atlantic.In this region both data types, surface phosphate as well as POC, decrease from April until September.Exp A has no way to reproduce such behavior, and produces quite unrealistic phosphate and POC values during the entire period.In Exp B the linkage between phosphate and POC is less strict because of the additional independent parameter sets p α and p γ that affect the surface POC distribution but not phosphate.Thus Exp B has more flexibility than Exp A. The inherent limitation of Exp A is also reflected in the final misfit values for POC and surface phosphate, which both are about a factor of 3 larger than for Exp B.

Misfit analysis
The mean and root-mean-square (RMS) misfit values of POC and phosphate are listed in Table 4. Consistent with the smaller cost function values discussed above, Exp B also exhibits smaller POC mean and RMS misfits than Exp A. Both Exp A and Exp B show an underestimation of POC concentrations, which in the case of Exp B is quite small.
The phosphate misfits are calculated separately for the three depth ranges also used in the cost function terms: surface layer (< 60 m), twilight zone (60 to 400 m), and deep water (> 400 m).In general, Exp B exhibits smaller mean and RMS phosphate misfits in the surface layer than Exp A. This is consistent with the smaller value of the phosphate (< 60 m) cost function term of Exp B described above.The mean phosphate misfit of Exp B in the twilight zone is slightly larger than for Exp A. In the deep ocean both experiments produce about the same misfits.Overall, Exp B performs better than Exp A, especially in the surface layer.

Surface phosphate and POC analysis
Figure 7 shows the model-simulated phosphate fields for June and December.Overall, both experiments reproduce the main features in the observations (see Fig. 1) well.We find low concentrations in the subtropical gyres and high concentrations in the subpolar regions as well as in equatorial and coastal upwelling regions.In Exp A, there are occasional small negative values in parts of the Atlantic subtropical gyres.These are unrealistic values, caused in the model by insufficient phosphate resupply to the ocean surface.In Exp B, these unrealistic values are not observed.
When comparing more closely, the model-data phosphate misfits for Exp A (Fig. 7b and d) reveal significant and large-scale systematic offsets.Exp A overestimates surface phosphate in the entire north Pacific, in the southeast Pacific upwelling region and at high latitudes in the South and North Atlantic.Phosphate underestimation in Exp A occurs at low latitudes, especially in the Atlantic.This general misfit pattern is found for June and December, and seems to prevail throughout the year.In Exp B the magnitude of surface phosphate misfits is generally smaller than for Exp A, but large-scale systematic features are also observed.
Figure 8 shows the model-simulated POC fields for June and December.Both Exp A and Exp B capture the main structure of satellite POC data: low concentrations in the subtropical gyres, and high concentrations in the coastal regions.However, significant systematic misfits are present both in Exp A and Exp B. In Exp A, the POC concentrations are underestimated in most regions in June and December.They are overestimated in the Southern Ocean during In Exp B, the extent of overestimation and underestimation is smaller compared to Exp A, but large-scale systematic patterns are also found.
The agreement between simulated and observed POC values for both experiments is best at low and poorest at high latitudes (Fig. 8b, d, f, and h).The large POC misfits in subpolar and polar regions may again be caused by model shortcomings.However, it is also possible that the POC data quality is poorer in these regions because of scarcity of satellite optical data due to light limitation and frequent cloud coverage.In addition, because of the high logistical cost, direct, shipboard measurements of in situ POC concentrations required for the calibration of the POC retrieval algorithms are rare (Stramski et al., 2008), and it is quite possible that the satellite-based POC data still contain systematic biases at these high latitudes.
As an example of the temporal evolution of surface phosphate and POC concentrations over the course of an entire year, Fig. 6 shows data values and model simulations for the location in the central North Atlantic marked in Fig. 3.This location was chosen because of relatively large seasonal changes of both the dissolved nutrient as well as POC concentrations (see Figs. 1 and 2).At this location surface phosphate concentrations are highest during winter and spring (Fig. 6a), when the MLD is deepest (Fig. 6d).During late spring, when the upper water column warms and MLD decreases, phosphate concentrations start to decline as biological production is starting to use up nutrients.The phosphate drawdown continues until September, when finally concentrations begin to increase again towards winter values.During winter and spring, when phosphate is highest, POC concentrations are lowest.At the chosen location, POC concentrations start rising in early spring, reaching maximum values in April and exhibiting a steady decline during most of the summer season.This decline occurs while biological production is high, as indicated by the decreasing phosphate concentrations.
Figure 6b and c show that Exp B is able to reproduce the observed trends in phosphate and POC, while results of Exp A are far from the data and essentially show very small seasonal amplitude in both cases.As the main reason for this failure of Exp A we see the rigidity of its POC source/sink setup (Eqs.4 and 5), which strictly links phosphate drawdown with large POC source term and increasing POC concentrations.Simultaneous decreases of both phosphate and POC such as in Fig. 6 cannot be reproduced by Exp A by design.

POC export analysis
Figure 9 shows the downward POC export flux at the base of the euphotic zone (133 m depth) of Exp A and Exp B for June and December.Both Exp A and Exp B show large seasonal changes of the POC export, especially in the high-export regions.In June, high POC exports are found at high latitudes in the North Atlantic, the North Pacific, and near the east coast of Africa.In December, the highest exports occur in the Southern Ocean, the coastal region of Africa, and the southern subpolar regions.The seasonal evolution of high-export regions matches the seasonal changes of high productivity as seen, for instance, in satellite chlorophyll maps.The annually averaged global POC export of Exp A is 9.9 Gt C yr −1 , while for Exp B we obtain an about 25 % higher flux of 12.3 Gt C yr −1 .These values are within the wide range of literature values between 11 and 22 Gt C yr −1 (Laws et al., 2000;Schlitzer, 2000;Eppley and Peterson, 1979).

Identifiability of independent parameters
The notion of identifiability addresses the question of whether it is at all possible to obtain unique solutions of the inverse problem for unknown parameters of interest in a meteorological/oceanic model from data collected in the spatial and temporal domains (Navon, 1998).In order to address the identifiability of our adjoint model, we conducted Experiment Exp B1 has to be considered an extreme case with all independent parameters close to zero.Any resulting structures in the final solutions are developed during the data assimilation process and are not already included in the initial fields.
Results from these runs are listed in Table 5.While starting with wildly different initial cost function values (not shown), all sensitivity runs show similar final values for the total cost function as well as for individual terms, all differing by less than 5 %.The export-related independent parameters p e in the three sensitivity runs are surprisingly similar given the vastly different initial fields, especially for B1 with its spatially constant near-zero values.The mean values of p 2 e and the associated global carbon exports agree within 2 % (Table 5) and the spatial fields all reveal nearly identical export patterns (not shown).The situation is different for the two other groups of independent parameters p α and p γ that determine sources and sinks of POC.These two parameter groups appear to be positively correlated, and similar POC distributions can obviously be produced (see POC misfit term in Table 5) with both POC sources (p α ) and sinks (p γ ) small or large.The observed spread in the three runs is almost 100 %.Overall, we find a situation that is quite common in optimization problems, where some parameters are tightly constrained (export parameters) while others vary largely.The fact that Exp B1, which starts from a "state of ignorance", still produces a carbon export field nearly identical to the other two runs shows that indeed the water column nutrient and satellite-based surface POC data contain sufficient information for estimating the marine carbon export flux.

Summary and conclusions
In this study, we use a medium-resolution coupled biogeochemical-physical ocean model and assimilate satellite POC as well as water column nutrient data to estimate the seasonal evolution of carbon export in the global ocean.As an extension of previous work, the addition of satellite data with their good spatial and temporal coverage allows for estimating the temporal variation of the global carbon export fields.The extended model simulates surface water POC concentrations in addition to nutrient concentrations throughout the water column.Both simulated fields are compared with  respective observations, and model misfits are accumulated in the cost function of the model.The adjoint method is applied to drive the model to the satellite POC as well as water column nutrient data and to optimize the export production values.Experiments are done using two different POP production scenarios (C / P Redfield ratio is used to convert between POC and POP), one with a POP source term being proportional to export production (Exp A), while in the second case the POP source is proportional to independently obtained new primary production (Exp B).
The model results show that the adjoint method worked well and that significant decreases of the cost function were achieved.Final (optimal) phosphate and POP model distributions agree much better with observations than initially, and the model has driven the model simulations closer to observations.The integrated carbon export in the two experiments amounts to 9.9 (Exp A) and 12.3 Gt C yr −1 (Exp B).
In the surface layer, we allow for material exchanges between dissolved phosphate and particulate phosphorous POP in a way that makes the two essentially anticorrelated.When POP builds up, phosphate is consumed, leading to a decrease in phosphate concentration.On the other hand, when POP is remineralized, dissolved phosphate is released, and its concentrations rise.Such an anticorrelation between the two parameters is actually also found in the data of water column phosphate and satellite POC in many regions.However, in some locations and months, the observations show a positive correlation between the two, with both parameters simultaneously increasing or decreasing for some period of time.This behavior cannot be reproduced by the model because of a setup implying anticorrelation.In such regions, simulated C D and C P values often differ much from observations, and the model has no way of improving both tracers at the same time because of the inherent anticorrelation built into the model.This is especially the case in Exp A with its strict linkage of POC export and POC production.In Exp B, the influence of this design limitation is reduced by the two additional sets of independent parameters.
The existence of significant and systematic differences between model and observations in the final, optimal solutions strongly suggests that the treatment of POP budgets and the coupling with dissolved nutrients is overly simplistic and unrealistic in the present setup.This is especially true for Exp A, which exhibits relatively small improvements in the POP fields.While the present study has shown that in principle the adjoint method can be applied for determination of time-varying export flux fields using satellite and water column data, the present results have to be considered preliminary and more refined model setups coupling the dissolved, and particulate phases of nutrients and carbon are needed.

Gradient test
Before the actual model runs, a gradient test (Navon et al., 1992) was performed to check the correctness of the adjoint model.The mathematics is based on Taylor expansion of the cost function.The cost function at a perturbed point p + αdp can be written as The quantity R defined as should vary quadratically with α for small enough αdp.The test involves repeated calculation of R for a sequence of α values starting with a very small value (e.g., 10 −5 ) and increasing α by one order of magnitude each step.In the case of a valid adjoint and correct gradient, R will increase by two orders of magnitude step by step over a wide range of α values.An example of such a successful gradient test is shown in Table A1.The gradient test is done for each independent parameter individually.For each type of independent parameter, we chose an arbitrary position to apply a change of independent parameter by dp, normally 0.5.This kind of gradient test should be done for each type of independent parameter.For example, in Exp B, it should be done three times, just changing one position in the type of p e , with the other independent parameters kept unchanged.Only if all gradient tests are successful do we start the production run.

Fig. 3 .
Fig. 3. Variable resolution model grid.The green circle indicates the position of the data in Fig. 6.

Fig. 4 .
Fig. 4. Schematic diagram of tracers and processes included in model.The two simulated tracers phosphate and particular organic phosphorous are represented by blue squares and colored dots, with symbol size reflecting typical concentrations.The size of the black arrow indicates decreasing export flux with depth.Phosphate is transported by the ocean circulation (two big blue arrows) and also influenced by the seasonal changes of mixed layer depth.
Dd is taken from the WOA09 phosphate standard error or 0.03 µmol kg −1 , whichever is greater.σ C Pd is taken as 30 % of the satellite POC or 13 mg m −3 , whichever is greater.p c , p w , p e , p n , and p s are center, west, east, north, and south values,respectively.

Fig. 6 .
Fig. 6.Comparison of monthly phosphate and POC values in the central North Atlantic (position marked by green circle in Fig. 3).(a) WOA09 phosphate and satellite POC, (b) model-simulated phosphate for Exp A and Exp B together with WOA09 phosphate, (c) model-simulated POC for Exp A and Exp B together with satellite POC, and (d) monthly mixed layer depth (MLD, Monterey and Levitus, 1997).

Fig. 7 .
Fig. 7. Model-simulated surface water phosphate and model-data misfits for June and December.The left-hand column (a-d) is for Exp A, and the right-hand column (e-h) for Exp B.

Fig. 8 .
Fig. 8. Model-simulated surface water POC and model-data misfits for June and December.The left-hand column (a-d) is for Exp A, and the right-hand column (e-h) for Exp B.

Fig. 9 .
Fig. 9. Optimal POC export fields for June and December for Exp A (a and b) and Exp B (c and d).

Table 1 .
The parameters used in the model.
Schlitzer (2007)D budgets exist for all depths, while C P budgets only exist in the surface layer, the number of C D parameters is much larger than the number of C P parameters.The initial field of C D is taken fromSchlitzer (2007), while the C P field is initialized with zero uniformly.The time evolution of C D and C P is calculated by running a 10 yr model simulation.The model values for the last year are compared with phosphate and POP data month by month, and the observed model-data misfits are fed into the adjoint model described below.

Table 2 .
Description of individual terms in the cost function F .

Table 3 .
Values of total cost function and individual terms (units of 10 7 ) for Exp A and Exp B.

Table 4 .
Mean and RMS of phosphate and POC model-data misfits for Exp A and Exp B.

Table 5 .
Summary of Exp B sensitivity runs starting at different initial independent parameters.

Table A1 .
Example results of a successful gradient test.