Geoscientific Model Development Simulation of variability in atmospheric carbon dioxide using a global coupled Eulerian – Lagrangian transport model

This study assesses the advantages of using a coupled atmospheric-tracer transport model, comprising a global Eulerian model and a global Lagrangian particle dispersion model, to improve the reproducibility of tracer-gas variations affected by the near-field surface emissions and transport around observation sites. The ability to resolve variability in atmospheric composition on an hourly time-scale and a spatial scale of several kilometers would be beneficial for analyzing data from continuous ground-based monitoring and from upcoming space-based observations. The coupled model yields an increase in the horizontal resolution of transport and fluxes, and has been tested in regional-scale studies of atmospheric chemistry. By applying the Lagrangian component to the global domain, we extend this approach to the global scale, thereby enabling computationally efficient global inverse modeling and data assimilation. To validate the coupled model, we compare model-simulated CO 2 concentrations with continuous observations at three sites: two operated by the National Oceanic and Atmospheric Administration, USA, and one operated by the National Institute for Environmental Studies, Japan. As the goal of this study is limited to introducing the new modeling approach, we selected a transport simulation at these three sites to demonstrate how the model may perform at various geographical areas. The coupled model provides improved agreement between modeled and observed CO 2 concentrations in comparison to the Eulerian model. In an area where variability in CO2 concentration is dominated by a fossil fuel signal, the correlation coefficient between modeled and observed concentrations increases by between 0.05 to 0.1 from the original values of 0.5–0.6 achieved with the Eulerian model. Correspondence to: Y. Koyama (koyama.yuji@nies.go.jp)


Introduction
This study reports on the development and validation of a coupled atmospheric-tracer transport model, which comprises a global Eulerian grid model combined with a Lagrangian particle dispersion model (LPDM).In principle, LPDMs such as FLEXPART (Stohl et al., 1998) and STILT (Lin et al., 2003) can simulate highly resolved observation footprints, and also allow for a degree of flexibility in related decision-making because they calculate the trajectory of each particle individually.In addition, LPDMs overcome the problem of numerical diffusion that is associated with Eulerian models.LPDMs therefore avoid transport errors caused by numerical diffusion and non-monotonic behavior, which is inevitable in higher-order schemes (Godunov, 1959).Thompson (1971) developed an early model of particle diffusion that provided the foundation for subsequent LPDMs.Although LPDMs have been under development for several decades, they have only recently become popular as a tool for investigating the transport processes associated with atmospheric events.This slow adoption of LPDMs has occurred because such models require abundant computer resources: to accurately describe a multi-day transport episode requires simulations of the separate trajectories of a large number of particles over multiple days.Since the 1990s, many studies have, however, employed LPDMs (e.g., Zannetti, 1992;Uliasz, 1994;Rodean, 1996;Wilson and Sawford, 1996).
More recently, LPDMs have been used for various research applications at a regional scale (Lin et al., 2003(Lin et al., , 2006(Lin et al., , 2007;;Paris et al., 2010;Warneke et al., 2009).When applying an LPDM for global analysis, a long-term integration is required (see Stohl et al., 2010), although this approach is not as computationally efficient as that of Eulerian models.This problem may be overcome by using a coupled model as the resolution of the flux data for an LPDM can be selected Y. Koyama et al.: Simulation of variability in atmospheric carbon dioxide independently from the resolution of the simulation grid for a Eulerian model and its flux data.Several studies have investigated the coupling of Lagrangian and Eulerian models at a regional scale (Vermeulen et al., 2006;Trusilova et al., 2010).The benefits of coupled transport modeling at regional and global scales for locations downwind of highly heterogeneous surface fluxes can be fully realized by using emission inventories at resolutions of either 10 km (Olivier et al., 2008;Gurney et al., 2009) or 1 km (Kannari et al., 2007;Oda and Maksyutov, 2011).Benefits are also seen when coupled transport models are used for coastal background monitoring sites that are influenced by the large variations in composition of air arriving from the land and from the ocean (Ramonet and Monfray, 1996).Hence, a coupled model can calculate tracer-gas concentrations at any given global location, with high spatial-resolution surface fluxes requiring less computation time than conventional Eulerian models.However, the computational efficiency of a coupled model may be reduced in the case of simulations using very large observational datasets, such as those from satellite observations.
In this study, we expand the application of a coupled model to the global scale by coupling a global LPDM and a Eulerian model at a time boundary, as opposed to the domain boundary used in regional modeling (Trusilova et al., 2010).Hence, the spatial locations of coupling are not limited to a certain regional domain.Also, a global coupled model setup makes it possible to implement a single-stage inversion and data assimilation schemes, as opposed to the two-stage approach proposed for nested model setups (Peylin et al., 2005;Rödenbeck et al., 2009).Combining Eulerian and Lagrangian models offer the opportunity of constructing a cost-efficient, high-resolution surface flux data assimilation system.The high-resolution Lagrangian model requires just one run for every observation (e.g., Stohl et al., 2010) regardless of number of iterations, whereas the Eulerian component has lower resolution and requires forward and adjoint model runs for every iteration before convergence (Chevallier et al., 2007).
To develop the coupled model, we employ the National Institute for Environmental Studies transport model (NIES TM) (Maksyutov et al., 2008) as a Eulerian model, and use FLEXPART (version 3.2) (Stohl et al., 1998) as an LPDM.To validate the coupled model, we compare modeled CO 2 concentrations with the National Oceanic and Atmospheric Administration's (NOAA) continuous measurements at Barrow, USA (BRW; 71.32 • N, 156.61 • W) and at Samoa (SMO; 14.25 • S, 170.56 • E) (Thoning et al., 2007).We use data that have passed through three quality-control flags, as adopted by NOAA/ESRL (Earth System Research Laboratory).We also compare modeled CO 2 concentrations with continuous observations obtained by NIES at Hateruma station (HAT; 24.05 • N, 123.80 • E).This station is situated on the eastern end of Hateruma Island (Mukai et al., 2001), which is a small island (12.5 km 2 ) located at the southwestern end of the Japanese Archipelago, 220 km east of Taiwan.This site is characterized by southeasterly winds during the summer and northwesterly winds during the winter, and wintertime pollution events are common as a result of emissions transported from the Asian continent.
The main purpose of this study and of coupling a LPDM with a global Eulerian model is to understand whether tracer variations are dominated by the surface flux around the selected sites.In this instance, the reproducibility of short-term synoptic-scale variability is more important than contributions from large-scale latitudinal gradients provided by the Eulerian model alone.In the Northern Hemisphere in particular, variations in CO 2 concentrations are strongly affected by large-scale variability, with seasonal variations primarily relating to respiration and photosynthesis processes.It is therefore useful to subtract the seasonal influence when assessing the capabilities of the coupled model.Synoptic-scale variability was also extracted from simulated and observed CO 2 concentrations by subtracting the running average, and correlation coefficients and variance ratios were then calculated between simulated and observed data.

Materials and methods
In the coupled tracer-transport algorithm employed in this study, CO 2 transport is divided into two time periods.During the first stage, CO 2 concentrations are simulated by the Eulerian model on a global scale at a medium resolution.The LPDM is then used in the second stage to simulate transport with corresponding fluxes at a higher resolution.The optimal duration of the second stage -the period of backward plume transport -is determined here empirically for the Hateruma site as 7 days.Results presented by Gloor et al. (2001) suggest that a shorter period may suffice for other sites.To test the sensitivity of the coupled model to the duration of the Lagrangian simulation, we calculated CO 2 concentrations at Hateruma for a period of a year using a double backward transport period of 14 days.Comparison of the 7-and 14-day simulations revealed no significant difference in peak shape or amplitude of CO 2 concentrations, although a shorter duration leads to a decrease in peak amplitudes.Hence, we employ 7 days as the backward transport period.Concentration differences were within 2 ppm, and the averaged absolute value of the differences was 0.52 ppm, which is a modest discrepancy comparable with the representation error in observations.
Simulated concentrations are predicted in four steps, as follows.(1) FLEXPART simulates 7-day back-plume transport, which represents the time reversed trajectory of 10 000 air particles released from the location of interest.The actual transport time varies between 6.875 and 7 days, depending on the time of particle release, because coupling with the global Eulerian model occurs at 3-h intervals and the release of 10 000 particles takes 3 h.(2) The obtained particle distributions are combined with surface CO 2 fluxes at the location of each particle in order to calculate the contribution of the surface fluxes, during the 7-day simulation period, to the observed concentration ( C). (3) The CO 2 concentrations simulated by the NIES TM at the location of each particle, at the moment of coupling, is averaged over the total number of particles to yield an initial CO 2 concentration for the Eulerian stage (C init ).This initial concentration includes the contribution of surface fluxes before the start of the Lagrangian stage.(4) The initial CO 2 concentration (C init ) and the change in concentration along the trajectories of the Lagrangian stage ( C), as contributed by surface fluxes, are then added together to give the predicted concentration (C cpld ) at the observation time and location.
Consequently, the CO 2 concentration predicted by the coupled model, C cpld , is expressed as C can be calculated based on trajectory statistics (Seibert and Frank, 2004), and is expressed as where M o is the total mass of air particles released from the monitoring site; T b and T e are the beginning and end times of the back-plume trajectory calculations, respectively; and m(r,t) is the particle density field, which can be expressed via individual particle positions, as follows: where N is the number of particles, µ i = M 0 N and is i-th particle mass, r = (x,y,z) is the position vector, r i is an i-th particle position at time t, and δ(r) is a Dirac delta function.
Assuming the surface flux, F , is injected homogeneously into a near-surface layer of thickness z s , we can approximate the rate of change of the near-surface CO 2 concentration contributed by a surface flux, δC/δt, as follows: Here, z s is assumed to be 300 m, which is approximate to the lowest model level of the NIES TM and is selected to balance the exact approximation of the transport equation.This is achieved with small z s but requires a very large number of particles for transport.Hence, z s is given the largest possible value without compromising the results for a scenario of a well-mixed boundary layer during the day (assuming that observations of stable tracers that are frequently sampled during the day will be analyzed).Therefore, the selected value of z s offers a rough approximation of the minimum value of the daytime mixed-layer height.From Eqs. ( 2), (3), and (4), C can be calculated as follows: A time interval, dt, of 1 h is assigned for the purposes of this study.
C init is calculated as follows: where C 3-D is the CO 2 concentration in each simulation grid cell of the NIES TM.Therefore, using Eq. ( 3) for the particle density field, m(r,T e ), both C and C init (T e ) can be calculated without determining explicitly the global particle density field in the Lagrangian stage (FLEXPART), which is a convenient property when the surface flux, F , is given at a very high resolution.
Forward calculations can also be performed using an LPDM alone, by using Eq. ( 5) as a crude representation of C 3-D , as tested by Stohl et al. (2009).To reproduce seasonal variations, it is necessary to obtain transport data covering a period of at least 3 months (Stohl et al., 2009).However, computing a 10 000 particle plume transport over 3 months with FLEXPART for each 3-hourly observation is computationally expensive, with the result that such calculations are inefficient.In contrast, the coupled model requires only a short period of back-plume trajectory calculations to reproduce seasonal variations.
Here we use the coupled model to simulate variations in CO 2 concentrations at Barrow (USA), Samoa and Hateruma (Japan).To achieve this, the coupled model is used to simulate the transport of the biospheric CO 2 tracer, the oceanic CO 2 tracer, and the anthropogenic CO 2 tracer.Biospheric CO 2 is simulated using CASA model fluxes (Randerson et al., 1997), and the ocean CO 2 tracer is based on air-sea fluxes reported by Takahashi et al. (2002).These monthly fluxes are converted into daily fluxes by linear interpolation.For the anthropogenic CO 2 tracer, we employ the flux of fossil fuel emissions produced by the Emissions Database for Global Atmospheric Research (EDGAR), version 4 (EDGAR, 2009).Fossil fuel emissions are kept constant throughout the simulation.All fluxes used in this study have a horizontal resolution of 1.0 • × 1.0 • .FLEXPART is forced with 6-hourly analysis data provided by the Global Forecast System (GFS) model of NOAA's National Center for Environmental Prediction (NCEP).GFS data have a horizontal resolution of 1.0 • × 1.0 • and 26 pressure levels.The NIES TM is driven by 12-hourly re-analysis data provided by NCEP's National Centre for Atmospheric Research (NCAR) Re-analysis Project (Kalnay et al., 1996), and is run with a horizontal resolution of 2.5 • × 2.5 • and 15 sigma levels.In the simulation of the NIES TM, the initial concentration of all grids is set to zero and the first 2 yr of simulation are used as spin-up time.
To deseasonalize CO 2 concentrations, we employ a running average method, as follows: www.geosci-model-dev.net/4/317/2011/Geosci.Model Dev., 4, 317-324, 2011 where C i is the CO 2 concentration corresponding to the sampling date and time (i), n is the time window for averaging (set to 90 days), and Ci is the deseasonalized CO 2 concentration.

Results and discussion
Figure As shown in Fig. 2, the coupled model performs better than the NIES TM in simulating deseasonalized CO 2 concentrations, particularly at Hateruma.Figure 3 shows the absolute value of the difference between deseasonalized modelsimulated and observed CO 2 concentrations (| CO 2 |) at Hateruma.This figure illustrates how the coupled model is generally closer to observations than is the NIES TM, showing an improvement in the fit between models during shortterm, high-concentration events.The root mean square errors (RMSEs) of | CO 2 | for the NIES TM and the coupled model are 2.21 and 1.81, respectively.
These RMSEs are confirmed by the correlation coefficients and variance ratios between deseasonalized simulations and observations, calculated using annual data for the period 2002-2004 (Tables 1 and 2, respectively).The variance ratio is the standard deviation of the model-simulated CO 2 concentrations, normalized by the standard deviation of observed CO 2 concentrations.The use of correlation coefficients and variance ratios for model validation is a common method (e.g., the TransCom Project; Patra et al., 2008).It is possible to obtain better variance ratios and correlation coefficients if they are calculated using data for winter periods, as opposed to annual data.This is due to the fact that numerous peaks in observed CO 2 concentrations occur during the winter at Hateruma, as the region is influenced by contaminated air masses from the Asian continent and by the monsoon season.
Table 1 lists the correlation coefficients between simulated (NIES TM and coupled model) and observed CO 2 concentrations.Correlation coefficients determined from the Barrow data are similar for both the NIES TM and the coupled model.However, for Samoa and Hateruma, the coupled model is in better agreement with observations, as the cor- relation coefficients for the coupled model are >0.1 higher than those determined from the NIES TM for all years (Table 1).This improvement in correlation coefficients can be considered significant, given the typical variability in correlation coefficients for transport models (e.g., Patra et al., 2008), which range between 0.2 and 0.8, although the higher values in this range are rarely obtained.
To demonstrate the high reproducibility of observed CO 2 concentrations during winter periods in comparison with other times of the year, correlation coefficients are calculated using winter data as well as annual data (Table 1).As a result, the correlation coefficient between model and observations increases by between 0.05 and 0.1 (from annual values of 0.5-0.6). Figure 4 shows a scatter plot of the relationship between model-simulated and observed CO 2 concentrations at the synoptic scale for the period 2002-2004 at Hateruma.Linear regressions in the figure are based on the least squares fitting technique.Most of the data are within ±4 ppm of the linear regression, and the coupled model provides a better fit (and a higher correlation coefficient) than does the NIES TM (Table 1).
Table 2 lists the variance ratios between simulated and observed CO 2 concentrations, and highlights how variance ratios calculated for the coupled model are closer to unity than those determined for the NIES TM (with the exception of Samoa).This finding indicates that coupling with LPDMs helps to compensate for the smearing of concentration fields,  which arises as a result of the numerical diffusion intrinsic to Eulerian models.Variance ratios for winter periods only are calculated for Hateruma, and are more realistic than annual variance ratios or model-observation correlation coefficients.
Use of the coupled model is not, however, equally beneficial for the three sites, which were selected to represent different meteorological conditions.This is because the correlation coefficients between modeled and observed time-series vary spatially, as has been examined in previous studies using Eulerian models.For example, Patra et al. (2008) highlight difficulties in simulating the Barrow site.Remote locations such as Samoa tend to be characterized by lower correlation coefficients and variance in CO 2 concentrations because these regions are less affected by pollution emissions.The lack of observations at remote sites makes it difficult to correct for errors associated with, for example, wind speed and direction, and hence leads to a degradation in the quality of the simulation.Thus, simulations for remote/marine sites tend to be poor in comparison with simulations for sites situated closer to strong emission sources (e.g., Fig. 6d in Patra et al., 2008).Therefore, it is not surprising that the correlation coefficients and variance ratios determined in this study are not drastically improved by introducing a more accurate transport method, despite the use of wind fields of similar quality.

Conclusions
We To investigate the ability of the coupled model statistically, the model-observation correlation coefficients and variance ratios were calculated for three sites over 3 yr: 2002-2004.In the case of Hateruma, the correlation coefficient of the coupled model exceeds that of the NIES TM by more than 0.1 in all years.The advantage of the coupled model is seen clearly when reproducing pollutant events observed during winter at Hateruma, showing how the coupled model offers improved simulations for sites influenced by contaminated air.
There are, however, some instances where the advantage of the coupled model over Eulerian models is less certain.For example, the correlation coefficients of the coupled model and the NIES TM are similar for the Barrow site, and the variance ratios of the NIES TM are closer to observations than those of the coupled model at Barrow and Samoa.However, the reproducibility of observations is not only dependent on the transport model employed, but also on the quality of surface flux data.Another important issue is the quality of data for wind and vertical mixing parameters, which are used for the forward simulations in both the Eulerian and Lagrangian stages of the coupled model.The misfit between modeled and observed data can be ascribed to numerical diffusion in the Eulerian model, whereas the major issue for the Lagrangian model is the imperfection of the wind field, which is a primary limiting factor in terms of the accuracy of forward and inverse modeling.
Medium-resolution Eulerian models can only resolve concentration variations at close to daily time-scales.Thus, considerable useful information at shorter time-scales is lost, including tracer correlations.In contrast, the coupled model is able to resolve concentration variations at hourly time-scales or less.Consequently, the coupled model performs better than a global Eulerian model alone for simulating synopticscale variability.Hence, the coupled model can be used for analyses variations in CO 2 that occur as a result of synoptic-scale meteorological phenomena, without the need for nested regional-scale modeling.
The use of a coupled model also enables the calculation of surface flux contributions to tracer concentrations at a resolution higher than that of the meteorological data used in an LPDM.This is because the spatial scale of heterogeneity in surface fluxes is often much smaller than the correlation radius of meteorological data.Accordingly, surface flux data with a very high resolution (e.g., several kilometers) can be used for forward transport modeling without compromising the effective resolution of tracer transport.In the case of mesoscale circulations, it is necessary to obtain sufficiently fine resolution wind data; thus, application of large-scale interpolated winds may not be appropriate.

Figure 3 .Fig. 3 .
Figure 3. Absolute value of the difference between deseasonalized model-simulated and observed CO2 concentrations at Hateruma.Blue lines: NIES TM minus observations; red lines: coupled model minus observations.Root mean square error (RMSE) is shown in the figure.Blue: NIES TM; red: coupled model.

Figure 4 .Fig. 4 .
Figure 4. Scatter plot showing the relationship between model-simulated and observed CO2 concentrations at a synoptic scale for the period 2002-2004 at Hateruma.Blue circles: NIES TM versus observations; red circles: coupled model versus observations.

Table 1 .
Correlation coefficients between simulated and observed deseasonalized CO 2 concentrations (Winter: from the beginning of January to the end of March, and from the beginning of November to the end of December).

Table 2 .
Variance ratios between simulated and observed deseasonalized CO 2 concentrations (Winter: from the beginning of January to the end of March, and from the beginning of November to the end of December).