Biogenic CO₂ flux uncertainty: numerical experiments and validation over south-eastern South America

2.1. Model simulations and description of the cases

The WRF-Chem v.4.4.2 model was used to simulate the atmospheric dynamics over the MPRArg region with two nested domains. The resolutions of the outer and inner domains were 12 and 4 km in the horizontal direction respectively; thus, no deep moist convection scheme was required for the inner domain (Fig. 1). Eight different combinations of WRF ABL, cloud microphysics and radiative parameterisations were run for three cases with different meteorological conditions. The reanalysis from the European Centre for Medium-Range Weather Forecasts (ERA5; Hersbach et al. 2023a, 2023b) available every 3 h was used as input and boundary conditions, and the WRF outputs were saved every hour. The most relevant parameters and configurations of the WRF simulations are shown in Table 1.

Fig. 1.

Outer (a) and inner (b, c) domains of the WRF simulations. The symbols indicate the location of the surface synoptic (crosses and circles) and upper-level (diamonds) weather stations used to validate the simulations. The greyscale in (a, b) indicates altitude and in (c) the dominant vegetation class. Measurements include 2-m air temperature (T2), shortwave incoming radiation at the surface (SWDOWN) and precipitaion (PP).

Each WRF simulation lasted for 30 h and was initialised at 00:00 hours UTC (Coordinated Universal Time) on 19 November, 24 November and 17 November 2018 (local time GMT–3); these cases were named sunny, cloudy and rainy respectively, owing to their dominant synoptic situations. Fig. 2 shows the 2-m air temperature (T2) diel cycles at five surface weather stations, 925-hPa geopotential height from ERA5 and GOES-16 (Geostationary Operational Environmental Satellite) albedo from Band 2 centred at 0.64 µm for the three selected events. The sunny case had mainly clear skies over the whole domain, except for a small cloudy area that lasted a few hours in the south-western part of the domain. By contrast, the domain was mostly covered by thin, high-level clouds during the cloudy case, but no significant precipitation was recorded during this event. The rainy case was associated with a cold front moving from south-west to north-east across the domain. This cold front triggered organised deep moist convection in the form of a mesoscale convective system at the north-east corner of the study region during the local afternoon, and the station-averaged recorded precipitation was 12.2 mm day⁻¹ and the maximum was 57.0 mm day⁻¹. These three cases represent increasing levels of complexity and presumably increasing uncertainty associated with the computation of surface CO₂ fluxes.

Fig. 2.

T2 station data during the simulated cases for the (a) sunny, (d) cloudy and (g) rainy cases in Gualeguaychú, Rosario, Marcos Juarez, Río Cuarto and Reconquista stations (solid, dashed, dot-dashed, densely dotted and loosely dotted lines respectively). Their geographical coordinates ae −33.00, −32.92, −32.70, 33.12 and −29.18°S; −58.62, −60.78, −62.15, 64.23 and −59.70°W respectively. ERA5 geopotential height at 925 hPa in geopotential metres (mgp) from 16:00 hours UTC (13:00 hours local time) for the (b) sunny, (e) cloudy and (h) rainy cases. GOES-16 albedo from Band 2 centred at 0.64 µm at 16:00 hours UTC for the (c) sunny, (f) cloudy and (i) rainy cases.

To evaluate the sensitivity of CO₂ flux computation, we performed simulations over the three selected cases using different combinations of the parametrisation schemes of short-wave radiation, ABL, microphysics and soil model. The Control experiment of WRF was designed by selecting the same set of parameterisations used by the National Meteorological Service of Argentina (Servicio Meteorológico Nacional, SMN) for their 2-day forecasts (Dillon et al. 2020). To test the sensitivity of the input variables to the model configuration, different runs were carried out changing only one of the schemes used in the Control setting (Table 2). Currently, over 70 physical parameterisations are available in WRF. A full analysis of all these would increase the complexity of the study beyond its scope. Therefore, we limited the number of WRF scenarios to eight. The selected schemes have been used widely by other researchers in the region with positive results (e.g. Ruiz et al. 2010; Ferreyra et al. 2016; Müller et al. 2016; Dillon et al. 2020, 2021; Alvarez Imaz et al. 2021; Casaretto et al. 2022; Suárez et al. 2022; Luque et al. 2024; Merino 2024). A brief description of the eight WRF configurations is detailed in the following subsections.

2.2. Validation data

Conventional weather stations from the SMN and automatic weather stations from the Instituto Nacional de Tecnología Agropecuaria (INTA) inside the MPRArg were used to validate WRF-simulated T2, incoming shortwave radiation at the surface (SWDOWN) and precipitation (PP; Fig. 1). Among all available stations within the model domain, a subset was selected seeking a more homogeneous spatial distribution (i.e. reducing the station density in highly populated areas). The characteristics and locations of the surface weather stations used in the validation of the simulations are presented in the Supplementary material (Supplementary Tables S1–S3).

The selected dates allowed us to validate the results with measurements from RELAMPAGO-CACTI (Nesbitt et al. 2021; Varble et al. 2021) and Mar Chiquita (Bautista et al. 2023) micrometeorological field campaigns. The performance of the WRF simulations in representing ABL structure was validated using atmospheric soundings from the RELAMPAGO-CACTI campaign. The dataset consisted of 5, 5 and 13 sounding profiles respectively launched during the sunny, cloudy and rainy cases at Villa María del Río Seco, Córdoba, Ezeiza and Santa Rosa (Fig. 1). All radiosonde launching times were matched to the closest WRF output hour.

WRF-simulated SWDOWN was validated using data from the INTA automatic surface weather station network. This dataset was complemented with extraordinary observations from the flux station measurements from RELAMPAGO-CACTI (Nesbitt et al. 2021; National Center for Atmospheric Research Earth Observing Laboratory In situ Sensing Facility and Oncley 2021; Varble et al. 2021) and from the Mar Chiquita flux station (Bautista et al. 2023). It is important to highlight that, even with the extended database, the stations measuring radiation were limited to the centre and south-east regions of the inner domain; thus, SWDOWN validation was restricted to that area.

In addition to in situ measurements, hourly outputs from ERA5 were used to compare T2, SWDOWN, PP and ABL structure with the WRF simulations and with the surface and sounding measurements.

2.3 Processing and validation

The model was evaluated using the bias, Root Mean Squared Error (RMSE) and Pearson correlation index (ρ), which are defined as follows:

where X_modi and X_measi are the respective modelled and observed variables, and N is the number of available observations.

The indices from Eqn 1–3 were used to quantify and characterise the uncertainties of the model. In particular, non-zero biases are systematic errors that lead to over- and under-estimations by the model. However, this index can be zero even when there are large discrepancies between the model and the observations, as long as the mean is accurately estimated. In those cases, the RMSE comes in handy because large errors result in a high RMSE. Nevertheless, neither bias nor RMSE gives insights into the diel cycle of the variables. This aspect can be elucidated by examining ρ. A ρ close to 1 or −1 signifies that the two variables are respectively nearly in phase or counterphase. Conversely, a ρ value close to zero suggests that they share neither a phase nor a counterphase relationship.

The observation data were interpolated in time to the WRF output time using a nearest neighbour approach. This applies to INTA, RELAMPAGO-CACTI and Mar Chiquita datasets, whose sampling frequency is of the order of a few minutes. The nearest-neighbour approach was used to interpolate WRF model data to the location of the weather stations. In this way, the grid point closest to the observation location was selected in the WRF and ERA5 domains. The only exception was WRF PP, where the four grid points closest to the station’s location were averaged to smooth out the small-scale spatial variability from the model output. However, radiosonde data were compared with simulated data at the ERA5 vertical levels. WRF and radiosonde data were interpolated to ERA5 vertical levels using the same nearest-neighbour approach that was used to interpolate the horizontal grid.

2.4. CO₂ flux uncertainties

We calculated uncertainties in the CO₂ flux estimations for two R_eco models and one GPP model: the R_eco and GPP models developed by Mahadevan et al. (2008) and already implemented in VPRM-WRF (R_ecoM and GPP_M), and a modified version of R_ecoM model developed by Gourdji et al. (2022) (R_ecoG). They are represented by Eqn 4–7:

In Eqn 4–7, where T is the near-surface temperature (in this study, we use T2), β, α, λ, PAR0, T_min, T_max and T_opt are parameters that represent the baseline ecosystem respiration, the temperature sensitivity of ecosystem respiration, the maximum quantum yield, the half-saturation value of photosynthetic active radiation, and the minimum, maximum and optimal temperatures for photosynthesis respectively, whose value is defined in the VPRM-Europe table (Ahmadov et al. 2009). All these parameters depend on the type of vegetation. However, β₁, α₁, α₂, γ, k₁, k₂ and k₃ are parameters that quantify the influence of each term on R_eco whose values are taken from Gourdji et al. (2022); EVI, P_scale, W_scale and T_scale are parameters ranging from 0 to 1 respectively representing greenness, leaf age, water stress and influence of T2 on vegetation photosynthesis. Absolute errors for all the models can be calculated after expanding Eqn 4–7 by the Gaussian Error Propagation method, resulting in the following error expressions:

Eqn 8–11 represent the absolute errors of R_ecoM due to T2 (Eqn 8), GPP_M due to SWDOWN (Eqn 9), T_scale due to T2 (Eqn 10), and R_ecoG due to T2 (Eqn 11). A_min, A_opt and A_max are respectively T − T_min, T − T_opt and T − T_max. We used the T2 RMSE obtained during the validation as dT and the SWDOWN RMSE as dSW. Then, we calculated relative errors for the main vegetation types over the inner domain, which are crops, grass and shrub. We fixed the parameters EVI, P_scale, W_scale and T_scale to 1 when calculating the relative errors of the R_ecoG model to get its upper bound. Regarding relative GPP_M errors due to T2, the parameters EVI, P_scale and W_scale were not relevant because they cancelled each other out, which made their results equal to relative errors in T_scale.

3.1. Overall performance of the simulations

As shown in Fig. 3, the simulation for the sunny case shows almost no quantity of condensates (QC) and no precipitation. The simulation for the cloudy case shows higher condensate contents but almost no precipitation, and the simulation for the rainy case shows substantially higher condensate concentrations as well as higher precipitation values. YSU and ACM2 had the lowest QC during the cloudy and rainy cases, although the lowest rainy case PP occurred in WDM6 and Thompson. However, in the cloudy and rainy cases, ERA5 had QC higher than the WRF simulations that used WSM6 as their microphysics scheme, but lower than WDM6 and Thompson. Regardless, the resulting ERA5 precipitation in the rainy case was the highest.

Fig. 3.

(a)Domain and temporally averaged QC, and (b) domain-averaged PP for each case and model configuration. Solid lines are Control scenario values for reference.

Front areas denote strong gradients of equivalent potential temperature (θ_e) due to large temperature and humidity contrasts. In each simulation made during the rainy case, the θ_e = 330 K line, which approximately follows the location of the cold front in this case, had a similar shape (Fig. 4a). In particular, in YSU and ACM2 runs, the front was located ~0.5° to the south, suggesting a slower front displacement than in the other scenarios. In the case of ERA5, its θ_e = 330 K line was ~1° further to the south between 64 and 60°W and ~0.5° furhter to the north between 60 and 57°W compared with the simulations (Fig. 4b).

Fig. 4.

(a) θ_e (greyscale) at 925 hPa for the Control configuration during the rainy case at 15:00 hours UTC. The solid line represents the θ_e = 330 K contour. (b) θ_e = 330 K contour for each WRF model configuration and for ERA5 during the rainy case at 15:00 hours UTC.

These results indicate that the WRF model adequately captured the main characteristics of these three events in terms of cloud formation and precipitation. The comparison also shows that differences among WRF simulations are smaller than the difference between WRF and ERA5. This is expected behaviour because ERA5 data are based on a global model with lower effective resolution in which processes such as deep convective clouds are parametrised.

3.2. Impact of different parametrisation schemes on VPRM input variables

3.2.1. 2-m air temperature

In this subsection, we present the validation results for T2. We averaged the metrics from Eqn 1–3 over all the stations for each run. All configurations had ρ greater than 0.87, RMSE lower than 2.5°C and negative bias, except for YSU and ACM2, which had positive biases during the cloudy (both) and rainy (ACM2) cases (Fig. 5). ACM2 had the best performance of the WRF configurations during the sunny case. Dudhia and WDM6 had the lowest RMSE in the cloudy case, and ACM2 outperformed all the other configurations during the rainy case, even though the smallest bias was in YSU. However, ERA5 was the best model in the sunny and rainy cases, although it was slightly worse than Dudhia in the cloudy case.

Fig. 5.

RMSE vs bias in degrees Celsius for T2 and the mean ρ calculated over the hourly T2 time series for each run during the (a) sunny, (b) cloudy and (c) rainy cases. Solid lines represent the RMSE and bias from the Control case, the dashed line represents bias = 0°C and the symbol size is the mean ρ.

A full description of the spatial and temporal distribution of these errors is presented in the Supplementary material (Supplementary Fig. S1–S4). One advantage of high-resolution models is in the representation processes associated with complex terrain. To investigate if this has a particular impact on our simulations, Fig. 6 shows the performance of WRF simulations at stations located over 200 m in altitude and below that value. The 200-m threshold was reached after binning the stations by height into five categories, with 20% in each group. All the Student’s t-tests had P values higher than 0.15 except for WDM6 in the rainy case and for ERA5 in all cases (Fig. 6). This implies that ERA5 estimated T2 for stations located above 200 m worse than for those situated below 200 m, whereas all the WRF configurations but WDM6 during the rainy case performed similarly regardless of this threshold. Therefore, complex topography added uncertainty to ERA5 estimates that did not reduce the performance of the WRF configurations. The complete statistical results are presented in the Supplementary material (Supplementary Table S4).

Fig. 6.

Boxplots calculated over the hourly T2 time series for each run during the (a) sunny, (b) cloudy and (c) rainy cases in stations with altitude below or above 200 m. The boxes and whiskers indicate the 25th–75th and the 5th–95th percentiles respectively. White rectangles and triangles are the medians and means respectively. Horizontal lines represent the means from the Control scenario. Single and double asterisks indicate P-values from the mean t-test smaller than 0.15 or 0.01 respectively.

During the rainy case, the precipitation systems are associated with enhanced mesoscale variability in surface variables such as T2. This enhanced variability at smaller spatio-temporal scales is also linked with increased uncertainty in the simulations, which can result in higher uncertainty in the simulated CO₂ fluxes. One feature frequently connected with mesoscale convective systems is the cold pool. All simulations in the rainy case show a large and intense cold pool over the north-eastern part of the domain at 15:00 hours UTC (Fig. 7). This cold pool is well ahead of the lower-temperature region that is observed in the south-western part of the domain and that is associated with the cold front.

Fig. 7.

T2 for each model configuration (a–h) and for ERA5 (i) during the rainy case at 15:00 hours UTC. The WRF model configurations were (a) Control, (b) Dudhia, (c) Goddard, (d) Noah, (e) WDM6, (f) Thompson, (g) YSU and (h) ACM2.

The cold pool varied its intensity depending on the radiation scheme used, but had the same location in all of them (Fig. 7a–c). However, in Thompson and WDM6, the cold pool extended to the south, merging with the colder air behind the cold front (Fig. 7e–f). The distance between the cold pool and the front was longest in ACM2 (Fig. 7h). Its position was very similar among the runs, where the most important differences were in Noah and YSU as it was located slightly to the north and south in these scenarios respectively (Fig 7d, g). ERA5 had the cold pool with the lowest intensity and split into two smaller pools in the west–east direction (Fig. 7i). This behaviour is expected considering that deep moist convection is not explicitly accounted for in ERA5 so that cold pool generation and maintenance processes are only partially represented.

3.2.2. Incoming shortwave radiation at the surface

In this subsection, we present the validation results for the 24-h integrated SWDOWN and for its diel cycle. Fig. 8 shows the RMSE and bias for the daily-averaged SWDOWN and the correlation coefficient computed from the hourly SWDOWN values. We averaged their values over all the stations for each run. The differences among the simulations were dominated by the radiation scheme in the sunny case, but the other parameterisations were also important in the other two cases. Error metrics were highest during the rainy case and smallest during the sunny case, whereas ρ had the opposite behaviour (Fig. 8). The only exception was WDM6 with a higher RMSE and smaller ρ in the cloudy case; this is consistent with the relatively worse performance of this configuration for T2 in this same case. In the sunny and rainy cases, Dudhia outperformed all the WRF simulations, but YSU had the best RMSE during the cloudy case. It is interesting to note that although Dudhia performs best in the sunny case in terms of shortwave radiation, it produces the worst results in the same case in terms of T2. This indicates that producing the right radiative forcing may not lead to the optimal results in terms of temperature, probably owing to model errors in terms of other processes such as turbulent mixing within the ABL or land surface processes.

Fig. 8.

RMSE vs bias for SWDOWN (MJ m⁻² day⁻¹) calculated over the daily time series, and the mean ρ calculated over the hourly time series for each different WRF configurations and for the ERA5 during the (a) sunny, (b) cloudy and (c) rainy cases. Solid lines represent the RMSE and bias from the Control case, the dashed line indicates when the bias = 0 MJ m⁻² day⁻¹, and the circle size is the mean ρ.

All the models simulated SWDOWN with similar biases during the sunny case, but there were large differences between the models in the cloudy and rainy cases. For instance, for the station with the largest differences across the configurations, located at −31.85°S, 60.54°W, its bias respectively ranged from −20.6 to 1.9 MJ m⁻² day⁻¹ and from −4.0 to 12.2 MJ m⁻² day⁻¹ in the cloudy and rainy cases. The SWDOWN diel cycle was similar to measured radiation in the sunny case, but the differences in its patterns were larger during the other two cases. A detailed presentation of the figures that illustrate these examples is given in the Supplementary material (Supplementary Fig. S5–S8).

3.2.3. Total precipitation

The simulated PP was validated by applying Eqns 1, 2 to 24-h accumulated precipitation. As this variable has a higher small-scale spatial variability than T2 and SWDOWN, WRF-simulated PP was interpolated to the station location using the average of the closest four grid points to each station.

In the sunny and cloudy cases, PP errors were below 0.5 mm day⁻¹; thus, the analysis was focused on the rainy case. In this case, all the simulations had a positive bias except WDM6, Thompson and ACM2, which respectively had negative, negative and zero biases (Fig. 9). However, the RMSE of PP was below 12.7 mm for all the simulations, which is of the same order of magnitude as the average of 12.2 mm measured across all the stations. The lowest RMSE of the WRF simulations was achieved in Goddard, with a value of 7.5 mm day⁻¹. WDM6 and Thompson produced the largest difference in terms of RMSE with respect to the Control run. This is expected because the cloud microphysics parameterisation is the most closely related to rain. However, other schemes such as Noah (land surface processes) and Dudhia (radiation processes) produce similar impacts in terms of RMSE. This behaviour may be related to the complex dynamics associated with precipitating systems, which can be significantly affected by the ABL structure, which in turn depends on surface fluxes and radiative forcing. Moreover, the chaotic behaviour of the atmosphere is exacerbated in these situations owing to the existence of rapid instabilities that accelerate error growth.

Fig. 9.

Bias and RMSE (mm day⁻¹) during the rainy case. Solid lines are the Control scenario values for reference.

Regarding the spatial distribution of the errors, all the WRF simulations had a positive bias at most stations but a negative bias over the ones located near the north-east corner of the domain, where a pre-frontal mesoscale convective system took place. The largest differences among the WRF simulations were observed in this region. The spatial distribution of PP errors is presented in the Supplementary material (Supplementary Fig. S9–S10).

3.2.4. Atmospheric boundary layer structure

To validate how WRF simulations represent the ABL structure, we applied Eqn 1 and 2 to compare model outputs and observations vertically interpolated to a common set of 10 vertical levels between 975 and 750 hPa, evenly spaced 25 hPa apart. Error metrics were then averaged over all the available soundings. The resulting potential temperature RMSE and bias within the ABL were respectively below 1.3, 1.2 and 2.8°C during the sunny, cloudy and rainy cases. ACM2 had the best metrics among the WRF simulations, but ERA5 outperformed all of them. The only exceptions were the cloudy case bias where YSU was slightly better than ACM2 and the rainy case RMSE where ACM2 was better than ERA5.

Fig. 10 shows the vertical profiles of RMSE and bias averaged over all available soundings. Potential temperature bias and RMSE below 900 hPa were respectively between −3.5 and 3.0°C, and below 6.0°C. They were larger below 950 hPa, except in the rainy case, where the highest errors were between 950 and 850 hPa. The three cases had cold biases for all the WRF simulations and vertical levels, except for YSU and ACM2, which were hotter near the surface. These results are consistent with the warm T2 bias presented earlier for these two configurations. The WRF simulations had the worst performance in the rainy case, but ERA5 had the worst score in the sunny case. However, dew point temperatures showed the opposite behaviour. They were better below 950 hPa than in the upper layers, and they had the lowest errors in the rainy case. Readers interested in the individual soundings are referred to the Supplementary material (Supplementary Fig. S11–S13).

Fig. 10.

Potential temperature (black) and dew point temperature (grey) bias (left), and RMSE (right) averaged over all the soundings from the (a, b) sunny, (c, d) cloudy and (e, f) rainy cases.

3.3. Impact of VPRM input variable uncertainties on the CO₂ flux simulation

We propagated the T2 and SWDOWN RMSEs estimated in the previous sections into the R_eco and GPP models to estimate their impact on CO₂ flux simulation. In the case of the daily accumulated SWDOWN, we divided its RMSE by 13.5 h to convert it to the appropriate units required by WRF-VPRM (W m⁻² s⁻¹). We chose 13.5 h because it is close to the number of daytime hours for the domains and cases simulated. In the same way, we analysed the relative errors due to T2 and SWDOWN within the relevant ranges of these variables for the domain and cases simulated, which were ranged from 8 to 38°C and from 50 to 2500 W m⁻². The following subsections illustrate R_ecoM, R_ecoG and GPP_M errors through Fig. 11 and 12.

Fig. 11.

Relative errors estimated for GPP_M and for the linear and polynomial R_eco models due to errors in T2 and as a function of T2 for each run. Solid lines are the medians of the configurations and the hatched areas are bounded by their minimum and maximum errors. The upper (a, b, c), middle (d, e, f) and lower rows (g, h, i) correspond to the sunny, cloudy and rainy cases respectively. The left (a, d, g), middle (b, e, h) and right columns (c, f, i) correspond to the crop, grass and shrub vegetation types respectively.

Fig. 12.

As in Fig. 11, but for the relative errors in GPP_M due to the error in SWDOWN and as a function of SWDOWN. The upper (a, b, c), middle (d, e, f) and lower rows (g, h, i) correspond to the sunny, cloudy, and rainy cases respectively. The left (a, d, g), middle (b, e, h) and right columns (c, f, i) correspond to the crop, grass and shrub vegetation types respectively.

4.1. Sensitivity of the most relevant variables involved in the CO₂ flux simulation to the choice of model configuration

Hourly T2 errors from the WRF simulations made in this study were smaller than 2.5°C in all cases (Fig. 7). On average, ERA5 performed better than WRF with errors below 1.9°C, but its performance was worse than WRF for stations with complex topography. Over the MPRArg, Ruiz et al. (2010) tested different schemes in lower-resolution forecasts (40 km) and obtained a spatially averaged RMSE between 2.4 and 3.1°C in the hourly T2 prediction for all schemes evaluated. Another study from the same region, Müller et al. (2016), evaluated 15-km resolution WRF forecasts for a 2-year period between 2012 and 2014. They reported absolute mean errors of daily mean temperatures up to 2.7°C. In the same way, Dillon et al. (2021) validated 10-km resolution WRF forecasts with different microphysics and ABL schemes during 44 days from 2018 over the MPRArg. Their averaged RMSE and bias over the 36-h forecasts were respectively below 3.0 and 1.0°C. In summary, other studies from the region had errors larger than in our study. Given that other sensitivity analyses conducted in this region and documented in the scientific literature exhibit lower resolution, a comparison will be drawn with simulations from different regions such as the United States, China and Greece (Giannaros et al. 2013; Wyszogrodzki et al. 2013; Ju et al. 2022). These authors performed simulations with respective horizontal grid spacings of 4, 3 and 2 km, and the maximum reported error in T2 was 3.1°C.

Daily integrated SWDOWN errors were respectively smaller than 2.2, 6.7 and 7.8 MJ m⁻² day⁻¹ during sunny, cloudy and rainy situations (Fig. 8). Over the southern portion of the domain, monthly mean SWDOWN ranges between 8.0 MJ m⁻² day⁻¹ in winter and 30.0 MJ m⁻² day⁻¹ in summer, and has been increasing over the past decades (Fernández et al. 2021). Numerical estimations for this variable over the region are scarce and involve models different from WRF. Among them, Pessacg et al. (2014) validated several regional models with a grid spacing of 50 km using satellite data. They found biases up to 1.7 MJ m⁻² day⁻¹ for the period 1990–2008. In other regions, there are numerical studies focusing on cloudy cases. For instance, 1.3-km resolution simulations of the AROME model over France during 2020 produced a maximum RMSE of 4.8 MJ m⁻² day⁻¹ in cases where there were clouds in the model and in the observations (Magnaldo et al. 2023). As another example, 4- and 1.5-km resolution runs of the Met Office Unified Model (MetUM) during 2017–18 over South Africa had respective errors of 4.6 and 5.2 MJ m⁻² day⁻¹ (Mendes et al. 2023). SWDOWN estimations by WRF have been validated in regions different from the MPRArg. Examples of these studies involving high-resolution WRF simulations include a 5-km resolution experiment over the Iberian Peninsula during the period 2000–2009. It produced daily SWDOWN RMSEs up to 3.1 MJ m⁻² day⁻¹ (Perdigão et al. 2017). Also, a validation of 3-km resolution simulations over north-eastern Germany in 87 frontal passages resulted in a maximum RMSE of 9.3 MJ m⁻² day⁻¹ (Mierzwiak et al. 2023). In all the numerical experiments mentioned, we assumed an average of 12 daytime hours to extrapolate their results from watts per square metre (W m⁻²) to megajoules per square metre per day (MJ m⁻² day⁻¹). In consequence, from the information stated before, SWDOWN errors from the present study share the same order of magnitude as other simulations and expand previous works by providing high-resolution WRF simulations over the region.

Accumulated PP forecasts inside the MPRArg have been studied using WRF in Ruiz et al. (2010), Müller et al. (2016), Casaretto et al. (2022) and Dillon et al. (2021). In our rainy case, PP averaged 12.2 mm across all stations, and its bias and RMSE were respectively bounded between −2.5 and 2.7 mm day⁻¹ and between 7.5 and 12.7 mm day⁻¹ (Fig. 9). This range is smaller than the bias range found in Ruiz et al. (2010), which was between −10.0 and 15.0 mm day⁻¹ in the MPRArg during summer (December, January and February) with a 40-km resolution model, and higher than the 5.0 mm day⁻¹ from Müller et al. (2016) with a 15-km resolution model in a 2-year period between 2012 and 2014. In the same way, using 4- and 3-km resolution models, Casaretto et al. (2022) obtained PP biases between −5.0 and 5.0 mm day⁻¹ during 43 days of November–December in a smaller region with complex terrain inside the MPRArg, with an average of 3.3 mm day⁻¹. However, Dillon et al. (2021) compared 10-km resolution WRF forecasts with PP satellite estimates over 44 days inside the MPRArg during 2018. They found the average bias of precipitation events to be below 3.0 mm grid⁻¹. Even though PP is particularly sensitive to the method of validation chosen, like using nearest-neighbour interpolation with one grid point or averaging a few grid points surrounding a station, daily PP from our simulations is in agreement with other assessments from the same region.

Temperature biases and RMSEs within the ABL were respectively between −3.5 and 3.0°C, and below 6.0°C (Fig. 10). Over the region, the ABL structure has been analysed in Ruiz et al. (2010) for YSU and MYJ schemes using a 40-km resolution model and soundings from Resistencia and Santiago del Estero stations (27.45°S, 59.05°W; 27.77°S, 64.30°W). In Resistencia, the biases were bounded by 2.0°C and they were cold near the surface and warm at higher levels. However, the bias at Santiago del Estero was cold, reaching up to 3.0°C within the ABL, with YSU errors being warmer at lower heights, as happened in our work. The ABL structure showed greater sensitivity to the choice of both ABL and surface layer schemes in the rainy case. Therefore, warmer, wetter and more unstable ABLs in YSU and ACM2 were generated. The ABL structure has been studied with higher-resolution models in other regions of the world. Nevertheless, our RMSE values were higher than those calculated by Coniglio et al. (2013) where the RMSE for the ABL potential temperature was below 1.0°C for a 4-km resolution model over the United States. Similarly, their ABL structure sensitivity was smaller than in our case. In their experiments, the errors were almost constant above 100 m, which is similar behaviour to that found for the sunny case in the present study. This is an expected result because Coniglio et al. (2013) considered only the soundings obtained under clear-sky conditions. Our study expands their analysis by showing that potential temperature errors can show a different vertical structure within the ABL under synoptically active conditions.

Each target variable was more sensitive to a particular physical process depending on the synoptic situation. For instance, the ABL–surface layer and radiation schemes analysed had the largest T2 and SWDOWN bias spread in most cases. This result highlights the importance of choosing appropriate ABL–surface layer and radiation schemes when analysing CO₂ fluxes in simulations. Selecting ill-performing parameterisations for these processes can increase the CO₂ flux errors. For that reason, it is important for future studies to carefully consider their selected schemes. Specifically, based on our results, we recommend including a set of ABL parameterisations with local and non-local mixing schemes. Future research can enhance the robustness of this analysis by incorporating additional case studies, extending time periods and using independent validation data. It would also be beneficial to evaluate the performance of similar analyses using different combinations of schemes than those selected here. Attention should be paid to the ABL–surface layer combination, as Maroneze et al. (2021) demonstrated that they significantly affect the simulation of the stable boundary layer.

4.2. Uncertainties of CO₂ flux simulations

In this study, we investigated the impact of temperature and SWDOWN errors on CO₂ fluxes. The R_eco and GPP models from the original VPRM – R_ecoM and GPP_M – had larger relative errors when either T2 or SWDOWN approached zero, which are common daytime conditions in the anticyclonic flow behind a cold front. The former was due to R_eco and GPP being smaller at colder temperatures when the vegetation is close to its T_min. In a similar way, as plants decrease their photosynthesis rate when SWDOWN is scarce, GPP_M relative errors are larger at sunrise or sunset. As a consequence, when both effects are considered, the highest relative errors are expected during sunrise and just behind the front zone, where temperature and incoming radiation are low. This outcome can be mitigated if a respiration model with low uncertainties at cold temperatures is used, like R_ecoG, instead of R_ecoM.

Li et al. (2020) performed a sensitivity analysis of the original VPRM to its parameters. They found the highest NEE uncertainties to be up to 60%, occurring near the local solar noon when SWDOWN was highest. Our study complements their results by pointing out that uncertainties in the environmental variables simulated by a numerical model increase the NEE relative errors during the morning and evening, accounting for ~20% of the cases assessed in their work (Fig. 11 and 12). Furthermore, at least over the MPRArg, the CO₂ flux uncertainty will be larger in areas with complex topography during rainy cases because T2 and SWDOWN errors were larger. Therefore, when both sources of uncertainty are taken into account, daytime relative errors may be greater than 60%.

Even though the simulations from our study were similar, slight differences will change the fluxes. For instance, as the nocturnal cooling rate in ERA5 was slower than in the WRF simulations in most locations, nighttime R_eco fluxes calculated with ERA5 will be higher, leading to a more positive CO₂ daily budget. By contrast, the nighttime carbon balance will be more negative in simulations that use Noah as the land surface scheme because its nocturnal cooling rate was faster.

Regarding the effect of precipitation on near-surface temperatures, there were runs with larger cold pools, like Thompson and WDM6, that result in smaller R_eco fluxes. Moreover, the position of the cold pool may also influence R_eco in the presence of high spatial variability in the vegetation type. For example, the cold pool simulated in the Control setting covered a smaller percentage of grass than in YSU because it was slightly shifted to the west. Therefore, differences in the location of the cold pools result in uncertainties in the NEE fluxes. Mislocation of cold pools is not the objective of our study and should be studied in depth in future work.

It is important to highlight that most of the runs that used MYJ as the ABL scheme underestimated T2 and overestimated SWDOWN. This will result in R_eco and GPP fluxes smaller and larger than the real ones respectively, and the combined effect over NEE will depend on the specific values of both variables. This suggests that ABL schemes with non-local mixing, like YSU or ACM2, are more appropriate for the region. Even though these schemes had positive T2 biases during the cloudy and rainy cases, they would result in lower R_eco and GPP errors.

The CO₂ flux uncertainty depended on the dominant synoptic situation of the case. In this way, GPP errors in cloudy weather were slightly larger than in fair (sunny) weather but much smaller than in rainy conditions. This implies that the uncertainty added by the front and the mesoscale convective system was more important than that produced by clouds alone. In consequence, future GPP estimates should consider each particular synoptic situation separately and then adjust its uncertainty level accordingly.

The Intergovernmental Panel on Climate Change (IPCC) includes CO₂ atmospheric concentration estimates for the near, mid and long-term futures based on global climate models able to account for CO₂ emissions from biogenic sources (Canadell et al. 2021; Friedlingstein et al. 2022). These models utilise flux inventories as prior information through a Bayesian approach. The priors can be based on observations, simulations from models like VPRM or a combination of both. As we have shown here, the uncertainty in the CO₂ fluxes can be very high and may depend on the synoptic situation. This fact added to other sources of error not analysed in this study, like model selection, parameters and resolution, contribute to the total uncertainty. Therefore, our results strongly suggest a more extensive investigation of the uncertainties associated with the representation of physical processes in the global and regional climate models, in particular those that are directly associated with the variables that modulate surface CO₂ emissions from biogenic sources. We recommend that future studies employ global and/or regional climate models capable of accounting for CO₂ emissions from biogenic sources to conduct sensitivity analyses similar to ours, while carefully considering all sources of error, particularly the uncertainty associated with the parametrisation of unresolved scale processes in the CO₂ flux.

4.3. Most appropriate WRF schemes for simulating CO₂ fluxes on each dominant synoptic situation

The performance of the WRF parameterisations tested in this study characterised the uncertainty associated with the selection of different parameterisations on atmospheric variables, essential for computing surface CO₂ fluxes over the MPRArg. Overall, all performed better than ERA5 in stations surrounded by complex topography. Among the radiation schemes, RRTMG_SW and Goddard outperformed Dudhia in the estimation of near-surface temperatures, but the latter produced more accurate estimates of shortwave radiation fluxes during the three synoptic situations simulated. The Noah MP land surface model produced better simulated near-surface temperatures and shortwave radiative fluxes than Noah. Also, Thompson microphysics parameterisations were better than WDM6 in the simulation of near-surface temperature. WSM6 and WDM6 were the best microphysics schemes to simulate shortwave radiation fluxes during the cloudy and rainy cases respectively. Non-local ABL schemes, YSU and ACM2, surpassed MYJ when simulating near-surface temperatures in the sunny and rainy cases. The opposite occurred in the cloudy case and with shortwave radiation fluxes in all situations. All in all, even when the sensitivity to a change in the parameterisation chosen was very small in many cases, it was possible to provide insights about an appropriate set of WRF schemes for each particular synoptic situation that can be used in future study cases.

Overall, relative NEE errors due to T2 were similar in all the configurations and cases tested, with the exception of WDM6 in the cloudy case, which propagated larger errors. When considering relative NEE errors due to SWDOWN, Dudhia outperformed the other configurations in the sunny and rainy cases, whereas the same happened for ACM2 in the cloudy case. Relative errors were smaller when T2 was close to T_opt and when SWDOWN was high; thus, all the configurations tested in this analysis will have a similar performance in study cases with these atmospheric conditions, like, for example, temperate and sunny cases. However, the spread of relative errors due to T2 and SWDOWN was larger when T2 was close to T_max or T_min and when SWDOWN was low respectively. These conditions happen at sunrise and sunset and just behind the front zone. Therefore, it is possible to get smaller relative NEE errors in study cases involving frontal passages when choosing Dudhia or ACM2 over the other schemes tested.