Analysis of total column CO 2 and CH 4 measurements in Berlin with WRF-GHG

Though they cover less than 3 % of the global land area, urban areas are responsible for over 70 % of the global greenhouse gas (GHG) emissions and contain 55 % of the global population. A quantitative tracking of GHG emissions in urban areas is therefore of great importance, with the aim of accurately assessing the amount of emissions and identifying the emission sources. The Weather Research and Forecasting model (WRF) coupled with GHG modules (WRFGHG) developed for mesoscale atmospheric GHG transport can predict column-averaged abundances of CO2 and CH4 (XCO2 and XCH4). In this study, we use WRF-GHG to model the Berlin area at a high spatial resolution of 1 km. The simulated wind and concentration fields were compared with the measurements from a campaign performed around Berlin in 2014 (Hase et al., 2015). The measured and simulated wind fields mostly demonstrate good agreement. The simulated XCO2 shows quite similar trends with the measurement but with approximately 1 ppm bias, while a bias in the simulated XCH4 of around 2.7 % is found. The bias could potentially be the result of relatively high background concentrations, the errors at the tropopause height, etc. We find that an analysis using differential column methodology (DCM) works well for the XCH4 comparison, as corresponding background biases are then canceled out. From the tracer analysis, we find that the enhancement of XCH4 is highly dependent on human activities. The XCO2 enhancement in the vicinity of Berlin is dominated by anthropogenic behavior rather than biogenic activities. We conclude that DCM is an effective method for comparing models to observations independently of biases caused, e.g., by initial conditions. It allows us to use our high-resolution WRF-GHG model to detect and understand major sources of GHG emissions in urban areas.

Abstract. Though they cover less than 3 % of the global land area, urban areas are responsible for over 70 % of the global greenhouse gas (GHG) emissions and contain 55 % of the global population. A quantitative tracking of GHG emissions in urban areas is therefore of great importance, with the aim of accurately assessing the amount of emissions and identifying the emission sources. The Weather Research and Forecasting model (WRF) coupled with GHG modules (WRF-GHG) developed for mesoscale atmospheric GHG transport can predict column-averaged abundances of CO 2 and CH 4 (XCO 2 and XCH 4 ). In this study, we use WRF-GHG to model the Berlin area at a high spatial resolution of 1 km. The simulated wind and concentration fields were compared with the measurements from a campaign performed around Berlin in 2014 . The measured and simulated wind fields mostly demonstrate good agreement. The simulated XCO 2 shows quite similar trends with the measurement but with approximately 1 ppm bias, while a bias in the simulated XCH 4 of around 2.7 % is found. The bias could potentially be the result of relatively high background concentrations, the errors at the tropopause height, etc. We find that an analysis using differential column methodology (DCM) works well for the XCH 4 comparison, as corresponding background biases are then canceled out. From the tracer analysis, we find that the enhancement of XCH 4 is highly dependent on human activities. The XCO 2 enhancement in the vicinity of Berlin is dominated by anthropogenic behavior rather than biogenic activities. We conclude that DCM is an effective method for comparing models to observations independently of biases caused, e.g., by initial conditions. It allows us to use our high-resolution WRF-GHG model to detect and understand major sources of GHG emissions in urban areas.

Introduction
The share of greenhouse gas (GHG) emissions released from urban areas has continued to increase as a result of urbanization (IEA, 2008;Kennedy et al., 2009;Parshall et al., 2010;IPCC, 2014). At present 55 % of the global population resides in urban areas (UNDESA, 2014), a number that is projected to rise to 68 % by 2050 (UNDESA, 2018). Meanwhile urban areas cover less than 3 % of the land surface worldwide (Wu et al., 2016) but consume over 66 % of the world's energy (Fragkias et al., 2013) and generate more than 70 % of anthropogenic GHG emissions (Hopkins et al., 2016). Carbon dioxide (CO 2 ) emissions from energy use in cities are estimated to comprise more than 75 % of the global energy-related CO 2 , with a rise of 1.8 % yr −1 projected under business-as-usual scenarios between 2006and 2030(IEA, 2009. Methane (CH 4 ) emissions from energy, waste, agriculture and transportation in urban areas make up approximately 21 % of the global CH 4 emissions (Marcotullio et al., 2013;Hopkins et al., 2016). As emission hotspots, urban areas therefore play a vital role in GHG mitigation. It is crucial to find appropriate methods for understanding and projecting the effects of GHG emissions on urban areas and for formulating mitigation strategies.
There are two methods for the quantitative analysis of GHG emissions: the bottom-up approach and the top-down approach Caulton et al., 2014;Newman et al., 2016). The bottom-up approach calculates emissions based on activity data (i.e., a quantitative measure of the activity that can emit GHGs) and emission factors . This approach has some uncertainty, e.g., on the national fossil-fuel CO 2 emission estimates, ranging from a few percent (e.g., 3 %-5 % for the US) to a maximum of over 50 % for countries with fewer resources for data collection and a poor statistical framework (Andres et al., 2012). The considerable uncertainties are caused by the large variability in source-specific and country-specific emission factors and the incomplete understanding of emission processes (Montzka et al., 2011;Bergamaschi et al., 2015). These uncertainties grow larger at subnational scales, when estimating the disaggregation of the national annual totals in space and time. The top-down approach can not only provide estimated global fluxes but also verify the consistency and assess the uncertainties of bottom-up emission inventories (Wunch et al., 2009;Montzka et al., 2011;Bergamaschi et al., 2018). However, it is hard to quantify the statistical errors attached to both atmospheric observations and prior knowledge about the distribution of emissions and sinks (Cressot et al., 2014). McKain et al. (2012) suggested that column measurements can provide a promising route to improving the detection of CO 2 emitted from major source regions, possibly avoiding extensive surface measurements near such regions. Such measurements, i.e., measurements of concentration averaged over a column of air, are performed to help to disentangle the effects of atmospheric mixing from the surface exchange (Wunch et al., 2011) and decrease the biases associated with estimates of carbon sources and sinks in atmospheric inversions (Olsen and Randerson, 2004). Compared to surface values, urban enhancements in columns are less sensitive to boundary-layer heights (Wunch et al., 2011;McKain et al., 2012;Kivi and Heikkinen, 2016), and column observations have the potential to mitigate mixing height errors in an atmospheric inversion system . Atmospheric GHG column measurements combined with inverse models are thus a promising method for analyzing GHG emissions and can be used to analyze their spatial and temporal variability (Ohyama et al., 2009;Pillai et al., 2011;Ostler et al., 2016;Kivi and Heikkinen, 2016).
In order to focus the top-down approach on concentration differences caused by local and regional emission sources, and in particular to quantify urban emissions, the differential column methodology (DCM) was proposed. It evalu-ates differences between column measurements at different sites. Chen et al. (2016) applied the DCM using the compact Fourier-transform spectrometers (FTSs) EM27/SUN (Bruker Optik, Germany) and demonstrated the capability of differential column measurements for determining urban and local emissions in combination with column models. Citywide GHG column measurement campaigns have been carried out, e.g., in Boston (Chen et al., 2014), Indianapolis , San Francisco, Berlin  and Munich . However, only a few studies have combined differential column measurements with highresolution models. Toja-Silva et al. (2017) simulated the column data at upwind and downwind sites of a gas-fired power plant in Munich using the computational fluid dynamic (CFD) model and compared them with the column measurements. Viatte et al. (2017) quantified CH 4 emissions from the largest dairies in the southern California region, using four EM27/SUNs in combination with the Weather Research and Forecasting model (WRF) in the large-eddy simulation mode. Vogel et al. (2019) deployed five EM27/SUN spectrometers in the Paris metropolitan area and analyzed the data with the atmospheric transport model framework CHIMERE-CAMS.
This paper carries out a quantitative analysis of GHG for the Berlin area in combination with DCM. We utilize the mesoscale WRF model (Skamarock et al., 2008) coupled with GHG modules (WRF-GHG; Beck et al., 2011) at a high resolution of 1 km. The aim is to assess the precision of WRF-GHG and to provide insights on how to detect and understand sources of GHGs (CO 2 and CH 4 ) within urban areas. WRF is a numerical weather prediction system and can be used for both atmospheric research and operational forecasting on a mesoscale range from tens of meters to thousands of kilometers (e.g., Chen et al., 2011). To produce high-resolution regional simulations of atmospheric CH 4 passive tracer transport, WRF was coupled with the Vegetation Photosynthesis and Respiration module (WRF-VPRM; Ahmadov et al., 2007). WRF-VPRM has been widely employed in several studies in which both the generally good agreement of the simulations with measurements and model biases were assessed in detail (Ahmadov et al., 2009;Pillai et al., 2011;Pillai et al., 2012;Kretschmer et al., 2012). Biogenic carbon fluxes given by the VPRM model tend to underestimate urban ecosystem carbon exchange, owing to the incomplete understanding of urban vegetation and to conditions related to urban heat islands and altered urban phenology (Hardiman et al., 2017). WRF-VPRM was later extended to WRF-GHG (Beck et al., 2011), which can simulate the regional passive tracer transport for GHGs (CH 4 , CO 2 and carbon monoxide -CO). Relatively few studies using WRF-GHG have been published as of yet. Pillai et al. (2016) utilized a Bayesian inversion approach based on WRF-GHG at a high spatial resolution of 10 km for Berlin to obtain anthropogenic CO 2 emissions and to quantify the uncertainties in retrieved anthropogenic emissions related to instruments (e.g., CarbonSat) and modeling errors. An observation system simulation experiment was studied in Pillai et al. (2016) based on synthetic data rather than on real observations, as in our study. In the present paper, our focus is on a highresolution (1 km) study of both CO 2 and CH 4 in Berlin and assessing the performance of WRF-GHG through comparing the simulated wind and concentration fields to observations from wind stations and ground-based solar-viewing spectrometers. Then DCM is tested as a proper approach for model analysis, which can cancel out the bias from initialization conditions and highlight regional emission tracers. The simulation workflow is also adapted to this purpose where needed. This study is the fundamental study of the WRF-GHG mesoscale modeling framework in urban areas.
The total annual CO 2 emissions of Berlin (21.3 million t in 2010) approximately correspond to those of Croatia, Jordan or the Dominican Republic (Reusswig and Lass, 2014). With its strong regulatory influence as a state within Germany, and having a strongly supportive policy, Berlin has already transformed itself into a climate-friendly city in which CO 2 emissions have been reduced by a third compared with 1990 levels, aiming for carbon neutrality by 2050 (Homann, 2018). Berlin therefore needs to assess and identify the emission sources accurately at the current stage to provide solid scientific support for the selection of mitigation options. Additionally, Berlin is an ideal pilot case for developing and testing simulations because the city is relatively isolated from other large cities with high emissions, such that anthropogenic GHG anomalies around Berlin can confidently be attributed to the city itself.
The major goals of our work in this context are (1) to simulate high-resolution (1 km) CO 2 and CH 4 concentrations for Berlin using WRF-GHG, attributing the changes in concentrations to different emission processes, (2) to compare the simulation outputs with the observations from a column measurement network in Berlin , assessing the precision of WRF-GHG, and (3) to use DCM in the simulation analysis, testing the feasibly of this approach. The structure of this paper is as follows: the model with its domain and external data sources are described in Sect. 2. A comparison analysis for wind fields and concentration fields is presented in Sect. 3, and CO 2 and CH 4 concentrations related to different processes (e.g., the anthropogenic component) are discussed. DCM, for the comparison of concentration fields and the tracer analysis, is presented and discussed in Sect. 4. Section 5 provides the discussion and summary of this study.

WRF-GHG modeling system
As mentioned in Sect. 1, we use the WRF model version 3.2 coupled with GHG modules to quantify the uptake and emission of atmospheric GHGs around Berlin at a high resolution of 1 km. WRF follows the fully compressible nonhydrostatic Euler equations (Skamarock et al., 2005(Skamarock et al., , 2008 and is based on the actual meteorological data in this case study. The meteorological initial conditions and lateral boundary conditions were taken from the Global Forecast System (GFS) model reanalysis in which in situ measurements and satellite observations were assimilated. Tracers in WRF-GHG are transported online in a passive way, i.e., without any chemical loss or production, when the tracer transport option is used (Ahmadov et al., 2007;Beck et al., 2011). As shown in Fig. 1, three domains are set up here, whose dimensions are 70 × 50 horizontal grid points with a spacing of 9 km for the coarsest domain (d01), 3 km for the middle domain (d02) and 1 km for the innermost domain (d03). WRF uses a terrain-following hydrostatic pressure vertical coordinate (Skamarock et al., 2008). In our case, 26 vertical levels are defined from the surface up to 50 hPa, 14 of which are in the lowest 2 km of the atmosphere. The innermost domain, d03, envelops all five measurement sites (see Sect. 3.1) to assess the simulation by comparing with the measured data. Berlin lies in the North European Plain on flat land (crossed by northward-flowing watercourses), which avoids the vertical interpolation problems caused by topography differences (Fig. 1). The Lambert conformal conic (LCC) projection is selected as a map projection. The simulated time span is from 18:00 UTC on 30 June to 00:00 UTC on 11 July in 2014. The description of the workflow for running WRF-GHG can be found in Appendix A.
The meteorological fields are obtained from the Global Forecast System (GFS) model at a horizontal resolution of 0.5 • , with 64 vertical layers and a temporal resolution of 3 h (as available via the NOAA's National Center for Environmental Information; https://www.ncdc.noaa.gov/, last access: 22 July 2019). The GFS uses hydrostatic equations for the prediction of atmospheric conditions, and its output includes large amounts of atmospheric and land-soil variables, wind fields, temperature, precipitation and soil moisture, etc. The initial and lateral boundary conditions for our WRF-GHG concentration fields are implemented using Copernicus Atmosphere Monitoring Service (CAMS) data (Agusti-Panareda et al., 2017). CAMS provides the estimated mixing ratios of CO 2 and CH 4 , with a spatial resolution of 0.8 • on 137 vertical levels and with a temporal resolution of 6 h (as available via https://atmosphere.copernicus.eu, last access: 22 July 2019).
The simulation of CO 2 and CH 4 fluxes with different emission tracers in WRF-GHG is based on flux models and emission inventories which are either already implemented inside the model modules (online calculation) or constitute external datasets (offline calculation). The flux values from external emission inventories are converted into atmospheric concentrations and added to the corresponding tracer variables. In combination with the background concentration fields for CO 2 and CH 4 that refer to the CO 2 and CH 4 values without any sources and sinks in the targeted domain, the tracer contributions are summed up to obtain the total con- centrations: CO 2,total = CO 2,bgd + CO 2,VPRM + CO 2,anthro + CO 2 , CH 4,total = CH 4,bgd + CH 4,anthro + CH 4,soil + CH 4 , where CO 2,total and CH 4,total represent the total CO 2 and CH 4 , CO 2,bgd and CH 4,bgd are the background CO 2 and CH 4 , CO 2,anthro and CH 4,anthro stand for the changes in CO 2 from the anthropogenic emissions, CO 2,VPRM is the change in CO 2 from the biogenic activities and CH 4,soil is the change in CH 4 from soil uptake, and CO 2 and CH 4 are the tiny computational errors for CO 2 and CH 4 that are described in detail in Appendix B. In the transport process, the relationship shown in Eq. (1) holds for each vertical level. The biogenic CO 2 emission is calculated online using VPRM (Mahadevan et al., 2008), in which the hourly Net Ecosystem Exchange (NEE) of CO 2 reflects the biospheric fluxes between the terrestrial biosphere and the atmosphere, estimated by the sum of gross ecosystem exchange (GEE) and respiration. VPRM in WRF-GHG calculates biogenic fluxes initialized by vegetation indices (land surface water index -LSWI, enhanced vegetation index -EVI, etc.) from the MODIS satellite (as available via https://modis.gsfc.nasa. gov/, last access: 22 July 2019). The reflectance data from the SYNMAP vegetation classification at a resolution of 1 km and 8 d from the MODIS satellite at 0.5-1 km spatial resolution (depending on the wavelength band) are aggregated onto the LCC projection within the VPRM preprocessor. Then, the data, including the high-solution vegetation indices at a resolution of 1 km, are available on the model domains.
We use the external dataset Emission Database for Global Atmospheric Research version 4.1 (EDGAR V.4.1) for the anthropogenic fluxes in our study. EDGAR V.4.1 provides annually varying global anthropogenic GHG emissions and air pollutants at a spatial resolution of 0.1 • Janssens-Maenhout et al., 2015), whose source sectors include industrial processes, on-road and off-road sources in transport, large-scale biomass burning, and other anthropogenic sources (Saikawa et al., 2017). Here we apply time factors for seasonal, weekly, daily and diurnal variations defined by the time profiles published on the EDGAR website (http://themasites.pbl.nl/tridion/en/themasites/edgar/ documentation/content/Temporal-variation.html, last access: 22 July 2019); however, considerable uncertainties are to be expected in applying these time factors. This temporal variation set is derived based on western European data such that the representativity for other European countries and even other world regions may be quite poor. The coarse emission fluxes used for the initialization of the anthropogenic tracer in WRF-GHG can cause problems when locating emission points within the high-resolution model grid and can weaken the impact from the real high-emission hotspots in the fine domain of our study. The chemical sink for atmospheric CH 4 (e.g., photochemistry in the stratosphere) can be ignored in the model, owing to its relatively long lifespan (9.5±1.3 year, Holmes, 2018), the small-scale domains and the limited simulation period (10 d) in our case.
3 Model analysis and model-measurement comparison

Description of measurement sites
The measurement campaign used for comparison with WRF-GHG in this paper was performed from 23 June to 11 July 2014 in Berlin using five spectrometers . It allows us to both test the precision of WRF-GHG (Sect. 3) and verify differential column methodology (DCM) as our analytic methodology (Sect. 4). In their measurement campaign, Hase et al. (2015) used five portable Bruker EM27/SUN FTSs for atmospheric measurements based on solar absorption spectroscopy. Five sampling stations around Berlin were set up, four of which (Mahlsdorf, Heiligensee, Lindenberg and Lichtenrade) were roughly situated along a circle with a radius of 12 km around the center of Berlin. Another sampling site was closer to the city center and located inside the Berlin motorway ring at Charlottenburg (Fig. 6). Detailed information on this measurement campaign is given in Hase et al. (2015), and Frey et al. (2015) provide additional details on the calibration of the spectrometers, precision and instrument-to-instrument biases.

Comparison of wind fields at 10 m
Winds have a strong impact on the vertical mixing of GHGs and a direct influence on their atmospheric transport patterns. Hence, we firstly compare the wind speeds and wind directions obtained from WRF-GHG to the measurements such that deviations between the simulated and measured wind fields are assessed. The wind measurements are not exactly co-located with the spectrometers mentioned in Sect. 3.1, but are rather located at three sampling sites (Tegel, Schönefeld and Tempelhof) and measure at a height of 10 m above the ground. The simulated wind speed at 10 m (ws 10 m ) and wind direction at 10 m (wd 10 m ) are calculated following the equations ws 10 m = u 2 10 m + v 2 10 m , where u 10 m and v 10 m are the components of the horizontal wind, towards the east and north, respectively, which can be obtained from WRF-GHG output files. Figure 2 shows the comparisons of wind speeds (Fig. 2a) and wind directions (Fig. 2b) between simulations and observations at 10 m from 1 to 10 July and the modelmeasurement differences. EM27/SUN only operates in the daytime when there is sufficient sunlight; the detailed description of the instrument can be found in Gisi et al. (2012), Frey et al. (2015) and Vogel et al. (2019). The instrumental working periods are marked by gray shaded boxes in Fig. 2. The measured (dashed lines) and simulated (solid) wind speeds (Fig. 2a) at 10 m show similar trends and demonstrate relatively good agreement over the 10 d time series, with a root-mean-square error (RMSE) of 0.9247 m s −1 . Large uncertainties in wind speeds are found to appear always with the lower wind speeds, mostly at night. In terms of wind directions at 10 m, we observe that the simulated wind directions show similar but slightly underestimated fluctuations (Fig. 2b), which result in an RMSE of 60.8328 • . Larger uncertainties in wind directions always exist during the low wind speed periods (Fig. 2a, b). During the instrumental working period (within the daytime), the simulations fit better with the measurements with relatively lower RMSEs of 0.6928 m s −1 for wind speeds and 41.4793 • for wind directions. We find that the measured wind fields (both wind speeds and wind directions) have more fluctuations compared to the simulations. This could be caused by really fast wind changes which the model, simulating a somewhat idealized environment, is not able to capture. To be specific, local turbulence given by urban canopy, buildings, etc., is not represented well in the model.

Comparison of pressure-weighted column-averaged concentrations
In the following, we use the measured concentration fields to compare with the simulated fields. An FTS EM27/SUN can measure the column-integrated amount of a tracer through the atmospheric column with excellent precision, yielding the column-averaged dry-air mole fractions (DMFs) of the target gases Hedelius et al., 2016). The measured DMFs of CO 2 and CH 4 are denoted by XCO 2 and XCH 4 . Hase et al. (2015) used constant a priori profile shapes in the retrievals of measurements. When comparing remote sensing observations to model data (or also datasets from different remote sensing instruments to one another), limitations of the instruments in reconstructing the actual atmospheric state need to be taken into account. In general, this requires the a priori profile that is used for the retrieval and the averaging kernel matrix, which specifies the loss of vertical resolution (fine vertical details of the actual trace gas profile cannot be resolved) and limited sensitivity (e.g., Rodgers and Connor, 2003). In the case of EM27/SUN, the spectrometers used in the network offer only a low spectral resolution of 0.5 cm −1 . Therefore, performing a simple least-squares fit by scaling retrieval of the a priori profile is generally appropriate. In this case, there is no need to specify a full averaging kernel matrix; instead, the specification of a total column sensitivity is sufficient. The total column sensitivity is a vector (being a function of altitude), which specifies to which degree an excess partial column superimposed on the actual profile at a certain input altitude is reflected in the retrieved total column amount. This sensitivity vector is a function of a solar zenith angle (SZA; and ground pressure), mainly due to the fact that the observed signal levels in different channels building the spectral scene used for the retrieval are shaped by a mixture of weaker and stronger absorptions. (If all spectral lines in the spectral scene are optically thin and too narrow to be resolved by the spectral measurement, the sensitivity would approach unity throughout.) In order to ensure measurement quality and enough sample points for further concentration comparisons, we select five measurement dates (1, 3, 4, 6 and 10 July) with relatively good measurement qualities (from fair, "++", to very good, "++++") based on Hase et al. (2015). The pressuredependent column sensitivities for CO 2 (Fig. 3b) and CH 4 (Fig. 3c) are derived from measurements performed in Lindenberg on 4 July (the best-quality day in terms of measurements). Details about the measurements can be found in Hase et al. (2015) and Frey et al. (2015). The shape and values of the column sensitivities from Lindenberg in Berlin closely resemble the results of Hedelius et al. (2016) in Pasadena. As depicted in Fig. 3a, the SZAs are almost identical for each day in our study (at each hour), rendering the shape of column sensitivities (at a specific hour of the day) practically independent of the measurement date. The column sensitivities for 4 July (Fig. 3b, c) are taken as a basis for our smoothing process below. The a priori CO 2 and CH 4 profiles have been taken from the Whole Atmosphere Community Climate Model (WACCM) version 6. A smoothed profile for a target gas G is then obtained as in Eq. (3) in Vogel et al. (2019), where G is the modeled profile from WRF-GHG, I is the identity matrix, K is a diagonal matrix containing the averaging kernel and G p is the a priori profile. In order to compare the simulated smoothed concentration fields with the observations, the simulated smoothed pressure-weighted column-averaged concentration for a target gas G (XG) is calculated as Here, p i is proportional to the differences of the pressure values P (i) at the bottom and P (i + 1) at the top of the ith vertical grid cell, P top and P sf represent the hydrostatic pressures at the top and at the surface of the model domain, and G s (i) stands for the simulated concentration of the target gas G at the ith vertical level. In Figs. D1 and D2 of Appendix D, we compare the simulated XCO 2 and XCH 4 with and without smoothing. The simulated concentrations are only slightly enlarged after smoothing, at approximately 1-2 ppm for XCO 2 and 2 ppb for XCH 4 , while the variations are mostly unchanged. Compared to the period with lower SZAs (at noon), the smoothed values in the morning and afternoon with higher SZAs hold relatively larger enlargements. Figure 4a shows the measured and smoothed modeled variations in XCO 2 and XCH 4 for these 5 d. Compared to the measurements, the smoothed simulated pressure-weighted column-averaged concentrations for CO 2 (XCO 2 ) show quite similar trends but with approximately 1-2 ppm bias, indicated by an RMSE of 1.2534 ppm. The simulated XCO 2 values are overestimated for 1, 3 and 4 July, while on 6 and 10 July, the model is underestimated, which could be the result of uncertainties from the coarse anthropogenic surface emission fluxes, background concentrations from CAMS (Sembhi et al., 2015) and the ignorance of the influence from the line of sight of the sun. Figure 4b shows the comparison of the pressure-weighted column-averaged concentrations for CH 4 (XCH 4 ) between observations and smoothed simulations on the five selected dates (1, 3, 4, 6 and 10 July). We find that there is an approximate offset of 50-60 ppb between observations and models (RMSE is 58.1082 ppb). The simulated XCH 4 is around 1860 ppb while the measured value is around 1810 ppb, which is comparable to the values (1790-1810 ppb) observed at two Total Carbon Column Observing Network (TC-CON) measurement sites in June and July 2014 in Bremen in Germany  and Bialystok in Poland . This bias of the simulated XCH 4 seems to be constant (around 2.7 %) each day. Thus, we introduce an offset applied to all sites for each simulation date to compare the model and the measured data, effectively removing the bias, which we attribute to too high a back-ground XCH 4 . The daily offset is assumed to be the difference between the smoothed simulated and measured daily mean XCH 4 . After applying the daily offset, the measured XCH 4 shows a somewhat better agreement and a similar trend but with larger variability compared to the simulation (RMSE is 3.1690 ppb). The smaller variations from the simulation results can, for example, be caused by the error from the spatio-temporal treatment of emission maps, underestimated emissions from anthropogenic activities, the coarse wind data and/or the smoothing of actual extreme values in the simulation.
A major offset in modeled CH 4 concentration fields could potentially be attributed to the errors in the troposphere height and a general offset from CAMS. In the CH 4 vertical concentration profile, we find that the typical sharp decrease occurs at the tropopause height. Tukiainen et al. (2016) also find the similar sharp decrease when using the AirCore to retrieve atmospheric CH 4 profiles in Finland. During the simulation, the background concentration values of CAMS are directly fitted to the WRF pressure axis without considering the actual tropopause height; thus this could cause some error. An illustration of the vertical distribution for CH 4 is provided in Appendix C. In contrast, the CO 2 vertical distribution shows decrease that is quite flat with the increase in pressure, and there is no need to consider the tropopause height during the grid treatment in the vertical layer. In terms of CAMS, the reports from Monitoring Atmospheric Composition and Climate (MACC) stated that CAMS has a bias and RMSE (approximately 50 ppb) in each part of the world, compared to the Integrated Carbon Observation System (ICOS) obser- vations in 2017 (Basart et al., 2017). Galkowski et al. (2019) also mentioned one CH 4 offset (approximately 30 ppb within troposphere) when initializing the concentration fields using CAMS. Apart from these two major potential reasons for the bias, the influence from the inaccurate simulated planetary boundary layers and the shape of the constant a priori profile used for the retrievals could both potentially contribute to the discrepancies for the concentration fields. Due to the lack of fine measured vertical concentration profiles, it is not easy to quantify these errors and attribute these potential reasons to this 2.7 % error quantitatively. Thus, a DCM-based analysis is presented in Sect. 4, aiming at eliminating the bias from these relatively high initialization values for CH 4 and making it easier to assess WRF-GHG results with respect to the measurements.

Contributions of different sources and sinks to the total signal: individual emission tracers
As described in Sect. 2, the various flux models implemented in WRF-GHG are advected as separate tracers, making it possible to distinguish the signals in concentration space for different source and sink categories for CO 2 and CH 4 (Beck et al., 2011). Berlin is located in an area of low-lying, marshy woodlands with a mainly flat topography (Kindler et al., 2018). There is no wetland in Berlin according to the MODIS Land Cover Map (Friedl et al., 2010). The land cov-ered by forests, green and open spaces (e.g., farmlands, parks and allotment gardens) accounts for 35 % of the total area in Berlin (SenStadtH, 2016). Additionally, 11 power plants are currently being operated in Berlin, 8 of which have a capacity of over 100 MW (Fraunhofer-Gesellschaft, 2018). In accordance with the geographical characteristics of the district and potential emission sources in Berlin, we focus on understanding the major emissions caused by vegetation photosynthesis and respiration (XCO 2,VPRM ) as well as anthropogenic activities (XCO 2,anthro ) for CO 2 and by soil uptake (XCH 4,soil ) as well as human activities (XCH 4,anthro ) for CH 4 . As an instructive example of an analysis involving these tracers, we look at the diurnal cycle of contributions from the different tracers mentioned above in Charlottenburg (Fig. 5).
The mean values, averaged over 9 d (from 2 to 10 July), as well as a 95 % confidential interval calculated in the averaging process are shown in Fig. 5. Figure 5a clearly shows a decline during the day and a rise at night in the XCO 2 enhancement over the background (blue: XCO 2,total -XCO 2,bgd ), with a maximum decrease over the course of the day of around 2 ppm. The XCO 2 enhancement over the background reaches its daily peak during morning rush hour (07:00 UTC). The morning peak corresponds to XCO 2 changes from human activities, depicted by the black line from 04:00 to 07:00 UTC (marked by a red square in Fig. 5a). Before the evening rush hour (16:00 UTC), XCO 2 over the background then decreases, owning to biogenic uptake. Beginning in the evening, values increase again. The fluctuation in the evening (17:00-19:00 UTC) is dominated by XCO 2 enhancements from human activities, while the substantial rise from 19:00 UTC onward is generated by the VPRM tracer, specifically the accumulation of the vegetation respiration in the evening.
XCO 2 is weaker compared to the strong biogenic uptake. To further highlight the role of anthropogenic activities in XCO 2 changes within the urban area, DCM is applied in Sect. 4. More specifically, we will use downwind-minusupwind column differences of CO 2 ( XCO 2 ) to describe the XCO 2 enhancement over an upwind site, as the difference between the downwind and upwind sites can be attributed to urban emissions.
Turning to XCH 4 in Fig. 5b, we plot the variations in the mean hourly contributions from the anthropogenic (black line: XCH 4,anthro ) and soil uptake tracer (blue line: XCH 4,soil ) in Charlottenburg. The contributions by anthropogenic activities fluctuate slightly around 2 ppb in the morning and at noon; then a peak occurs at the start of the evening rush hour (16:00 UTC). After 18:00 UTC, values clearly decrease, reaching approximately 2 ppb. From 21:00 UTC, XCH 4 stabilizes, exhibiting only moderate fluctuations. The XCH 4 enhancement above the background (green: XCH 4,total −XCH 4,bgd ) depends largely on the XCH 4 contributions by human activities. The changes in concen-trations caused by the soil uptake tracer (blue), whose values fluctuate between 0.001 and 0.01 ppb, have almost no influence on the variation in the XCH 4 enhancement over the background in the urban area.
4 Model analysis using differential column methodology

Comparison of differential column concentrations
The DCM can be employed to detect and estimate local emission sources within an area, based on calculated concentration differences between downwind and upwind sites . The difference ( XG) of a specific gas G in column-averaged DMFs across the downwind and upwind sites is defined as where XG downwind and XG upwind represent the columnaverage DMFs at the downwind and upwind sites. In this study, DCM is applied to measurements and models in the spirit of a post-processing analysis. This approach is not only useful for canceling out the bias of the simulated XCH 4 (see Sect. 3.3) but also for assessing the role of anthropogenic activities in XCO 2 changes more appropriately.
A necessary prerequisite for DCM is distinguishing the upwind and downwind sites among all five sampling sites. Wind direction thus plays a pivotal role in the calculation of the downwind-minus-upwind column differences. In this study, the hourly simulated vertically averaged wind directions are assumed as a standard to classify the sites into downwind and upwind sites. The tracer transport calculations in the first few hours are not stable in WRF-GHG. Thus, we select 3, 4, 6 and 10 July as our targeted dates. Table 1 shows the daily averaged wind directions with standard derivations and the details on the downwind and upwind sites for these four target dates. West wind is the prevailing wind direction on 3 July. That is to say, Mahlsdorf and Lindenberg are downwind sites, and the upwind sites corresponding to these are Charlottenburg and Heiligensee, described in Eq. (6). The wind on 10 July is northeasterly, and the combination of downwind and upwind sites are selected to be opposite of the ones on 3 July, see Eq. (8). The prevailing winds on 4 and 6 July are easterly. The upwind site is Lichtenrade, and the corresponding downwind sites are Heiligensee and Lindenberg, see Eq. (7). Based on the selection of downwind and upwind sites shown in Table 1 and Eq. (5), differential column concentrations ( XCH 4 ) are, therefore, calculated as Western wind (3 July) : Northern wind (4 and 6 July) : Northeastern wind (10 July) : Figure 7 depicts the variations in the wind fields (wind speeds and wind directions) and XCH 4 (corresponding to Eqs. 6, 7 and 8) on 3, 4, 6 and 10 July. As depicted in the Fig. 7ad, the hourly vertically averaged simulated wind speeds and directions at downwind and upwind sites are homogeneous. Thus, it is reasonable to use the daily mean wind directions as  Figure 7. Modeled wind fields for downwind (blue lines) and upwind (red lines) sites (a-d), and downwind-minus-upwind differential evaluation for measured (blue) and simulated (black lines) XCH 4 (e-h) on 3, 4, 6 and 10 July 2014. Based on the selection of downwind and upwind sites in Table 1, XCH 4 is calculated using Eqs. (6), (7) and (8), depicted by blue lines for measurements and black lines for simulations. The black error bars in (e-h) are the standard derivations of the minute values of the hourly mean.
the standard for the selection of downwind and upwind sites. The general trends in the simulated XCH 4 values, shown in Fig. 7e-h, seem to be roughly reproduced by the observations but slightly overestimated, with an RMSE of 1.3895 ppb. Yet DCM as presented here has the potential to highlight the role of anthropogenic activities, which we demonstrate, applying it to CO 2 tracers in the simulation. Thus, the analysis on anthropogenic and biogenic tracers for CO 2 will be especially prominent here. As described above, we continue to take 3, 4, 6 and 10 July as examples (see Fig. 8a-d).
The variations in XCO 2 (corresponding to Eqs. 6, 7 and 8) on 3, 4, 6 and 10 July are shown. In contrast to the variations in XCO 2 values (Sect. 3.4; Fig. 5a), the simulated XCO 2 (Fig. 8a-d, blue lines) is not so much influenced by the XCO 2 changes from the VPRM tracer ( Fig. 8a-d, green) but more closely follows the XCO 2 changes from anthropogenic activities (Fig. 8a-d, red). With DCM, the role of human activities in XCO 2 changes is highlighted, and the strong effect from the biogenic component is canceled out. The XCO 2 measurements ( Fig. 8a-d, black) show similar trends as the simulation with an RMSE of 0.2973 ppm.
To further understand the differences of XCO 2 and XCH 4 between measurements and simulations (see Fig. 7e-h and Fig. 8a-d), the comparison of hourly mean XCO 2 and XCH 4 values for these four targeted dates is illustrated in the right column of Fig. 8. Due to the restriction of measured wind information, we illustrate the differences of simulated and measured wind directions at 10 m (i.e., Fig. 2b) with respect to the hourly mean XCO 2 and XCH 4 . We find that the real hourly mean XCO 2 and XCH 4 values are generally higher than the simulated values. Extreme points are colored by red and blue in the right column of Fig. 8e-f, standing for large differences between measured and simulated wind directions at 10 m. We see that a large difference of wind directions is a necessary but insufficient condition for the bias of XCO 2 and XCH 4 be- Figure 8. Measured (black lines) and simulated (blue lines) XCO 2 on 3, 4, 6 and 10 July 2014, and comparison of hourly mean XCO 2 and XCH 4 for these 4 d. The XCO 2 , calculated using Eqs. (6), (7) and (8), are depicted by blue lines in (a-d). The red and green lines show the variation in the differences between downwind and upwind sites in XCO 2 changes from anthropogenic and biogenic activities, respectively. The points in (e-f) are coded by the difference of the simulated and measured wind directions at 10 m. The black error bars in (a-d) are the standard derivations of the minute values of the hourly mean. tween measurements and simulations. In future studies, this is suggested as something to be verified further.
We conclude that DCM, as applied in this plot, reduces the model bias caused by the simulation initialization but introduces unpleasant effects which may be attributed to errors in the assumed or simulated wind directions.

Comparison between differential column
concentrations and modeling results after the elimination of wind influence As described in Sect. 4.1, the wind direction impacts the distinction between downwind and upwind sites for DCM. Devising meaningful and accurate recipes for determining the wind directions is not easy, sometimes resulting in mixedquality results (of Sect. 4.1). Our simulated output provides the hourly wind and concentration fields. The instruments measure the concentration value every minute . We simply assume the wind direction to be a constant value within 1 h (the hourly vertically averaged values) in our calculation also when it comes to selecting upwind and downwind sites. This may create inaccuracies in the calculation of the measured XCH 4 . In this section, we test replacing the upwind values in DCM by an all-site mean to provide a potential solution for the elimination of such problems while still applying the DCM. The mean of the column-averaged DMFs over all sampling sites (XG specific site ) is assumed to be the back-ground concentration within the entire urban region, replacing the XCH 4 at the upwind site. The differences between the specific site and the mean of all the sites for each gas G ( XG specific site ) is then evaluated, i.e., XG specific site = XG specific site − XG all sites , where XG specific site is the column-averaged DMF at the respective sampling site. We now test this form of DCM for the same four targeted dates (3, 4, 6 and 10 July). The distance between any two sampling sites is around 25 km. The general trends of the simulated (Fig. 9, blue lines) and measured (Fig. 9, black lines) XCH 4 values appear to be more similar with an RMSE of 0.6698 ppb compared to the comparison of XCH 4 in Fig. 7e-h (RMSE of 1.3895 ppb). The modelmeasurement bias can be caused by underestimated emissions from anthropogenic activities, the smoothing of actual extreme values in the simulation and the ignorance of the line of the sun sight for the simulation. The variations in the XCH 4 at the five different sampling sites on the same day are similar (Fig. 9), but the measurements show more extreme values (e.g., 4 July) compared to the simulations. A further analysis in a future study is suggested to provide deeper insight into site-specific transport characteristics.
As a final point in our analysis, we focus on simulated XCO 2 values for these four target dates (Fig. 10). The XCO 2 values (blue line) on 3, 4, 6 and 10 July in five

Discussion and conclusion
We used WRF-GHG to quantitatively simulate the uptake, emission and transport of CO 2 and CH 4 for Berlin with a high resolution of 1 km. The simulated wind and concentration fields were compared with observations from 2014. Then, differential column methodology (DCM) was utilized as a post-processing method for the XCH 4 comparison and the XCO 2 tracer analysis. The measured and simulated wind fields at 10 m mostly demonstrate good agreement, but with slight errors in the wind directions. The simulated pressure vertical profile and the averaging kernel from the solar-viewing spectrometer (EM27/SUN) are used to obtain the smoothed pressureweighted average concentration for further comparisons. The simulated XCO 2 concentrations actually reproduce the observations well, but with approximately 1-2 ppm bias, which can be attributed to the coarse emission inventory, background concentrations from CAMS and the ignorance of the line of the sun sight for the simulation. Compared to the measured XCH 4 , some deviations can clearly be noted in the simulated XCH 4 , mostly caused by the relatively high background concentration fields and the errors at the tropopause height. We discussed the diurnal variation in concentration components corresponding to the major emission tracers for both CO 2 and CH 4 . The biogenic component plays a pivotal role in the variations in XCO 2 . The impact from anthropogenic emission sources is somewhat weak compared to this, while the XCH 4 enhancement is dominated by human activities.
We then concentrated on using DCM for focusing our analysis on relevant CO 2 and CH 4 contributions from the urban area. DCM highlights that the enhancement of XCO 2 over the background within the inner Berlin urban area is mostly caused by anthropogenic activities. In DCM, wind direction plays a vital role in defining the upwind and downwind sites, which directly influence the calculation of differential column concentrations. In the CO 2 tracer analysis, it turns out that XCO 2 , the difference with respect to a mean value instead of a specific upwind site, exhibits a more visible and clearer trend, which proves that the CO 2 enhancement is dominated by anthropogenic activities within the urban area. We conclude that DCM, when applied with care, helps in highlighting the relevant emission sources. Similarly, for XCH 4 , DCM eliminates the bias of the simulated values. Furthermore, when XCH 4 values suffer from in- Figure 10. XCO 2 (blue lines for simulations and black for measurements) for five sampling sites (i.e., the difference between XCO 2 at the site and the mean XCO 2 of five sampling sites): Charlottenburg (a: Char), Heiligensee (b: Heili), Lindenberg (c: Lind), Lichtenrade (d: Licht) and Mahlsdorf (e: Mahls). We furthermore show the differences in the simulated XCO 2 changes from biogenic (green lines) and anthropogenic (red lines) activities. The black error bars in each subplot are the standard derivations of the minute values of the hourly mean. consistent wind directions, we consider XCH 4 to be a useful quantity for analysis.
An analysis of XCO 2 in the Paris hotspot region was carried out by Vogel et al. (2019). Some of their results can be compared to the conclusions we drew in this paper. In Vogel et al. (2019), the modeled XCO 2 was calculated based on the chemistry transport model CHIMERE (2 km) and flux framework CAMS (15 km), with hourly anthropogenic emissions from the IER (Institut für Energiewirtschaft und Rationelle Energieanwendung; University of Stuttgart, Germany) and EDGAR emission inventories and the natural fluxes prescribed by the CTESSEL model (Sect. 2 in Vogel et al., 2019). When comparing results from our simulation, the diurnal variation in the XCO 2 enhancement over the background (Sect. 3.4 and Fig. 5a of our paper) is comparable to the findings of Vogel et al. (2019). For the analysis on the comparison of XCH 4 between simulations and measurements in Sect. 4.1, we found that negative column concentration differences between downwind and upwind sites appear for some periods, owing to the variation in wind directions that causes the conversion of upwind and downwind sites, which was also mentioned for the XCO 2 analysis in Vogel et al. (2019). Based on the CHIMERE-CAMS modeling framework, they showed that the strong decrease in XCO 2 during daytime can be linked to net ecosystem exchange, while a significant enhancement compared to the background is caused by XCO 2 from fossil-fuel emissions, but this is often compensated by net ecosystem exchange. We utilized DCM to bring out the role of anthropogenic activities within urban areas (see the XCO 2 tracer analysis in Sect. 4 of our paper).
We conclude that WRF-GHG is a suitable model for precise GHG transport analysis in urban areas, especially when combined with DCM. DCM is not only useful for the direct evaluation of measurements but also helps us to understand the results of tracer transport models, canceling out the bias caused by initialization conditions, for example, and highlighting regional emission sources. This case is a fundamental study for the WRF-GHG mesoscale modeling framework. Emission flux estimations using WRF-GHG would be our further target to be demonstrated for the case of Munich. This Munich case is combined with the first worldwide permanent column measurement network designed in Munich. Various emission tracers will be run for this case in which more emission tracers (e.g., biogenic emissions from wetland for XCH 4 , traffic emission and strong point source emissions in urban areas) are being separated and analyzed using the longer time period of available measurements.
In future work, we suggest running WRF-GHG for more urban areas such that, for example, different transport, more emission tracers, topography, emission scenarios and the quantification of model errors can be studied. The influence from the line of the sun sight should be taken into account, and the relative sensitivity analysis is suggested. The WRF-GHG mesoscale simulation framework may also be combined with microscale atmospheric transport models to simu-late crucial details of emission sources and transport patterns precisely, with the aim of tracing urban GHG emissions. A further promising direction for future studies may be the application of DCM and model-based analysis to satellite measurements to assess gradients across column concentrations with a dense spatial sampling.
Data availability. The simulation data that support the findings of this study are available on request from the corresponding author. The measurement data are available at https://doi.org/10.5194/amt-8-3059-2015 .

Appendix A: WRF-GHG running process
A detailed description on how to run WRF-GHG is provided in Beck et al. (2011), and thus, only the initialization process for our study in particular is summarized here. One daily simulation with WRF-GHG is normally performed for a 30 h time period, including a 6 h spin-up for the meteorology from 18:00 to 24:00 UTC of the previous day and a 24 h simulation of the tracer transport on the actual simulation day (Beck et al., 2011).
As for the boundary conditions, a small constant offset needs to be added into the WRF boundary files for the biospheric CO 2 and the soil sink CH 4 tracers at the start of each run because these tracers can result in a net sink. When the concentrations become negative, the advected tracer fields will "disappear", as the WRF code does not allow tracers with negative values. An offset applied in the initialization process helps to avoid this problem and later is subtracted in the post-processing. As for the initial conditions, the meteorological conditions are initialized with external data sources (GFS in our model) each day to update the WRF meteorological fields properly. The tracers for the total and background CO 2 and CH 4 flux fields are initialized only once, at the first day of the simulation period, using CAMS as an external data source. Furthermore, the lateral boundary conditions of the outer domain d01 are also initialized by the CAMS. Then, for the other days within the simulation period, these tracers for the total and background CO 2 and CH 4 fluxes are directly taken from the final WRF output at 24:00 UTC of the previous day to make the entire simulation continuous. The CO 2 tracer for VPRM and the CH 4 tracer for soil uptake are also initialized with a constant offset to avoid the appearance of negative values caused, for example, by the vegetation respiration (Beck et al., 2011). In terms of the other flux tracers, the tracer variables are initialized each day, using external data sources to provide the updated emission data for each tracer.

Appendix B: Model systematic equation errors for Eq. (1)
In the passive tracer transport simulation, the total concentration of each GHG is represented as a separate tracer, giving redundant information (with respect to the sum of all tracers for each GHG) and allowing for consistency checks. A variety of flux models and emission inventories implemented in the modules of WRF-GHG are used for the estimation of GHG fluxes. The flux values from external emission inventories are gridded and absorbed into the model. In the transport process, the relationship among the changes in concentrations from different emission tracers, the total and background concentrations (Eq. 1) should then be satisfied, ideally with CO 2 and CH 4 computational errors during the simulation process being zero. Nonzero values of CO 2 and CH 4 reflect the limited precision of the tracer transport calculation in WRF-GHG. Figure B1. The mean values (solid lines) and the 95 % confidence intervals of the computational error CO 2 (a) and CH 4 (b). CO 2 and CH 4 are calculated using Eq. (1). Figure B1 thus shows the mean values (solid lines) and the 95 % confidence intervals of CO 2 and CH 4 . As depicted in the figure, CO 2 ranges from −0.005 to 0.01 ppm, while CH 4 is in the range of −0.01 to 0.02 ppb. Divided by typical absolute values of the concentrations from different flux processes for XCO 2 (around 1 ppm) and XCH 4 (around 2-3 ppb) depicted in Fig. 4, the relative computational error is found to be ∼ 1 % for both CO 2 and CH 4 .
These tiny computational errors can be caused by the slight non-linearity of the advection scheme used in the WRF-GHG model, which makes the sum of the concentrations in CO 2 and CH 4 from all individual flux tracers not exactly equal to the concentration from the sum tracer, representing the total sum of all fluxes related to different processes.
Appendix D: Accounting for instrumental limitations in comparison of measured to simulated XCO 2 and XCH 4 Figure D1. Comparison of XCO 2 from WRF-GHG with and without smoothing (using our column sensitivities for EM27/SUN) for the first four simulated dates. The five colors stand for the concentrations from five sample sites. Dotted lines with the crosses represent the XCO 2 without smoothing, while solid lines with the circles stand for the smoothed values. Figure D2. Comparison of XCH 4 from WRF-GHG with and without smoothing (using our column sensitivities for EM27/SUN) for the first four simulated dates. The five colors stand for the concentrations from five sample sites. Dotted lines with the crosses represent the XCH 4 without smoothing, while solid lines with the circles stand for the smoothed values. 11298 X. Zhao et al.: Analysis of total column CO 2 and CH 4 measurements in Berlin with WRF-GHG Appendix E: The vertical wind profiles for wind speeds and wind directions Figure E1. The vertical distribution of wind fields (wind speeds and wind directions) on 3 July (a, b) and 4 July (c, d) in Tegel. The colors from black to blue represent the time from morning to evening.
Author contributions. XZ performed the simulations, with the support and guidance of JM, CG, JC and SH. JM provided the CAMS fields for the initialization. JC supplied the anthropogenic emission source, and CG offered the VPRM used for the simulations. MF and FH provided the measurement data for Berlin in 2014 and fruitful discussions related to the measurements. SH provided the guidance related to the running of the simulations in the Linux cluster. XZ, JC and SH designed the computational framework. XZ and JC performed the analysis of the results. XZ wrote the paper, with input from all authors. All authors provided critical feedback and helped shape the research, analysis and paper. istry for the biogenically related CH 4 flux estimates. The a priori concentration profiles from the Whole Atmosphere Community Climate Model (WACCM) were provided by James W. Hannigan (NCAR). Jia Chen is partly supported by the Technical University of Munich Institute for Advanced Study, funded by the German Excellence Initiative and the European Union Seventh Framework Programme under grant agreement no. 291763. The simulations presented in this work have been run on the Linux cluster (CooLMUC-2) of the Leibniz Supercomputing Centre (LRZ; Garching). We acknowledge support by the ACROSS research infrastructure of the Helmholtz Association.
Financial support. The article processing charges for this openaccess publication were covered by the Max Planck Society.
Review statement. This paper was edited by Stefano Galmarini and reviewed by two anonymous referees.