Greenhouse gas network design using backward Lagrangian particle dispersion modelling – Part 1: Methodology and Australian test case

. This paper describes the generation of optimal atmospheric measurement networks for determining carbon dioxide ﬂuxes over Australia using inverse methods. A Lagrangian particle dispersion model is used in reverse mode together with a Bayesian inverse modelling framework to calculate the relationship between weekly surface ﬂuxes, comprising contributions from the biosphere and fossil fuel combustion, and hourly concentration observations for the Australian continent. Meteorological driving ﬁelds are provided by the regional version of the Australian Community Climate and Earth System Simulator (ACCESS) at 12 km resolution at an hourly timescale. Prior uncertainties are derived on a weekly timescale for biosphere ﬂuxes and fossil fuel emissions from high-resolution model runs using the Community Atmosphere Biosphere Land Exchange (CA-BLE) model and the Fossil Fuel Data Assimilation System (FFDAS) respectively. The inﬂuence from outside the modelled domain is investigated, but proves to be negligible for the network design. Existing ground-based measurement stations in Australia are assessed in terms of their ability to constrain local ﬂux estimates from the land. We ﬁnd that the six stations that are currently operational are already able to reduce the uncertainties on surface ﬂux estimates by about 30 %. A candidate list of 59 stations is generated based on logistic constraints and an incremental optimisation scheme is used to extend the network of existing stations. In order to achieve an uncertainty reduction of about 50 %, we need to double the number of measurement stations in Australia. Assuming equal data uncertainties for all sites, new stations would be mainly located in the northern and eastern part of the continent.


Introduction
Inverse modelling has been used extensively over the last 2 decades in carbon cycle research to estimate surface fluxes of carbon dioxide (CO 2 ) on multiple temporal and spatial scales (i.e. Enting and Mansbridge, 1989;Rayner et al., 1999;Rödenbeck et al., 2003;Chevallier et al., 2010), employing mainly flask and in situ data. These observations include measurements from surface stations, tall towers, air planes and ships. More recently total column data (i.e. from the Total Carbon Column Observing Network (TCCON), Wunsch et al., 2011) have also been included in inversion studies (Chevallier et al., 2011). The main focus in most studies has been on deriving CO 2 fluxes from atmospheric CO 2 concentration observations through the inversion of an atmospheric transport model at the global scale for large land regions or a coarse grid (Peylin et al., 2013).

9364
T. Ziehn et al.: Greenhouse gas network design over 100 stations using mainly flask samples. Flask measurements can be made with high accuracy and precision, and the GLOBALVIEW product states uncertainties of 0.5 to 1 ppm depending on the station's location. However, flask samples are usually only provided weekly or fortnightly, which results in a poor temporal resolution and sampling selected for background conditions. This is partly compensated for by continuous in situ measurements, which are becoming increasingly available from a number of stations worldwide. Nevertheless, the sampling network is still too sparse with many gaps (i.e. in the tropics) to derive CO 2 sources and sinks at a local scale due to the under-determined nature of the inverse problem (i.e. number of sites is smaller than the number of grid cells) (Kaminski et al., 1999). Many inversion studies therefore focus on the estimation of fluxes for large regions of the continental or subcontinental scale. A unique solution of the inverse problem can be obtained by including prior information on the CO 2 surface fluxes, which can be derived for example from high-resolution model simulations that include the terrestrial biosphere and ocean fluxes (Kaminski et al., 1999).
Another issue that arises from the existing network of sampling stations is that they are mainly located at remote sites away from strong sources and sinks so that they can sample clean (well-mixed) air. Key stations such as Mauna Loa, Hawaii (Keeling et al., 1976), or Cape Grim, Australia, provide valuable long-term time series of atmospheric CO 2 concentrations, which are crucial for monitoring global atmospheric trends (Francey et al., 2013), but they are not ideally placed to detect local changes. In fact, their focus is sampling under baseline conditions, so that any influence from local land sources is minimised.
In order to derive reliable estimates of CO 2 sources and sinks at a local scale, the existing network of CO 2 measurement stations needs to be extended. Network design studies focus on optimal extensions of the existing network by considering potential stations where no data are available yet. The network design is usually performed in two steps: (1) running an atmospheric transport model for a given network in backward mode to evaluate the cost function for the surface flux inversion and (2) running an optimisation algorithm to minimise the cost function for the network design.
One of the first network design studies for CO 2 was performed by Rayner et al. (1996), where they used Bayesian synthesis inversion and simulated annealing to optimise the location of atmospheric CO 2 and δ 13 C measurements to constrain the global carbon budget. The network was optimised for the uncertainty variance in global ocean uptake using the GISS tracer transport model at a very coarse resolution (24×36 grid points). Rayner et al. (1996) also added one station at a time successively placed at every model grid point to identify the global minimum for one extra station. This is known as incremental optimisation, described in detail by Patra and Maksyutov (2002), who preferred this method over simulated annealing for its computational efficiency and to provide a continuous evolution of the observation network. Patra and Maksyutov (2002) demonstrated that both methods perform equally well using a semi-Lagrangian model at the global scale with the resolution set to 2.5 • ×2.5 • . Flux uncertainties were calculated for 11 ocean and 11 land regions, using a base network of 115 stations from the GLOBALVIEW data set and a list of 446 pre-selected potential stations. Incremental optimisation was also used by Law et al. (2004) to identify where sites are best located to minimise the uncertainties on annual mean flux estimates for 12 subregions of Australia. The inversions were performed using response functions for 116 regions globally, with the focus on Australia, represented by 44 grid points which were treated as potential new locations for the network extension.
In this study we aim to improve the methodology used in Law et al. (2004) in order to assess how the existing network of CO 2 -observing stations in Australia can be extended to minimise uncertainties in CO 2 flux estimates for Australia. The novelty of our approach is to derive the atmospheric transport matrix, which is required to relate surface fluxes to concentration observations, from a Lagrangian particle dispersion model (LPDM) run in backward mode. Running the model in backward mode is more efficient in the network design case because the number of sources exceeds the number of receptors by far. In addition to the transport matrix, we require only the error statistics of the data but not their actual values, in order to calculate the optimal network. This allows us to extend or create a network of stations where no data are available. However, as with all network design methods based on inversion modelling, our approach is dependent on specific choices made in the set-up of the estimation problem such as the resolution at which fluxes are estimated or how the error statistics are represented. The error statistics are usually provided in the form of a covariance matrix, which is difficult to obtain, but it has a large impact on the network design (Rayner, 2004).
When applying our new approach to Australia as a test case, we use a list of candidate stations instead of evaluating optimal locations on a regular grid. This is more efficient than treating every grid point as a potential location. It also has the advantage that we can easily take existing infrastructure into account, which will consequently result in a more realistic and cost-effective network extension. We explore the regular grid approach in the companion paper . Although we aim to design a cost-effective network by preselecting potential stations, we do not intend to perform a comprehensive economic evaluation of those stations. This would require specific information with regard to actual costs in setting up measurement equipment and in maintaining a site. Costs may also differ greatly between different sites, and a thorough cost analysis would be required, which is beyond the scope of this paper. However, the approach that we introduce here for the network design is generic, allowing for the optimisation of a number of properties of the network (including cost efficiency). It is also possible to implement the network design in two stages: (1) perform a general search based on a regular grid and (2) perform a specific search accounting for the costs associated with setting up new sites.
In contrast to many previous studies that mainly used flask measurements from GLOBALVIEW, we consider continuous measurements at an hourly timescale for existing and potential stations. This allows us to derive CO 2 fluxes at a high spatial and temporal resolution.
This paper (Part 1) develops the generic framework for the network design and introduces the Lagrangian particle dispersion model which we run in backward mode to obtain the source-receptor (s-r) relationship. We then apply this concept to the Australian continent as a test case. In a first step, we evaluate the existing network of CO 2 ground-based measurement stations in terms of its ability to provide reliable flux estimates. In a second step, we demonstrate how the existing network in Australia can be extended in an optimal way.
A companion paper (Part 2; Nickless et al., 2014) focuses mainly on sensitivity analysis of parameters and choices necessary for running the optimal network design and their consequences on the results. This will be demonstrated for a South African test case, where the optimal network is created on a regular grid using continuous measurements from five new instruments.

Methodology
The network design is based on a combination of Bayesian inverse modelling methodology applied to an atmospheric transport model and the optimisation of a cost function. For a given network the cost function is calculated from the posterior statistics of an inversion to infer CO 2 surface fluxes from CO 2 concentration measurements. The optimisation will then find the optimal network by minimising the cost function through altering the given network by adding or removing new stations.
The observed concentration (c) at a particular station at a particular time can be expressed as the sum of different contributions: where c f is the contribution due to surface fluxes within the modelled domain, c b the contribution from outside the region of interest (boundary inflow) and c i the contribution from the initial conditions. For the network design the initial conditions are neglected because they are very well constrained by the observations and their contribution to the flux uncertainty is therefore thought to be small. The contribution from the boundaries has to be assessed, and if the influence on the flux uncertainties is not negligible, then the boundary conditions have to be included in the network design process as well.

Surface flux inversion
We use a Bayesian synthesis inversion scheme (Tarantola, 1987;Enting, 2002) which allows us to infer CO 2 surface fluxes from CO 2 measurements. A simple linear expression can be used to model the relationship between the surface fluxes and concentrations: where c mod is the vector of the modelled concentrations and f the vector of the (unknown) surface fluxes. T is the transport or sensitivity matrix which needs to be determined. At this stage we do not include any influence from outside the domain and therefore assume that the modelled concentrations (c mod ) are equal to the concentrations (c mod f ) which are derived from the surface fluxes only. If we assume a Gaussian error distribution for the surface fluxes and concentrations, we can obtain the maximum likelihood estimate for f by minimising the cost function: where C c is the error covariance matrix of the observations, vector f 0 contains prior flux estimates, vector f represents predicted fluxes and C f 0 is the prior error covariance matrix of the surface fluxes. The cost function therefore ensures that we simultaneously minimise the mismatch between modelled concentrations and measurements and the mismatch between prior flux estimates and predicted fluxes. The solution of the optimisation problem expressed through the cost function in Eq. (3) provides optimal surface fluxes based on the observations provided and also posterior uncertainties for the CO 2 fluxes expressed through the posterior covariance matrix C f . For the network design approach we are only interested in the latter, because our aim is to find a network (set of observations) that minimises the CO 2 flux uncertainties. The posterior covariance matrix can be calculated by either of the two equivalent expressions (Tarantola, 1987): As noted by Hardt and Scherbaum (1994), the calculation of the posterior flux uncertainties does not depend on a particular value of the surface fluxes or concentration observations. It only depends on the transport model, the prior flux uncertainties and observational uncertainties. This has the advantage that we can evaluate potential stations for which we do not have real observations yet and without the need to generate synthetic data.

Lagrangian particle dispersion model (LPDM)
The relationship between surface fluxes and atmospheric concentrations is embodied in the transport matrix T. A common approach in deriving T is to use a Lagrangian stochastic particle dispersion model. The conventional approach is to run the model in forward mode, where particles are released at the surface (source) and tracked until they have passed the measurement station (receptor), which means that all particles need to be tracked even if they do not pass through the receptor. However, if the number of sources exceeds the number of receptors, then it is more efficient to run the Lagrangian model in reverse or backward mode, where the particles are released at the receptor and tracked backwards in time to any potential surface source. The source-receptor relationship can then be used to derive the transport matrix T.
Here we use a LPDM (Uliasz, 1994), which we run in reverse mode for each potential and existing measurement station we would like to include in the network design process. Particles are released (from the known or proposed measurement height) every 20 s for a total of 4 weeks for different seasons of the year, and the particles position is recorded at 15 min intervals. Particles that are near the surface are counted for each grid cell to determine the surface influence or sensitivity. This can be used to generate a footprint for each station, which shows the area of influence, and also to calculate the s-r relationship, which forms the transport matrix T. Here, we follow Seibert and Frank (2004) to derive the elements of that matrix.
According to Seibert and Frank (2004) the s-r relationship for a point source (one grid cell source) is given as whereχ is a mass mixing ratio (receptor) andq in is a mass flux density (source). The abbreviation µ tot stands for the total initial mass released at the receptor in a time interval, c in is the mass concentration and ρ in is the air density. The index in indicates the ith grid element and the nth time interval of length T . The overbar indicates temporal averaging over the time interval T , and V i is the volume of grid element i. The LPDM output does not provide mass concentrations (i.e. c in ) for a grid cell, but the number of particles near the surface. However, the number of particles in a given grid cell is directly related to their contribution in mass. Therefore, we do not need to assign a mass to the particles, but instead we can express the mass of the particles as a function of the number of particles, and after cancelling terms the particle count can be used directly to express the source-receptor relationship in the following way: with N in the number of particles in a grid element (source) at each time interval T and N tot the total number of particles released during a time interval. Using the number of particles instead of their mass concentration in Eq. (6), we get Note that our s-r relationship now becomes independent of the grid cell volume V i . The density of air can be calculated as where P is the pressure difference in the surface layer and g is the gravity of Earth. We also apply a conversion from mass mixing ratio to volume mixing ratio. This is simply done by multiplying with the ratio of the molecular mass of air to the molecular mass of carbon, which is our quantity of interest. The elements of the matrix T are now calculated as withχ expressed in ppm. For the network design we are interested in weekly fluxes of carbon divided into day-and night-time contributions, which reflects the way a flux inversion is usually done. This means that we have to provide the particle count N in as the sum over 1 week ( T = 1 week for day and night, divided at 06:00 and 18:00 Australian Eastern Time). Therefore, the mass flux densityq in in Eq. (11) has units of kg C m −2 week −1 (day/night). We set the surface layer height to 50 m, which corresponds to approximately 600 Pa ( P ). If we consider well-mixed conditions, then the s-r relationship should be independent of the thickness of the surface layer as long as the layer is not too deep (Seibert and Frank, 2004). This is further investigated in the companion paper .

Influence from outside the modelled domain
The inflow from the boundary can affect the concentrations measured at a certain point. These so-called boundary effects can be included in our modelling approach in two different ways: (a) if they are significant, then we have to explicitly solve for them; (b) if they are small enough, we can treat them as contribution to noise. If we decide to solve for the boundary concentrations on top of all the surface fluxes, then the modelled concentrations are given as where c mod b is the modelled contribution from the boundaries.
The contribution from fluxes outside the modelled domain can be treated via their effect on boundary concentrations (c B ). In order to assess the influence of the boundary concentrations on the observed concentrations (c) we need to determine the strength of the connection between the two. This can be done by calculating the Jacobian, which provides the sensitivities of observed concentrations to boundary concentrations. The boundary contribution can then be written as where M B is the Jacobian. Depending on the elements of M B , we might need to include the boundary conditions in the network design. The elements of the Jacobian for the boundary conditions can be calculated by accounting for the number of particles that disappear from the model domain during the simulation. LPDM can be set up to write out the location and time when particles leave the domain, and one can decide on a spatial and temporal resolution (Lauvaux et al., 2012). Here, we consider four boundaries (north, south, east and west), and we calculate the sensitivity of hourly observed concentrations to weekly boundary concentrations. In this way the Jacobian (M B ) for each site has 32 columns (4 boundaries × 4 weeks × 2 (day/night)) and 672 rows (hourly observations over 4 weeks), with its elements calculated as where N B is the number of particles leaving the domain at one of the four boundaries during 1 week (day/night) and N tot is the total number of particles released during 1 h. Ideally, we need to calculate M B for each station and then use a criterion to assess whether or not the boundary conditions affect that station. Note that we are neglecting the influence of the top boundary of the domain on the observations. This is likely to be both small and homogeneous (hence indistinguishable from the initial condition). We can use the following simple test to assess the effect of the boundary concentrations on the network design: where C I is the identity matrix. The diagonal elements of C b provide us with the uncertainty contribution of the boundary concentrations to the uncertainty of the observations. If they are small compared to the assumed observational uncertainty, then the uncertainty contribution of the boundary concentrations can also be considered small and we do not need to include them explicitly in the network design process. This means that we could use Eq.
(2) again and treat the boundary effects as contribution to noise instead.
The reason for assessing the effects of the boundary conditions first instead of using Eq. (12) as a standard case is that, if we need to include them in the inversion, then we would have to solve for the boundary concentration (c B ) in addition to all the surface fluxes. This not only means that we need to combine the transport matrix (T) and the Jacobian (M B ) into a new expanded transport matrix which has a much larger dimension, but we would also need to provide prior estimates and uncertainties for the boundary concentrations. These are hard to assign. The optimal network should seek to reduce the uncertainty of the surface fluxes, and the improvement of the tracer transport in the global circulation models should be left as a separate problem. Since it is quite challenging to provide sensible estimates for the prior uncertainties of the boundary concentrations, we would like to include them only if required (i.e. if it changes the outcome of the network design).

Network design for Australia
For the network design we run LPDM in backward mode for each station that we would like to include in this study. We start by assessing the stations in the base network in terms of their ability to reduce the uncertainties on net CO 2 flux estimates. We then add new stations from the candidate list to the base network using an optimisation scheme. Finally, we compare this optimal network with a network that was designed from scratch (i.e. we assume no existing stations).

Ground-based measurement stations
Australia has nine established ground-based measurement stations (see Fig. 1a and Table 1) run by CSIRO or the University of Wollongong. We exclude Cape Ferguson and Otway because they currently provide only flask data, and we also exclude Tumbarumba because it is not operational at this time. The remaining six stations provide continuous CO 2 measurements and form our base network. From the location of the six stations ( Fig. 1a) it is obvious that Australia as a whole is not very well covered since the site locations were not determined with the goal of estimating Australian CO 2 fluxes. Rather the sites consist of (a) Global Atmosphere Watch (GAW) locations, focused on measuring baseline air; (b) TCCON locations; (c) locations of the institutions running the sites; and (d) locations linked to specific projects. In order to estimate CO 2 fluxes from the terrestrial biosphere, we require stations that are able to pick up the signal from local sources. For the existing network this will depend on the wind direction, and we will show later to what degree the base network is already able to reduce the uncertainties on Australian CO 2 flux estimates.
In order to improve the accuracy of CO 2 flux estimates for Australia, we need to add new stations to the base network. The optimal location of new stations is determined by minimising a cost function (see Sect. 2.4.3) which is calculated  Table 1 for existing sites and Table 2 for potential sites. Existing stations that are not included in the base network are marked in light blue. for a number of potential locations. There are several ways of setting up a list of potential stations or candidate stations. The simplest way is to assign the stations according to a regular grid. However, this might lead to a very large number of potential stations, of which many may be located in inaccessible areas.
To design a more realistic and cost-efficient network, we can include logistic constraints such as the availability of supporting infrastructure as a limiting factor in the selection of stations for the candidate list. For example, one could use the location of airports or wind farms in Australia. There is also a large number of telecommunication towers along main roads which could potentially be used. Here, we use the location of the Australian Bureau of Meteorology weather watch radar stations (NRL, 2014) as potential stations. This guarantees that all stations are accessible by road and maintained. The list of all 59 potential stations can be found in Table 2, with their location shown in Fig. 1b.

Driving data and prior uncertainties
LPDM requires meteorological driving fields, which are provided in this study by the regional version of the Australian Community Climate and Earth System Simulator (ACCESS-R) (NMOC, 2013) at 12 km resolution for the Australian region at an hourly timescale. Driving data include the 3-D wind field, temperature and turbulent kinetic energy (TKE) at 39 vertical levels up to 18 km in height as well as surface pressure. These fields are provided for one example month (4 weeks) for Southern Hemisphere (SH) winter (July) and summer (January).
We also need to derive prior surface flux uncertainties for Australia and an estimate of the observational uncertainties (i.e. accuracy of concentration measurements). In terms of the prior surface flux uncertainties we consider contributions from the biosphere and from fossil fuel combustion. We assume that sources change every week, and flux uncertainties are therefore calculated on a weekly basis.
The biosphere flux uncertainties (expressed as the standard deviation) are estimated using the following simple Atmos. Chem. Phys., 14, 9363-9378 relationship (Chevallier et al., 2010): where NEP is the net ecosystem productivity (net carbon flux) and NPP the net primary productivity. NPP is derived for the Australian continent from BIOS2 model simulations (Haverd et al., 2013) at a daily timescale ( Fig. 2a and c). BIOS2 is a modelling framework that uses the Community Atmosphere Biosphere Land Exchange (CABLE) model (Wang et al., 2010) at 5 km resolution (0.05 • × 0.05 • ). We then aggregate the high-resolution fluxes to the resolution that we use for the network design (1.8 • × 1.8 • ) and estimate the uncertainties for NEP according to Eq. (16) for each week divided into day-and night-time ( Fig. 2b and d).
Fossil fuel uncertainties are derived from the Fossil Fuel Data Assimilation System (FFDAS) (Rayner et al., 2010;Asefi-Najafabady et al., 2014). We use 10 realisations from FFDAS version II at 0.1 • × 0.1 • , aggregate them to our network design resolution and then calculate the uncertainties from the 10 realisations. Due to the fact that fossil fuel fluxes are derived on the basis of power plant locations and night lights, they are very localised and vary a lot in magnitude (Fig. 3a). When we aggregate those high-resolution fluxes to our 1.8 • × 1.8 • network design resolution, we "smooth out" most of the very large fluxes (Fig. 3c). Consequently, the variation between the 10 realisations of the aggregated fluxes also becomes smoother, which leads to only small uncertainties (Fig. 3d). As a result, fossil fuel flux uncertainties are much smaller (< 0.3 g C m −2 week −1 ) than the uncertainties from the biosphere fluxes, and their influence will also be small. Figure 3b shows the uncertainties for the 10 realisations based on the original 0.1 • × 0.1 • resolution, which are much larger for individual grid cells than the uncertainties calculated for the aggregated fluxes. If we performed inversions for only a small region of Australia using a much higher resolution, then the fossil fuel uncertainties would become much more important and, depending on the resolution, they might even dominate the overall surface flux uncertainties. However, due to computational limitations we decided not to increase the resolution for the network design in this study. The influence of the spatial surface flux resolution on the outcome of the optimal network design is investigated in Part 2 .
Finally, we estimate the prior error covariance matrix of the land surface fluxes as where vector l f contains the land fractions, vector b σ 2 the variance for the biosphere fluxes and vector u σ 2 the variance for the fossil fuel emissions for each grid cell and each week (separated into day-and night-time). The operator "diag" returns a diagonal matrix with the vector elements as the diagonal, which means that we assume no correlations among different fluxes. The effect of correlation length between different fluxes is investigated in Part 2 . Multiplying by the land fractions guarantees that the prior uncertainties for coastal grid cells are scaled accordingly and ocean-only grid cells are set to 0. This is important because in the network design we want to focus on the reduction of uncertainty for the land fluxes only. Ocean flux uncertainties are not considered in this study because they are usually by a factor of 10 per unit area smaller than the land flux uncertainties (Chevallier, 2007). Due to the small amount of ocean grid cells in our modelled domain, the impact of the ocean flux uncertainties is expected to be small. Nevertheless, the contribution of the ocean fluxes to the posterior covariance matrix and the optimal location of stations is investigated in Part 2  for a South African test case.
Observational uncertainties are set to 2 ppm for all existing and potential stations (except in one sensitivity test). Again, the uncertainties are specified in terms of their standard deviation, and we assume no correlations among the uncertainties of different observations. In this way C c also becomes a diagonal matrix.

Cost function for the network design
We must optimise some scalar quantity derived from the posterior covariance. Rayner et al. (1996) noted the sensitivity of the optimal network to this choice. Common options are the average uncertainty of individual fluxes (the trace of the covariance, cost function J Ct ) or the uncertainty of the integrated flux (the sum of all elements, cost function J Ce ): where n is the number of elements in the diagonal of the matrix C f . In this study we use Eq. (19) because our focus is on the uncertainty reduction of the total flux estimate. The impact on the optimal network by using one cost function over the other is investigated in Part 2 . If we start the optimisation from the base network, then the transport matrix always includes the s-r relationship for our six stations in the base network (see Table 1). We can then add the s-r relationship for the three remaining existing stations and/or for the stations from the candidate list and construct C f . In order to find the set of stations that minimises our cost function, we apply the incremental optimisation, where we add only one station at a time from the candidate list to the base network and calculate C f . We choose the station that gives us the smallest cost function value, add it to the network and also remove it from the candidate list. We than repeat the process until our optimal network has reached a certain maximum size or the candidate list is empty.
We assume that observations (CO 2 concentration measurements) will be available from all stations from the candidate list and the base network at an hourly timescale. The s-r relationship we calculate with LPDM therefore represents the sensitivity of hourly observations to weekly fluxes.
We evaluate the different networks in terms of the uncertainty reduction defined as whereĴ Ce is the optimal cost function value and J Ce prior the cost function value based on the prior uncertainties. Instead of J Ce prior we could also use J Ce base , which is the cost function value for the base network. The footprint is the sum over the influence functions for 1 month and shows the number of particles that are in touch with the surface.

Results and discussion
After running LPDM for all existing and potential stations, we calculate the influence function or sensitivity matrix for each of the stations. We can also sum over the influence functions for the whole month, and this provides us with the surface footprint for each station. Figure 4a and b present the footprint for Cape Grim as an example of an existing station from the base network for July and January respectively. It can be seen that the area that is observed by Cape Grim differs by a large amount between the two seasons. In SH summer, Cape Grim samples mainly clean air coming from the Southern Ocean. The influence from the land is very small. However, in SH winter the dominant wind direction varies, and Cape Grim is sampling air that may also be influenced by surface fluxes from the south-eastern part of Australia. The surface footprint for a potential station in Alice Springs is presented in Fig. 4c and d for both seasons. Due to its central location, a station in Alice Springs would be able to detect the influence of potential surface fluxes from a large part of the Australian continent. However, from the surface footprint alone we cannot estimate how much a station at Alice Springs would help us to reduce the uncertainties on net CO 2 fluxes.
We then use Eq. (15) to decide whether or not we have to include the boundary conditions in our inversions explicitly or if we can treat them as contribution to noise. In order to do this, we investigate existing and potential stations close to the north, south, east and west coast of Australia (i.e. Darwin, Aspendale, Arcturus and Geraldton respectively). All diagonal elements of C b turn out to be small for those stations, which means that the uncertainty contribution of the boundary concentrations to the uncertainty of the observations can be considered negligible. Therefore, we decided not to include the boundary concentrations in the network design process. Note that this would change for a smaller domain or one where the large-scale concentrations were more uncertain.

Base network
Currently, there are six ground-based measurement stations in Australia that measure CO 2 continuously. As discussed earlier, some of these stations were designed to measure well-mixed air (i.e. Cape Grim) or background concentrations for detecting fugitive emissions (i.e. Arcturus). However, the surface footprint of Cape Grim, for example, indicates that our existing stations are also able to pick up the influence from the land depending on the dominant wind direction. Here, we test how useful our existing stations are Table 3. Ranking and uncertainty reduction (U R ) for the existing stations in the base network in terms of their ability to reduce the uncertainties on CO 2 flux estimates for two seasons (SH summer and winter), represented by January and July. The data uncertainty for all stations is set to 2 ppm. The station number is provided in brackets.

Rank Station July
U R Station January U R Station July + January U R  Table 4. Ranking and uncertainty reduction (U R ) for the existing stations in the base network in terms of their ability to reduce the uncertainties on CO 2 flux estimates for two seasons (SH summer and winter), represented by January and July. The data uncertainty is set to 1 ppm for Cape Grim, to 3 ppm for Aspendale and Wollongong and to 2 ppm for all remaining stations. The station number is provided in brackets.
Rank Station July U R Station January U R Station July + January U R in terms of estimating CO 2 fluxes from CO 2 concentration measurements. In a first experiment, we assume that all stations provide the same quality of measurements, and we set the data uncertainty to 2 ppm for each station. In a second experiment, we assign a lower uncertainty to measurements from Cape Grim and a higher uncertainty to measurements from Aspendale and Wollongong. Table 3 shows the ranking and uncertainty reduction for the first experiment for all stations in the base network for the two seasons individually and together. We use J Ce as a cost function, which means that we include all elements of the posterior covariance matrix. Incremental optimisation is used to determine the ranking of the stations and their overall contribution to the uncertainty reduction. We start with an empty network and then add the station which provides the greatest reduction in uncertainty. We repeat this until all six stations have been added to the network, which then forms the base network.
The results vary for the two seasons. We get a larger reduction of uncertainty in SH winter (July) than in SH summer (January) due to the difference in dominant wind direction and due to the fact that the prior biosphere flux uncertainties are also larger in July (see Fig. 2). However, for both seasons individually and together, Darwin, Wollongong, Arcturus and Aspendale rank as the four most important stations in our base network. These stations are already able to reduce the uncertainties on CO 2 flux estimates by more than 27 %. Gunn Point and Cape Grim are the least important stations. However, this is misleading if one is only assessing the ranking as presented in Table 3. In fact, Gunn Point alone provides about the same reduction in uncertainty as Darwin (12.47 % vs. 12.81 %) because these two stations are located very close together. Due to the fact that the uncertainty reduction for Darwin is slightly larger than the one obtained from Gunn Point, Darwin is added first to the network which makes Gunn Point "redundant". Cape Grim on the other hand provides the smallest reduction in uncertainty even when assessed on its own in an empty network, because Cape Grim samples the "cleanest" air of all the six stations.
The ranking of the existing stations in the base network also depends on the observational error assigned to those stations. As has been noted elsewhere (e.g. Rayner et al., 2010), this uncertainty includes not only the error in the actual measurement but the difficulty in simulating it within the model used in the inversion. In a sensitivity experiment, we set the the observational uncertainty for Cape Grim to 1 ppm, because this is Australia's primary ground-based measurement station and we expect a high accuracy for the data. We increase the observational uncertainty for Aspendale and Wollongong to 3 ppm, because these two stations are located close to large sources of fossil fuel emissions. We keep the uncertainty at 2 ppm for all remaining stations. The new ranking of the existing stations can be found in Table 4. It can be seen that Cape Grim now becomes one of the most important stations in the base network. In contrast, Aspendale and Wollongong, which were ranked high in the first experiment, become less important. This highlights the sensitivity of the network design to the observational uncertainty assigned to Table 5. Ranking and uncertainty reduction (U R ) for the new stations added to the base network for two seasons (SH summer and winter), represented by January and July. The data uncertainty for all stations is set to 2 ppm. The station number is provided in brackets.

Rank Station July
U R Station January U R Station July + January U R The small uncertainty reduction that we achieve in SH summer suggests that the current network is not suitable for estimating biosphere fluxes for Australia for that season. Overall, the six existing ground-based measurement stations are able to reduce the uncertainties on CO 2 flux estimates for Australia by nearly 30 % in both experiments for the two seasons together. This is an interesting result since most stations were not primarily designed to measure the contributions from land fluxes.

Extended network
We extend the base network by one station at a time using the incremental optimisation for each season individually and for both seasons together. We add a total of six new stations from the candidate list, and the results are presented in Table 5. The results show different network extensions for the two seasons, with only one station (Moree) in common. The base network already provides a substantial reduction in flux uncertainties for the SH winter season (43 %), and the six new stations allow for a further 18 %. Stations are mainly added in the northern part of the Australian continent (see Fig. 5a). In the SH summer season, the base network can only provide an uncertainty reduction of 10 %, but with the new stations added we can achieve an additional 30 %. The six new stations are mainly added in the north-eastern part of Australia (see Fig. 5b), filling the gaps between existing stations.
If we focus on the results for the network considering both seasons together, the first four stations added to the base network are the same as for the SH summer case, albeit in a slightly different order. This is in agreement with the fact that new stations are able to provide a greater reduction in uncertainty for the SH summer than for the SH winter. The two additional stations (Longreach and Tennant Creek) in the extended network for both seasons are located in the centraleastern part of the country (see Fig. 5c).
Adding six new stations to the base network would lead to a doubling in the number of ground-based measurement stations in Australia, and we would be able to achieve a reduction on the prior uncertainties of CO 2 flux estimates of more than 50 %. It is worth pointing out that Tumbarumba (ranked third in the optimisation; see Table 5) is actually not a new station from the candidate list. Tumbarumba is an established station that is currently not operational (see Table 1). Tumbarumba used to measure CO 2 continuously, and the network design indicates that it has great value for estimating local land fluxes.

New network
Another interesting scenario is to perform the network design by starting from an empty network (i.e. assume we do not have existing ground-based measurement stations). In this way we will be able to create the most efficient network and be able to compare it with the extended network from the previous section that is based on our existing stations.
The ranking of the stations added to the network by the incremental optimisation is shown in Table 6. Again, the networks vary depending on the season due to the difference in wind direction and prior biosphere flux uncertainties, but also show some similarities. Mornington Island and Moree, for example, are the highest-ranked stations for each season individually and combined. In the SH summer season Tumbarumba turns out to be even more important. Stations at Mornington Island and Moree would be able to reduce the CO 2 flux uncertainties by more than 27 % (summer and winter). In comparison with the base network, we can see that these two stations alone would be able to provide the same reduction in uncertainty as the four highest-ranked existing stations (Darwin, Wollongong, Aspendale and Arcturus) together. If we wanted to achieve the same performance as the extended network consisting of 12 stations (6 existing stations plus 6 new stations), we would only require 9 stations when designing the network from scratch. The distribution of the stations in the new network (which is not based on the existing stations) does not look that much different from the extended network (see Fig. 5d-f); however the optimisation has more freedom to place stations closer to regions where prior uncertainties are largest.

Comparison with previous study
In a previous study, Law et al. (2004) found that, if we wanted to erect a new measurement station in Australia, it should ideally be located in the north-west or central part of the continent. This is in contrast to our results, which suggest that it would be most beneficial to add a new station in the north-eastern part of Australia. Both studies use a similar metric for the optimal network design, namely the reduction in flux uncertainty over the Australian continent, but there are also many differences in the methodology and set-up between the two studies. For example, Law et al. (2004) use response functions for the inversion and divide Australia into only 12 subregions for which average annual mean uncertainties are calculated. The base network in Law et al. (2004) consists of only two stations (Cape Grim and Cape Ferguson). In this study the base network comprises six stations, with two stations located in the north of the continent. However, even if we start the network design from the same base network as in Law et al. (2004), the first station added to the network will still be located in the north-east. The largest impact on the difference in the new stations' locations might be due to the formulation of the prior flux uncertainties. In Law et al. (2004) the prior flux uncertainties are either assumed to be constant (i.e. set to 1 kg C m −2 yr −1 ) for all regions or Table 6. Ranking and uncertainty reduction (U R ) for the new stations starting from an empty network for two seasons (SH summer and winter), represented by January and July. The data uncertainty for all stations is set to 2 ppm. The station number is provided in brackets.

Rank Station July
U R Station January U R Station July + January U R In both cases, the largest prior flux uncertainties are assigned to the north-west or central part of Australia. It is therefore not surprising that a new station would then also be located in the same region. In our study, the prior flux uncertainties are scaled with NPP. The largest uncertainties can be found in the productive north-eastern part of the continent, and this is where we would add the first new station. Again, this highlights that the location chosen for new stations (or network) critically depends on the prior knowledge (i.e. prior flux uncertainties) provided.

Logistic constraints for potential stations
Potential stations in this study were selected based on existing locations of stations in the Australian Bureau of Meteorology National Radar Loop (NRL, 2014). Although this ensures that all potential sites are accessible and maintained, we cannot differentiate between the sites in terms of actual costs associated in setting up the equipment to measure CO 2 concentrations. For example, some sites may require the erection of a tall tower, whereas other sites may already have one that could potentially be used for additional equipment. Costs for maintaining a certain site may also differ by a large amount due to the site's distance away from the nearest major town or airport or its location being offshore. It is very challenging to include all this information, which might not even be available at the time, into the cost function. Weights need to be assigned for all penalty terms, and that may require tuning, which is time consuming. One way to circumvent this problem would be to change a station's observational uncertainty as a proxy for logistical issues. For example, we could use a smaller observational uncertainty for a station that is easily accessible and cost efficient to run, and a higher observational uncertainty for a station that is located offshore and does not provide a tall tower already. In this way the cost function does not need to be changed and we can use all available information for a potential site by changing only one quantity. However, if information with regards to costs in setting up and maintaining a site is available, then this should be included in the cost function so that one can account for the exact economic costs. This will be investigated in a future study.

Summary and conclusions
Running a Lagrangian particle dispersion model in reverse mode provides an efficient way of obtaining the relationship between concentration observations (receptor) and ground fluxes (source). Here, we used LPDM and the Bayesian framework to obtain the transport matrix (source-receptor relationship) for existing ground-based measurement stations and a list of candidate stations. An incremental optimisation scheme was then used to design an optimal network of ground-based measurement stations for Australia. Existing stations were assessed and ranked in terms of their ability to reduce CO 2 flux uncertainties for the whole continent. New observational networks were designed based on existing stations and starting from an empty network. We found that the influence from outside the domain (boundary concentrations) has only a small impact on the network design, and we therefore did not include uncertainties related to boundary concentrations in this study. In addition to biosphere flux uncertainties, which were derived from high-resolution BIOS2 model runs, we also considered uncertainties for fossil fuel emissions. These uncertainties were derived from 10 realisations from FFDAS at 0.1 • resolution. However, the fossil fuel fluxes are very localised with a range of many orders of magnitude. When we aggregated those fluxes to the 1.8 • resolution used for the network design, we smoothed out most of the large fluxes and their variability across the 10 realisations.
The assessment of existing ground-based measurement stations in Australia showed that they would be able to reduce surface flux uncertainties by about 30 %, which indicates the value of making in situ measurements (taken from all wind directions) at sites that are designed primarily for baseline measurements, such as Cape Grim.
If we want to halve the uncertainties on Australian flux estimates, we need to double the number of existing stations that are currently operational. Assuming all sites give measurements of the same quality and can be modelled equally well, the two most important new stations would be located in the north (Mornington Island) and in the east (Moree) of the continent. This also shows that new stations do not necessarily need to be located far inland in order to pick up the influence from local sources. In fact, a new station at Mornington Island would be located offshore, but still be able to observe the influence from the Australian biosphere. It is also worth noting that Tumbarumba, an existing station that is currently not operational, has great value for the estimation of local sources. In terms of costs associated with erecting a new site, it might be more efficient to put an already existing site back into operation.
Although we included logistic constraints in setting up the candidate list of stations, we did not include actual costs (i.e. to erect the station or for maintenance) in the network design. A relative measure of the ease of taking and modelling measurements at any given location can be accounted for in the network design through using a variable data uncertainty across measurement locations. This would be a valuable extension to this study, given the sensitivity found in the one test case undertaken in which data uncertainties were modified. We also plan to extend the study to optimise a network for estimating both CO 2 and CH 4 fluxes since some in situ instruments are designed to measure both species.