A model for mechanistic and system assessments of biochar effects on soils and crops and trade‐offs

We developed a biochar model within the Agricultural Production Systems sIMulator (APSIM) software that integrates biochar knowledge and enables simulation of biochar effects within cropping systems. The model has algorithms that mechanistically connect biochar to soil organic carbon (SOC), soil water, bulk density (BD), pH, cation exchange capacity, and organic and mineral nitrogen. Soil moisture (SW)–temperature–nitrogen limitations on the rate of biochar decomposition were included as well as biochar‐induced priming effect on SOC mineralization. The model has 10 parameters that capture the diversity of biochar types, 15 parameters that address biochar‐soil interactions and 4 constants. The range of values and their sensitivity is reported. The biochar model was connected to APSIM's maize and wheat crop models to investigate long‐term (30 years) biochar effects on US maize and Australia wheat in various soils. Results from this sensitivity analysis showed that the effect of biochar was the largest in a sandy soil (Australian wheat) and the smallest in clay loam soil (US maize). On average across cropping systems and soils the order of sensitivity and the magnitude of the response of biochar to various soil‐plant processes was (from high to low): SOC (11% to 86%) > N2O emissions (−10% to 43%43%) > plant available water content (0.6% to 12.9%) > BD (−6.5% to −1.7%) > pH (−0.8% to 6.3%) > net N mineralization (−19% to 10%) > CO2 emissions (−2.0% to 4.3%) > water filled pore space (−3.7% to 3.4%) > grain yield (−3.3% to 1.8%) > biomass (−1.6% to 1.4%). Our analysis showed that biochar has a larger impact on environmental outcomes rather than agricultural production. The mechanistic model has the potential to optimize biochar application strategies to enhance environmental and agronomic outcomes but more work is needed to fill knowledge gaps identified in this work.


Introduction
Driven by the potential synergy of simultaneously sequestering carbon, enhancing soil quality and producing renewable energy (Laird, 2008), research on the pyrolysis-biochar platform is growing rapidly and various initially small scale commercialization efforts are underway. Pyrolysis, the heating of biomass in the absence of oxygen yields bio-oil and syngas, which are bioenergy raw materials, and a solid biochar residue, which has potential use as a soil amendment (Laird et al., 2009). The physical and chemical properties of biochar, however, vary substantially depending on the nature of the biomass feedstock and the thermochemical processing conditions . Such variation influences the stability of biochar carbon in soil environments, hence carbon sequestration potential, as well as soil quality and agronomic responses to soil biochar amendments. Furthermore, complex interactions between soil type, biochar type, climate, and management practices can substantially influence crop response to biochar applications. Thus a systems level understanding of soil-crop-climate-biochar interactions is needed to predict a priori whether the application of a particular biochar on a particular soil will result in positive, negative or neutral agronomic and environmental responses.
To date, the majority of biochar studies have been carried out as short-term lab incubation or greenhouse pot experiments focusing on soil quality responses but neglecting crop and climate interactions. Results from these studies have demonstrated both desirable and nondesirable effects of various biochar types on soil properties and processes and plant growth. Biochar has been found to increase plant available water content (PAWC), which is promising, given that soil water is often the main limiting factor of crop production (Basso et al., 2013). Most studies report a decrease in soil bulk density (BD) as a result of biochar application, which is another desirable characteristic as low BD enhances soil aeration and ultimately the ability of plant roots to explore and fully use the soil volume (Laird et al., 2010). The increase of soil organic carbon (SOC) as a result of biochar application is well documented in the literature (Brockhoff et al., 2010;Haefelea et al., 2011;Rogovska et al., 2014). However, biochar carbon is qualitatively different from biogenic SOC and thus is not expected to have the same influence on nutrient cycling and other soil processes as biogenic SOC. The effect of biochar additions on soil carbon and nutrient cycling is still under investigation and the experimental findings are inconclusive, as both increased and decreased greenhouse gas (GHG) emissions and N mineralization rates have been reported (Spokas & Reicosky, 2009;Clough et al., 2010;Bruun et al., 2011a;Strubel et al., 2011;Ameloot et al., 2013;Cayuela et al., 2014Cayuela et al., , 2015Gai et al., 2014;Case et al., 2015). Direct effects of biochar on soil cation exchange capacity (CEC), pH, PAWC, BD, and nutrient adsorption may have cascading influences on carbon and nutrient cycling. Furthermore, there is evidence that biochar amendments induce neutral, positive and negative priming effects on the turnover of SOC and the mineralization of crop residues (Cross & Sohi, 2011;Zimmerman et al., 2011;Jones et al., 2012;Spokas et al., 2012;Wang et al., 2015).
In contrast to incubation or pot studies, field trials examining the effect of biochar on intensively managed and fertile agricultural soils are currently rare and the effects on soil properties and processes are commonly less pronounced or not evident in these studies (Major et al., 2010;Haefelea et al., 2011;Guerena et al., 2013;Jones et al., 2012;Rogovska et al., 2014). Using a metaanalysis, Jeffery et al. (2011) found that biochar, on average, increased crop yield by 15% in pot studies and by 5% in field studies. The variation in yield response to biochar was high and included some negative responses, which reflects the complexity that exists in agroecosystems.
The complexity in biochar research increases greatly when different types of biochar (feedstock source x production method 9 peak pyrolysis temperature) are considered, which result in substantial variation in physical and chemical properties. Thus, determining where and what type of biochar to apply to optimize positive agronomic and environmental outcomes is very challenging. Yet this is essential for the strategic targeting of biochar applications and will add economic value to the pyrolysis-biochar platform. To address this, biochar experimental findings must be integrated and re-analyzed under a common framework. Meta-analysis studies (Jeffery et al. 2011;Biederman & Harpole, 2013;Crane-Droesch et al., 2013;Cayuela et al., 2014;Liu et al., 2015;Wang et al., 2015) are the first steps in the integration of biochar information, but they lack the power of prediction that a mechanistic approach can offer. Our effort is to link the individual biochar effects in a way that incorporates the complex interactions and feedback mechanisms. Jeffery et al. (2015) argued that biochar research needs a systems framework to allow effective comparisons between experiments, across scales and locations, to guide policy and recommendations concerning biochar application to soil. We hypothesize that agricultural systems models, which have a long history and have made significant contributions in agricultural research, can be used to quantify the complexities involved in the biochar-soil-plant-climate-management systems, and that doing so will help identify critical knowledge gaps, stimulate new strategic biochar research and accelerate commercialization.
At present, modeling of biochar is still in its infancy, and published efforts have focused on a particular aspect of the system without addressing concurrent effects of biochar on soil water, nitrogen, carbon and their interactions (Woolf et al., 2010;Foereid et al., 2011;Woolf & Lehmann, 2012;Lychuk et al., 2015;Lim et al., 2016). Our goal is to advance biochar modeling by accounting for concurrent biochar effects on soil carbon, water, nitrogen, and other soil properties and their dynamic interactions. Our specific objectives were: (i) develop a conceptual mechanistic biochar model by linking biochar effects to the soil-plant-atmosphere system within the Agricultural Production Systems sIMulator (APSIM) software platform (Holzworth et al., 2014); (ii) provide a range of parameter values that address diversity among different biochar types and conduct a sensitivity analysis of those parameters; (iii) use the model to investigate long-term biochar effects on US maize and Australia wheat cropping systems; (iv) conduct a preliminary evaluation of the modeling approach by comparing predicted and measured soil and crop responses to biochar applications, and (v) identify knowledge gaps and biochar research priorities.

Biochar model development
We synthesized biochar information and programmed the biochar model within APSIM (Holzworth et al., 2014). APSIM is a worldwide leading agricultural systems modeling platform, which is frequently updated, freely available, and contains more than 30 crop models and soil/environmental models along with flexible specification of farm management. Hence it can assist investigation of biochar effects on different cropping systems growing in different environments (soil 9 climate combinations). Because we developed and run the biochar model within APSIM we ensured that the biochar model is compatible with the soil process models already developed within APSIM.

Background on APSIM
Agricultural Production Systems sIMulator is a field scale model that operates on a daily time step (www. apsim.info). Its modular design allows users to select a crop or a combination of crops in rotation, cultivars, soils, climates, and management practices. In turn the model provides outputs for numerous soil-plant variables; including crop growth processes, soil water, soil temperature, nitrogen and carbon cycling, GHG emissions and residue dynamics (Probert et al., 1998;Thorburn et al., 2001Thorburn et al., , 2010Keating et al., 2003;Holzworth et al., 2014). The major processes simulated by the model and their feedback mechanisms are shown in Fig. 1. The model uses a multi-layer soil profile and soil processes are simulated in each layer separately. SOC is divided into three pools, each with a different decomposition rate constant; including BIOM (microbial biomass pool with a fast decomposition rate), HUM (humic pool with an intermediate decomposition rate) and INERT (inert pool which does not decompose). Fresh organic matter (FOM) is a separate pool, which contains roots from the previous crop and crop residues if they have been incorporated into the soil by tillage. The FOM is also divided into three sub-pools reflecting fast (FPOOL1), slow (FPOOL2) and medium (FPOOL3) decomposition rates. Both FOM and BIOM are fast decomposing pools, but they have different carbon to nitrogen (C : N) ratios as they represent different organic materials. The C : N ratio of BIOM is fixed at 8, while the C : N ratio of FOM varies depending on the source of residue/roots.   underlined theory, the equations and the parameter values developed for the biochar model (Table 1). The development of each biochar parameter was accompanied by a sensitivity analysis.

Biochar decomposition
The starting point in the development of the biochar model was to define a new biochar C pool for APSIM that contains both labile and recalcitrant components as commonly described in the literature (Woolf et al., 2010). The daily rate of biochar decomposition was modeled using a modified double exponential decay function (Archontoulis & Miguez, 2015) that accounts for water, temperature and N limitations: where, dlt BC is the amount of biochar C that decomposes every day (kg C ha À1 day À1 ), BC o is the amount of biochar applied (kg dm ha À1 ), f loss is the fraction of biochar that is lost during application (0 to 1), f carbon is the fraction of C in biochar (0 to 1), f labile is the fraction of the labile biochar pool (0 to 1), MRT 1 and MRT 2 are the mean residence times for the labile and recalcitrant biochar pools, respectively, and w f , t f and n f are SW, temperature and N modifiers (0 to 1) for the rate of biochar decomposition. Given the lack of specific experimental data, we used modifiers similar to that used by APSIM for other organic materials including soil organic matter, residue and manure (Probert et al., 1998Thorburn et al., 2001Thorburn et al., , 2010. The equations and the parameters are provided below. Table 1 provides literature values for f carbon , f labile , MRT 1 and MRT 2 . Their values vary substantially among biochar types in relation to feedstock source and production method, and are therefore used in the model to characterize biochar diversity. These parameters are meaningful, measurable and user defined. The f carbon is determined by thermal combustion and may vary from 0.1 to 0.9 (Spokas & Reicosky, 2009;Atkinson et al., 2010;Brewer et al., 2011;Lehmann et al., 2011;Gai et al., 2014). The f labile is determined by proximate analysis {[(total C À fixed C)/total C] using ASTM standard method D 1762-84, 2007} and typically varies from 0.03 to 0.30 (Woolf et al., 2010;Wang et al., 2015). In general, f labile decreases with increasing pyrolysis temperature (Bruun et al., 2011b), but is also influenced by feedstock properties such as lignin content, and the thermochemical process. The MRT 1 is based on incubation studies and varies from 0.3 to 25 years (Woolf et al., 2010;Wang et al., 2015). The MRT 2 can be estimated based on H : C elemental ratio determined by thermal combustion analysis and varies from 50 to 50 000 years (Woolf et al., 2010;Woolf & Lehmann, 2012). The parameter f loss is used to correct for losses that may occur during biochar application in the field.
Equation (1) can take very different shapes depending on the parameter values used and the prevailing soilweather conditions (Fig. 2). A sensitivity analysis of BC o , f carbon , f labile, MRT 1 and MRT 2 parameters showed that the MRT 2 has substantially less impact on the rate of biochar decomposition than the other parameters (Fig. 3a). The above parameters set the potential decomposition rate for biochar (Fig. 4). In practice, the prevailing soil-crop-weather conditions determine the actual decomposition rate, which is invariably lower. For example, Fig. 4 shows that the actual decomposition rate was 15% of the potential on an annual basis, with large day-to-day fluctuations. This is because soil water, temperature and N are rarely at optimum levels in the field; by contrast with laboratory incubations, which are commonly run under optimum or near optimum conditions for biochar decomposition. As it is difficult to maintain optimum levels in the long term, this explains why the reported MRT 1 values in the literature declines as the experimentation duration increases (Wang et al., 2015).
To account for limitations on the decomposition rate of biochar we used APSIM's temperature and water modifiers on the rate of SOC decomposition (Probert et al., 1998). For the N modifier (n f ), we developed a new function similar to that currently used by APSIM for residue decomposition: where the cnrf_bc parameter determines the slope of the curve, Soil BCL is the labile biochar in the soil, N avail is the sum of inorganic N in the soil and N released during biochar decomposition (see below), and the opt_bc parameter determines the change from immobilization to mineralization. The default cnrf_bc and opt_bc values were assumed similar to APSIM residue (Table 1), but more research is needed to verify the values of these constants. The sensitivity of the N modifier to these parameters is presented in Fig. 3c. The term Soil BCL /N avail changes dynamically in relation to the amount of biochar in the soil and the soil mineral N available for biochar decomposition and in general this is highly variable (Fig. 4). The modifiers w f , t f and n f , though generic, are biochar specific because different biochar types affect the SW and N availability differently.

Partitioning of the decomposed biochar carbon
The gradual degradation of biochar (Eqn 1) emits some biochar C to the atmosphere, and transfers biochar C to the microbial biomass (BIOM) and humic (HUM) pools (Kramer et al., 2004;Brodowski et al., 2007;Kuzyakov et al., 2009;Woolf & Lehmann, 2012). We modeled the partitioning of decomposed biochar C by defining an efficiency coefficient (ef_bc) that represents the fraction of decomposed biochar C that is retained in the soil and coefficients to distribute this retained C between the BIOM and HUM pools ( Fig. 1):  Daily rate of N mineralized or immobilized during biochar decomposition (kg N ha À1 day À1 ) Soil BC Labile and recalcitrant biochar remaining in the soil (kg C ha À1 ) Soil BCL Labile biochar remaining in the soil (kg C ha À1 ) w f, t f , n f Water, temperature and nitrogen modifiers on the rate of biochar C decomposition ( †FOM: fresh organic matter pool (mainly roots and residues). ‡BIOM: microbial soil organic matter pool (the labile portion of the soil organic matter). §HUM: humic soil organic matter pool (the bulk of the soil organic matter). ¶LL: soil water drained lower limit or permanent wilting point. kDUL: soil water drained upper limit or field capacity. **BD: soil bulk density. where: where, dlt_bc_CO 2 is the daily loss of biochar C to the atmosphere (kg C ha À1 day À1 ), dlt_bc_biom and dlt_bc _hum are the daily flux of biochar C (kg C ha À1 day À1 ) to BIOM and HUM, respectively, and fr_bcbiom is the fraction of retained biochar C that ends up in BIOM. The default value for ef_bc (Table 1) was assumed equal to the efficiency defined for humus decomposition in APSIM (Probert et al., 1998) as has been assumed in other biochar models (Woolf & Lehmann, 2012). The default value for fr_bcbiom was set so a small portion (0.05 9 0.4 = 2%; Table 1) of the decomposed biochar C was partitioned to BIOM. This assumption is reasonable given experimental findings that the size of the microbial pool increases following biochar application (Clough et al., 2010;Bruun et al., 2012;Zhang et al., 2014;Liu et al., 2015). The parameters described in this section have either not been experimentally measured or have been measured for only a limited range of biochar-soil-weather conditions (e.g. using isotopically labeled biochars, Zimmerman et al., 2011). Figure 3b shows the sensitivity of the rate of CO 2 evolution to BC o (amount applied) and ef_bc (C retention; Table 1).

Biochar priming effects
In addition to the changes in the size of the SOC pools (BIOM and HUM) imposed by biochar application and decomposition, some biochar types can alter the turnover rates of the biogenic SOC pools (Fig. 1). This is known as the priming effect. Woolf & Lehmann (2012) investigated priming effects on the labile FOM, while Cross & Sohi (2011), Zimmerman et al. (2011 and Wang et al. (2015) studied priming effects on SOC. Using this information, we developed functions in the biochar model to address both positive (accelerates decomposition rate) and negative (decelerates decomposition rates) priming effects. For the FOM pool, we defined a positive priming coefficient (P FOM ) as a function of Soil BC similar to Woolf & Lehmann (2012): where, X bc is the new dynamic decomposition rate (day À1 ) for each of the three APSIM FOM pools, X is APSIM's default decomposition rate of each FOM pool (day À1 ), P FOM is the positive priming coefficient (m 2 kg À1 C) that is assumed to be the same for all the FOM pools, Soil BC is the biochar remaining in the soil (kg C ha À1 ), and 10 000 is used for unit conversion. Biochar negative priming effects on FOM were modeled by increasing the C transfer from FOM to the more stable HUM pool and simultaneously decreasing the fraction of FOM C that is respired; in other words the internal cycling of carbon was increased similarly to Woolf & Lehmann (2012): where, ef_fombc is the modified ef_fom parameter that determines the amount of decomposed C retained in the system (0 to 1 value), fr_fom_biombc is the modified fr_-fom_biom parameter that determines the C flow from the FOM to BIOM pool (0 to 1) (note: the remaining C is allocated to HUM), and the parameters P e and P f are the negative priming coefficients (m 2 kg À1 C). Biochar induced positive and negative priming effects on SOC pools (BIOM and HUM) were computed in a similar way. All of the priming coefficients can be derived from lab incubation studies such as that of Kuzyakov et al. (2009) where isotopically labeled biochar was used to distinguish between CO 2 emission from biochar and SOC decomposition. Woolf & Lehmann (2012) reported that P FOM (FOM positive priming; Table 1) varies from 0% to 13.8% and the negative priming or stabilization coefficient varies from 0% to 3.6%. Estimates for P BIOM and P HUM (BIOM and HUM positive priming's  Tables 2 and 3). coefficients) were difficult to extract from the literature because of the challenge in separating priming effects from the contribution of the labile biochar on C miner-alization (Cross & Sohi, 2011). Zimmerman et al. (2011) studied C mineralization rates for 20 biochar types on five soils, and found that the magnitude of priming var- Fig. 3 Sensitivity analysis of biochar model parameter. The output variables shown in each panel are: dltbc: daily rate of biochar C decomposition; dlt_bc_co2: daily rate of biochar C lost to the environment; dlt_nbc: daily rate of nitrogen mineralization/immobilization due to biochar decomposition; soil CEC BC: soil cation exchange capacity after biochar application; soil pH BC: soil pH after biochar application; NH4 ads: soil NH4 adsorption due to biochar application; dlt_DUL: daily rate of change in drained upper limit; dlt_BD: daily rate of change in soil bulk density due to biochar application. The values in parentheses next to each parameter are the default values used in the sensitivity analysis. For symbols explanation and units see Table 1. ied from À52% to +89% at the end of the first year. Wang et al. (2015) in a review study reported lower biochar priming values on SOC (from À8% to +20%) and that the direction and the magnitude of the priming depended on various factors including biochar characteristics, soil properties and experimental duration (biochar age). The priming coefficients in our biochar model are linked to Soil BC (Eqns 4-6), which is determined by biochar characteristics (MRT 1 , MRT 2 , BC o , f carbon , and f labile ) and soil properties (w f , t f and n f , see soil hydrology and BD later). In contrast, our model does not include yet an age factor because of the lack of specific data to support modeling. Nevertheless, potential users can set different priming values in the model (Table 1) and test hypotheses related to the magnitude of the priming effect on soil processes and crop yields (see Results section).

Biochar effects on soil N mineralization and immobilization
During the decomposition of the biochar (Eqn 1), N is mineralized or immobilized depending primarily on the biochar C : N ratio, which is an input to the model (Table 1). Mineralization or immobilization of mineral N is determined as the balance between the release of N during biochar decomposition and the immobilization of mineral N during the formation of microbial biomass and new humic materials (Probert et al., 1998). The N needed for biochar decomposition is calculated as: where, dlt_ nbc_need is the daily amount of N required for biochar decomposition (kg N ha À1 day À1 ), dlt_bc_biom and dlt_bc_hum is the daily biochar C that goes to BIOM and HUM pools, respectively (kg C ha À1 day À1 ), CN_biom is the C/N ratio of the BIOM pool (default = 8; APSIM version 7.7), and CN_hum is the C/N ratio of the HUM pool (around 11-14; user defined). The N demand may be satisfied by the N that is released during the decomposition of the biochar, which is estimated as: where, dlt_nbc_released is the mineralized N during biochar decomposition (kg N ha À1 day À1 ) and CN_bc is the CN ratio of the biochar (Table 1; user defined). The C : N ratio of various biochar types varies from 7 to 769 in the literature (Atkinson et al., 2010, and it is typically measured by thermal combustion. We believe that the C : N ratio of the labile fraction of biochar may be different from the C : N ratio of the recalcitrant fraction, but due to the lack of specific data and difficulties in measuring individual biochar pool C: N ratios, we currently use the C : N ratio of the whole biochar. This is another gap in biochar research. Finally, the net rate of N mineralization or immobilization during biochar decomposition was calculated as: where, dlt_nbc is the N mineralized or immobilized during biochar decomposition (kg N ha À1 day À1 ). For example, in the case study shown in Fig. 4 approximately 2 kg N ha À1 year À1 was immobilized during biochar decomposition. The mineralized or immobilized N is added or subtracted to/from the NH 4 pool on a daily basis, and this can greatly affect subsequent N cycling processes such as nitrification and denitrification as well as plant N uptake (Fig. 1). Our assumption here Fig. 4 Simulated daily decomposition rate of biochar with their associated N immobilization (left panels) and their stress modifiers used to reduce the potential decomposition rate to actual (right panel; 1 means no stress; 0 means full stress). The biochar parameters values used are shown in Table 1 under bicohar type 1. The simulation refers to a maize crop growing in Ames, Iowa in a loamy soil (see Tables 2 and 3). The sharp change in the nitrogen stress factor on day of the year 125 is due to 180 kg N ha À1 fertilizer application.
is that a biochar with a CN ratio above 25 (parameter opt_bc; Table 1) will immobilize N, but the rate of N immobilization is generally very low because of the slow rate of biochar decomposition (Fig. 4). A sensitivity analysis of the parameters affecting Eqn (9) showed that all are important (Fig. 3c) and can affect the process in very different directions. This varying sensitivity highlights the complexity that exists in the soil-biochar system.

Biochar effects on soil CEC
The biochar model has functions to estimate soil CEC with (Soil CECBC ) and without (Soil CEC ) biochar present. When Soil CEC data are not available, the model estimates initial Soil CEC from SOC and noncarbonate clay% (n clay ) using a set of equations developed by Seybold et al. (2005) for different soil orders: The model then simulates the effect of biochar on soil CEC. Literature reports on the impact of biochar amendments on soil CEC are inconsistent with some reports showing an increase in CEC (Laird et al., 2010;Van Zwieten et al., 2010;Peng et al., 2011) and others showing no effect (Haefelea et al., 2011;Rogovska et al., 2014). Differences are most likely related to the biochar type, the amount added, the native soil CEC, and biochar age. Over time, in soil environments surfaces of biochar particles are oxidized forming various organic functional groups; dissolved organic compounds may be adsorbed on, or leached from, biochar surfaces; microorganisms may colonize biochar surfaces; carbonates and other ash components may be solubilized exposing new biochar surfaces; and the biochar equilibrates with the ambient soil pH. All of these aging processes influence the net impact of biochar amendments on soil CEC. We used a user defined effective biochar CEC (BC ECEC ; Table 1), which is assessed as the net change in CEC per kg of soil approximately 1 year after a biochar addition relative to the mass of biochar in 1 kg of the soil. The BC ECEC approach integrates all of the above aging effects without identifying the specific mechanisms. The default value in the model for BC ECEC is 187 cmolc kg À1 (Laird et al., 2010) while the range varies between 0.5 and 68 for fresh and 96 to 235 cmol kg À1 for aged biochars, respectively (Cheng et al., 2008;Gai et al., 2014). The impact of a biochar on soil CEC is calculated as: where Soil CEC and Soil CECBC are the CEC of the soil before and after a biochar addition, Mass fr is the mass fraction of biochar in the soil (g g À1 ) estimated from layer depth, soil BD and biochar in the soil [Mass fr = Soil BC /(f carbon ÁBD o Ád layer Á10 000)], and BC ECEC is the effective CEC of the biochar (cmolc kg À1 ). Sensitivity analysis showed that BC ECEC has a relatively small impact on Soil CECBC (Fig. 3d).

Biochar effects on soil pH
Most biochars are alkaline due to the presence of both carbonate rich inorganic ash, which is mixed with the biochar, and the presence of organic functional groups on biochar surfaces, both of which can react with protons (Yuan et al., 2011;Fidel, 2012;). Hence the impact of a biochar addition on soil pH is similar to the addition of agricultural lime. We used a user defined lime value (BC LV ; Table 1) to represent the amount of acid that can be neutralized by the biochar. Typically, biochar lime values range from 20 to 200 cmol c kg À1 . By contrast, pure CaCO 3 has a lime value of 2000 cmol c kg À1 . To estimate the impact of a biochar addition (or an agricultural lime addition) on soil pH we used a sigmoidal function (Nelson & Su, 2010) assuming that the buffering capacity of soil is proportional to the soil CEC: where Soil pH and Soil pHBC are the soil pH before and after biochar application, Mass fr is the mass fraction of biochar in the soil (g g À1 ), BC LV is the biochar lime value (cmolc kg À1 ), U pH and L pH are the upper and lower limit of soil pH (here set at 8.3 and 3.5 assuming the agronomic pH range is limited by the precipitation of carbonates and solubility of aluminum hydroxide, respectively), and P 1 is a fitted constant that adjusts the buffering capacity of soils (here 10). Among these parameters, the U pH and the initial soil pH have the greater impact on Soil pHBC (Fig. 3e). The soil pH changes rapidly after a biochar or lime addition or when the soil is amended with ammonium fertilizer otherwise the daily change in soil pH should be negligible. Equation (12) can result in different shapes in response to input parameters (Table 1)

Biochar effects on NH 4 adsorption and desorption
The adsorption of the NH 4 due to biochar addition is considered to be mainly due to the CEC. We quantified the effect of biochar on NH 4 adsorption using the Langmuir model, which is frequently used for describing sorption isotherms (Gai et al., 2014): where NH 4ads is the adsorbed NH 4 from the soil solution (mg L À1 ), NH 4 is the concentration of the ammonium in the soil solution (mg L À1 ), K ads refers to adsorption capacity (L mg À1 ) and is defined by the user in the biochar model (Table 1; Fig. 5), Soil CECBC and Soil CEC is the soil CEC before and after biochar application. Literature values for K ads vary from 0.0026 to 0.125 L mg À1 (Table 1; Ding et al., 2010;Gai et al., 2014). Equation (13) is very sensitive to the K ads parameter (Fig. 3f). Equation (13) shows that the ability of the soil to adsorb NH 4 due to the biochar addition is proportional to CEC, which is a function of the added biochar in the model (see Eqn 11). A CEC of 1 cmol kg À1 is equivalent to a maximum sorption of 0.14 mg NH 4 g À1 biochar, which means that biochar's maximum adsorption capacity is from 0.07 to 9.5 mg NH 4 g À1 biochar assuming fresh biochar CEC values from 0.5 to 68 (Gai et al., 2014;Jassal et al., 2015). However, recent literature findings have shown even higher maximum adsorption capacities, i.e. 1.1-28.4 (Jassal et al., 2015) and 2.6-15.4 (Gai et al., 2014), which led to a hypothesis that most of the sorption must be due to physical entrapment of NH 4 in (fresh) biochar pores (Spokas et al., 2011;Jassal et al., 2015). The above values clearly show that the ability of the biochar to adsorb NH 4 is not universal, but depends on several factors including biochar type, age, and soil type (Ding et al., 2010;Nelissen et al., 2012;Yao et al., 2012;Gai et al., 2014;Sarkhot et al., 2014;Jassal et al., 2015). Equation (13) can be extended further to account for more mechanisms or variables such as O%, C%, O/C ratio, and H/C ratio, which have been found to correlate well with the rate of NH 4 adsorption (Gai et al., 2014).
In contrast to NH 4 adsorption, much less information is available about the release of the adsorbed NH 4 to the soil solution. If the soil-biochar system adsorbs NH 4 , then theoretically, at some point in time it should release it to the soil solution until equilibrium is achieved. Based on this assumption, we modeled NH 4 desorption similarly to adsorption (see Eqn 13). The rate of desorption is driven by a user defined K des parameter (Table 1; Fig. 5). Figure 5 illustrates K ads and K des effects on simulated NH 4 over time. The value of K des is largely unknown and we used in the biochar model similar values as the K ads , despite the fact that some short duration experiments have shown lower values. Clearly more research is needed in this area given the importance of NH 4 in N cycling (Sarkhot et al., 2014;Jassal et al., 2015).

Biochar effects on soil water and BD
The effect of biochar on APSIM soil water retention parameters [lower limit (LL, mm mm À1 ), drained upper limit (DUL, mm mm À1 ), and saturation point (SAT, mm mm À1 )], as well as soil BD (in g cm À3 ) was modelled by incorporating relationships between SOC and these parameters from the Saxton & Rawls (2006) pedotransfer functions into the biochar model (Fig. 1). Inputs for the model are clay and sand content of the soil and the rate of change in SOC due to biochar addition. The rate of change in SOC can be positive as a result of biochar addition or negative when SOC declines over time. Clay and sand contents mediate the rate of change. Outputs from the model are the rate of change in LL, DUL and BD on a daily basis. The rate of change in BD is used to compute SAT from the change in total porosity which is a function of BD [SAT = (1 À BD/2.65) 9 0.95]. The value of 2.65 is the particle density of the mineral soils (g cm À3 ) and 0.95 reflects the air filled porosity (Dalgliesh & Foale, 2005).
The Saxton & Rawls (2006) equations are driven by the amount of SOC and have not been tested for Fig. 5 Simulated soil NH 4 -N (mg L À1 ; 0-10 cm depth) response to adsorption (K ads ) and desorption (K des ) coefficients. Biochar type 1 parameter values were used in this simulation (see Table 1) expect for the K ads and K des parameters, which values are shown in the panel. The simulation refers to a maize crop growing in Ames, Iowa in a loamy soil (see Tables 2 and  3). The sharp increase in NH 4 -N on day of the year 123 is due to application of 100 kg N ha À1 urea-N fertilizer in this example. biochar. We hypothesized that the effect of biochar on soil hydrology will be somewhat different from that of SOC as biochar C is qualitatively different from SOC. Literature supports this hypothesis as it shows that different types of biochars have very different surface pore spaces (0.2 to 8.2 nm; Gai et al., 2014;Lim et al., 2016) and variable effects on DUL, LL, and BD, which cannot be explained by the amount of biochar applied alone (Brockhoff et al., 2010;Laird et al., 2010;Liu et al., 2012;Basso et al., 2013;Cornelissen et al., 2013;Herath et al., 2013). These studies showed that the LL is rarely changed due to biochar addition, DUL and SAT can be changed but the rate of change depends on the initial SOC before application and on biochar type. To account for the literature findings and to add further flexibility in the model, we developed and incorporated into the biochar model 0-1 quality modifiers on the rate of change in LL, DUL and BD (Q LL , Q DUL , and Q BD , respectively). The Q LL was set to zero (Table 1) because biochar does not change LL, in contrast to Saxton & Rawls (2006) equations that show a change in LL. The Q DUL was set as a function of SOC [Q DUL = exp (ÀK DUL ÁSOC)]. The K DUL parameter is used to correct for the negligible effect of biochar on DUL in soils rich with organic C. The Q BD was set as function of SOC similar to Q DUL but with a different input parameter (see K BD in Table 1) because the effect of biochar on SAT is usually larger than that on DUL (e.g. Liu et al., 2012). The parameters Q LL , K DUL and K BD are user defined (Table 1). By setting Q LL to 1 and K DUL and K BD to zero, the quality modifiers are cancelled, and the equations take the original Saxton & Rawls (2006) form. A sensitivity analysis of the biochar hydrological parameters is presented in Fig. 3g.
The biochar model is capable of representing the effect of soil type (via texture), the physical aspects of biochar (amount added) as well as quality aspects of biochar (pore space) on soil water retention and BD. As SOC increases due to biochar addition, DUL and consequently PAWC increase and the magnitude is larger in coarse textured than in fine textured soils. Figure 6 shows simulated PAWC response to biochar additions in soils with different levels of organic C. The magnitude of the simulated changes in PAWC is consistent with literature findings (Brockhoff et al., 2010;Laird et al., 2010;Liu et al., 2012;Basso et al., 2013;Cornelissen et al., 2013;Herath et al., 2013).
For the biochar effect on BD we considered a tillage factor in addition to the effect of C. This is because most laboratory incubation studies have shown a decreased BD in biochar amended soils relative to controls (Laird et al., 2010;Lim et al., 2016), but results from field studies are more equivocal with most studies either showing no significant effect or a decrease in BD after biochar application (Haefelea et al., 2011;Rogovska et al., 2014). Under field conditions biochar is commonly incorporated by tillage and it can be difficult to distinguish the biochar effect from the tillage effect on soil BD. To address this challenge we used the law of maximum absolute change between the two processes: biochar effect on BD (through the Saxton and Rawls pedotransfer functions) and tillage effect on BD (through Eqn 14). To simulate tillage effects on BD we used the following equation (Andales et al., 2000): where BD till is the calculated BD (g cm À3 ) that changes over time. BD o is the user input initial BD (g cm À3 ). BD a is the BD just after tillage operation (this value is defined by the user and varies with tillage equipment; Table 1). The AS parameter reflects soil aggregate stability and equals 5Á(1-0.205 Á SOC), and C rain is the cumulative amount of precipitation since the latest tillage event (see Andales et al., 2000). Figure 3h illustrates the sensitivity of BD related parameters on Eqn (14). In general, the biochar model decreases BD after tillage to BD a ( Fig. 7; Table 1). Then, with precipitation as the main driver, BD returns to the initial condition or to the new condition set by the added biochar C (Fig. 7). Literature findings have shown an average reduction in soil BD of 11.1 AE 8.6% due to biochar addition (Lim et al., 2016). The reduction in BD due to biochar in this example was 5.8% ( Fig. 7; loamy soil). The final BD estimate from the biochar model is used to compute SAT on a daily basis. It should be mentioned that the dynamic changes in LL, DUL, and SAT have a substantial impact on the simulation process because they are linked to water modifiers that can alter C and N cycling and plant growth Fig. 6 Simulated changes in plant available water content (PAWC = DUL À LL) by the biochar model as affected by the application rate and soil organic carbon. A biochar type 1 parameter values were used in this simulation (see values in Table 1; expect the BCo which varied) and a loamy soil (see Table 2, expect the initial soil organic carbon which varied). (Figs 1 and 4). As the biochar model affects both DUL and SAT, the SW will change as well. This in turn affects the water filled pore space (WFPS = SW/SAT), which is one of the main drivers for N 2 O emissions in the model. This shows the feedback mechanisms that exist in the model (Fig. 1). It would be possible to incorporate biochar effects on saturated hydraulic conductivity (K sat ) in the model through the Saxton and Rawls equations, but it was not finally included in this version for two reasons: (i) literature findings are still contradictory (see Table 3 in Lim et al., 2016) and (ii) a preliminary sensitivity analysis of K sat (i.e. 100 mm day À1 AE 50% change) had a very minor impact on the simulation output (data not shown).

Summary of biochar parameters
We classified the biochar model parameters in three major categories (Table 1): (i) those that are highly dependent on biochar type (production method 9 feedstock source 9 production temperature) and management (application amount, incorporation depth); (ii) those that are determined by soil and biochar interactions (effects on soil hydrology, priming); and (iii) finally those that are assumed to be constants in the model. All of the above parameters are important and affect various soil processes and properties differently (Figs 3 and 8). Parameters for the first category can be extracted from various literature resources or measured as outlined above or can be taken from the Biochar Engineering web tool (http:// spark.rstudio.com/veromora/BiocharEng/; Morales et al., 2015) or other biochar databases (http:// biochar.ucdavis.edu/). For the second and third category, literature information is incomplete and specific experiments are needed to fill modeling knowledge gaps (see Discussion later).

Preliminary model testing
We used experimental data from Rogovska et al. (2014) to test various aspects of the biochar model. Briefly, in this experiment a woodchip biochar (see parameter values in Table 1) was applied at six rates (range from 0 to 98 Mg ha À1 ) to a maize field located in central Iowa. For more details we refer to Rogovska et al. (2014). Our simulation was for a continuous maize crop (Table 2) growing in a Clarion soil (Table 3) in Ames, IA.

Long term sensitivity analysis on systems performance
A sensitivity analysis of the individual biochar model parameters on key soil processes was presented above (Fig. 3). However, this type of sensitivity analysis (Fig. 3, change one parameter at a time while keeping all others constants) has some important limitations when the results are to be scaled up to the final product (i.e. grain yield; Fig. 1). This is because the results of the sensitivity analysis are dependent on the soil type used, cropping system, climate, management, and biochar parameters. The fact that there is not yet a biochar type that can serve as a reference in the sensitivity analysis, makes the initial choice of the biochar parameters very critical.
Given all these issues we performed a second sensitivity analysis in which we explored the effect of different biochars (named BC1, BC2, BC3; Table 1) on diverse cropping systems (crop, management, climate and soils; Tables 2 and 3). In this analysis we output several production and environmental variables in the longterm (30 years) to get a complete picture of systems performance. The results are presented in Fig. 8. On average, across five cropping systems and three biochar types, the order of sensitivity of various soil-plant  Table 1 (biochat type 1). The simulations refer to a maize crop growing in Ames, Iowa in a loamy soil (see Tables 2 and 3). processes and properties to biochar application was (from high to low): SOC > N 2 O emissions > PAWC > soil BD > soil pH > soil net N mineral-ization > CO 2 emissions > WFPS > grain yield > crop biomass. The direction and the magnitude of the change (compared to the control) was highly variable among  Tables 2 and 3) as affected by three different biochars (see Table 1 for the parameter values). The duration of each simulation was 30 years and the results averaged across this period (expect the value of SOC) and compared to the control simulation (without biochar) for each cropping system. A negative value indicates that the biochar reduced the output variable (yield, CO 2 emissions, etc.) and vice versa. The values in each panel and cropping system indicate the long term value of the control simulation. cropping systems and biochar types: SOC varied from 11% to 86%, N 2 O emissions from À10% to 43%, PAWC from 0.6% to 12.9%, soil BD from À6.5% to À1.7%, soil pH from À0.8% to 6.3%, soil net N mineralization from À19% to 10%, CO 2 emissions from À2.0% to 4.3%, WFPS from À3.7% to 3.4%, grain yield from À3.3% to 1.8%, and crop biomass from À1.6% to 1.4%. The effect of biochar on various variables was largest in the loam sandy soil (Australian wheat) and smallest in clay loam soil (US maize) ( Table 2, Fig. 8).
Among the tested variables, SOC, PAWC and BD showed a monotonic increase or decrease while all the other variables showed both increased and decreased trends over time (Figs 8 and 9). Our results indicated that the effect of biochar on crop production variables was much smaller compared to the effect on environmental variables such N 2 O and CO 2 emissions. The reason might be the feedback mechanisms that exist in biochar-soil-crop-atmosphere systems and are integrated into our model. These mechanisms may substantially counterbalance biochar effects on individual processes or the biochar effects are not large enough to affect crop yields. For example, an increase in SW that favors crop growth might trigger NO 3 -N loss via denitrification or leaching, which can limit crop growth. Given that crop yield is the final product of several interrelated soil-plant processes (Fig. 1), our model analysis suggests that the effect of biochar diminishes all the way to the final product, grain yield. This result also shows that it is possible for biochar to decrease GHG emissions without negatively impacting crop yields which is important (see the magnitude of the biochar response: À10% to 43% for N 2 O vs. À3.3% to 1.8% for grain yield). This can be achieved by identifying appropriate biochar types within a cropping system and environment. The biochar model coupled with appropriate optimization functions can assist in this task.
The magnitude of the obtained results in Fig. 8 is consistent with literature findings on biochar and reveal for the first time the trade-offs that exist in the biocharsoil-crop-atmosphere system. For example, the application of BC2 resulted in higher grain yield and biomass production and at the same time in higher GHG emissions. That was caused by the higher PAWC and soil N availability (see PAWC and N mineralization in Fig. 8). The decrease in WFPS was not enough to reduce N 2 O emissions. In contrast, the BC3 produced the opposite responses compared to BC2. It decreased both grain yield and GHG emissions variables. Literature studies have shown decreased and increased N 2 O and CO 2 emissions due to biochar application (Clough et al., 2010;Bruun et al., 2011a;Case et al., 2012Case et al., , 2014Case et al., , 2015Ameloot et al., 2013;Cayuela et al., 2014Cayuela et al., , 2015Liu et al., 2015). Our study adds significant value to the biochar research because it reveals for the first time the reasons responsible for increased or decreased GHG emissions and grain yields (see parameters in Table 1, mechanisms in Fig. 1 and the magnitude of the response in Fig. 8). Figure 8 presents the average effect of the biochars across 30 years of cropping. It should be mentioned that the response of biochar on production and environmental variables was not consistent across years, but variable with positive, neutral and negative responses over time (Fig. 9). This is consistent with observations and highlights the importance of climate and the carry-over effects on the biochar response. This also verifies the need for long term field assessments of biochar before it can be recommended for large-scale applications (Jones et al., 2012;Nelissen et al., 2012). Modeling can assist in this task by extending knowledge, but long-term field trials are ultimately needed to calibrate and validate modeling.
Preliminary biochar model evaluation Figure 10 shows the goodness of fit of the biochar model to several soil-plant variables obtained under field conditions (Rogovska et al., 2014). The calibrated biochar parameter values used in this simulation are shown in Table 1. Overall, the model matched experimental observations with a mean relative absolute error from À0.4% to 13.1%, which is very good. The prediction of soil BD and pH was the most accurate, while the ability of the model to predict more complex variables such as SOC was lower (Fig. 10). In this particular set of data (year: 2012, location: Ames, IA), the model predicted that corn grain yields were limited by SW, which agrees with the experimental observations (Rogovska et al., 2014). The model predicted that biochar application positively influenced SW (Fig. 10f), but this increase was not enough to compensate for the extremely dry conditions that occurred in 2012.

Discussion
The potential of the biochar model In a recent state-of-the-art biochar review study, Jeffery et al. (2015) argued for the need of a mechanistic understanding of the effects of biochar application in order to quantify and predict benefits and trade-offs. In this paper we presented the first complete biochar simulation model that integrates and synthesizes broad biochar knowledge and enables simulation of biochar effects on diverse cropping systems growing in different environments. We also demonstrated the potential of the model to assist system level understanding and assessments of biochar and their trade-offs (Figs 8 and  10). The biochar model can assist biochar research in various ways. It can be applied to analyze data from the vast lab incubation studies or field experiments and to extrapolate the results beyond the study period in time and space in order to evaluate long-term consequences, which are definitely needed. Once the model fits the experimental observations (calibration task), it can provide information on soil-plant processes that have not been measured experimentally and therefore allow for a system analysis of the biochar effects. The model can be used to test or generate hypotheses to inform experimental research or to perform experiments that are expensive, time consuming and difficult to implement in real-life situations. For example, what is the optimum biochar rate and type to maximize yields and minimize GHG emissions from a specific cropping system in a particular location? The model can provide an answer to this question relatively quickly and inexpensively. By running a factorial simulation with multiple biochar rates, multiple biochar types within a rate, different management practices (N application rate, cultivars, planting dates, tillage type and time) for each biochar type, and considering 30 years of weather variability, the system can be optimized and the win-win combination can be found theoretically. Then the best scenario can be tested experimentally to verify model predictions. The biochar model can be coupled with other advanced simulation platforms such as pSIMS (Elliot et al., 2014) to assist regional-to global-scale assessments of biochar.

Current status and future steps
As this was the first attempt to develop a generic tool for system assessments, we initially focused on developing algorithms and predicting the direction and the magnitude of the biochar responses (Figs 4-8). Our work on model conceptualization and development is driven by diverse literature information and not from a site-crop-biochar specific experimental dataset. Usually modelers build their models using a particular experimental dataset and use another data set for validation (Lychuk et al., 2015). In contrast, we developed a comprehensive theoretical framework by integrating a wide range of biochar research findings and synthesized it from a systems perspective and performed a Fig. 9 Simulated temporal variability of biochar response (BC1 ; Table 1) to grain yield, annual cumulative CO 2 emission and annual cumulative N 2 O emissions as compared to the control (no biochar). The simulation refers to US maize growing on a loamy soil (see Tables 2 and 3). A negative value indicates that the biochar reduced the output variable (yield, CO 2 emissions, etc.) and vice versa.
preliminary testing of the model. Comprehensive model calibration and validation is the next step. Given the diversity in biochar types and their complex interactions with crops, soils and climate, the task of calibrating and validating the biochar model must be iterative over time and include multi-locations, multi-species and multibiochar type trials. In that respect, the APSIM modeling platform becomes even more important as it is frequently updated, freely available in the public domain and contains several crop models and soil/environmental models along with flexible managers to assist investigation of biochar effects on diverse cropping systems in different environments. The current version of the biochar model captures the most important biochar effects on soils and extrapolates the effects to crop yields by elaborating mechanistic crop models available in APSIM. Regarding the N cycling which is of particular importance, the biochar model affects N cycling via the following direct and indirect pathways (Fig. 1): 1. increase in the size of the SOC pools that can stimulate higher N mineralization rates (Eqn 2) 2. changes in the turnover rates of FOM and SOC pools due to priming (Eqns 4-6) 3. changes in CEC and pH that affect urea hydrolysis and nitrification rates (Eqns 11 and 12); 4. changes in NH 4 pool due to mineralization/immobilization of N caused by the decomposition of biochar (Eqns 2, 7-9); 5. changes in NH 4 pool due to adsorption and desorption caused by biochar CEC (Eqns 11, 13); 6. changes in field capacity and SW which affect water stress factors and therefore indirectly affect all of the N transformation processes (Fig. 4) and plant growth (Fig. 8); 7. changes in BD, which affect SAT and therefore WFPS. WFPS, in part, drives N 2 O emissions together with NO 3 , temperature and C in APSIM (Thorburn et al., 2010). Our next steps in biochar modeling are: (i) finalize on-going calibration tasks (e.g. Fig. 10) as well as expand model calibration tasks to diverse environments, cropping systems and biochar types; (ii) develop a database with calibrated biochar parameters (Table 1) that will allow subsequent system analyses of biochar effects by the users; (iii) add new mechanisms and improve existing ones in the model. For instance, quantification and modeling of biochar effects on root  Table 1) and the measurements were obtained in 2012. The simulation refers to a corn crop (Table 2) growing in a Clarion soil (Table 3) in Ames, IA. functioning, seedling emergence, possible allelopathy effects, and distinction between fresh and aged biochar effects are some of the future challenges, among others.

Knowledge gaps and research priorities to support mechanistic modeling
In this study we encountered several knowledge gaps that would benefit from further research. We propose the following for prioritization: (i) NH 4 absorption and release patterns in the long term. The size of the NH 4 pool over time is the key to achieve production and environmental benefits (Fig. 1); (ii) experimental verification of the Saxton & Rawls (2006) pedotransfer functions for biochar response because correct prediction of soil hydrology affects the overall system substantially ( Fig. 1); (iii) estimation of the C : N ratio of the labile and recalcitrant biochar separately because its value strongly affects the rate of N immobilization (Eqn 2); (iv) development of biochar specific temperature and SW modifiers; (v) experimental verification of the biochar C partitioning coefficients used in the model (Eqn 3; Table 1); (vi) quantification of the duration of priming, i.e. does the priming effect last for days, weeks, months, years and how is this related to biochar type, soil, and management?; (vii) distinction of fresh from aged biochars; (viii) the effects of fertilizers and residue management on pH (Yuan et al., 2011) needs to be addressed for a better representation of soil processes in the model.
In summary, this study set the basis for a mechanistic and system assessment of biochar effects on soils and crops and their trade-offs. We envision that the biochar model will shed light on critical questions and generate hypotheses for experimental research. The model has a great potential to assist agronomists, soil scientists, economists, and policy makers to design sustainable and profitable biochar based cropping systems.