Relocating croplands could drastically reduce the environmental impacts of global food production
We use the notation in Table 1.Table 1 Notation used in the description of the optimisation framework.Full size tableCurrent crop production and areas, P
i(x), H
i(x)We used 5-arc-minute maps of the fresh-weight production Pi(x) (Mg year−1) and cropping area Hi(x) (ha) of 25 major crops (Table 2) in the year 201037. These represent the most recent spatially explicit and crop-specific global data75. Separate maps were available for irrigated and rainfed croplands, allowing us to estimate the worldwide proportion of irrigated areas as 21% of all croplands.Table 2 Crops included in the analysis.Full size tableAgro-ecologically attainable yields ({widehat{Y}}_{i}(x))
We used 5-arc-minute maps of the agro-ecologically attainable dry-weight yield (Mg ha −1 year−1) of the same 25 crops on worldwide potential growing areas (Supplementary Movie 3) from the GAEZ v4 model, which incorporates thermal, moisture, agro-climatic, soil, and terrain conditions42. These yield estimates were derived based on the assumption of rainfed water supply (i.e., without additional irrigation) and are available for current climatic conditions and, assuming a CO2 fertilisation effect, for four future (2071–2100 period) climate scenarios corresponding to representative concentration pathways (RCPs) 2.6, 4.5, 6.0, and 8.576 simulated by the HadGEM2-ES model77. Potential rainfed yield estimates for current climatic conditions were available for a low- and a high-input crop management level, representing, respectively, subsistence-based organic farming systems and advanced, fully mechanised production using high-yielding crop varieties and optimum fertiliser and pesticide application42. We additionally considered potential yields representing a medium-input management scenario, given by the mean of the relevant low- and high-input yields. Future potential yields were available only for the high-input management level. Thus, we considered a total of 175 (=25 × 3 present + 25 × 4 future) potential yield maps. Potential dry-weight yields were converted to fresh-weight yields, ({widehat{Y}}_{i}(x)), using crop-specific conversion factors42,78.Both current and future potential rainfed yields from GAEZ v4 were simulated based on daily weather data, and therefore account for short-term events such as frost days, heat waves, and wet and dry spells42. However, the estimates represent averages of annual yields across 30-year periods; thus, whilst the need for irrigation on cropping areas identified in our approach during particularly dry years may in principle be obviated by suitable storage of crop production79, in practice, ad hoc irrigation may be an economically desirable measure to maintain productivity during times of drought, which are projected to increase in different geographic regions due to climate change80,81.Carbon impact C
i(x)Following an earlier approach8, the carbon impact of crop production, Ci(x), in a 5-arc-minute grid cell was estimated as the difference between the potential natural carbon stocks and the cropland-specific carbon stocks, each given by the sum of the relevant vegetation- and soil-specific carbon. The change in vegetation carbon stocks resulting from land conversion is given by the difference between carbon stored in the potential natural vegetation, available as a 5-arc-minute global map8 (Supplementary Fig. 1a), and carbon stored in the crops, for which we used available estimates8,78. Regarding soil, spatially explicit global estimates of soil organic carbon (SOC) changes from land cover change are not available. We therefore chose a simple approach, consistent with estimates across large spatial scales, rather than a complex spatially explicit model for which, given the limited empirical data, robust predictions across and beyond currently cultivated areas would be difficult to achieve. Following an earlier approach8, and supported by empirical meta-analyses82,83,84,85,86, we assumed that the conversion of natural habitat to cropland results in a 25% reduction of the potential natural SOC. For the latter, we used a 5-arc-minute global map of pre-agricultural SOC stocks7 (Supplementary Fig. 1b). Thus, the total local carbon impact (Mg C ha−1) of the production of crop i in the grid cell x was estimated as$${C}_{i}(x)={{C}}_{{{{{{rm{potential}}}}}},{{{{{rm{vegetation}}}}}}}(x)+0.25cdot {C}_{{{{{{rm{potential}}}}}},{{{{{rm{SOC}}}}}}}(x)-{C}_{{{{{{rm{crop}}}}}}}(i)$$
(1)
where ({{C}}_{{{{{{rm{potential}}}}}},{{{{{rm{vegetation}}}}}}}(x)) and ({C}_{{{{{{rm{potential}}}}}},{{{{{rm{SOC}}}}}}}(x)) denote the potential natural carbon stocks in the vegetation and the soil in x, respectively, and ({C}_{{{{{{rm{crop}}}}}}}(i)) denotes the carbon stocks of crop i (all in Mg C ha−1). By design, the approach allows us to estimate the carbon impact of the conversion of natural habitat to cropland regardless of whether an area is currently cultivated or not.In our analysis, we did not consider greenhouse gas emissions from sources other than from land use change, including nitrous emissions from fertilised soils and methane emissions from rice paddies87. In contrast to the one-off land use change emissions considered here, those are ongoing emissions that incur continually in the production process. We would assume that the magnitude of these emissions in a scenario of redistribution of agricultural areas, in which the total production of each crop remains constant, is roughly similar to that associated with the current distribution of areas. We also did not consider emissions associated with transport; however, these have been shown to be small compared to other food chain emissions88 and poorly correlated with the distance travelled by agricultural products89.Biodiversity impact B
i(x)Analogous to our approach for carbon, we estimated the biodiversity impact of crop production, Bi(x), in a 5-arc-minute grid cell as the difference between the local biodiversity associated with the natural habitat and that associated with cropland. For our main analysis, we quantified local biodiversity in terms of range rarity (given by the sum of inverse species range sizes; see below) of mammals, birds, and amphibians. Range rarity has been advocated as a biodiversity measure particularly relevant to conservation planning in general39,90,91,92,93 and the protection of endemic species in particular39. In a supplementary analysis, we additionally considered biodiversity in terms of species richness.We used 5-arc-minute global maps of the range rarity and species richness of mammals, birds, and amphibians under potential natural vegetation (Supplementary Fig. 1c, d) and under cropland land cover94. The methodology used to generate these data38 combines species-specific extents of occurrence (spatial envelopes of species’ outermost geographic limits40) and habitat preferences (lists of land cover categories in which species can live95), both available for all mammals, birds, and amphibians96,97, with a global map of potential natural biomes44 in order to estimate which species would be present in a grid cell for natural habitat conditions. Incorporating information on species’ ability to live in croplands, included in the habitat preferences, allows for determining the species that would, and those that would not, tolerate a local conversion of natural habitat to cropland. The species richness impact of crop production in a grid cell is then obtained as the number of species estimated to be locally lost when natural habitat is converted to cropland. Instead of weighing all species equally, the range rarity impact in a grid cell is calculated as the sum of the inverse potential natural range sizes of the species locally lost when natural habitat is converted; thus, increased weight is attributed to range-restricted species, which tend to be at higher extinction risk40,41.As in the case of carbon, the approach allows us to estimate the biodiversity impact of crop production in both currently cultivated and uncultivated areas.Land potentially available for agriculture, V(x)We defined the area V(x) (ha) potentially available for crop production in a given grid cell x, as the area not currently covered by water bodies42, land unsuitable due to soil and terrain constraints42, built-up land (urban areas, infrastructure, roads)1, pasture lands1, crops not considered in our analysis37, or protected areas42 (Supplementary Fig. 1e). In the scenario of a partial relocation of crop production, in which a proportion of existing croplands is not moved, the relevant retained areas are additionally subtracted from the potentially available area, as described further below.Optimal transnational relocationWe first consider the scenario in which all current croplands are relocated across national borders based on current climate (Fig. 3a, dark blue line). For each crop i and each grid cell x, we determined the local (i.e., grid-cell-specific) area ({widehat{H}}_{i}(x)) (ha) on which crop i is grown in cell x so that the total production of each crop i equals the current production and the environmental impact is minimal. Denoting by$${bar{P}}_{i}={sum }_{x}{P}_{i}(x)$$
(2)
the current global production of crop i, any solution ({widehat{H}}_{i}(x)) must satisfy the equality constraints$${sum }_{x}{widehat{H}}_{i}(x)cdot {widehat{Y}}_{i}(x)={bar{P}}_{{{{{{rm{i}}}}}}},{{{{{rm{for}}}}}}quad{{{{{rm{each}}}}}},{{{{{rm{crop}}}}}},i$$
(3)
requiring the total production of each individual crop after relocation to be equal to the current one. A solution must also satisfy the inequality constraints$${sum }_{i}{widehat{H}}_{i}(x)le V(x),{{{{{rm{for}}}}}}quad{{{{{rm{each}}}}}},{{{{{rm{grid}}}}}},{{{{{rm{cell}}}}}},x,,$$
(4)
ensuring that the local sum of cropping areas is not larger than the locally available area V(x) (see above). Given these constraints, we can identify the global configuration of croplands that minimises the associated total carbon or biodiversity impact by minimising the objective function$${sum }_{x}{widehat{H}}_{i}(x)cdot {C}_{i}(x)to ,{{min }}quad{{{{{rm{or}}}}}}quad{sum }_{x}{widehat{H}}_{i}(x)cdot {B}_{i}(x)to ,{{min }}$$
(5)
respectively. More generally, we can minimise a combined carbon and biodiversity impact measure, and examine potential trade-offs between minimising each of the two impacts, by considering the weighted objective function$${sum }_{x}{widehat{H}}_{i}(x)cdot (alpha cdot {C}_{i}(x)+(1-alpha )cdot {B}_{i}(x))to ,{{min }}$$
(6)
where the weighting parameter α ranges between 0 and 1.Considering all crops across all grid cells, we denote by$$bar{C}={sum }_{i}{sum }_{x}{H}_{i}(x)cdot {C}_{i}(x)$$
(7)
the global carbon impact associated with the current distribution of croplands, and by$$hat{C}(alpha )={sum }_{i}{sum }_{x}{hat{H}}_{i}(x)cdot {C}_{i}(x)$$
(8)
the global carbon impact associated with the optimal distribution ({{{widehat{H}}_{i}(x)}}_{i,x}(={{{widehat{H}}_{i}^{alpha }(x)}}_{i,x})) of croplands for some carbon-biodiversity weighting (alpha in [0,1]). The relative change between the current and the optimal carbon impact is then given by$$hat{c}(alpha )=100 % cdot frac{hat{C}(alpha )-bar{C}}{bar{C}}$$
(9)
Using analogous notation, the relative change between the current and the optimal global biodiversity impact across all crops and grid cells is given by$$widehat{b}(alpha )=100 % cdot frac{widehat{B}(alpha )-bar{B}}{bar{B}}$$
(10)
The dark blue line in Fig. 3a visualises (widehat{c}(alpha )) and (widehat{b}(alpha )) for the full range of carbon-biodiversity weightings (alpha in [0,1]), each of which corresponds to a specific optimal distribution ({{{widehat{H}}_{i}(x)}}_{i,x}) of croplands. We defined an optimal weighting ({alpha }_{{{{{{rm{opt}}}}}}}), meant to represent a scenario in which the trade-off between minimising the total carbon impact and minimising the total biodiversity impact is as small as possible. Such a weighting is necessarily subjective; here, we defined it as$${alpha }_{{{{{{rm{opt}}}}}}}={{arg }},{{{min }}}_{alpha in [0,1]}left|begin{array}{ll}frac{frac{partial {hat{c}}(alpha)} {partial {hat{b}}(alpha)}}{hat{c}(alpha)} cdot frac{frac{partial {hat{b}}(alpha)} {partial {hat{c}}(alpha)}}{hat{b}(alpha)}end{array}right|$$
(11)
Each of the two factors on the right-hand side represents the relative rate of change in the reduction of one impact type with respect to the change in the reduction of the other one as α varies. Thus, αopt represents the weighting at which neither impact type can be further reduced by varying α without increasing the relative impact of the other by at least the same amount. Scenarios based on this optimal weighting are shown in Figs. 1, 2, and Supplementary Figs. 3–6, and are represented by the black markers in Fig. 3.Our approach does not account for multiple cropping; i.e., part of a grid cell is not allocated to more than one crop, and the assumed annual yield is based on a single harvest. Allowing for multiple crops to be successively planted in the same location during a growing period would increase the dimensionality of the optimisation problem substantially. However, given that only 5% of current global rainfed areas are under multiple cropping98, this is likely not a strong limitation of our rainfed-based analysis. As a result of this approach, our results may even slightly underestimate local crop production potential and therefore global impact reduction potentials.Optimal national relocationIn the case of areas being relocated within national borders, the mathematical framework is identical with the exception that the sum over relevant grid cells x in Eqs. (2) and (4) is taken over the cells that define the given country of interest, instead of the whole world. In this way, the total production of each crop within each country for optimally distributed croplands is the same as for current areas. The optimisation problem is then solved independently for each country.Optimal partial relocationWhen (either for national or transnational relocation) only a certain proportion (lambda in [0,1]) of the production of each crop (of a country or the world) is being relocated rather than the total production, Eq. (3) changes to$$mathop{sum}limits_{x}{widehat{H}}_{i}(x)cdot {widehat{Y}}_{i}(x)=lambda cdot {bar{P}}_{i},{{{{{rm{for}}}}}},{{{{{rm{each}}}}}},{{{{{rm{crop}}}}}},i,.$$
(12)
In addition, the area potentially available for new croplands, V(x), (see above) is reduced by the area that remains occupied by current croplands accounting for the proportion ((1-lambda )) of production that is not being relocated. We denote by ({H}_{i}^{lambda }(x)) the area that continues to be used for the production of crop i in grid cell x in the scenario where the proportion λ of the production is being optimally redistributed. In particular, ({H}_{i}^{0}(x)={H}_{i}(x)) and ({H}_{i}^{1}(x)=0) for all i and x. For a given carbon-biodiversity weighting (alpha in [0,1]) in Eq. (6), ({H}_{i}^{lambda }(x)) is calculated as follows. First, all grid cells in which crop i is currently grown are ordered according to their agro-environmental efficiency, i.e., the grid-cell-specific ratio between the environmental impact attributed to the production of the crop and the local production,$${E}_{i}^{alpha }(x)=frac{{H}_{i}(x)cdot (alpha cdot {C}_{i}(x)+(1-alpha )cdot {B}_{i}(x))}{{P}_{i}(x)}.$$
(13)
Let ({x}_{1}(={x}_{1}(i,alpha ))) denote the index of the grid cell in which crop i is currently grow for which ({E}_{i}^{alpha }) is smallest among all grid cells in which the crop is grown. Then let x2 be the index for which ({E}_{i}^{alpha }) is second smallest (or equal to the smallest), and so on. Thus, the vector (({x}_{1},{x}_{2},{x}_{3},ldots )) contains all indices of grid cells where crop i is currently grown in descending order of agro-environmental efficiency. The area ({H}_{i}^{lambda }({x}_{n})) retained in some grid cell ({x}_{n}) is then given by$${H}_{i}^{lambda }({x}_{n})=left{begin{array}{ll}{H}_{i}({x}_{n}) & {{{{{rm{if}}}}}};mathop{sum }limits_{m=1}^{n}{P}_{i}({x}_{m})le (1-lambda )cdot {bar{P}}_{i}\ 0, & hskip-7.5pc{{{{{rm{else}}}}}}end{array}right.$$
(14)
Thus, cropping areas in a grid cell ({x}_{n}) are retained if they are amongst the most agro-environmentally efficient ones of crop i on which the combined production does not exceed ((1-lambda )cdot {bar{P}}_{i}) (which is not being relocated). Growing areas in the remaining, less agro-environmental efficient grid cells are abandoned and become potentially available for other relocated crops. Note that ({H}_{i}^{lambda }) depends on the weighting α of carbon against biodiversity impacts. Finally, instead of Eq. (4), we have, in the case of the partial relocation of the proportion λ of the total production,$$mathop{sum}limits_{i}{widehat{H}}_{i}(x)le V(x)-{H}_{i}^{lambda }(x)quad{{{{{rm{for}}}}}},{{{{{rm{each}}}}}},{{{{{rm{grid}}}}}},{{{{{rm{cell}}}}}},x,.$$
(15)
Solving the optimisation problemAll datasets needed in the optimisation (i.e., (A(x)), ({P}_{i}(x)), ({H}_{i}(x)), ({C}_{i}(x)), ({B}_{i}(x)), ({widehat{Y}}_{i}(x)), (V(x))) are available at a 5 arc-minute (0.083°) resolution; however, computational constraints required us to upscale these to a 20-arc-minute grid (0.33°) spatial grid. At this resolution, Eq. (6) defines a 1.12 × 106-dimensional linear optimisation problem in the scenario of across-border relocation. The high dimensionality of the problem is in part due to the requirement in Eq. (3) that the individual production level of each crop is maintained. Requiring instead that, for example, only the total caloric production is maintained31,99 reduces Eq. (6) to a 1-dimensional problem. However, in such a scenario, the production of individual crops, and therefore of macro- and micronutrients, would generally be very different from current levels, implicitly assuming potentially drastic dietary shifts that may not be nutritionally or culturally realistic.The optimisation problem in Eq. (6) was solved using the dual-simplex algorithm in the function linprog of the Matlab R2021b Optimization Toolbox100 for a termination tolerance on the dual feasibility of 10−7 and a feasibility tolerance for constraints of 10−4.In the case of a transnational relocation of crop production, the algorithm always converged to the optimal solution, i.e., for all crop management levels, climate scenarios, and proportions of production that were being relocated. For the relocation within national borders, this was not always the case. This is because some countries produce small quantities of crops which, according to the GAEZ v4 potential yield estimates, could not be grown in the relevant quantities anywhere in the country under natural climatic conditions and for rainfed water supply; these crops likely require greenhouse cultivation or irrigation can therefore not be successfully relocated within our framework. Across all countries, this was the case for production occurring on 0.6% of all croplands. When this was the case for a certain country and crop, we excluded the crop from the optimisation routine, and a country’s total carbon and biodiversity impacts were calculated as the sum of the impacts of optimally relocated crops plus the current impacts of non-relocatable crops.This issue is linked to why determining the optimal distribution of croplands within national borders is not a well-defined problem for future climatic conditions. Under current climatic conditions, if a crop cannot be relocated within our framework, then its current distribution offers a fall-back solution that provides the current production level and allows us to quantify environmental impacts. Different climatic conditions in the future mean that the production of a crop across current growing locations will not be the same as it is today, and therefore the fall-back solution available for the present is no longer available, so that a consistent quantification of the environmental impacts of a non-relocatable crop is not possible.Carbon and biodiversity recovery trajectoriesOur analysis in Supplementary Fig. 6 requires spatially explicit estimates of the carbon recovery trajectory on abandoned croplands. Whilst carbon and biodiversity regeneration have been shown to follow certain general patterns, recovery is context-specific (Supplementary Note 1) in that, depending on local conditions, the regeneration in a specific location can take place at slower or faster speeds than would typically be the case in the broader ecoregion. Here, we assumed that these caveats can be accommodated by using conservative estimates of recovery times and by assuming that local factors will average out at the spatial resolution of our analysis. The carbon recovery times assumed here are based on ecosystem-specific estimates of the time required for abandoned agricultural areas to retain pre-disturbance carbon stocks82. Aiming for a conservative approach, we assumed carbon recovery times equal to at least three times these estimates, rounded up to the nearest quarter century (Table 3). Independent empirical estimates from specific sites and from meta-analyses are well within these time scales (Supplementary Note 1).Table 3 Assumed times required for carbon stocks on abandoned cropland to reach pre-disturbance levels.Full size tableApplying the values in Table 3 to a global map of potential natural biomes44 provides a map of carbon recovery times. We assumed a square root-shaped carbon recovery trajectory across these regeneration periods101; similar trajectories, sometimes modelled by faster-converging exponential functions, have been identified in other studies25,27,30,102,103,104,105. Thus, the carbon stocks in an area of a grid cell x previously used to grow crop i were assumed to regenerate according to the function$$left{begin{array}{ll}{{C}}_{{{{{{rm{agricultural}}}}}}}(x)+sqrt{frac{t}{{{T}}_{{{{{{rm{carbon}}}}}}}(x)}}cdot ({{C}}_{{{{{{rm{potential}}}}}}}(x)-{{C}}_{{{{{{rm{agricultural}}}}}}}(x)) & {{{{{rm{if}}}}}},t ; < ; {{T}}_{{{{{{rm{carbon}}}}}}}\ hskip14.7pc{{C}}_{{{{{{rm{potential}}}}}}}(x) & {{{{{rm{if}}}}}},tge {{T}}_{{{{{{rm{carbon}}}}}}}end{array}right.$$
(16)
where, using the same notation as further above$${{C}}_{{{{{{rm{potential}}}}}}}(x) ={{C}}_{{{{{{rm{potential}}}}}},{{{{{rm{vegetation}}}}}}}(x)+{{C}}_{{{{{{rm{potential}}}}}},{{{{{rm{SOC}}}}}}}(x)\ {{C}}_{{{{{{rm{agricultural}}}}}}}(x) ={{C}}_{i}(x)+0.75cdot {{C}}_{{{{{{rm{potential}}}}}},{{{{{rm{SOC}}}}}}}(x)$$
(17)
Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More
