Relocating croplands could drastically reduce the environmental impacts of global food production

We use the notation in Table 1.

Current crop production and areas, P
_i(x), H
_i(x)

We used 5-arc-minute maps of the fresh-weight production P_i(x) (Mg year⁻¹) and cropping area H_i(x) (ha) of 25 major crops (Table 2) in the year 2010³⁷. These represent the most recent spatially explicit and crop-specific global data⁷⁵. Separate maps were available for irrigated and rainfed croplands, allowing us to estimate the worldwide proportion of irrigated areas as 21% of all croplands.

Agro-ecologically attainable yields (widehatY_i(x))

We used 5-arc-minute maps of the agro-ecologically attainable dry-weight yield (Mg ha ⁻¹ year⁻¹) 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 conditions⁴². 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 CO₂ fertilisation effect, for four future (2071–2100 period) climate scenarios corresponding to representative concentration pathways (RCPs) 2.6, 4.5, 6.0, and 8.5⁷⁶ simulated by the HadGEM2-ES model⁷⁷. 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 application⁴². 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, (widehatY_i(x)), using crop-specific conversion factors^42,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 spells⁴². 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 production⁷⁹, 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 change^80,81.

Carbon impact C
_i(x)

Following an earlier approach⁸, the carbon impact of crop production, C_i(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 map⁸ (Supplementary Fig. 1a), and carbon stored in the crops, for which we used available estimates^8,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 approach⁸, and supported by empirical meta-analyses^{82,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 stocks⁷ (Supplementary Fig. 1b). Thus, the total local carbon impact (Mg C ha⁻¹) of the production of crop i in the grid cell x was estimated as

where (C_{rmpotential,{{rmvegetation}}}(x)) and (C_{rmpotential,rmSOC}(x)) denote the potential natural carbon stocks in the vegetation and the soil in x, respectively, and (C_{rmcrop}(i)) denotes the carbon stocks of crop i (all in Mg C ha⁻¹). 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 paddies⁸⁷. 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 emissions⁸⁸ and poorly correlated with the distance travelled by agricultural products⁸⁹.

Biodiversity impact B
_i(x)

Analogous to our approach for carbon, we estimated the biodiversity impact of crop production, B_i(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 general^{39,90,91,92,93} and the protection of endemic species in particular³⁹. 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 cover⁹⁴. The methodology used to generate these data³⁸ combines species-specific extents of occurrence (spatial envelopes of species’ outermost geographic limits⁴⁰) and habitat preferences (lists of land cover categories in which species can live⁹⁵), both available for all mammals, birds, and amphibians^96,97, with a global map of potential natural biomes⁴⁴ 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 risk^40,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 bodies⁴², land unsuitable due to soil and terrain constraints⁴², built-up land (urban areas, infrastructure, roads)¹, pasture lands¹, crops not considered in our analysis³⁷, or protected areas⁴² (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 relocation

We 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 (widehatH_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

the current global production of crop i, any solution (widehatH_i(x)) must satisfy the equality constraints

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

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

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

where the weighting parameter α ranges between 0 and 1.

Considering all crops across all grid cells, we denote by

the global carbon impact associated with the current distribution of croplands, and by

the global carbon impact associated with the optimal distribution (\widehatH_i(x)_i,x(=\widehatH_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

Using analogous notation, the relative change between the current and the optimal global biodiversity impact across all crops and grid cells is given by

The dark blue line in Fig. 3a visualises (widehatc(alpha )) and (widehatb(alpha )) for the full range of carbon-biodiversity weightings (alpha in [0,1]), each of which corresponds to a specific optimal distribution (\widehatH_i(x)_i,x) of croplands. We defined an optimal weighting (alpha _rmopt), 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

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 cropping⁹⁸, 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 relocation

In 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 relocation

When (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

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,

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 x₂ 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

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 barP_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,

Solving the optimisation problem

All datasets needed in the optimisation (i.e., (A(x)), (P_i(x)), (H_i(x)), (C_i(x)), (B_i(x)), (widehatY_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 × 10⁶-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 maintained^31,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 Toolbox¹⁰⁰ for a termination tolerance on the dual feasibility of 10⁻⁷ and a feasibility tolerance for constraints of 10⁻⁴.

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 trajectories

Our 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 stocks⁸². 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).

Applying the values in Table 3 to a global map of potential natural biomes⁴⁴ provides a map of carbon recovery times. We assumed a square root-shaped carbon recovery trajectory across these regeneration periods¹⁰¹; similar trajectories, sometimes modelled by faster-converging exponential functions, have been identified in other studies^{25,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

where, using the same notation as further above

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.