Jukuri, open repository of the Natural Resources Institute Finland (Luke) All material supplied via Jukuri is protected by copyright and other intellectual property rights. Duplication or sale, in electronic or print form, of any part of the repository collections is prohibited. Making electronic or print copies of the material is permitted only for your own personal use or for educational purposes. For other purposes, this article may be used in accordance with the publisher’s terms. There may be differences between this version and the publisher’s version. You are advised to cite the publisher’s version. This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Matti Sihvonen, Jussi Lintunen, Elena Valkama & Kari Hyytiäinen Title: Management of legacy nutrient stores through nitrogen and phosphorus fertilization, catch crops, and gypsum treatment Year: 2020 Version: Published version Copyright: The Author(s) 2020 Rights: CC BY 4.0 Rights url: http://creativecommons.org/licenses/by/4.0/ Please cite the original version: Sihvonen M, Lintunen J, Valkama E, Hyytiäinen K.Management of legacy nutrient stores through nitrogen and phosphorus fertilization,catch crops, and gypsum treatment. Natural Resource Modeling. 2020;33:e12289. https://doi.org/10.1111/nrm.12289 Natural Resource Modeling. 2020;33:e12289. wileyonlinelibrary.com/journal/nrm | 1 of 40 https://doi.org/10.1111/nrm.12289 Received: 15 January 2020 | Accepted: 29 September 2020 DOI: 10.1111/nrm.12289 Management of legacynutrient stores through nitrogen and phosphorus fertilization, catch crops, and gypsum treatment Matti Sihvonen1 | Jussi Lintunen2 | Elena Valkama2 | Kari Hyytiäinen1 1University of Helsinki, Helsinki, Finland 2The Natural Resources Institute Finland (Luke), Helsinki, Finland Corresspondence Matti Sihvonen, University of Helsinki, Helsinki, Finland, Email: matti.sihvonen@helsinki.fi Abstract We develop a modeling framework, based on discrete‐time dynamic optimization, to study the effect of legacy nutrient stores and soil nutrient dynamics on optimal nutrient management and agri‐environmental policy in crop production. Three alternative measures are studied to reduce nutrient loss: reduced fertilization, nonlegume catch crop cultivation and gypsum amendment. According to our results, landowner can bring down excessively high phosphorus stocks causing environmental damage within decades, by si- multaneous optimization of the nitrogen and phosphorus fertilizers on the economic basis of profit maximization. Our results suggest that ni- trogen loss abatement with catch crops is socially optimal, whereas the use of gypsum is well justi- fied as a temporary measure on soils with high soil phosphorus levels. A dynamic tax‐subsidy‐scheme, which takes into account the current soil nutrient levels and field attributes such as soil texture, can - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. Natural Resource Modeling published by Wiley Periodicals LLC enforce the socially optimal outcome. The welfare losses of the static steady‐state‐based tax‐subsidy‐ schemes are increasing functions of the legacy nutrient stores and soil's ability to hold nutrients. Recommendations for Resource Managers • We develop a modeling framework to study the effect of the legacy nutrient stores and the soil nutrient dynamics on the optimal nutrient man- agement and agri‐environmental policy in crop production. • Nonlegume catch crop cultivation is a socially op- timal long‐term measure for nitrogen loss abate- ment, whereas phosphorus loss abatement with gypsum is socially optimal on soils with high soil phosphorus levels. • A dynamic tax‐subsidy‐scheme, which is adjusted annually according to the soil nutrient stocks, leads to social optimum. Although this can be difficult to implement in practice, it can be useful in the derivation of the simpler, static tax‐subsidy‐ schemes. • If a gypsum subsidy is paid for those years, where the soil P level is above the threshold level for the gypsum application, the welfare loss of the static steady‐state‐based tax‐subsidy‐schemes is al- most zero. • Simultaneous adjustment of the N and P fertilizer rates enables the use of simple, static and soil‐ texture‐ignorant tax‐subsidy schemes, without a notable social welfare loss KEYWORDS dynamic optimization, fertilizer management, gypsum, legacy nutrient store, nonlegume catch crop, nutrient pollution, tax‐subsidy‐scheme 2 of 40 | Natural Resource Modeling SIHVONEN ET AL. 1 | INTRODUCTION Efficient management of nitrogen (N) and phosphorus (P) is essential in crop production because those are the two most important nutrients for crop growth. If leaked to aquatic ecosystems, however, N and P nutrients are a major cause of eutrophication (Carpenter et al., 1998; Vitousek et al., 1997). Eutrophication is the primary global water quality issue for most of the freshwater and coastal marine ecosystems (Smith & Schindler, 2009). Agri‐environmental programs have been successful at reducing nutrient losses at a field, but significant water‐quality improvements in the receiving water ecosystems remain unobserved. The major reason is the continuing release of nutrients from legacy nutrient stores, which have accumulated in agricultural soils over decades of fertilizer and manure applications (Sharpley et al., 2013; Stackpoole et al., 2019; Van Meter et al., 2018; Whiters et al., 2014). Legacy stores can be sufficient for maintaining eutrophication in the receiving water ecosystems for decades after fertilization is reduced or stopped (Carpenter, 2005; Watkins, 1998). However, soils with high nutrient legacies are fertile and produce high yields (Barrow, 1980; Dodd & Mallarino, 2005). Therefore, a soil nutrient stock is an indicator of both the productive capacity of fields and the nutrient load on water bodies. In Finland, farmers used high P fertilizer rates in 1940–1980 to increase soil P stocks of the fields because the native P reserves of Finnish soils are poorly available for the plants (Saarela, 2002; Tilman et al., 2002). Soil P stocks have been rising until the present day, and there are still many fields with high soil P stocks (Figure 1). About half of agricultural soils have such high soil P levels that a yield response to P fertilization is unlikely (Ylivainio et al., 2014). Moreover, about 60% of the anthropogenic P and 50% of the N loading in Finland originates from agricultural land (Vuorenmaa et al., 2002). The total N and P loads from Finland to Baltic Sea were 82,500 and 3,255 tons in 2014 (HELCOM, 2018). If the marginal cost of N load is 6.6 €/kg (Gren & Folmer, 2003), the value of the load would be 545 million €. Similarly, the value of the P load from Finland to Baltic Sea would be 155 million €. Despite various mitigation policies, riverine nutrient export from agricultural lands remains high (Räike et al., 2019). Past research has shown that a long‐term approach is necessary to capture the effect of legacy nutrient stores on a socially optimal nutrient management (Whiters et al., 2014). Also, accounting for the externalities of crop cultivation requires a combination of policy instruments and abatement measures, because no single solution is likely to suffice (Pretty et al., 2001). To this end, in this study we focus on a socially optimal joint long‐term management of P and N nutrients together with nonlegume catch crop cultivation and gypsum abatement in spring Coarse soils 0 20 40 60 80 Soil P level (mg/l) 0 2000 4000 6000 8000 10000 Fr eq ue nc y Clay soils 0 20 40 60 80 Soil P level (mg/l) 0 1000 2000 3000 (a) (b) FIGURE 1 Distribution of the soil P stocks on Finnish agricultural soils (based on the initial data describing soil P levels in field plots used in the VEMALA model (Huttunen et al., 2016) SIHVONEN ET AL. Natural Resource Modeling | 3 of 40 barley production on heterogeneous agricultural land. Heterogeneity arises from the different legacy nutrient stores and soil textures of the field hectare. The economic theory of nonpoint‐source water pollution control dates back to Griffin and Bromley (1982), Shortle and Dunn (1986), and Segerson (1988). Xepapadeas (1992) proposed a general dynamic framework for the determination of the intertemporal incentive schemes for diffusion pollution. Schnitkey and Miranda (1993) conducted a steady‐state analysis for controlling P runoff from farms combining livestock and crop production. Larson et al. (1996) developed a welfare‐theoretic approach for comparing second‐best single‐tax policies using a crop yield and pollution production model. Goetz (1997) analyzed optimal private and social soil use paths when farmer's annual decisions consist of crop mix and use of inputs, while considering the nonlinear dynamic soil losses. Yadav (1997) studied optimal N and P fertilizer use in a dynamic setting with groundwater contamination externalities. Berntsen et al. (2003) studied effects of the three different nitrogen taxation scenarios for four farm types. Martıńez and Albiac (2004) studied optimal input taxes on N fertilizer and irrigation water. Iho (2010) and Iho and Laukkanen (2012a) studied an input tax on P fertilizer or soil P stock in a situation where planner's choices included P fertilizer rate and the width of a vegetative filter strip, whereas Iho and Laukkanen (2012b) considered a similar situation where the abatement measure was gypsum amendment. Gren et al. (2013) and Ahlvik et al. (2014) studied cost‐effective spatial and dynamic allocations of N and P loss abatement for reaching environmental targets in Baltic Sea. Despite the extensive literature of a non‐point source pollution control, economic literature on a long‐term management of both N and P fertilizers and optimal environmental policy accounting for the nutrient carryover information is limited. In particular, there appears to be no study where both soil P and N are state variables in an economic analysis. We address the following questions: (a) how socially and privately optimal fertilization patterns differ from each other on soils with different legacy nutrient stores? (b) What kind of tax‐subsidy‐schemes can encourage producers to engage in a socially optimal nutrient management and abatement? (c) How the soil nutrient dynamics affect optimal fertilization patterns, abatement, and agri‐environmental tax‐subsidy‐schemes? (d) What are the welfare gains of the increasingly detailed tax‐subsidy‐schemes? To answer these questions we apply discrete‐time dynamic optimization. Quantitative assessment is based on the models and data de- scribing barley (Hordeum vulgare L.) cultivation in southern Finland on coarse and clay textured soils. Our modelling approach extends the approach by Jomini et al. (1991), Bäckman and Lansink (2005), and Lambert et al. (2007) who studied economically optimal P and N rates (in slightly different settings) while considering P carryover effects, by including N carryover, nutrient loss processes, environmental policy, and several nonlinear interactions among the model components in the analysis. We also extend the approach by Iho (2010), Iho and Laukkanen (2012a, 2012b) by including N in the analysis. More specifically, we focus on a situation where a decision maker has access to three alternative measures to control leakage of N and P: application of nonlegume catch crops, gypsum amendment, and reduced fertilization inputs. 2 | THE MODEL 2.1 | The setup We consider crop production on an average representative hectare. A time step of the model is 1 year, denoted by t {0,1,2, …}∈ . Each year, nitrogen and phosphorus fertilization inputs (kg/ha), Nt and Pt, respectively, are used to supplement the natural supplies of soil nutrients to 4 of 40 | Natural Resource Modeling SIHVONEN ET AL. increase a barley yield (kg/ha), yt . Phosphorus accumulates in the soil over time if P fertilizer inputs exceed P removal. Each period of the planning horizon is characterized by a certain level of a plant‐available soil P stock (mg/L), denoted by xt. Plant‐available soil P stock comprises inorganic P dissolved in water/soil solution that is readily available for plant uptake.1 P car- ryover equation describes the evolution of xt from 1 year to another: x x x P y= + ϑ( , , ),t t t t t+1 (1) where transition function, ϑ, defines the rate of change of xt. The rate of P accumulation is an increasing function of the P fertilization, and a decreasing function of the annual crop yields because of the P crop‐uptake effect. Each period of the planning horizon is also characterized by a certain level of potentially plant available soil N (kg/ha), that is, N in mineralized form, nt (cf., Sela et al., 2017; Thomas, 2003; Watkins, 1998). We define the N carryover equation as follows: ( )n n ϕ n N y e= + , , , ,t t t t t tN+1 (2) where a transition function, ϕ, defines the rate of change of nt. This function is an increasing function of N fertilizer rate, a decreasing function of the annual yield because of the N crop‐ uptake effect, and a decreasing function of the annual N leaching to waterbodies, etN . Annual N and P fertilization and soil nutrient stocks have negative environmental con- sequences in the form of N and P losses (kg·ha−1·yr−1) to water ecosystems, etN and etP, re- spectively. Landowner can reduce the nutrient losses by reducing fertilization or by using abatement measures. The cultivation of nonlegume (non‐N‐fixing) catch crops reduces effec- tively N losses (Aronsson et al., 2016).2 Each year a landowner decides the share of a hectare that is devoted to catch crop cultivation, φ [0,1]tcc ∈ ; we assume that the choice of catch crop is annual. The current knowledge regarding the effect of catch crops on P losses is conflicting (Liu et al., 2014) and is omitted here. Gypsum treatment, denoted by φ [0,1]tgyp ∈ , is an alternative measure for reducing P loss that has gained widespread interest (Kleinman et al., 2019). Fur- thermore, φtgyp has no effect on soil P (Ekholm et al., 2012). We assume also that there is no direct interaction between φtgyp and φtcc. A one‐time gypsum application remains effective 5 years (Ollikainen et al., 2018). In the analytical model, however, we assume that a planner makes the decision annually (to ease the notation). Annual P load, etP is an increasing function of both xt and Pt, because etP consists of two forms: the leaching of dissolved reactive phos- phorus (DRP) and the loss of particulate phosphorus (PP). DRP is mostly determined by xt (McDowell & Sharpley, 2001; Peltovuori, 1999; Schnitkey & Miranda, 1993; Sharara et al., 2017). Annual N load, etN is also an increasing function of both nt and Nt (Van Meter et al., 1998; Van Meter et al., 2018). We describe N and P loss processes by the following functions: ( )e e N φ n= , , ,tN N t t tcc (3) ( )e e P φ x= , , .tP P t t tgyp (4) 1There are several P pools in soil: (a) plant‐available P in soil solution, (b) sorbed P, and (c) mineral P. In here, soil P stock refers to soil test P (STP), which provides an indication of the level of soil P that is potentially available to a plant. Typically, STP is extracted from the soil with acid ammonium acetate, pH 4.65 (Vuorinen & Mäkitie, 1955). 2Catch crops retain nutrients in the root zone, whereas cover crops prevent soil erosion and minimize the risk of surface runoff through the improved infiltration. SIHVONEN ET AL. Natural Resource Modeling | 5 of 40 According to the optimum law of Liebscher, annual P fertilizer response is a decreasing function of soil P. Thus, there is a trade‐off between the immediate benefit of P fertilization and the benefit of soil P accumulation. One can consider plant available soil N and added N fertilizer as substitutes: total available N to crop in a given period is the sum of applied N and soil N (Park et al., 2007; Stauber et al., 1975). Soil N and N fertilizer are, however, imperfect substitutes because their effect on a yield is different (Stoecker & Onken, 1989). We also assume that there is no cross‐effect between N and P responses, because the functional form in the numerical model excludes the interaction of the fertilizer inputs. We assume that catch crop does not improve yield through enhancing soil fertility and crop N supply.3 Instead, catch crop reduces yield due to the competition between a main crop and a catch crop for various factors of growth (Valkama et al., 2015). Consequently, catch crops have a positive indirect effect on soil P via the reduction in crop‐uptake (Lemola et al., 2000). We assume that gypsum application has no effect on barley yield (Iho & Laukkanen, 2012b). Hence, the annual yield is ( )y f P N φ x n= , , , , ,t t t tcc t t (5) where f is a concave yield response function satisfying the law of diminishing marginal returns. The model specification resembles a Cobb‐Douglas Mitschelich‐Baule type production function, which is flexible regarding the extent of isoquant convexity, and it imposes plateau growth (Yadav, 1997). 2.2 | Problem of the social planner Social planner maximizes the net present value (NPV) of the revenues from crop production minus the monetary value of the environmental damage caused by nutrient pollution by choosing optimal levels of P and N fertilization and catch crop and gypsum application areas, Pt, N ,t φ ,tcc and φtgyp, respectively. Note that φtcc and φtgyp are continuous variables, but due to model structure they turn up as corner solutions: either they are applied on the whole field hectare or not. The optimization problem is { }{ } β p y p P p N p φ p φ μ e μ e x n x P N φ φ max − + + + − − Subject to equations (1) and (2), and are given, , , 0, {0, 1}, and {0, 1} P N φ φ x n t T t y t P t N t cc t cc gyp t gyp P t P N t N t t t t cc t gyp , , , , , =0 0 0 t t t cc t gyp t t+1 +1 ⎡⎣ ⎤⎦∑ ≥ ∈ ∈ (6) where β ρ= (1 + )−1 is a social discount factor, and ρ 0≥ is a discount rate. We follow a common tradition in economic research of using an infinite‐horizon problem formulation: T → ∞. A soil management problem considering long‐term biophysical and economic con- sequences is typically formalized as a profit‐maximization model with an infinite time horizon (Hediger, 2003).4 Constant exogenous prices for barley yield, P fertilizer, and N fertilizer (€/kg) 3Even nonlegume catch crops may improve biological N fixation, and therefore affect the N use efficiency (Lemola et al., 2000). 4Stationary problems with an infinite time horizon lead to simplified optimal stationary strategies because there is no need to specify what happens after the terminal state. This approach is appropriate if the distant future has no significant influence on the optimal path for the near future (Sydsæter et al., 2005). In here, discounting takes care of the influence of the distant future on the decisions in the beginning of the planning horizon. 6 of 40 | Natural Resource Modeling SIHVONEN ET AL. are denoted by py, pP, and pN , respectively. Constant exogenous costs for catch crop and gypsum treatment (€/ha) are denoted by pcc and pgyp. Parameters μP and μN are the constant exogenous marginal damage costs of P and N losses, respectively. The marginal damage costs reflect the average monetized marginal change in the benefits that people derive from the aquatic ecosystems when nutrient emissions decrease marginally. These costs are necessary for monetarizing the damages of N and P losses. 2.3 | Optimality condition for the annual input decisions Society's Lagrangean is ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) L β p f P N φ x n p P p N p φ p φ μ e P φ x μ e N φ n βλ x x P f P N φ x n x βλ n ϕ n N f P N φ x n e n γ φ γ φ γ φ γ φ = , , , , − + + + − , , − , , + + ϑ , , , , , , − + + , , , , , , , − + + + 1 − + 1 − . t t y t t t cc t t P t N t cc t cc gyp t gyp P P t t gyp t N N t t cc t t x t t t t t t cc t t t t n t t t t t t cc t t t N t t cc t cc t gyp t t cc t cc t gyp t gyp =0 +1 +1 +1 +1 ,0 ,0 gyp ,1 ,1 ⎜ ⎟ ⎛ ⎝ ⎡⎣ ⎤⎦ ⎡ ⎣⎢ ⎤ ⎦⎥ ⎡ ⎣⎢ ⎤ ⎦⎥ ⎞ ⎠ ∑∞ (7) Here λtx+1 and λtn+1 are the shadow prices of the plant available soil P and N stocks, re- spectively. Lagrange multipliers γti0 and γti1 are respectively connected to the lower and upper limits of catch crop and gypsum decisions, φti, where i cc gyp{ , }∈ . We obtain the first‐order conditions (multiplied with β t− ) for an interior maximum by differentiating the Lagrangean with respect to the input choices and the soil P and N stocks, because, by the assumption, the yield response function is concave. Annual P fertilizer rate is at its optimal level when the marginal value product of P fertilizer and the discounted marginal value product of soil P and N carryovers is equal to the marginal cost of a P fertilizer and the externality of the P fertilization: P p f β λ f λ ϕ f p μ e: + (ϑ + ϑ ) + = + .t y P t tx P t f t P t tn f t P t P P P tP, +1 , , , +1 , , ,⎡⎣ ⎤⎦ (8) Here subscripts denote partial derivatives. Thus, the optimal P rate is an increasing function of a P yield response and a decreasing function of the effect of P fertilization on P loss and the marginal damage cost of P loss. Optimal P rate is an increasing function of the shadow price of a soil P stock. Instead, optimal P rate is a decreasing function of the shadow price of a soil N stock, because the effect of a P rate on a soil N stock is negative through the increased N yield uptake. The condition for the optimal N fertilizer rate is N p f β λ f λ ϕ ϕ f p μ e: + ϑ + ( + ) = + .t y N t tx f t N t tn N t f t N t N N N tN, +1 , , +1 , , , ,⎡⎣ ⎤⎦ (9) In this case, the discounted marginal value product of a soil P carryover is negative, because N fertilization has a negative effect on a soil P stock by increased P yield uptake. SIHVONEN ET AL. Natural Resource Modeling | 7 of 40 Also, in this case the discounted marginal value product of a soil N carryover is positive. The condition for catch crop is similar to those for the N and P fertilizer rates, with the exception that a yield reduction implies negative marginal value product, and the abate- ment measure's externality effect is negative also: φ μ e β λ f λ ϕ f γ p γ p f:− + ϑ + + = + − .tcc N φ tN tx f t φ t tn f t φ t tcc cc tcc y φ t, +1 , , +1 , , ,0 ,1 ,cc cc cc cc⎡⎣ ⎤⎦ (10) This implies that optimal catch crop cultivation decision is an increasing function of the impact of a catch crop on N loss, and the marginal damage cost of N loss, and a decreasing function of the impact of a catch crop on an annual yield. Moreover, the effect of a catch crop on P and N carryovers is positive because of the reduced crop‐uptake effects.5 Gypsum application is at its optimal level when the negative externality of a gypsum application is equal to the marginal cost of the gypsum application: φ μ e γ p γ:− + = + .tgyp P φ tP tgyp gyp tgyp, ,0 ,1gyp (11) Gypsum is applied only when a soil P level is above a certain soil texture‐specific threshold value, denoted by xˆgyp. When x x> ˆt gyp, it holds that gypsum price is lower than its marginal environment benefit p μ e<φ P φ tP,gyp , and gypsum is applied. The condition for the optimal plant available soil P stock equalizes the sum of the effects of soil P on P and N carryovers and the marginal value product of a soil P stock with the current shadow price of the soil P stock and the marginal damage cost of a soil P stock: x p f β λ f λ ϕ f μ e λ: + (1 + ϑ + ϑ ) + = + .t y x t tx x t f t x t tn f t x t P x tP tx, +1 , , , +1 , , ,⎡⎣ ⎤⎦ (12) Similarly, the optimality condition for the plant available soil N stock is ( )n p f β λ ϕ ϕ f ϕ e λ f μ e λ: + 1 + + + + ϑ = + .t y n t tn n t f t n t e t n tN tx f t n t N n tN tn, +1 , , , , , +1 , , ,N⎡⎣ ⎤⎦ (13) We also have transversality conditions to guarantee that asymptotically the shadow prices of the soil stocks are zero: β λ xlim = 0 t t t x t+1→∞ and β λ nlim = 0 t t t n t+1→∞ . 2.4 | Steady state If the system is in a steady state, the optimal input decisions remain unchanged from one period to another, and therefore we can drop time indexes from Equations (8)–(13). In such a case, we can derive the symmetric steady‐state soil P and N shadow prices λx and λn (see A1 for the detailed derivation): 5If we would consider a legume catch crop, instead of a nonlegume one, there would also be a direct positive effect on the plant available N stock. 8 of 40 | Natural Resource Modeling SIHVONEN ET AL. ( ) λ p f μ e ϕ f ρ f = − + − (ϑ + ϑ ) − x y x P x P p f μ e ρ ϕ ϕ f ϕ e f x ρ x f x f ϕ f ρ ϕ ϕ f ϕ e − − ( + + ) 1 1 + ϑ − ( + + ) y n N n N n f n eN n N f n f x n f n eN n N ⎡ ⎣⎢ ⎤ ⎦⎥ (14) ( ) λ p f μ e f ρ ϕ ϕ f ϕ e = − + ϑ − + + − ,n y n N n N p f μ e ρ f f n ρ n f n e n N ϕ f f ρ f − − (ϑ + ϑ ) 1 1 + ϑ − (ϑ + ϑ ) y x P x P x f x N f x f n x f x ⎡ ⎣⎢ ⎤ ⎦⎥ (15) where p fy x is the marginal value product of a soil P stock, p fy n is the marginal value product of a soil N stock, μ eP xP is the marginal effect of the soil P stock on P loss externality, and μ eN nN is the marginal effect of the soil N stock on N loss externality. The third term in equation (14), ( )ϕ f p f μ e ρ ϕ ϕ f ϕ e f x − − + + y n N n N n f n eN n N , captures the marginal effect of the soil P stock on the social net benefit of the soil N stock, which is scaled down with the difference between the discount rate and the marginal effect of the soil N stock on the N carryover. Note that the term ϕ ff x captures the negative effect of the soil P stock on N carryover (through the crop uptake). Thus, if the net effect of the soil N stock is positive, the third term in Equation (14) decreases the steady‐state shadow value of the soil P stock. The interpretation is analogous for the third term in Equation (15). Hence, Equations (14) and (15) show that in a steady state social planner's shadow prices of the soil nutrient stocks are defined by the net periodic marginal social net benefits of the soil P and N stocks. These net benefits are divided with the so‐called effective discount factors. The effectiveness means that when discounting the perpetual effect of small addition of x or n, we do not use mere ρ, but the discount rate is adjusted with the decay rates of P and N stocks, that is, fϑ + ϑx f x and ϕ ϕ f ϕ e+ +n f n e nNN , and the interaction of the stocks, that is, P stock affects the decay of N stock and vice versa. In the case of P stock, we have soil N effects, f ϕ fϑf n f x, and their perpetual value is scaled down with the term ρ ϕ ϕ f ϕ e− ( + + )n f n e nNN , which again is the discount rate that is adjusted with the decay rate of a N stock. The interpretation is analogous for the effective discount rate of the soil N stock. Hence, these third terms are related to fact that P change affects N stock, which again affects Pstock. Note, that the factor ρ1/(1 + ) (or β) stems from the discrete time approach. This is a technical correction, which would not exist in a continuous time model. Note that if ϑ < 0x and ϕ < 0n , the discounting is in fact stronger than without the decay. In other words, the effective discount rate is typically greater than the actual discount rate. This is natural as the effect declines as the stock effect vanishes due to the decay. Soil P increases social welfare by increasing the yields, but it also has a negative effect on social welfare through increased P loss. Furthermore, soil P has a negative effect on a soil N stock via increased N yield‐uptake. Analogously, soil N has a positive effect on yield, negative environ- mental effect of increased N loss, and a negative effect on a soil P stock. Therefore, assuming that the effective discount rates are positive, the shadow prices are positive (negative) if the benefit from increased profits of a crop yield are greater (smaller) than the environmental externalities and the effects on the other nutrient stock. When λ > (<)0, the NPV is an in- creasing (decreasing) function of the nutrient stock. The result implies that the social value of a soil nutrient stock depends on the relative size of the long‐term costs and benefits of the soil nutrient stock. SIHVONEN ET AL. Natural Resource Modeling | 9 of 40 3 | ENVIRONMENTAL POLICY 3.1 | First‐best regulation We make a standard assumption that under the conditions of no regulations, landowners maximize the NPV of the revenues from crop production, and the environmental impacts of agricultural practices do not influence their decisions (Goetz & Zilberman, 2000; Tilman et al., 2002).6 There are two basic strategies for internalizing the externalities of food produc- tion: regulation and market incentives. We study agri‐environmental tax‐subsidy‐schemes, because they allow producers to respond to changes in economic conditions. A planner cannot tax nonpoint source emissions from agricultural land because those are difficult to measure and monitor. Instead, a regulator can alter measurable management choices causing the ex- ternalities. In addition, in practice it could be difficult to set a tax on soil nutrients. Thus, we focus on fertilizer input taxes, and subsidies on the abatement measures. However, a regulator needs to take into account the fact that P and N fertilization will cause lasting gradual emissions through the soil nutrient stocks. Therefore, a regulator needs to adjust fertilizer taxes and abatement subsidies. Under the conditions of perfect information, this tax‐subsidy‐scheme is a first‐best policy instrument resulting in the socially optimal nutrient pollution levels. To derive the tax‐subsidy‐scheme, we must first define the problem of a landowner under the environmental policy. Regulated private producer maximizes the NPV of annual net returns from crop production: { } ( ) ( ) ( ) ( ) { } β p y p τ P p τ N p σ φ p σ φ x n x P N φ φ max − + + + + − + − , Subject to equations (1) and (2), and are given, , , 0, {0, 1}, and {0, 1} P N φ φ x n t T t y t P t P t N t N t cc t cc t cc gyp t gyp t gyp t t t t cc t gyp , , , , , =0 0 0 t t t cc t gyp t t+1 +1 ⎡⎣ ⎤⎦ ∑ ≥ ∈ ∈ (16) Again, for simplicity, the time horizon is set to infinite: T → ∞.7 The first‐order conditions for the annual optimal decisions in the regulated private optimum are similar to those in the social optimum, with the exception that the respective taxes replace the marginal damages of the inputs (see A2). By equalizing the first‐order conditions for the social and regulated private optimums, and solving the equations for the respected taxes and subsidies, we obtain the first‐ best tax‐subsidy‐scheme (see A3 for details): τ μ e β f ϕ f= + (ϑ + ϑ )Δ + Δ ,tP P P tP P t f t P t tx f t P t tn, , , , +1 , , +1⎡⎣ ⎤⎦ (17) τ μ e β f ϕ ϕ f= + ϑ Δ + ( + )Δ ,tN N N tN f t N t tx N t f t N t tn, , , +1 , , , +1⎡⎣ ⎤⎦ (18) 6This is a simplification. In reality, a private producer may consider also the environmental effects of the production, at least to some extent, in addition to profits from the crop production. 7If the farmers are landowners, they can be long‐term profit maximizers, because the farmer can sell the land at value that corresponds to its future expected productive value. 10 of 40 | Natural Resource Modeling SIHVONEN ET AL. σ μ e β f ϕ f= − − ϑ Δ + Δ ,tcc N φ tN f t φ t tx f t φ t tn, , , +1 , , +1cc cc cc⎡⎣ ⎤⎦ (19) σ μ e= − .t P φ tPgyp ,gyp (20) where μ β a e μ β b eΔ = + ,tx P s s s x t s P N s s s n t s N +1 =0 , +1+ =1 , +1+∑ ∑ ∞ ∞ (21) μ β c e μ β d eΔ = + ,tn N s s s n t s N P s s s x t s P +1 =0 , +1+ =1 , +1+∑ ∑ ∞ ∞ (22) where parameters as, bs, cs, and ds are defined recursively as follows: a = 10 , a = ΘtP1 +1, and b = ΘtPN1 +1, and the rest are recursively defined by relations a a b= Θ + Θs s t sP s t sNP−1 + −1 + and b b a= Θ + Θs s t sN s t sPN−1 + −1 + . Similarly for nitrogen, c = 10 , c = ΘtN1 +1, and d = ΘtNP1 +1, and the rest are recursively defined by relations c c d= Θ + Θs s t sN s t sPN−1 + −1 + and d d c= Θ + Θs s t sP s t sNP−1 + −1 + . The transition factors are defined as fΘ 1 + ϑ + ϑtP x t f t x t+1 , +1 , +1 , +1≔ (the effect of the soil P stock on P carryover), ϕ fΘtPN f t x t+1 , +1 , +1≔ (the effect of the soil P stock on N carryover), ϕ ϕ f ϕ eΘ 1 + + +tN n t f t n t e t n tN+1 , +1 , +1 , +1 , +1 , +1N≔ (the effect of the soil N stock on N carryover), and fΘ ϑtNP f t n t+1 , +1 , +1≔ (the effect of the soil N stock on P carryover). Here, Δtx+1 and Δtn+1 define the social costs of phosphorus and nitrogen stocks, that is, the externalities of the legacy nutrient stocks from the current period to infinity. The regulator adjusts the taxes and the subsidies by their effect on soil N and P, based on decisions’ impact on future soil N and P losses. Equations (17)–(20) show that the first‐best taxes and subsidies take into account the direct emissions, discounted next‐period soil nutrient emissions, and the marginal costs of the environmental externality of the soil nutrient stocks from the current period onwards. The influence of the soil nutrient externalities decreases in time, because of discounting. Soil nutrient externalities are weighted with the effect of a given input on a soil nutrient carryover. Note that Equation (20) describing the subsidy for gypsum has only one term because gypsum does not affect the soil nutrient dynamics in the model. Equation (20) shows that the producer uses gypsum only if the gypsum subsidy exceeds the price of gypsum. A first‐best tax‐subsidy‐scheme yields zero social welfare loss. 3.2 | Static tax‐subsidy‐schemes The first‐best tax‐subsidy‐scheme requires information about the soil nutrient levels and the soil texture of the target field hectare because nutrient losses depend on a soil texture.8 Since the optimal tax is time‐dependent, it requires constant changes. This might not be feasible due to higher administrational costs. Also, calculating the future externalities of the legacy nutrient 8This information is available in Finland for those farms that are part of the Finnish Agri‐Environmental Program, FAEP (over 90% of the farms). FAEP requires farmers to test their fields for soil fertility. Farmers must renew the soil fertility tests in every 5 years. Samples for the tests must be taken for every field parcel that is larger than 0.5 hectare. Peltovuori (1999). The soil fertility tests are necessary for the precision farming (Iho & Laukkanen, 2012a). However, the soil tests in Finland typically gather information about soil P levels, while ignoring soil N levels, since soil N information is considered less relevant (Peltovuori, 1999). SIHVONEN ET AL. Natural Resource Modeling | 11 of 40 stores may be difficult. Thus, a planner may have to rely on suboptimal tax‐subsidy‐schemes approximating the first best policy. Here we focus on static schemes based on the steady state equilibrium. We examine a scheme that can differentiate between soil textures and a scheme that cannot. In a steady state, the tax‐subsidy‐scheme is time‐invariant. Therefore, we can drop time indexes from Equations (17)–(20). This leads to steady‐state‐based taxes and subsidies (see A4 for the detailed derivation): ( )τ μ e f ϕ f= + Φ ϑ + ϑ + Φ ,P P PP P P f P N f P⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ (23) ( )τ μ e f ϕ ϕ f= + Φ ϑ + Φ + ,N N NN P f N N N f N⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ (24) σ μ e f ϕ f= + Φ ϑ + Φ ,cc N φN P f φ N f φ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎cc cc cc (25) σ μ e= ,P φPgyp ⁎gyp (26) where ( ) ( ) ( ) μ e ϕ f ρ f Φ = + − ϑ + ϑ − ,P P x P μ e ρ ϕ ϕ f ϕ e f x x f x f ϕ f ρ ϕ ϕ f ϕ e ⁎ − + + ⁎ ⁎ ⁎ ⁎ ⁎ ϑ − + + N n N n f n eN n N f n f x n f n eN n N ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ (27) ( ) ( ) ( ) μ e f ρ ϕ ϕ f ϕ e Φ = + ϑ − + + − .N N n N μ e ρ f f n n f n e n N ϕ f f ρ f ⁎ − ϑ + ϑ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ϑ − ϑ + ϑ P x P x f x N f x f n x f x ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ ⁎ (28) Here the terms ΦP and ΦN are the differences between private and social shadow prices of the soil P and N stocks in the steady state, respectively. Equations (23)–(26) show that the socially optimal steady‐state tax‐subsidy‐scheme captures the direct effect of the input choices on the associated nutrient loss and the indirect effects on soil nutrient accumu- lation/depletion. The latter effect is scaled with the difference between private and social shadow prices of the nutrient stocks in the steady state. The steady‐state‐based tax‐ subsidy‐scheme is suboptimal for those periods where the soil nutrient stocks differ from their steady‐state levels. Therefore, it generates social welfare loss. However, the design of the steady‐state‐based policy instruments requires less information than the first‐best policy instrument: the only information needed is the soil texture of a target field hectare. Hence, administration costs should be lower, which compensates the environmental welfare losses. Finally, we consider also a situation where a regulator does not utilize the information of the soil texture of a target hectare either. Such a scheme can be determined, for example as an weighted average of the soil texture‐specific steady‐state taxes and subsidies: τ ω τ=i jn j ji∑ , where ωj is a share of a given soil texture j of the whole spectrum in the target area. This tax‐ subsidy‐scheme requires only information on the distribution of soil textures in the target area. Hence, the information requirements and the resulting administration costs are even lower. 12 of 40 | Natural Resource Modeling SIHVONEN ET AL. However, the environmental welfare loss of this tax‐subsidy‐scheme is necessarily greater than that of the policy instrument utilizing the soil texture information. 4 | QUANTITATIVE ASSESSMENT 4.1 | The setup Our numerical demonstration uses data and models estimated and validated for spring barley cultivation on clay and coarse soil textures in Finland (see A5). These two soil textures are essentially different in their characteristics: clay soils have greater ability to hold nutrients compared to that on coarse soils (Goetz, 1997; Hazelton & Murphy, 2007). This holds particularly on Finnish young glacial clay soils (Iho & Laukkanen, 2012b). Cultivation technology is conventional tillage. Table 1 shows the price parameters. Applied social discount rate is 3%. The choice of the discount rate is crucial for a long‐term analysis (Weitzman, 2001, 1998), and therefore we studied the sensitivity of the results to this choice (see S1). 4.2 | Soil P and N paths, annual input choices, yields, and nutrient losses in social and private optimums We first study the development of soil P and N stocks in the social and private optimums. The reduction time paths for high legacy P and N stores, x0 and n0 are almost the same in the private and social optimums on both soil textures (Figure 2). In the private optimum, high legacy stores TABLE 1 Price parameter estimates applied in the empirical model Parameter Estimate Source Barley price 0.135 €/kg Kotimaan hinnat (www.vyr.fi) N fertilizer cost 0.91 €/kg Luke (2015) P fertilizer costb 1.99 €/kg Luke (2015) Gypsum amendment costb 44 €·ha−1·yr−1 Ollikainen et al. (2018) Catch crop cultivation cost 33.5 €·ha−1·yr−1 COWI (2007) Marginal damage of N loss 6.6 €/kg Gren and Folmer (2003) Marginal damage of P lossc 47 €/kg Gren and Folmer (2003) aThe assumption of the constant P price may not be reasonable because rock phosphate deposits are finite, and the depletion of those deposits may become an issue in the future (Berntsen et al., 2003). We studied the effect of the increasing P price on the results (S2). bThe annual cost is obtained by dividing total cost, 220 €/ha by five because gypsum is effective for 5 years after the application. cWe applied a Redfield ratio (C:N:P = 106:16:1) to transfer marginal damage estimate for N‐loss, μN to marginal damage estimate for P‐loss, μP (Redfield, 1934). When the ratio is viewed in grams, following weights are applied: C‐weight= 12 g/mol, N‐weight= 14 g/mol and P‐weight= 31 g/mol. When all the weights in the ratio are divided by 31, the ratio becomes 41:7.2:1, and the obtained estimate for μP is 47 €/kg (7.2 times 6.6 €/kg). SIHVONEN ET AL. Natural Resource Modeling | 13 of 40 are decreased, because it is not economically efficient to maintain such high nutrient stores.9 In the social optimum, high legacy stores are decreased also because of the high nutrient loss externalities. When x0 is low, private producer increases soil P level faster and to a higher level than what is a social optimal path. However, on clay soils the difference between privately and socially optimal time paths of the soil N stores is practically negligible. Since the run‐off externalities have no influence on the private input decisions, the optimal annual fertilization rates are notably higher in private optimum than in social optimum (S3). Hence, the fertilization differences are explained by the high environmental damage cost es- timates. These damage costs, however, are prominently uncertain (Keeler et al., 2016; Keiser et al., 2019), and the results are sensitive to these cost estimates (S4 and S5). The socially and privately optimal N rates are higher on clay soils than on coarse soils, because of a higher N yield response on clay soils, and because N loss is higher on coarse soils. Instead, socially and privately optimal P rate is higher on coarse soils than on clay soils, because of the higher P yield response on coarse soils. In addition, high P rate on coarse soils is needed to maintain a relatively high soil P level due to coarse soil's low ability to hold nutrients. 0 20 40 60 80 100 120 0 10 20 30 40 So ilP (m g/ l) Coarse soils 0 20 40 60 80 100 120 0 10 20 30 40 Clay soils Social, Soil P(1): 1 mg/l Social, Soil P(1): 36 mg/l Private, Soil P(1): 1 mg/l Private, Soil P(1): 36 mg/l 0 20 40 60 80 100 120 Time (years) 0 100 200 300 400 500 So ilN (kg /ha ) Coarse soils 0 20 40 60 80 100 120 Time (years) 0 100 200 300 400 500 Clay soils Social, Soil N(1): 0.5*SS-level Social, Soil N(1): 1.5*SS-level Private, Soil N(1): 0.5*SS-level Private, Soil N(1): 1.5*SS-level (a) (b) (c) (d) FIGURE 2 Socially and privately optimal soil P and N paths for low and high legacy nutrient stores. SS‐level stands for the steady‐state N stock, which was 250 kg/ha on clay soils and 100 kg/ha on coarse soils. In (a and b) the initial soil N stock is at its steady‐state level, and in (c and d) the initial soil P stock is at its steady‐state level, which is 3.23 and 8.5 mg/L on clay and coarse soils, respectively 9This suggests that the existence of the high legacy nutrient stores, particularly high P stores, is a result of the past inefficient P management (Svanbäck et al., 2019). There are many uncertain and stochastic external factors, however, which may have led to decisions causing too high legacy stores. In some cases, it may be reasonable for a risk averse farmer to accumulate soil nutrients to the levels, which are suboptimal according to a deterministic model. In addition, in the case of livestock production, the manure of the farm animals is spread on the field. Typically, transporting excess manure away from the field is too expensive, which may lead to increasing nutrient stocks. 14 of 40 | Natural Resource Modeling SIHVONEN ET AL. In the social optimum, a catch crop is cultivated throughout the planning horizon regardless of the legacy nutrient stores on both soil textures. This is because the negative effect of a catch crop on N loss enables higher N rates in a social optimum compared to the situation where catch crop cultivation is not possible. In addition, catch crop cultivation cost, 33.85 €/ha is below the threshold level where cultivation becomes optimal: 85 and 41 €/ha on coarse and clay soils, respectively. Catch crop is not cultivated in a private optimum, because private producer ignores the benefits of reduced N losses. Gypsum is applied in social optimum only temporarily on soils with high legacy P stores (S6). The threshold level of a soil P stock is lower on clay soils (16 mg/L) compared to that on coarse soils (32 mg/L). Gypsum is applied longer on clay soils than on coarse soils, because the threshold level is reached faster on coarse soils with low ability to hold nutrients. Gypsum is not applied in the private optimum. The lower fertilizer rates in the social optimum lead to lower annual yields compared to those in the private optimum (S7). In addition, annual nutrient losses are higher in the private optimum than in the social optimum (S7). When a legacy P stock, x0, is high, P loss is higher on clay soils than on coarse soils, because annual erosion is higher on clay soils (see A5). Annual N loss, on the other hand, is higher on coarse soils than on clay soils, in line with the previous literature (Coote & Ramsey, 1983; Nichols, 1984). 4.3 | Social NPVs Social NPV is, of course, higher in the social optimum than in the private optimum, because of the high social costs of the nutrient losses (Figure 3). The social welfare loss of the private optimum is greater on coarse soils than on clay soils, because of the high N losses on coarse soils. On both soils, a social welfare loss is an increasing function of the legacy nutrient stores. Catch crop has a substantial effect on social NPV on coarse soils (Figure 3a). Instead, catch crop's effect is relatively minor on clay soils (Figure 3b). The welfare effect of gypsum is only observed on soils with x0 above the threshold for gypsum application. Overall, environmental policies are necessary to reduce privately optimal fertilizer rates and to induce the private producers to use the catch crop and gypsum, when socially optimal. The social NPV is an increasing function of x0 on coarse soils and a decreasing function of x0 on clay soils (Figure 3a,b). On coarse soils, the benefit of increased soil fertility outweighs the damage of increased P loss, implying that the steady‐state shadow price of a soil P stock as given by Equation (14), is positive (75.5 €/ha).10 On clay soils the P loss externalities of an increased soil P stock outweighs the associated benefits from a crop production, and the steady‐state shadow price of a soil P stock is negative (−209 €/ha). Hence, the shadow price of a soil P stock is positive (negative) on soils with low (high) ability to hold nutrients. However, the steady‐state level of a soil P stock on clay soils is above zero, because crops need a certain amount of P to grow. On coarse soils, the steady‐state level of a soil P stock on coarse soils is moderate, because maintaining higher soil P stock generates costs in terms of both input costs and externality costs. Analogously, a social NPV is a decreasing function of n0 on coarse textured soils, where the steady‐state shadow price of a soil N stock as given by Equation (15), is negative (−0.30 €/ha), whereas a social NPV is an increasing function ofn0 on clay soils, where the steady‐state shadow price of a soil N stock is positive (0.03 €/ha). However, when catch crop 10We obtained the numeric values for the shadow prices by inserting steady‐state levels for the input choices and the soil nutrient stocks to the formulas shown in Equations (14) and (15). SIHVONEN ET AL. Natural Resource Modeling | 15 of 40 cultivation is not possible, the shadow price of a soil N stock is negative also on clay soils, and therefore a social NPV decreases when n0 increases. 4.4 | Policy instruments: Tax‐subsidy‐schemes As explained in Section 3.2, a social planner can implement first‐best tax‐subsidy‐schemes if she/he has information about the legacy nutrient stores and the soil textures of the target area, and capability to implement annually changing taxes and subsidies. In the case of a low legacy P store, x0, a socially optimal annual P‐tax, as described by Equation (17), is initially high, particularly on clay soils, because in such a case privately optimal P rate is high (Figure 4b). At the same time, privately optimal N rate is low, and therefore a socially optimal N‐tax, as described by Equation (18), is below its steady‐state level (Figure 4a). A social planner uses also nonlegume catch crops and gypsum to reduce nutrient losses. In this case, an optimal catch crop subsidy, as described by Equation (19), is initially below its steady‐state level because N loss is low in the private optimum (Figure 4c). At the same time, optimal gypsum subsidy, as described by Equation (20) is low, because gypsum is applied only on soils with high legacy P stores (Figure 4d). In the case of high x0, privately optimal N rate is initially high. As a result, annual N‐loss is high in the private optimum. Therefore, an optimal annual N‐tax is above its steady‐state level before the steady state is reached (Figure 4a). An optimal catch crop subsidy is also high because of high N loss (Figure 4c). Instead, in this case the privately optimal P rate is low, because plant obtains most of the needed P from the soil P stock. Therefore, an optimal annual P‐tax is below its steady‐state level (Figure 4b). However, on coarse soils, P‐tax is high also with 0 20 40 60 Initial soil P level (mg/l) 0 2000 4000 6000 8000 Coarse soils 0 20 40 60 Initial soil P level (mg/l) 0 2000 4000 6000 8000 Clay soils Social Private Social, no CC Social, no GYP Social, no CC & GYP 0 100 200 300 400 500 Initial soil N level (kg/ha) 0 2000 4000 6000 8000 Coarse soils 0 100 200 300 400 500 Initial soil N level (kg/ha) 0 2000 4000 6000 8000 Clay soils (a) (b) (c) (d) FIGURE 3 Social NPV as a function of the legacy nutrient stores for the socially and privately optimal input choices over the planning horizon, and for the cases where CC is not available, GYP is not available, and where neither measure is available. Note that in (c and d) the curves associated with the gypsum are missing, because initial soil P level is at its steady‐state level, which is below the threshold for gypsum application. CC, catch crop; GYP, gypsum; NPV, net present value 16 of 40 | Natural Resource Modeling SIHVONEN ET AL. 0 20 40 60 80 100 120 Time (year) 0 0.2 0.4 0.6 (a) 0 20 40 60 80 100 120 Time (year) 0 1 2 3 (b) Coarse, Soil P(1): 1 mg/l Coarse, Soil P(1): 48 mg/l Clay, Soil P(1): 1 mg/l Clay, Soil P(1): 48 mg/l 0 20 40 60 80 100 120 40 60 80 100 (c) 0 20 40 60 80 100 120 0 20 40 60 80 100 (d) 0 50 100 150 200 Time (year) 0.2 0.3 0.4 0.5 (e) 0 50 100 150 200 Time (year) 0 0.5 1 1.5 2 2.5 (f) Coarse, Soil N(1): 0.5*SS-level Coarse, Soil N(1): 1.5*SS-level Clay, Soil N(1): 0.5*SS-level Clay, Soil N(1): 1.5*SS-level 0 50 100 150 200 40 60 80 100 (g) 0 50 100 150 200 0 5 10 15 20 25 (h) FIGURE 4 First‐best tax‐subsidy‐schemes for low and high legacy nutrient stocks. In (a–d) initial soil N level is at its SS level (260 kg N ha−1 on clay soils and 100 kg N ha−1 on coarse soils), and in (e–h) initial soil P level is at its steady‐state level (3.23mg P L−1 on clay soils and 8.5 mg P L−1 on coarse soils). SS, steady‐state SIHVONEN ET AL. Natural Resource Modeling | 17 of 40 high x0, because P loss is resulting from the high optimal soil P stock. Therefore, P‐tax paths on coarse soils are almost identical despite the different legacy P stores (Figure 4b). If the N legacy store,n0, is low (high), N loss is below its steady‐state level, and therefore a N tax is also below (above) its steady‐state level (Figure 4e). Similarly, an optimal catch crop subsidy is above (below) its steady‐state level on soils with high (low) n0 (Figure 4g). The effect of n0 on P tax and gypsum subsidy is almost negligible, because the N legacy store's effect on P losses is almost non- existent. Independently of the legacy nutrient stores, N‐tax and catch crop subsidy are higher on coarse soils than on clay soils, because of the higher N losses. Similarly, P‐tax is higher and gypsum subsidy is lower on clay soils because of the higher P losses compared to those on coarse soils. In practice, a social planner may not be able to use an annually changing tax‐subsidy‐scheme. Instead, a planner may rely on simpler static policy instruments that are fixed to a certain level for a given number of years, such as EU programs. Steady‐state‐based static taxes and subsidies are rea- sonable approximations of the first best instruments because the system converges asymptotically to the steady state independently of the legacy nutrient stores. On coarse soils, a steady‐state‐based N‐tax is higher than that on clay soils because a socially optimal steady‐state N rate is higher on clay soils than on coarse soils (Table 2). Instead, on coarse soils a P‐tax is relatively low because the socially optimal steady‐state P rate is higher on coarse soils than on clay soils. Catch crop subsidy is higher on coarse soils than on clay soils because the benefit of catch crop cultivation on coarse soils is higher. Steady‐state‐based gypsum subsidy is lower than the marginal cost of the gypsum application on both soil textures. Consequently, gypsum is not applied in the regulated private optimum. Note, however, that the first‐best gypsum subsidy is above the marginal cost for those periods where the soil P is above the threshold level for the gypsum application (S6), and therefore, gypsum is used also in the regulated private optimum (Figure 4d). If the regulator cannot implement soil texture‐specific policies, she/he needs to simplify the tax‐subsidy‐scheme further. A soil texture‐ignorant static N‐tax is higher (lower) than the soil texture‐specific steady‐state N‐tax on clay (coarse) soils (Table 2). Instead, a soil texture‐ignorant P‐tax is lower (higher) than the soil texture specific steady‐state P‐tax on clay (coarse) soils. 4.5 | Welfare losses of the suboptimal policies The static steady‐state‐based tax‐subsidy‐scheme yields almost the same NPV as the first‐best instrument, unless the soil of a target area has a high ability to hold nutrients (clay soils) and TABLE 2 Steady‐state‐based tax‐subsidy‐schemes Steady‐state‐based soil texture‐ specific tax‐subsidy‐scheme Steady‐state‐based soil texture‐ ignorant tax‐subsidy‐schemea Soil texture Coarse Clay Coarse/Clay N‐tax (€/kg) 0.46 0.26 0.36 P‐tax (€/kg) 0.63 1.98 1.31 Catch crop subsidy (€/ha) 97.7 55.4 76.55 Gypsum subsidy (€/ha) 17.1 11.5 14.3 aA tax‐subsidy‐scheme which ignores the soil textures of the target area is based on the average of the soil texture‐specific steady‐state‐based taxes and subsidies. 18 of 40 | Natural Resource Modeling SIHVONEN ET AL. the legacy P store, x0 is considerably above its steady‐state level, x ⁎ (Figure 5). This is because on coarse soils the steady state is reached fast in the beginning of a planning horizon, in line with the findings of the previous literature (Jomini et al., 1991; Kennedy et al., 1973). If x0 is far from x ⁎, the static tax‐subsidy‐scheme generates welfare losses from a long period in the be- ginning of a planning horizon. Therefore, welfare loss is an increasing function of the difference between x0 and x ⁎. The main reason for the welfare loss is that gypsum is not used in the static tax‐subsidy‐schemes, because the subsidies are lower than the marginal cost of a gypsum application. Gypsum has a significant effect on a social NPV, especially on clay soils, when x0 is considerably above x ⁎, because it takes a longer time to reach the threshold soil P level, under which gypsum is no longer optimal, because of the naturally slower soil nutrient dynamics on clay soils than on coarse soils. Thus, the differences in the use of gypsum mostly explains the difference between the social NPVs associated with the first‐best tax‐subsidy‐schemes and the static tax‐subsidy‐schemes. For this reason, the welfare losses of the static tax‐subsidy‐schemes are almost nonexistent for different N legacy stores, if the initial soil P level is below the threshold level for the gypsum application (Figure 5c,d). Moreover, the increase in social welfare that is obtained by differentiating the static tax‐ subsidy‐scheme for different soil textures is almost negligible (Figure 5). This is because the soil‐texture‐ignorant fertilizer taxes for N and P change to opposite directions compared to the soil‐texture specific taxes. As a result, both N and P fertilizer rates adjust simultaneously to changed fertilizer taxes in such a way that the yields are almost the same in cases where soil‐ texture‐specific and soil‐texture‐ignorant tax‐subsidy‐schemes are used. Moreover, the NPVs are almost the same in these cases, independently of the legacy P stores, because in either of the cases gypsum is used. On clay soils, less N and more P is used when a soil‐texture‐ignorant scheme is used, compared to the situation where a soil‐texture‐specific scheme is used. This is 0 20 40 60 Initial soil P level (mg/l) 3000 4000 5000 6000 7000 Coarse soils 0 20 40 60 Initial soil P level (mg/l) 3000 4000 5000 6000 7000 Clay soils 0 100 200 300 400 500 Initial soil N level (kg/ha) 5000 6000 7000 Coarse soils 0 100 200 300 400 500 Initial soil N level (kg/ha) 5000 6000 7000 Clay soils First-best Static soil texture-specific Static soil texture-ignorant (a) (b) (c) (d) FIGURE 5 Social NPVs obtained by applying different tax‐subsidy‐schemes. When the x‐axis is the soil P level, the soil N level is fixed to its steady‐state level (260 kg N ha−1 on clay soils and 100 kg N ha−1 on coarse soils). Similarly, when the x‐axis is the soil N level, the soil P level is fixed to its steady‐state level (3.23mg P L−1 on clay soils and 8.5 mg P L−1 on coarse soils). NPV, net present value SIHVONEN ET AL. Natural Resource Modeling | 19 of 40 because soil‐texture‐ignorant N tax is higher, and P tax is lower compared to the corresponding taxes in the soil‐texture‐specific scheme. This leads to higher soil P levels and lower soil N levels in the case of a soil‐texture‐ignorant scheme, compared to those for a soil‐texture‐specific scheme (Figure 6b,c). On coarse soils soil P levels are lower and soil N levels higher in the case of a soil‐texture‐ignorant scheme, compared to those for a soil‐texture‐specific scheme (Figure 6f,g). However, the effect of the legacy P and N stocks on optimal P and N paths, respectively, is notably stronger than the effect of the different tax‐subsidy‐schemes (Figure 6a,d,e,h). These results highlight the benefits of the simultaneous adjustment of N and P nutrients. 5 | DISCUSSION AND CONCLUSIONS We developed a modeling framework, based on a discrete‐time dynamic optimization, for studying the effect of legacy nutrient stores and soil nutrient dynamics on an optimal nutrient management and agri‐environmental policy in crop production. The model suggests that high legacy nutrient stores are reduced to their steady‐state levels at similar rates in both the social and the private optima, if N and P fertilizer rates are simultaneously optimized on the economic basis of profit maximization. This is in line with the suggestion that a producer can recover the legacy P stores to reduce the reliance on synthetic fertilizers (Sharpley et al., 2013; Svanbäck et al., 2019). Moreover, according to our results, if the benefit from increased profits of a crop yield is greater (smaller) than the environmental externalities of the soil nutrient stocks, a social NPV is an increasing (decreasing) function of the legacy nutrient stock. Hence, the high legacy nutrient stock is not automatically a bad thing. Our results suggest that N loss abatement with catch crops is socially optimal. However, the welfare effect was notable only on coarse soils with high annual N loss. This supports a previous suggestion that there is a need for increased use of catch crops in Finland (Aronsson et al., 2016). The use of gypsum is well justified as a temporary measure on coarse soils with high soil P levels (cf., Iho and Laukkanen, 2012b). The optimal gypsum use could be more extensive if the model would consider its positive crop yield (Chen et al., 2008) and soil P development (Anderson et al., 1995) effects. Note also that abatement parameterization of gypsum was based on clay soil measurements (Ekholm et al., 2012), so the gypsum's impact on coarse soils may be less accurate in the model. Despite the similar reduction paths of the legacy nutrient stores in the private and social optimums, the private optimum was notably worse in terms of social welfare due to its high level of nutrient losses. Hence, a model was extended to study an optimal environmental policy. A dynamic tax‐subsidy‐scheme, which takes into account the current soil nutrient levels and field attributes such as soil texture, can enforce socially optimal outcome (cf. Xepapadeas, 1992). The practicality of such a policy instrument, however, depends on the local conditions, avail- ability of information, program administration costs, legal limitations, and the abilities of a regulatory agency (Horan & Ribaudo, 1999; Larson et al., 1996). However, it can be useful in the derivation of the simpler, static tax‐subsidy‐schemes, and for studying the welfare losses of the suboptimal policies. To examine the welfare losses of administratively more simple, suboptimal policies (Chukwuemeka, 2015; Claassen & Horan, 2001; Kolstad, 1987), we applied static policy schemes that were based on steady‐state shadow prices. These suboptimal policies did not induce gypsum use, which was the main source of their welfare losses. The soil texture could be ignored without a notable welfare loss because the joint optimization of N and P fertilizers could 20 of 40 | Natural Resource Modeling SIHVONEN ET AL. 0 20 40 60 80 100 Time (years) 0 10 20 30 40 So il P (m g/ l) Clay soils, high and low legacy soil P stock First-best policy First-best policy Static, soil-texture -specific policy Static, soil-texture -specific policy Static, soil-texture -ignorant policy Static, soil-texture -ignorant policy 0 20 40 60 80 100 Time (years) 240 245 250 255 260 So ilN (kg /h a) Clay soils, high and low legacy soil P stock 0 20 40 60 80 100 3 3.2 3.4 3.6 So ilP (m g/ l) Clay soils, high and low legacy soil N stock 0 20 40 60 80 100 100 200 300 400 So ilN (kg /h a) Clay soils, high and low legacy soil N stock 0 20 40 60 80 100 0 10 20 30 40 So il P (m g/ l) Coarse soils, high and low legacy soil Pstock 0 20 40 60 80 100 90 100 110 120 So ilN (kg /h a) Coarse soils, high and low legacy soil P stock 0 20 40 60 80 100 Time (years) 7 7.5 8 8.5 9 So ilP (m g/ l) Coarse soils, high and low legacy soil N stock 0 20 40 60 80 100 Time (years) 50 100 150 So ilN (kg /h a) Coarse soils, high and low legacy soil N stock (a) (b) (c) (d) (e) (f) (g) (h) FIGURE 6 Optimal paths for soil nutrient stocks under different tax‐subsidy‐schemes for high and low nutrient legacy stores. High legacy P store is 36 mg/L and low legacy P store is 1 mg/L, on both soil texture. High legacy N store is 50% above the steady‐state level and low legacy N store is 50% below the steady‐state level (steady‐state level is 260 kg N ha−1 on clay soils and 100 kg N ha−1 on coarse soils) SIHVONEN ET AL. Natural Resource Modeling | 21 of 40 compensate the opposite imperfect incentives of the scheme (cf., Nkonya & Featherstone, 2000; Schlegel et al., 1996). Hence, the simultaneous adjustment of the N and P fertilizer rates enables the use of simple, static and soil‐texture‐ignorant tax‐subsidy schemes. The low welfare losses in the case of low legacy nutrient stocks and high discount factor suggest that unified tax‐subsidy‐schemes may well be cost‐effective (see S1). Unified here means that a policy instrument is the same despite the initial soil nutrient levels or the soil textures of the site. This is in line with the previous literature (Helfand & House, 1995; Lötjönen et al., 2020; Shortle et al., 1998). Based on our results we can conclude that a mixed policy instrument, where the static tax‐subsidy‐schemes are used with a dynamic subsidy for gypsum application would lead to outcome with almost‐zero welfare loss. Our results show that the welfare losses of the static steady‐state‐based tax‐subsidy‐schemes are increasing functions of the legacy nutrient stores and soil's ability to hold nutrients, because the higher is this ability, the longer it takes to reach the steady state. The derived steady‐state‐based fertilizer taxes are in the same order of magnitude than the corresponding taxes applied in Finland and in EU (Hildén et al., 2007; Nemming & Hansen, 2015). However, relatively few countries have used fertilizer taxes to control nutrient loading from agricultural land: Finland (years 1976–1994), Holland, Sweden, Norway, Austria, and Denmark (Hildén et al., 2007) and some American states (Weersink et al., 1998). Economic incentive‐based policies are challenging to use when the source of the pollution and the effect of the abatement practices on the pollution are difficult or costly to monitor (Griffin & Bromley, 1982; Segerson, 1988; Shortle & Dunn, 1986). Instead, policymakers typically set environmental objectives and standards; command and control has dominated environmental policy in most countries (Weersink et al., 1998). In EU, the most important policies aiming to reduce nonpoint source pollution have been direct controls such as the Drinking Water Directive and the Nitrate Directive (Hanley, 2001). In addition, subsidies to abatement measures have been more common than taxes to polluting inputs. Here, the derived subsidy for catch crop cultivation on coarse soils corresponds to the current subsidy in Finland (100 €/ha), whereas the derived catch crop subsidy on clay soils corresponds to the subsidy before the year 2015 (52 €/ha; Aronsson et al., 2016). There are at least three natural ways to extend the model. First, and most importantly, large nutrient legacy stores are often in areas with high livestock density, and the spatial differ- entiation of crop and livestock production is often the underlying cause of the high nutrient legacies (Svanbäck et al., 2019). We could extend the model to include manure as an input, but then we should also consider a spatial dimension, because in such a case transport costs become essential future of the problem (Lötjönen et al., 2020; Sharara et al., 2017). Second, we ignored the impacts of climate change, which may have significant influence on agricultural pro- ductivity and the externalities on the long run (Gornall et al., 2010). Third, models including multiple crops are more accurate compared to the monoculture models, since the majority of the European farmers typically use crop rotation patterns (Thomas, 2003). Crop rotation sys- tems can reduce the dependence on fertilizer inputs through internal nutrient recycling, be beneficial for soil fertility and farm profitability, and reduce environmental externalities of the production (Guertal et al., 1997; Ikerd, 1991; Vandermeer et al., 1998). AUTHOR CONTRIBUTIONS Matti Sihvonen and Kari Hyytiäinen developed the original research idea. Matti Sihvonen carried out the coding, performed all the theoretical derivations and the numerical optimization calculations, simulations and estimations, produced the results (the figures and the tables), and wrote the manuscript. Jussi Lintunen provided assistance in the mathematical formulation of the equations, commented the manuscript, and wrote selected, short sections of the manuscript. Kari Hyytiäinen 22 of 40 | Natural Resource Modeling SIHVONEN ET AL. commented the manuscript, and wrote selected, short sections of the manuscript. Elena Valkama provided information about the nitrogen loss functions and commented the manuscript. ORCID Matti Sihvonen http://orcid.org/0000-0002-3603-1760 Jussi Lintunen https://orcid.org/0000-0002-4635-634X Elena Valkama https://orcid.org/0000-0002-8337-8070 Kari Hyytiäinen https://orcid.org/0000-0002-4366-0186 REFERENCES Ahlvik, L., Ekholm, P., Hyytiäinen, K., & Pitkänen, H. (2014). An economic‐ecological model to evaluate impacts of nutrient abatement in the Baltic Sea. Environmental Modelling & Software, 55, 164–175. Anderson, D. L., Tuovinen, O. H., Faber, A., & Ostrokowski, I. (1995). Use of soil amendments to reduce soluble phosphorus in dairy soils. Ecological Engineering, 5, 229–246. Aronsson, H., Hansen, E. M., Thomsen, I. K., Liu, J., Øgaard, A. F., Känkänen, H., & Ulén, B. (2016). The ability of cover crops to reduce nitrogen and phosphorus losses from arable land in southern Scandinavia and Finland. Journal of Soil and Water Conservation, 71, 41–55. Barrow, N. J. (1980). Evaluation and utilization of residual phosphorus in soils. In F. E. Khasawneh, E. C., Sample, E. J., Kamprath (Eds), The Role of Phosphorus in Agriculture, 1976 (pp. 333–359). American Society of Agronomy Crop Science Society of America Soil Science Society of America Madison. Bechmann, M., Eggestad, H. O., & Vagstad, N. (1998). Nitrogen balances and leaching in four agricultural catchments in southeastern Norway. Environmental Pollution, 102(Sl), 493–499. Berntsen, J., Petersen, B. M., Jacobsen, B. H., Olesen, J. E., & Hutchings, N. J. (2003). Evaluating nitrogen taxation scenarios using the dynamic whole farm simulation model FASSET. Agricultural Systems, 76, 817–839. Bäckman, S., & Lansink, A. O. (2005). Crop and soil specific N and P efficiency and productivity in Finland. Agricultural and Food Science, 14, 262–276. Carpenter, S., Caraco, C. N. F., Correll, D. L., Howarth, R. W., Sharpley, A. N., & Smith, V. H. (1998). Nonpoint pollution of surface waters with phosphorus and nitrogen. Ecological Applications, 8, 559–568. Carpenter, S. R. (2005). Eutrophication of aquatic ecosystems: Biostability and soil phosphorus. Proceedings of the National Academy of Sciences of the United States of America, 102(29), 10002–10005. Chen, L., Dick, W. A., & Kost, D. (2008). Flue gas desulfurization products as sulfur sources for corn. Soil Science Society of American Journal, 72, 1464–1470. Chukwuemeka, O. (2015). Managing an accumulative inorganic pollutant: an optimal tax prescription for the social planner. International Journal of Economics, Commerce and Management, 3(8), 1–14. Claassen, R., & Horan, R. (2001). Uniform and non‐uniform second‐best input taxes—‐The significance of market price effects on efficiency and equity. Environmental and Resource Economics, 19, 1–22. Coote, D. R., & Ramsey, J. F. (1983). Quantification of the effects of over 35 years of intensive cultivation on four soils. Canadian Journal of Soil Science, 63(1), 1–14. COWI. (2007). Economic analysis of the BSAP with focus on eutrophication. Final report. 112 pp. Dodd, J. R., & Mallarino, A. (2005). Soil‐test phosphorus and crop grain yield responses to long‐term phosphorus fertilization for corn‐soybean rotations. Soil Science Society of America Journal, 69, 1118–1128. Ekholm, P., Turtola, E., Grönroos, J., Seuri, P., & Ylivainio, K. (2005). Phosphorus loss from different farming systems estimated from soil surface phosphorus balance. Agriculture, Ecosystems and Environment, 110, 226–278. Ekholm, P., Valkama, P., Jaakkola, E., Kiirikki, M., Lahti, K., & Pietola, L. (2012). Gypsum amendment of soils reduces phosphorus losses in an agricultural catchment. Agricultural and Food Science, 21, 279–291. Goetz, R. (1997). Diversification in agricultural production: A dynamic model of optimal cropping to manage soil erosion. American Journal of Agricultural Economics, 79(2), 341–356. Goetz, R., & Zilberman, D. (2000). The dynamics of spatial pollution: The case of phosphorus runoff from agricultural land. Journal of Economic Dynamics & Control, 24, 143–163. SIHVONEN ET AL. Natural Resource Modeling | 23 of 40 Gornall, J., Betts, R., Burke, E., Clark, R., Camp, J., Willett, K., & Wiltshire, A. (2010). Implications of climate change for agricultural productivity in the early twenty‐first century. Philosophical Transactions of the Royal Society of London, Series B, Biological Sciences, 365(1554), 2973–2989. Gren, I. M., & Folmer, H. (2003). Cooperation with respect to cleaning of an international water body with stochastic environmental damage: The case of the Baltic Sea. Ecological Economics, 47, 33–42. Gren, I. M., Savchuk, O. P., & Jansson, T. (2013). Cost‐effective spatial and dynamic management of a eutrophied Baltic Sea. Marine Resource Economics, 28, 263–284. Griffin, R. C., & Bromley, D. W. (1982). Agricultural runoff as a nonpoint externality: A theoretical development. American Journal of Agricultural Economics, 64, 547–552. Guertal, E., Bauske, M., & Edwards, J. (1997). Crop rotation effects on sweet potato yield and quality. Journal of Production Agriculture, 10(1), 70–73. Hanley, N. (2001). Policy on agricultural pollution in the European Union. In J. S. Shortle & D. Abler (Eds.), Environmental policies for agricultural pollution control (p. 224). Cabi Publishing. Hazelton, P. A., & Murphy, B. W. (2007). Interpreting soil test results: What do all the numbers mean? (2nd ed., p. 152). Collingwood, VIC: CSIRO Publishing. Hediger, W. (2003). Sustainable farm income in the presence of soil erosion: An agricultural Hartwick rule. Ecological Economics, 45, 221–236. HELCOM (2018). Sources and pathways of nutrients to the Baltic Sea. In Baltic Sea Environment Proceedings, 153. Helfand, G. E., & House, B. W. (1995). Regulating nonpoint source pollution under heterogeneous conditions. American Journal of Agricultural Economics, 77, 1024–1032. Hildén, M., Huhtala, A., Koikkalainen, K., Ojanen, M., Grönroos, J., Helin, J., Isolahti, M., Kaljonen, M., Kangas, A., Känkänen, H., Puustinen, M., Salo, T., Turtola, E., & Uusitalo, R. (2007). Verotukseen perustuva ohjaus maatalouden ravinnepäästöjen rajoittamisessa. Ympäristöministeriön raportteja, 15. Retrieved from https:// julkaisut.valtioneuvosto.fi/bitstream/handle/10138/41355/YMra_15_Verotukseen_perustuva_ohjaus_maat- alouden_ravinnepaastojen_rajoittamisessa.pdf?sequence=2 Horan, R. D., & Ribaudo, M. O. (1999). Policy objectives and economic incentives for controlling agricultural sources of nonpoint pollution. Journal of the American Water Resources Association, 35(5), 1023–1035. Huttunen, I., Huttunen, M., Piirainen, V., Korppoo, M., Lepistö, A., Räike, A., Tattari, S., & Vehviläinen, B. (2016). A national scale nutrient loading model for Finnish watersheds – VEMALA. Environmental Modelling and Assessment, 21(1), 83–109. Hyytiäinen, K., Niemi, J., Koikkalainen, K., Palosuo, T., & Salo, T. (2011). Adaptive optimization of crop production and nitrogen leaching abatement under yield uncertainty. Agricultural Systems, 104, 634–644. Iho, A. (2010). Spatially optimal steady‐state phosphorus policies in crop production. European Review of Agricultural Economics, 37(2), 187–208. Iho, A., & Laukkanen, M. (2012a). Precision phosphorus management and agricultural phosphorus loading. Ecological Economics, 77, 91–102. Iho, A., & Laukkanen, M. (2012b). Gypsum amendment as a means to reduce agricultural phosphorus loading: An economic appraisal. Agricultural and Food Science, 21, 307–324. Ikerd, J. (1991). A decision support system for sustainable farming. Northeastern Journal of Agricultural Economics, 20(1), 109–113. Jansson, P. E. (2012). Coup model: Model use, calibration, and validation. Transactions of the ASABE, 55(4), 1337–1346. Jomini, P. A., Deusin, R. R., & Lowenberg‐DeBoer, J. (1991). Modelling stochastic crop response to fertilization when carry‐over matters. Agricultural Economics, 6, 97–113. Keeler, B. L., Gourevitch, J. D., Polansky, S., Isbell, F., Tessum, C. W., Hill, J. D., & Marshall, J. D. (2016). The social costs of nitrogen. Science Advances, 2, 1–9. Keiser, D. A., Catherine, L. K., & Shapiro, J. S. (2019). The low but uncertain measured benefits of US water quality policy. Proceedings of the National Academy of Sciences of the United States of America, 116(12), 5262–5269. Kennedy, J. O. S., Whan, I. F., Jackson, R., & Dillon, J. L. (1973). Optimal fertilizer carry‐over and crop recycling policies for a tropical grain crop. Australian Journal of Agricultural Economics, 147, 104–113. 24 of 40 | Natural Resource Modeling SIHVONEN ET AL. Kleinman, P. J. A., Fanelli, R. M., Hirsch, R. M., Buda, A. R., Easton, Z. M., Wainger, L. A., … Shenk, G. W. (2019). Phosphorus and the Chesapeake Bay: Lingering issues and emerging concerns for agriculture. Journal of Environmental Quality, 48(5), 1191–1203. Kolstad, C. D. (1987). Uniformity versus differentiation in regulating externalities. Journal of Environmental Economics and Management, 14, 386–399. Kotimaan hinnat (2020). Markkinatietoa kotimaasta, EU:sta ja maailmalta. Retrieved from https://www.vyr.fi/ fin/markkinatietoa/kotimaan-hinnat/ Lambert, D. M., Lowenberg‐DeBoer, J., & Malzer, G. (2007). Managing phosphorus soil dynamics over space and time. Agricultural Economics, 37, 43–53. Larson, D. M., Helfand, G. E., & House, B. W. (1996). Second‐best tax policies to reduce nonpoint source pollution. American Journal of Agricultural Economics, 78, 1108–1117. Lemola, R., Turtola, E., & Christian, E. (2000). Undersowing Italian ryegrass diminishes nitrogen leaching from spring barley. Agricultural and Food Science in Finland, 9, 201–215. Liu, J., Ulen, B., Bergkvist, G., & Aronsson, H. (2014). Freezing‐thawing effects on phosphorus leaching from catch crops. Nutrient Cycling in Agroecosystems, 99, 17–30. Luke. (2015). Kasper: Ajankohtaista tietoa pelto‐ja puutarhaviljelystä sekä kasvinsuojelusta: Fosforilaskuri. Retrieved from https://portal.mtt.fi/portal/page/portal/kasper/pelto/peltopalvelut/fosforilaskuri Lötjönen, S., Temmes, E., & Ollikainen, M. (2020). Dairy farm management when nutrient runoff and climate emissions count. American Journal of Agricultural Economics. https://doi.org/10.1002/ajae.12003 Martıńez, Y., & Albiac, J. (2004). Agricultural pollution control under Spanish and European environmental policies. Water Resources Research, 40, 1–12. McDowell, R. W., & Sharpley, A. N. (2001). Approximating phosphorus release from soils to surface runoff and surface drainage. Journal of Environmental Quality, 30, 508–520. Van Meter, K. J., Van Cappellen, P., & Basu, N. B. (2018). Legacy nitrogen may prevent achievement of water quality goals in the Gulf of Mexico. Science, 360, 427–430. Nemming, A., & Hansen, R. V. (2015). Fertilizer accounts in Denmark. Ministry of Food, Agriculture and Fisheries of Denmark. HELCOM workshop, Oldenburg 6. Maj. Retrieved from https://helcom.fi/media/ documents/01_HELCOM-Workshop_Denmark.pdf Nichols, J. D. (1984). Relation of organic carbon to soil properties and climate in the southern Great Plains. Soil Science Society of America Journal, 48, 1382–1384. Nkonya, E. M., & Featherstone, A. M. (2000). Determining socially optimal nitrogen application rates using a delayed response model: The case of irrigated corn in western Kansas. Journal of Agricultural and Resource Economics, 25(2), 453–467. Ollikainen, M., Ekholm, P., Punttila, E., Ala‐Harja, V., Riihimäki, J., Puroila, S., Kosenius, A.‐K., & Iho, A. (2018). Peltojen kipsikäsittely maatalouden vesiensuojelukeinona. Retrieved from https://blogs.helsinki.fi/ save-kipsihanke/files/2018/11/SAVE-hankkeen-tietopaketti-kipsik%C3%A4sittelyst%C3%A4.pdf Park, S. C., Stoecker, A. L., Hattey, J. A., & Turner, J. C. (2007). Long‐term profitability of animal manure using optimal nitrogen application rate. 2007 Annual Meeting, Southern Agricultural Economics Association, Mobile, AL, February 4–7. Retrieved from http://ageconsearch.umn.edu/record/34830/files/sp07pa04.pdf Peltovuori, T. (1999). Precision of commercial soil testing practice for phosphorus fertilizer recommendations in Finland. Agricultural and Food Science, 8(3), 299–308. Pretty, J., Brett, C., Gee, D., Hine, R., Mason, C., Morison, J., Rayment, M., Van Der Bijl, G., & Dobbs, T. (2001). Policy challenges and priorities for internalizing the externalities of modern agriculture. Journal of Environmental and Management, 44(2), 263–283. Rankinen, K., Salo, T., Granlund, K., & Rita, H. (2007). Simulated nitrogen leaching, nitrogen mass field balances and their correlation on four farms in south‐western Finland during the period 2000–2005. Agricultural and food science, 16, 387–406. Redfield, A. C. (1934). On the proportions of organic derivatives in seawater and their relation to the composition of plankton. In R. J. Daniel (Ed.), James Johnston Memorial Volume (pp. 176–192). University Press of Liverpool. Räike, A., Taskinen, A., & Knuuttila, S. (2019). Nutrient export from Finnish rivers into the Baltic Sea has not decreased despite water protection measures. Ambio, 49, 460–474. SIHVONEN ET AL. Natural Resource Modeling | 25 of 40 Saarela, I. (2002). Phosphorus in Finnish soils in the 1900s with particular reference to the acid ammonium acetate soil test. Agricultural and Food Science, 11(4), 257–271. Salo, T. J., Palosuo, T., Kersebaum, K. C., Nendel, C., Angulo, C., Ewert, F., … Rötter, R. P. (2016). Comparing the performance of 11 crop simulation models in predicting yield response to nitrogen fertilization. The Journal of Agricultural Science, 154, 1218–1240. Schlegel, A. J., Dhuyvetter, K. C., & Havlin, J. L. (1996). Economic and environmental impacts of long‐term nitrogen and phosphorus fertilization. Journal of Production Agriculture, 9(1), 114–118. Schnitkey, G. D., & Miranda, M. J. (1993). The Impact of pollution control on livestock‐crop producers. Journal of Agricultural and Resource Economics, 18(1), 25–36. Segerson, K. (1988). Uncertainty and incentives for nonpoint pollution control. Journal of Environmental Economics and Management, 15, 87–88. Sela, S., van Es, H. M., Moebius‐Clune, B. N., Marjerison, R., Moebius‐Clune, D., Schindelbeck, R., … Young, E. (2017). Dynamic model improves agronomic and environmental outcomes for maize nitrogen management over static approach. Journal of Environmental Quality, 46, 311–319. Sharara, M., Sampat, A., Good, L. W., Smith, A. S., Porter, P., Zavala, V. M., Larson, R., & Runge, T. (2017). Spatially explicit methodology for coordinated manure management in shared watersheds. Journal of Environmental Management, 192, 48–56. Sharpley, A., Jarvie, H. P., Buda, A., May, L., Spears, B., & Kleinman, P. (2013). Phosphorus legacy: Overcoming the effects of past management practices to mitigate future water quality impairment. Journal of Environmental Quality, 42, 1308–1326. Shortle, J., & Dunn, J. (1986). The relative efficiency of agricultural source water pollution control policies. American Journal of Agricultural Economics, 68, 668–677. Shortle, J. S., Horan, R. D., & Abler, D. G. (1998). Research issues in nonpoint pollution control. Environmental and Resource Economics, 11(314), 571–585. Sihvonen, M., Hyytiäinen, K., Valkama, E., & Turtola, E. (2018). Phosphorus and nitrogen yield response models for dynamic bio‐economic optimization: An empirical approach. Agronomy, 8, 41. Smith, V. H., & Schindler, D. W. (2009). Eutrophication science: Where do we go from here? Trends in Ecology & Evolution, 24(4), 201–207. Stackpoole, S. M., Stets, E. G., & Sprague, L. A. (2019). Variable impacts of contemporary versus legacy agricultural phosphorus on US river water quality. Proceedings of the National Academy of Sciences of the United States of America, 116(41), 20562–20567. Stauber, M. S., Burt, O. R., & Linse, F. (1975). An economic evaluation of nitrogen fertilization when carry‐over is significant. American Journal of Agricultural Economics, 57, 463–47. Stoecker, A. L., & Onken, A. B. (1989). Optimal fertilizer nitrogen and residual nitrate‐nitrogen levels for irrigated corn and effects of nitrogen limitation: An economic analysis. Journal of Production Agriculture, 2(4), 309–317. Svanbäck, A., McCrackin, M. L., Swaney, D. P., Linefur, H., Gustafsson, B. G., Howarth, R. W., & Humborg, C. (2019). Reducing agricultural nutrient surpluses in a large catchment–Links tolivestock density. Science of the Total Environment, 648, 1549–1559. Sydsæter, K., Hammond, P., Seierstad, A., & Strom, A. (2005). Further Mathematics for Economic Analysis (2nd ed.). Pearson Education Limited. Thomas, A. (2003). A dynamic model of on‐farm integrated nitrogen management. European Journal of Agricultural Economics, 30(4), 439–460. Tilman, D., Cassman, K. G., Matson, P. A., Naylor, R., & Polansky, S. (2002). Agricultural sustainability and intensive production practices. Insight review articles. Nature, 418(8), 671–677. Turtola, E., Salo, T., Miettinen, A., Iho, A., Valkama, E., Rankinen, K., Virkajärvi, P., Tuomisto, J., Sipilä, A. Muurinen, S., Turakainen, M., Lemola, R., Jauhiainen, L., Uusitalo, R., Grönroos, J., Myllys, M., Heikkinen, J., Merilaita, S., Bernal, J. C., … Jaakkola, M. (2017). Hyötyä taseista – Ravinnetaseiden tulkinta ympäristön ja viljelyn hyödyksi. Luonnonvara‐ ja biotalouden tutkimus, 15/2017. Uusitalo, R., Hyväluoma, J., Valkama, E., Ketoja, E., Vaahtoranta, A., Virkajärvi, P., Grönroos, J., Lemola, R., Ylivainio, K., Rasa, K., & Turtola, E. (2016). A simple dynamic model of soil test phosphorus responses to phosphorus balances. Journal of Environmental Quality, 45(3), 977–983. 26 of 40 | Natural Resource Modeling SIHVONEN ET AL. Uusitalo, R., & Jansson, H. (2002). Dissolved reactive phosphorus in runoff assessed by soil extraction with acetate buffer. Agricultural and Food Science, 11, 343–353. Valkama, E., Lemola, R., Känkänen, H., & Turtola, E. (2015). Meta‐analysis of the effects of undersown catch crops on nitrogen leaching loss and grain yields in the Nordic countries. Agriculture, Ecosystems and Environment, 203, 93–101. Vandermeer, J., van Noordwijk, M., Anderson, J., Ong, C., & Perfecto, I. (1998). Global change and multi‐species agroecosystems: Concepts and issues. Agriculture, Ecosystems & Environment, 67(1), 1–22. Vitousek, P. M., Mooney, H. A., Lubchenco, J., & Melillo, J. M. (1997). Human domination of earth's ecosystems. Science, 277, 494–499. Vuorenmaa, J., Rekolainen, S., Kenttämies, K., & Kauppila, P. (2002). Losses of nitrogen and phosphorus from agricultural and forest areas in Finland during the 1980s and 1990s. Environmental Monitoring and Assessment, 76, 213–248. Vuorinen, J., & Mäkitie, O. (1955). The method of soil testing in use in Finland. Agrogeolical Publications, 63, 1–44. Watkins, K. B. (1998). Y‐c., Lu, W‐y., Huang. economic and environmental feasibility of variable rate nitrogen fertilizer application with carryover effects. Journal of Agriculture and Resource Economics, 23(2), 401–426. Weersink, A., Livernois, J., Shogren, J. F., & Shortle, J. S. (1998). Economic instruments and environmental policy in agriculture. Canadian Public Policy—Analyse de Politiques, XXIV(3), 309–327. Weitzman, M. L. (1998). Why the far‐distant future should be discounted at its lowest possible rate. Journal of Environmental Economics and Management, 36, 201–208. Weitzman, M. L. (2001). Gamma discounting. American Economic Review, 91, 260–271. Whiters, P. J. A., Neal, C., Jarvie, H. P., & Doody, D. G. (2014). Agriculture and eutrophication: Where do we go from here? Sustainability, 6, 5853–5875. Xepapadeas, A. P. (1992). Environmental policy design and dynamic nonpoint‐source pollution. Journal of Environmental Economics and Management, 23, 22–39. Yadav, S. N. (1997). Dynamic optimization of nitrogen use when groundwater contamination is internalized at the standard in the long run. American Journal of Agricultural Economics, 79, 931–945. Ylivainio, K., Sarvi, M., Lemola, R., Uusitalo, R., & Turtola, E. (2014). Regional P stocks in soil and in animal manure as compared to P requirement of plants in Finland: Baltic Forum for Innovative Technologies for Sustainable Manure Management. WP4 Standardisation of manure types with focus on phosphorus. MTT REPORT, 124. Retrieved from http://www.mtt.fi/mttraportti/pdf/mttraportti124.pdf SUPPORTING INFORMATION Additional supporting information may be found online in the Supporting Information section. How to cite this article: Sihvonen M, Lintunen J, Valkama E, Hyytiäinen K. Management of legacy nutrient stores through nitrogen and phosphorus fertilization, catch crops, and gypsum treatment. Natural Resource Modeling. 2020;33:e12289. https://doi.org/10.1111/nrm.12289 SIHVONEN ET AL. Natural Resource Modeling | 27 of 40 APPENDIX A: A1: Derivation of the steady‐state shadow prices of the soil p and n stocks If the system is in a steady state, the optimal input decisions remain unchanged from one period to another, and therefore we can drop time indexes from Equations (8)–(13). In such a case we can solve Equation (12) for the steady‐state soil P shadow price λx. First, we drop the time indexes to obtain the form p f β λ f λ ϕ f λ μ e+ [ (1 + ϑ + ϑ ) + ] = + .y x x x f x n f x x P xP (A1) Second, we solve the equation for λx: λ p f μ e λ ϕ f ρ f= − + ( − (ϑ + ϑ )) . x y x P x P ρ n f x ρ x f x 1 1 + 1 1 + (A2) Note that the equation still includes the unknown λn. Similarly, we can solve Equation (13) for the steady‐state soil N shadow price λn. First, we drop the time indexes to obtain the form )(p f β λ ϕ ϕ f ϕ e λ f λ μ e+ 1 + + + + ϑ = + .y n n n f n e nN x f n n N nNN⎡⎣ ⎤⎦ (A3) Second, we solve the equation for λn ( )( )λ p f μ e λ f ρ ϕ ϕ f ϕ e = − + ϑ − + + .n y n N n N ρ x f n ρ n f n e n N 1 1 + 1 1 + N (A4) Now we have a system of equations that we solve for λx and λn. We obtain the expressions for the steady‐state shadow prices: λ p f μ e ϕ f ρ f = − + ( − (ϑ + ϑ )) − andx y x P x P p f μ e ρ ϕ ϕ f ϕ e f x ρ x f x f ϕ f ρ ϕ ϕ f ϕ e − − ( + + ) 1 1 + ϑ − ( + + ) y n N n N n f n eN n N f n f x n f n eN n N ⎡ ⎣⎢ ⎤ ⎦⎥ (A5) ( ) λ p f μ e f ρ ϕ ϕ f ϕ e = − + ϑ − + + − .n y n N n N p f μ e ρ f f n ρ n f n e n N ϕ f f ρ f − − (ϑ + ϑ ) 1 1 + ϑ − (ϑ + ϑ ) y x P x P x f x N f x f n x f x ⎡ ⎣⎢ ⎤ ⎦⎥ (A6) A2: Optimality conditions for the regulated private optimum To obtain the expressions for the annual taxes and the subsidies, we must first derive the first‐ order conditions for the regulated private optimum. Regulated private producer maximizes NPV of the annual net returns from the crop production: 28 of 40 | Natural Resource Modeling SIHVONEN ET AL. { } ( ) ( ) ( ) ( ) ( ) { } β p y p τ P p τ N p σ φ p σ φ x x x P y n n ϕ n N y e x n max − + + + + − + − , Subject to = + ϑ( , , ) = + , , , and is given P N φ φ x n t T t y t P t P t N t N t cc t cc t cc gyp t gyp t gyp t t t t t t t t t t t N , , , , , =0 +1 +1 0 0 t t t cc t gyp t t+1 +1 ⎡⎣ ⎤⎦ ∑ (A7) Lagrangean of the problem is ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) L β p f P N φ x n p τ P p τ N p σ φ p σ φ βη x x P f P N φ x n x βη n ϕ n N f P N φ x n e n γ φ γ φ γ φ γ φ = , , , , − + + + + − + − + + ϑ , , , , , , − + + , , , , , , , − + + + 1 − + 1 − t t y t t t cc t t P t P t N t N t cc t cc t cc gyp t gyp t gyp t x t t t t t t cc t t t t n t t t t t t cc t t t N t t cc t cc t gyp t gyp t cc t cc t gyp t gyp =0 +1 +1 +1 +1 ,0 ,0 ,1 ,1 ⎜ ⎟ ⎛ ⎝ ⎡⎣ ⎤⎦ ⎡⎣⎢ ⎤ ⎦⎥ ⎡ ⎣⎢ ⎤ ⎦⎥ ⎞ ⎠ ∑∞ (A8) Note that here we denote the private shadow prices of the soil P and N stocks by ηx and ηn, respectively. Taking a derivative w.r.t. Nt, dividing the equation with βt and rearranging gives: ( )N p f β η f η ϕ ϕ f p τ: + ϑ + ( + ) = + .t y N t tx f t N t tn N t f t N t N tN, +1 , , +1 , , , (A9) Similarly, taking a derivative w.r.t. Pt, dividing the equation with βt and rearranging gives: P p f βη f βη ϕ f p τ: + (ϑ + ϑ ) + = + .t y P t tx P t f t P t tn f t P t P tP, +1 , , , +1 , , (A10) The derivative of the Lagrangean w.r.t. φtcc, dividing the equation with βt and rearranging gives: ( )φ σ p f β η f η ϕ f γ p γ: + + ϑ + + = + .tcc tcc y φ t tx f t φ t tn f t φ t tcc cc tcc, +1 , , +1 , , ,0 ,1cc cc cc (A11) The derivative of the Lagrangean w.r.t. φtgyp, dividing the equation with βt and rearranging gives: φ σ γ p γ: + = + .tgyp tgyp tgyp gyp tgyp,0 ,1 (A12) The derivative of the Lagrangean w.r.t. xt, dividing the equation with βt and rearranging gives: ( )x p f β η f η ϕ f η: + (1 + ϑ + ϑ ) + = .t y x t tx x t f t x t tn f t x t tx, +1 , , , +1 , , (A13) SIHVONEN ET AL. Natural Resource Modeling | 29 of 40 And finally, the derivative of the Lagrangean w.r.t. nt, dividing the equation with βt and rearranging gives: ( )( )n p f β η ϕ ϕ f ϕ e η f η: + 1 + + + + ϑ = .t y n t tn n t f t n t e t n tN tx f t n t tn, +1 , , , , +1 , ,N (A14) A3: Derivation of the first‐best taxes/subsidies To obtain the expressions for the first‐best taxes, we equalize the conditions for the social and regulated private optimal N and P fertilizer rates, catch crop cultivation and gypsum application. Let us start with the conditions for the P fertilization: p f β λ f λ ϕ f p μ e p f β η f η ϕ f p τ + [ (ϑ + ϑ ) + ] − − = + [ (ϑ + ϑ ) + ] − − y P t t x P t f t P t t n f t P t P P P t P y P t t x P t f t P t t n f t P t P t P , +1 , , , +1 , , , , +1 , , , +1 , , (A15) Cancelling common terms and rearranging gives: ( ) ( )τ μ e β η λ f η λ ϕ f= + − (ϑ + ϑ ) + −tP P P tP tx tx P t f t P t tn tn f t P t, +1 +1 , , , +1 +1 , ,⎡⎣ ⎤⎦ (A16) Let us then focus on the conditions for the N fertilization: p f β λ f λ ϕ ϕ f p μ e p f β η f η ϕ ϕ f p τ + ϑ + ( + ) − − = + ϑ + ( + ) − − y N t t x f t N t t n N t f t N t N N N t N y N t t x f t N t t n N t f t N t N t N , +1 , , +1 , , , , , +1 , , +1 , , , ⎡⎣ ⎤⎦ ⎡⎣ ⎤⎦ (A17) Cancelling common terms and rearranging gives: ( ) ( )τ μ e β η λ f η λ ϕ ϕ f= + − ϑ + − ( + )tN N N tN tx tx f t N t tn tn N t f t N t, +1 +1 , , +1 +1 , , ,⎡⎣ ⎤⎦ (A18) Similarly combining the conditions for the catch crop cultivation and gypsum application gives: ( ) ( )σ μ e β η λ f η λ ϕ f= − − − ϑ + −tcc N φ tN tx tx f t φ t tn tn f t φ t, +1 +1 , , +1 +1 , ,cc cc cc⎡⎣ ⎤⎦ (A19) and σ μ e= − .tgyp P φ tP ,gyp (A20) To obtain the explicit expressions for the taxes and the subsidies, however, we must de- termine the expressions for the differences of the shadow prices η λ−tx tx+1 +1 and η λ−tn tn+1 +1. First, we need the expressions for the forwarded shadow prices. Let us first focus on the conditions for the socially optimal P and N stocks. We solve the conditions for λtx and λtn, and then we forward the obtained expressions by 1 year to obtain 30 of 40 | Natural Resource Modeling SIHVONEN ET AL. λ p f μ e β λ f λ ϕ f= − + (1 + ϑ + ϑ ) +tx y x t P x tP tx x t f t x t tn f t x t+1 , +1 , +1 +2 , +1 , +1 , +1 +2 , +1 , +1⎡⎣ ⎤⎦ (A21) ( )λ p f β λ ϕ ϕ f ϕ e λ f μ e = + 1 + + + + ϑ − t n y n t t n n t f t n t e t n t N t x f t n t N n t N +1 , +1 +2 , +1 , +1 , +1 , +1 , +1 +2 , +1 , +1 , +1 N ⎡⎣ ⎤⎦ (A22) Let us then focus on the conditions for the optimal P and N stocks in the regulated private optimum. We solve the conditions for ηtx and ηtn, and then we forward the obtained expressions by 1 year to obtain η p f β η f η ϕ f= + (1 + ϑ + ϑ ) +tx y x t tx x t f t x t tn f t x t+1 , +1 +2 , +1 , +1 , +1 +2 , +1 , +1⎡⎣ ⎤⎦ (A23) ( )η p f β η ϕ ϕ f ϕ e η f= + 1 + + + + ϑtn y n t tn n t f t n t e t n tN tx f t n t+1 , +1 +2 , +1 , +1 , +1 , +1 , +1 +2 , +1 , +1N⎡⎣ ⎤⎦ (A24) Combining the obtain expressions yields ( ) ( ) η λ μ e β η λ f η λ ϕ f − = + − (1 + ϑ + ϑ ) + − t x t x P x t P t x t x x t f t x t t n t n f t x t +1 +1 , +1 +2 +2 , +1 , +1 , +1 +2 +2 , +1 , +1 ⎡⎣ ⎤⎦ (A25) ( )( ) ( ) η λ μ e β η λ ϕ ϕ f ϕ e η λ f − = + − 1 + + + + − ϑ t n t n N n t N t x t x n t f t n t e t n t N t n t n f t n t +1 +1 , +1 +2 +2 , +1 , +1 , +1 , +1 , +1 +2 +2 , +1 , +1 N ⎡⎣ ⎤⎦ (A26) Thus, we may rewrite the tax‐subsidy‐scheme as follows: τ μ e β f ϕ f= + (ϑ + ϑ )Δ + ΔtP P P tP P t f t P t tx f t P t tn, , , , +1 , , +1⎡⎣ ⎤⎦ (A27) τ μ e β f ϕ ϕ f= + ϑ Δ + ( + )ΔtN N N tN f t N t tx N t f t N t tn, , , +1 , , , +1⎡⎣ ⎤⎦ (A28) σ μ e β f ϕ f= − − ϑ Δ + Δtcc N φ tN f t φ t tx f t φ t tn, , , +1 , , +1cc cc cc⎡⎣ ⎤⎦ (A29) σ μ e= −tgyp P φ tP ,gyp (A30) where μ e β f ϕ fΔ = + (1 + ϑ + ϑ )Δ + Δtx P x tP x t f t x t tx f t x t tn+1 , +1 , +1 , +1 , +1 +2 , +1 , +1 +2⎡⎣ ⎤⎦ (A31) ( )μ e β ϕ ϕ f ϕ e fΔ = + 1 + + + Δ + ϑ Δtn N n tN n t f t n t e t n tN tn f t n t tx+1 , +1 , +1 , +1 , +1 , +1 , +1 +2 , +1 , +1 +2N⎡⎣ ⎤⎦(A32) We rewrite these expressions: μ e β βΔ = + Θ Δ + Θ Δtx P x tP tP tx tPN tn+1 , +1 +1 +2 +1 +2 (A33) SIHVONEN ET AL. Natural Resource Modeling | 31 of 40 μ e β βΔ = + Θ Δ + Θ Δtn N n tN tN tn tNP tx+1 , +1 +1 +2 +1 +2 (A34) where f ϕ f ϕ ϕ f ϕ e f Θ 1 + ϑ + ϑ Θ Θ 1 + + + Θ ϑ t P x t f t x t t PN f t x t t N n t f t n t e t n t N t NP f t n t +1 , +1 , +1 , +1 +1 , +1 , +1 +1 , +1 , +1 , +1 , +1 , +1 +1 , +1 , +1N ≔ ≔ ≔ ≔ We then try to obtain expressions for Δtx+1 and Δtn+1, which does not include terms Δtx+2 and Δtn+2, by using recursion. We focus only on Δtx+1, because Δtn+1 is obtained immediately by the symmetry. Forwarding the expressions (A33) and (A34) we obtain: μ e β βΔ = + Θ Δ + Θ Δtx P x tP tP tx tPN tn+2 , +2 +2 +3 +2 +3 (A35) μ e β βΔ = + Θ Δ + Θ Δtn N n tN tN tn tNP tx+2 , +2 +2 +3 +2 +3 (A36) Plugging (A35) back to (A33) we obtain ( ) ( ) μ e βμ e βμ e β β Δ = + Θ + Θ + Θ Θ + Θ Θ Δ + Θ Θ + Θ Θ Δ t x P x t P P t P x t P N t PN n t N t P t P t PN t NP t x t P t PN t PN t N t n +1 , +1 +1 , +2 +1 , +2 2 +1 +2 +1 +2 +3 2 +1 +2 +1 +2 +3 (A37) Forwarding (A33) by 2 years we obtain μ e β βΔ = + Θ Δ + Θ Δtx P x tP tP tx tPN tn+3 , +3 +3 +4 +3 +4 (A38) Plugging (A38) to (A37) we obtain ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) μ β e β e β e μ β e β e β β β β Δ = + Θ + Θ Θ + Θ Θ + Θ + Θ Θ + Θ Θ + Θ Θ + Θ Θ Θ Δ + Θ Θ + Θ Θ Θ Δ + Θ Θ + Θ Θ Θ Δ + Θ Θ + Θ Θ Θ Δ t x P x t P t P x t P t P t P t PN t NP x t P N t PN n t N t P t PN t PN t N n t N t P t P t PN t NP t P t x t P t P t PN t NP t PN t n t P t PN t PN t N t N t n t P t PN t PN t N t NP t x +1 0 , +1 1 +1 , +2 2 +1 +2 +1 +2 , +3 1 +1 , +2 2 +1 +2 +1 +2 , +3 3 +1 +2 +1 +2 +3 +4 3 +1 +2 +1 +2 +3 +4 3 +1 +2 +1 +2 +3 +4 3 +1 +2 +1 +2 +3 +4 (A39) Now, focusing on the five first terms, we can observe the pattern. Using the sum operators, we obtain following sum of two sums μ β a e μ β b eΔ = +tx P s s s x t s P N s s s n t s N +1 =0 , +1+ =1 , +1+∑ ∑ ∞ ∞ (A40) where a = 10 , a = ΘtP1 +1, and b = ΘtPN1 +1, and the rest are recursively defined by relations a a b= Θ + Θs s t sP s t sNP−1 + −1 + 32 of 40 | Natural Resource Modeling SIHVONEN ET AL. b b a= Θ + Θs s t sN s t sPN−1 + −1 + From symmetry, we can deduce that μ β c e μ β d eΔ = +tn N s s s n t s N P s s s x t s P +1 =0 , +1+ =1 , +1+∑ ∑ ∞ ∞ (A41) where c = 10 , c = ΘtN1 +1, and d = ΘtNP1 +1, and the rest are recursively defined by relations c c d= Θ + Θs s t sN s t sPN−1 + −1 + d d c= Θ + Θs s t sP s t sNP−1 + −1 + Thus, the optimal tax‐subsidy‐scheme is τ μ e β f ϕ f= + (ϑ + ϑ )Δ + ΔtP P P tP P t f t P t tx f t P t tn, , , , +1 , , +1⎡⎣ ⎤⎦ (A42) τ μ e β f ϕ ϕ f= + ϑ Δ + ( + )ΔtN N N tN f t N t tx N t f t N t tn, , , +1 , , , +1⎡⎣ ⎤⎦ (A43) σ μ e β f ϕ f= − − ϑ Δ + Δtcc N φ tN f t φ t tx f t φ t tn, , , +1 , , +1cc cc cc⎡⎣ ⎤⎦ (A44) σ μ e= −tgyp P φ tP ,gyp (A45) where the social costs of phosphorus and nitrogen stocks are defined by the (A40) and (A41), respectively. A4: Static tax‐subsidy‐schemes We can derive the steady‐state tax‐subsidy‐scheme by equalizing the socially optimal steady‐ state optimality conditions: P p f β λ f λ ϕ f p μ e: + [ (ϑ + ϑ ) + ] − − = 0,y P x P f P n f P P P PP (A46) N p f β λ f λ ϕ ϕ f p μ e: + [ ϑ + ( + )] − − = 0,y N x f N n N f N N N NN (A47) φ μ e β λ f λ ϕ f γ p γ p f:− + [ ϑ + ] + − − + = 0,cc N φN x f φ n f φ cc cc cc y φ,0 ,1cc cc cc cc (A48) φ μ e γ p γ:− + − − = 0,gyp P φP gyp gyp gyp,0 ,1gyp (A49) where λ p f μ e ϕ f ρ f = − + ( − (ϑ + ϑ )) − andx y x P x P p f μ e ρ ϕ ϕ f ϕ e f x ρ x f x f ϕ f ρ ϕ ϕ f ϕ e − − ( + + ) 1 1 + ϑ − ( + + ) y n N n N n f n eN n N f n f x n f n eN n N ⎡ ⎣⎢ ⎤ ⎦⎥ (A50) SIHVONEN ET AL. Natural Resource Modeling | 33 of 40 ( ) λ p f μ e f ρ ϕ ϕ f ϕ e = − + ϑ − + + − ,n y n N n N p f μ e ρ f f n ρ n f n e n N ϕ f f ρ f − − (ϑ + ϑ ) 1 1 + ϑ − (ϑ + ϑ ) y x P x P x f x N f x f n x f x ⎡ ⎣⎢ ⎤ ⎦⎥ (A51) and the conditions determining the regulated private optimum: P p f βη f βη ϕ f p τ: + (ϑ + ϑ ) + − − = 0,y P x P f P n f P P P (A52) N p f β η f η ϕ ϕ f p τ: + ( ϑ + ( + )) − − = 0,y N x f N n N f N N N (A53) φ p f β η f η ϕ f γ p σ γ: + ( ϑ + ) + − − − = 0,cc y φ x f φ n f φ cc cc cc cc,0 ,1cc cc cc (A54) φ γ p σ γ: − − − = 0,gyp gyp gyp gyp gyp,0 ,1 (A55) where ( ) ( ) η p f ϕ f ρ f = + ( − (ϑ + ϑ )) − andx y x p f ρ ϕ ϕ f ϕ e f x ρ x f x f ϕ f ρ ϕ ϕ f ϕ e − + + 1 1 + ϑ − + + y n n f n eN n N f n f x n f n eN n N ⎡ ⎣⎢ ⎤ ⎦⎥ (A56) ( ) η p f f ρ ϕ ϕ f ϕ e = + ϑ − + + − .n y n p f ρ f f n ρ n f n e n N ϕ f f ρ f − (ϑ + ϑ ) 1 1 + ϑ − (ϑ + ϑ ) y x x f x N f x f n x f x ⎡ ⎣⎢ ⎤ ⎦⎥ (A57) Next, we simply equalize these optimality conditions and cancel the common terms from both sides of the equations and then we solve the obtained equations for the taxes/subsidies: τ μ e β η λ f β η λ ϕ f= + [ − ](ϑ + ϑ ) + [ − ] ,P P PP x x P f P n n f P (A58) τ μ e β η λ f β η λ ϕ ϕ f= + [ − ]ϑ + [ − ]( + ),N N NN x x f N n n N f N (A59) σ μ e β η λ f β η λ ϕ f= + [ − ]ϑ + [ − ] ,N φN x x f φ n n f φcc cc cc cc (A60) σ μ e= ,P φPgyp gyp (A61) where η λ μ e ϕ f ρ f − = + ( − (ϑ + ϑ ) − andx x P x P μ e ρ ϕ ϕ f ϕ e f x ρ x f x f ϕ f ρ ϕ ϕ f ϕ e − ( + + ) 1 1 + ϑ − ( + + ) N n N n f n eN n N f n f x n f n eN n N⎟ ⎡ ⎣⎢ ⎞ ⎠ ⎤ ⎦⎥ (A62) η λ μ e f ρ ϕ ϕ f ϕ e − = + ϑ − ( + + ) − .n n N n N μ e ρ f f n ρ n f n e n N ϕ f f ρ f − (ϑ + ϑ ) 1 1 + ϑ − (ϑ + ϑ ) P x P x f x N f x f n x f x ⎡ ⎣⎢ ⎤ ⎦⎥ (A63) 34 of 40 | Natural Resource Modeling SIHVONEN ET AL. A5: Functions used in the empirical application Nonlinear yield response function in equation (5) for coarse soils is (Sihvonen et al., 2018): y f P N x x θ x θ Pθ x θ θ N θ Y= ( , , ) = ( (1 + )) ( 1 + exp( ) + 1) + (log( + 1)) 1 + ,t t t t t t t t t N 1,1 −2 2,1 2,2 2 3,1 3,2 2 3,3 0 2 (A64) where yt is the average yield (kg·ha−1·yr−1), f =t average annual yield response on coarse soils, x =t soil P (mg/L1), P =t P fertilizer rate (kg·ha−1·yr−1), N N=t fertilizer rate (kg·ha−1·yr−1), and θ = 0.01331,1 , θ = 0.07972,1 , θ = 0.00952,2 , θ = 0.7273,1 , θ = 0.02693,2 , and θ = 3.653⁎10 ,3,3 −8 YN0 = 2,148 (N control yield (kg/ha), that is, a yield without added N fertilizer). Applied yield response function in Equation (5) for clay soils is (Sihvonen et al., 2018): y f P N x θ x θ P θθ x ω θ N θ Y= ( , , ) = exp( log( + 1) ) ( + 1 + exp( ) + ) + 1 ( + 1) ,t t t t t t t P t N 1,1 2 1/2 2,1 2,2 2,3 ,min 3,1 3,2 0 2 (A65) where θ = 0.3171,1 , θ = 65.801,2 , θ = 0.03722,1 , θ = 0.012,2 , θ = 0.07993,1 , and θ = 0.0377,3,2 Y =N0 2,466, ωP,min = 0.96. The effect of nonlegume catch crop, φtcc on the yield is included in Equation (5) as follows: ( )y f P N φ x φ f P N x= , , , = 1 − 0.03 ( , , ),t t t tcc t tcc t t t⎡⎣ ⎤⎦ (A66) Thus, yield is 3% lower when nonlegume catch crop is cultivated compared to the situation where catch crop is not cultivated (Valkama et al., 2015). Applied soil P transition function in Equation (1) is (Uusitalo et al., 2016): ς ς P ς P x ς xϑ = + + − ,t t t t t1 2 bal, 3 bal, 4 (A67) where ϑ =t P soil transition (mg/L), P =tbal, P balance (kg·ha−1·yr−1), x =t soil P (mg/L), and ς = 0.1945061 , ς = 7.0375 102 −3∗ , ς = 0.41898 103 −3∗ , ς = −0.0289274 , on coarse soils and ς = 0.127381 , ς = 6.7745 102 −3∗ , ς = 0.3225 103 −3∗ , ς = −0.0235914 , on clay soils. Thus, the state Equation (1) takes the following form: x x ς ς P ς P x ς x= + + + −t t t t t t+1 1 2 bal, 3 bal, 4 . For P tbal, a following function is applied (Iho & Laukkanen, 2012a): P P β x β y= − ( log( ) + ) ,bal t t t t, 1 2 (A68) where P =tbal, P balance (kg·ha−1·yr−1), x =t soil P (mg/L), y =t yield (kg·ha−1·yr−1), and β = 0.000841 , β = 0.0001862 . There is no state equation for the soil N for Finnish agricultural soils. Moreover, it is clear that we cannot directly couple existing N transition functions with the functions in our model, because those are estimated and validated for entirely different kind of environmental condi- tions. Therefore, following Hyytiäinen et al. (2011), we used the Coupled heat and mass transfer model for the soil‐plant‐atmosphere system (COUP; Jansson, 2012) to simulate data for the estimation of the soil N transition function. The simulations were done for inorganic nitrogen (N) rates ranging from 0 to 200 kg·ha−1·yr−1 with 25 kg/ha steps. We repeated the simulations for nine initial states: initial total amount of N in the whole soil profile was changed from −50% SIHVONEN ET AL. Natural Resource Modeling | 35 of 40 to 50% of the baseline level with 25% step while holding initial organic C/N ratio fixed. Fur- thermore, initial organic C/N ratios were changed using −50%, −25%, 0%, +50%, and +100% of the baseline while holding initial organic N concentration fixed. The baseline N and C stocks were 1.33e+4 kg/ha and 1.44e+5 kg/ha and the baseline C/N ratio was 10.7. The simulations where done for both coarse and clay soils. The total number of the simulations was 162. CoupModel was parameterized for crop production in southern Finland (Rankinen et al., 2007; Salo et al., 2016). The meteorological data (for the period 1980–2011) were from the Jokioinen Observatory of Finnish Meteorological Institute (60.81°N, 23.50 °E, and altitude 104m). We converted daily simulation output data to annual data. Most of the N in the soil stock is in the form that is not directly usable for the plants. To get some kind of an estimation of the soil N that is potentially usable for the plants, we estimated first a simple model for the both soil textures: y θ n=N t t0 1 0 (A69) where yN t0 is the annual N control yield, that is a yield without added N, in N fertilizer experiments, and nt0 is the annual soil N in the experiments where no N fertilizer was given. The results are shown in Table A1. Thus, approximately 9.3% and 8.2% of the annual soil N stock is in plant available form on clay and coarse soils, respectively. The difference between current and the next period plant available soil N (kg·ha−1·yr−1), n n−t t+1 is the soil N transition, or N carryover, ϕt. In the estimation of the N carryover equation, the explained variable, ϕt is multiplied with the ratio 0.093 and 0.082 on clay and coarse soils, and so is nt in the right‐hand side of the equation. The obtained equation for clay soils was ϕ β n β N β e β y β c β= + + ln ( ) + + ln( ) + deposition ,t t t tN t t t1 2 3 4 5 6 (A70) and for coarse soils ϕ β n β N β e β y β= + + + + deposition ,t t t tN t t1 2 3 4 5 (A71) TABLE A1 Estimation results for the parameters of Equation (A69) Parameter Estimate SE t value Pr(>|t|) Clay soils θ1 0.09339 0.00196 47.6 <2.00E−16 *** Adjusted R2: 0.8935 F‐statistic: 2266 on 1 and 269 DF, p value: <2.2e−16 Coarse soils θ1 0.081997 0.001676 48.92 <2.00E−16 *** Adjusted R2: 0.8986 F‐statistic: 2394 on 1 and 269 DF, p value: <2.2e−16 36 of 40 | Natural Resource Modeling SIHVONEN ET AL. where nt is the annual soil N stock (kg·ha−1·yr−1), Nt is the annual N fertilizer input (kg·ha−1·yr−1), etN is the annual N loss (kg·ha−1·yr−1), yt is the annual yield (kg·ha−1·yr−1). On coarse soils the annual yield was replaced with a logarithm of the yield due to the nonlinear relationship. We also added two factors into the equation that are not considered in the analytical model: soil C stock (kg·ha−1·yr−1), ct and atmospheric N deposit (kg·ha−1·yr−1), depositiont, because those had a significant effect on the N carry over in the data. We ignored the soil C stock in the N transition equation on coarse soils because it did not improve the fit, but instead made the fit much worse. However, we treated both factors as constants fixed to their average level in the optimization. Soil C dynamics are important, but beyond the scope of this study. One more thing to consider in the estimation was the fact that the N loss is also a function of the soil N. Thus, N loss and soil N transition function needed to be estimated simultaneously by using the three state least squares (3SLS) method. We estimated the following N loss function for both soil textures: TABLE A2 3SLS estimation results for the parameters of the soil N carryover equations Parameter Estimate SE t value Pr(>|t|) Clay soils N transition equation: ϕ β n β N β e β y β c β= + + ln ( ) + + ln( ) + depositiont t t tN t t t1 2 3 4 5 6 β1 −3.21E−03 1.40E−04 −22.892 <2.22E−16 *** β2 1.99E−02 1.22E−03 16.32 <2.22E−16 *** β3 −6.17E−01 9.46E−03 −6.523 8.33E−11 *** β4 −6.25E−05 6.35E−05 −9.849 <2.22E−16 *** β5 2.19E−01 3.54E−02 6.192 6.94E−10 *** β6 1.15E−01 7.11E−02 1.621 0.105 Adjusted R2: 0.30 N loss equation: e e=tN a a N a ntotalrunoff + +t t t1 2 3 a1 4.21E−03 6.19E−05 67.99 <2.22E−16 *** a2 5.31E−02 1.48E−04 35.94 <2.22E−16 *** a3 7.08E−04 1.91E−05 37.11 <2.22E−16 *** Adjusted R2: 0.34 Coarse soils N transition equation: ϕ β n β N β e β y β= + + + + depositiont t t tN t t1 2 3 4 5 β1 −3.84E−03 2.08E−04 −18.47 <2.22E−16 *** β2 2.43E−02 1.64E−03 14.76 <2.22E−16 *** β3 −1.86E−01 1.80E−03 −10.34 <2.22E−16 *** β4 −6.77E−04 7.42E−05 −9.12 <2.22E−16 *** β5 3.68E−01 5.13E−02 7.185 8.90E−13 *** Adjusted R2: 0.22 N loss equation: e e=tN a a a N a n+ totalrunoff + +t t t1 2 3 4 a1 1.73 5.18E−02 33.35 <2.22E−16 *** a2 2.58E−03 7.42E−05 34.80 <2.22E−16 *** a3 4.64E−03 1.40E−04 33.22 <2.22E−16 *** a4 5.67E−04 3.36E−05 16.91 <2.22E−16 *** Adjusted R2: 0.46 SIHVONEN ET AL. Natural Resource Modeling | 37 of 40 e e=tN a a N a ntotalrunoff + +t t t1 2 3 (A72) However, the equation for the N loss on coarse soils also included a constant term. Estimation results are shown in Table A2. The effect of the current plant available soil N on the N carryover is negative. Overall goodness was evaluated with McEloroy's R2, which was 0.36 for coarse soils models and 0.32 for clay soils models. The adjusted R2 values of the individual models are relatively low, particularly for the N transitions function on coarse soils. If the models would have been estimated sepa- rately by OLS, the R2 values would have been considerably higher, particularly for the models on clay soils (Table A3). Also, the differences in the parameter estimates and TABLE A3 OLS estimation results for the parameters of the soil N carryover equations Parameter Estimate SE t value Pr(>|t|) Clay soils N transition equation: ϕ β n β N β e β y β c β= + + ln ( ) + + ln( ) + depositiont t t tN t t t1 2 3 4 5 6 β1 −3.13E−03 1.40E−04 −22.32 <2.22E−16 *** β2 2.05E−02 1.22E−03 16.74 <2.22E−16 *** β3 −7.17E−01 9.47E−03 −7.57 5.17E−14 *** β4 −6.34E−05 6.35E−05 −9.98 <2.22E−16 *** β5 2.21E−01 3.54E−02 6.233 5.40E−10 *** β6 1.51E−01 7.11E−02 2.121 0.034 * Adjusted R2: 0.65 N loss equation: e e=tN a a N a ntotalrunoff + +t t t1 2 3 a1 4.20E−03 6.19E−05 67.97 <2.22E−16 *** a2 5.31E−02 1.48E−04 35.94 <2.22E−16 *** a3 7.08E−04 1.91E−05 37.12 <2.22E−16 *** Adjusted R2: 0.97 Coarse soils N transition equation: ϕ β n β N β e β y β= + + + + depositiont t t tN t t1 2 3 4 5 β1 −3.47E−03 1.95E−04 −17.76 <2.22E−16 *** β2 2.61E‐02 1.65E−03 15.84 <2.22E−16 *** β3 −2.57E−02 1.80E−03 −14.36 <2.22E−16 *** β4 −6.89E−04 7.44E−05 −9.25 <2.22E−16 *** β5 3.26E−01 4.55E−02 7.17 1E−12 *** Adjusted R2: 0.62 N loss equation: e e=tN a a a N a n+ totalrunoff + +t t t1 2 3 4 a1 1.76 5.21E−02 33.83 <2.22E−16 *** a2 2.55E−03 7.45E−05 34.25 <2.22E−16 *** a3 4.61E−03 1.40E−04 32.95 <2.22E−16 *** a4 5.53E−04 3.36E−05 16.44 <2.22E−16 *** Adjusted R2: 0.46 38 of 40 | Natural Resource Modeling SIHVONEN ET AL. the standard errors between the models estimated by 3SLS and OLS are minor, or almost non‐ existent. The biggest difference is the reduction in the estimate for the effect of the N loss on N transition on coarse soils. Hence, we conclude that the models are rather robust. The N loss functions were scaled to match the recent estimates of N‐loss from coarse and clay soil (Turtola et al., 2017). The applied yield response functions for N response, that is, the third elements in Equations (A64) and (A65) are functions of annual N fertilization rates. Thus, we had to add the effect of the plant available soil N to these functions. We assumed that the annual total N input is equal to the plant N uptake. We obtained the estimates for the N uptake equation again by estimation. Thus, we regressed the N yield uptake on the N fertilization and soil N on both soil textures: N a N a nuptake = +t t t1 2 (A73) Estimation results suggest that 1% increase in N fertilization increases N uptake by 71% and 1% increase in plant available soil N stock increases N uptake by 2.4% on clay soils (Table A4). On coarse soils the corresponding increases are 63% for the N fertilization and 2% for the soil N. Variation in the data is explained by the simple linear models. For DRP‐loss in equation (4), we applied a function from Uusitalo and Jansson (2002) with runoff parameters are from Ekholm et al. (2005): e ζ θ x θ P θ= [run off ( ( + )) − ]/100,tDRP DRP t t t1 2 3 (A74) where e =tDRP DRP‐loss (kg·ha−1·yr−1), x =t soil P (mg/L−1), P =t P fertilizer rate (kg·ha−1·yr−1), runoff = 270t on clay soils, and runoff = 230t on coarse soils, θ = 0.0211 , θ = 0.012 , θ = 0.0153 . For PP‐loss in Equation (3), we again applied a function from Uusitalo and Jansson (2002) with runoff parameters from Ekholm et al. (2005): e ζ θ x θ P θ= erosion { ln[ + ] − }10 ,tPP PP t t t1 2 3 −6 (A75) TABLE A4 Estimation results for the parameters of the N yield uptake function Parameter Estimate SE t value Pr(>|t|) Clay soils a1 0.713 0.0102 69.98 <2.00E−16 *** a2 0.024 0.0009 27.13 <2.00E−16 *** Adjusted R2: 0.91 F‐statistic: 1.216E+4 on 2 and 2428 DF, p value: <2.2e−16 Coarse soils a1 0.63 0.0095 65.94 <2.00E−16 *** a2 0.02 0.0010 20.99 <2.00E−16 *** Adjusted R2: 0.8996 F‐statistic: 1.089E+4 on 2 and 2428 DF, p value: <2.2e−16 SIHVONEN ET AL. Natural Resource Modeling | 39 of 40 where e =tPP PP‐loss (kg·ha−1·yr−1), x =t soil P (mg/L−1), P =t P fertilizer rate (kg·ha−1·yr−1), and erosion = 1200t on clay soils, erosion = 500t on coarse soils, θ = 2501 , θ = 0.012 , θ = 1503 . The effect of gypsum on the P loss is included in Equation (3) as follows: ( ) ( )e δ φ e P x δ φ e P x= 1 − 1 − ( , ) + 1 − 1 − ( , )tP PPgyp tgyp PP t t DPRgyp tgyp DPR t t⎡⎣ ⎤⎦ ⎡⎣ ⎤⎦ (A76) where δ [0,1]PPgyp ∈ is PP‐loss reduction effect of the catch crop cultivation, and δ [0,1]DPRgyp ∈ is DRP‐loss reduction effect of the catch crop cultivation. We use the estimate of 0.36 for δPPgyp and 0.67 for δDPRgyp (Iho & Laukkanen, 2012a). Thus, the applied P loss function is ( ) ( ) e δ φ ζ θ x θ P θ δ φ ζ θ x θ P θ = 1 − 1 − [runoff ( ( + )) − ] 100 + 1 − 1 − erosion { ln[ + ] − } 10 . t P PP gyp t gyp DRP t t t DPR gyp t gyp PP t t t 1 2 3 1 2 3 6 ⎡⎣ ⎤⎦ ⎡⎣ ⎤⎦ (A77) 40 of 40 | Natural Resource Modeling SIHVONEN ET AL.