Connecting earth observation anomalies to farmer surveys for monitoring impacts of agricultural drought on rainfed rice yields in Nigeria

Agricultural drought threatens rainfed rice production in Nigeria, where smallholder farmers depend on rainfall and have limited capacity to buffer climate shocks. While meteorological drought indices such as the SPI and the SPEI are widely used in national early warning systems, their ability to capture the impacts of droughts on rainfed rice yields at the smallholder field-level remains uncertain. This study evaluates the added value of earth observation (EO)-derived vegetation and soil moisture anomalies for monitoring and predicting drought impacts on rainfed rice yields in Nigeria. Satellite-based Normalized Difference Vegetation Index anomalies (NDVIA) and Soil Water Index anomalies (SWIA) were derived using a zonal clustering and thresholding approach and combined with farmer survey and yield data collected from 146 rainfed rice farmers across four major rice-growing states between 2019 and 2024. Multivariate regression models were used to assess the relationships between EO indicator anomalies and annual yield changes, and the effects of different zonal clustering and anomaly thresholds on anomaly calculation were evaluated. Results show that SPI and SPEI explain a substantial share of yield variability in some years, particularly when droughts coincide with sensitive phenological stages. However, EO-based anomaly indicators, especially SWIA (maximum improved R² = 0.25), provide complementary information and significantly improve yield predictions in years when meteorological indices alone perform poorly. The timing of anomalies relative to rice phenology was critical, with droughts during panicle initiation having the largest yield impacts. Integrating EO-based vegetation and soil moisture anomaly indicators with existing meteorological indices can contribute to the monitoring of agricultural droughts and improve the operational relevance of early warning systems for rainfed rice farmers in Nigeria.

1 Introduction

Droughts have serious impacts on farmers worldwide, resulting in yield losses that threaten farmer incomes and food security for local and global populations (Masih et al., 2014). Meteorological droughts are characterized by a water deficit caused by an imbalance between precipitation and evaporation. The resulting decrease in soil moisture reduces water availability for plants, causing agricultural droughts that damage crops and negatively affect yields and farmer incomes, especially for smallholder farmers without the means or infrastructure for irrigation (Liu et al., 2016). A large proportion of the population in Nigeria relies on subsistence farming (Onyeneke, 2021) and is vulnerable to the impacts of changing climate conditions on rainfall patterns. Rainfall patterns in Nigeria are projected to change drastically, with increasing drought intensity despite increasing rainfall, driven in large part by greater variability during the rainy season, which brings more extreme conditions (Ogunrinde et al., 2019, 2021). Climate modelling has shown an increasing frequency and duration of droughts in Nigeria, threatening the availability of water resources for agricultural use (Ogunrinde et al., 2022).

Rice is a staple crop in Nigeria, with consumer demand growing rapidly as domestic consumption patterns change. On the production side, rice has become a leading cereal crop in terms of value for farmers and is central to food security policy in the country (Gyimah-Brempong et al., 2016). Rainfed rice systems are crucial for rice production in Nigeria, with rainfed lowland ecologies accounting for 47 % of total production area, and rainfed upland ecologies representing over 30 % of rice production by land area (over 560 000 ha) (Alagbo et al., 2022). Despite large areas of rice production, Nigerian rice production accounts for 50 %–60 % of rice consumed in the country, with the remaining deficit made up by foreign rice imports (FAO in Nigeria, 2026; Onafeso et al., 2016; Udemezue, 2018). Rice production in Nigeria is dominated by smallholder farmers (approximately 90 %) (Chiaka et al., 2022), who cultivate small plots of land ranging from 0.2 to 3 ha. Rice farmers in the country remain vulnerable to extreme weather events (such as droughts), which can result in lower yields, stunted rice crops, and poor economic outcomes for farmers (Ajetomobi et al., 2011; Mbah et al., 2016; Nwalieji and Uzuegbunam, 2012; Ojo et al., 2022). Methods to monitor and predict drought impacts on rice yields are important at both the policy and farmer level to ensure food security and safeguard local rice supply against climate change.

Historically, several droughts have occurred in Nigeria, especially in the arid Northern states that include the Lake Chad region. Past drought events that significantly affected the population include the severe droughts of 1971–1973, 1983–1984, and 1991–1995 (Federal Ministry of Environment, 2018), which led to widespread crop failures, livestock deaths, depleted groundwater levels, and famine. These effects on agricultural production prompted the Federal Government of Nigeria to establish institutional frameworks and programs to minimize drought-related vulnerabilities. The Nigerian Meteorological Agency (NiMet) is responsible for analysing meteorological data, incorporating both satellite and weather station data, and delivering advice to farmers based on seasonal predictions, particularly for seasonal rainfall prediction and crop performance (Ofuoku and Obıazı, 2021). Nonetheless, current techniques are limited (Oloyede et al., 2023) and mainly rely on observations from sparse ground weather stations and long-term climate modelling (Akinyemi et al., 2020; Nigerian Meteorological Agency, 2025). AGRHYMET, a regional institution that monitors droughts in West Africa, regularly provides weather bulletins, crop monitoring analyses based on simulation models, and vegetation indices covering West Africa (including Nigeria) that track the growth and condition of vegetation (Traore et al., 2014). Earlier research has modelled and forecast rainfall patterns across West Africa, including Nigeria, incorporating both regional and global climate models (Olaniyan et al., 2022; Onyejuruwa et al., 2025). Long-term and seasonal predictions support farmers in seasonal planning; however, acute weather phenomena such as droughts and dry spells are difficult to predict, and farmers need early warning systems (EWS) to respond to such events. This is particularly relevant in West Africa, where rainfed agriculture is crucial for food security and is increasingly threatened by inconsistent rainfall patterns, droughts during crucial crop growth stages, and changes in the frequency and intensity of dry spells during the rainy season (Tefera et al., 2024). While many Earth Observation (EO)-based EWS exist for drought and vegetation condition at the global scale (Gommes et al., 2017; Rojas, 2021; Van Hoolst et al., 2016), they are often aimed at informing policy rather than providing farmers with clear warnings and guidelines during the growing season – even when tailored to the national scale (van Ginkel and Biradar, 2021; Vos et al., 2023). It is also unclear to what extent existing EWS are accessible to farmers or whether they are effective at predicting the impacts of drought events on farmers' harvest yields and income.

Currently, techniques used by NiMet to provide early warnings for farmers incorporate EO to monitor meteorological droughts. Meteorological drought indices are a common approach to understanding characteristics such as drought severity and are often linked to insurance schemes for smallholder farmers. The Standardized Precipitation Index (SPI) is one such index and has been used to model precipitation deficits in Nigeria, however researchers note that drought management responses remain reactive rather than proactive (Okpara et al., 2017) – a challenge for all systems that rely on monitoring alone. SPI is currently the primary drought index used by NiMet in quarterly bulletins analysing meteorological trends across Nigeria (Nigerian Meteorological Agency, 2025). SPI is also highlighted by the World Meteorological Organization as a useful indicator for drought monitoring (WMO and GWP, 2016). While the SPI is based only on precipitation data, the Standardized Precipitation and Evapotranspiration Index (SPEI) considers water balance deficit through also incorporating evapotranspiration from plants and soils. SPEI has been used in the past to assess climatic impact on rice yields in the Gambia, demonstrating that increasing temperatures and lower precipitations are correlated with significant decreases in rice yields (Jabbi et al., 2021). Notably, SPI and SPEI correspond well in humid regions, however in semi-arid and arid regions the signals diverge, with SPEI detecting droughts more frequently and earlier than SPI (Mwinjuma et al., 2026). This is particularly relevant when studying agricultural trends in countries such as Nigeria, which contain gradients of aridity along which the same crop (in this case rice) is grown. However, these indices ultimately detect meteorological droughts, which depend on precipitation and evaporation patterns, while agricultural droughts link the impacts of meteorological droughts to adverse impacts on crop productivity due to a deficit in root zone soil moisture (De Vos et al., 2025; Hazaymeh et al., 2016). The impacts of meteorological droughts on crops can also be varied, especially in complex and heterogeneous environments found in African agricultural landscapes (Guo et al., 2023).

Assessing the impact of agricultural droughts involves monitoring the evolution of crop conditions over time in response to meteorological droughts. However, the vegetation response to drought conditions is not always easy to predict from EO data, and may follow complex, localized patterns depending on the type of vegetation and the region (Henchiri et al., 2020). EO can be used to monitor vegetation dynamics using satellite-derived vegetation indices. The Normalized Difference Vegetation Index (NDVI) is a common vegetation index calculated from the red and near-infrared bands of multispectral satellite images, and NDVI anomalies (NDVIA) can be used to assess the degree of plant stress due to a variety of factors, including drought. NDVI has been employed in the past to develop an open-source EWS for drought monitoring in Nigeria (Adedeji et al., 2020). One advantage of NDVI as a drought indicator is the long historical timeseries available from past and present EO missions, allowing for long-term statistical analysis and robust anomaly calculation. Indices based on the analysis of NDVIA over a historical time period, or during an individual growing season, have been implemented as indicators of vegetation condition and to assess drought risks across large regions (Kogan, 1990, 1995; Traore et al., 2014). Changes in NDVI or NDVIA during multiple stages of rice development are associated with impacts on grain yield (Ali et al., 2014; Phyu et al., 2020). However, it remains unclear how changes in vegetation indices translate into actionable information for farmers about potential drought impacts on their rice crops in rainfed rice systems. This is particularly the case for NDVIA-based indicators, which inherently involve methodological choices when defining an anomaly threshold relative to a spatiotemporal baseline.

While NDVI can be used to monitor vegetation conditions, soil moisture can also serve as an effective proxy for agricultural drought. The Soil Water Index (SWI) estimates the soil moisture profile from EO data and can be used for representing the root-zone soil moisture available at different depths (Bauer-Marshallinger, 2018; Paulik et al., 2014). A maximum depth of 60 cm has been used in the past to describe the rooting depth of common rainfed rice varieties grown in Africa (Kijoji et al., 2014). Soil moisture anomalies calculated from EO data have been used in past studies (sometimes in conjunction with NDVI) to identify drought events and monitor drought impacts on vegetation and agriculture (De Vos et al., 2025; Tian et al., 2022; Zribi et al., 2021). Due to the freshwater demands of rice crops, soil moisture is a key variable to monitor for understanding the yield impacts of anomalies in soil water content (Zampieri et al., 2018). Soil moisture time series have been linked to rice yields from national-scale farmer surveys in Madagascar (Rigden et al., 2022), however a link with localised surveys with smallholder farmers is missing. In particular, an operational gap exists between the EO indicator-based monitoring of droughts and how they affect rainfed rice yields for smallholder farmers most affected by such droughts. Our research aims to bridge this gap by contributing to a better understanding of how EO indicator-based drought analysis can be tailored to increase relevance for smallholder farmers and ultimately improve the operational effectiveness of EWS systems employed by NiMet.

The objective of this research is, therefore, to assess the added value of EO indicator anomalies for predicting the impacts of agricultural droughts on rice yields in Nigeria in addition to the indicators already monitored by NiMet. We use a statistical approach to calculate pixel-level anomaly values from existing EO indicators to monitor vegetation condition and soil moisture (Adefisan and Abatan, 2015). Furthermore, we assess the impacts of key parameters in calibrating EO anomaly indicators, namely the selection of the optimal number of agricultural zones and the upper threshold for identifying negative drought impacts. Finally, we compare this approach with one based solely on meteorological indices to assess whether EO anomaly indicators improve monitoring of drought impacts on rainfed rice production. We validate EO anomaly indicator-based drought predictions using farmer survey data, a comprehensive analysis of zonal calibration and threshold sensitivity for calculating EO indicator anomalies, and a comparison of EO-anomaly indicators with traditional drought indices for predicting drought impacts on rice yields. We use the results of these analyses to understand the extent to which EO can be used as an input to improve EWS for predicting drought impacts on rainfed rice farmers in Nigeria, with the intention of informing policy and supporting the development of improved EWS systems for farmers in the region.

2 Methods

To evaluate how well traditional drought indices and the proposed NDVIA approach monitor drought impacts on Nigeria's rainfed rice systems, we combine satellite- and climate-derived indicators with plot-level survey, GPS, and yield data from major rice-growing states (Fig. 1). We then compare index signals with reported drought impacts and model rice yields using these inputs to assess early-warning potential and quantify drought–yield relationships.

Figure 1 Overview of methodology followed in this study.

2.1 Study area

With two main seasons (dry and wet), climate conditions are generally favourable for agriculture in Nigeria. The Tropical Continental airmass drives the dusty Harmattan winds from the Sahara Desert during the dry season, while the Tropical Maritime airmass from the Atlantic Ocean has a significant impact during the rainy season (Eludoyin et al., 2014). Total annual rainfall ranges from 1500 to 4000 mm in the South and 300 to 1000 mm in the North, with droughts occurring mostly at the start of the wet season due to the impenetrability of the moisture laden airmasses deep into the hinterlands (Federal Ministry of Environment, 2018).

Figure 2 Location within Nigeria of the states of interest and local government areas (LGAs) surveyed in this study. Rainfed rice production percentages obtained from (Mba et al., 2021; Omoare and Oyediran, 2020) and average rice yields obtained from (National Bureau of Statistics, 2019). Local ecological zones adapted from (Adenle et al., 2020; States in Nigeria by Ecological Zone, 2026; FORMECU, 1998). Additional details in Table A1.

Nigeria consists of a gradient of agroecological zones that run from south to north from the Atlantic coast to the arid Sahel plains (Fig. 2) (Adenle et al., 2020; FORMECU, 1998). The humid tropical forest zone of the South, with its bimodal rainfall pattern and rich sandy loam-clay soils, supports rainfed rice production, while the Sudan-Sahel zones of the arid North rely solely on intensive irrigation for rice cultivation due to the below-average rainfall and sandy soil with low water-holding capacity (Nigeria – Global yield gap atlas, 2026). Rainfed rice production in Nigeria is categorized into two major systems: rainfed upland and rainfed lowland (referred to as fadama in northern Nigeria – rice cultivation in low-lying valleys and floodplains). The upland system which accounts for 30 % to 35 % of the total rice cultivated area in Nigeria (Erenstein et al., 2003) is characterized by low fertility, well-drained soils that are entirely dependent on rainfall for moisture, leading to poor yields in the absence of secondary water sources. The lowland rice system, the dominant system in the study area, covers nearly 50 % of the total rice area (Erenstein et al., 2003). Yields in the lowland are substantially higher because the rice fields are constantly flooded or waterlogged during the rains, with water supplied from multiple sources, including direct rainfall.

The rice calendar of Nigeria is complex and varies across the different rice production systems. Generally, in the South, land preparations commence in March to April with the first rains, with sowing or transplanting occurring between April and June, and subsequent harvest from September-October, depending on variety. With a shorter wet season, land preparations in the North occur between May and June, with sowing in June and July, followed by harvests between October and November.

Fieldwork was conducted in March and April 2025 with the purpose of administering surveys to rainfed rice farmers in the target states, as well as collecting GPS coordinates of their rice plots. In total, 10 local government areas across four states were visited, where farmers from rainfed rice-farming communities responded to predefined survey questions through face-to-face interviews. After reviewing the surveys for completion and accuracy, farmers went with trained enumerators to identify their rice plots and collect GPS coordinates. In total, surveys were administered to 226 rainfed rice farmers, with 154 rice plots identified by GPS location. Information collected in the surveys included data on predominant soil type, rice varieties, timing of transplanting and harvesting, and the yield impacts of past droughts experienced over the years 2019–2024. Yield data were collected according to farmers' recollections of annual harvests, converted to metric tonnes, and tonnes per hectare were calculated per farmer plot using GPS-delineated polygons. The field data were used to assess and validate the drought analysis approaches implemented in this study.

2.2 Satellite data

Following the approach of De Vos et al. (2025) NDVI time series were extracted using openEO from two harmonized Copernicus Land Monitoring Service products for Nigeria. The CGLS NDVI v3 product, generated from SPOT/VEGETATION and PROBA-V sensors at a 1 km spatial resolution (Smets et al., 2020), was obtained in 10 d intervals (dekads) covering dates from 1999 until 2020, and from 2020 until 2024 the CGLS NDVI v2 was obtained at 300 m spatial resolution (Swinnen et al., 2022), generated from the Sentinel-3 OLCI sensor was used with the same time intervals and resampled to 1 km spatial resolution using nearest neighbour resampling in openEO in order to match the spatial resolution of the v3 product. Invalid pixels (e.g., due to cloud coverage) were filtered out. SWI data was available through the CGLS SWI v4 product, generated from the H29 Metop ASCAT 12.5 km SSM Near Real Time product, available at a 12.5 km spatial resolution at 10 d intervals (Raml and Bauer-Marshallinger, 2025). The lower resolution data was upscaled using bilinear interpolation to match the 1 km resolution of the NDVI products.

Next, in addition to per-pixel descriptive statistics (long-term average, standard deviation, 10th, 50th, and 90th percentiles), a harmonic analysis was conducted to model the multiannual dynamics of seasonal vegetation changes. This harmonic analysis was conducted because vegetation dynamics can exhibit periodicity that varies with vegetation type, agroecological zoning, and management practices. By using sinusoidal waves to describe those patterns, parameters such as the amplitude and phase of the sinusoidal waves can be used to differentiate between different crop types or systems (Jung and Chang, 2015). A least-squares approach was used to find the coefficients of sine and cosine harmonic waves for the first three harmonics (Eq. 1) to represent the maximum number of potential growing seasons in the study area (with two potential rice seasons and an additional fallow/vegetable season) and to align the wavelength with the duration of a rice growing season (approximately 100 d). Pixel-level sine/cosine coefficients (a_f, b_f), amplitude (A_f) (Eq. 2), and phase (ϕ_f) (Eq. 3) were also used as inputs for the zonal clustering approach.

$$\begin{array}{}\text{(1)}& \mathrm{NDVI}\left(t\right)=a_{\mathrm{0}}+\sum _{f=\mathrm{1}}^{\mathrm{3}}\left(a_f\cos \left(\mathrm{2}\mathit{\pi }ft\right)+b_f\sin \left(\mathrm{2}\mathit{\pi }ft\right)\right)\end{array}$$

$$\begin{array}{}\text{(2)}& {\displaystyle }{A}_f=\sqrt{a_f^{\mathrm{2}}+b_f^{\mathrm{2}}}\text{(3)}& {\displaystyle }{\mathit{\varphi }}_f=\mathrm{atan}\mathrm{2}\left(-b_f,a_f\right)=\mathrm{2}\arctan \left({\displaystyle \frac{\sqrt{a_f^{\mathrm{2}}+b_f^{\mathrm{2}}}-a_f}{b_f}}\right)\end{array}$$

Due to the different agroecological zones in Nigeria and considering the differences in agricultural management practices, access to water, and crops grown across different communities, a spatial zoning approach was used to identify drought patterns as they occur within individual regions that share similar characteristics (Adefisan and Abatan, 2015; Page et al., 2017; Perera et al., 2022). To this end, zonal clustering based on the KMeans algorithm was chosen as an appropriate approach to separate pixels into separate zones based on the similarity of their NDVI values, harmonics, and location (latitude and longitude included as features). The number of zones used in clustering was tested to understand the impacts of zoning on anomaly indicator values. This approach has been used in the past, especially in the crop insurance industry, where zoning is important for accurately modelling climate risk based on agricultural production patterns and satellite-derived environmental data (Kirchner et al., 2025). The zone delineations calculated from historical NDVI data were also used to delineate zones for processing SWI data (Bauer-Marshallinger, 2018; Paulik et al., 2014) in the following steps. The SWI band “SWI_040” was used, representing a T latency factor of 40, which roughly represents the water present in soils up to a depth of 60 cm (Kijoji et al., 2014; Raml and Bauer-Marshallinger, 2025).

2.3 Zonal adjustment and anomaly calculations

Reflecting differences in environmental and management variables across agroecosystems in Nigeria, zonal dekadal statistics were calculated for both EO indicators (SWI and NDVI), summarising the two historical datasets per dekad and per pixel for all available years, including mean, standard deviation, and percentile values. The dekad pixel mean was adjusted by the respective zonal median to obtain an adjustment factor, which was then applied to the EO indicator data at the pixel scale to reduce outlier values in each zone. However, the total number of zones had a large impact on both the calculation of EO indicator anomaly values as well as the final results. Therefore, we tested a range of zones (5–120) in 10-zone increments to assess the impact of different zone distributions (and as a result, different zonal statistics) on yield model results (Fig. 3).

Figure 3 Different spatial zones over study area showing increasing complexity of zone distribution with (a) 10 zones, (b) 60 zones, and (c) 120 zones. The pixel zone clusters were used to calculate zonal adjustment factors for both NDVI and SWI anomalies.

The dekadal anomaly calculation is similar to the Vegetation Condition Index (VCI) proposed by (Kogan, 1990), however rather than identifying minimum and maximum EO indicator values over the historical time period, we calculate pixel-level statistical thresholds at low historical EO indicator values. Equation (4) is used to normalize EO indicator values to upper and lower thresholds to values between 0 and 100. Zero values indicate no EO indicator anomalies, while non-zero values scale linearly to 100, which indicates an EO indicator anomaly below the lower threshold value.

$$\begin{array}{}\text{(4)}& {\mathrm{EOIA}}_d={\displaystyle \frac{{\mathrm{EOI}}_d-{\mathrm{EOI}}_{\mathrm{upper}}}{{\mathrm{EOI}}_{\mathrm{lower}}-{\mathrm{EOI}}_{\mathrm{upper}}}}\end{array}$$

Where for a single pixel, the EO indicator anomaly (EOIA) value for a dekad d is determined by a ratio of differences between the EO indicator value for that dekad (EOI_d) and the lower (5th percentile) and upper (e.g. zonal 15th percentile) threshold values (EOI_lower, EOI_upper). These threshold values were used by De Vos et al. (2025), however we also tested a range of EO indicator upper threshold values up to 50 % as higher upper threshold values by design will result in more pixels identified as drought, ultimately influencing model outcomes. This range of values was tested to identify the upper threshold value leading to the best yield model results, while the lower threshold was maintained at 5 %. Using the zone distributions generated from NDVI data, we applied Eq. (4) to calculate anomaly values for both SWI and NDVI.

2.4 Multivariate regression analysis

To understand how both traditional meteorological indices and EO indicator anomaly variables could be used to monitor drought impacts on farmers' rice yields, available meteorological and vegetation indices, together with field data, were used as input variables to fit multivariate regression models. Field data (soil type, rice variety) were converted from categorical to dummy variables for use as inputs. For input variables derived from meteorological indices and remote sensing data, summary statistics (median, minimum, maximum, standard deviation) were used to summarize dekad values for the entirety of the rice growing season. SPI/SPEI variables were first used to fit regression models with all aggregation time windows (1-, 3-, 6-, and 12-month periods), but models were also fit to the index of each time window to understand whether the preceding precipitation and evapotranspiration trends were important for predicting yield changes. This approach was also aimed at understanding whether different time windows might be more correlated to yield changes in different years due to different drought conditions and patterns. Additionally, for EO indicator anomaly variables, the number of dekads with non-zero anomaly values during the growing season was used as an input variable. Inputs were grouped by category to assess the additional benefits of each set of input variables in improving model accuracy (Table 1). Annual changes in rice yields were calculated per plot using the annual yields reported by farmers in the field surveys. Yield data were cleaned by removing null yields and yield change values outside the 5th/95th percentile range for each year to prevent outliers from affecting model accuracy.

Table 1 Input data for multivariate regression analysis.

To assess whether NDVI/SWI patterns across distinct parts of the rice growth cycle were correlated with rice yield changes, we also calculated anomaly variable statistics over separate phenological stages, rather than over the entire growing period. Due to the short duration of the panicle initiation growth stage (approximately three dekads), the rice growing season was only divided into two phenological periods (vegetative and maturity), separated at the beginning of the panicle initiation period. Phenological summary statistics (median, minimum, maximum, standard deviation) summarized dekad values over each phenological period. These summary statistics were also used as model inputs and model results were compared with those obtained with NDVIA/SWIA summary statistics calculated over the entire growing season. Multivariate regression models were trained, and accuracy metrics (adjusted R² score [adj-R²], RMSE, model p-value) were iteratively calculated for EO indicator anomaly values calculated in a range of 5 to 120 zones, and 15 % to 50 % upper threshold values. The other input variables remained the same across all regression models. Regression models were fit to yield data for 2020–2024, and model adj-R² scores were used to evaluate impacts of zoning and upper threshold values on model results.

3 Results

3.1 Field survey results

Droughts in rainfed rice plots were mostly reported by farmers at the height of the rice growing season, between June and September. From the 5 years surveyed, 2021 stands out as a year with less farmer-reported droughts throughout all growing season months. In contrast, there were many more droughts reported by farmers in 2024 than in all previous years – particularly in the months of July and August, there were more than double the number of plots with drought reported than in any of the previous years (Fig. 4). The increased number of farmers reporting drought impacts on their rainfed rice plots in 2024 underscores the importance of improving climate-resilient agricultural practices, scaling up local irrigation infrastructure, and expanding EWS and localized climate advisories to help farmers adapt in real time to changing climate conditions. Due to differences in transplant time and rice varieties used by farmers, varieties reach different stages of their development cycle at different times in the growing season. Though farmer responses were not divided on an annual basis, 62 % of the surveyed farmers had experienced drought stress during the panicle initiation stage of their rice crops, suggesting large negative impacts of drought on their final yields (Fig. 4). Cleaned rice yield changes generally followed normal distributions centred around no annual change (Fig. A1). However, 2020 and 2024 yield changes were skewed toward negative yield changes, indicating most farmers' yields reduced those years, while 2021 yield changes skewed positive, likely due to the yield losses the year before. 2024 also had the highest variability in yield changes, with some reported yield losses even below yield values that would be possible in rainfed rice fields. Yield change distribution was skewed slightly negative in 2022 as well, indicating a cluster of farmers that experienced yield losses, while 2023 yield changes were evenly distributed between yield losses and gains.

Figure 4 Months with drought(s) reported by farmers during field surveys (top), and rice crop development stages during which farmers experienced drought stress on their crops (bottom).

3.2 Anomaly calculation parameters

The number of dekads with non-zero NDVIA and SWIA values during the entire rice growing season at each surveyed plot location increased with higher upper threshold values. This was expected, as an increased upper threshold decreases the sensitivity of anomaly calculations, thereby identifying a broader range of EO indicator values as anomalies. Notably, an increased number of zones decreased the total number of non-zero NDVIA values, while an increased number of zones led to slightly more non-zero SWIA values, particularly at larger upper threshold values (Fig. 5). Initial principal component analyses showed that NDVI statistical values were distributed in a gradient across the first two principal components, however a silhouette score analysis identified slight elbow patterns around the 20–30 zone and 50–60 zone marks (Fig. A2). Furthermore, there were differences in the number of non-zero NDVIA and SWIA values per year, with the lowest number of non-zero NDVIA and SWIA values in 2021 and 2023, respectively. Both indices identified 2020 as the year with the largest number of non-zero anomaly values. This does not necessarily align with farmers' self-reported annual drought impacts, which were lowest in 2021 and reached a peak in 2024. However, as we count all dekads during a growing season that had non-zero NDVIA/SWIA values, some of the counted drought dekads did occur during parts of the growing season that would have had limited impact on rice yields.

Figure 5 Number of dekads with non-zero NDVIA (top) and SWIA (bottom) anomaly values for tested combinations of zones and upper threshold values.

Evaluating the spatial distribution of NDVIA with different combinations of zones and upper threshold values, clear patterns emerge as well during identified drought moments (Fig. 6). At a very low zone number, pixel values of the very arid north of Nigeria are included in summary statistics for some of the northern zones, which means that desert-like seasonal patterns influence the calculation of both vegetation seasonality and anomaly values. By mixing desert and vegetation pixels in the same zonal clusters, the algorithm reduces the specificity of the anomaly predictions, leading to higher and more generalised NDVIA values. In contrast, once zonal complexity is more accurately captured (such as above 40 zones in Fig. 6), there is not much difference in anomaly values with an increasing number of zones. Increasing upper threshold values increase the number of non-zero NDVIA values, and increase the NDVIA values themselves, regardless of the number of zones. For SWIA results, similar patterns emerge, with increasing upper threshold values increasing the extent of pixels with non-zero SWIA values. Increasing the zone number also identifies more areas with soil moisture anomalies, particularly in the more arid north of the country (Fig. A3). From the spatial distribution of anomaly pixels, there is a clear minimum of zones necessary for more accurate drought predictions, however an optimal number of zones is difficult to identify.

Figure 6 NDVIA values for different combinations of zones and upper threshold values capturing anomalies in vegetation greenness during the dekad 11 June 2020.

3.3 Multivariate regression analysis

All available variables as described in Table 1 were used as inputs to multivariate regression models to predict drought impacts on annual yield changes per rice plot. Year-over-year changes in rice yields (as reported in farmer surveys) were used as the dependent variable, while assuming that fields and management practices remained largely similar. Adjusted R² scores were used to account for the increasing number of input variables in the different models (Table 2). Of the farmer survey variables, neither soil type nor rice variety had any correlation with annual yield changes, and these survey variables were removed from subsequent models. The meteorological indices had the highest significant adj-R² scores and were statistically significant for the 2024, 2020, and 2021 models. Both SPI and SPEI had notably higher scores than other variables when used as the sole model inputs, however they had intermediate scores when both indices were used as model inputs. The models incorporating all variables had lower adj-R² scores relative to the meteorological variable models, suggesting that both collinearity and the impacts of the non-significant survey variables were detrimental to model accuracy. While the “All years” models had lower adj-R² scores, there was still some predictive power for the multi-year models incorporating only meteorological variables. Notably, RMSE values differed between years, with 2024 models showing higher RMSE values, larger than typical rice yield values, while RMSE values for 2020–2023 were low in comparison.

Table 2 Adjusted R² scores and RMSE (in brackets) for annual and all year regression models trained with individual data categories, and all categories except anomaly variables (NDVIA, SWIA). Values in bold represent statistically significant model adjusted R² scores (p > 0.05).

Table 3 Adjusted R² scores and RMSE (in brackets) for annual and all year regression models trained with individual aggregations of SPI and SPEI input variables. Bolded values represent statistically significant model adjusted R² scores (p < 0.05).

As SPI and SPEI are calculated over different time windows (of the preceding 1-, 3-, 6-, and 12-month periods), we also compared the input variables from each index as separate inputs to multivariate regression models (Table 3). All SPI and SPEI aggregations had significant correlations with yield change values for the years 2020, 2021, and 2024, as well as for the “All years” models. Only SPI6 had a significant correlation with 2022 yield changes. SPEI generally performed better than SPI, due to the inclusion of evapotranspiration, which is a proxy for crop water demand. There was little difference in performance between the short-term and long-term indices, with the majority of the difference attributable to models trained on yield data from different years. All input variables' adj-R² scores were both very low and not significant for the yield changes observed in 2023. The low adj-R² scores of the “All years” model suggest there is high interannual variability, leading to model instability across years, preventing model generalization to other years. Furthermore, high RMSE values for 2024 suggest that while the models for that year have higher R² scores, the spread of their predicted yield changes is very large, beyond the expected rice yield values for the region. This would affect the applicability of such a model for actual yield predictions and may have been due to large yield change values in the 2024 training data.

Anomaly variables calculated with incremental combinations of zone numbers and upper threshold values were used as inputs together with SPI and SPEI to multivariate regression models. Model adj-R² scores, RMSE, and p-values were calculated for each model, and the improvement in adj-R² compared to the highest yearly adj-R² model as described in Table 2, was visualized (Fig. 7).

Figure 7 Improvement in model adjusted R² score relative to best annual model trained without anomaly variables. Improvement in model adjusted R² shown for NDVIA (top), SWIA (second from top), combined (third from top), and only NDVIA/SWIA (bottom) for models with different combinations of zones and upper threshold values.

The addition of EO indicator anomaly variables improves the adj-R² scores of the 2020 model more than for the other years. SWIA data was associated with increased the model adj-R² more than NDVIA data, but the model incorporating both anomaly variables and meteorological variables (combined model) had the highest performance. Similarly, SWIA data improved model adj-R² scores for the 2022 model, as well as slightly improving the 2023 scores. Notably, the models incorporating SWIA performed worse than those using only meteorological variables in 2021 and 2024. NDVIA data did not improve performance for the models of any other year and performed much more poorly than the meteorological-only model in 2021. For the combined model incorporating meteorological and anomaly variables, we observed the greatest improvement in adj-R² scores for 2020, suggesting that both NDVIA and SWIA provided additional, complementary information that improved the regression results. Some of the 2022 models also improved adj-R² scores more than models with a single-anomaly input. When meteorological variables were removed and only the EO indicator anomalies were used, model performance decreased, most notably in 2024. Improved model performance in 2020 persisted (albeit at lower levels), suggesting that both NDVIA and SWIA were indeed effective at improving model performance that year. Furthermore, for models with anomaly variables calculated with less than 30 zones, there is a banding pattern with identical model results. This is likely due to the uneven distribution of clusters in the KMeans algorithm, which despite creating up to 30 zones, tended to assign new cluster labels to smaller and more heterogenous zones outside the study areas. This resulted in identical zone delineation for low numbers of zones, and anomaly calculations thus also included a wider variety of pixels, resulting in less identified anomalies and many invalid (zero) anomaly values as model inputs. Zone complexity in the study areas tends to increase when clustering into more than 30 zones.

Figure 8 Improvement in phenological model adjusted R² score relative to best annual model trained without anomaly variables. Improvement in model adjusted R² shown for NDVIA (top), SWIA (second from top), combined (third from top), and only NDVIA/SWIA (bottom) for phenological models with different combinations of zones and upper threshold values.

Finally, model performance was also assessed with phenological anomaly variables, for which statistics were calculated separately during the vegetative and maturity phenological stages of the rice growth season, based on average transplant and harvest dates for farmer plots in each state (Fig. 8). Looking at the phenological anomaly variables, NDVIA increases model performance in 2020, 2022, and 2023. Notable improvements are seen for 2023, a year that otherwise shows very low model performance. Phenological SWIA data increases model performance more than seasonal SWIA data, and there are also notable improvements for 2021, but some zone/threshold combinations reduce model performance, particularly in 2023 and 2024. The combined models do show some improvements, particularly below a 25 % upper threshold, and around 30–40 zones. However, the 2022 models show greatly reduced performance relative to the seasonal anomaly variable models. When using only the anomaly variables and no meteorological variables, both the 2020 and 2024 models show improvement in model performance, while 2023 model performance is greatly reduced. This suggests that while phenological NDVIA variables do improve model performance for 2023, they are complementary to meteorological variables rather than effective completely on their own. Some banding patterns were observed at low threshold values for phenological SWIA values in 2024, likely resulting from too strict an upper threshold in that year failing to identify dekads with anomalies and resulting in zero-value model inputs. Similar banding patterns were observed at low zone numbers as described previously, likely due to the same effect. Although there is no clear “best” combination or application framework, these results suggest that anomaly-based indicators provide complementary information, particularly in years with weak meteorological-yield relationships.

Generally, while the optimal number of zones or upper threshold value was unclear, some clear patterns were observed. Low zone complexity (under 30 zones) showed a banding pattern for both NDVIA and SWIA variables, potentially indicating that rice plots are contained within limited zones with few differences, leading to similar NDVIA/SWIA values and lower adj-R² scores due to similar model inputs. Increased zone complexity improved model performance for models with only SWIA and meteorological variables, particularly when using seasonal rather than phenological stage statistics. However, a slightly higher model performance was observed around 50-70 zones for the seasonal SWIA and combined (2020 and 2022) and phenological NDVIA (2023) model, and around 70–100 zones for the phenological SWIA (2021) and NDVIA/SWIA models (2020 and 2024). This suggests that a zone complexity of around 70 zones could provide higher correlations between anomaly variables and rice yield changes across multiple years, however it is important to note that there is no clear “optimal” zone number, and the optimal zone range does vary year to year. While increasing the upper threshold percentage did increase the number of dekads identified as anomalies by both NDVIA and SWIA (Fig. 5), it did not necessarily improve the correlation between anomaly variables and rice yield change. For both seasonal and phenological NDVIA and SWIA variables in 2020, reduced model improvement scores are seen at upper threshold values above 35 %. For numerous combinations of input variables and years, upper threshold values of 15 % or 20 % show the highest model improvement scores. This may be because at these threshold values, the NDVI/SWI value that is identified as an “anomaly” is limited, resulting in a more specific definition of drought dekads, and leading to higher correlations and greater explainability of annual yield change values. Model errors vary across years and model types (Fig. A4), with the highest RMSE values observed for models using 2024 data and for NDVIA/SWIA (both seasonal and phenological). High RMSE values for 2024 data may be explained by the greater variance of yield data for that year, despite relatively high adj-R² scores and good explainability by meteorological variables alone.

3.4 Modelling of individual years and regions

The predictive capacity of the multivariate regression models to estimate rice yield changes from meteorological variable and anomaly indicator inputs depended on the year, and generally followed the patterns observed in Figs. 7–8. Figure 9 visualizes model predictions for a multivariate regression model fit on both meteorological variables and anomaly indicator values calculated with 70 zones and an upper threshold value of 20 %. Other combinations of zones and upper threshold values showed a similar grouping of results. While the mean values of both predicted and reported yields remained similar, in the years 2020, 2021, and 2023, the models tended to overestimate positive yield changes. In both years 2022 and 2024 yield change predictions were more evenly distributed, but with higher errors for negative yield changes, while positive yield changes (increasing yields year-over-year) followed a more clustered 1 : 1 relationship. Figure 9 suggests that the model predictive capacity was better at more neutral yield changes (closer to 0 t ha⁻¹), while for points towards the extremes of increasing or decreasing annual yield changes, the model was less capable of predicting accurate yield change values.

Figure 9 Annual rice yield changes predicted by model including meteorological variables (SPEI/SPI), and anomaly indicators (NDVIA, SWIA) calculated with 70 zones and an upper threshold value of 20 % versus annual rice yield change values as reported by farmers, per individual rice field. Red line indicates 1 : 1 reported vs predicted relationship.

Clear differences were observed at the state level in the performance of NDVIA/SWIA variables in describing annual changes in rainfed rice yields during drought periods. In the majority of surveyed locations, droughts reported by farmers aligned with NDVIA and SWIA results, but it was evident that the timing of the drought played an important role in yield impacts. For Niger state, the state with the highest number of surveyed rice plots, 2024 had the highest number of drought months reported by surveyed farmers (Fig. 10). During the 2024 growing season, both SPI3 and SPEI3 decreased around April, after transplanting started. This decrease corresponds to some high NDVIA values, indicating a negative anomaly in vegetation greenness over the rice plot pixels, consistent with the effects of an agricultural drought. The large decrease in SPI3/SPEI3 values continues throughout the rice tillering stage and into the panicle initiation stage in August, where both indices remain below −2.5, only increasing slightly during the flowering and grain-filling stages. During this period, SWIA anomaly values remain mostly constant throughout the panicle initiation and flowering stages. This indicates a deficit of soil moisture that can affect rainfed rice crops, especially during key stages for drought impacts on rice yields. Ultimately, the drought effects in 2024 are already quite well delineated by meteorological variables, while SWIA aligns well with farmer-reported drought periods. Rice yields in Niger state were relatively stable in 2023, with some yield increases, while SPI3/SPEI3 values were higher during all stages of rice growth. SWIA and NDVIA also showed very few anomaly values. SPI3/SPEI3 values were negative during the panicle initiation and maturity stages of 2022, and both NDVIA and SWIA values showed higher peaks, despite stable yields relative to the year prior. 2021 stands out as an outlier, with high NDVIA and SWIA values, as well as decreased SPI/SPEI values indicating agricultural drought during all stages of rice development, despite farmers reporting yield increases during that year. Finally, 2020 shows relatively high SPI3/SPEI3 values, but high SWIA values during the maturity stage align with yield decreases reported by farmers. This may in part explain the high improvement in model performance by adding SWIA variables for 2020 model inputs.

Figure 10 Summary graph comparing meteorological variables (SPI3, SPEI3; left), anomaly variables (NDVIA, SWIA; centre left), farmer-reported drought periods (centre-right), and annual/5-year mean yields (right) for Niger state. Anomaly variable values are derived from 70 zone and 20 % upper threshold equation parameters. State-specific rice growth calendar shows different stages of rice development according to field survey results: vegetative growth/tillering (green), panicle initiation (grey), and flowering/grain-filling (yellow). Statistically significant differences (Student's t-test, p < 0.05) between annual and 5-year mean yields is shown via letters (a, b) above boxplots.

High NDVIA values during the tillering phase, decreasing towards panicle initiation as SWIA values increase, were observed for all four states, especially during the years 2020, 2021, and 2024 (Figs. A5–A7). Nasarawa state notably showed much higher NDVIA/SWIA values, aligning well with both farmer-reported drought dekads and reductions in yield for almost all years. This state also had high SWIA values throughout almost the entire rice-growing season, even at reduced upper threshold values (under 30 %). The high SWIA values observed in this state are also decoupled from SPI3/SPEI3 patterns, suggesting that SWIA captures changes in soil water availability that are not easily explained by changes in precipitation/evaporation patterns. At the same time, these increases in SWIA align well with farmer-reported drought dekads. Ebonyi state shows stable SPI3/SPEI3 values during the panicle initiation and maturity stages of 2022 and 2023, which also correspond to NDVIA/SWIA zero values throughout the entire growing season for both years. This also aligns with almost no yield changes reported by farmers during those years. 2024 in Ebonyi state had a small drought early in the rice growing season, which is captured both by the meteorological indices as well as anomaly variables. However, this drought did not have a large effect on yield values. However, Ebonyi state had high NDVIA values early in the 2020 rice growing season, which could have contributed to the improvement in model performance for that year. Ogun state had fluctuating yield change values with many outliers which do not align well with neither SPI3/SPEI3 nor NDVIA/SWIA dynamics. This may have been due to a shifting cultivation system practiced in the state, which suggests that farmers may not have retained the same rice plots year-on-year. However, it's important to note that there were only a few plots in this state, suggesting too low a sample size to draw accurate links between drought detection indices and farmer-reported rice yields at the state level.

4 Discussion

Rainfed rice production in the study areas faces clear challenges due to drought, as evidenced by farmers' self-reported yield losses due to droughts over the duration of the study period. Existing meteorological indices are not always sufficient for monitoring drought impacts on rice yields, and regression models exhibit low predictive power for yield changes in years and locations with complex drought patterns. The four states surveyed in this study also have different management practices, environmental conditions, and sometimes even different rice varieties grown by farmers. FARO 44, the most common variety grown by surveyed farmers, is drought-resistant (Afiukwa et al., 2016), demonstrating that farmers are already adapting to changing climate conditions. However, both FARO 44 and OFADA, the second most common rice variety grown by surveyed farmers, show negative impacts on grain formation and grain yields when droughts occur during the panicle initiation and grain-filling stages (Akpoilih and Dada, 2025). Furthermore, survey results collected in this study showed that nearly 99 % of surveyed farmers use alert systems to receive early warnings of droughts, with the majority relying on radios and extension officers as their source of early warnings. These warnings are important for farmers to be able to react to coming droughts, either by accessing financial resources to absorb the economic shock, implementing basic irrigation if possible, or changing rice varieties or planting schedules; despite most smallholder rice farming households maintaining a low level of resilience to water insecurity shocks (Kim et al., 2017; Sanusi and Dries, 2025). The developmental stage of a rice crop at which agricultural drought occurs has different impacts on final yields, with drought stress throughout the panicle initiation stage resulting in severely reduced yields (Boonjung and Fukai, 1996). Furthermore, farmers' survey responses show that the majority (58 %) experienced major yield losses due to droughts, specifically during the panicle initiation stage of rice crop development. The reported patterns of yield changes (Fig. A1) and reported drought experiences are largely supported by the anomaly results in this study, specifically with SWIA capturing soil water anomalies during the later stages of the rice growing season and correlating well to negative yield changes, suggesting that monitoring seasonality of drought events is important for understanding their impacts on rice yields.

Annual temporal differences in NDVIA/SWIA patterns over the surveyed rice plots align with NiMet seasonal climate predictions and farmer bulletins, which show a variety of different rainfall patterns and crop impacts throughout the study areas. In 2020, there was a mid-season cessation of rains for 5–6 weeks in July and August across both Niger and Nasarawa, and a shorter dry spell in Ebonyi state (NAERLS and FMARD, 2020). The timing of these droughts directly overlapped with the reproductive and maturity stages of the rice crops in those states. These reported droughts do align with non-zero NDVIA and SWIA values in Nasarawa and Ogun states (Figs. A6, A7), as well as later season SWIA in Niger state (Fig. 8), potentially due to lag-time effect (Guo et al., 2026). Higher rainfall was reported across the country in 2021, although Niger and Nasarawa states experienced moderate dry spells (NAERLS and FMARD, 2021). These are also clearly visible in the mid-season increase in SWIA for both states. Ebonyi state (Fig. A6) also shows a peak in SWIA values, however this occurs earlier in the season and does not have a significant impact on yields. SWIA values also peak later in the rice growing season in Ogun state, also without a significant impact on yields. Severe floods occurred in 2022 from June to October, particularly affecting Nasarawa and Ebonyi states, with more minor impacts on Niger state (National Bureau of Statistics, 2023). This aligns with reduced anomaly values in Nasarawa state, and no observed anomaly values in Ebonyi state. NiMet bulletins and a post-season survey show that 2023 was an uneventful year in terms of drought, with normal rainfall generally, and some floods in Ebonyi state (NAERLS and FMARD, 2023; NiMet, 2023). Rice plot pixels in all states (except Nasarawa) show drastically reduced anomaly values in 2023 and almost no impact on rice yields. Finally, a nationwide dry spell was reported earlier in the rice season in 2024 with moderate to severe conditions in all states in this study (NiMet, 2024). However, this was later followed by excessive rainfall, sometimes resulting in flooding that affected rice fields (NAERLS and FMARD, 2024). For all four states, this is the year in which SPI3/SPEI3 values most closely follow the reported drought occurrences, whereas anomaly variables do not. Lower SPI3/SPEI3 values (particularly during the panicle initiation stage) are associated with significantly reduced rice yields in Ogun and Niger states, while the other two states show earlier drops in precipitation that do not align with reduced rice yields, further supporting the argument that the panicle initiation stage is an effective moment to use drought indices for monitoring drought impacts on rice yields (Boonjung and Fukai, 1996).

The drought patterns reported by NiMet in the study region may help to understand the multivariate regression model results for both the seasonal and phenological anomaly variables and rice yields. Despite reported drought impacts in 2020, SPI3/SPEI3 remain high, except in Ogun state. In contrast, SWIA/NDVIA peak, particularly in Niger and Nasarawa states. Phenological NDVIA variables are the only inputs that improve model performance in 2023, which could be partially explained by non-zero NDVIA values during the vegetative and panicle initiation stages in Nasarawa state, whereas there are very few anomaly values in other states. While it was difficult to understand exactly why phenological NDVIA improved model performance for some years, these results align with Zhang et al. (2019), who found that NDVI during early and mid-growth stages could be used to predict rice grain yields. Notably, Nasarawa state has the highest NDVIA/SWIA values for almost every year, and even the least sensitive zone/upper threshold combinations showed high anomaly values for Nasarawa. Lay et al. (2026) note that Nasarawa state is a particularly drought-prone region in North-Central Nigeria, which could explain the especially high anomaly values. Regional differences that were difficult to control for may account for the discrepancy in model adj-R² scores between years – especially as in Nasarawa, high SWIA values seem disconnected from the relatively high SPI/SPEI values during the tillering and panicle initiation stages. In Ogun state, decreased yields are particularly noticeable in 2020 and 2024, which have high SWIA values during the flowering/grain filling stage. Individual points also had large impacts on model adj-R² scores for some years, for example due to high yield variance in 2024 (which also explains higher RMSE values for that year; Fig. A1). Additionally, larger field sizes and higher yields were recorded for multiple surveyed locations in Nasarawa state, along with clay soil types that indicate a tendency towards water retention and external water supply to supplement rainfall. During fieldwork, it was also noted that some farmers (particularly in Ogun state) cultivated rice over two seasons and counted total rice yields, while others had external water sources that were difficult to document during fieldwork. These factors may have impacted the effectiveness of both meteorological and anomaly variables in predicting annual rice yield changes.

The ideal number of zones and threshold values for anomaly variable calculations had a large impact on both the calculation of EO indicator anomaly values, as well as the results of the multivariate regression analysis. Zones were calculated using seasonal NDVI patterns across the country, which did not show clear separation between zones, making it challenging to cluster pixels, as there was no “best” solution for the optimal number of clusters beyond basic aridity-based zoning. There is a minimum number of zones necessary for more accurate drought predictions, but this is due to the KMeans clustering algorithm which assigns smaller zones to more heterogenous pixels in the south, while grouping the majority of northern pixels together for lower zone numbers (Fig. 3). Once this minimum zone number (around 30–40 zones) is surpassed, drought pattern predictions remain largely similar for increasing numbers of zones. This is evident in Figs. 6 and A3, where anomaly values during high NDVIA/SWIA events in 2020 don't noticeably change beyond 40 zones. This aligns with Adefisan and Abatan (2015), who created a zoning system for Nigeria based on diverse drought and rainfall patterns, ultimately resulting in a three-zone system largely based on the rainfall gradient, similar to the spatial distribution of zone clusters beyond a zone number of 30. This finding goes against the use of 100+ zones for the determination of anomaly values, however previous studies have focused on index insurance applications, which have different end goals. While Kirchner et al. (2025) found that increased number of index insurance zones could provide more tailored insurance for farmers, they caution that too many zones can increase the variance within each zone and necessitate more field data to reduce insurance risk. Increasing threshold values increased the number of dekads with drought anomalies for both SWIA and NDVIA, although the greatest change in the number of drought dekads occurred around the 30 % threshold. Below that threshold however, drought dekad detection was more specific to drought as it related to yield changes for rainfed rice farmers and led to better model performance with the addition of anomaly variables. This is also observed in Figs. 6 and A3, where increasing upper threshold values do identify more anomaly pixels, but these anomalies are not necessarily more specific to the occurring drought events that cause said anomalies. This aligns with Bayissa et al. (2019), who defined drought intensity according to percentage deviation from combined drought index values, with 30 % and 20 % rankings corresponding to abnormally dry and moderate drought categories, respectively. This is relevant for applications such as index-based crop insurance, where area-specific risk and index thresholds can determine farmer payouts in the event of droughts (Bayissa et al., 2019). It further aligns with the previous use of anomaly variables for monitoring agricultural drought (Eze et al., 2020).

Both meteorological indices used, as well as the anomalies calculated from NDVI and SWI data, had large degrees of collinearity, ultimately affecting model accuracy. Initial principal component analyses showed similar patterns across non-anomaly variables, with SPEI more effectively capturing drought impacts on yield changes, aligning with previous studies (Tirivarombo et al., 2018). Significant adj-R² scores for meteorological indices align well with the successful use of SPI/SPEI to predict drought impacts on cereal yields in Africa as observed in previous studies (Makuya et al., 2024; Tirivarombo et al., 2018). Our results show that, despite collinearity among individual variables and across variable sources, the addition of both NDVIA and SWIA improved model accuracy in predicting yield change, especially in years when meteorological variables alone were insufficient. Despite the model improvement, the adj-R² values are still relatively low, indicating that there are more complex factors leading to rice yield changes than those that can be modelled with soil moisture and vegetation anomalies and meteorological variables alone. Anomaly variables do not add much yield explainability for years when meteorological variables already provide high model adj-R² scores (such as 2024). But for years when yield change is less accurately explained by meteorological variables, such as 2023, anomaly variables are valuable for identifying periods of drought and understanding the impacts of those drought periods on rice yields. These results align with previous research: Winkler et al. (2017) demonstrate that vegetation indices, in particular VCI, are able to detect much smaller-scale agricultural droughts and demonstrate different patterns of drought detection relative to purely meteorological indices. Similarly, we understand that the spatial resolution of NDVI/SWI data used in this study encompasses more than just the individual rice plots, or rice-producing regions. Multiple land-cover and land-use types exist within each pixel, ultimately influencing the spectral signal that indicates vegetation stress. Because small-scale irrigation or field improvement techniques (such as bunding) may be present at the field level, it is a challenge to fully disentangle the rice crop stress signal and thus the impacts of drought events on rice yields. Furthermore, Boschetti et al. (2013) note the disconnection between rainfall and vegetation greenness patterns in West Africa is often driven by human activity, which was difficult to account for entirely within this study.

Temporal patterns of EO indicator anomalies and their relationship to changes in rainfed rice yield differed across the states investigated in this study (Figs. A5–A7). Many years and states show NDVIA values peaking during the tillering stage of the rice season and decreasing as the rice crops approach the panicle initiation stage. The vegetation greenness anomalies during this period seem to be linked to changes in precipitation: SPI3 values were also low during each increased NDVIA period, however, SWI anomalies are only identified later in the season. This could be due to differences in the response time lag of the two anomaly variables to precipitation changes. The use of the SWI layer with a T latency factor of 40 d (Raml and Bauer-Marshallinger, 2018) further supports that SWIA values describe the impacts of low precipitation earlier in the season on later-season soil moisture, while NDVIA values may be more immediately responsive to changes in precipitation. Guo et al. (2026) found that the water use efficiency of African agroecosystems was particularly susceptible to deep soil moisture, and both NDVI- and SWI-derived indices had differing lag-times to meteorological drought (as determined by SPEI). In particular, the deep soil moisture (at 40–100 cm, including the SWI depth used in this study), was responsible for the most variation in agroecosystem water use efficiency, supporting our findings regarding SWIA as an effective variable to monitor agricultural drought impacts on rice yields. Lag time between soil moisture and SPI varies by SPI term, with SPI-1 and -3 showing little to no lag, whereas long-term SPI layers exhibit lags of multiple dekads (Cammalleri et al., 2024). Conversely, in the Horn of Africa, Odongo et al. (2023) found that arable land had a shorter (1 to 4 months) accumulation period during which the effects of meteorological drought would impact soil moisture. This may help explain the late-season increase in SWIA values observed in many states, leading to eventual yield impacts. The early- and mid-season non-zero NDVIA values observed in this study also align with Medida et al. (2023), who further note that soil texture was an important factor in the correlation between precipitation and vegetation stress. The aggregation of anomaly values into different phenological stages throughout the growing season did in fact improve model performance, particularly for NDVIA in 2023 and for SWIA in 2021. This suggests that EO indicator anomaly values at specific points in the growing season can be used for monitoring the impacts of drought on rice yields.

The variability in crop management practices among individual farmers across states, as well as the heterogeneity of crop and land cover for individual pixels, presented serious challenges to the development of a “one-fits-all” model, as evidenced by the lower performance of multiyear models. Furthermore, because drought occurrence at different stages of the rice growth cycle has different impacts on yields, it is also important to understand why NDVIA and SWIA exhibit different temporal patterns in identifying drought moments throughout the growing season. There are also challenges with using farmer-reported data on droughts, especially since the survey data were collected in 2025, potentially introducing recall bias in farmers' estimates of their rice yields and the impacts of past droughts dating back to 2019. This can be seen in the mismatch between drought as identified by SPI/SPEI patterns and farmer-reported drought dekads. Another limitation to this study was a spatial resolution mismatch between EO indices (1 km resolution), and farmer plots (which were rarely larger than one hectare), resulting in an EO signal that comes from more than just rice fields, and may include forests, arid regions, and urban areas. The two meteorological indices are also available at a much lower spatial resolution (approx. 28 km), which is a limitation for accurately matching the complexity between precipitation/evapotranspiration trends and their impacts on vegetation greenness and soil moisture. These spatial resolution mismatches are further complicated by the heterogeneity of smallholder farmer plots in much of sub-Saharan Africa. However, the lower spatial resolution is more relevant for drought monitoring at the landscape level and allows for some disregard of the impacts of management practices or field-level differences. The regression analyses between observations at both landscape and field scales thus also account for soil moisture or vegetation greenness in other land cover types, in addition to rainfed rice fields. This may have affected model performance, especially in years when farmer-reported rice yields did not align with broader landscape-level patterns. The use of linear regression models may also have limited the accuracy of the results, as in some years it is possible that the anomaly-yield change relationship was non-linear or too complex to be explained purely by anomaly variable statistics. Previous studies have identified better relationships between indices such as NDVI and rice crop yields using more complex models (Zhang et al., 2019).

Future research can help unpack these questions, while implementing anomaly-based drought indices in existing EWS can help assess their effectiveness in providing rice farmers with accurate and timely warnings, especially if anomaly indicators are integrated with seasonal forecasting approaches to deliver early drought predictions (De Vos et al., 2025). Future research could also investigate further the influence of rice variety selection for drought resilience of farmers in the region, as differences in growth duration, yields, and tolerance to water insecurity can have large impacts on farmers' economic outcomes (Kim et al., 2017; Omoare and Oyediran, 2020; Sanusi and Dries, 2025). The temporal patterns and correlations between EO indicator anomalies at different phenological stages and rice yields suggest that these indicators could be effective indices for monitoring how drought timing will affect farmers, although more research is needed to understand how regional differences affect yield model performance. Integrating yield modelling to understand the yield impacts of agricultural droughts could be a powerful next step for NiMet, allowing for more specific predictions of crop-level drought impacts. Using only NDVI and SWI in this study did not allow us to account for the utility of other EO indicators, such as Land Surface Temperature, which have been used to monitor drought impacts in past studies (Bayissa et al., 2019; Henchiri et al., 2020) and may be relevant for any follow-up research. Additional thermal data sources such as ECOSTRESS and novel indicators such as the Cooling Efficiency Factor Index (CEFI) also offer promising avenues for future research that could expand this framework to also monitor thermal anomalies (Wu et al., 2022; Zampieri et al., 2025). Additional studies could identify other anomaly indices that may be more discriminative in their identification of agricultural droughts that affect rice yields and investigate how remote sensing products with increased spatial resolution (for example Sentinel-2) could be integrated into this approach to provide analysis of drought impacts on rainfed rice production at the field level.

5 Conclusions

This study has demonstrated the added value of indices based on vegetation and soil moisture anomalies for monitoring drought impacts on rainfed rice yields in Nigeria and has connected temporal patterns in meteorological variables to farmer-reported yields. The use of field surveys distributed across multiple rice-growing states of Nigeria to train yield models was a unique application and gave relevant insights into how both meteorological and anomaly variables can be used to monitor yield changes due to drought. Although it remains unclear why some discrepancies persist between farmer-reported drought, annual changes in rice yields, and EO indicator-derived drought periods, anomaly indices provide an additional indicator of agricultural droughts with higher spatial resolution than traditional meteorological indices. Zonal clustering maximizes performance around 70 zones but should always be tailored to the specific area of study. Upper threshold values affect the number of drought pixels observed, but a lower threshold increases model specificity. The temporal differences in detected NDVIA/SWIA anomalies suggest that more complex temporal drought patterns are at play, and the discrepancy between non-zero NDVIA/SWIA values and reduced SPI/SPEI values suggests that the anomaly indices provide information on vegetation and soil conditions that are not covered by current standard drought EWS in Nigeria. SWIA is more robust than NDVIA at improving model performance across multiple low-skill years, however optimal effects can be seen where the use of both anomaly indicators improves yield prediction beyond the capacity of meteorological indices alone. Using anomaly indicators calculated at specific phenological stages, rather than over the entire growing season, improves performance, especially for NDVIA, suggesting that monitoring vegetation greenness at these stages can be an effective way to understand drought impacts on crop yields. Despite the need for further research to understand the reasons for interannual instability of using EO anomaly indicators for predicting on-farm rice yields, it is clear that EO anomaly indicators can improve model performance under specific conditions relative to meteorological indices and could be an effective tool for NiMet to improve existing EWS systems. The connection of field survey results to remote sensing anomaly indicators in this study demonstrates the effectiveness of such indicators for identifying agricultural drought impacts on rainfed rice yields in such complex, heterogeneous, and smallholder-dominated landscapes, but additional work could collect field-level data over multiple seasons rather than relying on farmers' recall of past season harvests. Further research and continued efforts to advance EWS employed by NiMet, including the integration of anomaly-based indicators, can ensure that farmers are better informed about the likelihood of upcoming droughts, enabling them to respond or prepare in any possible way. Implementation of the open-source anomaly indicator datasets could be pursued through future projects and through integration of open-source methods developed in this study into existing drought monitoring frameworks.

Appendix A

Table A1 Individual communities within local government areas and states where field surveys were conducted during this study.

Figure A1 Distribution of annual rice yield changes (left to right: 2020, 2021, 2022, 2023, 2024), filtered to remove top and bottom 5 % of observations, used to fit multivariate regression models in this study. Yield change values were derived from farmer surveys collected in the field.

Figure A2 Analysis into ideal number of clusters for Kmeans clustering algorithm showing gradient pattern between first two principal components (left) and continuously decreasing gradient of silhouette scores for 10 % pixel sample (right).

Figure A3 SWIA values for different combinations of zones and upper threshold values capturing anomalies in soil moisture during the dekad 11 August 2020.

Figure A4 Mean RMSE (t ha⁻¹) values for multivariate regression models using NDVIA, SWIA, and both as input variables for individual year models.

Figure A5 Summary graph comparing meteorological variables (SPI3, SPEI3; left), anomaly variables (NDVIA, SWIA; centre left), farmer-reported drought periods (centre-right), and annual/5-year mean yields (right) for Ebonyi state. Anomaly variable values are derived from 70 zone and 20 % upper threshold equation parameters. State-specific rice growth calendar shows different stages of rice development according to field survey results: vegetative growth/tillering (green), panicle initiation (grey), and flowering/grain-filling (yellow). Statistically significant difference (Student's t-test, p < 0.05) between annual and 5-year mean yields is shown via letters (a, b) above boxplots.

Figure A6 Summary graph comparing meteorological variables (SPI3, SPEI3; left), anomaly variables (NDVIA, SWIA; centre left), farmer-reported drought periods (centre-right), and annual/5-year mean yields (right) for Ebonyi state. Anomaly variable values are derived from 70 zone and 20 % upper threshold equation parameters. State-specific rice growth calendar shows different stages of rice development according to field survey results: vegetative growth/tillering (green), panicle initiation (grey), and flowering/grain-filling (yellow). Statistically significant difference (Student's t-test, p < 0.05) between annual and 5-year mean yields is shown via letters (a, b) above boxplots.

Figure A7 Summary graph comparing meteorological variables (SPI3, SPEI3; left), anomaly variables (NDVIA, SWIA; centre left), farmer-reported drought periods (centre-right), and annual/5-year mean yields (right) for Ebonyi state. Anomaly variable values are derived from 70 zone and 20 % upper threshold equation parameters. State-specific rice growth calendar shows different stages of rice development according to field survey results: vegetative growth/tillering (green), panicle initiation (grey), and flowering/grain-filling (yellow). Statistically significant difference (Student's t-test, p < 0.05) between annual and 5-year mean yields is shown via letters (a, b) above boxplots.