Assessing Effects of Agronomic Nitrogen Management on Crop Nitrogen Use and Nitrogen Losses in the Western Canadian Prairies

Effective agronomic nitrogen management strategies ensure optimum productivity, reduce nitrogen losses, and enhance economic profitability and environmental quality. Farmers in western Canada make key decisions on formulation, rate, timing, and placement of fertilizer nitrogen that are suitable for soils, weather, and farming operations within which they operate. Suitability of agronomic nitrogen management options are assessed by estimates from linear interpolations and extrapolations of temporally and spatially discrete field-plot measurements of nitrogen responses. Such estimates do not account for non-linear and offsetting biogeochemical feedbacks of nitrogen cycles and cannot provide comprehensive nitrogen budgets for alternative nitrogen management options. These limitations can be overcome by using process-based agro-ecosystem models that adequately simulate basic processes of nitrogen biogeochemical cycles and are rigorously tested against site observations. Ecosys is a process-based ecosystem model that successfully simulated the biogeochemical feedbacks among nitrogen, carbon, and phosphorus cycles across different agro-ecosystems. This study deployed ecosys to generate spatially and temporally continuous estimates to assess crop nitrogen use and agronomic nitrogen losses from the crop fields across Alberta for alternative nitrogen fertilizer management scenarios. The study simulated effects of four nitrogen management scenarios: fall banded urea, fall banded ESN (Environmentally Smart Nitrogen), spring banded urea, and spring banded ESN on nitrogen recovery and losses from barley fields on mid-slope landforms. These simulations were done at township grids of ~10 km × 10 km over 2011–2015 utilizing provincial soil and climate datasets. Modeled annual N2O, N2, and NH3 emissions, and nitrogen losses in surface runoff and sub-surface discharge were lower by about 25, 30, 70, and 40%, respectively, with spring banding than in fall banding across Alberta. Modeled barley yields and grain nitrogen uptake were similar in spring and fall banding, indicating agro-economic and environmental sustainability advantage of spring banding in Alberta. These modeled estimates were consistent with estimates based on plot and laboratory research for Alberta and similar prairie conditions. This study pioneered a methodology of process-based agroecosystem modeling, which is replicable and scalable to assess cumulative impacts of alternative agronomic nitrogen management options on crop production and the environment on provincial, regional, federal, continental, and global scales.


INTRODUCTION
Sustainable fertilizer use management is key for optimizing crop production. Fertilizer nitrogen (N) is a major input for crop production across western Canadian prairies (Shen et al., 2019a,b). About 75% of total Canadian fertilizer N application takes place across the three prairie provinces of Alberta, Saskatchewan, and Manitoba (Statistics Canada, 2016). Alberta contains one third of the agricultural land area in Canada and encompasses a wide range of agro-climatic and soil conditions (Statistics Canada, 2016). The agricultural zone (i.e., white zone) of Alberta encompasses about 42% of total agricultural lands of Canada, which includes about 21 million ha of farmlands (Statistics Canada, 2016). However, nitrogen use efficiency (NUE) on Alberta farmlands is low, as only about 30-50% of the total applied fertilizer N is recovered in crops (Cassman et al., 2002;Janzen et al., 2003). Nitrogen fertilizers not taken up by the crop can either be immobilized by microbes or vulnerable to losses as ammonia (NH 3 ), nitrate (NO 3 − ), nitrous oxide (N 2 O), and di-nitrogen (N 2 ) through volatilization, leaching, and denitrification processes, respectively (Janzen et al., 2003;Qiao et al., 2015;Shrestha et al., 2018;Cui and Wang, 2019;Shen et al., 2019a,b). Losses of N can contaminate water bodies, release potent greenhouse gas (GHG) (e.g., N 2 O) and create particulate aerosols (e.g., NH 3 ) that impact human health (Forster et al., 2007;Qiao et al., 2015;Shrestha et al., 2018;Cui and Wang, 2019;Shen et al., 2019a,b). Minimizing N losses from agro-ecosystems could thus provide opportunities for reducing the environmental footprint of crop production in western Canada (Shen et al., 2019a,b). In addition, reduced fertilizer N loss can enhance fertilizer NUE, which minimizes fertilizer requirements and saves input costs (Grant and Wu, 2008). Optimal fertilizer N management ensures high crop production for human and livestock consumption and is the foundation of value-added sustainable agricultural products. Efficient fertilizer N use would increase consumers' confidence in crop production sustainability and food security and hence would support an increased market access for the producers (Urso and Gilbertson, 2018). Management of the balance between fertilizer N losses and recovery in agriculture is also very critical since humanity has already exceeded the safety threshold of biogeochemical N cycling (Rockström et al., 2009). The need to reduce N losses and increase NUE is one of the top 10 global environmental priorities (UNEP, 2014), which is also a part of the sustainable development goals and aligns with the growing popularity of the cyclical economy awareness of industry and governments (MacArthur Foundation, 2019).
Optimization of agronomic N fertilizer application depends on agro-climatic, soil, crop, nutrient management, and economic variables (Snyder, 2017). In western Canada, farmers' agronomic nutrient planning includes making key decisions on fertilizer products, and timing, placement, and rates of application based on operational variables such as product availability, time, budget, labor, and equipment (Grant and Wu, 2008). The majority of N fertilizers used in western Canada are ammoniumbased (Statistics Canada, 2016). Although the majority of N fertilizers in western Canada are applied in the spring, fall application is sometimes preferred to distribute the workload and take advantage of lower fertilizer prices in the fall (Statistics Canada, 2016). However, fall applied fertilizer N is prone to losses, especially during the subsequent spring thaw before the crop can utilize it. After application, N fertilizers undergo hydrolysis and release ammonium (NH 4 + ), which is further oxidized to nitrate (NO 3 − ) by a microbial process called nitrification (Butterbach-Bahl et al., 2013;Bhanja et al., 2019;Shen et al., 2019a,b;Li et al., 2020). Ammonium (NH 4 + ) and nitrate (NO 3 − ) are the forms of N utilized by crops. However, NO 3 − is highly mobile through the soil and hence more prone to losses through runoff and leaching through the soil especially during spring snowmelt, which can cause surface water and groundwater pollution. NO 3 − can also be denitrified to produce potent GHG N 2 O and inert N 2 gases under saturated conditions (Butterbach-Bahl et al., 2013;Bhanja et al., 2019;Shen et al., 2019a,b;Li et al., 2020). These N losses can be minimized with a spring application that can optimize fertilizer NUE and reduce loss concerns. For instance, N 2 O emissions from prairie crop fields can be reduced by up to 30% depending on weather and soil type by avoiding spring thaw following a fall N fertilizer application (Dunmola et al., 2010;Glenn et al., 2012;Li et al., 2012Li et al., , 2016Maas et al., 2013).
Fertilizer application placement also plays an important role in crop N use efficiency and agronomic N losses. Although most of the N fertilizers in western Canada are applied in some forms of in-soil banding, broadcasting is still a dominant placement method in large areas of pastures or forages and in split or inseason applications (Grant and Wu, 2008). However, surface broadcasting can be very inefficient agronomically, economically, and environmentally since it can cause up to 50% more N losses through NH 3 volatilization compared to banding (Sheppard et al., 2010). Various enhanced efficiency products such as nitrification and urease inhibitors, and coated urea, have been developed with the objective of reducing N losses and to improve crop N use efficiency (Li et al., 2020). For instance, Environmentally Smart N (ESN, Nutrien) is a polymer-coated urea that is designed to slow down the rate of N release to better match the crop demand than conventional urea, which would improve crop N use and reduce N losses (Cahill et al., 2010;Gao et al., 2015). Selecting the right combinations of fertilizer products, rate, application timing, and placement, i.e., the 4R (Right Source @ Right Rate, Right Time, and Right Place) nutrient stewardship, can thus be an economically viable and environmentally sustainable strategy in western Canadian prairies (Malhi et al., 2001;Grant et al., 2002).
Effectiveness of 4R options for optimizing crop N use efficiency and minimizing N losses varies with variations in soils, landforms, and weather. Suitability assessment of a given combination of 4R for a given weather, landform, and soil condition is a prerequisite before a farmer makes a decision on adopting it. A provincial or prairie-wide numerical inventory of the fate of N applied in the agro-ecosystems would identify dominant regional N loss pathways and direct farmers to the most beneficial N management practices across various landforms, soils, and weather conditions (Dimitrov and Wang, 2019). The fertilizer industry can also use these assessments to identify opportunities of developing and commercializing enhanced efficiency fertilizer products. These estimates can also be scaled-up and displayed spatially across the province to show locations of "hot-spots, " which would support governments in targeting and designing appropriate incentive programs and policies. Currently, the effects of various 4R N management combinations on the fate of applied N in agro-ecosystems are being assessed based on site measurements, which are usually temporally and spatially discrete (e.g., Rawluk et al., 2001;Asgedom et al., 2014;Gao et al., 2015). The field measurements of different forms of N losses require linear temporal interpolations and spatial extrapolations for these assessments, which impart substantial uncertainties into these evaluations, since N transformation processes are highly nonlinear and often involve offsetting mechanisms (Flesch et al., 2018). Moreover, the site measurements are limited to fewer soil types, weather conditions, and management options for logistic reasons. Often different field projects are intended to evaluate different pathways of N transformations, which makes construction of a comprehensive budget for the fate of applied N into the agro-ecosystems very difficult. Consequently, a comprehensive inventory with temporally and spatially continuous estimates of crop N uses and fertilizer N loss under various agronomic N management does not exist to date for any of the prairie provinces.
Process-based mechanistic agro-ecosystem models can provide spatially and temporally scaled-up numerical estimates of different N pools under different N management scenarios. However, such a process model has to be built upon siteindependent algorithms from independent research, which can then be applied to various soil, weather, and agronomic management conditions without site-specific calibration of the model algorithms. The model should be able to reproduce a particular field condition from site-specific model inputs on soils, weather, and land, crop, and nutrient managements rather than tweaking the model codes for each scenario depending on the training dataset. The model outputs of different N pools have to be rigorously validated against site measurements under varying soils, weather, and agronomic management to evaluate model precision. Ecosys is such an ecosystem model, which successfully simulated soil-plant-atmosphere N continuum across different agro-ecosystems within and outside western Canada (Grant and Pattey, 1999Grant, 2001;Grant et al., 2006Grant et al., , 2016Metivier et al., 2009). Building upon those field-level validation studies, this study aimed at deploying the ecosys model spatially to generate provincial estimates of crop N use and agronomic N losses for Alberta crop fields under alternative N fertilizer management scenarios. In this study, we describe ecosys simulations to assess the effects of four selected N management scenarios: fall banded urea, fall banded ESN, spring banded urea, and spring banded ESN on yield, N uptake, and N losses across dryland barley fields of Alberta. We then corroborate the modeled results against available data and literature values to examine the adequacy of the simulated results in describing spatial distribution of the pathways of movement of applied N within Alberta agro-ecosystems. These selected scenarios serve as prototypes for establishing a scaling-up methodology, which numerically estimate the fate of fertilizer N under various N management scenarios.

Model Description
Ecosys is a process-based, hourly time-step, terrestrial, ecosystem model where transformations, transport, and exchanges of N within the modeled ecosystem are simulated in conjunction with those of carbon (C) and phosphorus (P) in a comprehensive modeling scheme, in which they are coupled with soil water, heat, and solute transport (Grant, 2001). Nitrogen transformation processes in ecosys are predominantly governed by coupled reduction-oxidation reactions, which result in microbial and root energy yields, decomposition, and growth, and hence drive N mineralization-immobilization, nitrification, and denitrification. Root and mycorrhizal N uptake occur through ion exchange, radial diffusion, and convection, which affects rubisco (ribulose-1,5-bisphosphate carboxylase/oxygenase) activation that drives modeled crop productivity, growth, and yields in ecosys. Rubisco activation in ecosys can also be affected by functions of water, temperature, and oxygen stresses, and availability of other nutrients such as P. Gaseous and aqueous N losses through NH 3 gas and dissolved organic and inorganic N are also modeled in ecosys. Nitrogen availability in the modeled soil solution is also affected by simulated adsorption and desorption of NH 4 + between soil solution and clay surfaces.
Nitrogen inputs to a modeled ecosystem in ecosys include various chemical fertilizer formulations, manure, organic amendments, atmospheric deposition, and biological N 2 . Formulation, timing, placement, and rates of a fertilizer application event are explicitly defined by model inputs. Granular NH 3 -based fertilizers (e.g., urea) undergo hydrolysis, which controls the rate of N release from fertilizer granules. The hydrolysis process in ecosys is a function of soil moisture and temperature that is mediated by microbial activity. The rate of hydrolysis is calculated from a specific rate constant multiplied by total heterotrophic microbial activity, the urea concentration relative to its Michaelis-Menten (MM) constant, and an Arrhenius function of soil temperature, and can be reduced by urease inhibition. The urease inhibition is calculated from another MM function of the aqueous concentration of total active heterotrophic activity that serves as a proxy of urease activity. For ESN, the specific rate constant is one-fourth of the rate constant for urea until 10% of the total applied N is hydrolyzed, after which the rate constants for both ESN and urea become the same. This algorithm simulates an initial lag in N release from ESN as opposed to urea, which approaches a sigmoidal N release response for ESN (Cahill et al., 2010). All of the above algorithms are parameterized from the kinetics and equilibria of complex biogeochemical and eco-physiological processes reported by independent research on ecosystem functioning within a broad scope of spatial scale. These algorithms in ecosys thus do not require calibration, training, or parameterizing for each unique space-time scenario. Instead, a modeled agro-ecosystem in ecosys is simulated from site-specific model inputs of weather, soil, and agronomic management data. A more detailed description of ecosys algorithms representing N transformation, transport, and exchange including all key equations, variable definitions, parameters, and references can be found in Grant and Pattey (2003), Grant et al. (2006), Metivier et al. (2009), andGrant et al. (2002).

Methodology and Model Inputs
The N transformation, transport, and exchange algorithms in ecosys were used to derive numerical estimates of effects of N fertilizer timing and products on crop N use and agronomic N losses in the western Canadian province of Alberta. Alberta is the fourth largest province of Canada, which extends between 49 • −60 • N and 110 • −114 • W occupying an area of 661,848 km² (Figure 1). The southern part of Alberta has a semi-arid climate (Köppen climate classification BSk) whereas central and northern Alberta experience humid continental climate (Köppen climate classification Dfb). A total of four parallel sets of model simulations were set up to simulate four selected N management scenarios: fall banded urea, fall banded ESN, spring banded urea, and spring banded ESN. Modeled outputs of N recovery in yield and N uptake of dryland/rain-fed barley; gaseous N losses as N 2 O, NH 3 , and N 2 ; and aqueous N losses through surface runoff and subsurface discharge from those barley fields were used to determine the effects of the selected N management scenarios on the fate of N in Alberta crop fields. Each of these four simulations had a total of 3,063 township scale (∼10 km × 10 km) spatially explicit grid cells that spread across the agricultural areas of Alberta (Figure 2; Table 1). Each of the township grid cells was divided into four landforms: top-, mid-, and foot-slopes, and depressional areas for simulations, depending on slope classification for different landforms (MacMillan and Pettapiece, 2000). Since mid-slope landforms comprise the largest proportion (40%) of the total arable lands in Alberta, we simulated the mid-slope land form elements in this study (Table 1). However, the other three landforms will eventually be simulated in future phases of this modeling project to generate more comprehensive scenarios of fate of N in prairie agroecosystems. Model inputs of the soil properties represented the key characteristics of four major agricultural soil groups of Alberta, i.e., Brown Chernozems (Aridic Borolls), Dark Brown Chernozems (Typic Borolls), Black Chernozems (Udic Borolls), and Dark Gray and Gray Chernozems and Luvisols (Boralfs and Mollic Cryoboralfs) (Figure 2; Table 1). The soil properties were derived from the most frequently occurring soil (modal or dominant) profile in each slope position for each township (Protz et al., 1968). The modal soil profiles were selected from the lands falling within Land Suitability Rating System's rating of 2-4 to represent Alberta's best arable lands for spring seeded small grains (Bock et al., 2018;LSRS, 2019).
Each grid cell was seeded with a barley plant functional type (PFT) ( Table 1). Barley was selected as the model crop for this study since it is the third largest crop grown in Canada, with over 90% of production located in the three Prairie Provinces of Alberta, Saskatchewan, and Manitoba (Statistics Canada, 2016). The barley PFT was built by customizing the wheat PFT in Grant et al. (2011) for eco-physiological adaptation for barley crops in Alberta. The ecological adaptation in the barley PFT were represented by adjusting the crop climate zone adaptation to the appropriate Köppen climate zones for different regions in Alberta. The physiological adaptation in the barley PFT was represented by raising the fractions of leaf proteins in rubisco and in mesophyll chlorophyll by 40%, which would simulate higher productivity and hence more rapid accumulation of grain biomasses in barley than in the wheat PFT in Grant et al. (2011). These adjustments to the barley PFT in ecosys with respect to the wheat PFT were made based on relative performance between ecosys simulated barley and wheat yields, biomass growth, and nitrogen uptake, which were rigorously tested by Grant et al. (2020) against Alberta field data. Typical soil, crop, and nutrient management practices and recommended N fertilizer rates across different regions of Alberta were used as model inputs (Figure 2; Table 1). The ratios of organic C to N, and N to P, in each soil layer were assumed as 10-1, which were typical to agricultural soils in Alberta. Sustained grain removal would create P limitation in the modeled crop. To eliminate P limitation, a phosphate fertilizer at a rate of 2.5 kg P ha −1 year −1 was applied to each grid cell in each simulation along with banded N fertilizers ( Table 1). The spinup runs ensured that the modeled ecosystems attained mass and energy balances to represent stable site conditions. All the simulations started with spin-ups from 2001 to 2010, which then extended to simulation runs from 2011 to 2015 using gridded, real-time, daily, weather data, i.e., maximum and minimum air temperature, incoming shortwave radiation, precipitation, wind speed, and relative humidity (Table 1) (ACIS, 2019). Since the ecosys model is an hourly time-step model, the model inputs of daily weather variables were first scaled down to hourly, to be implemented as hourly vertical model boundary conditions, to drive hourly model calculations. This temporal downscaling of the weather variables from daily to hourly was done internally inside the ecosys weather sub-model. The daily incoming shortwave radiation was downscaled to hourly values using a sinusoidal curve for radiation based on day length (Figure 3). The maximum and minimum daily air temperatures were used to drive a sinusoidal curve that calculated hourly temperature, so that the minimum temperature was reached at dawn and the maximum temperature was reached at 3 h after solar noon (Figure 3). Relative humidity was downscaled from daily to hourly by using a similar sinusoidal curve as the temperature (Figure 3). The daily precipitation was equally distributed to each hour in a day, and the average daily wind speed was used for each hour of the day as model upper boundary conditions. The temporal downscaling of air temperature and radiation would enable simulation of diurnal variations in N, C, heat, and water balance in the model. However, the temporal downscaling procedure of daily weather variables to hourly values, as described above, could still miss a sudden drop or rise of temperature from 1 h to another in a day, or a large precipitation event that occurred in some hours of a day. Lack of model inputs for these episodic events may affect the model's capability of accurately simulating sudden flushes of N losses (e.g., N 2 O). Model vertical boundary conditions, as described by the model inputs of air temperature and precipitation above, adequately represented drier growing seasons in the southern part of the province and long harsh winters, with relatively shorter and wetter growing seasons in the north. This is apparent in the monthly mean temperature and precipitation distribution across Alberta over the period of 5 simulation years from 2011 to 2015 (Figure 4).

Validation and Analyses of Modeled Outputs
Enhanced efficiency N fertilizers like ESN differ from conventional urea in their N release patterns from fertilizer granules. The granules of urea are coated with a polymer FIGURE 2 | (A) Spatial distribution of the soil groups comprising the agricultural zone of Alberta, which represents the study area for ecosys simulation. (B) Spatial distribution of top soil (0-15 cm) organic carbon contents, and (C) pH of the agricultural zone of Alberta that was simulated in this study. (D) Spatial distribution of typical fertilizer nitrogen application rates used in ecosys simulations. The white areas on the map were not modeled based on the selection criteria as described in section Methodology and Model Inputs. The thick dashed lines on the maps represent Alberta's provincial boundary, and the thin black lines within the maps demarcate soil group areas.
to slow down N release rates in ESN fertilizers. This slower release is designed to better match crop N demand to enhance crop N uptake and minimize agronomic N losses. Simulated N release patterns for ESN vs. urea were compared against laboratory data to examine how well ecosys simulated the observed differences between ESN and urea in their N release patterns. For this purpose, daily modeled N releases were accumulated and averaged as percentages of total applied N, for all four simulations, for all years, across all soils. The N release percentages of urea vs. ESN were then plotted against thermal time expressed as degree days. The thermal time was cumulative of daily accumulated modeled hourly soil temperatures above 0 • C at the depth of N fertilizer banding. The N release to thermal time relationship curves were then compared with similar curves constructed by data from a laboratory experiment.
The laboratory experiment was performed with a commercial top soil "Greensmix" of a sandy loam texture with a pH of 6.6 (Dowbenko, personal communication). Two separate sets of measurements were performed to account for N release from urea vs. ESN granules with a gradual increment and decline of temperature at the rate of about 5 • C/week. Each of these sets of measurements was performed under two moisture levels: 50 and 75% of the field capacity (Dowbenko, personal communication). During these laboratory experiments, the soil samples were maintained in sample pots within a growth chamber with designated temperature controls. While preparing the experiment pots, the fertilizer samples were evenly dispersed and covered with a 0.6-cm-thick soil layer. The two stated levels of moisture contents in the pots were checked daily and maintained throughout the experiment period (3-4 months). The rate of fertilizer release for each sample pot was measured once every week. For this purpose, the soil layer at the top of the fertilizer layer was removed very carefully to prevent any damage to the fertilizer granules. The loose soil particles were washed out of the granules by using a gentle stream of deionized water. The granules were then analyzed for N concentration by using colorimetry in aqueous solutions.
Seasonal and interannual variations in weather can significantly affect yield and crop N uptake. Modeled annual barley grain yields at typical grain moisture contents (13.5%) were averaged for each township over the simulation period (2011)(2012)(2013)(2014)(2015) to include effects of a range of weather conditions. The averaged modeled barley yields and annual grain N uptake were reported as rates in kg ha −1 year −1 . The rates of modeled barley yield and grain N uptake were then compared for the four scenarios to assess the effects of N fertilizer timing and products Model inputs of weather data Gridded daily weather data-maximum and minimum air temperature, incoming shortwave radiation, relative humidity, precipitation, and wind speed (ACIS, 2019) Implemented weather data into the model Daily incoming shortwave radiation (MJ m −2 d −1 ) downscaled to hourly radiation (W m −2 ) Daily maximum and minimum air temperatures ( • C) downscaled to hourly air temperature ( • C) Daily average relative humidity (%) downscaled to hourly relative humidity (%) Daily average wind speed measured at 10 m (km h −1 ) downscaled to hourly average wind speed (m s −1 ) Daily precipitation (mm day −1 ) were equally redistributed as hourly precipitation ( Organic carbon-to-nitrogen ratio 10:1 for each vertical soil layer in each grid cell, which is typical for Alberta agricultural soils Organic nitrogen-to-phosphorus ratio 10:1 for each vertical soil layer in each grid cell, which is typical for Alberta agricultural soils on the recovery of applied N fertilizer. These comparisons were performed spatially at township scales in maps and also by comparing soil group averages in bar charts. Modeled barley grain yields were validated against Agriculture Financial Services Corporation (AFSC) of Alberta data across agricultural areas of Alberta for the simulation period (2011)(2012)(2013)(2014)(2015). In Alberta, farmers using crop insurance have to report annual crop yields to be in compliance with the AFSC for crop insurance purposes (AFSC, 2019). The AFSC compiles and publishes the reported annual yields for each crop based on an area weighted averaging for a total of 22 risk zones across the agricultural areas of Alberta (AFSC, 2019). In this study, modeled barley grain yields were averaged for all modeled townships that fell within each of the AFSC agricultural risk zones, for all four N management scenarios, over the simulation period to facilitate comparison against AFSC compiled observed barley yields, averaged over the same time period (2011)(2012)(2013)(2014)(2015). The comparison between modeled and observed (AFSC data) barley yields would provide a measure of model accuracy in simulating geo-spatial variations in barley yields across Alberta over the simulation period. Model accuracy was evaluated by geo-spatial Pearson's correlation, slope, intercept, and root mean square for errors (RMSE) of linear regression of modeled vs. observed yields. This test would provide details on model accuracy and uncertainties in simulating regional and provincial scale crop yields and N uptake.
Modeled hourly outputs for different forms of agronomic N losses were accumulated annually as rates in kg N ha −1 year −1 and were also averaged for each grid cell in each simulation over the simulation period. These N losses were mapped and compared to facilitate township scale spatial comparisons among the four N management scenarios. The rates of different forms FIGURE 3 | Example representations of temporal downscaling of daily weather data to hourly weather variables to be implemented as hourly model inputs for vertical boundary conditions over two selected township grids using sinusoidal functions. IR sw , incoming shortwave radiation; T, air temperature; T min , daily minimum air temperature; T max , daily maximum air temperature; and RH, relative humidity. The thick dashed lines on the maps represent Alberta's provincial boundary, and the thin black lines within the maps demarcate soil group areas.
of N losses were also averaged and scaled up to the soil group levels to facilitate regional comparisons among different N management scenarios. While averaging by soil groups over the simulation period (2011)(2012)(2013)(2014)(2015), standard deviations were illustrated and coefficients of variations were reported to demonstrate the spatiotemporal variations of the modeled N recovery and losses due to variations in soils and weather.
Percent changes in modeled yields and in the key components of modeled N budget for a change in N fertilizer timing or product were listed to facilitate a comprehensive summary of the N management scenario analyses. Modeled trends, and magnitudes and ranges of modeled values for yields, grain N uptake, and various forms of N losses, were also compared with field observations and estimates from available published research for Alberta or similar prairie conditions.

Modeling Nitrogen Release From Urea vs. ESN
Modeled ESN and urea differed from each other in their average N release patterns, which was corroborated well by the observed laboratory results (Figure 5). The rate of modeled average N release from urea initially increased rapidly with increasing thermal time, after which it plateaued (Figure 5). Averaged modeled N release rate from ESN was slower than that of urea at lower thermal time after which it became very close to that of urea, approaching a sigmoidal N release pattern (Figure 5). Modeled N release rates for both fertilizer products for a given exposure to a thermal time varied [Coefficient of variation (CV) up to ±25%] due to variations in soils and weather across Alberta agricultural areas (spatial distribution not shown) (Figure 5). The observed laboratory results for urea showed a similar increase with increasing thermal time at a gradually decreasing rate as modeled (Figure 5). Observed N release rate for ESN followed an initial lag similar to the N release pattern as modeled (Figure 5).

Modeling the Effects of Agronomic Nitrogen Management on Nitrogen Recovery in Barley Grain
Effects of agronomic N management on the recovery of applied N were assessed by the variations of modeled estimates of barley grain yields and N uptake with variations in N application timing and products. Geo-spatial variations in modeled barley yields across Alberta agricultural areas during 2011-2015 corroborated well against observed (AFSC data) barley yields as indicated by a strong geo-spatial correlation between modeled and observed yields (Figure 6). However, a slope of 1.1 and an intercept of 538 kg ha −1 year −1 from a simple linear regression of modeled vs. observed barley yields, meant the modeled barley yields were larger than the observed yields (Figure 6). A smaller RMSE of 215 kg ha −1 year −1 from a simple linear regression of modeled vs. observed barley yields, however, showed lower model uncertainties in predicting geo-spatial variations in barley yields across Alberta during the simulation period (Figure 6). Modeled barley yields varied across the province by region, with the drier brown soils having the lowest average yields and the wetter black soils having the highest average yields (Figures 5, 7). Overall, the modeled barley yields remained mostly unaffected, either by a variation in the timing of application (fall vs. spring) or by a variation in N fertilizer products (urea vs. ESN) (Figure 7; Table 2). A change in the application timing from fall to spring resulted in only about 2% overall increase in modeled barley yields across Alberta during 2011-2015 ( Table 2). Contrary to the expectation that ESN would produce higher yields, modeled ESN application indicated no significant yield effect when compared to modeled urea application across the province (Table 2). However, there were localized effects of fertilizer timing and products on modeled barley yields, which was revealed at the township-scale spatial distribution of the modeled barley yields (Figure 7). For instance, there was about 20% reduction in modeled barley yields in some parts of the southeast dark graygray soil zone caused by a change from fall to spring application (Figure 7). These reductions were greater in ESN than in urea (Figure 7). In contrast, there were increases in modeled barley yields in some parts of the northwest dark gray-gray soil zone and in the southern dark brown zone, for a change from fall to spring application (Figure 7).
Modeled barley grain quality, as represented by grain N uptake (or content), showed a similar spatial pattern as the modeled grain yields, with brown soils having the lowest, and black soils having the highest average grain N uptake (Figure 8). Like modeled barley grain yields, grain N uptake did not show any regional or province-wide consistently discernible effects of N fertilizer timing or products (Figure 8; Table 2). Modeled grain N content was higher by about 2% in the spring application over the fall application and was down by only 1%, while ESN was applied instead of urea across the province (Figure 8; Table 2). However, the modeled barley grain N uptake also showed similar localized effects of N fertilizer timing and products, as did the modeled grain yields. Modeled grain N uptake declined by about 10-15% from fall to spring applications in parts of the southeast dark gray-gray soil zone and increased in parts of the northwest dark gray-gray and southern dark brown soil zones (Figure 8).

Modeling Effects of Agronomic Nitrogen Management on Nitrogen Losses
Variations in N application timing (fall vs. spring) and products (urea vs. ESN) had different effects on modeled N losses in the forms of N 2 O, N 2 , and NH 3 gases, and N losses in surface runoff and sub-surface discharge. Average modeled annual soil N 2 O emissions varied from 0.68 (CV ± 20%) to 1.88 (CV ± 47%) kg N ha −1 year −1 across the soil groups in all scenarios (Figure 9). Modeled average annual N 2 O emission varied among soil groups. Black and dark gray-gray soils had higher modeled N 2 O emissions than the brown and dark brown soils (Figure 9). Modeled N 2 O also varied substantially within each soil group area (Figure 9). Modeled N 2 O emissions were smaller in spring banding than in fall banding, irrespective of fertilizer products (Figure 9). On average, modeled N 2 O emissions from spring banding was 24% less than that from fall banding throughout dryland barley fields on mid-slope landforms across Alberta during 2011-2015 (Figure 9; Table 2). Annual reductions in N 2 O emissions with a change from fall to spring banding was almost double in black and dark gray-gray soils than those in brown and dark brown soils (Figure 9; Table 2). However, there was no discernible soil-group-wide or province-wide difference in modeled N 2 O emissions for variations in N products from urea to ESN (Figure 9; Table 2).
Complete denitrification simulated agronomically inconsequential N losses in the form of N 2 , averages of which ranged between 0.63 (CV ± 48%) and 1.20 (CV ± 56%) kg N ha −1 year −1 across the soil groups in all scenarios (Figure 10). A change from fall to spring banding reduced N 2 -N losses by about 32% across the province (Figure 10; Table 2). Variations in N products (urea vs. ESN) did not simulate any discernible change in N 2 -N losses (Figure 10; Table 2).
Modeled average volatilization of NH 3 ranged from a consumption of 0.19 (CV ± 150%) kg N ha −1 year −1 to an emission of 1.18 (CV ± 138%) kg N ha −1 year −1 across the soil groups in all scenarios (Figure 11). On an average, dark gray-gray soils consumed NH 3 from the air, and the other soils emitted NH 3 (Figure 11). Spring banding of both fertilizers showed a reduction of NH 3 emissions by about 67% across brown, dark brown, and black soil zones (Figure 11; Table 2). However, only a 5% reduction in NH 3 -N loss was simulated for a change from urea to ESN across these soils (Figure 11; Table 2). Effects of N fertilizer timing and product on NH 3 -N were not relevant for dark gray-gray soils since all the modeled fertilizer management scenarios simulated average consumptions of NH 3 (Figure 11; Table 2).
Modeled average dissolved organic and inorganic N in surface runoff and sub-surface discharge ranged from 0.32 (CV ± 90%) to 1.15 (CV ± 66%) kg N ha −1 year −1 across the soil groups in all fertilizer management scenarios (Figure 12). Nitrogen losses in surface runoff and sub-surface discharge varied regionally with black and dark gray-gray soils producing higher losses than the other two soils (Figure 12). However, there were large variations in modeled N losses in surface runoff and sub-surface discharge within each soil group area (Figure 12). Overall, N losses in surface runoff and sub-surface discharge were 37% less in spring banding than in fall banding across the province during 2011-2015 (Figure 12). Variations in N products from urea to ESN did not apparently affect N losses in surface runoff and sub-surface discharge (Figure 12; Table 2).

DISCUSSION
The geo-spatial variations of modeled dryland barley grain yields and N uptake across different soil groups and climates in Alberta during 2011-2015 matched reasonably well against observed variations in barley yields and N uptake reported across Alberta (Figures 6-8) (McKenzie et al., 2004;Anbessa and Juskiw, 2012;Perrott, 2016). However, modeled barley yields were about 25% larger than the AFSC reported observed yields during this period across Alberta (Figure 6). We provide the following explanations for such deviation. The observed yields published by AFSC included reduction in yields due to hail events, insect and pest damages, and lodging that were not simulated in the modeled scenarios (AFSC, 2019). Moreover, the AFSC data included reported yields from landforms that were not modeled in this study (e.g., top-and foot-slope positions and depressional areas) and from fields where lower than recommended N fertilizers may have been applied, which further contributed to the modeled vs. observed grain yield divergence. Model inputs for P fertilizer ( Table 1) would have alleviated any possible P limitation to the simulated barley yield, which could have also contributed to the larger modeled yields vs. observed AFSC data. Barley yields in top-slope positions could be suppressed by low moisture availability. Yields in footslope and depressional area could also be affected by excessive moisture and lodging. These phenomena were not accounted for in these simulations since these simulations included only welldrained, mid-slope landforms, which tend to be higher yielding. Although the modeled dryland barley grain yields were larger than the AFSC data, the modeled grain yields and grain N uptake were well within the ranges of long-term experimental data on barley grain yields (i.e., 4,300-6,900 kg ha −1 year −1 ) and N uptake (i.e., 81-131 kg N ha −1 year −1 ) across brown, dark brown, black, and gray soils of southern and central Alberta (Figures 7, 8) (McKenzie et al., 2004;Anbessa and Juskiw, 2012;Perrott, 2016).
Variations in barley grain yields and N uptake across soil groups were modeled predominantly by adequate simulation of moisture and N availability to the modeled crop (Figures 6-8). Brown soils in southern Alberta have the lowest soil organic matter levels and received the lowest fertilizer rates and the lowest mid-growing season (July) precipitation during the simulation period (Figures 2, 3). Low soil organic matter in drier soils provided less substrate for N mineralization in the model, which, along with low fertilizer inputs, caused low available N for modeled crop growth and uptake. Low moisture availability also simulated crop water stress and, hence, a further decline in modeled crop growth and N uptake. Lower N and moisture availability caused lower modeled barley grain yields and N uptake in brown soils compared to the other soils that received relatively higher moisture and N inputs and had higher N mineralization from higher organic matter (Figures 2, 3, 6-8).
Modeled spring banding produced about 5% higher barley grain yields than fall banding in the dark brown soils, with no apparent change in yields between fall and spring banding in other soils (Figure 7; Table 2). This is supported by Malhi et al. (1992), whose plot research found that spring banding of urea fertilizer produced about 8% higher barley grain yield than late fall banding across Alberta. Modeled ESN demonstrated an initial lag response in N release compared to urea, which was corroborated well by independent laboratory studies ( Figure 5) (Cahill et al., 2010). Despite the slower release, modeled ESN produced similar barley yields and N uptake to those under urea application (Figures 7, 8; Table 2). Gao et al. (2015) also found no significant change in spring wheat yields and grain N uptake between banded urea and ESN in two black chernozemic soils in Canadian prairies. Although initial slow release of ESN compared to urea was modeled and measured in laboratory studies (Figure 5) (Cahill et al., 2010), there was no field data available to validate the differences in modeled release pattern of N between urea and ESN. While the validation of modeled N release pattern of urea vs. ESN against the laboratory studies facilitated a comparative validation of modeled N release pattern between urea vs. ESN, it is still very important to corroborate the modeled N release pattern against field data where various factors like weather, soil temperature, moisture, land management, and crop uptake interact frequently. Validation of modeled N release against field data in future studies would further reduce the uncertainties in modeling the effects of conventional vs. controlled release N fertilizers and further improve our predictive capacity on the fate of these N management practices.
N fertilizer timing and products affected modeled barley yields and grain uptake differently based on soils and weather conditions. For instance, modeled grain yields and N uptake were lower in spring application than in fall application in some parts of southeast dark gray-gray soil group area (Figures 7,  8). The reduction in grain yields and N uptake from fall to spring were higher with ESN than with urea (Figures 7, 8). This was predominantly modeled from slower crop N availability governed by slower rates of urea hydrolysis as limited by lower microbial activity due to less soil organic C (Figure 2). Less over winter and spring precipitation over that area also increased aqueous concentration of microbial biomass, which further inhibited urease activity and hence rates of urea hydrolysis in spring banding scenarios that caused slower N release and 2 | Effects of nitrogen fertilizer application timing (fall vs. spring) and products (urea vs. ESN) on modeled annual dryland or rain-fed barley grain yields, modeled annual grain N uptake, and different forms of modeled annual nitrogen losses averaged over the simulation period (2011)(2012)(2013)(2014)(2015) for all the townships across mid-slope landforms of Alberta under simulated barley cultivation.

Change in application timing or product
Percent (%) change in modeled outputs for a change in application timing or product (a negative value means percent reduction and a positive value means percent increase for an associated change in nitrogen application timing or product) reduced modeled grain yields and N uptake compared to those in fall banding scenarios (Figures 4, 7, 8). Slower N release in ESN than in urea further hindered modeled grain yields and N uptake in spring banded ESN in those areas (Figures 5, 7, 8). Asgedom et al. (2014) also found that slower release of N from banded ESN reduced spring wheat and rapeseed yields progressively over 2 years in a black chernozemic soil in the prairies. The simulations in this study, however, did not include any blend of ESN and urea, which is becoming a farm practice in some areas of Alberta to overcome the early N needs and to reduce the higher cost of ESN.
On the contrary to the southeast dark gray-gray soil group areas, modeled barley grain yields and N uptake in spring banding were higher than the fall banding in some parts of northwest dark gray-gray and southern dark brown soil group areas (Figures 7, 8). These increases in grain yields and N uptake from fall to spring banding were modeled from rapid N release from spring application as facilitated by adequate microbial activity from higher soil organic C and higher urease activity due to lower aqueous microbial concentration caused by adequate over winter and spring precipitation in those areas (Figures 2, 3, 7, 8).
Modeled N 2 O-N losses were well within the range of estimates (i.e., 0.1-3.0 kg N ha −1 year −1 ) based on linear temporal interpolations of measured N 2 O fluxes in no-till fields across different soils and climates of Alberta under spring wheat and barley cultivation (Figure 9) (Lemke et al., 1998(Lemke et al., , 1999Rochette et al., 2008;Soon et al., 2011). Higher N 2 O emissions from black and dark gray-gray soils compared to brown and dark brown soils were simulated predominantly from larger NO 3 − accumulation from greater N inputs and N mineralization in wetter soils that enhanced the denitrification process in saturated soils (Figure 9). Like the modeled trend, increased N 2 O-N losses from fertilized crop fields were also reported for Alberta and globally with the increase in available substrate (NO 3 − ) for denitrification (Rochette et al., 2008;Shcherbak et al., 2014;Chai, 2017). Estimates from periodic field measurements of soil N 2 O fluxes showed larger N 2 O fluxes in black and gray soils than in brown and dark brown soils, which corroborates very well with the modeled geo-spatial distribution of N 2 O emissions from dryland barley fields on mid-slope landforms of Alberta across different soils and climates (Lemke et al., 1999;Rochette et al., 2008). Although modeled N release was slower in ESN than in urea, N 2 O emissions did not differ markedly in ESN simulations than in urea simulations for both application timing (Figures 5, 9; Table 2). Annual estimates based on periodic field measurements of N 2 O fluxes also showed that for recommended application rates of N fertilizer, ESN did not show any significant reduction in N 2 O-N losses from spring wheat, barley, and canola fields across Alberta and other prairie provinces as compared to the conventional urea (Li et al., 2012(Li et al., , 2016Gao et al., 2015). However, estimated annual N 2 O-N losses were 15-25% larger in urea than in ESN when application rates were 1.5 times higher than the recommended rates, or when there was considerable seeding delay in the spring that created excessive NO 3 − accumulation from the spring banded N fertilizer (Li et al., 2012(Li et al., , 2016Gao et al., 2015). Annual modeled N 2 O emissions were predominantly contributed by large flushes of N 2 O fluxes during spring thaw when inadequate O 2 supply forced modeled microbes to reduce NO 3 − as alternate electron acceptors. Consequently, larger soil N 2 O fluxes (by up to 30%) from fall banding than from spring banding were modeled for both urea and ESN applications (Figure 9; Table 2). Estimates based on field measurements also showed that up to 30% of the total annual N 2 O emissions could be contributed by the large flushes of soil N 2 O fluxes during winter and spring thaw from fertilized crop fields across Canadian prairies (Dunmola  Glenn et al., 2012;Maas et al., 2013). Some estimates based on field data showed that fall banding could cause up to 50% greater N 2 O-N losses than spring banded urea and ESN from wheat-barley-canola systems in dark gray soils of Alberta (Soon et al., 2011). From a plot-based periodic flux measurement study, Hao et al. (2001) estimated about 60% greater N 2 O emissions with fall broadcasting of ammonium nitrate fertilizer followed by tillage than with similar spring applications. Modeled N 2 O emissions varied spatially with NO 3 − availability and degree of soil saturation (Figure 9). Higher NO 3 − availability in wetter soils produced large modeled N 2 O emissions of up to 4 kg N ha −1 year −1 in some parts of northern and central Alberta (Figure 9). Wetter soils were simulated in those areas during the simulation period (2011-2015) since water inputs through precipitation were greater than vertical water losses through evaporation in the modeled landscapes representing those areas. Based on periodically measured N 2 O fluxes, Nyborg et al. (1997) estimated that the N 2 O-N losses can be up to 3.5 kg N ha −1 within only an 11-day period of spring thaw, further indicating potentials of greater N 2 O losses due to NO 3 − accumulation during wetter spring. Current field data-based N 2 O-N loss estimates are predominantly based on linear interpolations of periodic flux measurements, which may miss the episodic flushes of N 2 O and, hence, may underestimate annual N 2 O losses. Flesch et al. (2018) further emphasizes the importance of higher temporal resolution measurements in estimating annual N 2 O-N losses while measuring 6-hourly N 2 O emissions from barley fields on gray luvisolic soils of central Alberta by using micrometeorological techniques. This study revealed that N 2 O-N emissions could accumulate up to 5.3 kg N ha −1 from barley fields over only a month during spring thaw, which were larger than the current annual estimates for western Canadian prairies (Flesch et al., 2018). However, maintaining long-term high-resolution N 2 O measurements and replicating them spatially are highly demanding of time, technology, and cost. Some of these limitations could be overcome by supplementing current estimates with the modeled annual N 2 O estimates in this study that were derived by accumulating hourly calculated outputs (section Validation and Analyses of Modeled Outputs). Di-nitrogen (N 2 ) emission is usually overlooked but can be a very significant form of N loss from agro-ecosystems that receive large N inputs and are under prolonged saturation (Zistl-Schlingmann et al., 2019). Modeled N 2 -N emissions were simulated predominantly from complete denitrification under saturated soil conditions during spring thaw, which eventually contributed to greater modeled N 2 -N losses in fall than spring banded urea and ESN (Figure 10; Table 2). Delayed release in N from ESN than urea did not produce a lower N 2 emission (Figures 5,  10; Table 2). Modeled N 2 -N was not large enough to be agronomically and economically significant (Figure 10; Table 2).
Modeled NH 3 emissions were higher in fall vs. spring banded urea and ESN in brown, dark brown, and black soils across Alberta (Figure 11; Table 2). The higher NH 3 volatilization in fall was simulated from higher availability of NH 4 + for volatilization in fall banding simulations resulting from over winter and early spring hydrolyses. Generally moist soils during and shortly after spring application increased NH 4 + solubility in the model and, hence, also reduced modeled NH 3 volatilization from spring banding than in fall banding. Reduced rate of urea hydrolysis in ESN caused about 5-8% reduction in NH 3 volatilization compared to urea across dark brown and black soils in Alberta (Figure 11; Table 2) (Rawluk et al., 2001). Overall, modeled NH 3 emissions across southern and central Alberta were lower than the estimate of about 2.5 kg NH 3 -N ha −1 year −1 for western Canadian wheat fields in an emission factor based monthly NH 3 emission modeling (Sheppard et al., 2010). Simulated NH 3 emissions were also affected by soil pH. For instance, black and dark brown soils that had pH between 7 and 8 shifted the modeled chemical equilibrium in such a way that NH 4 + solubility was reduced, which ultimately enhanced NH 3 volatilization (Figures 2, 11) (Sommer and Ersbøll, 1996;Bouwman et al., 2002;Havlin et al., 2013;Grant et al., 2016). In contrast, most of the dark gray-gray soils had pH under 6, which enhanced modeled NH 4 + solubility and hence caused net NH 3 consumption instead of emission (Figures 2, 11) (Sommer and Ersbøll, 1996;Bouwman et al., 2002;Havlin et al., 2013;Grant et al., 2016). Nitrogen losses in surface runoff and sub-surface discharge from crop fields in western Canadian prairies are primarily snowmelt driven (Casson et al., 2008;Tiessen et al., 2010). Dissolved inorganic N from fall banded urea and ESN were transported along with the surface runoff and sub-surface discharge waters during snowmelt and spring thaw, which produced larger modeled N-runoff in fall vs. spring banding across Alberta over the simulation period (Figure 12; Table 2). Less runoff, combined with lower NO 3 − accumulation from lower N inputs and less N mineralization, caused lower modeled N-runoff fluxes in brown soils than the other soils (Figures 2,  12) (Casson et al., 2008). Modeled N-runoff fluxes were overall smaller compared to the estimates of 7.2-11.7 kg N-runoff losses ha −1 year −1 from eastern Canadian crop fields that received higher precipitation and N inputs (De Jong et al., 2009). However, in western Canada, N loss through NO 3 − leaching was estimated to vary from 1.5 to 4.5 kg N ha −1 year −1 for various N application rates in various rotations in long-term (over 30 years) research plots (Campbell et al., 1994(Campbell et al., , 2006. These leaching losses in western Canada were negligible in continuously cropped and fall seeded fields but were very large in excessively fertilized fields and fields under summer fallows, which favored NO 3 − accumulation and soil moisture buildup that are precursors of NO 3 − leaching (Campbell et al., 1994(Campbell et al., , 2006. Modeled NO 3 − losses through subsurface drainage in this study could be a proxy of NO 3 − leaching since these NO 3 − were transported out of the modeled root zone along with subsurface drainage through lateral model boundaries. Assuming modeled NO 3 − loss through sub-surface drainage as a proxy of N leaching, modeled losses of dissolved organic and inorganic N from continuously cropped fields through surface runoff and subsurface discharge in this study can be considered very conservative when compared to the aqueous N loss estimates through leaching from long-term field studies in western Canada (Figure 12) (Campbell et al., 1994(Campbell et al., , 2006. Although the N loss estimates in this modeling are within the range of most field-plot based estimates, the modeled estimates can be considered conservative since the simulation did not include foot-slope or depressional areas where N losses can be 2-to 3-fold higher than those from mid-slope landforms (Izaurralde et al., 2004). Extending these simulations to remaining landforms in interconnected transects of top-, mid-, and foot-slopes, and depressional areas would thus provide more comprehensive estimates of agronomic N losses and crop N use from various N management scenarios. Reproducing these simulations for highly fertilized crop (e.g., canola) fields, irrigated lands, and for agronomic management such as pulses in rotations, and residue retention, would also include simulations of extreme N losses and, hence, would provide better approximations of the fate of N from various agronomic N management. Besides, the simulations in this study were performed over large spatial extents of each grid cell (∼10 km × 10 km) using only four selected N management scenarios, extrapolated weather data, and dominant soil properties. Given the variabilities in soils, weather, crop, land use, and management practices within each grid cell of this size, the soil, weather, and N management practices used as inputs for these simulations may not always adequately represent the conditions in a field or a farm within a grid. So, the estimates from this modeling study can only be a first approximation of crop N recovery and agronomic N losses at a field or a farm scale. However, finer resolution modeling can be performed at a field or farm scale by selecting locations of interests from the modeled landscape and providing inputs to the model for soil, weather, crop, and management practices to adequately represent a field or a farm within the simulations (Table 1).

CONCLUSIONS
The process-based modeling analyses in this study indicated that the spring application of fertilizer could be an optimal N fertilizer application timing for Alberta farmers in reducing N losses while not compromising agronomic and economic returns in dryland barley cultivation across mid-slope landforms. Effectiveness of the spring banding in optimizing benefits and minimizing N losses, however, would be dependent on variations in soils, climate, and rates of N inputs. The modeled results, however, did not show discernible differences in barley N use or agronomic N losses from Alberta barley fields on the mid-slope landscapes due to a difference in N fertilizer products between urea and ESN.
Resilience of any agronomic N management option in terms of long-term sustainability and profitability is a key to successful farming operations. The desire to maximize production rather than optimize it may end up with N fertilizer application rates beyond economic profitability or environmental sustainability. This study opens up windows of opportunities for assessing the potential impacts of increasing N fertilizer application rates on agronomic N loss and crop yields and N uptake in Alberta. Such a study can also be used in assessing topographic influence on variable N fertilizer rates in precision farming to optimize crop productivity and minimize agricultural N losses. Application timing and placement, which are considered less efficient and more prone to losses, such as fall application and surface broadcast, are sometimes preferred for operational reasons. "What if " scenario analyses based on such modeling would provide the farmers with options of working with different enhanced efficiency N fertilizer products such as coated urea, urease inhibitors, nitrification inhibitors, or any combination of these technologies to reduce N losses and optimize production while operating within their operational limitations in terms of timing and placement.
This modeling approach can be used to identify "hot spots" or sensitive areas that are more prone to N losses. Policy makers can use this information to formulate applicable and sustainable policies and to devise incentive plans for promoting environmental stewardship in farming operations. The fertilizer industry can get valuable first-hand forecast to formulate and commercialize suitable products for profitable and sustainable agri-business. This study pioneered a methodology to assess the suitability of 4R nutrient stewardship options for sustainable crop production across a broad area of about 21 million ha. The simulations are also scalable to regional, federal, continental, and global scales by feeding the model with soil, climate, and management data appropriate to the scales (Table 1). Therefore, this study has important practical application, replicability, and validity in contributing to the existing knowledge pool of agricultural nutrient management science.

DATA AVAILABILITY STATEMENT
The original sources and links of all data used as model inputs and for model validation are appropriately cited in the article.

AUTHOR CONTRIBUTIONS
SM led the study and contributed to the formulation of the original idea, extracting data, designing modeling experiment, performing simulations, analyzing data, and writing the manuscript. DS, DH, LK, and DI contributed to the design and implementation of data extraction procedure for the modeling experiment and to manuscript writing. TG and JI contributed to the design of modeled output visualizations and manuscript writing. RG is the original developer of the model ecosys used in this study and also contributed to the modeling experiment and manuscript writing. All authors contributed to the article and approved the submitted version.

FUNDING
Funding for this project was provided by Growing Forward 2 Research and Innovation Internal Initiatives of federalprovincial partnership between Canada and Alberta. Computing support for this project was provided by Compute Canada WestGrid. Organization and analyses of model inputs and modeled outputs were performed by using R 3.6.3 packages of dplyr and tidyr. All the graphs and the maps (except Figure 1) were created by using R 3.6.3 packages of ggplot2 and rgdal.