Technical Note : Comparison of ensemble Kalman filter and variational approaches for CO 2 data assimilation

Data assimilation (DA) approaches, including variational and the ensemble Kalman filter methods, provide a computationally efficient framework for solving the CO 2 source–sink estimation problem. Unlike DA applications for weather prediction and constituent assimilation, however, the advantages and disadvantages of DA approaches for CO 2 flux estimation have not been extensively explored. In this study, we compare and assess estimates from two advanced DA approaches (an ensemble square root filter and a variational technique) using a batch inverse modeling setup as a benchmark, within the context of a simple one-dimensional advection–diffusion prototypical inverse problem that has been designed to capture the nuances of a real CO 2 flux estimation problem. Experiments are designed to identify the impact of the observational density, heterogeneity, and uncertainty, as well as operational constraints (i.e., ensemble size, number of descent iterations) on the DA estimates relative to the estimates from a batch inverse modeling scheme. No dynamical model is explicitly specified for the DA approaches to keep the problem setup analogous to a typical real CO2 flux estimation problem. Results demonstrate that the performance of the DA approaches depends on a complex interplay between the measurement network and the operational constraints. Overall, the variational approach (contingent on the availability of an adjoint transport model) more reliably captures the large-scale source–sink patterns. Conversely, the ensemble square root filter provides more realistic uncertainty estimates. Selection of one approach over the other must therefore be guided by the carbon science questions being asked and the operational constraints under which the approaches are being applied.


Introduction
Data assimilation (DA) is best known as a tool in numerical weather prediction (NWP; e.g., Swinbank, 2010) and has been applied to analyze complex data sets and estimate parameters in a variety of fields, including atmospheric constituent (e.g., Lahoz and Errera, 2010;Elbern et al., 2010), oceanographic (e.g., Haines, 2010), and land surface (e.g., Reichle, 2008;Houser et al., 2010) assimilation problems.In all such applications, a DA system aims to optimally combine the information from available observations with a prior model estimate (or the background derived from a model forecast) based on their respective uncertainty estimates.
DA approaches for estimating CO 2 fluxes aim to constrain the spatial and temporal distributions of CO 2 sources and sinks by integrating atmospheric, terrestrial and oceanic data together into a common analysis framework.CO 2 -DA applications (e.g., ensemble-filter-based approaches; Peters et al., 2005;Feng et al., 2009;Miyazaki et al., 2011;Chatterjee et al., 2012;Kang et al., 2012); variational based approaches (Rayner et al., 2005;Chevallier et al., 2005;Rödenbeck, 2005;Baker et al., 2006), and hybrid ensemble-variational approaches such as the maximum likelihood ensemble filter (Zupanski et al., 2007;Lokupitiya et al., 2008) have been in vogue for nearly a decade and are viewed as an alternative to more traditional batch inverse modeling schemes.Unlike these newer DA applications, which use a combination of numerical approximations and time-stepping approaches, the batch schemes directly solve the linear system of equations relating the fluxes and the atmospheric CO 2 observations in a single step.The DA approaches are attractive because of their computational efficiency (e.g., Rayner, 2010), but the impact of their underlying numerical approximations on the final estimates and their associated uncertainties is often unclear.Chatterjee et al. (2012) pointed out fundamental differences between the carbon flux estimation (i.e., the inverse framework) and the NWP/constituent (i.e., assimilation framework) problems -namely that (a) the performance of the DA approaches are not necessarily equivalent for the two frameworks and (b) that only under specific inversion scenarios are the DA approaches able to perform optimally.Differences between the two frameworks are mainly driven by the ill-conditioned and highly diffusive nature of the flux estimation problem, as well as the absence of an explicit dynamical model that can evolve a set of estimated fluxes forward in time.By propagating the state vector between different assimilation time steps, a dynamical model directly contributes to the growth of the eigenvalue spectrum of the state covariance matrix in certain preferred directions and the decay in others (Bengtsson et al., 2003;Furrer and Bengtsson, 2007).For the CO 2 inverse problem, however, the dynamics are embedded within the atmospheric transport and are not used explicitly to inform the temporal evolution of the state vector.A few recent studies have attempted to use an explicit dynamical flux model (e.g., Kuppel et al., 2013), but they note that model errors significantly reduce the performance of the inversion in terms of the quality of the estimated fluxes.Currently, in most carbon flux estimation studies, dynamical models are not used, resulting in a loss of potentially valuable information to the DA system.The absence of this information, coupled with the availability of only sparse observational data sets, may result in the DA approaches performing suboptimally.
The authors are not aware of any study specifically related to the CO 2 flux estimation problem that attempts to evaluate the relative performance of DA techniques.This is unlike the weather forecasting community, where several studies have evaluated the strengths and weaknesses of ensemble and variational approaches for different weather-related applications ranging from simple to chaotic nonlinear systems (e.g., Lorenc, 2003;Caya et al., 2005;Fertig et al., 2007;Kalnay et al., 2007;Liu et al., 2008;Whitaker et al., 2009;Buehner et al., 2010a, b;Jardak et al., 2010;Zhang et al., 2011; also see the special collection of papers on intercomparison at http:// journals.ametsoc.org/page/Ensemble_Kalman_Filter).Apart from NWP-related comparison studies, DA approaches have also been intercompared for chemical (e.g., Carmichael et al., 2008) and constituent (e.g., ozone - Wu et al., 2008) assimilation problems.These comparison studies cannot be used as a baseline, however, because of differences between the flux estimation and the NWP/constituent DA frameworks, as stated earlier.
The main purpose of this work is thus to fill this gap and build on the existing body of intercomparison studies from the perspective of the CO 2 flux estimation problem.Specifi-cally, this study aims to answer the following two questions: (1) What is the relative performance of two state-of-the-art DA approaches (ensemble square root filter, EnSRF (e.g., Whitaker and Hamill, 2002), and 4-dimensional variational, 4D-VAR (e.g., Talagrand, 2010) for solving the CO 2 inverse problem, and (2) how well can the DA approaches reproduce the flux estimates and associated uncertainties from a batch inverse modeling (BIM) scheme?
To facilitate the intercomparison, we consider here a onedimensional (1-D) passive tracer transport problem.Similar to previous studies (e.g., Liu and Rabier, 2002;Park and Kalnay, 2004), the 1-D framework allows us flexibility in setting up the problem because multiple experiments can be simulated in a computationally efficient manner.The low computational cost associated with the 1-D problem enables the implementation of a BIM approach in addition to the DA approaches.The DA estimates are thus compared both to the true signal and to the BIM estimates in order to isolate the degradation due to the underlying numerical approximations in the DA approaches.This study assesses whether these approximations limit the ability of the examined DA approaches to be used as suitable long-term replacements for the BIM approach under different inversion conditions.
When designing the 1-D problem, we focus on a framework that allows us to examine an underdetermined and finescale CO 2 flux estimation problem.This setup is necessary to mimic the challenges of a true CO 2 flux estimation problem in which atmospheric mixing coupled with the sparseness of available observations results in the inverse problem being highly underdetermined and ill posed.The underdetermined nature of the problem is accentuated by the need for estimating CO 2 fluxes at fine spatial and temporal scales, which is necessary to not only avoid spatiotemporal aggregation errors (e.g., Kaminski et al., 2001;Gourdji et al., 2012) but also to improve the understanding of the fine-scale processes driving the carbon cycle.This paradigm shift has brought about a computational bottleneck in solving the BIM problem, which requires the atmospheric transport model to be run either once per estimated flux region/period combination or once per observation if an adjoint to the transport model is available.This in turn has prompted the use of computationally efficient alternatives, such as DA approaches, in which the number of atmospheric transport model runs is proportional to the number of ensemble members (in the ensemble approach) or the number of descent iterations (in the variational approach), both of which are typically set to be orders of magnitude lower than the number of estimated parameters or available observations.Analogous to a real CO 2 flux estimation problem, no dynamical model is explicitly specified for solving the 1-D problem.
The experiments are specifically targeted to evaluate the impact of three factors on the two DA approaches: (a) the impact of the observational density and homogeneity, (b) the impact of the model-data mismatch covariance, and (c) the impact of the operational parameters of the DA system (i.e., ensemble size, number of iterations).While examining the first two factors, issues of sampling and convergence error are minimized by specifying a large number of ensemble members and descent iterations for the EnSRF and the 4D-VAR 1 approaches, respectively.More realistic operational constraints are subsequently imposed in a latter set of experiments to not only evaluate the fundamental differences between the two DA approaches but also the effect of the compromises necessary to make the algorithms practically applicable.This study represents the first comparison of the EnSRF and the 4D-VAR approaches for a flux estimation problem, and is expected to guide the development of future intercomparison experiments with real data, including satellite observations of atmospheric CO 2 .

Estimation methods
In a Bayesian framework, prior information and likelihood are expressed as probability density functions or pdfs (e.g., Enting, 2002).If the pdfs can be approximated as Gaussian, then maximizing the posterior probability of the state is equivalent to finding the minimum of a quadratic objective function: where s is a m × 1 vector of the discretized state (e.g., CO 2 flux), z is the n × 1 vector of observations (e.g., CO 2 observations), h is a forward model that is often a combination of an observation operator and an atmospheric transport model with embedded dynamics, R is the n × n model-data mismatch covariance matrix, s b is the m × 1 prior estimate of the state, Q b is the m × m error covariance matrix of the prior estimate s b , and the superscript T denotes the matrix transpose operation.Note that in the case of the atmospheric CO 2 inverse problem, h is linear and typically represented via its matrix form H (i.e., sensitivity matrix with dimensions n × m, or Jacobian matrix), which captures the sensitivity of the observations z to the fluxes s (i.e., H i,j = ∂z i ∂s j ).The inverse problem as formulated via Eq.( 1) determines a leastsquares fit of the state estimate to the observations and the 1 Typically in the DA community the term 4D-VAR is used to represent the three-dimensional space plus time.In this study, the variational approach is applied to a one-dimensional space plus time, which may suggest that the term "1 + 1D-VAR" may be more appropriate.Within the geophysical community, however, the term 1-D in 1D-VAR specifically refers to the vertical column, and is quite popular for radio occultation data (e.g., Eyre et al., 1993;Poli et al., 2002), total column water vapor (e.g., Marécal and Mahfouf, 2002;Bauer et al., 2006), cloud (e.g., Janiskova et al., 2012) assimilation etc.Since in the current study 1-D refers to a single dimension along the horizontal space and not necessarily the vertical column, we persist with usage of the term 4D-VAR rather than the term 1D-VAR.
prior, with the ultimate aim of estimating the posterior pdf of the true state.
The estimate of the pdf of the true state is obtained here via three approaches -BIM, EnSRF and 4D-VAR.For Gaussian pdfs and a linear model, BIM analytically solves the linear system of equations resulting from the minimization outlined in Eq. ( 1) to obtain the parameters of the posterior pdf.4D-VAR and EnSRF provide numerical approximations of this solution, and if perfectly implemented, they will yield the same solution, which in practice is never feasible due to operational constraints.The best estimate for EnSRF is defined as the mean estimate across an ensemble of source-sink realizations, and the posterior uncertainties are defined based on the ensemble spread at each location/time.The best estimate for 4D-VAR is defined as the maximum likelihood estimate that minimizes the misfit between a prior guess of the sourcesink estimate and the observations that are available over a given assimilation window.The posterior analysis covariances for 4D-VAR need to be calculated indirectly, however, for example via a Monte Carlo algorithm (e.g., Chevallier et al., 2007).In the Monte Carlo framework, different 4D-VAR estimates based on perturbations of the observations (using the specified model-data mismatch statistics), perturbations of the prior (using the specified prior error statistics), and combinations thereof.In practice, for a real CO 2 flux inversion problem it may not be computationally feasible to run 4D-VAR with a large number of perturbations.The computational efficiency of the one-dimensional framework in this study allows us to be generous with the number of perturbations (25 total) to assess the quality of the Monte Carloalgorithm-based error statistics relative to the uncertainty estimates obtained directly from the other two approaches.
A review of these approaches and their underlying mathematical framework is available in Appendix A, along with an exposition of the algorithmic choices necessary to adapt the DA approaches for solving the CO 2 flux estimation problem.

Problem description
A 1-D advection-diffusion problem of a passive tracer is used to emulate the CO 2 flux estimation problem.In the 1-D illustration, the passive tracer represents atmospheric CO 2 .Tracer fluxes get released from a series of locations over a finite duration and get transported by a tracer transport model that emulates the physics of advection and diffusion.No sink is specified, and there is therefore a gradual buildup of the passive tracer within the domain.Observations of the tracer are obtained at various locations and times within the domain.The locations and times of the observations as well as their precision can be regulated to simulate a variety of inversion scenarios.The inverse problem involves using the noisy tracer observations along with the transport information to infer the original tracer fluxes.
In the following description, the units of mass, length and time are reported as [M], [L], and [T ] to keep the problem The tracer observations (z) and the tracer fluxes (s) are related in the following fashion: where H is the sensitivity matrix that is generated using a 1-D tracer transport model as where x r and t r are the tracer flux release locations and times, x o and t o are the tracer observation locations and times, erfc represents the complementary error function.The tracer transport model embedded in Eq. ( 4) assumes conservation of mass and is based on a well-known one-dimensional analytical solution for a conservative tracer with a continuous source under steady-state conditions (e.g., Ogata and Banks, 1961;Runkel, 1996).The tracer observations obtained at a particular time step are sensitive to the tracer flux released at multiple previous time steps.Given that the total length of the domain is 300 [L] and the advection velocity is 50 [L T −1 ], the maximum residence time of the tracer within the domain is approximately 6 [T ].Based on the form of Eq. ( 2), however, the majority of fluxes occur at x r >50 [L] (Fig. 1b), and the typical residence time is therefore <5 [T ].This means that an observation taken at time t o [T ] provides information about the tracer flux approximately up to time t o − 5 [T ].In all subsequent experiments, the lag window size for the DA approaches is thus set to 5 [T ], such that this "long" lag window allows nearly all of the flux influence on observations to be represented within the DA approaches.The finite lag window also recognizes the operational limitations associated with the implementation of DA approaches in a real-world setting, and we foresee this as a realistic scenario for future global inversions (see Appendix A for a more detailed discussion).
The final piece of information necessary for setting up the inverse problem is the prior estimate (s b ) of the tracer flux and its error characteristics (Fig. 1b).s b is chosen here to be constant across all time periods: Its error covariance matrix Q b is based on an exponential decay model in space, with a correlation length (3l . The 1-D framework was designed to capture most of the characteristic features of the CO 2 flux estimation problem.For a real CO 2 inversion, units are most typically [µmol (m 2 s) −1 ] for s and σ Q , [ppm] for z and σ R , [ppm µmol −1 (m 2 s −1 )] for H, and [km] for l Q .

Experiments
Experiments are designed to explore the impact of three factors on the ability of the DA approaches to solve the inverse problem: (a) the observational density and homogeneity, (b) the model-data mismatch covariance, and (c) the operational constraints of the DA system (i.e., ensemble size, number of descent iterations).In all the experiments, the size of the state vector or the total number of fluxes to be inferred is 10 500 × 1 (i.e., 300 locations × 35 times).
The first set of experiments (Table 1 -experiments A through C) aims to investigate the effect of the density and spatiotemporal homogeneity of the observational network.Three different observational networks are designed (Fig. 2).In the first network configuration (denoted as REF -"reference observational set" as outlined in Sect.2.2), observations are obtained throughout the domain (x o = 1, 2,. . ., 299, 300 [L]) and for all 35 measurement times (t o = 1.5, 2.5,. . ., 34.5, 35.5 [T ]) (Fig. 2a).The total number of observations available is thus 10 500 (i.e., 300 locations × 35 times).In the second network configuration (denoted as HM -"homogeneous"), observations are obtained at 25 equally spaced locations within the 1-D domain (x o = 10, 22, 34. . . , 298 [L]) for all 35 time periods (t o = 1.5, 2.5,. . ., 34.5, 35.5 [T ]) (Fig. 2b).The total number of observations is thus reduced to 875 (i.e., 25 locations × 35 times).In the final configuration (denoted as HT -"heterogeneous"), observations are taken at 25 randomly selected locations for each measurement time (t o = 1.5, 2.5,. . ., 34.5, 35.5 [T ]) (Fig. 2c), and these locations vary from one time to the next.Similarly to HM, the total number of observations in HT is 875 (i.e., 25 locations × 35 times), but the observations are neither uniform in space nor consistent in time.Note that unlike REF, both HM and HT represent underdetermined inversion problems where the total number of observations is substantially lower than the number of unknowns in the state space to be estimated.In reality, the HT network configuration scheme is the closest to current CO 2 monitoring networks where different monitoring locations (ground based or remote sensing) can come online and go offline over different periods.
For all the three network configurations (experiments A through C), random errors with a variance of 10 [M 2 L −2 ] are added to the observations to represent measurement, transport, aggregation, and representation errors encountered in real applications.The model-data mismatch covariance matrix R (Eq. 1) has this same variance as its diagonal elements.In contrast to the prior error covariance, the errors  have no spatial or temporal correlation.Finally, all three approaches use the same numerical realization of errors, thus ensuring that they are solving the same inverse problem.
The second set of experiments (Table 1 -experiments AR through CR) examines the effect of the model-data mismatch variance on the best estimates and their associated uncertainties.For all the network configurations, the variance of the random errors is increased to 400 [M 2 L −2 ] with the diagonal values of the model-data mismatch covariance matrix R increased accordingly.All other parameters are kept the same as in the first set of experiments.
The third set of experiments (Table 1 -experiments AO through CO) explores the impact of operational constraints, which are always an important consideration in implementing a DA system.To minimize numerical approximations and avoid sampling or convergence errors, the ensemble size (for EnSRF) and the number of descent iterations (for 4D-VAR) for the first two sets of experiments (Table 1 -experiments A-C and experiments AR-CR) are set to 1000 and 250, respectively.The number of descent iterations is prescribed to be lower than the number of ensemble members, keeping in mind that 4D-VAR typically requires more model integrations (i.e., both forward and adjoint model run) than EnSRF.Given that it is not feasible to either run a large number of ensemble members or specify a large number of descent iterations for real atmospheric applications, these numbers are reduced to 100 ensemble members for EnSRF and 25 descent iterations for 4D-VAR in the third set of experiments.The noise added to the observations is kept the same as in the first set of experiments, namely 10 [M 2 L −2 ], to allow for a direct comparison with experiments A-C.

Results
In the sections that follow, results from the nine experiments are interpreted at both the native and aggregated spatial scales, and estimates from the EnSRF and the 4D-VAR approaches are compared both to the truth and to the estimates from the BIM approach.Taylor diagrams (Taylor, 2001) are used to assess the root-mean-square difference (RMSD) and the correlation coefficient (CC) between the flux estimates and the truth, as well as the standard deviation (SD) of the flux estimates and the truth.These metrics are calculated across 30 time periods (t r = 6, 7,. . ., 34, 35 [T ]) to be representative of the overall experiment after discarding the first five time periods as spin-up.

Impact of observational density and homogeneity
For the REF network (experiment A), all three approaches perform well in recovering the true flux (e.g., Fig. 3a), and in fitting the observations within the specified model-data mismatch errors (results not shown).For the sample time period presented in Fig. 3a, both the 4D-VAR and the En-SRF estimates capture the flux profile, including its large and small peaks.These results are typical of the performance of the three approaches across other estimation times.The performance across the full examined time period is summarized in Fig. 4a, where all three approaches show a high CC (∼0.97), low RMSD (∼0.3 [M L −1 T −1 ]), and standard deviations (∼1.5 [M L −1 T −1 ]) that are similar to that of the true fluxes.
The performance of all three approaches degrades as the observational density and homogeneity decrease in going from experiments A to C. This is evident by looking at Fig. 3b and c, where the estimates fail to capture the smaller  1.
double peak around 100-200 [L], and the Taylor diagrams in Fig. 4b and c show a corresponding drop in CC and an increase in RMSD.In general, for observations with spatially uncorrelated model-data mismatch errors such as those used here, decreasing the observational density is expected to decrease the analysis accuracy.The response of the two DA approaches mirrors the BIM approach in such cases, including the inference of an incorrect flux pattern for the HT network around 100 [L] in Fig. 3c.This result indicates that in the absence of operational constraints, best estimates from the DA approaches are consistent with the BIM estimate even for an underdetermined inverse problem.
In terms of the recovered posterior uncertainty estimates, the EnSRF uncertainty estimates are more consistent with the BIM uncertainty estimates relative to 4D-VAR (Fig. 3ac).For the REF observational configuration, the average ra-tios of the predicted posterior uncertainty of the individual flux estimates in EnSRF (σ ŝEnSRF ) and 4D-VAR (σ ŝ4D−VAR ) to those from BIM (σ ŝBIM ) are approximately 0.98 and 0.84, respectively; that is, on average, EnSRF and 4D-VAR underestimate the posterior uncertainties by 2 and 16 %, respectively, relative to BIM (Fig. 5a).As the observational density changes, EnSRF overestimates the uncertainty by 2 and 6 % for HM and HT, respectively, while 4D-VAR underestimates the uncertainty substantially, by 22 and 18 %, respectively (Fig. 5b, c).For all cases, the 4D-VAR uncertainties for individual locations/times over-or underestimates the BIM estimates even more substantially than evidenced by the average statistics, however, as seen by the spread in the histograms in Fig. 5.
In the EnSRF framework, the uncertainties are directly related to the ensemble spread.In the absence of a dynamical  1.For each experiment, statistics are calculated between the estimates and the true fluxes across all locations and all 30 time periods and are represented on a Taylor diagram.For each Taylor diagram, the true flux is represented by a point along the abscissa corresponding to the standard deviation of the true fluxes ("Truth").All other points ("BIM", "EnSRF", "4D-VAR"), which represent the estimated fluxes, are positioned such that their standard deviation is the radial distance from the origin, the correlation coefficient between the estimates and the truth is the cosine of the azimuthal angle, and the root-mean-square difference (RMSD) between the estimates and the truth is the distance to the observed point.In the limit of perfect agreement, these other points would coincide with "Truth" (i.e., RMSD = 0, CC = 1, and SD of the estimates would be the same as that of the truth).model, there is little source of variability for the ensemble to maintain a consistent spread.As observations are assimilated, the ensemble members tend to collapse to the ensemble mean and the adaptive inflation (see Appendix A.3) has to compensate for this degeneracy by inflating the ensemble spread.In the HT case, however, the inflation technique has a delayed response in adjusting to the heterogeneity in the observational network, as different observation locations come into and out of the network.For the adaptive inflation component to function well, we find it beneficial to have a consistent set of observations to maintain a reasonable ensemble spread (e.g., Chatterjee et al., 2012).It is worthwhile to mention here that the magnitudes of the inflation factors are very small in experiments A-C.This is not surprising given that a large number of ensemble members have been specified and the sampling error is hence quite low.
For 4D-VAR the posterior uncertainties are obtained via the Monte Carlo technique.For the time period shown in Fig. 3a-c, the 4D-VAR uncertainties not only underestimate the BIM uncertainties on average but they are also a lot more variable (i.e., too high or too low) for individual fluxes relative to BIM.We believe this heterogeneity in the posterior uncertainty estimates to be a result of the number of perturbations specified for the Monte Carlo technique.A sensitivity Fig. 5. Histogram showing the ratio of the estimated posterior uncertainties from the EnSRF (σ ŝEnSRF ) and the 4D-VAR (σ ŝ4D-VAR ) approaches to the posterior uncertainties from the BIM (σ ŝBIM ) approach.The ratios are calculated for each estimated flux time/location over all 30 time periods to be representative of the overall experiments outlined in Table 1.In the ideal case -that is, if the posterior uncertainty from a DA approach is equal to the posterior uncertainty from BIM -the histogram would be a single line centered around 1.
test in which the total number of perturbations is increased (or decreased) indicated that the heterogeneity of the uncertainty estimates obtained via the Monte Carlo technique decreased (or increased) correspondingly.When averaged over all the time periods (Fig. 5a-c), the posterior uncertainty estimates clearly underestimate the BIM uncertainty.Even then, the large spread in the histograms shown in Fig. 5a-c reinforces our earlier conclusions that the uncertainty estimates for individual fluxes over/underestimate the BIM uncertainty estimates substantially.
Overall, we find that both 4D-VAR and EnSRF can reproduce the performance of BIM in terms of the best estimates of the fluxes for all three observational network configurations.Even though small discrepancies are noticeable, the impact of a sparse and/or heterogeneous observational net-work is similar for the DA approaches compared to the BIM approach.Both the DA approaches have some difficulty in reproducing the BIM posterior uncertainty estimates, albeit for different reasons.In the absence of any operational constraints, however, EnSRF provides more realistic and useful uncertainty bounds than 4D-VAR.

Impact of model-data mismatch covariance
For all the network configurations, the quality of the estimates degrades when a higher model-data mismatch error is prescribed (comparing Fig. 3ar-cr with Fig. 3a-c), although the heterogeneous network estimates show the most pronounced degradation.An increase in σ 2 R leads to higher uncertainties for the estimates, indicative of the decreased confidence in the analysis.Analogous to the first set of experiments, EnSRF and 4D-VAR best estimates respond similarly to BIM when the model-data mismatch covariance changes, and both the approaches track the BIM best estimates quite well for all the three experiments.Figure 4ar-cr confirm that the best estimates from all the three approaches have a lower CC, higher RMSD, and lower SD when compared to Fig. 4a-c.
The standard deviation of the flux estimates change considerably between Fig. 4ar (∼1.45-1.5 [M L −1 T −1 ]) and 4cr (∼0.94-1.00[M L −1 T −1 ]) for the three approaches.Increasing the σ 2 R to 400 [M 2 L −2 ] results in the analysis rejecting the information from the observations and giving more weight to the prior, yielding overly smooth a posteriori estimates.A typical example of this is seen by comparing the estimated peak around 50-100 [L] in Fig. 3c and cr.Estimates in both these panels are based on the same observational network but the estimates in Fig. 3cr do not capture the amplitude of the two large peaks in the true flux signal.
For experiments AR-CR, the posterior uncertainty estimates for all the three approaches are higher compared to experiments A-C, as expected due to the higher prescribed model-data mismatch error.Similarly to experiments A-C, the 4D-VAR uncertainty estimates for individual locations/times are too variable relative to BIM (Fig. 5).Averaged over time and space, the 4D-VAR uncertainty estimates underestimate the BIM uncertainty estimates by approximately 25 % (Fig. 5ar-cr).Thus, even though the 4D-VAR uncertainty estimates for experiments AR-CR are higher than the corresponding uncertainty estimates for experiments A-C, they fail to capture the full magnitude of the BIM uncertainty estimates.This makes intuitive sense due to the indirect approach adopted for generating the 4D-VAR uncertainty estimates.Conversely, as the observational network becomes sparser and more heterogeneous, the EnSRF slightly overestimates the BIM average uncertainties by 3 % (HM; Fig. 5br) and 5 % (HT; Fig. 5cr), while it underestimates the uncertainty by only 1 % for the reference network (Fig. 5ar).The EnSRF uncertainty estimates for individual locations/times are more closely distributed around the BIM estimates (Fig. 5).The better performance of EnSRF in terms of the uncertainty estimation can be directly related to the ensemble spread.Relative to experiments A-C, when the prescribed model-data mismatch error is high in experiments AR-CR, the initial ensemble spread is reduced by a lower amount as observations are now being given less weight and hence have lower impact on the ensemble spread.Consequently, the ensemble members maintain a large spread throughout the analysis and results in large posterior uncertainty estimates that are more realistic relative to 4D-VAR.
Experiments AR-CR reconfirm that in the absence of operational limitations, an increase (or decrease) in the modeldata mismatch covariance does not affect the ability of the DA approaches to reproduce the BIM best estimates.Similar to experiments A-C, in terms of the posterior uncertainty estimates, however, (a) both the DA approaches are less skilled at reproducing the uncertainty estimates from BIM, and (b) the uncertainty estimates from 4D-VAR severely underestimate the BIM uncertainty estimates and are less realistic than the EnSRF uncertainty estimates.

Impact of operational constraints
Operational constraints hinder the performance of the DA approaches, and the impact is further intensified as the observational network becomes more heterogeneous.
For 4D-VAR, an inadequate number of iterations may lead to a failure to find the minimum of the quadratic objective function (convergence results not shown here).When the observational network is heterogeneous, the minimization has even more difficulty in finding the path towards the minimum.Thus, comparing Fig. 4ao-co, the 4D-VAR estimates diverge from the BIM estimates for the HT network configuration.In general, we find that for the HT network, 4D-VAR needs approximately 50 iterations to converge completely for the case studies presented here.Conversely, for the REF and the HM network, 4D-VAR requires only approximately 20 to 30 iterations to reach full convergence.For all the three experiments, however, the value of the objective function is reduced relative to that for the prior fluxes, indicating an improvement over the prior estimates.As pointed out by Rödenbeck (2005), the minimization determines the large-scale gradient in the initial iterations, while in subsequent iterations fine-scale tuning is performed to capture the optimum.By artificially limiting the number of iterations in experiments AO-CO, the ability of 4D-VAR to make small-scale adjustments is hindered, which manifests itself clearly in experiment CO (panel co in Figs. 3, 4 and 5).
For EnSRF, the degradation is attributable to sampling error caused by the limited ensemble size.This reduces the estimation accuracy (both flux estimates and their uncertainties) and makes the filter sensitive to the observational density.Note that a Schur-based localization scheme was implemented for EnSRF (see Appendix A3).Since the localization length scale is dependent on the ensemble size, when an ensemble size of 1000 was used (experiments A-C, or AR-CR), a long localization length scale of 90-120 [L] can be used.The localization length scale is determined subjectively based on sensitivity tests, and hence a range of values (i.e., 90-120 [L]) is acceptable within which the EnSRF estimates are not contaminated by spurious noise.Reducing the ensemble size to 100 requires the use of a shorter localization length scale.It was found beneficial to have different length scales for the different observational networks, namely 10-30 [L] for the  for the sparser networks.Specifying a longer localization length scale than 30 [L] for the REF network led to a divergence of the En-SRF system.In this case, the spurious noise in the ensemble outweighs the positive impact of the observations.The complex interplay between the ensemble size and the observational density makes it difficult to identify a mathematical or Atmos.Chem.Phys., 13, 11643-11660, 2013 www.atmos-chem-phys.net/13/11643/2013/physical basis for selecting an appropriate localization length scale.We refer the reader to Chatterjee et al. (2012), as well as the sensitivity tests presented therein, for a more detailed discussion of the role of localization for the CO 2 source-sink estimation problem.Furthermore, when operational constraints are imposed, the posterior flux uncertainty estimates obtained using both EnSRF and 4D-VAR are unable to reproduce those from BIM.The degree to which posterior uncertainties from both DA approaches over/underestimate the BIM uncertainties for individual locations and times increases.This is clearly evident in Figs. 3 and 5, panels ao-co.Averaged over all time periods and locations, the posterior uncertainty estimates from EnSRF are still closer to the BIM uncertainty estimates, being within 10 % of the averaged BIM uncertainties for experiments AO, BO, and CO, whereas the average 4D-VAR uncertainties underestimate the BIM uncertainties by up to 25 %.The degradation in the uncertainty estimates provided by 4D-VAR and EnSRF is due to different reasons.While for EnSRF the large sampling error plays a dominant role, for 4D-VAR the perturbations in the Monte Carlo technique are unable to capture the true range of the posterior uncertainties.The impact for 4D-VAR is accentuated for experiments BO and CO, where the sparse network exacerbates the need for more iterations to reach convergence.
An important caveat here is that the results for both the DA approaches could potentially be improved through further tuning of each algorithm.For example, the implementation of more sophisticated algorithms to precondition and obtain faster convergence, or stronger localization schemes to dampen the spurious noise in the ensemble members, might provide slightly different responses and reduce the error incurred due to the numerical approximations.In spite of having state-of-the-art algorithms, however, once the underlying numerical approximations come into play, (a) EnSRF fails to reproduce the BIM best estimates, with the EnSRF performance decreasing as the observational network becomes sparser and more heterogeneous, and (b) 4D-VAR also fails to reproduce the BIM best estimates when the observational network is heterogeneous but still performs better than En-SRF.The better performance of 4D-VAR in terms of the flux estimates is offset by the fact that EnSRF provides more realistic uncertainty bounds on the recovered flux estimates.
The DA approaches are particularly sensitive to the information flow from the observations because of the lack of an explicit dynamical model.As discussed in Sect. 1, a dynamical model adds to the information content of the system, albeit at the cost of additional model errors that must be taken into account.Identifying/developing an appropriate dynamical model relevant to the CO 2 flux estimation problem, and subsequently repeating experiments such as those presented here, may further inform the assessment of the interplay between the operational constraints and the observational network.

Examining results at aggregated spatial and temporal scales
In the previous sections, the performance of the DA and the BIM approaches were analyzed and reported at the native estimation scales (both in space and time).In contemporary CO 2 inversion studies, the fine-scale estimates and their uncertainties are typically averaged a posteriori in space and/or in time to coarser scales (e.g., daily grid-scale fluxes averaged to monthly continental scales) for additional interpretation.
To represent this process, the domain used here is divided into two subregions, 1-150 [L] and 151-300 [L], and the true flux and the flux estimates are aggregated over each of these areas and examined across time.This is qualitatively analogous to aggregating fluxes a posteriori to "large regions" (e.g., biomes, continents) within the inversion domain.Because the true flux differs between these two subregions, and these differences themselves vary in time, it is possible to examine the ability of the various approaches to capture these spatiotemporal variations.Figure 6 presents the comparison at the aggregated scales in the form of a Taylor diagram.For all the 9 experiments, 4D-VAR is able to match the temporal variation of the spatially aggregated BIM estimates better than EnSRF.Even when the number of descent iterations is reduced (Fig. 6ao-co), the differences between the BIM and the 4D-VAR best estimates are negligible at the examined spatially aggregated scales.Comparing Figs.4co and 6co demonstrates that the differences observed at fine scales are substantially reduced when aggregating the estimates to a coarser resolution.
Because there is no explicit dynamical model to evolve the information between assimilation time periods, the EnSRF estimates are always contaminated with small-scale sampling errors, but these errors partially cancel out when the estimates are spatially aggregated.Overall, the EnSRF estimates still, however, exhibit more spurious variability relative to the 4D-VAR estimates.Especially when the number of ensemble members is reduced (Fig. 6ao-co), both the sampling error as well as the observational density and heterogeneity start to play a role, leading to less reliable estimates at aggregated scales.For example, in Fig. 6ao and 6co, the CC between the spatially aggregated EnSRF and the true flux estimates drops from 0.99 to 0.90, and the RMSD increases from 0.03 to 0.21 [M L −1 T −1 ].The standard deviation of these spatially aggregated EnSRF estimates also increases from 0.37 to 0.40 [M L −1 T −1 ], leading to an overestimation relative to the true flux, which has a standard deviation of 0.36 Analysis at aggregated scales demonstrates that when operational constraints are not imposed, both the DA approaches provide aggregated estimates that are close to the aggregated BIM estimates.Even when operational constraints are imposed, the aggregated flux estimates for En-SRF and 4D-VAR (Fig. 6co) have substantially higher CC A. Chatterjee and A. M. Michalak: Technical Note: Inter-comparison of EnKF and 4D-VAR for CO 2 -DA Fig. 6.Performance of BIM, EnSRF, and 4D-VAR at spatially aggregated scales for the different experiments outlined in Table 1.The details of the plot are as described in the caption of Fig. 4. and lower RMSD than the corresponding estimates at the fine scale (Fig. 4co).This is encouraging from the perspective of a real CO 2 flux estimation problem, as it implies that the DA approaches may provide reliable flux estimates at aggregated scales, even under circumstances when their performance at fine scales is compromised.

Discussion
Overall, the choice between 4D-VAR and EnSRF approaches for the CO 2 flux estimation problem should be based on the carbon science questions being targeted, as well as the tradeoff between the impact of incomplete convergence of the minimization algorithm (for 4D-VAR) and the impact of sampling error (for EnSRF) on the estimated fluxes and their uncertainties.
When an adjoint for the atmospheric transport model is available, and if the best estimates are the primary target, 4D-VAR outperforms EnSRF, especially when the measurement network is heterogeneous.With a small number of iterations, 4D-VAR may not converge to the minimum, but it still reliably captures the majority of the large-scale features, with the final estimates close to those expected from a BIM solution.The main disadvantages of the 4D-VAR approach, on the other hand, are its more cumbersome implementation and its inability to provide explicit estimates of the analysis error, which can instead be obtained indirectly via the Monte Carlo technique.In addition, the 4D-VAR approach may potentially be more advantageous due to its ability to more easily account for correlated model-data mismatch errors (O.Talagrand, personal communication, 2012).Although this specific issue has not been examined in this study, it is worth keeping in mind with the increasing use of satellite-based CO 2 measurements.Recent work by Brankart et al. (2009) has demonstrated techniques to cope with such correlations in an ensemble filter setting as well.Such techniques, however, are numerically efficient only for certain types of error correlation structures.
In the absence of correlated model-data mismatch errors, however, the serial EnSRF is easier to implement than a 4D-VAR system, and does not require the development and maintenance of an adjoint model.Due to restrictions on the size of the ensemble, however, it is necessary to adapt and tune ancillary algorithms such as localization and inflation, which improve the ensemble approximations to the full-rank Kalman filter.Unlike adaptive inflation, the choices to be made in implementing a localization scheme remain highly subjective.Experiments in this study show, for example, that the localization length scale is dependent on both the ensemble size and the observational density.Will increasing volumes of observations push us towards specifying shorter localization length scales?If so, what is the limit beyond which decreasing the localization length scale may actually degrade the analysis?It is necessary to identify more rigorously a basis for selecting the localization parameters (e.g., Anderson, 2012) or adaptive approaches that may be less sensitive to variations in the observational network.
One critical requirement of an inverse problem is to obtain reliable second-order statistical moments for the estimated system states.Even though we have demonstrated that posterior error statistics can be obtained for both EnSRF (directly) and 4D-VAR (indirectly via a Monte Carlo technique), the EnSRF uncertainty estimates reproduce the BIM uncertainty estimates consistently better.The indirect techniques necessary for obtaining posterior error statistics for 4D-VAR has important caveats, as demonstrated by the large over/underestimation of uncertainties for individual flux locations and times, as well as the underestimation of the posterior uncertainties relative to BIM when averaged across all locations and times.Therefore, EnSRF is more desirable for attribution purposes, wherein source-sink estimates with realistic confidence bounds can be used to gain a better understanding of the mechanistic processes driving the carbon cycle or to reconcile estimates from top-down and bottom-up approaches.
With both 4D-VAR and EnSRF, there is a direct tradeoff between computational savings and estimation accuracy, which is intensified when solving an underdetermined problem with a heterogeneous observational network.For largescale flux estimation problems, operational constraints will always exist, as will scarce and inconsistent observations, transport model biases and uncertainties (thus further limiting the use of available observations), etc.The HT scheme with a limited number of ensemble members/descent iterations (panel co in Figs. 3, 4 and 5) serves as the closest analogue to a real inversion problem.Even if we account for the increase in remote-sensing measurements of CO 2 , the observational network is going to be a complex hybrid between the REF and the HT scheme.In this scenario, the accuracy and precision of either of the DA approaches will be compromised relative to BIM, as demonstrated in this study.

Conclusions
We present a comparative assessment of two advanced DA approaches with the BIM approach for the atmospheric CO 2 inversion problem.The performance of the DA approaches is found to depend on a complex interplay between the underlying numerical approximations and the information available from the observations.Overall, the 4D-VAR scheme is found to be more robust for obtaining the best estimates, while the EnSRF scheme has the advantage of providing direct estimates of analysis error that are more realistic than those operationally feasible for 4D-VAR.The relative performance of 4D-VAR and EnSRF best estimates, when a large and homogeneous set of observations is available, is consistent with the conclusions obtained from intercomparison studies carried out by other DA communities.The sensitivity of the approaches to the observational scheme in the absence of an explicit dynamical model, and specifically for solving an underdetermined inverse problem, however, had not previously been thoroughly explored.Beyond CO 2 source-sink estimation problems, the conclusions of this study are therefore also relevant to other DA problems where a dynamical model is lacking.
The sensitivity experiments demonstrate that when a large number of ensemble members or descent iterations is specified, the best estimates obtained from state-of-the-art implementations of the 4D-VAR and the EnSRF approaches are similar to those from BIM, irrespective of the observational characteristics.Even under these optimal conditions, however, the uncertainty estimates from 4D-VAR are unable to reproduce those from BIM.When operational constraints are imposed, both the characteristics of the observational network and the numerical approximations play a greater role in differentiating the performance of the two DA approaches.Because such operational constraints are always present for real CO 2 source-sink estimation problems, the choice of an approach should be based on (a) the carbon science questions being targeted (i.e., whether the science calls only for best estimates of fluxes, or for flux estimates with realistic uncertainties) and (b) the inversion conditions under which they are being applied (i.e., characteristics of the observational network, such as data density and heterogeneity).
Recently, 4D-VAR and EnSRF have begun to influence each other's development, and EnVar DA approaches (Lorenc, 2013) may take center stage for solving the CO 2 flux problem in the future.In addition, the emerging focus on the assimilation of satellite CO 2 observations will require additional intercomparisons that go beyond the scope of the work presented here, including the consideration of temporal error correlations in the atmospheric CO 2 observations.In the BIM approach (e.g., Enting, 2002), the analytical solution for the a posteriori estimate and the associated covariances of the objective function (Eq. 1) are given by where ŝa is the posterior best estimate of the state and Q a is the a posteriori covariance of the recovered best estimate.
The diagonal elements of Q a represent the predicted error variance (σ 2 ŝ ) of individual elements in ŝa .As stated in Sect. 1, for CO 2 inversion studies, the generation of the matrix H requires an atmospheric transport model to be run either once per estimated element of the state vector, or once per observation.The large number of model runs ultimately makes the BIM approach computationally intractable for solving very large-scale problems.

A2 Four-dimensional variational (4D-VAR)
In 4D-VAR, the fluxes ŝa that minimize the objective function in Eq. ( 1) are sought iteratively by minimizing the misfit between a feasible state trajectory and the observations that are available over a given assimilation window.The overall approximation lies in the fact that the minimization can be stopped by artificially limiting the number of iterations or by requiring that the norm of the gradient decreases by a predefined amount during each iteration.Most minimization schemes rely on the availability of the gradient of the objective function with respect to the state (or control vector in 4D-VAR terms), Instead of analytically calculating the gradient, the adjoint of the forward transport model is used to compute the term σ 2 ŝ directly, which is then added to ŝ in Eq. (A4) above.While a variety of minimization schemes may be used (e.g., conjugate gradient or BFGS, e.g., Nocedal and Wright, 2006), in this study, we use a Lanczos minimizer algorithm (e.g., Nocedal and Wright, 2006), which produces similar results as the conjugate-gradient technique in terms of the reduction of the cost function and gradient norm but converges substantially faster.Additionally, a new variable ( ) is defined to precondition the minimization as where now becomes the control variable with respect to which the objective function is minimized instead of s directly.The optimal preconditioning matrix to reduce the number of iterations required to solve the minimization problem is the inverse Hessian (Axelsson and Barker, 2001).For a real high-dimensional CO 2 inversion problem, however, the size and structure of the prior covariance matrix may constitute a significant impediment for calculating the inverse Hessian.Any algebraic manipulations, such as taking inverses or calculating the square roots for preconditioning purposes, become operationally cumbersome and thus matrix manipulation techniques (e.g., Yadav and Michalak, 2013) are necessary to sidestep these operational challenges.
The main caveat with the variational approaches, such as 4D-VAR, is that a direct estimate of the analysis error is not available (no clear analogue of Eq.A2).Mathematically this can be obtained from the inverse of the Hessian (e.g., Le Dimet et al., 2002;Rödenbeck, 2005;Meirink et al., 2008) but operational challenges restrict the calculation and storage of the Hessian for high-dimensional problems.Recent applications of 4D-VAR for NWP problems have shown that computationally efficient alternatives do exist (e.g., Cheng et al., 2010;Gejadze et al., 2013).Although some of these techniques are more suited for nonlinear problems (e.g., Gejadze et al., 2013), they retain potential applicability to the CO 2 flux estimation problem as well.In this study, we use a Monte Carlo technique (e.g., Chevallier et al., 2007) where both the observations and the prior are perturbed multiple times with the prespecified model-data mismatch error statistics and prior error statistics, respectively.
4D-VAR is implemented with successive overlapping windows to (a) account for the long residence times of CO 2 in the atmosphere and (b) to avoid the operational cost associated with calculating the inverse of the full prior covariance matrix Q b in Eqs.(A4) and (A5).In order to determine the length of the moving window, however, one needs to keep in mind two primary factors: (a) observations typically inform fluxes over a finite period of time (e.g., 4-to 6-month time frame based on the analysis reported in Bruhwiler et al., 2005), and (b) temporal correlation in the prior flux errors (e.g., for the prior land flux errors Chevallier et al. (2012) has suggested the temporal correlation to be strongly positive for lags <85 days and mildly positive for lags >275 days).Given a sufficiently long window the 4D-VAR approach using a moving window will emulate the 4D-VAR approach using a single window, while at the same time providing substantial computational savings.In the absence of a dynamical model, this implementation of 4D-VAR becomes similar to the FGAT-3DVAR (Massart et al., 2010) variant occasionally used within the NWP community.

A3 Ensemble square root filter (EnSRF)
In the ensemble filter approach, the key innovation is to work in a reduced subspace of the error covariance matrices.
Observations are assimilated to update the ensemble representation of the error covariance matrices.The optimal analysis states and an estimate of the analyses error are determined in a similar fashion to Eqs. (A1) and (A2), but the calculation of Q b H T and HQ b H T is approximated by running transport model using the ensemble members directly: While several variants of the ensemble filter approach exist, here we have used a serial ensemble square root filter (Whitaker and Hamill, 2002) implemented in a fixed-lag smoother form (e.g., Whitaker and Compo, 2002;Chatterjee et al., 2012).Similarly to an ensemble square root filter, the ensemble smoother uses Monte Carlo estimates of the error covariances to compute a Kalman smoother gain matrix.This is applied iteratively to a time series of observations, where the analysis at the first time step is equivalent to an ensemble square root filter analysis as it only utilizes observations taken up to and including the analysis time.All subsequent time steps utilize observations taken a number of observing times past the analysis time.The localization scheme (i.e., to cut down spurious noise in the ensemble members) is based on Houtekamer and Mitchell (2001) using a fifth-order Gaspari-Cohn function (Gaspari and Cohn, 1999), while the adaptive inflation algorithm (i.e., to counter spurious variance deficiency among the ensemble members) is based on Anderson (2009).Implementation of these algorithms within an ensemble smoother framework is described in further detail in Chatterjee et al. (2012).Despite these ancillary algorithms, the overall implementation of EnSRF is quite simple and computationally efficient.

A4 Sensitivity tests and additional approaches considered
The setup of 4D-VAR (with overlapping time windows) and EnSRF (expressed as a fixed-lag smoother) reflect state-ofthe-art implementations of these two DA approaches, keeping in mind the nature of the atmospheric CO 2 process, the fact that CO 2 observations only provide significant information about fluxes over a finite preceding time window, and the operational constraints for estimating CO 2 fluxes at high spatiotemporal scales (e.g., spatial ∼1 • , temporal ∼daily).
Note that we have implemented both approaches with a single long window as a sensitivity test.Other than an increase in the computational cost, the conclusions regarding the performance of the approaches relative to the BIM approach are the same as those reported in the manuscript.We have also tested other varieties of the ensemble filter (EnKF with perturbed observations (Evensen, 2003) and the variational approach (PSAS; Courtier, 1997) and found that the overall conclusions from the presented experiments remain con-sistent across these algorithmic choices.We encourage the reader to look at Lorenc (2003) and Nichols (2010) (and references therein) for a more detailed discussion on these DA approaches.

Fig. 1 .
Fig. 1.(a) Spatiotemporal variability of the tracer flux.(b) True flux profile for a particular time period corresponding to the dashed white line in (a).Also shown in (b) is the prior estimate of the tracer flux profile with its ±2σ Q prior uncertainty (dashed lines).

Fig. 2 .
Fig. 2. Observations of the tracer obtained from the three network configurations -REF (a), HM (b), and HT (c).Note that going from the REF to the HM and the HT networks, the total number of observations decreases by a factor of 12, whereas in going from the HM to the HT network, the observational network becomes more heterogeneous in space and time.

Fig. 3 .
Fig. 3. Example of estimated tracer fluxes (lines) and associated ± 2σ ŝ uncertainties (shaded areas) for the different approaches assessed in this study.All values are shown for the 25th time period, which is representative of the observed performance over other time periods.The panel titles correspond to the different experiments outlined in Table1.

Fig. 4 .
Fig.4.Performance of the BIM, the EnSRF, and the 4D-VAR approaches for the different experiments outlined in Table1.For each experiment, statistics are calculated between the estimates and the true fluxes across all locations and all 30 time periods and are represented on a Taylor diagram.For each Taylor diagram, the true flux is represented by a point along the abscissa corresponding to the standard deviation of the true fluxes ("Truth").All other points ("BIM", "EnSRF", "4D-VAR"), which represent the estimated fluxes, are positioned such that their standard deviation is the radial distance from the origin, the correlation coefficient between the estimates and the truth is the cosine of the azimuthal angle, and the root-mean-square difference (RMSD) between the estimates and the truth is the distance to the observed point.In the limit of perfect agreement, these other points would coincide with "Truth" (i.e., RMSD = 0, CC = 1, and SD of the estimates would be the same as that of the truth).

Table 1 .
Summary of the experiments outlined in Sect.2.3.The following parameters are held constant for all the experiments in this study: s