Land Suitability for Cocoa Cultivation in Peru: AHP and MaxEnt Modeling in a GIS Environment

Abstract

Peru is one of the world’s leading exporters of cocoa beans, which directly impacts the household economy of millions of small farmers. Currently, the expansion and modernization of the cocoa-growing area require the zoning of the territory with suitable biophysical and infrastructural conditions to facilitate optimizing productivity factors. Therefore, we analyzed land suitability for cocoa (Theobroma cacao L.) production on the Peruvian mainland as a support measure for sustainable agriculture. To this end, the climatological, edaphological, orographic, and socioeconomic criteria determining sustainable cocoa cultivation were identified and mapped. Three modeling approaches (Analytic Hierarchy Process—AHP, Maximum Entropy—MaxEnt, and AHP—MaxEnt combined) were further used to hierarchize the importance of the criteria and to model the potential territory for sustainable cocoa cultivation. In all three modeling approaches, climatological criteria stood out among the five most important criteria. Elevation (orographic criteria) is also featured in this group. On the other hand, San Martin and Amazonas emerged as the five regions with the largest area ‘Highly suitable’ for cocoa cultivation in all three modeling approaches, followed by Loreto, Ucayali, Madre de Dios, Cusco, Junín, and Puno, which alternated according to modeling approach. From most to least restrictive, the AHP, MaxEnt, and AHP–MaxEnt modeling approaches indicate that 1.5%, 5.3%, and 23.0% of the Peruvian territory is ‘Highly suitable’ for cocoa cultivation, respectively.

1. Introduction

Cocoa (Theobroma cacao L.) is grown from 100 to 1400 m a.s.l., in landscapes ranging from mountains to alluvial plains with dry and pre-humid environments [1]. That is, in multiple edaphic, physiographic and climatic conditions, which originate a wide range of agroecological environments that respond differentially to technological recommendations and crop management options. Peru is considered one of the main producers and suppliers of fine cocoa, which in turn is the world’s second-largest producer of organic cocoa, with 48.6% on exports of cocoa beans, of which 20% has the organic and fair-trade certification. It is also the world’s eighth largest producer of cocoa beans, accounting for 1.7% of world cocoa bean production. Cocoa is also the second most important alternative to illegal crops, after coffee, which highlights its growing importance [2].

Nowadays, the expansion and modernization of the cocoa-growing area, under new production strategies and criteria of competitiveness and sustainability, require the zoning of the territory with appropriate biophysical conditions (climate, soil, orography) and infrastructure (accessibility, nearby populations, etc.), so as to facilitate the optimization in production [3]. Agricultural zoning for cocoa cultivation becomes particularly important in Peru as it is the fifth country with the largest number of cultivated hectares of cocoa in Latin America and the Caribbean [2].

In that sense, the integration of remote sensing (RS), Geographic Information Systems (GIS) and Multi-Criteria Evaluation (MCE) techniques is a support system for decision-making problems, such as agro-zoning, which include a large set of factors or constraints [4]. This system has been widely applied to model land suitability for sustainable agriculture and other environmental and socioeconomic sciences. For planning sustainable agriculture, it was applied, for instance, for coffee [5,6,7], rice [8], cocoa [9,10,11,12], and combined crops [13,14,15], among other crops.

Although there are different MCE techniques, the most common approach is to estimate the weight of importance of each criterion using expert opinions [16] by the Analytical Hierarchy Process (AHP) [17]. In Peru, previous studies have been reported using this technique on specific crops such as Malus domestica in the Mala Valley (Lima Region) [18] and Coffea arabica [19] and Solanum tuberosum [20] in Amazonas Region. However, expert opinions have a strong influence on the results, and the weights could be biased [21]. This concern is addressed by applying machine learning techniques to estimate the weight of each criterion. In this regard, Species Distribution Models (SDM), such as the Maximum Entropy (MaxEnt) approach [22], apart from modeling the potential distribution of suitable habitat for a given species, provide valuable information on the relative contribution of the criteria used to the model. This relative contribution may then be used as a weight of importance in the MCE [23].

In the literature, MCE and SDM have been integrated to determine suitable areas for sustainable aquaculture [24], priority areas for archaeological site protection [25], priority areas for species conservation [26] and vulnerable areas within a protected area [23]. However, there is no evidence of integrated modeling (AHP and MaxEnt) on land suitability for sustainable cocoa or other agricultural crop farming [27]. Therefore, in this study, land suitability for sustainable cocoa cultivation in Peru was determined using AHP-only, MaxEnt-only and combined AHP–MaxEnt modeling. To this end, three specific objectives were implemented: (i) to map the key criteria for sustainable cocoa cultivation in Peru, (ii) to rank the importance of the key criteria for sustainable cocoa cultivation and (iii) to model the current potential territory for sustainable cocoa cultivation in Peru.

2. Materials and Methods

2.1. Study Area and Methodology Framework

Peru is located in the intertropical zone of South America, with an altitudinal gradient from sea level to 6768 m a.s.l. (Figure 1). Peru covers an area of approximately 1,285,215 km², making it the twentieth-largest country on Earth and the third-largest in South America. It has an enormous landscape diversity due to its geographic conditions, which in turn gives it a great diversity of natural resources and agroecosystems. Peru, with a population of 31,237,385 inhabitants, a population density of 24.3 inhabitants/km² and an annual growth rate of 1.07%, is the fifth most populous country in South America [28]. Of the 24 regions of Peru, 12 have georeferenced cocoa records (Figure 1).

Figure 2 shows the methodological process to determine the suitability of the territory for cocoa cultivation in Peru. Three modeling approaches were used: (i) AHP, (ii) MaxEnt, and (iii) AHP—MaxEnt combination. The three approaches were worked in the GIS environment with spatial data, expert opinion and spatial statistics.

2.2. Identification and Mapping of Criteria for Cocoa Cultivation

Different criteria were used according to their availability in national and international spatial databases and the requirements of each modeling approach [27] (

Table 1. Criteria used for each modeling approach.

Criteria (1st Hierarchy AHP)	Subcriteria (2nd Hierarchy AHP)	Model in Which It Was Used			Spatial Database
		AHP	MaxEnt	AHP–MaxEnt
Climatological	Bio1: Mean annual temperature	x	x	x	[40]
	Bio2–Bio11: bioclimatics derived from temperature		x		[40]
	Bio12: Mean annual rainfall	x	x	x	[40]
	Bio13–Bio19: bioclimatics derived from precipitation		x		[40]
	Annual mean max temperature	x		x	[40]
	Annual mean min temperature	x		x	[40]
	Number of dry months	x		x	[40]
	Relative humidity	x	x	x	[41]
	Solar radiation		x		[40]
Edaphological at 0.30 m	pH in H2O	x	x	x	[42]
	Coarse fragment content	x	x	x	[42]
	Organic carbon	x	x	x	[42]
	Texture	x		x	[42]
	Total nitrogen	x	x	x	[42]
	CEC—Cation exchange capacity	x	x	x	[42]
	Bulk density		x		[42]
	Proportion of sand particles		x		[42]
	Proportion of silt particles		x		[42]
	Proportion of clay particles		x		[42]
Orographic	Elevation	x	x	x	CGIAR
	Slope	x	x	x	CGIAR
	Aspect	x	x	x	CGIAR
Socioeconomical	LULC—Land Use and Land Cover	x			[43,44]
	Distance to urban centers	x			[45]
	Distance to roads	x			[46]
	Distance to rivers	x			[45]
	Distance to protected natural areas	x			[47]

). Namely, the AHP approach requires knowledge of specific ranges of crop suitability for each criterion. Meanwhile, the MaxEnt approach does not need these ranges, and it is possible to use commonly unknown criteria, such as bioclimatic criteria. For AHP modeling, a set of 20 criteria that condition/favors cocoa cultivation were identified based on previous studies on agricultural land suitability [1,9,10,11,12,29] and technical manuals [30,31,32,33,34] on cocoa cultivation. For the MaxEnt modeling, 33 environmental criteria were established based on previous cocoa potential distribution modeling studies [35,36,37,38,39]. Of all these criteria, for AHP–MaxEnt modeling, socioeconomic criteria were excepted.

Spatial layers of precipitation, temperature and solar radiation, with 30″ spatial resolution, were obtained from WorldClim 2.1 [40]. A dry month was considered as a month where twice the monthly mean temperature was lower than the monthly precipitation, according to the Gaussen xerothermal index [48]. Monthly point data of relative humidity, with 10′ spatial resolution, were obtained from the Climatic Research Unit [41]. Nine interpolation techniques (inverse distance weighted, natural neighbor, spline: regularized and tension, ordinary kriging: spherical, circular, gaussian, linear and exponential) were used to generate continuous relative humidity maps (250 m spatial resolution) in ArcGIS 10.5 [20]. The best interpolation technique for each month was determined based on four statistics (coefficient of determination, mean bias error, mean absolute bias error, root mean square error and t-Student [49]), with 23% (1030) of the point data. The spline tension (11 months) and ordinary linear kriging (March) techniques performed best [50].

Soil physicochemical properties, with 250 m spatial resolution, were obtained from SoilGrids 2.0 [42]. Orographic variables were derived from the 250 m spatial resolution Digital Elevation Model, downloaded from the CGIAR Consortium for Spatial Information (www.srtm.csi.cgiar.org/; accessed on 17 January 2022). The Land Use and Land Cover (LULC) base map was obtained from the Copernicus Global Land Service-Land Cover (CGLS-LC100)-Collection 3-2019 at 100 m spatial resolution [43]. In this map, land uses (urban area, agricultural area and secondary vegetation) from the National Ecosystem Map of Peru [44,51], LULC maps from Ecological and Economic Zoning studies of 14/24 regions (Amazonas [52], Ayacucho, Cajamarca, Cusco, Huancavelica, Huánuco, Junín, Lambayeque, Madre de Dios, Piura, Puno, San Martin, Tacna and Ucayali) and from a local LULC map (province of Rodríguez de Mendoza [53]) were incorporated.

Urban polygons were extracted from the final LULC map, and population centers (points) were obtained from the Ministry of Education [45]. We used the road network from the Ministry of Transport and Communications [46] and the rivers from the 341 national charts of the National Geographic Institute [45]. Then, distances to roads, rivers and towns were calculated using Euclidean distance. We also used the protected areas and their buffer zones updated to 2022 by the National Service of Natural Areas Protected by the State [47].

In summary, 42 base layers were prepared in a raster model, with one thematic map for each sub-criterion of spatial suitability. These were standardized at a spatial resolution of 250 m and in the WGS84 geographic coordinate system.

2.3. Modelling Approach with the Analytical Hierarchical Process—AHP

2.3.1. Construction of Hierarchies and Thresholds of Criteria Suitability

In the AHP, the problem/objective is hierarchically structured into different levels comprising a predefined number of elements [54]. A hierarchy was constructed consisting of 20 sub-criteria (2nd hierarchy), grouped within four criteria (1st hierarchy) (Table 1). The subcriteria were reclassified according to thresholds of the suitability of the territory for cocoa cultivation (3rd hierarchy,

Table 2. Suitability thresholds of key criteria for cocoa cultivation in Peru, AHP model.

Criteria/Subcriteria		Land Suitability Classes (3rd Hierarchy AHP)
		Highly Suitable	Moderately Suitable	Marginally Suitable	Not Suitable	Adapted From
Climatological
Mean annual temperature (°C)		25–28	22–25/28–32	20–22/32–35	<20/>35	[10,29,30,31]
Annual mean min temperature (°C)		18–21	15–18/>21	12–15	<12	[11,12,32]
Annual mean max temperature (°C)		28–30	30–32/25–28	>32/22–25	<22	[11,12,32]
Mean annual rainfall (mm)		1600–2500	2500–3500/1400–1600	1200–1400/3500–4400	<1200/>4400	[30,31]
Number of dry months		0–2	3	4	>4	[12,32]
Relative humidity (%)		70–80	80–85/60–70	85–90/50–60	>90/<50	[11,33,34]
Edaphological at 0.30 m
pH in H2O		6–7	5–6/7–7.6	4.2–5/7.6–8.2	<4.2/>8.2	[29,30,31,32,34]
Texture 1		SiCL, CL, SiL	L, SCL, SC	Si, SL, C	LS, S, SiC	[30,31,32]
Coarse fragment content (%)		<15	15–35	35–55	>55	[10,30,31]
Organic carbon (%)		>1.5	0.8–1.5	<0.8	–	[30,31,32]
CEC (cmol/kg)		>24	20–24	16–20	<16	[10]
Total nitrogen (%)		>0.18	0.15–0.18	0.1–0.15	<0.1	[10]
Orographic
Elevation (m asl)		400–800	0–400/800–1200	1200–1600	>1600	[1,34]
Slope (%)		<8	8–16	16–30	>30	[29,30,31,32]
Aspect		N, NE, NW, Flat	W, E	SE, SW	S	[19,56]
Socioeconomical
LULC	CGLS-LC100 2	40	20	30	0, 50–90, >100	[19]
	Ecosystems of Peru	Agricultural area	–	Secondary vegetation	Urban/built	[20]
	Agricultural Map of Peru	Agriculture	–	–	–
	ZEE	Agriculture	–	Cattle raising	Urban/built	[20]
	Global urban borders	–	–	–	Urban/built
Distance to roads (km)	National-axis	0–6	6–9	9–12	>12	[19,20]
	Departmental	0–4	4–8	8–10	>10	[19,20]
	Local	0–2	2–4	4–8	>8	[19,20]
Distance to rivers (km)		0–0.5	0.5–2	2–5	>5	[19]
Distance to urban centers (km)	Urban areas	0–3	3–6	6–10	>10	[19]
	Population centers	0–1	1–3	3–5	>5	[19]
Distance to protected natural areas		Out	–	Buffer zone	Within	[19]

¹ S: Sand, LS: Loamy sand, SL: Sandy loam, L: Loam, SiL: Silt loam, Si: Silt, CL: Clay loam, SCL: Sandy clay loam, SiCL: Silty clay loam, SC: Sandy clay, SiC: Silty clay, C: Clay. ² CGLS-LC100 [43]: 0—No data, 20—Shrubs, 30—Herbaceous vegetation, 40—Cropland, 50—Urban/built up, 60—Bare/sparse vegetation, 70—Snow and ice, 80—water bodies, 90—Herbaceous wetland, and >100—all the forests.

). The commonly used approach to classifying land suitability thresholds is “FAO: A framework for land evaluation” [55]: Highly suitable, Moderately suitable, Marginally suitable, Currently unsuitable and Permanently unsuitable. In this study, as in other studies [18,19,20], the last two levels were combined since it is difficult to establish internal limits for these two levels.

2.3.2. Determination of Importance Weights of Criteria

The initial development of the first and second hierarchies required the construction of Pairwise Comparison Matrices (PCM), where cocoa experts compared one criterion against the others (pairwise) and established a degree of importance between them [7]. This section is not dealt with here but has been extracted and is discussed in the companion article to this one [27].

2.3.3. AHP Sub-Model Generation and AHP Suitability Modeling

The final development of the first and second hierarchies consisted of integrating the re-classified thematic maps (3rd hierarchy based on Table 2), according to the hierarchical group, by weighted superposition [19,20,57]. The resulting suitability depended on the reclassified map pixel score and the sub-criterion importance weight calculated by PCM. The integration of sub-criteria generated the climatological, edaphological, orographic and socioeconomic suitability sub-models, and the integration of these sub-models generated the final suitability model.

2.4. Modelling Approach with Maximum Entropy—MaxEnt

2.4.1. Georeferenced Cocoa Records

Georeferenced records were obtained from iNaturalist (www.inaturalist.org/observations; accessed on 9 January 2022), TROPICOS Missouri Botanical Garden (www.tropicos.org; accessed on 9 January 2022) and GBIF Global Biodiversity Information Facility (www.gbif.org/; accessed on 9 January 2022) through three QGIS 3.10 plugins (GBIF occurrences, Species Explorer and Natusfera) [35,36]. These were complemented with georeferenced records of native organic cocoa [58]. To remove spatial sampling bias and improve model performance [59], georeferenced records were filtered to a 250 m grid (equal to the spatial resolution of the criteria). The spatial filter reduced the georeferenced records from 546 to 196 (Figure 1).

2.4.2. Selection of Environmental Criteria

Collinearity between criteria causes overfitting problems, increases uncertainty, and decreases the statistical power of the model [60]. Therefore, using the ‘removeCollinearity’ function of the ‘virtualspecies’ package [61] in R 3.6, (i) Pearson’s correlation coefficients between criteria were calculated, from which (ii) a distance matrix was calculated, which in turn was used to (iii) construct a hierarchical cluster dendrogram. The criteria were grouped according to an r ≥ 0.7. This threshold is an acceptable measure to minimize multicollinearity of the adjusted models [60].

In order to select one important criterion per cluster, we ran a preliminary MaxEnt model (the setup is explained in Section 2.4.3) using all criteria, then we selected the criterion with the best performance in the Jackknife test [62] (i.e., the smallest difference in regularized training gains obtained from a model generated with all criteria except the criterion of interest, and a model generated with just the criterion of interest [63]). It was thus selected the following criteria, three orographic (elevation, slope and aspect), three bioclimatic (Bio04—Seasonality of temperature, Bio12—Annual precipitation, Bio19—Precipitation of the coldest quarter), and seven edaphological (CEC, Organic carbon, Bulk density, Total nitrogen and Coarse fragments, silt and, sand contents).

2.4.3. Modelling the Potential Distribution

The cocoa potential distribution model was generated by the Maximum Entropy principle algorithm [22], implemented in MaxEnt 3.4.4 (https://biodiversityinformatics.amnh.org/open_source/maxent/; accessed on 21 February 2022). 75% and 25% of the georeferenced records (randomly selected) were used for training and validation of each model, respectively [22]. The algorithm was run using 100 replicates over 1000 iterations with different random partitions (Bootstrap method), a convergence threshold of 0.00001 and 10000 maximum background points [63,64]. Other default settings were kept, as MaxEnt is able to select the appropriate function for the number of samples used for a model [60].

Model performance was evaluated by the Area Under the Curve (AUC), calculated from the Receiver Operating Characteristic [22]. Five levels of performance were differentiated according to the AUC [65]: excellent (>0.9), good (0.8–0.9), accepted (0.7–0.8), poor (0.6–0.7) and invalid (<0.6). The Cloglog output format of the model generated a map of continuous probability values for the potential cocoa distribution ranging from 0 to 1 [66]. These were reclassified into four ranges [63,64]: ‘Highly suitable’ (>0.6), ‘Moderately suitable’ (0.4–0.6) and ‘Marginally suitable’ (0.2–0.4) potential distribution, as well as Not suitable distribution (<0.2).

2.5. AHP–MaxEnt Modeling Approach

Reclassified thematic maps (based on Table 2) were integrated by weighted overlay [19,20,57]. The resulting suitability depended on the reclassified map pixel score and the sub-criterion importance weight. This weight, unlike the AHP model (Section 2.3.3), was not obtained by expert PCM (Section 2.3.2). For this model, a MaxEnt model was generated (the modeling setup was explained in Section 2.4.3), including criteria for this model (Table 1), obtaining the contribution percentage to the model. Then, this contribution percentage was assumed as the importance weight [23]. The integration of sub-criteria by weighted overlay [19,20,57] generated the final land suitability model for cocoa cultivation.

3. Results

3.1. Model Based on the Analytical Hierarchical Process—AHP

3.1.1. Suitability Map of Subcriteria

Figure 3 and Appendix A,

Table A1. Suitability area of subcriteria for cocoa cultivation in Peru, AHP model.

Criteria/Subcriteria	Highly Suitable		Moderately Suitable		Marginally Suitable		Not Suitable
	km2	%	km2	%	km2	%	km2	%
Climatological
Mean annual temperature	495,696.05	38.5	211,336.61	16.4	61,939.16	4.8	501,135.87	38.9
Annual mean min temperature	323,401.20	25.1	433,151.18	33.6	95,775.96	7.4	417,779.31	32.4
Annual mean max temperature	124,463.94	9.7	610,457.37	47.4	177,066.49	13.7	358,119.86	27.8
Mean annual rainfall	402,160.95	31.2	275,838.57	21.4	46,432.18	3.6	545,675.95	42.3
Number of dry months	970,538.09	75.3	39,724.27	3.1	38,972.61	3.0	220,872.67	17.1
Relative humidity	422,733.62	32.8	604,382.27	46.9	236,748.40	18.4	6319.39	0.5
Edaphological at 0.30 m
pH in H2O	284,809.84	22.1	333,425.91	25.9	615,584.49	47.8	36,363.58	2.8
Texture	463,633.52	36.0	534,435.04	41.5	250,084.29	19.4	22,031.38	1.7
Coarse fragment content	812,754.45	63.1	456,692.84	35.4	736.40	0.1	0.0	0.0
Organic carbon	1,178,174.42	91.4	78,533.32	6.1	13,476.01	1.0	0.0	0.0
CEC	190,558.84	14.8	221,422.78	17.2	418,014.65	32.4	440,187.65	34.2
Total nitrogen	1,025,291.05	79.6	144,885.17	11.2	79,769.70	6.2	20,238.07	1.6
Orographic
Elevation	112,490.36	8.7	662,751.90	51.4	56,421.60	4.4	438,500.97	34.0
Slope	814,094.16	63.2	210,170.69	16.3	199,388.82	15.5	46,510.98	3.6
Aspect	464,149.07	36.0	327,161.32	25.4	320,736.76	24.9	158,118.11	12.3
Socioeconomical
LULC	73,007.50	5.7	117,842.57	9.1	254,814.97	19.8	824,499.95	64.0
Distance to urban centers	259,032.66	20.1	323,477.46	25.1	160,753.82	12.5	526,900.84	40.9
Distance to roads	465,633.65	36.1	153,914.95	11.9	87,072.91	6.8	563,562.33	43.7
Distance to rivers	615,475.80	47.8	437,734.36	34.0	148,793.14	11.5	68,180.82	5.3
Distance to protected natural areas	911,644.79	70.7	0.00	0.0	133,880.37	10.4	224,658.54	17.4

show the reclassified maps and areas according to suitability thresholds (Table 2) of the climatological, edaphological, orographic and socioeconomic subcriteria. The subcriteria with the largest ‘Highly suitable’ area with respect to their criteria group are the number of dry months (970,538.09 km², 75.3%), organic carbon (1,178,174.42 km², 91.4%), slope (814,094.16 km², 63.2%) and protected areas (911,644.79 km², 70.7%). While those with the highest ‘Not suitable’ area are annual precipitation (545,675.95 km², 42.3%), CEC (440,187.65 km², 34.2%), elevation (438,500.97 km², 34.0%), and land cover and land use (824,499.95 km², 64.0%). In all maps, 1.4% (18,477.30 km²) of the Peruvian territory was discriminated from the analysis, corresponding to a mask of main water bodies, glaciers and urban areas.

3.1.2. Submodels and Land Suitability Models

In the AHP model, climatological (35.7%) and edaphological (29.1%) are the most important criteria, followed by socioeconomic (18.2%) and orographic (17.0%) (Companion article to this one [27]). On the other hand, the subcriteria, annual precipitation, CEC, elevation and distance to the water network scored the highest weighting with respect to their group of criteria (

Table 3. Weight of importance (%) and relative contribution of subcriteria to land suitability modeling for cocoa cultivation.

AHP Model 1,2			MaxEnt Model 1			AHP–MaxEnt Model 1
Bio12: Mean annual rainfall	9.9		Bio19: Precipitation of coldest quarter	41.4		Annual mean min temperature	40.7
Elevation	9.7		Elevation	23.6		Number of dry months	14.2
CEC	7.1		Bio12: Mean annual rainfall	6.0		Elevation	10.4
Relative humidity	6.7		Bio4: Temperature seasonality	5.7		Relative humidity	7.4
Texture	6.0		Slope	4.2		pH in H2O	5.6
Number of dry months	5.8		Proportion of silt particles	3.0		Bio1: Mean annual temperature	4.2
Total nitrogen	5.5		Coarse fragment content	3.0		Bio12: Mean annual rainfall	3.3
Distance to rivers	5.0		Aspect	3.0		Aspect	2.5
Bio1: Mean annual temperature	4.6		Organic carbon	2.6		CEC	2.1
Annual mean min temperature	4.6		CEC	2.5		Slope	2.1
Slope	4.6		Bulk density	2.0		Texture	2.0
pH in H2O	4.5		Proportion of sand particles	1.5		Organic carbon	1.7
Annual mean max temperature	4.1		Total nitrogen	1.5		Coarse fragment content	1.5
Organic carbon	4.1		Bio1: Mean annual temperature	*		Total nitrogen	1.3
LULC	4.1		Bio2, Bio3, Bio5–Bio11, Bio13–Bio18	*		Annual mean max temperature	1.1
Distance to roads	3.6		Relative humidity	*
Distance to PNA	3.4		Solar radiation	*
Aspect	2.7		pH in H2O	*
Distance to urban centers	2.1		Proportion of clay particles	*
Coarse fragment content	2.0

¹Italics = Climatological; Bold = Edaphological; Underlined = Orographic; Normal = Socioeconomical. ² Adapted from [27]. * Initially considered but removed from final suitability modeling by modeling approach.

). With the weighted overlay of sub-criteria, suitability submodels were generated for each hierarchical group. Indeed, climatological (363,379.34 km², 28.2%) and edaphological (290,845.16 km², 22.6%) are the sub-models with the highest ‘Highly suitable’ areas for cocoa cultivation (Figure 4).

A weighted overlay of sub-models generated an AHP model of land suitability for cocoa cultivation in Peru (Figure 5a). Here, 1.5% (19,437.63 km²), 80.6% (1,038,036.17 km²), 16.5% (211,982.87 km²), and 0.05% (630.04 km²) showed, respectively ‘Highly suitable’, ‘Moderately suitable’, ‘Marginally suitable’ and ‘Not suitable’ territory for cocoa cultivation (Appendix A,

Table A2. Land suitability for cocoa cultivation in Peruvian regions, AHP model.

Regions	Highly Suitable		Moderately Suitable		Marginally Suitable		Not Suitable		Non-Classified
	km2	%	km2	%	km2	%	km2	%	km2	%
Amazonas	2627.36	6.7	34,035.21	86.6	2366.54	6.0	0.00	0.0	277.35	0.7
Ancash	0.00	0.0	24,959.06	69.4	9987.54	27.8	0.00	0.0	1015.65	2.8
Apurímac	0.00	0.0	14,305.92	67.8	6729.24	31.9	0.00	0.0	78.99	0.4
Arequipa	0.00	0.0	21,537.55	34.0	40,332.23	63.8	496.06	0.8	890.04	1.4
Ayacucho	12.82	0.0	27,761.70	63.8	15,582.39	35.8	3.90	0.0	143.01	0.3
Cajamarca	255.28	0.8	26,147.56	79.1	6502.76	19.7	0.00	0.0	139.07	0.4
Callao	0.00	0.0	51.47	36.4	7.45	5.3	0.00	0.0	82.49	58.3
Cusco	2351.81	3.3	52,497.77	72.8	16,216.07	22.5	0.00	0.0	1010.49	1.4
Huancavelica	0.00	0.0	16,789.43	76.1	5165.07	23.4	0.00	0.0	110.54	0.5
Huánuco	1928.73	5.2	31,166.61	83.8	3852.00	10.4	0.00	0.0	253.19	0.7
Ica	0.00	0.0	8653.51	41.0	12,269.83	58.2	0.79	0.0	156.63	0.7
Junín	2128.12	4.8	36,898.31	83.9	4512.84	10.3	0.00	0.0	458.03	1.0
La Libertad	0.00	0.0	15,465.58	61.1	9700.39	38.3	0.00	0.0	130.00	0.5
Lambayeque	0.00	0.0	7683.17	53.6	6359.99	44.3	0.00	0.0	299.15	2.1
Lima	0.00	0.0	20,490.05	58.6	13,274.72	37.9	0.00	0.0	1225.23	3.5
Loreto	708.86	0.2	367,167.06	97.9	331.10	0.1	0.00	0.0	6908.97	1.8
Madre De Dios	1340.31	1.6	82,533.63	97.0	640.25	0.8	0.00	0.0	531.67	0.6
Moquegua	0.00	0.0	3993.39	25.3	11,551.11	73.1	82.89	0.5	179.92	1.1
Pasco	646.19	2.7	21,202.58	87.9	2092.50	8.7	0.00	0.0	172.68	0.7
Piura	0.00	0.0	20,674.83	57.3	13,867.65	38.5	0.00	0.0	1522.58	4.2
Puno	4.57	0.0	53,726.91	79.1	13,068.23	19.2	0.00	0.0	1163.11	1.7
San Martin	4732.75	9.3	40,866.60	80.2	5099.64	10.0	0.00	0.0	262.27	0.5
Tacna	0.00	0.0	3835.20	23.8	11,919.89	74.1	46.4	0.3	281.58	1.8
Tumbes	0.00	0.0	4169.10	88.9	410.28	8.7	0.00	0.0	110.89	2.4
Ucayali	2700.82	2.6	101,423.98	96.3	143.18	0.1	0.00	0.0	1073.78	1.0
Perú	19,437.63	1.5	1,038,036.17	80.6	211,982.87	16.5	630.04	0.05	18,477.30	1.4

). Regarding regions, San Martin (4732.75 km²), Ucayali (2700.82 km²), Amazonas (2627.36 km²), Cusco (2351.81 km²), Junín (2128.12 km²), Huánuco (1928.73 km²), and Madre De Dios (1340.31 km²) have the largest cocoa-growing areas with ‘Highly suitable’ land, on the contrary, Loreto (708.86 km²), Pasco (646.19 km²), Cajamarca (255.28 km²), Ayacucho (12.82 km²), and Puno (4.57 km²) have the smallest areas (

Table 4. ‘Highly suitable’ area (km²) of land suitability for cocoa cultivation in Peru, based on regions and modeling approaches.

AHP Model			MaxEnt Model			AHP–MaxEnt Model
San Martin	4732.75		Madre De Dios	26,285.36		Ucayali	92,224.56
Ucayali	2700.82		Loreto	19,827.18		Madre De Dios	75,261.35
Amazonas	2627.36		San Martin	9829.90		Loreto	60,526.92
Cusco	2351.81		Amazonas	4107.07		San Martin	18,009.91
Junín	2128.12		Puno	3149.39		Amazonas	11,913.73
Huánuco	1928.73		Cajamarca	1694.55		Cusco	10,873.61
Madre De Dios	1340.31		Junín	774.32		Huánuco	9342.31
Loreto	708.86		Cusco	751.75		Pasco	7173.36
Pasco	646.19		Ucayali	485.92		Puno	5435.40
Cajamarca	255.28		Pasco	450.47		Junín	5109.47
Ayacucho	12.82		Huánuco	174.66		Cajamarca	543.78
Puno	4.57		Tumbes	103.31		Ayacucho	131.31
Ancash	0		Lima	51.86		Ancash	0
Apurímac	0		Ancash	39.34		Apurímac	0
Arequipa	0		La Libertad	25.58		Arequipa	0
Callao	0		Piura	19.48		Callao	0
Huancavelica	0		Ayacucho	16.90		Huancavelica	0
Ica	0		Callao	0.12		Ica	0
La Libertad	0		Arequipa	0.05		La Libertad	0
Lambayeque	0		Apurímac	0		Lambayeque	0
Lima	0		Huancavelica	0		Lima	0
Moquegua	0		Ica	0		Moquegua	0
Piura	0		Lambayeque	0		Piura	0
Tacna	0		Moquegua	0		Tacna	0
Tumbes	0		Tacna	0		Tumbes	0

).

3.2. Maximum Entropy Model—Maxent

3.2.1. Model Performance and Importance of Sub-Criteria

The average AUC over 100 MaxEnt replicates 0.916, with a standard deviation of 0.008, indicating an excellent predictive performance of the model. According to the Jackknife test of variable importance, the Elevation is the environmental variable with the highest gain when used in isolation; it, therefore, seems to have the most useful information itself. Elevation, when omitted, is also the environmental variable decreasing the gain the most, and therefore appears to have the most information missing in the other variables. It was found that 76.7% of the potential cocoa distribution is driven by four environmental variables, namely Bio19–Precipitation of the coldest quarter (41.4), Elevation (23.6%), Bio12–Annual precipitation (6.0%) and Bio04–Seasonality of temperature (5.7%) (Table 3). At the same time, Sand Content (1.5%) and Total Nitrogen (1.5%) contributed the least.

3.2.2. Potential Distribution

Areas of ‘high’ potential distribution probability for cocoa were identified mainly in the lowlands of the Peruvian Amazon (Figure 5b). Areas of ‘high’, ‘moderate’, ‘low’ and ‘no potential’ cocoa distribution cover 5.3% (67,787.22 km²), 7.2% (92,791.09 km²), 20.3% (261,335.27 km²) and 65.8% (848,173.42 km²) of Peru’s territory, respectively (Appendix A,

Table A3. Potential distribution of cocoa in Peruvian regions, MaxEnt model.

Regions	Highly Suitable		Moderately Suitable		Marginally Suitable		Not Suitable		Non-Classified
	km2	%	km2	%	km2	%	km2	%	km2	%
Amazonas	4107.07	10.4	8039.43	20.5	11,974.36	30.5	14,908.26	37.9	277.35	0.7
Ancash	39.34	0.1	80.84	0.2	231.55	0.6	34,594.89	96.2	1015.63	2.8
Apurímac	0.00	0.0	0.00	0.0	0.00	0.0	21,035.16	99.6	78.99	0.4
Arequipa	0.05	0.0	0.45	0.0	5.34	0.0	62,360.05	98.6	890.00	1.4
Ayacucho	16.90	0.0	161.15	0.4	516.10	1.2	42,666.66	98.1	143.02	0.3
Cajamarca	1694.55	5.1	1000.71	3.0	2436.93	7.4	27,773.40	84.0	139.07	0.4
Callao	0.12	0.1	0.06	0.0	1.27	0.9	57.46	40.6	82.49	58.3
Cusco	751.75	1.0	3591.87	5.0	6685.78	9.3	60,036.25	83.3	1010.48	1.4
Huancavelica	0.00	0.0	0.00	0.0	4.8	0.0	21,949.70	99.5	110.54	0.5
Huánuco	174.66	0.5	1651.59	4.4	7839.97	21.1	27,281.12	73.3	253.18	0.7
Ica	0.00	0.0	0.23	0.0	3.36	0.0	20,920.58	99.2	156.60	0.7
Junín	774.32	1.8	3557.1	8.1	6441.02	14.6	32,766.83	74.5	458.02	1.0
La Libertad	25.58	0.1	88.28	0.3	267.96	1.1	24,784.18	98.0	129.98	0.5
Lambayeque	0.00	0.0	9.32	0.1	444.94	3.1	13,588.92	94.7	299.13	2.1
Lima	51.86	0.1	115.09	0.3	311.33	0.9	33,286.51	95.1	1225.20	3.5
Loreto	19,827.18	5.3	31,559.59	8.4	140,822.04	37.5	175,998.12	46.9	6909.06	1.8
Madre De Dios	26,285.36	30.9	23,172.08	27.2	22,018.65	25.9	13,038.19	15.3	531.59	0.6
Moquegua	0.00	0.0	0.00	0.0	0.76	0.0	15,626.63	98.9	179.92	1.1
Pasco	450.47	1.9	2859.58	11.9	6003.89	24.9	14,627.35	60.7	172.67	0.7
Piura	19.48	0.1	132.9	0.4	1079.19	3.0	33,310.93	92.4	1522.57	4.2
Puno	3149.39	4.6	1572.31	2.3	2668.18	3.9	59,409.85	87.4	1163.09	1.7
San Martin	9829.90	19.3	9625.48	18.9	13,939.97	27.4	17,303.63	34.0	262.27	0.5
Tacna	0.00	0.0	0.00	0.0	0.12	0.0	15,801.38	98.2	281.57	1.8
Tumbes	103.31	2.2	110.1	2.3	382.53	8.2	3983.41	84.9	110.91	2.4
Ucayali	485.92	0.5	5462.94	5.2	37,255.24	35.4	61,063.97	58.0	1073.69	1.0
Perú	67,787.22	5.3	92,791.09	7.2	261,335.27	20.3	848,173.42	65.8	18,477.03	1.4

). On the regional distribution, Madre De Dios (26,285.36 km²), Loreto (19,827.18 km²), San Martin (9829.90 km²), Amazonas (4107.07 km²), Puno (3149.39 km²), and Cajamarca (1694.55 km²) have the largest areas with ‘high’ potential distribution for cocoa cultivation (Table 4).

3.3. Model Based on AHP–MaxEnt

3.3.1. Importance and/or Weights of Subcriteria

The average AUC for the 100 MaxEnt replicates 0.920, and the standard deviation is 0.007, suggesting an excellent predictive performance of the model. It was found that 72.7% of the potential cocoa distribution is driven by four environmental variables, namely, Mean annual minimum temperature (40.7%), Number of dry months (14.2%), Elevation (10.4%) and Relative humidity (7.4%) (Table 3). While Coarse Fragment Content (1.5%), Total Nitrogen (1.3%), and Mean Annual Maximum Temperature (1.1%) contributed the least. The subcriteria, mean annual minimum temperature (40.7%), pH in H₂O (5.6%), and elevation (10.4%) had the highest weighting within their criteria group; not-withstanding, mean annual maximum temperature (1.1%), Total nitrogen (1.3%), and Terrain slope (2.1%) were the least weighted.

3.3.2. Land Suitability Model

With the weighted overlay of subcriteria, it was generated the land suitability model for cocoa cultivation in Peru (Figure 5c). In Peru, 23.0% (296,545.69 km²), 37.4% (482,489.88 km²), 35.2% (453,379.97 km²), and 2.9% (37,671.17 km²) of the territory featured ‘Highly suitable’, ‘Moderately suitable’, ‘Marginally suitable’, and ‘Not suitable’, respectively (Appendix A,

Table A4. Land suitability for cocoa cultivation in Peruvian regions, AHP–MaxEnt model.

Regions	Highly Suitable		Moderately Suitable		Marginally Suitable		Not Suitable		Non-Classified
	km2	%	km2	%	km2	%	km2	%	km2	%
Amazonas	11,913.73	30.3	16,080.17	40.9	11,035.21	28.1	0.00	0.0	277.35	0.7
Ancash	0.00	0.0	3698.66	10.3	31,241.36	86.9	6.57	0.0	1015.65	2.8
Apurímac	0.00	0.0	28.43	0.1	20,915.26	99.1	91.48	0.4	78.99	0.4
Arequipa	0.00	0.0	2344.36	3.7	38,506.10	60.9	21,515.40	34	890.04	1.4
Ayacucho	131.31	0.3	1424.91	3.3	38,317.92	88.1	3486.68	8.0	143.01	0.3
Cajamarca	543.78	1.6	6555	19.8	25,801.48	78.1	5.33	0.0	139.07	0.4
Callao	0.00	0.0	33.31	23.6	25.61	18.1	0.00	0.0	82.49	58.3
Cusco	10,873.61	15.1	12,867.98	17.9	46,907.22	65.1	416.84	0.6	1010.49	1.4
Huancavelica	0.00	0.0	139.20	0.6	21,609.48	97.9	205.81	0.9	110.54	0.5
Huánuco	9342.31	25.1	8066.40	21.7	19,533.25	52.5	5.38	0.0	253.19	0.7
Ica	0.00	0.0	1985.98	9.4	18,077.74	85.8	860.42	4.1	156.63	0.7
Junín	5109.47	11.6	12,203.67	27.7	26,226.12	59.6	0.00	0.0	458.03	1.0
La Libertad	0.00	0.0	5857.69	23.2	19,301.92	76.3	6.36	0.0	130.0	0.5
Lambayeque	0.00	0.0	12,371.14	86.3	1670.14	11.6	1.89	0.0	299.15	2.1
Lima	0.00	0.0	1266.71	3.6	32,163.85	91.9	334.21	1.0	1225.23	3.5
Loreto	60,526.92	16.1	307,662.24	82.0	17.90	0.0	0.00	0.0	6908.94	1.8
Madre De Dios	75,261.35	88.5	8251.35	9.7	1001.44	1.2	0.00	0.0	531.72	0.6
Moquegua	0.00	0.0	218.60	1.4	9533.01	60.3	5875.78	37.2	179.92	1.1
Pasco	7173.36	29.7	7435.32	30.8	9332.53	38.7	0.06	0.0	172.68	0.7
Piura	0.00	0.0	26,559.96	73.6	7978.70	22.1	3.84	0.0	1522.58	4.2
Puno	5435.40	8.0	6332.45	9.3	54,856.01	80.7	175.86	0.3	1163.10	1.7
San Martin	18,009.91	35.3	23,806.63	46.7	8882.46	17.4	0.00	0.0	262.27	0.5
Tacna	0.00	0.0	739.35	4.6	10,382.87	64.6	4679.27	29.1	281.57	1.8
Tumbes	0.00	0.0	4538.01	96.8	41.36	0.9	0.00	0.0	110.89	2.4
Ucayali	92,224.56	87.5	12,022.37	11.4	21.05	0.0	0.00	0.0	1073.78	1.0
Perú	296,545.69	23	482,489.88	37.4	453,379.97	35.2	37,671.17	2.9	18,477.30	1.4

). Regionally, Ucayali (92,224.56 km²), Madre De Dios (75,261.35 km²), Loreto (60,526.92 km²), San Martin (18,009.91 km²), Amazonas (11,913.73 km²), Cusco (10,873.61 km²), and Huánuco (9342. 31 km²) have the largest areas with ‘Highly suitable’ land for cocoa cultivation, compared to Pasco (7173.36 km²), Puno (5435.40 km²), Junín (5109.47 km²), Cajamarca (543.78 km²), and Ayacucho (131.31 km²) showing the smallest ‘Highly suitable’ areas (Table 4).

4. Discussion

Cropland suitability analysis based on different modeling approaches has been well documented [67,68,69,70], and despite the potential gains achieved for crop zoning on an individual basis with MCE approaches such as AHP [19,20,71], and with SDM approaches such as MaxEnt [23,72], the integration of both models has recently become an important tool to enhance various factors of importance [73,74], regardless of individual suitability adjustment values. Therefore, for the first time in this research, high potential suitability cocoa lands in Peru are documented based on three models with (i) AHP, (ii) Maxent and (iii) AHP–MaxEnt approach.

There were 42 hierarchical key sub-criteria for sustainable cocoa cultivation in Peru, including 20, 33, and 15 for the AHP, MaxEnt, and AHP–MaxEnt modeling approaches, respectively. Although each model has different evaluation criteria [75], the three approaches showed similarities in their results regarding the most important criterion. Climatological criteria stood out in the top four positions of the most important criteria in all three modeling approaches. Elevation (orographic criterion) is also featured in this group. Differences in criteria used to respond to the need for input information in each approach. Namely, the EMC agro-zoning approaches, such as AHP, require input on specific ranges of crop suitability for each criterion. The commonly used classification guide for land suitability thresholds is the “FAO: Framework for Land Evaluation” [55], and ranges exist for a wide list of crops [30,31,32]. Meanwhile, machine learning modeling (~SDM) approaches such as MaxEnt has no need for these ranges, and it is possible to use a larger number of commonly unknown criteria, such as WorldClim’s bioclimatic criteria [40].

However, when using the MaxEnt approach, individually or in combination, it is advisable not to include socioeconomic variables, as we did in this study. Since this species distribution modeling is based on points of occurrence of the species in naturally suitable areas, with no human interaction. In common platforms (GBIF, iNaturalist, TROPICOS, speciesLink and others) for obtaining occurrence data for species modeling, there is no differentiation between wild and cultivated collections regarding crops [39]. Socioeconomic variables, though not supported by the MaxEnt model features, are still relevant and could be used as a restriction mask for MaxEnt and AHP–MaxEnt results.

It is assessed that the most suitable cocoa areas in Peru are mainly explained by the climatological criteria and elevation in the three approaches. Compared to previous studies on cocoa land suitability, using MCE [1,9,10,11,12,29] or SDM [35,36,37,38,39], this study included a greater number of sub-criteria (42 sub-criteria). On the one hand, this is because of the three approaches used, and on the other hand, in such studies, the large range of subcriteria depends on the study scope and spatial data availability [19,20]. In future studies, for example, economic (benefit-cost, productivity, crop rate of return or other [76]) and social (household skill level, labor availability, access to information, poverty rate or other [77]) subcriteria, not considered here for unavailable spatial data, may be included; and of course, crop risk maps such as disease [19,36] or cadmium (Cd) in the soil [78] may be incorporated, especially Cd due to its detrimental impact on cocoa [79]. However, an important issue when integrating more criteria is to consider that there will be a greater spatial heterogeneity of the data sources, which will ultimately influence the results [80].

The AHP approach determined that Bio12—Annual precipitation, elevation, and CEC are the top three sub-criteria predicting the model with 26.7% contribution. This approach allowed flexible decision-making by groups of sub-criteria, which can easily use both intangible and tangible variables in a systematic way [81]. In Peru, for example, this approach provides a structured and comparatively simple solution to multi-criteria decision-making problems on cocoa crop suitability. Despite the fact that the importance weights estimated by experts in the AHP approach may be influenced by respondents’ subjectivity [21], the findings demonstrate a much more homogeneous distribution of weights, and the importance is no longer concentrated on three or four criteria such as the other two approach-es using MaxEnt’s machine learning [22].

Since MaxEnt effectively addresses the suitability effects of variable environmental factors [82], it allowed us to quantitatively relate the model to potential areas [83], whereby 5.3% of the territory is highly suitable for cocoa cultivation based on 18 variables. This model showed that the combined contribution of the variables Bio19—Coolest quarter precipitation, elevation, and Bio12—Annual precipitation reached up to 71%.

In the third modeling approach (AHP–MaxEnt), the sub-criteria of minimum mean annual temperature, number of dry months, and elevation are matched by a 65.3% contribution to the iterative process with 15-variables model building. Thus, AHP–MaxEnt successfully addressed the uncertainty of expert opinion weighting bias when using the AHP-only model [23] while further boosting the climatic information provided by the MaxEnt model, thereby identifying up to 23% of highly suitable areas for cocoa cultivation in Peru.

The results of the combined AHP–MaxEnt approach show higher fractions of ‘Highly suitability’ areas, compared to the AHP and MaxEnt only, because assembly allowed adjusting the value of moderately suitable areas with respect to the subcriteria used in the AHP model. By assuming the expert criteria, the farmer can improve agricultural practices to enhance yields in highly suitable areas with techniques that improve the CEC and achieve values greater than 24 cmol/kg [30,31,32] in optimal cocoa altitudes between 400–800 [1,34], and with drainage practices, infiltration ditches or installation of technician irrigation if necessary to meet cocoa optimal water requirement between 1600–2500 mm per year [30,31].

In the three modeling approaches, San Martin and Amazonas were among the five regions with the largest ‘Highly Suitable’ area for cocoa cultivation, followed by Loreto, Ucayali, Madre de Dios, Cusco, Junín and Puno, having alternating positions according to the modeling approach. These regions have also recorded the highest production but have not the highest yields; even Cusco (366 kg/ha) and Amazonas (642 kg/ha) lie below the national yield average of 720 kg/ha [2]. Furthermore, despite being areas with lower Cd estimated in soil [78], forest losses also affect these regions, about 70% are patches of less than 5 ha (small-scale agriculture) [84].

According to the AHP approach, the regions have a significantly smaller ‘Highly Suitable’ area. It may be due to the socioeconomic criteria only used in this approach, restricting the naturally suitable areas (climate, edaphology and orography) as identified in the other approaches, but with no accessibility or infrastructure conditions. However, although limited spatial information may constrain crop-land suitability assessment, future modeling studies could include new variables that influence socioeconomic performance, such as farm size, expertise in cocoa, and partnership involvement in associations [3,85]. This will also allow taking advantage of the potential of the AHP approach to assigning weights/scores and add criteria and alternatives to important social, political, economic and technical variables, and a variety of objectives, criteria and alternatives [81], thus exploiting the potential of the ensemble models [3].

From most to least restrictive, the AHP, MaxEnt, and AHP–MaxEnt modeling approaches indicate that 1.5% (19,437.63 km²), 5.3% (67,787.22 km²), and 23.0% (296,545.69 km²) of the Peruvian territory is ‘Highly suitable’ for cocoa cultivation, respectively. However, the marked difference between these areas may be due to the different criteria and modeling contribution weights. Therefore, future studies would identify whether the high differences hold when using the same criteria for the three modeling approaches. The MaxEnt and AHP–MaxEnt approaches present a greater ‘Highly suitable’ range in the Amazonian regions of Peru. Namely, the centers of origin of cocoa are located in South America’s Amazon [86], and for these approaches, we used cocoa collection points of occurrence from that area. Notwithstanding, the AHP model discriminated the ‘Highly suitable’ areas because, in the Amazon, there are currently no conditions of accessibility or infrastructure for cocoa production (AHP—socioeconomic variables).

From an economic production approach, it is recommended to use the most restrictive model for the success of the crop. On the other hand, regarding the conservationist approach, for germplasm collection and/or genetic conservation purposes, it is suggested to use the least restrictive model in order to study a larger area and apply a precautionary principle [62]. Furthermore, using the area (4013.21 km²) identified as ‘Highly suitable’ in the three models is also recommended to ensure crop cultivation success.

The gap between the statistic of 1300 km² of cocoa cultivated area in Peru [2] and the estimated potential ‘Highly Suitable’ area in this study is significant. Yet, a national cocoa cultivated area map is needed to identify the regions with a spatial gap. In fact, 44.0% (8560.41 km²), 14.9% (10,070.01 km²), and 7.0% (20,759.83 km²) of the ‘Highly suitable’ area from the AHP, MaxEnt, and AHP–MaxEnt modeling approaches matches the national agricultural area map (11,6497.16 km²) [87], respectively. This suggests that currently, non-cocoa agricultural areas may also be reconverted to cocoa farms.

5. Conclusions

There were 42 hierarchical key sub-criteria for sustainable cocoa cultivation in Peru, including 20, 33, and 15 for the AHP, MaxEnt, and AHP–MaxEnt modeling approaches, respectively. This sub-criteria were grouped into climatological, edaphological, orographic, and socioeconomic criteria. Indeed, climatological criteria stood out among the top four most important criteria in the three modeling approaches. Elevation (orographic criterion) is also featured in this group. San Martin and Amazonas regions had the largest area ‘Highly suitable’ for cocoa cultivation among the top five regions, according to the three modeling approaches. These two regions were followed by Loreto, Ucayali, Madre de Dios, Cusco, Junín and Puno, which alternated depending on the modeling approach. From most to least restrictive, the AHP, MaxEnt, and AHP–MaxEnt modeling approaches report that 1.5%, 5.3%, and 23.0% of the Peruvian territory is ‘Highly suitable’ for cocoa cultivation, respectively.

The study will provide decision support for sustainable agricultural cocoa production in Peru, as well as an opportunity to improve agricultural planning by providing much-needed information to farmers and agricultural planners. The methodological approach used in this research integrates AHP and MaxEnt for land suitability analysis for cocoa cultivation, and it can definitely be applied to other cocoa-growing areas of the world, with the appropriate adjustments to local realities. This methodology can also be applied to other crops of nutritional, economic and environmental importance in Peru. The land suitability analysis identifies areas with suitable crop development, contributing in this way to not overexploiting soil resources and, consequently, practicing sustainable agriculture.

This study has shown that SDM (particularly MaxEnt) could be used together with MCE models (specifically AHP) in a complementary approach, providing a more robust method for land evaluation for agriculture. Additional case studies would be advantageous, and there is also the potential to explore other SDMs in addition to Maxent. The SDM provides additional information to support the MCE approach that would otherwise be difficult to acquire.