Recommended coupling to global meteorological fields for long-term tracer simulations with WRF-GHG

Atmospheric transport models are often used to simulate the distribution of greenhouse gases (GHGs). This can be in the context of forward modeling of tracer transport using surface–atmosphere fluxes or flux estimation through inverse modeling, whereby atmospheric tracer measurements are used in combination with simulated transport. In both of these contexts, transport errors can bias the results and should therefore be minimized.

Here, we analyze transport uncertainties in the commonly used Weather Research and Forecasting (WRF) model coupled with the greenhouse gas module (WRF-GHG), enabling passive tracer transport simulation of CO₂ and CH₄. As a mesoscale numerical weather prediction model, WRF's transport is constrained by global meteorological fields via initialization and at the lateral boundaries of the domain of interest. These global fields were generated by assimilating various meteorological data to increase the accuracy of modeled fields. However, in limited-domain models like WRF, the winds in the center of the domain can deviate considerably from these driving fields. As the accuracy of the wind speed and direction is critical to the prediction of tracer transport, maintaining a close link to the observations across the simulation domain is desired. On the other hand, a link that is too close to the global meteorological fields can degrade performance at smaller spatial scales that are better represented by the mesoscale model. In this work, we evaluated the performance of strategies for keeping WRF's meteorology compatible with meteorological observations. To avoid the complexity of assimilating meteorological observations directly, two main strategies of coupling WRF-GHG with ERA5 meteorological reanalysis data were tested over a 2-month-long simulation over the European domain: (a) restarting the model daily with fresh initial conditions (ICs) from ERA5 and (b) nudging the atmospheric winds, temperatures, and moisture to those of ERA5 continuously throughout the simulation period, using WRF's built-in four-dimensional data assimilation (FDDA) in grid-nudging mode.

Meteorological variables and simulated mole fractions of CO₂ and CH₄ were compared against observations to assess the performance of the different strategies. We also compared planetary boundary layer height (PBLH) with radiosonde-derived estimates. Either nudging or daily restarts similarly improved the meteorology and GHG transport in our simulations, with a small advantage of using both methods in combination. However, notable differences in soil moisture were found that accumulated over the course of the simulation when not using frequent restarts. The soil moisture drift had an impact on the simulated PBLH, presumably via changing the Bowen ratio. This is partially mitigated through nudging without requiring daily restarts, although not entirely alleviated. Soil moisture drift did not have a noticeable impact on GHG performance in our case, likely because it was dominated by other errors. However, since the PBLH is critical for accurately simulating GHG transport, we recommend transport model setups that tie soil moisture to observations. Our method of frequently re-initializing simulations with meteorological reanalysis fields proved suitable for this purpose.

1 Introduction

Quantification of carbon sources and sinks is an area of active scientific research with implications for global warming and climate change (IPCC AR6; Rama et al., 2022). In the context of the Paris Agreement, a dramatic reduction in carbon emissions is planned in order to limit the rising temperatures seen around the globe. Monitoring these emission reductions is key to the success of this global effort. One method to achieve this monitoring is via atmospheric inverse modeling, wherein atmospheric transport models are used to deduce anthropogenic greenhouse gas (GHG) emissions based on atmospheric measurements of these GHGs. Uncertainties within the inversion approach include prior emission fluxes, observation error (mainly from satellite-based measurements), and transport modeling errors (Feng et al., 2016). In the case of regional inversions, uncertainties in the lateral boundary conditions (LBCs) also need to be considered (Schuh et al., 2010; Feng et al., 2019a, b; Chen et al., 2019). Previous studies suggest that transport errors can have a considerable impact on simulated atmospheric GHG mole fractions (Lin, 2005; Díaz Isaac et al., 2014) and flux estimates (Baker et al., 2006; Lauvaux et al., 2009; Lauvaux and Davis, 2014). The spatial resolution of the model also plays a role in emission estimates, as pointed out in Feng et al. (2019b) and Gerken et al. (2021). An increase in resolution may be of benefit when compared with lower-resolution systems such as global models. In a study focusing on model uncertainties in regional atmospheric CO₂ simulations over North America, Feng et al. (2019a) found that transport uncertainties were as large as uncertainties due to biogenic fluxes in some seasons, which should be considered in the design and interpretation of inversion studies. Given these facts, efforts must be undertaken to reduce transport errors. As transport errors may accumulate following model initialization, this is especially important for longer simulations.

Transport models are driven by modeled meteorology fields. To reduce transport errors, one needs to use validated meteorological fields with high accuracy. The European Centre for Medium-Range Weather Forecasts (ECMWF) provides various datasets suitable for this purpose, such as reanalyses, analyses, and forecasts at different spatial resolutions, up to 0.08°×0.08° globally. In these meteorological fields, variables such as temperature, humidity, horizontal winds, and vertical mixing have a crucial impact on how accurate atmospheric transport models can be. Reanalysis data in particular, despite typically being produced at lower spatial resolution than forecast products, are informed by historical observations and are quality controlled, which is desirable for the task.

Despite the constraint from lateral boundary conditions, if a model is initialized with a unique initial condition (IC) and is allowed to run freely for periods of months or years, the simulated meteorology will deviate from the driving fields and therefore from reality. This is because what happens within the model domain depends on various components of the model's physics, land surface scheme, and parameterizations of subgrid-scale processes. To carry out atmospheric GHG transport simulations while maintaining a link to the reanalysis data, one of the simplest approaches is to frequently re-initialize the model with fresh meteorological initial conditions (Ahmadov et al., 2007, paragraph 25). This approach was adopted in previous regional studies using the Weather Research and Forecasting (WRF) mesoscale model (Skamarock et al., 2008) coupled with the greenhouse gas module (WRF-GHG; Beck et al., 2011). Specifically, regular re-initialization each day at 00:00 UTC, coupled with a 6 h meteorological spin-up starting at 18:00 UTC the previous day, was used in Pillai et al. (2011), Beck et al. (2013), Gałkowski et al. (2021), and others. This strategy is similar to the common practice in numerical weather prediction (NWP), where short-range forecasts are performed over limited periods due to the growth of the forecast error over time. Such an accumulation of errors might also be expected within the regional domain (Simmons et al., 1995; Molteni et al., 1996; DelSole and Hou, 1999; Danforth et al., 2007). Another method to avoid drifting from the driving fields within a regional model is to apply nudging inside the domain (Stauffer and Seaman, 1990), ensuring that the simulated meteorological fields do not deviate too far from the global fields while the simulation is in progress. Conveniently, grid nudging is one of the built-in four-dimensional data assimilation (FDDA) options within WRF (see Sect. 2.3).

Hence, some studies rely on nudging instead of frequent re-initialization of the model (Bullock et al., 2014; Spero et al., 2014; Markina et al., 2018; Zittis et al., 2018), e.g., in order to avoid discontinuity in the simulated transport (Lo et al., 2008; Vincent and Hahmann, 2015).

Few studies have assessed whether it is necessary to re-initialize periodically or how frequently such re-initializations should take place. Others have shown that, in long-term continuous simulations with a focus on meteorological variables (such as winds, temperature, pressure, and precipitation), these issues can be resolved when nudging is applied (Lo et al., 2008; Vincent and Hahmann, 2015), therefore suggesting that there is no need to do multiple short runs and that continuity of the fields is ensured. However, to our knowledge, no studies have investigated the impact on long-term transport of GHG tracers, including the benefits and drawbacks of combining frequent re-initialization and nudging for long-term transport of GHGs.

The goal of this study is to determine the optimal method for keeping long-term GHG tracer simulations close to real weather. To this end, we tested various strategies to keep WRF close to ERA5 meteorological fields, i.e., by re-initializing WRF daily, grid nudging, or both. This assessment focuses on the accuracy of meteorological parameters that are critical for GHG tracer transport, namely wind speed, wind direction, and planetary boundary layer height (PBLH), as well as soil moisture, surface temperature, and humidity, as they impact the PBLH. In addition, we simulate and evaluate transport of CO₂ and CH₄. The simulated values for both meteorological parameters and GHGs were compared to observations across the model domain. This paper is structured as follows: Sect. 2 provides detailed information about the model setup, including the products used for meteorological and chemical initial conditions (ICs) and lateral boundary conditions (LBCs), the technique of grid nudging, the frequent restart strategy, the experimental design, and the observational data used for model validation. Section 3 contains the results and statistical analysis. We discuss our main findings and compare them to similar studies in Sect. 4. Finally, a summary and conclusion is given in Sect. 5.

2 Data and methods

2.1 WRF-GHG

The core of our system is the Weather Research and Forecasting model (WRF) version 3.9.1.1 (Skamarock et al., 2008), run with the Advanced Research WRF (ARW) core. We use WRF-Chem, enabling the GHG option, which we will refer to as WRF-GHG in the text. The model uses fully compressible Eulerian non-hydrostatic equations on an Arakawa C-staggered grid, conserving mass, momentum, entropy, and scalars (Skamarock et al., 2008; Mahadevan et al., 2008). More than 1 decade ago, WRF was first used to simulate atmospheric CO₂ at the mesoscale (Ahmadov et al., 2007). Modules bundled with the WRF distribution, referred to as WRF-Chem (Grell et al., 2005), allow chemical processes and the transport of atmospheric tracers to be simulated, including gases and aerosols. Such a module was prepared by Beck et al. (2011) as the greenhouse gas module (WRF-GHG), which we use in our study to simulate CO₂ and CH₄ transport. Our simulations use the Noah land surface model (LSM) from Chen and Dudhia (2001), stochastic convective parameterization from Grell and Freitas (2014), and the MYNN2 planetary boundary layer (PBL) scheme (Nakanishi and Niino, 2006). In WRF-GHG, to simulate changes in mole fractions of CO₂ and CH₄, offline or online coupling to flux models can be used. These fluxes are associated with anthropogenic and biogenic flux components that are transported into the atmosphere to simulate the integrated signal of CO₂ and CH₄.

We compute biogenic CO₂ fluxes using the online Vegetation Photosynthesis and Respiration Model (VPRM; Mahadevan et al., 2008). VPRM utilizes remote-sensing products: the Enhanced Vegetation Index (EVI) and the Land Surface Water Index (LSWI) derived from reflectances measured by the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite (http://modis.gsfc.nasa.gov/, last access: December 2020). These indices are aggregated for various vegetation types and then projected onto the model domain at the spatial resolution of the transport model. They are then combined with model-simulated solar radiation and 2 m temperature, simulating the biogenic uptake and respiration of CO₂ at the resolution of the model time step. Therefore, biogenic CO₂ fluxes differ among the experiments.

Anthropogenic CO₂ and CH₄ fluxes were taken from emission inventories. CH₄ emissions are from the Emissions Database for Global Atmospheric Research (EDGAR) dataset, version 4.3.2 (Janssens-Maenhout et al., 2019), with a horizontal resolution of 0.1°×0.1°; we have used emissions from 2012 (latest available in that product). CO₂ emissions were taken from the European TNO-MACC-III inventory with a spatial resolution of 0.125°×0.0625°, which is an update of the earlier TNO-MACC-II dataset (Kuenen et al., 2014). Two different inventories were used for the two species because previous analysis (not shown) found that these inventories resulted in a better agreement with measurements. All the emission products were re-gridded and interpolated to match the spatial resolution of the WRF domain and disaggregated from the available annual emission data into hourly emissions, using country- and sector-specific temporal and vertical profiles of GHG emissions (Brunner et al., 2019).

2.2 Initial and boundary conditions

For meteorological initial and lateral boundary conditions, we use ERA5 reanalysis fields (horizontal winds, pressure, temperature, sea surface temperature etc.) from ECMWF (Hersbach et al., 2020). These 3D fields were downloaded at hourly resolution on a 0.25°×0.25° regular grid (approximately 31 km spatial resolution) on the 137-level ECMWF vertical grid using the Climate Data Store Application Program Interface (https://cds.climate.copernicus.eu/, last access: December 2020, European Reanalysis 5, 2020). The data were then interpolated to WRF grids using the WRF Preprocessing System (WPS) software, which is part of the regular WRF distribution.

Initial and boundary tracer conditions for CH₄ and CO₂ were taken from the Copernicus Atmosphere Monitoring Service (CAMS) GHG short-term forecast (experiment ID: gqpe, based on IFS cycle CY43R1; see Diamantakis and Agusti-Panareda, 2017, and Agusti-Panareda et al., 2017), henceforth referred to as the CAMS fields. This product benefits from the assimilation of satellite observations at the initialization stage (from TANSO-GOSAT for CO₂ and CH₄ and also from MetOp-IASI for CH₄). Gałkowski et al. (2021) showed small bias errors in this product for the free troposphere over Europe, especially for CH₄. For the full period of our simulation, we have used the first 24 h of the forecast, initialized at midnight at 3 h temporal resolution on the ECMWF L137 vertical grid. We performed a horizontal interpolation from the original TCo1279 Gaussian cubic octahedral grid (equivalent to approximately 9 km horizontal resolution) to an intermediate 0.125°×0.125° regular latitude–longitude grid. The 3D CH₄ and CO₂ fields were then interpolated onto the WRF grid.

2.3 Grid nudging

Nudging, also known as Newtonian relaxation, is a method of four-dimensional data assimilation (FDDA). This technique keeps the model close to an analysis field over the course of a simulation. Grid nudging gently forces the model simulation towards a series of physical reference states by adding a calculated relaxation term at every model grid point. It can be applied selectively, i.e., above a given model level, over a selected period, or throughout the simulation. This method provides a four-dimensional analysis that is moderately balanced dynamically with driving meteorological fields and preserves continuity while still allowing complex local topographical or convective variations that are resolved by the model subroutines.

These additional tendency terms are calculated by the difference in the model state and the (re-)analysis state in the nudged variable in addition to the normal tendency term originally derived by the model throughout the domain or at selected altitudes and at each time step, as shown in Eq. (1):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}=F\left(\mathit{\theta }\right)+G_{\mathit{\theta }}{W}_{\mathit{\theta }}({\widehat{\mathit{\theta }}}_{\mathrm{0}}-\mathit{\theta }),\end{array}$$

where θ represents the variable fields that are nudged (two horizontal wind components (u and v), temperature, and moisture), while $\widehat{{\mathit{\theta }}_{\mathrm{0}}}$ is the value of the corresponding variable to which the nudging relaxes the solution (in this case, the interpolated reanalysis value from ERA5). F(θ) is the tendency term obtained by the model's parameterization of physics, advection, etc.; G_θ is a timescale constant controlling the nudging strength (nudging coefficient); and W_θ is an additional spatiotemporal weight used to limit the effect of nudging. The nudging strength G_θ should be carefully selected such that modeled features are not overwritten.

For horizontal winds and temperature, G_θ was set to the default value of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹, as several studies found it to be acceptable (Spero et al., 2018). For moisture, a value of $\mathrm{4.5}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹ was used, following Spero et al. (2018), who found that the default value was too high. Nudging is turned off in the planetary boundary layer (PBL) so that local-scale features near the surface are allowed to develop within the regional model (Miguez-Macho et al., 2004; Lo et al., 2008; Bowden et al., 2012, 2013). Spero et al. (2014) found that capping nudging at the tropopause (while also restricting nudging below the PBL) improved the representation of radiation, clouds, and precipitation in WRF, which potentially affects the skill of simulating atmospheric CO₂, CH₄, and their fluxes (this is discussed in more detail in Sect. 4). Nudged simulations were all started from a single pre-simulation spin-up period (i.e., −6 h to 00:00 UTC), in which grid nudging was applied.

Other built-in nudging methods include spectral and observational nudging. Moreover, grid- and observational nudging can be used separately or in combination. For this study we focus primarily on grid nudging, but we also assessed spectral nudging. However, results from spectral nudging were not significantly different; thus it is not included in our analysis.

2.4 Experiment design

The model domain covers Europe, spanning from roughly 30° W to 55° E and 33 to 67° N, as shown in Fig. 1. The horizontal resolution of the grid is 5 km×5 km with 882×705 grid points. Simulations were performed using 60 vertical levels, with the model top at 50 hPa and 10 levels within the lowest 2 km. The internal time step of the model was 15 s, and we used instantaneous values stored hourly in our analysis. WRF-GHG simulations were performed for May and June 2018.

Figure 1 Single WRF domain over Europe with 5 km×5 km grid spatial resolution. Red circles with three-letter labels mark the locations of ICOS tall-tower sites for evaluating simulated GHG tracer mole fractions, blue diamonds show the IGRAv2 network of radiosondes for PBLH estimation, and yellow pentagons are the surface synoptic stations for assessing meteorological skills near the surface in Upper Silesia. Flight tracks of in situ CH₄ aircraft measurements, which were used for model–observation comparison in the Upper Silesian Coal Basin, are shown in the inset on the right.

Two main strategies for assuring the consistency of the WRF-GHG meteorological fields with those of ERA5 were tested. The first was to regularly (every 24 h) restart the model with fresh initial conditions from ERA5 (referred to as “DR” for “daily restarts”). The second strategy nudged the atmospheric winds, temperatures, and moisture to ERA5 values, continuously using the built-in four-dimensional data assimilation (FDDA) option (referred to as “GN” for “grid nudging”). In the DR simulations, the model is restarted each day at 00:00 UTC with ERA5 meteorological fields initialized at 18:00 UTC the previous day, following a spin-up period of 6 h to allow downscaling of the variables consistent with the WRF physics. CO₂ and CH₄ tracer fields at 00:00 UTC are then copied over from the previous day's simulation (Ahmadov et al., 2012).

Six different simulations were conducted using combinations and variants of these two main strategies, as shown in Table 1. A simulation with no nudging (NN) and no restarts (NR) served as the reference simulation (combined to give NN_NR). For simulations with grid nudging, we applied nudging only above the boundary layer, so that local features below were allowed to develop without interference by the coarser global analysis fields. We employed two methods for turning off nudging inside the PBL, nudging only above the PBL dynamically determined by the PBL scheme (simulations denoted simply as GN) and nudging only above a certain fixed level, to avoid uncertainties in predicting the PBLH (Díaz-Isaac et al., 2019). In this case, we chose 3 km or 700 hPa (at the 13th model level), such that the simulations were mostly above the boundary layer. These simulations have a suffix, resulting in GN_3km.

Table 1 Configuration of assimilating ERA5 for WRF simulations.

2.5 Validation data and methods

We compare modeled meteorological variables with observations from surface synoptic stations and radiosonde-derived PBLH estimates. Meteorological observations used for model validation are 10 m horizontal winds, 2 m temperature, and specific humidity; these were taken hourly from the NOAA Integrated Surface Database (ISD; https://registry.opendata.aws/noaa-isd, last access:: May 2021). A subset of the network in the Upper Silesian Coal Basin (USCB), Poland, was used for a total of 12 stations including two mountain sites (Zakopane and Kasprowy Wierch), as shown in Fig. 1.

Additionally, we evaluated our simulations against radiosonde data from the Integrated Global Radiosonde Archive (IGRA) Version 2 (Durre et al., 2016). This network consists of quality-controlled radiosonde observations of temperature, humidity, and wind at stations across all continents. There are 22 stations in the IGRA v2 network within our model domain, as shown in Fig. 1. Only data at 12:00 UTC were used for estimating the PBLH, derived using the bulk Richardson method (Vogelezang and Holtslag, 1996, Eq. 2) as described by Seidel et al. (2010, 2012). Thereafter, we analyzed simulated soil moisture (SMOIS), i.e., the differences among simulations, and demonstrate its sensitivity to the PBLH.

Apart from the meteorological evaluations, we also compared modeled GHG mixing ratios with two observation datasets. The first dataset comprised CO₂ and CH₄ measurements from 18 instrumented Integrated Carbon Observation System (ICOS; ICOS RI et al., 2022) tall towers. We used hourly measurements from the highest intake level of each tower, collected between 11:00 and 15:00 UTC to ensure a well-developed PBL. The second GHG dataset consisted of low-altitude aircraft GHG measurements that were obtained in May–June 2018 during the CoMet 1.0 campaign (Fiehn et al., 2020b) in the Upper Silesian Coal Basin (USCB) in southern Poland, where multiple coal mines operate. Due to the large number of mines characterized by high specific methane emissions (defined as the amount of methane emitted per 1 t of excavated coal; Swolkień et al., 2022) and continuously high levels of mining activity, this area is considered to be one of the major anthropogenic sources of CH₄ in Europe, responsible for emissions of approximately 475 kt CH₄ yr⁻¹ (CoMet ED v4.01; Swolkień et al., 2022). The aircraft, equipped with an in situ analyzer for greenhouse gases, flew upwind and downwind over the mine source cluster, capturing anthropogenic signals emitted from the mines for several hours on different days during May and June 2018 under clear weather conditions. We focus the analysis on CH₄ as an example of a strong GHG source in the near-field.

The latter two are independent observations to evaluate the performance of each model setup, whereas ERA5 has already assimilated data from various surface weather sites and radiosondes. Due to the availability of aircraft GHG data during this period, the comparisons against surface synoptic stations are exclusively focused on the Upper Silesian Coal Basin (USCB). The location of these measurements can also be seen in Fig. 1.

We evaluate our simulations based on how well they reproduce the observations described above. The performance metrics are mean error (ME), root-mean-square error (RMSE), and coefficient of determination (R²). In some cases, we show Taylor diagrams which provide additional information.

3 Results

3.1 Meteorological validation

In order to assess meteorological model performance, we compare the results with measurements at a subset of synoptic surface stations monitoring meteorological conditions over southern Poland. We evaluate the performance of the six scenarios mentioned in Table 1 by applying the statistical metrics on simulated and observed wind speed, wind direction, temperature, and specific humidity near the model surface, namely 10 m winds (U₁₀, V₁₀), 2 m temperature (T₂), and specific humidity (Q₂); ME and RMSE were computed using data from the 12 stations in the Upper Silesian region (see locations in Fig. 1). The results are shown in Figs. 2 and 3.

Figure 2 Box plots of statistical errors across 12 synoptic stations in Silesia analyzed hourly for May and June 2018. The whiskers extend a distance of 1.5 times the interquartile range. The top of the box represents the 25th percentile, the bottom of the box is the 75th percentile, the red line represents the median, and the white star indicates the mean. Outliers are plotted with the symbol “x”.

Figure 3 Box plots (same as in Fig. 2) of the coefficient of determination R² across 12 synoptic stations in Silesia analyzed hourly for May and June 2018.

With daily re-initialization introduced to the reference run (NN_NR), simulation NN_DR reduced the average RMSE for 10 m wind speed by 14.31 % (from 1.93 to 1.65 m s⁻¹; Fig. 4a), for T₂ by 38.16 % (from 3.04 to 1.88 K; Fig. 4c), and for Q₂ by 48.28 % (from 1.97 to $\mathrm{1.02}\times {\mathrm{10}}^{-\mathrm{3}}$ kg kg⁻¹; Fig. 4d). Grid nudging (GN_NR) achieved slightly better error reductions with 17.61 %, 43.09 %, and 52.79 % for 10 m wind speed, T₂, and Q₂ in RMSE (from 1.93 to 1.59 m s⁻¹, from 3.04 to 1.73 K, and from 1.97 to $\mathrm{0.93}\times {\mathrm{10}}^{-\mathrm{3}}$ kg kg⁻¹), respectively. A similar or slightly better error reduction than GN_NR relative to the reference run (NN_NR) is achieved with simulation GN_DR, with the average RMSE of 10 m wind speed falling from 1.93 to 1.58 m s⁻¹ (18.13 %) and the mean RMSE of T₂ and Q₂ dropping from 3.04 to 1.72 K (43.31 %) and from 1.97 to $\mathrm{0.92}\times {\mathrm{10}}^{-\mathrm{3}}$ kg kg⁻¹ (53.28 % error reduction), respectively.

Figure 4 Statistical metrics of model performance against radiosonde-derived estimates of the PBLH at 12:00 UTC for May and June 2018. The map shows the locations of the radiosonde stations used for validation. The Taylor diagram (a) provides a summary of the skill for each simulation, showing normalized standard deviation, RMSE, and correlation coefficient, averaged across all stations. Markers that lie closer to the red star (labeled as “Observed”) show better performance. Panels (b) and (c) are as in Figs. 2 and 3.

Simulation GN_3km_NR, while outperforming the reference run in most metrics, has poorer performance compared to NN_DR. This is especially the case in representing T₂ and Q₂: compared to the reference run, the mean RMSE is only reduced to 2.40 K (27.65 % larger error than NN_DR) and $\mathrm{1.54}\mathrm{kg}{\mathrm{kg}}^{-\mathrm{1}}\times {\mathrm{10}}^{-\mathrm{3}}$ (50.7 % larger error than NN_DR), respectively. A similar trend can be seen in the average ME and in the 10 m wind speed, T₂, and Q₂. We do not observe significant variability in wind direction performance across the six different simulations. Figures 2 and 3 show the performance metrics of the four variables that are assimilated from ERA5 using FDDA. Together with statistical measures found earlier and the highest R² score, GN_DR has the best performance in general among the scenarios, with GN_NR ranking a close second.

Figure 4 shows an overview of the PBLH performance, similar to the analysis of the synoptic surface stations. The largest error reduction is again seen in simulation GN_DR, dropping 24.64 % from 618.31 to 465.95 m in RMSE compared to the reference run, closely followed by GN_NR with 21.90 % (from 618.31 to 482.4). Similar to the evaluation with surface meteorological data, GN_3km_NR also shows the smallest reduction in RMSE and is slightly worse than NN_DR (524.93 and 502.80 m, respectively). This may be explained by the poor performance in simulating T₂ and Q₂, since both parameters drive the development of the PBL. All statistical results from this section are summarized in Table A1.

3.2 Evolution of soil moisture

We observe a divergence of SMOIS over time between scenarios with daily restarts (DR) and no restarts (NR) (Figs. 5 and 6). Soil moisture is modeled by the land surface model component of WRF (here: Noah). The key difference between the DR and NR simulations is that, in the DR scenarios, soil moisture is re-initialized every 24 h from ERA5 and thus remains close to the land surface model from the ECMWF Integrated Forecast System (IFS, HTESSEL ECMWF, 2016). In contrast, for the NR scenarios, SMOIS solely follows the course of the land surface model in WRF and, over time, drifts away from the ERA5 results (Fig. 6). In Sect. 4.3, we discuss the implications of divergent SMOIS.

Figure 5 Difference in soil moisture between the simulations (GN_3km_DR – GN_3km_NR) using grid nudging with daily restarts (GN_3km_DR) and using grid nudging with no restarts (GN_3km_NR) above 3 km.

Figure 6 Time series showing the drift in soil moisture between the simulations using grid nudging with daily restarts (GN_DR) and using grid nudging with no restarts (GN_NR). Within the model domain, six different locations/pixels are selected to demonstrate such discrepancies throughout the simulation period.

3.3 Evaluation of WRF-GHG tracer simulations against ICOS-ATM

Figure 7 summarizes the statistical evaluation of CO₂ and CH₄ against ICOS data. In Fig. 7a, although we observed a rather low ME for the reference scenario, we do see that the RMSE is higher than that for all five other scenarios. This is a sign of compensating errors of both signs. From the Taylor diagram for atmospheric CO₂ and CH₄ (Fig. 7c and f), we observed a cluster separated from the reference run in the direction of the observations (red star). This cluster consists of five simulations which all used either daily restarts and/or grid nudging. Interestingly, we do not see any significant distinction among these five scenarios. The simulations exhibit only small differences in the performance metrics for both atmospheric CO₂ and CH₄, as seen in the box plots of statistical errors and R² in Fig. 7b and e. By a very small margin, GN_DR yields the best performance and GN_3km_NR yields the worst, except for the reference run. This is the same result as in the meteorological evaluation in Sect. 3.1, albeit with smaller differences among all runs. For detailed statistical results for all six scenarios, readers are referred to Table A1.

Figure 7 Statistical overview of model performance against ICOS tall-tower observations of atmospheric CO₂ and CH₄ analyzed hourly from 11:00–15:00 UTC for May and June 2018. The first row (a–c) shows the evaluation for CO₂, and the second row (d–f) shows the evaluation for CH₄.

ICOS tall towers are mostly situated outside of urban areas where biogenic and background signals are dominant. In order to further investigate the model performance for GHG tracers, we shift our focus to an area with a stronger influence from anthropogenic emissions in the next section, where we evaluate the model skill using in situ aircraft measurements.

3.4 Evaluation against aircraft measurements

A quantitative statistical comparison between the data from all CoMet 1.0 flights against the different model setups is shown in Fig. 8. The reference scenario (NN_NR) performed the worst among the six scenarios, and no significant performance difference was observed among the other five setups, consistent with the comparison to tower measurements shown in Fig. 7. A detailed overview of the time series for all flights can be found in the Supplement (Figs. S1–S11).

Figure 8 Quantitative statistical metrics for all flights from the CoMet 1.0 campaign for different model scenarios analyzed for CH₄. A more intense color indicates a better score and vice versa. Flights that crossed so close to nearby point sources that we cannot represent them well are denoted with a single asterisk, and flights conducted outside the USCB are denoted with a double asterisk.

Considering the spatiotemporal resolution of the model (5 km, hourly) and the assigned emissions (10 km, hourly), our simulations face limitations in capturing fine features of the plume structure when the flight path is too close to the point source in certain cases. Specifically, flights such as 20180601a and 20180613a involve spiraling around point sources, with the horizontal distance from the source to the spiral loop being <5 km. As a result, our models cannot adequately represent peaks of CH₄ enhancements. Consequently, this led to poor statistical performance across these flights. Flights conducted outside the USCB region, such as 20180607b and 20180614a, fall under a different context for evaluating near-point source emissions. This sub-selection of flights is indicated with asterisks in Fig. 8, with a single asterisk indicating flights with large near-field sources and a double asterisk indicating flights that took place outside the USCB. Flights deemed suitable for model–data comparisons under the influence of a strong near-field source are 20180529a, 20180606a, 20180606b, and 20180611a, as these flights sampled their downwind wall relatively far downstream from individual sources under a well-developed PBL.

Among these selected flights, we found an interesting case where atmospheric transport was significantly improved by grid nudging. The improvement is visible in a comparison of NN_DR and GN_DR against aircraft measurements.

On 11 June (Fig. 9), we observed a noticeable difference of approximately 40 ppb at the upwind leg and the local background between NN_DR and GN_DR, with GN_DR showing a better match with observations on that day. The underlying mechanism driving this improvement in the simulation with grid nudging in addition to daily restarts is discussed in Sect. 4.4.

Figure 9 Flight tracks of the airborne platform colored based on either modeled CH₄ mole fractions (a NN_DR; c GN_DR) or in situ aircraft measurements (b) performed over the Upper Silesian Coal Basin. In situ observations have a high temporal resolution of 1 s spanning about 2.5 h (ca. 9000 time steps), and the model fields were stored hourly. Modeled values were extracted along the flight track from the nearest point in time and space. Also plotted are 2D planes of the modeled CH₄ at a model level at ca. 950 m a.s.l. to give an indication of the plume structure and the location of the point source.

4 Discussion

4.1 Performance in meteorology

The simulation GN_3km_NR performs poorly in comparison with the other simulation scenarios, based on the evaluation of skill in simulating T₂ and Q₂ (Sect. 3.1, Figs. 2 and 3c and d) and the R² of the PBLH (Fig. 4), whereas performance differences in wind speed and wind direction are smaller. The cause of the offset in T₂ and Q₂ relative to observations (see ME in Fig. 2c and d) could be a consequence of the discrepancy in SMOIS, making the atmosphere wetter or dryer (Fig. 10d and e). This may lead to different sensible and latent heat fluxes, which in turn could result in different temperature and humidity close to the surface. This interpretation is supported by significant dependencies of T₂ and Q₂ on SMOIS (Sect. 4.3.2). As there is no nudging below 3 km, T₂ and Q₂ are not being adjusted in simulation GN_3km_NR. Simulation GN_NR also lacks daily re-initializations, just like GN_3km_NR. Nevertheless, GN_NR does not seem to show any notable biases when compared against observations of T₂ and Q₂ (Figs. 2 and 3c and d). The only difference between experiments GN_3km_NR and GN_NR is the threshold altitude of grid nudging. GN_NR assimilates horizontal winds, temperature, and moisture from ERA5 above the model-simulated PBL, which is dynamically diagnosed and follows its diurnal cycle. Therefore, in GN_NR, each night the T and Q profiles are adjusted down to the (low) nocturnal PBL. While this is still above 2 m, T₂ and Q₂ are influenced by this adjustment the next day when turbulent mixing occurs. In this way, offsets in T₂ and Q₂ due to divergent SMOIS are moderately amended each day by grid nudging throughout the simulation period, similar to the effect of daily re-initializations (see Fig. S25 in the Supplement). However, in Sect. 4.3.2, we demonstrate that the divergence in SMOIS, which is still present in simulation GN_NR, has a residual influence on the convective environment, if to a lesser degree than in GN_3km_NR.

Figure 10 Linear dependencies among model quantities related to PBLH performance by comparing simulations GN_DR and NN_NR. Shown are the dependencies on SMOIS of the relative difference in the PBLH (a), the Bowen ratio (c), T₂ (d), and Q₂ (e). The correlation between the PBLH and the Bowen ratio is shown in panel (b). The model outputs were sampled at 12:00 UTC at the 22 radiosonde site locations as in Fig. 4. N_∩ shows the number of samples used in the statistics for each comparison. All the r values shown are significant, with p-values under 0.001. Plots for each location can be found in Figs. S26–S30.

The PBLH simulation in GN_3km_NR did not perform as well as the other nudging and/or restart runs. In Sect. 4.3.2, we demonstrate that the underlying reason is likely soil moisture drift. Simulations NN_DR and GN_3km_DR have almost identical skill in estimating PBLH. Thus, nudging only above 3 km has only a minimal impact in improving the meteorology within the boundary layer, whereas nudging dynamically above the model-simulated PBL shows better performance (seen when comparing GN_DR and GN_3km_DR). Scenario GN_NR performs reasonably well in representing PBLH, with the second-best R² value in Fig. 4 and thus with a small advantage over NN_DR. Finally, with the addition of daily restarts, the GN_DR scenario performs slightly better than GN_NR, having the lowest ME and RMSE and the highest R² among all experiments, albeit by a small margin. In Sect. 4.3.2, we show there is a residual influence of soil moisture drift on the PBLH, which likely explains the small performance advantage of GN_DR over GN_NR. In summary, either resetting the SMOIS regime periodically using ERA5 or mitigating the impact of SMOIS drift via grid nudging leads to a small improvement in the PBLH representation in WRF.

4.2 Performance in simulated greenhouse gases

In Sect. 3.3 and Fig. 7c and f, we can see improvements relative to the reference run for both CO₂ and CH₄ when the model was re-initialized daily and/or when grid nudging was employed; however, we could not clearly distinguish the five simulations that employ either strategy from one another based on these results. This includes the GN_3km_NR simulation, which performed slightly worse than the other nudging/restart simulations in terms of the PBLH. The reason for the similarity of the performances may be that other errors that are common to all simulations dominate over the impact of PBLH differences on the GHG simulation, e.g., transport errors, fluxes, and boundary conditions. This interpretation is supported by Feng et al. (2019a, b), who found that, despite contributions from transport and/or boundary conditions, the uncertainty in modeled atmospheric GHG mole fractions was primarily driven by the underlying fluxes, which in our case are fixed across all six scenarios. The evaluation with regard to the CoMet 1.0 campaign was very similar to the comparison with ICOS tower measurements; i.e., there was not much difference among the simulations. However, we uncovered a scenario that illustrates an improvement in long-range transport achieved through grid nudging. This case is discussed in detail in Sect. 4.4. Comprehensive flight comparisons are available in Figs. S1–S23 in the Supplement.

4.3 Impact of soil moisture drift on WRF performance

Lo et al. (2008) evaluated the skill of grid nudging with a continuous run and a weekly re-initialization run. Their results show that both simulations yield similar performances (seen also in our study (GN_NR vs. GN_DR), e.g., Figs. 2 and 3). They conclude that simulations should not be subdivided from a long simulation into shorter ones because soil parameters generally have a long memory. In the end, Lo et al. (2008) abandoned frequent restarts in favor of a continuous run but with nudging. Furthermore, Vincent and Hahmann (2015) stated that the disadvantage of frequent restarts is wasted computational power for the spin-up period and discontinuities between individual simulations. However, the approach of continuous runs with nudging overlooks the impact on SMOIS. The following sections explain our findings regarding the impacts of SMOIS on model performance, specifically humidity and the PBLH.

4.3.1 Impact of SMOIS on modeled humidity

For instance, both Bullock et al. (2014) and Zittis et al. (2018), who focus on the meteorological performance of WRF, encountered issues with surface-level water vapor being either too wet or too dry. These issues resemble what we have observed in the SMOIS differences between runs with and without restarts for different regions, as seen in Fig. 5. Therefore, SMOIS drift may explain the water vapor discrepancies that Bullock et al. (2014) and Zittis et al. (2018) observed.

Kim et al. (2020) found that, for fog simulation studies, further improvement was achieved when observed soil moisture information was utilized as an initial condition. This suggests that problems in near-surface humidity can be improved by frequent restarts, informed by historical observations from quality-controlled reanalysis fields.

4.3.2 Impact of SMOIS on modeled PBL height

In Sect. 3.1, we showed that both nudging down to the dynamically determined PBL and daily restarts improved the simulated PBLH compared to nudging only above 3 km and not employing daily restarts. We observed the same pattern in the performance of T₂ and Q₂. Here we show that these results can be explained by soil moisture drift in the WRF model and its effective mitigation by daily restarts or, to a slightly lesser extent, the mitigation of its impact on T₂ and Q₂ by grid nudging.

A positive bias in soil moisture increases humidity at the surface and therefore the amount of energy that is being used for evapotranspiration (latent heat). This energy is then not available in the form of sensible heat, reducing temperature. The difference may be quantified by the Bowen ratio, i.e., the ratio of sensible and latent heat flux. Benjamin et al. (2016) demonstrated that such a bias in humidity leads to a positive feedback affecting the development of the planetary boundary layer (PBL), linked with wet or dry bias. For example, large sensible heat flux triggers more turbulent convection, causing a deeper and drier PBL. The reverse is found for wetter regions. Therefore, a positive soil moisture bias results in a negative PBLH bias and vice versa, mediated by a change in the Bowen ratio. The resulting impact on the PBL is clearly not desirable for GHG tracer simulations.

To demonstrate the chain of effects that leads from soil moisture bias to PBLH bias in WRF, we compare the SMOIS, T₂, Q₂, Bowen ratio, and PBLH between simulation setups GN_DR and NN_NR, focusing on locations with radiosonde measurements at 12:00 UTC (Fig. 10). These simulations represent the best and worst PBLH performances, respectively, offering the most pronounced signals for our sensitivity analysis. We analyzed the difference in SMOIS against the relative difference in PBLH in Fig. 10a, where a negative slope (r<0) indicates that wetter soil leads to a lower PBLH. Positive correlations (r>0) between PBLH differences and Bowen ratio discrepancies are evident in (Fig. 10b), suggesting that higher convection corresponds to an increased PBLH. Conversely, Bowen ratio discrepancies negatively correlate with SMOIS divergence (Fig. 10c), implying that wetter soil is associated with reduced convection and vice versa. The influence of SMOIS on T₂ and Q₂ is also apparent (Fig. 10d and e), confirming that the indirect effect does indeed exist in our model when the SMOIS regime is distorted, albeit with a subtle impact (detailed analyses for each location can be found in Figs. S26–S30 in the Supplement). These results are consistent with the findings of Benjamin et al. (2016) and with the mechanism described above regarding how soil moisture differences have an impact on the modeled PBLH.

Consistent with the warm and dry bias reported by Benjamin et al. (2016), a bias in SMOIS will also have an effect on cloud cover and thus on the shortwave radiation reaching the surface. As this is a key driver of the online VPRM module, SMOIS drift may also impact the simulated biogenic CO₂ fluxes (not included in this study).

4.4 A case study demonstrating performance improvement by grid nudging

Our analyses of the CoMet 1.0 flights revealed a case where grid nudging significantly improved the model performance. In Sect. 3.4 we mentioned the notable contrast at the upwind leg between simulations NN_DR and GN_DR in Fig. 9. In simulation GN_DR, an enhancement of CH₄ is seen which is regional, not close to the point source. Figure 11 shows modeled CH₄ over time at the location from Fig. 9 (50.26° N, 18.47° E). We compare the free-tropospheric CH₄ of both simulations at 3500 m a.s.l. with observations when available. We see that the total CH₄ in simulation GN_DR agrees better with the observations than in NN_DR. The main contribution of the errors is not from the background but from transport of anthropogenic CH₄ emitted inside the modeled domain.

Figure 11 CH₄ (a NN_DR; b GN_DR) over time at Silesia (50.26° N, 18.47° E) in mole fraction (ppb). Background methane as CH₄_BCK (dashed lines) and total methane as CH₄_Sum, i.e., the sum of background and anthropogenic CH₄ (solid lines) at 3500 m a.s.l. The white stars show coincident aircraft-based in situ measurements in the free troposphere. Note that observed values are filtered based on model-simulated PBLH to extract data in the free troposphere; hence one data point is omitted for simulation NN_DR June-13-12:00:00 UTC.

We trace back in time and space to find out where the regional offset on 11 June originates and by what means. Figure 12 shows a series of snapshots of differences in simulated CH₄ and PBLH, which show that the creation of the regional offset began roughly 24 h backwards in time in northern Germany. There, the simulations show a disagreement in simulated PBLH, forming a regional enhancement of approximately 40 ppb difference in CH₄ within the atmosphere. This enhancement accumulated through time due to differences in PBLH and was transported southeast to finally reach Silesia at 12:00 UTC on 11 June, when the aircraft measurements were performed. This case demonstrates how simulated PBLH can have a critical impact on simulated GHG mixing ratios.

Figure 12 Difference in the evolution of atmospheric CH₄ at the second model level (top) and boundary layer height (bottom), with and without grid nudging (NN_DR and GN_DR). Columns show snapshots over time, 8-hourly from 12:00 UTC, 10 June, until 12:00 UTC, 11 June. The circled areas in the top row indicate the accumulated regional offset in methane being transported southeast from northern Germany to Upper Silesia. The offset originates from the difference in PBLH at 12:00 UTC on 10 June, also circled in the bottom row.

5 Summary and conclusions

Errors in atmospheric transport often limit the precision and accuracy of long-term modeling of atmospheric tracers, both forward in time and in inversions for estimating GHG sources and sinks. In order to reduce this error component, we have performed a sensitivity study to determine appropriate methodologies for using ERA5 reanalyses from ECMWF to drive high-resolution (5 km horizontally) simulations of WRF-GHG over Europe. Namely, we have focused on using the method of (1) restarting the model daily with fresh initial conditions, to maximize the consistency between WRF-simulated fields with ERA5, and/or (2) FDDA grid nudging throughout the modeled free troposphere. This applies an additional tendency term to the variables that are expected to be critical for transport (wind speed, wind direction, temperature, and moisture) above a given level at each grid cell to gently force the model state closer to that of ERA5. Note that our WRF-GHG experiment does not involve nesting, unlike most previous studies. Using one large domain allowed us to assess performance across a wide range of environments and make use of more data for model evaluation. Furthermore, in contrast to past studies that focused on meteorology, the performance differences in simulating passive atmospheric tracers was also considered.

Six different simulations with different configurations were conducted in order to assess the two main strategies outlined above, alone and in combination (see detailed description in Table 1). We found the following:

1. Applying either daily restarts, nudging above the PBL, or both considerably improved the performance in meteorology and in simulated CO₂ and CH₄ mixing ratios compared to a free run. A small advantage may be achieved by combining both daily restarts and nudging compared to employing only one of the two methods.

2. Without re-initialization and without nudging, the selected land surface model within WRF was unable to properly represent the hydrological cycle over longer simulation periods, causing soil moisture to drift away from observation-driven reanalysis fields, which led to a deterioration in surface temperature, surface humidity, and, to some extent, in the PBLH. Both frequent re-initializations and nudging down to the simulated PBLH alleviated the deterioration of these quantities.

3. Compared to the free-running reference simulation, grid nudging only above a fixed level of 3 km resulted in a considerably smaller performance improvement than dynamically nudging down to the model-simulated PBLH. The reason is that the latter method nudged atmospheric fields throughout nearly the entire vertical column during the nighttime, when PBLHs were low, thereby improving surface temperature and surface moisture in a manner similar to the daily restarts.

4. The modeled PBLH was sensitive to soil moisture drift. The fundamental mechanism is soil moisture's influence on the Bowen ratio, i.e., an impact on the sensible heat flux that drives the development of the PBL. Nudging surface temperature and moisture minimized the impact of SMOIS drift on the PBLH performance.

5. We identified two methods that effectively alleviate the impact of soil moisture drift on the modeled PBLH: restarting WRF daily and nudging down to the simulated PBL. Other methods, which we have not investigated, may be viable as well. These include employing a model that allows the nudging of soil moisture or surface T₂ and Q₂, such as surface nudging (a built-in option within WRF), soil moisture nudging in a PX land surface model (Pleim and Gilliam, 2009), and nudging SMOIS to satellite soil moisture data (Capecchi and Brocca, 2014). However, the former two methods may run into clashes between the two different LSMs (here, Noah for WRF and HTESSEL for ERA5) during the simulation, causing inconsistencies in soil moisture dynamics. The latter requires pre-processing satellite data and a dedicated sensitivity test for this method with our model setup. Hence, we prefer frequent restarts or nudging, as they are easier to implement.

6. Due to the impact of SMOIS drift on the PBLH when nudging only above 3 km without daily restarts (simulation GN_3km_NR) and to the importance of the PBLH for simulating GHG mixing ratios, we expected that daily restarts or nudging down to the dynamically determined PBL would improve GHG performance compared to GN_3km_NR. However, performance in GHG was very similar among those simulations. We conclude that the GHG performance in our simulations was dominated by errors other than the SMOIS drift.

7. SMOIS drift may, via its influence on precipitation, humidity, cloud formation, radiation, and near-surface temperature, disturb biospheric CO₂ fluxes simulated by the online coupled flux model VPRM. However, dedicated sensitivity tests with a longer simulation period would be required to assess this; hence it is outside the scope of this study. Nonetheless, we recommend addressing SMOIS drift in GHG transport model setup. This is successfully avoided by daily restarts to constrain the land–atmosphere exchange and convective environment for longer GHG tracer simulations.

Finally, we would like to note that the frequent re-initialization approach is not only suitable for WRF-GHG: it can also be applied to other mesoscale models simulating tracers, such as COSMO, CHIMERE, and ICON-ART. It ensures that models remain consistent with the reanalysis fields, especially with respect to the land–atmosphere exchange, which is associated with convection and consequently the tracer concentration distribution.

In summary, based on this study, we recommend the combination of grid nudging and frequent re-initialization of the meteorological reanalyses for tracer simulations over using either method alone, and we believe that the additional expense of computational time for spin-ups associated with daily restarts is time well spent.

Appendix A: Quantitative statistical metric scores

Table A1 Quantitative statistical metrics for Sect. 3.1 and 3.3.

Appendix B: Sites used for validation

Table B1 Meteorological sites used in this study.

Table B2 ICOS sites used in this study for evaluating CH₄ and CO₂.