Numerical simulation of atmospheric CO2 concentration and flux over the Korean Peninsula using WRF-VPRM model during Korus-AQ 2016 campaign

We conducted regional scale CO2 simulations using the Weather Research and Forecasting model (WRF) coupled with the Vegetation Photosynthesis and Respiration Model (VPRM). We contrasted simulated concentrations with column, ground and aircraft observations during the Korea-United States Air Quality (KORUS-AQ) 2016 field campaign. Overall, WRF-VPRM slightly underestimates CO2 concentrations at ground and column monitoring sites, but it significantly underestimates at an inland tower measurement site, especially within the stable (nocturnal) boundary layer in nighttime. The model successfully captures the airborne vertical profiles but showed a large offset within the planetary boundary layer (PBL) in the areas surrounding Seoul and around the Taeahn point source emissions in the west coastal area of the Korean Peninsula. A case study flight intended to capture Chinese influence observed no clear signals of long-range transport of CO2, due mainly to the much larger magnitude of background CO2 concentrations. The calculated Net Ecosystem Exchange (NEE) with flux measurements at a tower site in the South Korean Peninsula has also been evaluated comparing with CO2 flux measurements at a flux tower site, resulting in the underestimation by less than a factor of 1.


Introduction
More than half of the global population resides in urban areas, and by 2050, the urban population is expected to increase up to~69% of the world population [1]. As urbanization expands globally, the emissions of anthropogenic greenhouse gases (GHGs) have also increased. About 30%-40% of GHGs are emitted from urban areas, and these anthropogenic emissions are related to severe weather phenomena, such as the increase of average atmospheric PLOS ONE | https://doi.org/10.1371/journal.pone.0228106 January 24, 2020 1 / 21 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 temperature, the increase and/or decrease of precipitation, the higher occurrence, and intensity of serious weather phenomena [2]. Despite its low global warming potential relative to other GHGs, the relative abundance of carbon dioxide (CO 2 ) emissions make understanding its emission sources critical [3]. A wide variety of studies have been conducted to better understand the effects of CO 2 on the carbon cycle between the atmosphere and the Earth's surface, generally using either bottom-up or top-down approaches. A bottom-up approach involves processing information from surveys from individual base anthropogenic inventories, specifying industrial, residential, commercial and transport sectors. These individual inventories are then linked with carbon contents of each fuel type to produce a large, final emission inventory. Errors in this approach are about ±5% uncertainty at the global scale [4] and about 50-200% at the urban scale [5][6]. Upscaling of relations between biospheric CO 2 flux measurements and environmental factors, such as air temperature, radiation, water, and vegetation parameter, also belong to this approach [7,8]. In contrast, a top-down approach is an empirical downscaling method. For example, in inverse modeling, a predictive atmospheric transport model and a priori initial estimates of surface CO 2 fluxes are used to forecast observed concentrations and to estimate the terrestrial CO 2 exchange while maintaining consistency between observed and simulated CO 2 concentrations. The accuracy of this approach is limited by a priori initial estimates of CO 2 fluxes [9,10].
For decades, neighborhood-scale (microscale) flux measurement studies have been carried out to understand CO 2 interactions between the surface and the atmosphere. AmeriFlux and FluxNet have played an important role to quantify the Net Ecosystem Exchange (NEE) of CO 2 [11] in primarily biogenic regions, while efforts to conduct flux measurements in urban areas have been complicated by both anthropogenic emission sources and vegetative uptake [12][13][14]. Observed concentrations are also used to retrieve NEE via large (global) scale inverse modeling [15][16][17], However the lack of spatially explicit a priori flux estimates (i.e. representative errors) can cause large uncertainties for the production of posterior gridded emission data. Between large scale models and flux measurements, the coarse resolution of global inverse modeling makes it challenging to capture the fine resolution of neighborhood flux measurements. To fill the significant scale-gap between the microscale flux measurements and global inverse modeling, finer spatiotemporal resolution modeling is required. Mesoscale (regional) meteorological models can be used to establish the link across this scale gap. To this end, a vegetation photosynthesis and respiration model (VPRM) has been coupled in the Weather Research and Forecast (WRF) model [18] and has been built-in to WRF-Chem (hereafter WRF-VPRM) since version 3.4. This coupled model has been evaluated over various vegetation dominated areas, such as croplands and forests and complex mountainous terrain [19][20][21]. These studies reported significant improvements in capturing the spatiotemporal variation of observed CO 2 concentrations and fluxes not exhibited specifically in global models. Additionally, the coupled model has been evaluated in urban areas where CO 2 measurements are complicated by both man-made emission sources and vegetative contributions. Over the southern California region in the US, a sensitivity study focused on fossil fuel CO 2 emissions revealed that higher resolution input emissions resulted in a better improvement of CO 2 concentration simulation [22]. Another study over the same domain reported significant improvements in NEE simulation after optimization of VPRM biospheric parameters [23].
The rapid rate of development, and the associated increase in carbon emissions, over the last five decades in East Asia, including China, Taiwan, Japan, and Korea, make it an important target for understanding this coupling between global models and local measurements. Based on the Global Carbon Atlas database [24], the East Asia area is identified as the biggest CO 2 emission source: China ranked the 1 st largest emission country by 26.3% of the world total CO 2 emissions in 2016, followed by Japan (3.6%), South-and North-Korea (1.9%), and Taiwan (0.8%). Among fossil-fuel territories, coal (68.8%) was the most dominant emission source, followed by oil (17.0%), cement (8.9%) and gas (5.1%) [25,26]. The transportation and industrial processes primarily using fossil-fuel combustion become dominant anthropogenic CO 2 emission sources. The land-use changes also affect the CO 2 emissions via the interactions between the surface and the atmosphere: the East Asia uptake about-552 MtCO 2 , which is about 20% of the global emissions 2,488 MtCO 2 by land-use changes in 2010 [27,28].
In addition to their own domestic emissions, there are also diplomatic controversies among the countries in East Asia, because air pollutants can affect other countries regionally. Specifically, air pollutants emitted from China can move over the West Sea and have been shown to affect the air quality of other nations downwind, such as Taiwan, Japan, and Korea [29][30][31].
To investigate the factors contributing to poor air quality in Korea, the Korea-United States Air Quality Study (KORUS-AQ) campaign was executed by both Korea's National Institute of Environmental Research (NIER) and the United States National Aeronautics and Space Administration (NASA) in 2016 spring. The KORUS-AQ study conducted extensive observations of trace gases and aerosols in order to identify specific emission sources, and detailed information can be found at the KORUS-AQ white paper [32].
Among a wide variety of trace gases, in which most of the air pollutants are involved, we focused on CO 2 in this study. We calculated WRF-VPRM over the Korean Peninsula during the KORUS-AQ 2016 campaign period and evaluated the modeling performance comparing with the concentration and flux observations as well as the meteorology.

WRF-VPRM model
In this study, we used VPRM coupled in WRF-Chem (version 3.8.1); more detail can be found in [18]. We utilized two-way nesting with grids at 10 km and 3.3 km resolution (Fig 1) and 38 vertical layers (with 12 layers below 1.5 km) extending up to 100 hPa. Initial conditions (ICs) and boundary conditions (BCs) for meteorological fields for the WRF modeling were fed from the 3-hr GDAS/FNL 0.25 degree global tropospheric analyses and forecast grids data of National Centers for Environmental Prediction [33] and the 6-hr National Centers for Environmental Prediction (NCEP) SST dataset with 8-km horizontal resolution data [34]. The WRF Single Moment (WSM) 3-class microphysics scheme, New-Grell cumulus scheme (only for a coarse domain), and Rapid Radiative Transfer Model (RRTMG) short-and long-wave radiation scheme were used. The Yonsei University scheme [35] with the topo-wind and pblmix options were used for the Planetary Boundary Layer (PBL), which we determined based on sensitivity tests (S1 Table). The WRF four-dimensional data assimilation (FDDA) was applied using the grid nudging options. We conducted the numerical simulation during two periods: (1) from April 26 00:00 UTC to May 11 00:00 UTC, and (2) from May 12 00:00 UTC to June 10 00:00 UTC, with 5 days spin-up time for each simulation.
The Enhanced Vegetation Index (EVI) and Land Surface Water Index (LSWI) from the Moderate Resolution Imaging Spectroradiometer (MODIS) surface reflectance data are used to generate the scaled temperature (T scale ), physiological index (P scale ), and canopy water content (W scale ). These scaling factors and the VPRM parameters, including the maximum quantum yield (λ) and the half-saturation value of photosynthetically active radiation (PAR 0 ), are used to calculate the Gross Ecosystem Exchange (GEE). The simulated shortwave radiation (SW) is used for the photosynthetically active radiation (PAR), as SW is strongly correlated with PAR: SW � 0.505 × PAR (units: SW, W m -2 ; PAR, μmol m -2 s -1 ) [18]. The respiration rate (RESP) is calculated from the WRF's 2-m air temperature (T) along with first order linear parameters, including a slope (α), and an intercept (β). The summarized equations for the calculation of NEE (GEE + RESP) used in the VPRM module are below: The VPRM parameters should be optimized based on local CO 2 flux measurements for each representative land-use class. However, due to the lack of observations over the study domain during the study period, we used the default VPRM parameters.

Fossil fuel CO 2 emissions
In this study, we used the NOAA 2016 CarbonTracker 3D CO 2 mole fraction data (version CT2017 for CO 2 lateral boundary and initial conditions. The 3-hr CT2017 data, containing global background, photosynthesis/respiration by biosphere, fires, combustion of fossil fuels, and air-sea exchange, covers all globe up to the tropopause, with 2.5˚× 2.5˚spatial resolution [37]. To feed surface boundary conditions for CO 2 , hourly fossil fuel CO 2 emission data for the year 2015 from the Fossil-Fuel Data Assimilation System (FFDAS) [38], gridded at 10 km, were used. The FFDAS is a data assimilation system that estimates the fossil fuel CO 2 emissions at every grid cell by solving the Kaya identity, which expresses emissions as a product of population density, per capita economic activity, energy intensity of the economy, and carbon intensity of energy [39]. A series of spatially explicit observation data are used to constrain each factor in the Kaya identity. The FFDAS uses remote sensing-based nighttime lights data, gridded population sector-based fossil fuel CO 2 emissions from the International Energy Agency (IEA), and global power plant CO 2 emissions [6]. Natively, the FFDAS estimates fossil fuel CO 2 emissions at 0.1˚and annual resolutions over the globe. Nasaar et al. [40] derived scale factors that can apply to global emission data to represent the diurnal and weekly cycles of fossil fuel CO 2 emission variation at 0.25˚× 0.25˚resolution, which we used for this study.

Observation data
In order to better understand the factors contributing to poor air quality in Korea, the KORU-S-AQ campaign was performed by both Korea's National Institute of Environmental Research (NIER) and the United States National Aeronautics and Space Administration (NASA) from May to early June 2016. The KORUS-AQ study conducted extensively spatial and vertical observations of both trace gases and aerosols using ground sites, aircraft, and vessels in order to identify specific emission sources and support policymakers to establish air quality mitigation strategies. Detailed information is available on NASA's KORUS-AQ webpage [41].
The NASA DC-8 aircraft was instrumented to measure vertical profiles of meteorological variables and chemical species. Ambient CO 2 dry molar ratios were measured using a modified commercial non-dispersive infrared (NDIR) spectrometer (LI-COR 6252) with a combined precision and accuracy of 0.1 ppm (1σ) at a 1 Hz sampling rate. In addition, we used the 92 Automated Synoptic Observation System (ASOS) data and 6 rawinsonde data of the Korean Meteorological Administration (KMA) to evaluate the meteorological modeling performance. In order to estimate simulated column concentrations, we used the Total Carbon Column Observing Network (TCCON) data [42].
We also evaluated ground level CO 2 concentrations using a variety of measurements at different sites, summarized in Table 1. These include the World Meteorological Organization Global Atmosphere Watch (WMO GAW) data [43] at not only Anmyeondo (AN), Jeju Gosan (GO) and Uleungdo (UL) but also Boseong (BO). BO is a tall tower site which can measure vertical profiles of CO 2 at 60 m, 140 m and 300 m a.g.l., but only 60-m and 300-m data were available during the study period. The instrumentation is located in the middle of cropland, mainly consisting of rice paddies, that opens to a bay area to the south and is surrounded by hills in the other directions~2-3 km from the tower. Small residential houses and farming warehouses are located at the mountain base.

Results and discussion
The simulated results for the finest domain where, M i and O i indicate modeling results and observations, respectively, at each grid point (i). The time used in this manuscript is South Korean Local Standard Time (LST) unless otherwise specified.

Meteorological evaluation
We compared the simulated meteorological results from WRF-VPRM with 93 ground measurements at ASOS sites run by KMA. The observed, averaged 2-m temperature, which is an important meteorological factor for the calculation of vegetation respiration in the VPRM module, was 18.5˚C with the model underestimating by 1.0˚C (IOA = 0.94). The averaged 10-m wind speed was 2.1 m s -1 with the model underestimating by 0.2 m s -1 (IOA = 0.81). The 10-m wind directions ranged from 190-280˚, and the dominant wind direction was SW and W during the study period, except May 10 and 11 (NE wind direction). The simulated offset in wind direction was < 15˚in average. The vertical profiles of simulated meteorology were also evaluated with 6 rawinsonde observations [45]. The model underestimated the potential temperature by 0.1˚C and wind speeds by 0.4 m s -1 . The simulated PBL height (PBLH) is also critical, because trace gas concentrations are ultimately determined by the volume of PBL. PBLH was calculated using the bulk Richardson number [46]. WRF-VPRM successfully captured the diurnal variation in PBLH, which began to increase at~06:00 LST, reached a peak at~15:00 LST, and then decreased until2 0:00 LST. The model underestimated by~73 m in average.

Comparison with ground observations
Horizontal spatiotemporal distribution. Spatiotemporal variations in atmospheric CO 2 concentrations are typically associated with the spatial distribution of emission sources and biosphere combined with the temporal variations in PBLH. Fig 2 shows the horizontal distribution of land cover retrieved from SYNMAP and emissions from FFDAS in the study domain. The significant point sources are mainly located in the west side of mid-Korea, mostly in and near the Taeahn peninsula, in which Yeongheung and Boryeong power plants and Daesan and Dangjin industrial complexes exist. Gwangyang and Yeosu industrial areas are also identified as large point sources in the south. Among several major cities of Korea, Seoul Metropolitan Area and its vicinity (SMA, a capital area) showed the largest area emissions, followed by Ulsan (Korean largest petrochemical industrial city), and Busan (Korean largest harbor city). The ranges of emissions are reported by color code in S1 Fig.
The highest CO 2 concentrations appeared in the early morning (06:00-09:00 LST), associated with the lowest PBLH and the morning rush-hour vehicle emissions (not shown). Around noon, when the PBL fully developed and vegetation uptake thoroughly increased, the CO 2 concentration remained low until late afternoon. After that, the concentration began to gradually increase again by the start of evening rush-hour while PBLH decreased. The increased concentration remained high through the night until early next morning.
Column and ground concentrations. Before the further comparison of the simulated results with the ground and vertical observations, it is necessary to confirm the total mass consistency in the modeling domain. A total column concentration of CO 2 simulated can be compared with the observed column CO 2 concentration (xCO 2 ) data acquired from the TCCON site at AN. The averaged xCO 2 showed 405.4±0.7 ppm, and the vertical accumulated concentration of WRF-VPRM was underestimated by approximately 0.3 ppm, which is comparable to the resolution of TCCON (0.25 ppm).
At the Korean CO 2 background monitoring sites at AN, GO and UL observed mean concentrations of 414.1±2.9 ppm, 411.9±4.4 ppm, and 410.5±3.9 ppm, respectively, during the study period. Assuming that the UL mean is the background concentration of our study domain (Korea), it was about 2-7 ppm higher than the global CO 2 background measured at Mauna Loa (408±1 ppm) in May 2016 [47]. At each CO 2 monitoring site, WRF-VPRM underestimated the CO 2 concentration at AN by 2.4 ppm (RMSE = 3.6) and at GO by 0.4 ppm (RMSE = 4.7) and overestimated at UL by 0.7 ppm (RMSE = 4.6), on average. The larger MB at AN is likely from the relatively more sparse observations (only four days) compared with the other two sites.
Besides the uncertainties from simulated meteorology, the impact of the temporal scaling factors in the FFDAS emission data, which enables annual global emission data to have a diurnal and weekly cycle, can introduce uncertainty in our simulations. Based on sensitivity tests using GEOS-Chem [40], the impact of the scaling factor on the surface atmospheric CO 2 was negligible for areas far from sources and ranged approximately from 1.5 to 8 ppm over large urban areas. Its effect on the xCO 2 was as high as 0.1-0.5 ppm over the major urban areas. Judging on the locations of the three CO 2 monitoring sites, which are set up among the areas far-from and nearby sources, our simulation MB is tolerable.
The model captured the clear diurnal variations: dominant peaks appeared at 06:00-07:00 LST due to the CO 2 emitted on the early morning and/or carried-over from the previous day when PBLHs were still low. In spring, as the vegetation starts to grow, it is expected that the CO 2 concentration will gradually decrease as biospheric uptake increases. The variation of daily minimum CO 2 observations, which can minimize the influence of erratic anthropogenic CO 2 , can be a good tool to check the temporal trend. Fig 3 shows the 5-day moving average of observed and simulated daily minimum CO 2 concentrations, indicating that our model successfully captured the CO 2 decreasing trend as summer comes.
In addition to evaluation using ground CO 2 monitoring sites, we estimated the modeling performance with a tall tower observation at BO (Fig 1) at both 60 m and 300 m a.g.l. Note that vertical CO 2 concentrations provided by tower sites were only available at BO during the study period. Fig 4 shows the averaged diurnal variation of the observed and simulated CO 2 concentrations, temperature and winds at each level. The observations at 60 m and 300 m showed a peak at~06:00 LST and~09:00 LST, respectively, then both concentrations gradually decreased until~18:00 LST and then increased through the night. The concentrations at the lower level showed 3.8 ppm and 9.7 ppm higher than that at the upper level in day-and nighttime, respectively. The larger difference at night can be interpreted at the point of view of micrometeorology. At night, the 60 m level was within the stable (nocturnal) boundary layer, while the 300 m level was located in the residual layer. The nighttime inversion was observed in temperature as well, and the averaged PBLH calculated by WRF was~100m under the assumption that we can treat the nighttime PBLH as the top of the stable boundary layer. Interestingly, the inversion appeared unique to the study period, as it was not observed in the other spring (S2 Fig). During mid-day, in contrast, there was much less difference in concentrations (~2 ppm) because the two levels were together within the mixed layer. Overall, the observations showed a typical diurnal variation of CO 2 concentration within PBL.
Our model successfully captured the diurnal variation of CO 2 concentrations at 300m, and the simulated results underestimated by approximately 3.9 ppm (RMSE = 9.6) and 0.5 ppm (RMSE = 18.8) in day-and nighttime, respectively. At 60 m, on the other hand, the model failed to capture the amplitude of temporal variations and significantly underestimated nighttime concentrations by 10.3 ppm (RMSE = 14.3). The uncertainties of WRF-VPRM can be caused by vertical mixing, advection, initial fields, emissions and VPRM parameters, described in previous WRF-VPRM studies [19,21]. Considering the model's reliable simulation performance in micrometeorology, as previously described, the significant offset of CO 2 concentration at 60 m was likely due to the uncertainty of emissions around BO. The prevailing winds were from the south, where the local emission sources, such as fishery boats, vehicles, small factories and residential houses in the bay area, are located. These small emission sources are not mandatory to report, so they are not likely fully represented in the relatively coarser resolution of emission inventories. In sum, due to the underestimated nighttime emissions at each grid in model, the model failed to capture the accumulated concentrations within the stable (nocturnal) boundary layer, at least around BO during the study period.
We also evaluated the vertical gradient in CO 2 concentrations, using the daily minimum values at two levels, resulting in that the model underestimated by about a factor of 2, due to the significant underestimation of CO 2 concentration as described above. It is not clear whether this relatively poorer performance can be also shown in other areas, due to the sparser tower observations. Instead, aircraft observations are helpful for estimating the vertical gradient at higher levels.

Comparison with aircraft measurements
In this section, we evaluated WRF-VPRM's performance with the NASA DC-8 aircraft observations. The 4-D fields of simulated CO 2 concentrations were collocated with flight tracks in space and time at each hourly simulation output. The averaged offset of pressure between the simulated results and the DC-8 observations was~10±15 hPa in total 18-flight days. In each day of observations, DC-8 took off and flew along planned routes (Fig 1). For the clarity, we divided the aircraft observations into four groups, based on the land cover below the flight tracks and the type of emission sources: SMA, the West (Yellow) Sea, and two main domestic airways in Korea. The SMA represents air quality over the Korean capital city, in which anthropogenic combustion emissions are dominant. The flight tracks over the West (Yellow) Sea were intended to detect the long-range transport of air pollutants from China, where the DC-8 took meridional tracks over the ocean forming a "wall" between China and the Korean Peninsula. The two domestic airways, connecting Seoul to Jeju (J-route) and Busan (B-route), which are the most frequent observation routes that DC-8 took during the study period. The J- Numerical simulation of atmospheric CO 2 concentration and flux during Korus-AQ 2016 campaign route was designed to sample the air over local point sources, transported air from the West Sea, and croplands. In contrast, the B-route was designed to capture the biospheric activity mostly over mountains and the Busan urban emission.  (Fig 6 and Fig 7). The full set of time series for 18 individual flights is provided in S3 Fig and S4 Fig. During the study period, the aircraft flew over various land-cover and land-use types, maneuvering in and out of the PBL. Overall, as expected, the lower concentrations appeared at higher altitudes and/or over the ocean, while the higher concentrations emerged over urban areas and/or within the PBL. In the vertical profiles in Fig 7, observed CO 2 concentrations exhibited the maximum concentrations at < 1 km and gradually decreased along with the height up to 7-10 km, showing the minima at about 2-3 km. Overall, the averaged concentration observed was 410.2±5.2 ppm across all flights. WRF-VPRM underestimated concentrations by 1.7 ppm (RMSE = 6.9, IOA = 0.6), and the averaged CO 2 concentrations above the PBL ranged from 407.1 ppm to 409.7 ppm, while those within the PBL showed the largest biases mainly from SMA.
The Seoul metropolitan area (SMA). About 55% of the measurements were taken below 1.5 km by the DC-8 when it flew over the SMA. Within the PBL, this group showed the highest averaged concentration of 420.6 ppm compared with other groups. This high value is most likely to be driven by Seoul populated local emissions. WRF-VPRM captured the CO 2 vertical profile over the SMA like other groups, except near the ground (Fig 7). The simulated results underestimate concentrations by 7.3 ppm and 0.4 ppm within and above the PBL, respectively. The large bias in the PBL concentrations may be caused by the simulated deficiency from both/either the vertical gradient CO 2 near the ground and/or the incompletely established emissions over the SMA. We roughly assumed that the vertical gradient of CO 2 concentration within the PBL can be calculated by the difference in concentration between the ground and the 2-km level. The model underestimated the vertical gradient by~30%. This implies that densely populated emissions in Seoul are not fully resolved in our model and it causes the underestimation. Tang et al. [48] also reported a similar poor result in the same study domain and period, using a different model, the Copernicus Atmosphere Monitoring Service (CAMS), and pointed out the limit of the model's resolution.
The West (Yellow) Sea. The effects on CO 2 concentrations from long-range advection from the outside of the Korean Peninsula should be also considered. Legs of DC-8 flight tracks on May 4, 18, 25, 30, and 31 (LST) were designed to detect outflow plumes from China, flying along the "wall" tracks over the West Sea. We displayed the simulated results over East Asia on May 25, for example, in Fig 8. The polluted air plumes starting from highly populated cities, such as Shanghai, Tianjin, and Hangzhou in China, flew in westerlies and arrived to the "wall" within a day. The averaged observed CO 2 concentration was 410.9 ppm and was underestimated by WRF-VPRM by~2.3 ppm (RMSE = 5.5) for the five days of flight tracks along the wall. In order to quantify the magnitude of China's outflow, we conducted additional numerical experiments in which we used the same simulation configurations, described in Sect. 2, but enabling only China's FFDAS emissions in order to to detect the simulated advection of CO 2 originating from China. It was expected that the concentration difference between the two cases under westerly wind conditions would help quantify the effects of long-range CO 2 transport at the "wall". However, the averaged difference was~0.2 ppm (RMSE = 0.6; IOA = 0.99), which is too low to clearly quantify the effects of China's outflow compared with the atmospheric background CO 2 concentration. It also agrees to the results of Tang et al. [48].

Seoul-Jeju and Seoul-Busan airways.
Among several domestic commercial airways in Korea, two routes (J-and B-route) were chosen for the investigation by DC-8 aircraft observations. The J-route is close to the western coastal areas of the Korean Peninsula and has more anthropogenic emission sources and croplands, while the B-route is located over more forested inland areas. In May 2016, the J-and B-route had 8,176 and 1,646 flights in total, respectively, including both passenger and freight flights, based on statistics from Korea Airports Corporation [49]. However, the effects of this number of domestic aircraft on CO 2 concentrations is minor compared with the surface anthropogenic emissions below the flight tracks. In the South Korean Peninsula, there are a lot of huge point sources, such as fossil-fuel power plants and industrial areas. About 30% of the total 39 fossil-fuel power plants of Korea in 2016 were located in the western part of Korea, and they generated energy of~16 GW out of~32 GW of all Korean fossil-fuel power plants, reported by the electric power statistics information system [50]. Overall, the averaged concentrations on J-and B-route were 410.3 ppm and 408.1 ppm, respectively.
The effects of the point sources on the simulated results were also very apparent. On June 2 and 3 (Fig 5), for example, the flight tracks were designed to cover both J-and B-routes. DC-8 detected hot spots (up to~450 ppm) for two consecutive days, when the aircraft made a slight detour from the direct J-route to the northern tip of the Taeahn peninsula. The uncertainty caused by the point sources can be clearer by removing these hot spots in DC-8 observations. After removing the peaks, the MB in the simulated results was significantly reduced by a factor of~5. On June 5 (Fig 5), the flight track was specifically designed for a more intense survey of point sources around the Taeahn peninsula. Comparing with other aircraft observations, the simulated results failed to capture these huge peaks and significantly underestimated CO 2 by 4.6 ppm (RMSE = 11.5). The poorer statistics are due mainly to the scattered point sources that are not fully resolved and/or that are blended in coarser grids compared with the spatial scale of aircraft observations. That is, while the aircraft observed air pollutants at the subgrid scale, the model was not able to fully detect them within its lower resolution (3.3 km). Regional scale inverse modeling can help us better estimate the locations of these local sources using aircraft-based prior observations, as shown by Brioude et al. [51]. However, this is outside of the scope of this study.

CO 2 fluxes
Although CO 2 is chemically inert in the atmosphere, it continuously interacts with the biosphere, so CO 2 flux simulation is necessary to estimate atmospheric CO 2 . The interactions can be represented by NEE. VPRM provides the NEE as the sum of GEE and RESP (Sect. 2.2). The simulated NEE can be evaluated by comparing with measured CO 2 fluxes for each vegetation type. Under the assumption of that the simulated NEE equals to the measured CO 2 flux as suggested by Park et al. [23], we compared the observed flux data at the ground level with the NEE at BO during the study period (Fig 9). The simulated results showed the clear diurnal effects of vegetation: the daytime negative values due to photosynthesis and the nighttime positive values driven by biospheric respiration. The time series of NEE showed increased intensity in summer, which is inversely correlated with decreasing of CO 2 concentrations (Sect. 4.2). Median values of observed fluxes were -0.20 μmol m -2 s -1 and 0.84 μmol m -2 s -1 , and our model showed approximately +98% and -52% offset in day-and night-time, respectively.  The uncertainty in the simulated NEE can be influenced by several factors. First, errors in simulated 2-m air temperature can generate biased NEE, as pointed out in previous studies [18,20]. Our validation statistics showed the simulated 2-m temperature was quite close to observations, so the meteorological uncertainty can likely be neglected. Second, direct comparison of simulated NEEs with measured CO 2 fluxes can also cause large uncertainties, because the flux measurement naturally includes a CO 2 storage term within the canopy [52]. The bigger the storage term of CO 2 fluxes, the larger offset in NEEs. The canopy within the fetch of the flux site at BO dominantly consists of short rice fields, which would minimize the uncertainty from the storage term bias compared with tall and large-bodied tree canopies. Third, VPRM parameters can add uncertainty to the NEE. As described in Sect. 3, the VPRM parameters represent vegetation characteristics, such as light efficiency, photosynthetic active rate, and respiration, which should be optimized based on the flux measurement data. A few studies applying this optimization have been conducted: Jamroensan [20] improved CO 2 modeling performance by adding a new category of vegetation such as soybean over the Midwestern US; Park et al. [23] also improved the NEE simulation results, optimizing VPRM parameters over the Southern California basin.
The VPRM parameters should ideally be optimized at flux measurement sites in various vegetation types. We could not conduct this optimization due to the lack of flux measurement data during this study domain and period. In another study period, where the abundance of long-term (> 1yr) flux observations is available, future studies could suggest improved parameters for the Korean Peninsula.

Conclusions
Over a study domain, the South Korean Peninsula, we conducted numerical simulations of atmospheric CO 2 , using a coupled WRF-VPRM model. The model's performance was evaluated by comparing with observations acquired from ground CO 2 monitoring sites, a TCCON spectrometer, a CO 2 flux site, and NASA DC-8 aircraft observations during the KORUS-AQ 2016 field campaign. WRF-VPRM successfully captured the spatiotemporal variation of CO 2 concentrations over the South Korean Peninsula, strongly associated with the temporal variations of PBL and the spatial distribution of emission sources and biosphere. The simulated results were comparable with the ground CO 2 concentrations observed at two CO 2 monitoring sites, GO (MB~-0.4 ppm) and UL (MB~-0.7 ppm), and with the CO 2 column concentration from the TCCON spectrometer at AN (MB~-0.3 ppm). Beside the monitoring sites, we also compared simulations with tower observations at BO. The model significantly underestimated the concentration at 60 m at night and the vertical gradient of CO 2 Besides the monitoring sites, the model's performance was also evaluated at an inland tower site at BO. WRF-VPRM underestimated the CO 2 concentration by~3 ppm at 300 m and~9 ppm at 60 m. The largest bias occurred in nighttime at 60 m. During the study period, there was a nighttime inversion layer, the height of which was located between 60 m and 300 m. Considering the model's successful simulated results of the nighttime stable (nocturnal) boundary layer, the main culprit of the significant underestimation likely resides with the fact that the relatively coarser emission data FFDAS could not fully resolve the local sources around BO (Sect. 3.2).
At higher altitudes, the simulated results were evaluated with the NASA DC-8 aircraft observations, separated into four flight groups: SMA, the West Sea, and the J-and B-route. The model captured variations in CO 2 concentrations along the DC-8 flight tracks, depending on horizontal and vertical locations. However, it failed to capture large concentrations spikes appearing in the PBL, especially over the SMA where Seoul urban emissions are dominant, and over the Taeahn Peninsula, where power plants and industrial complexes are densely located. Overall, the western part of Korea contains the dominant emission sources in the study domain.
The influence of long-range transport of CO 2 from China was also investigated on May 25, when the DC-8 aircraft flew along the "wall" flight track over the West Sea. The model captured the time series of CO 2 concentrations. To investigate the effect of plume transported from China, we conducted additional simulations in which only Chinese emissions were established. Unexpectedly, the magnitude of concentrations transported from China was very low Numerical simulation of atmospheric CO 2 concentration and flux during Korus-AQ 2016 campaign compared to the background CO 2 concentration. Therefore, we were not able to quantify any significant CO 2 advection from China from this particular case.
In addition to the evaluation of CO 2 concentrations, we also compared the simulated NEE with measured CO 2 fluxes at a flux tower to evaluate the modeling performance for the CO 2 exchange between the biosphere and the atmosphere. Although the model captured the temporal variability in the fluxes: negative and positive values during day-and night-time, respectively, it showed a significant offset compared with the observed fluxes. VPRM parameters are variable season by season and region by region, so they should be optimized at each vegetation type. Further optimization of VPRM parameters could lead to better performance, as suggested in previous studies, but was not possible for this study due to the lack of available sufficient flux data at various vegetation types over a long time period (> 1 yr).
The output of this modeling work can be used to validate remote sensing instruments, such as GOSAT and OCO-2. The results are also useful to support environmental mitigation policy of the Ministry of Environment, especially for the fossil-fuel carbon emission controls, and/or biospheric management in Korean Government. Visualization: Soon-Young Park.

Supporting information
Writing -original draft: Changhyoun Park.