Crop pollination by wild bees

Coffee flowers for a few days after the first rains, and the bees that visit it nest somewhere else: in the ground of a pasture, in the hollow trunks of a forest, in dead wood. How many reach a farm depends on the land around it, which makes pollination a service of the landscape, not of the field. In this tutorial you will build InVEST's crop pollination model in Terraflux on Ribeirão Claro, with its cropland as coffee: you will describe three guilds of wild bees and the resources each land use offers them, map where they come from and where they forage, put sixteen coffee farms on the map and read their pollination and yield, see how the bees' flight distance and the crop's half saturation change the answer, and compare the scenario of 2040 with 2020. At the end you will know what the model's maps mean, where its numbers come from, and how not to over-read them.

The idea#

Most crops grown for fruit or seed yield more when animals pollinate them. Klein et al. (2007) found that 87 of the leading global food crops depend on animal pollination to some degree, while 28 do not; Ricketts et al. (2004) measured coffee yields 20% higher within about 1 km of a forest in Costa Rica, thanks to the bees that live in it. Wild bees are central-place foragers: they nest in one place and fly out to the flowers around it, so a farm gets bees only if nesting sites and flowers lie within flight range of it. Kremen et al. (2007) put this into a conceptual framework for services delivered by mobile organisms, and Lonsdorf et al. (2009) turned it into the first quantitative model that predicts the abundance of pollinators across a landscape from land use. That model is InVEST's crop pollination model.

The model works with guilds: groups of species that nest and forage alike (ground-nesting bees, bees nesting in cavities, bumble bees...), or with single species when they are known well. For each guild it asks four questions of every cell:

1. Can it nest here? Each land use offers some availability of each nesting substrate (bare ground, cavities in trees, stems...), and each guild uses some substrates.
2. What flowers are within its reach? Each land use has floral resources in each season, and each guild is more or less active in each season. The flowers around a cell are weighted by distance, with a kernel set by the guild's typical foraging distance.
3. How many bees does it send out? The pollinator supply: nesting sites times the flowers within reach.
4. How many bees forage here? The pollinator abundance: the bees nesting around a cell, shared out among the flowers they can reach.

The kernel \(e^{-D/\alpha}\) gives every cell a weight that falls with distance: at \(D = \alpha\) a flower counts 37% of one next to the nest, at \(2\alpha\) 14%. Terraflux, as InVEST's code, computes it in cells (\(\alpha\) divided by the cell size), cuts it at \(5\alpha\), where the weight is below 1%, and divides it by its sum.

Enlarge

A worked cell: take the stingless bees (nesting only in cavities, equally active in both seasons, so \(fa = 0.5\) in each; \(sa = 0.4\)) and a coffee cell whose flight range is half forest and half coffee, the flowers being the same everywhere in it. With this tutorial's table, forest offers \(FE = 0.5 \times 0.7 + 0.5 \times 0.5 = 0.6\) and coffee \(0.5 \times 0.6 + 0.5 \times 0.1 = 0.35\), so \(FR = 0.475\). The forest's cells supply \(PS = 0.475 \times 1.0 \times 0.4 = 0.19\) bees, the coffee's \(0.475 \times 0.05 \times 0.4 = 0.0095\): a mean of 0.0998 within range. In the wet season the coffee cell gets \(PA = (0.6 \times 0.5 / 0.475) \times 0.0998 = 0.063\). On the real map the stingless bees' mean abundance on the coffee within 60 m of forest is 0.065.

On the farms#

With farms, polygons of a crop, the model turns the abundance into pollination and yield. Each farm names the season its crop is pollinated, its half-saturation constant \(h\), the share of its pollination given by managed bees (hives), and how much of its yield depends on pollinators.

Enlarge

Farm 1 of this tutorial has a mean abundance \(PAT = 0.189\) in the wet season, \(h = 0.05\), no hives and \(\nu = 0.25\): \(FP = 0.189 \times 0.95 / (0.05 \times (1 - 0.378) + 0.189) = 0.816\), so \(Y_W = 0.25 \times 0.816 = 0.204\) and \(Y_T = 1 - 0.25 \times (1 - 0.816) = 0.954\). Without any bee it would still yield \(1 - 0.25 = 0.75\) of its potential: \(\nu\) caps what pollinators can add.

When to use it#

• Where do a crop's pollinators come from? The supply maps show which forest, pasture or field margins feed the farms, and so which patches are worth keeping for agriculture's sake.
• Scenarios. A plan that clears forest edges, plants a crop or restores riparian forest changes nesting and flowers at once; the model says what it does to the bees on the farms.
• Farm practices. Flower strips, hedgerows or nest boxes change a farm's own floral and nesting values, which the farm polygons carry.
• Beside the other services, as in the ideas of this series: carbon, habitat and water for the same scenario, and pollination for its agriculture.

Lonsdorf et al. (2009) tested the model against field data: coffee in Costa Rica, watermelon and sunflower in California, watermelon in New Jersey and Pennsylvania. In Costa Rica and California it explained up to 80% of the variance among farms in the bees observed; in New Jersey and Pennsylvania it did not predict them, probably because the land use map missed the small patches of flowers and nesting sites around the farms. Their sensitivity analyses found the predictions depend strongly on the resources inside the crop fields themselves. Ricketts et al. (2008) synthesised 23 studies of 16 crops: the visitation of native bees falls exponentially with the distance from natural habitat, to half at 0.6 km, more steeply in the tropics.

What it is not: it is an index of relative abundance, not a count of bees and not a population model (InVEST's guide says so plainly). It does not know the size of habitat patches, the dynamics of the populations, pesticides, or how effective each species is as a pollinator, and its yields are relative to a full pollination. Compare scenarios run with the same tables, not absolute numbers between studies. For the share of habitat around sites at several radii, see Sampling, time series and redundancy; for one species' distribution, a species distribution model.

What you need#

• A land use map (lulc_2020.tif of Ribeirão Claro, from ribeirao-claro.zip), projected in metres, with integer codes; lulc_2040.tif for the scenario.
• The guild table: a row per guild (or species) with the columns species, alpha (the typical foraging distance, metres), relative_abundance, foraging_activity_<season>_index (0 to 1) for each season and nesting_suitability_<substrate>_index (0 to 1) for each substrate. Here two seasons, dry (April to September) and wet (October to March), and two substrates, ground and cavity:

The three guilds are ground-nesting solitary bees, stingless bees nesting in tree cavities, and large carpenter bees boring into dead wood. In guilds.csv the columns have their full names (nesting_suitability_ground_index, foraging_activity_wet_index...).

• The biophysical table: a row per land use with lucode, nesting_<substrate>_availability_index for each substrate and floral_resources_<season>_index for each season, all from 0 to 1. Every class of the map needs a row. Coffee flowers after the first rains, eucalyptus in the dry season; pasture offers bare ground to nest in, forest the cavities:

• Farms (optional): polygons with the columns crop_type, half_sat, season (the season the crop is pollinated, named as in the tables), p_dep (the share of the yield that depends on pollinators, \(\nu\)), p_managed (the share of pollination given by managed bees, \(mp\)), fr_<season> and n_<substrate> (the farm's own flowers and nesting, which replace the land use's under the polygon). Here sixteen coffee farms, the patches of cropland of at least 13.5 ha that are cropland in 2020 and in 2040, all with half_sat 0.05, season wet, p_dep 0.25, and the cropland's flowers and nesting; farms 3 and 9 rent hives of honey bees (p_managed 0.3), the others have none.

Enlarge

The guild table, the biophysical table and the farms are in pollination-inputs.zip: unzip it beside the files of Ribeirão Claro. Every value is illustrative. In a real study the guilds and their foraging distances come from surveys of the bees visiting the crop and from the literature (Greenleaf et al. (2007) related the foraging distances of 62 bee species to their body size: larger bees fly farther), and the nesting and floral values from field surveys or experts, documented for every class.

Build the model#

Step 1: the land use and the guilds#

In the model file the guilds read as a list:

format = "terraflux/1"

[[nodes]]
id = "lulc"
op = "raster.read"
params = { path = "lulc_2020.tif" }
note = "Land use, 30 m, UTM 23S."

[[nodes]]
id = "biophysical"
op = "table.read"
params = { path = "biophysical.csv" }
note = "Nesting and flowers by class."

[[nodes]]
id = "pollination"
op = "services.pollination"
inputs = { lulc = "lulc.raster", biophysical = "biophysical.table" }
params = { guilds = [
  { species = "ground", nesting_suitability_ground_index = 1, nesting_suitability_cavity_index = 0, foraging_activity_dry_index = 0.5, foraging_activity_wet_index = 1, alpha = 400, relative_abundance = 0.4 },
  { species = "stingless", nesting_suitability_ground_index = 0, nesting_suitability_cavity_index = 1, foraging_activity_dry_index = 1, foraging_activity_wet_index = 1, alpha = 700, relative_abundance = 0.4 },
  { species = "carpenter", nesting_suitability_ground_index = 0, nesting_suitability_cavity_index = 1, foraging_activity_dry_index = 0.6, foraging_activity_wet_index = 1, alpha = 1000, relative_abundance = 0.2 },
] }

guilds = "guilds.csv" (InVEST's guild table) gives the same maps (checked for this tutorial). Like the path of raster.read, the file name is taken from the model's folder. The seasons and the substrates are read from the names of the columns, and the two tables must name the same ones: a biophysical table with a season the guild table lacks stops the node and says which.

Run it (Ctrl+Enter on the node). Its outputs, each a stack (several maps on one grid, each layer named):

• supply: the pollinator supply \(PS\) of each guild, its layers named as the guilds;
• abundance: the abundance \(PA\) of each guild in each season, its layers named <species>_<season> (ground_dry, ground_wet, ...);
• total_abundance: the abundance of all guilds, \(PAT\), by season, its layers named as the seasons.

Step 2: the farms#

Read farms.gpkg with vector.read and connect its vector to the input farms. The farms may be in another coordinate system: the node reprojects them to the land use's grid, as InVEST does. From now on the guilds come from the file.

format = "terraflux/1"

[[nodes]]
id = "lulc"
op = "raster.read"
params = { path = "lulc_2020.tif" }

[[nodes]]
id = "biophysical"
op = "table.read"
params = { path = "biophysical.csv" }

[[nodes]]
id = "farms"
op = "vector.read"
params = { path = "farms.gpkg" }
note = "Coffee: half_sat, season, p_dep, p_managed, fr_, n_."

[[nodes]]
id = "pollination"
op = "services.pollination"
inputs = { lulc = "lulc.raster", biophysical = "biophysical.table", farms = "farms.vector" }
params = { guilds = "guilds.csv" }

With farms the node also gives:

• farm_pollinators: \(FP\) on the farms' cells, each farm in its own season (no value elsewhere);
• total_yield and wild_yield: \(PYT\) and \(PYW\) on the farms' cells;
• farm_table: a row per farm, in the order of the polygons, with farm (1, 2, ...), p_abund (the mean abundance in the farm's season), y_tot (\(Y_T\)), pdep_y_w (\(\overline{PYW}\)) and y_wild (\(Y_W\)).

Step 3: 2020 and 2040 side by side#

For the scenario, add a second land use reader (lulc_2040.tif) and a second pollination node with the same table, guilds and farms, and write both farm tables with table.write:

format = "terraflux/1"

[[nodes]]
id = "lulc_2020"
op = "raster.read"
params = { path = "lulc_2020.tif" }

[[nodes]]
id = "lulc_2040"
op = "raster.read"
params = { path = "lulc_2040.tif" }

[[nodes]]
id = "biophysical"
op = "table.read"
params = { path = "biophysical.csv" }

[[nodes]]
id = "farms"
op = "vector.read"
params = { path = "farms.gpkg" }

[[nodes]]
id = "pollination_2020"
op = "services.pollination"
inputs = { lulc = "lulc_2020.raster", biophysical = "biophysical.table", farms = "farms.vector" }
params = { guilds = "guilds.csv" }

[[nodes]]
id = "pollination_2040"
op = "services.pollination"
inputs = { lulc = "lulc_2040.raster", biophysical = "biophysical.table", farms = "farms.vector" }
params = { guilds = "guilds.csv" }

[[nodes]]
id = "table_2020"
op = "table.write"
inputs = { table = "pollination_2020.farm_table" }
params = { path = "farms_2020.csv" }

[[nodes]]
id = "table_2040"
op = "table.write"
inputs = { table = "pollination_2040.farm_table" }
params = { path = "farms_2040.csv" }

From the command line, with the model saved as pollination.toml in the folder of the data:

terraflux run ./pollination.toml --out results

How close is this to InVEST? Terraflux's tests compare the operator, map by map, with landschaft's version of InVEST 3.20's model on five cases (several guilds, farms, a wide kernel...): every cell holds the same 32-bit number, and the farm table agrees to within \(10^{-12}\). landschaft's own tests compare those cases with InVEST's outputs, to the bit on Windows.

Read the results#

Where the bees come from#

The supply is the map of sources: where the nests are, weighted by the flowers within their reach. Its maximum is set by the relative abundance: the carpenter bees, with 0.2, reach at most 0.10, the others about 0.2.

Where they forage#

A row of the map shows the two steps at work:

Enlarge

For the crop, the season that matters is the one it flowers in: the wet season's map over the coffee. Its mean is 0.174 on the coffee, against 0.215 in the forest and 0.118 in pasture. How does it change with the distance to the forest?

This is the pattern Ricketts et al. (2008) found in the field, with a twist the model makes explicit: it holds for the guilds that nest in the forest, not for those that nest in the open. The gradient here is gentle because no coffee of Ribeirão Claro is more than 600 m from a forest (half of it is within 150 m).

The farms#

Enlarge

Three things to read in this table. First, the abundance varies more than the yield: p_abund goes from 0.145 to 0.189 (30%), y_wild from 0.191 to 0.204 (7%), because with \(h = 0.05\) every farm is on the flat part of the curve. Second, the distance to the nearest forest is not the whole story: farm 4, 160 m from the forest on average, has as many bees as farm 6 at 51 m (0.188 and 0.189), and more than farm 8 at 142 m (0.176); what counts is how much forest and pasture lies within the bees' range, not only the nearest edge. Third, managed bees mask the wild ones: on farms 3 and 9 the hives fill the yield whatever the wild bees do, so their y_wild does not tell what the wild bees could do.

Settings that change the answer#

The foraging distance#

\(\alpha\) is the number most worth testing: published foraging distances vary widely, and it sets how far the forest's bees reach into the fields. Move the stingless bees' \(\alpha\) from 100 m to 1,400 m, the other guilds as before:

Enlarge

PlayForaging distance of the stingless bees, alpha (m)

A longer flight distance does not create bees: it moves them. The coffee next to the forest gets slightly fewer as \(\alpha\) grows past 200 m, the coffee far from it many more. If your question is which farms lose when a forest goes, the answer depends on \(\alpha\); run the model with the range of distances the literature gives for your guilds.

The half saturation#

The farms' yields depend on \(h\) much more than on anything else, because it sets where on the curve the farms sit. The same farms with four values of half_sat:

InVEST's guide suggests 0.05 as a default, with a range of 0.05 to 0.1, and calls it a parameter to tune after a first run. Like the half-saturation constant of habitat quality, it does not change the ranking of the farms, but it changes every yield you report: choose it once, say how, and keep it for every scenario.

The scenario of 2040#

The scenario clears forest on the gentle edges, turns pasture into coffee, plants a new eucalyptus block in the north-east, grows the town and restores forest along the rivers. The farms stay as they are; their surroundings change.

2020

2040

Enlarge

Note what the mean distance of the coffee to the forest says: almost nothing (177 m against 176 m). The loss comes from the amount of forest within the bees' range, which a distance to the nearest edge does not measure and the model does.

Common pitfalls#

• Names with an underscore. Seasons and substrates are read from the column names (foraging_activity_wet_index), so their names cannot contain _: wet_season is not read as a season. Use wet, dry, ground, cavity.
• Tables that do not match. The guild table, the biophysical table and the farms must name the same seasons and substrates; a farm's season must be one of them. The node stops with a message that names the mismatch.
• A missing class. Every class of the land use needs a row of the biophysical table; a class without one stops the node.
• Geographic coordinates. \(\alpha\) is in metres and is divided by the cell size; a map in degrees makes it thousands of cells. Reproject the land use first (see Coordinate systems).
• The edge of the map. Near the edge, the flowers and the nests beyond it are not counted: the mean is over the part of the flight range inside the map. Include a margin of land around the farms of at least the largest \(\alpha\), better twice.
• Small resources the map misses. Field margins, hedgerows and gardens are often smaller than a cell and absent from a land use map, and Lonsdorf et al. (2009) found that they can decide the result. Give the farms their own fr_ and n_ values when you know them.
• Abundance read as bees. 0.17 is not 0.17 bees: the indices rank places and scenarios, they do not count insects.
• Yields read as harvests. \(Y_T = 0.95\) is 95% of the yield the farm would have with full pollination, not a number of bags of coffee.

Going further#

• Ecosystem service models: the ideas, part 1, for the landscape and the other services.
• Habitat quality and rarity for the forest the bees nest in, seen as habitat.
• Carbon storage and sequestration for the same scenario in tons of carbon.
• Several benefits and trade-offs to weigh pollination with other benefits in a restoration plan.
• Sampling, time series and redundancy for the landscape around field sites at several radii, to relate observed bees to the land around them.
• The operator's reference: services.pollination.

References#

• Greenleaf, S. S., Williams, N. M., Winfree, R., & Kremen, C. (2007). Bee foraging ranges and their relationship to body size. Oecologia, 153(3), 589–596. https://doi.org/10.1007/s00442-007-0752-9
• Klein, A.-M., Vaissière, B. E., Cane, J. H., Steffan-Dewenter, I., Cunningham, S. A., Kremen, C., & Tscharntke, T. (2007). Importance of pollinators in changing landscapes for world crops. Proceedings of the Royal Society B: Biological Sciences, 274(1608), 303–313. https://doi.org/10.1098/rspb.2006.3721
• Kremen, C., Williams, N. M., Aizen, M. A., Gemmill-Herren, B., LeBuhn, G., Minckley, R., Packer, L., Potts, S. G., Roulston, T., Steffan-Dewenter, I., Vázquez, D. P., Winfree, R., Adams, L., Crone, E. E., Greenleaf, S. S., Keitt, T. H., Klein, A.-M., Regetz, J., & Ricketts, T. H. (2007). Pollination and other ecosystem services produced by mobile organisms: a conceptual framework for the effects of land-use change. Ecology Letters, 10(4), 299–314. https://doi.org/10.1111/j.1461-0248.2007.01018.x
• Lonsdorf, E., Kremen, C., Ricketts, T., Winfree, R., Williams, N., & Greenleaf, S. (2009). Modelling pollination services across agricultural landscapes. Annals of Botany, 103(9), 1589–1600. https://doi.org/10.1093/aob/mcp069
• Natural Capital Alliance. (2026). InVEST User's Guide: Crop pollination (pollinator abundance). https://storage.googleapis.com/releases.naturalcapitalproject.org/invest-userguide/latest/en/croppollination.html
• Ricketts, T. H., Daily, G. C., Ehrlich, P. R., & Michener, C. D. (2004). Economic value of tropical forest to coffee production. Proceedings of the National Academy of Sciences, 101(34), 12579–12582. https://doi.org/10.1073/pnas.0405147101
• Ricketts, T. H., Regetz, J., Steffan-Dewenter, I., Cunningham, S. A., Kremen, C., Bogdanski, A., Gemmill-Herren, B., Greenleaf, S. S., Klein, A. M., Mayfield, M. M., Morandin, L. A., Ochieng', A., & Viana, B. F. (2008). Landscape effects on crop pollination services: are there general patterns? Ecology Letters, 11(5), 499–515. https://doi.org/10.1111/j.1461-0248.2008.01157.x