The DNDC model is a process-based biogeochemical model that simulates crop growth, soil C, and N cycling and associated biosphere–atmosphere exchange of C and N trace gases, as well as N leaching (Li et al., 1992a,b; Li, 2000; Giltrap et al., 2010; Deng et al., 2013). The DNDC model runs at daily/subdaily time intervals using initial conditions such as soil properties (bulk density, pH, texture, and soil organic carbon [SOC]) and meteorological drivers (daily air temperature, precipitation, and radiation), as well as agricultural management (planting/harvest date, tillage, and fertilization). Based on interacting simulations of six submodules (soil climate; plant growth, mineralization, fermentation, nitrification, and denitrification) and input information, DNDC predicts soil layer-specific environmental conditions of substrate C and N availability, soil temperature and moisture. In addition, it predicts oxygen concentrations for partitioning of anaerobic/aerobic microsites, which drive microbial C and N turnover processes of mineralization, nitrification and denitrification as well as associated losses of N₂O. The model can be run in either site or regional mode (Lugato et al., 2010; Giltrap et al., 2013).

The DNDC model was initially developed to predict soil C and N biogeochemistry in temperate agroecosystems (Li, 2000) but has largely been untested in tropical regions. Recent work has shown that the very deep soils (reaching 8 m or more) in the Mato Grosso region have high hydraulic conductivity and large pores despite fine soil textures (Scheffler et al., 2011; Dias et al., 2015). We developed new DNDC parameterizations to account for the hydraulic properties of tropical soils by adjusting soil hydraulic conductivity and soil porosity (Supplementary Table 1) and recalculating soil hydraulic parameters for all texture classes based on pedotransfer functions (PTFs) for Brazilian soils (Medrado and Lima, 2014) in equations for water-filled pore space (WFPS) and wilting point, respectively:

We compared both the standard and local parameterizations models against field observations (O'Connell, 2015; Jankowski et al., 2018) in western MT (Figure 1, Tanguro Ranch, 12°59′S, 52°23′W). Tanguro Ranch has a humid tropical climate (Köppen's Aw, Kottek et al., 2006) with a rainy season from October to April and a severely dry season from May to September. Mean annual precipitation averages about 1,900 mm and the mean annual temperature is 25°C. At Tanguro Ranch, soils are medium textured, highly weathered, base-poor ustic Oxisols (Latossolo vermelho-amarelo distrófico in the Brazilian Classification System). Soils and the depth of the groundwater table are deep (10–30 m), very well drained on plateaus, and have a mean soil texture of 55% sand, 2% silt and 43% clay. Soil mean pH was 3.9 under native forest and 6.0 under established croplands (Scheffler et al., 2011; Riskin et al., 2013).

Beginning in 2010, large areas of Tanguro Ranch were converted from single cropping of soybeans to double cropping of soybeans followed by maize in a single rainy season. At Tanguro, as in most of Mato Grosso, soybeans are planted at the start of the rainy season (Oct–Nov), fertilized with phosphorus (no nitrogen), and harvested in late January to early February. Maize is planted immediately following soybean harvest (late January) and harvested at the end of the rainy season (April). N fertilizer is typically applied to maize as urea or ammonium-nitrate at the time of seeding (~5 kg N ha⁻¹). Roughly 20–30 days after planting, ~80 kg N ha⁻¹ is broadcast as urea, for ~85 kg N ha⁻¹ total N fertilizer plus any residual N from soybeans of previous years.

To validate DNDC estimates of N₂O emissions, we used two distinct datasets of N₂O emissions from soybean-maize fields. First, to estimate emissions from the soy phase of fields planted in soybean-maize, we measured N₂O emissions from soybean-maize fields across Tanguro Ranch regularly from July 2013 through November 2014 (O'Connell, 2015). The precise amount and timing of N fertilization were not available for this dataset so only the measurements under soybean cropping were used for model validation, which involved no N fertilization. Emissions were measured using a closed chamber method (Venterea et al., 2005) where gas samples were drawn with a syringe four times after the chamber was closed. Nitrous oxide concentrations in the samples were determined using a gas chromatograph in the laboratory. The soil-to-atmosphere fluxes (dN₂O/dt) were calculated using a linear regression of N₂O concentrations against time. The fluxes were computed based on air temperature and air pressure at the time of the measurement and total volume and basal area of the chamber (O'Connell, 2015).

Next, to account for N fertilization effects on N₂O emissions from maize, we used data collected from a field-level fertilizer manipulation experiment during the maize-growing phase of the cropping cycle in a single field at Tanguro Ranch (Jankowski et al., 2018). This experiment measured N₂O emissions from established five fertilizer treatments (0, 80, 120, 160, and 200 kg N ha⁻¹), laid out within five replicate blocks that consisted of 33 m² plots (5.5 × 6 m) with 1 m buffers between plots. Maize was planted on 30 January 2015, the day after soybean harvest. All plots received 9 and 31 kg K ha⁻¹ and all N-fertilized treatments received 5 kg N ha⁻¹ as NH₄NO₃ with the maize seeds sowed by machine planters. A second dose of N was added by hand (urea) to the N-fertilized treatments 21 days after planting (concurrent with application date on the surrounding field) to create total N application rates of 80, 120, 160, and 200 kg N ha⁻¹. All treatments and controls had five replicates randomly assigned within blocks.

The plots were sampled for N₂O emissions with static chambers in the same manner as above and measurements were taken daily for 1 week following both planting and fertilization and weekly or bi-weekly until harvest on 2 June 2015. Fluxes were calculated based on O'Connell (2015). Data loggers were installed in each N fertilizer treatment to measure soil temperature and water content at a depth of 10 cm. Logger measurements were validated by field measurements of soil gravimetric moisture during gas flux measurements. We converted soil moisture content to percent of WFPS based on gravimetric moisture, bulk density and particle density. Maize biomass (dry weight from 10 maize plants per plot) was determined by harvest 125 days after seeding. Daily precipitation and average temperature were collected from a datalogger at the experimental plots (Jankowski et al., 2018).

The DNDC model performance was evaluated by the normalized root mean square prediction error (RMSE) and coefficient of determination (r²) based on the following equations (Moriasi et al., 2007):

where Xmea is the measured value and Xsim is the simulated value. X_mea is the average value of field measurements and X_sim is the average value of model simulations.

After the tests, the local parameterized DNDC model was used for a sensitivity test. DNDC was run for the same site but with varied climate, soil, and management conditions. The purpose of the sensitivity test was to identify the most sensitive factors that could influence the N₂O emissions from the target ecosystem.

After site validation, DNDC was used to simulate crop yields and N₂O emissions from single-cropping soybean and double-cropping soybean-maize for the entire state of Mato Grosso from 2001 to 2010 by compiling spatially distributed input data in a GIS database. To meet DNDC requirements for regional runs, the state was divided into 5 × 5 km grids of uniform agricultural management, soil, and climate. Each grid contained information about soil properties (clay content, initial SOC, pH, and bulk density), land use management (crops, tillage, fertilization, planting, and harvesting dates) and climate (precipitation, temperature, radiation, wind speed, and air humidity). Total N₂O emissions and crop yields for the entire state were a result of the sum of the results from all grids.

Spatial data on croplands were provided by Spera et al. (2016). This data set delineated cropland areas, single and double cropping and specific crop types within each cropping pattern. The data set was derived from time series analysis of MODIS EVI with 250 meter resolution and provided a very high degree of accuracy (Galford et al., 2008; Spera et al., 2014).

We assumed all single soybean and double soybean-maize cropping systems were, respectively, cultivated under conventional tillage (plowing soil at 0–20 cm seven days before the planting date) and no-tillage systems called “plantio direto” that only disturbs the soil in sowing or planting and is widely adopted in MT (Galford et al., 2010; Arvor et al., 2012). In addition, planting, and harvesting dates were assumed to be the same as those carried out by Tanguro Ranch for the 2013/2014 crop cultivation (Table 1). Soybean crops usually receive an input of 7 kg N ha⁻¹ in MT (Raucci et al., 2015). Soybeans typically do not require nitrogen additives, but most farmers in the region add this modest amount as “insurance” (Galford et al., 2010; Raucci et al., 2015).

Maize crops do require nitrogen fertilization even when planted immediately following soybeans. Studies show that optimal maize productivity in this region may require up to 120 kg N ha⁻¹ (Mar et al., 2003; Rezende Pereira et al., 2009), some of which may be provided by N fixed by soybeans grown as the first crop (Galford et al., 2010; Jankowski et al., 2018). At current productivity levels, maize crops in the Amazon and Cerrado regions need nitrogen fertilizer doses that range from 0 to 77 kg N ha⁻¹ (Mar et al., 2003; Cruz et al., 2005; Souza and Soratto, 2006). Applied fertilizer doses ranging from 0 to 10 kg N ha⁻¹ and 0 to 75 kg N ha⁻¹ were tested in the DNDC model as minimum and maximum values of N-fertilizer rates for soybean and maize, respectively (Table 1).

Soil organic C concentration, pH, clay content and bulk density were obtained using the MT soil map (IBGE, 2015), which was divided into three categories: oxisol, ultisol, and sandy. Physical and chemical properties of these soil types were extracted from a database of soil sampling carried out in several cropping systems in MT (Belizario, 2011) (Table 2).

Maximum and minimum air temperature, daily rainfall, wind speed, solar radiation and air humidity data from 2001 to 2014 were retrieved from the Brazilian National Institute of Meteorology (INMET, 2016). Because Mato Grosso had only 11 long-term meteorological stations over the study period (2001 to 2010) all retrieved products were spatially interpolated using the Voronoi diagram method to obtain 5 × 5 km gridded datasets (Supplementary Figures 1, 2).

Superimposing climate, soil, and management characteristics resulted in a total of 1,298 homogenous simulation grid cells for the entire decade of regional inventory. As the soil C pools were not necessarily in equilibrium at the start of each simulation, each year was run as a separate simulation. This was done rather than running a single 10-year simulation to prevent long-term changes in soil C from affecting the results (Giltrap et al., 2010). In addition, each grid was run in DNDC in “site” mode. This is preferable to “regional” mode because its use is more transparent and flexible (Perlman et al., 2013).

The error in the calculated N₂O emissions was taken as the sum of the variation in soil parameters and N fertilizer rate reported in the literature that addressed the region under analysis (Tables 1, 2). This produced the maximum error range for these two sources. The errors introduced by the assumptions made about farm management practices or the climate data sets were not included in this paper.

Measured N₂O emissions of all N-fertilizer treatments of 2015 field experiment (Figure 2; Jankowski et al., 2018) increased only following the second (broadcast) N-fertilization of maize. Emissions peaked on the third day after fertilizer application (from 0.12 ± 0.10 to 0.36 ± 0.28 mg m⁻² h⁻¹), but by the seventh day emissions were similar to the control. No significant emissions were detected from seeding N-fertilization, probably because of the low rate of application (5 kg N ha⁻¹).

Cumulative measured N₂O emissions from the N fertilization experiment ranged from 0.27 to 0.75 kg N ha⁻¹ and were consistent with findings that tropical soils tend to emit less N₂O compared with temperate soils fertilized at similar rates (Hickman et al., 2015; Meurer et al., 2016; Jankowski et al., 2018). Estimates of fertilizer N lost as N₂O-N during this period ranged from 0.10 to 0.24% (Table 3; Jankowski et al., 2018) and were well below the 1% Intergovernmental Panel on Climate Change default emission factor (IPCC, 2006). These results agreed with findings by Meurer et al. (2016), who reported that despite N fertilization increasing N₂O emissions for short periods (3 to 7 days after application), the annual responses of N₂O emissions to N application were negligible for application rates less than 100 kg N ha⁻¹.

In both the standard and local parameterized models, simulated N₂O emissions followed the same pattern of measured values for all N-fertilizer treatments (Figure 2). Despite good agreement with the timing of N₂O peak fluxes, the unmodified DNDC model overestimated the magnitude of emissions by up to 2.5 times until adjustments for soil physical properties were made (Figure 2). After parameterizing DNDC's hydrologic variables (see Soil Temperature and Moisture, Figure 3), modeled N₂O emissions agreed closely in timing and magnitude with field measurements (Figure 2). After model calibration, r² and the RMSE were substantially improved (Table 3). The values of RMSE between the standard and local parameterized models and observed N₂O emissions ranged between 1.9–56.2% and 2.0–4.4%, respectively. Improved correlations between measured and simulated values showed that the hydrological modifications improved simulations of N₂O emissions (Table 3).

Simulations revealed two important insights. First, the model captured the temporal pattern of N₂O emissions following both fertilization events. The inputs from biological N fixation by soybean and fertilizer application (5 kg N ha⁻¹) at the time of maize seeding produced almost no measurable N₂O emissions; both the standard and local parameterized models captured that dynamic. Next, the improvement in the model that occurred with calibration to soil conditions at Tanguro Ranch indicated that soil texture and physical properties were key factors for more accurately simulating N₂O fluxes.

There are remaining uncertainties in how well the model estimates fluxes with increasing fertilizer application. Compared with measured values, simulated accumulated N₂O emissions were close to the 1:1 line for the two lowest levels of fertilizer application but the model overestimated emissions for the higher fertilizer application levels (Table 3; Figure 4). Thus, the model was very accurate for the current regional application rates of ~70 kg/ha (Galford et al., 2010), but may overestimate emissions if application rates increase substantially. Regional crop types were mapped with 87% accuracy (Spera et al., 2016), introducing additional error to regional emissions estimates but perhaps countering slight over estimates of GHG emissions at high fertilization rates.

Field surveys of soil WFPS (O'Connell, 2015) was greater between October and April (rainy season) than between July and September (dry season). Accumulated rainfall averaged 1779 and 264 mm over these seasons, respectively. Soil temperatures ranged from 20 to 30°C but showed no distinct seasonal pattern (Figure 3). Within the corn portion of the 2015 growing season, the soil WFPS in the second experimental field was higher in February-March (rainy season) than April-May (dry season). Accumulated annual rainfall averaged 360 and 98 mm during these periods, respectively. Soil temperature averaged 27°C and differed slightly from January to June 2015 (Figure 3).

DNDC originally overestimated measured WFPS in the study system (Figure 3), howerver adjusting parameterizations based on infiltration measured from Tanguro Ranch improved the ability to simulate average WFPS and r² increased from 0.35 to 0.65 (Figure 3). The overestimation was consistent with a temperate-based parameterization for clay soils (55% clay in this case) as being very poorly drained. Many soils with high clay content in tropical croplands have increased permeability caused to the micro-aggregated structure of Oxisols (Tomasella and Hodnett, 1996, 2005).

Overall, model performance measures of simulated soil environmental conditions were comparable to other studies that simulated soil temperature and water dynamics (Li et al., 2006, 2012; Giltrap et al., 2010; Kim et al., 2014, 2015). Neither the original DNDC or standard and local parameterizations of DNDC captured diurnal temperature ranges observed in field data. We improved the ability of the DNDC to estimate WFPS, even though the soil temperature estimates were better than soil moisture estimates with the model adjustments for tropical conditions (Figure 3). Our study further confirms the finding that model performance in predictions of soil temperature is better than for soil moisture. Soil moisture simulations are more sensitive and dependent on a variety of inputs of physical characteristics (e.g., bulk density and saturated hydraulic conductivity) with high spatial variability that is challenging to capture. Both soil temperature and moisture deserve further refinement in future model development for Brazilian and other tropical cropland soils.

Simulated maize yields generally matched the field measurements of grain production over a range of N fertilization (Figure 5). Modeled soybean yield averaged 3.3 t DW ha⁻¹, which agreed well with Tanguro Ranch records and official statistics of MT state (CONAB, 2015) of 3.5 and 3.1 t DW ha⁻¹ for 2015, respectively. As the field measurements and simulations both indicated, crop yields did not significantly increase once nitrogen application rates exceeded ~100 kg N ha⁻¹ (Figure 5). This finding agreed with other nitrogen fertilizer trials for maize crops in Mato Grosso (Kappes, 2013; Kappes et al., 2013).

The sensitivity of DNDC simulated N₂O emissions to variations in input parameters has been widely tested at both site (Li et al., 2010) and regional scales (Lugato et al., 2010; Giltrap et al., 2013). These studies indicated that simulated N₂O emissions are most sensitive to input parameters of climate, soil texture, SOC, bulk density, pH and applied N rate. We examined effects on simulated N₂O emissions by varying each of these inputs, while keeping all other attributes constant.

Results from the sensitivity tests indicated that among all the environmental factors we considered, SOC content and bulk density had the largest effects on N₂O emission rates (Figure 6). When SOC increased from 1.8 to 2.3%, the annual N₂O emission rate increased from 0.5 to 0.7 kg N ha⁻¹y⁻¹. The modeled data indicated that higher SOC produced more dissolved organic carbon (DOC) and inorganic N (ammonium and nitrate) through decomposition that led to higher rates of nitrification and denitrification, the two processes that produce N₂O (Figure 6).

Soil texture and structure are relevant driving factors for N₂O emissions, mainly because they control nutrient availability and water balances. In general, soils with finer texture have a higher availability of N (Luizão et al., 1989) and consequently tend to emit higher amounts of N₂O than sandy soils (Matson et al., 1998). However, in our sensitivity analysis the model behaved just the opposite. This is likely a function of the local parameterization of the model that reflect rapid decrease of water content in fine soils, which makes them behave more like coarse-textured soils (Supplementary Table 1). As a result, the water holding capacity does not necessarily increase with an increase in clay content, and therefore nitrification is more likely to occur than denitrification. Tropical soils with high clay contents, formation of micro-aggregates, and high drainage can be expected to emit less N₂O than is reported for temperate soils (Meurer et al., 2016).

Sensitivity of simulated N₂O emissions to changes in bulk density was high. When soil density is doubled (from 0.8 to 1.6 g cm⁻³) N₂O emissions declined by 65% (Figure 6). This was a consequence of the implicit dependency of nitrification and denitrification on oxygen availability, and thus on the extent of anaerobic zones in the soil profile. An increase in bulk density will decrease total pore volume and thus reduce oxygen diffusion into the soil. This is notable following rainfall when the anaerobic volume fraction of the soil increases more slowly in soils with high bulk density (Stange et al., 2000; Kiese et al., 2005).

The simulated N₂O emissions were sensitive to both air temperature and precipitation (Figure 6) as these factors control the processes of decomposition, nitrification and denitrification in DNDC. In our case, the simulations suggested that higher precipitation induces higher rates of denitrification and higher temperatures stimulate microbial activity. However, the reverse effect can also occur. Higher precipitation may induce greater rates of nitrate leaching that compete with the denitrification process, together with elevated temperature that cause substantial reductions of soil moisture through increased transpiration. Nitrogen fertilizer application rates had a nearly linear effect on N₂O emissions. An increase in fertilizer application rate from 45 to 135 kg N ha⁻¹y⁻¹ increased the N₂O emission rate from 0.5 to 1.1 kg N ha⁻¹y⁻¹ (Figure 6).

The local parameterized DNDC captured the trends in observed data much better than the original DNDC (Table 3) but validation of models operating at daily or subdaily time steps remains difficult. DNDC was developed for more data-rich temperare agricultural systems, but in emerging tropical cropping systems there are far fewer field-based measurements of GHG emissions with experimental design of multiple levels of fertilization for crops like maize (e.g., Hickman et al., 2017; Jankowski et al., 2018). The remote nature of many tropical sites and the requirements for gas chromatography including limitations on accessibility and high costs of standard gases contribute to this data scarcity.

Model simulations showed that N₂O emissions from MT single soybean and double soybean-maizecropping systems increased almost fourfold during 2001–2010, from 1.1 ± 1.1 to 4.1 ± 3.2 Gg N–N₂O. Modeled average N₂O emissions ranged from 0.1 ± 0.1 to 0.7 ± 0.7 kg N₂O-N ha⁻¹ for single soybean and from 0.7 ± 0.7 to 1.2 ± 1.2 kg N₂O₋N ha⁻¹ for double crop soybean-maize areas over MT from 2001 to 2010 (Table 4; Figure 7). Higher average values for the double cropping system show that this recent change in cropping practice enhanced regional N₂O emissions.

Simulated regional production of maize and soybeans for 2010 varied from 5.2 ± 0.6 and 3.0 ± 0.0 t ha⁻¹ (Table 4; Figure 8), which agreed well with values of 4.0 and 3.1 t ha⁻¹, respectively, reported in national statistics (CONAB, 2015). Maize yields for 2001 to 2010 were probably underestimated because the model was parameterized with local data from a farm that has higher maize productivity when compared with the other simulated MT regions (Table 4; Figure 8).

The largest sensitivity in regional N₂O estimates arises from uncertainty in fertilizer application rates applied to soil and soil parameters, especially soil density (Table 2). In this region, soil organic carbon, clay and density potentially vary from 10% to 130%. Rates of nitrogen application varied from 0 to 85 kg N ha⁻¹ y⁻¹ (Tables 1, 2). Combined, we estimated uncertainties on the amount of nitrogen fertilizer applied to soil that ranged from 42% to 80% of total uncertainties, depending on the region. Clearly, reducing the uncertainties relative to N fertilizer is a priority for reducing uncertainties in N₂O emission estimates. In addition, the methodology applied in this work of using the maximum and minimum values of each soil property produces high uncertainties. If more data were available regarding the distribution of the soil properties within a soil class then a Monte Carlo method could be used to estimate the uncertainty in a probabilistic manner rather than using the extreme values (Giltrap et al., 2013).

Mato Grosso aims to increase maize production by 60% (to 38 million metric tons) to meet 2025 projected demands (IMEA, 2015). To accomplish this, the state is expected to expand its area planted to maize from 3 to 6 million ha over the next decade. With nearly three million ha currently under single soybean cultivation (plus three million ha double cropped with maize) and an additional five million ha expected to be added to support demand for soybeans that is estimated as requiring 14 million ha (IMEA, 2015). Our model results show that the area that currently exists as single-cropped soybean could easily accommodate double cropping systems to meet projected maize grain demand. Assuming no improvement in agronomic practices (and maintaining the current level of N fertilization), linearly extrapolating simulated results from 2001–2010 to 2025, to meet future demands of corn production by double cropping with corn three million ha of soybean fields currently single cropped would double N₂O emissions compared with 2010 levels (Figure 9).

We also evaluated whether increased N-fertilization application rates would increase maize yield sufficiently to meet maize demand. According to technical experts in the region, maize in MT is currently underfertilized with N (Kappes, 2013; Kappes et al., 2013). This is a consequence of the history of maize cultivation as a second crop starting in the 1990s, mostly for soil cover with promotion of the plantio direto system, which eventually evolved to an income generating crop. Under such conditions, maize benefited from soybeans that fixed and left behind modest amounts of N after harvest so farmers did not fertilize maize to the extent recommended. As incomes from maize crops increased and soybean prices fluctuated, fertilization of maize increased.

We hypothesized that farmers are now more likely to increase fertilization rates, and evaluated how increasing N application would impact maize productivity. We simulated N-fertilizer application levels for maize from 0 to 300 kg N ha⁻¹ and found that productivity peaked at 180 kg N ha⁻¹, after which point yield stabilized or declined slightly (Figure 10). Although this N application rate (300 kg N ha⁻¹) is almost six times greater than that currently average estimated for MT (Table 1), and therefore not likely to be applied, it is similar to maximum rates in other intensively fertilized regions. We found that a fertilizer application of 180 kg N ha⁻¹ produced 25% more maize compared with current practices utilizing 35 kg N ha⁻¹ (Figure 10). If this rate of fertilization were applied across all currently double cropped areas, however, N₂O emissions would increase by 2.6 times per ha on average (Figure 10).

Similar results have been found in experimental trials aiming at the evaluation of maize yield potential as a result of different N-fertilizer doses in Brazil. Mar et al. (2003); Jankowski et al. (2018) found maize productivity peaked at 120 kg N ha⁻¹, with productivity reductions when doses were applied above this rate. The authors argued that the over fertilization may unbalance the plant ability to uptake other nutrients as well as it enhances N losses through leaching and volatilization (Jankowski et al., 2018).

These results suggest that for MT to meet estimated maize demands in the next decade, producers will have two alternatives: expand soybean-maize double cropping intensification to an additional 3.0 million ha currently used for single soybean cropping, or further intensify soybean-maize double cropping on 2.3 million ha of single cropping soybean fields and increase N fertilization in total soybean-maize double cropping area by five times (180 kg N ha⁻¹). The second strategy would spare 0.7 million ha for other activities compared to the first alternative but may emit 35% more N₂O than the first option (Figure 9) and would certainly be more costly for farmers because of the increased fertilizer use. If the target is to meet future maize demands together with reduced N₂O emissions, our simulations suggest that to convert soybean areas currently producing a single crop to double cropping with maize, maintaining current N fertilizer application rates must be a priority over increasing the rate of N application. Furthermore, applying less N fertilizer has additional benefits, such as reducing N volatilization and leaching to water bodies that occur because of large and inefficient N fertilizer applications (Carpenter et al., 1998; Kim et al., 2015; Spera et al., 2016).