Spatiotemporal characterization and regional forecasting of the temperature and humidity index: an exploratory case study on dairy production in Villa Clara, Cuba

Full text

Text extracted from the article PDF for reading and search. Tables, figures and equations are in the PDF, which remains the authoritative version.

1. Introduction

Heat stress negatively affects the welfare, health, feed intake, milk production, and reproductive performance of dairy cattle [1- 3]. These effects are particularly significant in tropical and subtropical production systems, where high air temperatures and humidity can persist for extended periods and where on-farm access to environmental monitoring and cooling infrastructure may be uneven [4-6].

The temperature-humidity index (THI) is one of the most widely used indicators to describe the environmental heat load on livestock [7-9]. Originally developed to describe human thermal comfort (Thom, 1959) [10], it was later shown by Hahn & McQuigg (1967) [11] and Hahn & Osburn (1970) [12] to be clearly associated with losses in milk production; since then, it has been widely adopted in livestock bioclimatology (Nienaber & Hahn, 2007) [13]. However, THI thresholds are not biologically universal: their interpretation depends on genotype, milk production level, acclimatization, housing, solar radiation, wind speed, access to shade or cooling, and the duration of exposure (Khalifa, 2003) [14]. Consequently, fixed thresholds should be treated as screening values rather than universal physiological limits. Cuban dairy systems face the combined challenges of a tropical climate, distinct wet and dry seasons, a spatially uneven station network, and locally adapted yet heterogeneous cattle populations [7,8,15]. Long-term meteorological records can support the identification of recurring periods of heat exposure and the development of parsimonious statistical forecasting tools for livestock management [3,5]. At the same time, a gridded or regional product derived from a sparse station network represents an interpolated or aggregated estimate; its uncertainty must be quantified before it is used as information at the farm scale [1,16,17]. The study had three objectives: (i) to characterize the spatial and seasonal variability of the THI in Cuba using daily observations from 23 INSMET stations grouped into three geographic regions during the 1993–2013 period; (ii) to compare the out-of-sample descriptive accuracy of an STL+ETS model against automatic ARIMA models for forecasting the THI at weekly and monthly scales in each region; and (iii) to examine, as an exploratory analysis, whether monthly THI was associated with milk production indicators recorded at a representative farm (La Vitrina, Villa Clara). The third objective was considered exploratory because herd management, physiological, and health-related covariates were not available to perform a causal adjustment.

2. Materials and Methods

2.1. Study area and Station Network

The study covered the Cuban archipelago. Daily meteorological data from 23 INSMET surface stations for the 1993–2013 period (21 years) were grouped into three regions based on geographic coverage: Group ITH1 (Eastern), extending from Camagüey to Guantánamo; Group ITH2 (Central), spanning from Matanzas to Ciego de Ávila; and Group ITH3 (Western), extending from Pinar del Río to Isla de la Juventud.

Six stations were excluded from the final analysis due to sample sizes that were unequal compared to the rest of their group: 78361 (Group ITH1); 78308, 78335, and 78344 (Group ITH2); and 78315 and 78316 (Group ITH3). Station elevations ranged from 0.90 m (Cape San Antonio, Pinar del Río) and 1 130.00 m (Gran Piedra, Santiago de Cuba); the lowest recorded ITH values were concentrated at stations situated above approximately 400–750 m a.s.l. (e.g., Topes de Collantes, 767.33 m; Gran Piedra, 1130.00 m; Palenque de Yateras, 440.00 m).

The exploratory production case study utilized monthly records from the La Vitrina livestock enterprise in Villa Clara province, Cuba, which is representative of the region's tropical dairy systems.

2.2. Meteorological Observations and Quality Control

Air temperature and relative humidity data were obtained from the INSMET Climate Center. Quality control assessed spatial, temporal, and statistical continuity and consistency, including the individual review of suspected cases (outliers) and the calculation of statistics regarding central tendency, dispersion, and frequency. Data were organized into daily (365 values/station/year), weekly (52 values/station/year), monthly (12 values/station/year), and annual (21 data points) series for each of the 23 stations.

Missing values were identified and replaced with the series mean (mean imputation). The final descriptive dataset comprised 206 877 daily THI observations: 189 392 in the THI1 Group, 90 398 in the THI2 Group, and 98 087 in the THI3 Group.

2.3. Temperature and Humidity Index

THI (referred to as ITH in the source records) was calculated from dry-bulb air temperature and relative humidity following Thom (1959) and NRC (1971), as reported by Nienaber and Hahn (2007) [13]:

THI = 0.8 · Tdb + RH · (Tdb − 14.4) + 46.4 where Tdb is dry-bulb air temperature (°C) and RH is relative humidity expressed as a decimal fraction (0–1). Calculated THI values were classified using two complementary schemes reported in the thesis: (i) the four-category operational scale of Cartwright et al. (2023) [6] applied directly to THI values (Table 1a), used in Section 3.2 to describe the proportion of “alert,” “danger,” and “emergency” conditions; and (ii) the air-temperature-based thermal classification of Kulicov & Rudnev (1987) [18] (Table 1b), used to characterize the animals’ physiological state independently of the composite index.

The calculated THI values were classified using two complementary schemes reported in the thesis: (i) the four-category operational scale of Cartwright et al. (2023) [6], applied directly to the THI values (Table 1a) and used in Section 3.2 to describe the proportion of "alert," "danger," and "emergency" conditions; and (ii) the air temperature-based thermal classification of Kulicov and Rudnev (1987) [18] (Table 1b), used to characterize the animals' physiological state independently of the composite index.

Table 1: Operational classification of THI values (Cartwright et al., 2023) [6] THI range Classification ≤ 74 Normal 75–78 Alert 79–83 Danger ≥ 84 Emergency Daytime (Tfoto) and nighttime (Tnicto) temperatures, used to characterize diurnal thermal exposure, were calculated following Castillo & Castellví (2001) [19]:

Tfoto = Tx − (Tx − Tn) / 4 Tnicto = Tn + (Tx − Tn) / 4 where Tx and Tn are the daily maximum and minimum temperatures.

Independently of the THI value itself, ambient air temperature was also classified using the thermal-comfort ranges of Kulicov and Rudnev (1987) [18] (Table 2), which describe animal state as a function of temperature alone.

Table 2: Classification of air-temperature regime and associated animal thermal state (Kulicov & Rudnev, 1987) [18] Lower limit (°C) Upper limit (°C) Classification Animal state — 16–24 Ideal Best condition, maximum productivity 10–16 24–27 Favorable Active animals, satisfactory state, normal productivity 7–10 27–30 Not fully favorable Insignificant depression, mainly at midday — 30–35 Unfavorable Visible depression, reduced productivity — 35–40 Very unfavorable Significant depression, minimal productivity

2.4. National-Scale Gridded Thi Mapping (Short-Range

In parallel with the regional statistical analysis described below, the source methodology also describes a short-range (24, 48, 72, 96, and 120 h) operational THI mapping service, integrating six-hour forecast fields (12 UTC initialization) of the GFS (Global Forecast System) numerical weather-prediction model with the station-based THI/thermal-classification scheme, following the approach of Domínguez et al. (2017) [20]. This component used a national network of 68 INSMET stations, distinct from, and broader than, the 23- station regional network used for the time-series characterization and forecasting analysis reported in Sections 3.2–3.3.

Station temperature and specific humidity were reduced to mean sea level using the static-equation formulation of Matveev (1984) [21]:

where T0 is temperature reduced to sea level (K), Test is station-level temperature (K), hest is station elevation, and P1, P0 are station-level and sea-level-reduced pressure (hPa); S is specific humidity computed from station-level pressure and dew point, and S0 is specific humidity reduced to sea level. Reduced values were interpolated onto a regular 118 × 44- node grid (approximately 10-km spacing) covering the Cuban archipelago using Shepard’s inverse-distance-weighting method (power µ = 2) (Shepard, 1968) [22], then re-elevated to each grid node’s true height. For graphical display only, the 118 × 44 grid was further interpolated to a 718 × 324- pixel map using bicubic spline interpolation; this step did not generate additional independent observations and did not enter the statistical evaluation reported in Section 3.

2.5. Regional Time-Series Forecasting Models

For each of the three regions (Eastern, Central, and Western), the aggregated weekly and monthly ITH series were examined for stationarity using the Augmented Dickey– Fuller test and the KPSS test (Table 3). Both tests indicated that all six series (three regions × two-time scales) were stationary (Dickey–Fuller: p < 0.05 in all cases).

Table 3: Dickey–Fuller and KPSS stationarity tests by region and time scale Escale Group Dickey–Fuller Lag p KPSS KPSS truncation lag KPSS p Daily I −27.509 55 0.01 1.3020 95 0.01 Daily II −12.249 48 0.01 1.4743 78 0.01 Daily III −24.970 48 0.01 1.8322 77 0.01 Weekly I −8.6862 10 0.01 2.1199 7 0.01 Weekly II −8.7354 10 0.01 0.0760 7 0.01 Weekly III −8.7835 10 0.01 4.3920 7 0.01 Monthly I −5.4499 6 0.01 2.1262 3 0.01 Monthly II −8.6185 6 0.01 0.1005 3 0.01 Monthly III −5.1948 6 0.01 3.6182 3 0.01 Two forecasting approaches were compared: (i) an automatic seasonal ARIMA model, with orders (p,d,q) (P,D,Q) selected via the auto. Arima function in R, and (ii) an STLETS(M,N,N) model, combining Seasonal-Trend decomposition using LOESS with a multiplicative-error, no-trend, no-seasonality exponential-smoothing model applied to the seasonally adjusted series, following the Holt–Winters family of exponential-smoothing methods [...] Forecast accuracy was assessed using mean error (ME), root mean squared error (RMSE), mean absolute error (MAE), mean percentage error (MPE), mean absolute percentage error (MAPE), and mean absolute scaled error (MASE), following Hyndman & Koehler (2006) [23]. Model residuals were evaluated using ACF and PACF diagnostics.

where ℓt is the level at time t and ηt is a multiplicative error term, normally distributed with mean 0 and variance σ². Model residuals were evaluated using ACF and PACF diagnostics.

2.6. Exploratory Analysis of Dairy-Production Records

Monthly production and demographic records from the La Vitrina farm, Villa Clara, Cuba, number of milking cows, total milk production (L), total births, general mortality, and calf mortality, were compiled from the corresponding livestock indicator bulletins.

Associations between monthly THI and each production/demographic indicator were assessed using Pearson correlation coefficients (two-tailed significance). Losses attributable to theft or slaughter were excluded from the mortality counts.

3. Results

3.1. Data Completeness and Station Distribution

After quality control, the final data set comprised 206,877 daily THI observations from 23 stations (six additional stations excluded for unequal sample size, as noted in Section 2.1). Group-level descriptive statistics are show in table 4.

Table 4: Descriptive statistics of daily THI by region Groups n Range Minimum Maximum Mean SD ITH1 (Oriental) 189 392 30 55 85 75.39 3.784 ITH2 (Central) 90 398 40 43 83 73.85 4.621 ITH3 (Occidental) 98 087 31 51 83 74.10 4.521 Overall, the series ranged from 47 to 85 (mean = 74.63; n = 206,877). The absolute minimum (THI = 47) occurred on 10– 11 January 2010 at stations 78342 (Topes de Collantes, Sancti Spíritus, 767.33 m a.s.l.) and 78366 (Gran Piedra, Santiago de Cuba, 1,130.00 m a.s.l.), coinciding with the passage of the season’s tenth cold front under La Niña and negative North Atlantic Oscillation conditions. The absolute maximum (THI = 85) was sustained for three consecutive days (25–27 August 2010) at Punta Lucrecia, Holguín (4.00 m a.s.l.), associated with increased low- and mid-level humidity, strong afternoon heating, and afternoon sea-breeze convergence. Station elevation across the full 23-station network (see Section 2.1) ranged from 0.90 to 1,130.00 m.

3.2. Spatial and Seasonal Variability of THI

THI followed a unimodal annual cycle in all three regions, with maxima in July–August (values consistently above 75, classified as alert/danger for the animals) and minima from December to March (values below 74, classified as normal conditions). At the weekly scale, the highest values occurred around week 31 (late July–early August), consistent with the monthly pattern.

A spatial gradient of decreasing thermal comfort from west to east was observed: Group ITH1 (Oriental) showed the least favorable conditions, followed by Group ITH2 (Central), with Group ITH3 (Occidental) showing the most favorable conditions overall. This deterioration was more pronounced during July–August regardless of region. The most favorable stations in each group were located above approximately 400 m a.s.l. (e.g., Palenque de Yateras, 440.00 m, and Pinares de Mayarí, Group ITH1; Topes de Collantes, 767.33 m, Group ITH2), consistent with the expected effect of altitude and montane climate on temperature and humidity. Anomaly-detection (X-bar) analysis of daily THI series identified 12.5% (21,264 of 170,566 cases), 15.7%, and 21.6% (24,161 cases) of observations outside the established control limits for Groups ITH1, ITH2, and ITH3, respectively, predominantly below the lower control limit. Correspondingly, 87.5%, 84.3%, and 78.4% of observations for Groups ITH1, ITH2, and ITH3 fell within the “alert” or “danger” range for animal thermal comfort, indicating that unfavorable conditions predominate for most of the year across all three regions, with Group ITH3 (Occidental) showing the smallest proportion of such conditions. Weekly CUSUM analysis showed an increasing trend in cumulative anomalies for Group ITH1 (12.5% of points outside limits, mostly above the upper bound), no defined trend for Group ITH2 (no points outside limits), and both excess and deficit anomalies for Group ITH3, concentrated below the lower bound.

3.3. Regional THI Forecasting Performance

Table 4 summarizes forecast-accuracy statistics for the STL+ETS(M,N,N) and automatic-ARIMA models, by region and time scale.

3.4. Exploratory Association with Dairy-Production

Records At the La Vitrina farm, total milk production peaked during the rainy season (July–August) and was lowest during the dry season (February–April), despite THI being classified as “alert” or “danger” in both periods. Monthly THI was significantly and positively correlated with total milk production (r = 0.496, p = 0.043, n = 17) and with per-cow milk production (r = 0.489, p = 0.046, n = 17) (Table 5). No statistically significant correlations were found between THI and calf/heifer sales, births, general mortality, or calf mortality (all p > 0.05), although weak positive correlations were observed with sales-related variables (r = 0.272–0.293) and weak negative correlations with mortality and births (r = −0.126 to −0.327).

Table 5: Pearson correlations between monthly THI and production/demographic indicators, La Vitrina farm Variables r p (bilateral) n Potential litres per cow (LXVR) 0.278 0.280 17 Calf/heifer sales, real (VESTR) 0.293 0.253 17 Births, real (NACR) 0.103 0.704 16 General mortality, real (MTR) −0.327 0.201 17 Calf mortality, real (MCR) −0.015 0.955 16 Total milk production, real (LP) 0.496 0.043 17 Milk per cow, real (LR) 0.489 0.046 17 Milking cows (VOP) 0.154 0.555 17 Cows in herd (VOR) 0.199 0.443 17 Potential litres per cow, plan (LXVP) 0.272 0.291 17 Calf/heifer sales, plan (VESTP) 0.276 0.284 17 Births, plan (NACP) −0.126 0.643 16 General mortality, plan (MTP) −0.265 0.272 19 Calf mortality, plan (MCP) −0.051 0.851 16 Bold values are significant at p < 0.05. These results should be interpreted cautiously: the analysis used aggregate monthly farm-level records and did not include individual-animal physiological status, feeding, disease incidence, pasture availability, shade, or cooling interventions. It therefore describes temporal co-variation, not a causal effect of heat stress on milk production.

4. Discussion

This study documented a pronounced seasonal cycle and a west-to-east spatial gradient in THI across Cuba, and compared two forecasting approaches for regional THI series. The STL+ETS(M,N,N) model was more accurate than automatic ARIMA models by conventional scale-dependent and percentage error metrics (RMSE, MAE, MAPE, MASE) in the majority of the region–time-scale combinations examined, although it showed greater bias (ME) in several cases, a trade-off that should be weighed against the intended operational use of the forecasts (e.g., whether systematic bias or overall error magnitude is more consequential for livestock-management decisions).

The apparent paradox of milk production peaking during the period of least favorable thermal conditions (July–August) may reflect several non-mutually exclusive factors, which have also been discussed in the context of other tropical systems [5,7,23]. These include the genetic adaptation of locally raised cattle to heat and humidity; physiological adjustments (e.g., increased respiratory rate and peripheral blood flow) that partially offset productivity losses, and possibly the dominant factor, the greater availability of feed and pasture during the rainy season, which coincides with the period of highest heat load [2,3,9]. This aligns with the observation that the Temperature-Humidity Index (THI) showed a moderate yet significant association only with milk production variables, and no significant association with mortality, births, or animal sales; this suggests that factors other than thermal comfort (e.g., nutrition, management) are more decisive for those outcomes in this system [24,25,29]. The findings contrast with reports describing significant associations between the season of the year and calf mortality in Mexican livestock systems [26,28]; however, such relationships depend on the context and the system and should not be generalized without local validation [26,28,30]. Two components of the broader agrometeorological service described in the original methodology, the "leave-one-station-out" validation of spatial interpolation and the quantitative assessment of the skill of GFS/NCEP-based THI forecasts for the 24–120 h range, could not be evaluated. Given that the national grid (68 stations, 10 km nominal spacing) relies on a station network that is considerably sparser than the grid resolution, interpolation uncertainty likely varies by region and should not be assumed negligible; reporting this uncertainty is a priority before the gridded product is used for farm level decision making.

5. Conclusion

Daily THI in Cuba (1993–2013, 23 stations) exhibited a clear unimodal seasonal cycle, peaking in late July to early August, and a spatial gradient of decreasing thermal comfort from west to east, with the most favorable conditions recorded at mountain stations. At weekly and monthly scales, the STL+ETS(M,N,N) model outperformed automatic ARIMA models in terms of scale-dependent and percentage-based accuracy metrics across most region-scale combinations evaluated, albeit with higher bias in several instances. Thermal comfort is one of several contributing factors, rather than the dominant one regarding the production and demographic outcomes examined. Validating the spatial interpolation and the short-term forecasting service based on numerical weather prediction models remains a necessary step before both can be considered operationally validated.

Institutional Review Board Statement Not applicable.

Transparency The authors confirm that the manuscript is an honest, accurate, and transparent account of the study, and that any discrepancies from the study as planned have been explained. This study followed all ethical practices during writing.

Competing Interests The authors declare that they have no competing interests.

Author contributions Ismabel María Domínguez Hurtado: Conceptualization, Data curation, Formal analysis, Investigation, Visualization, writing – original draft, writing – review and editing. César Soto Valero: Formal analysis (statistical analysis). David del Valle Laveaga: Writing – review and editing, Funding acquisition. Rigoberto Fimia Duarte: Writing – review and editing. Neilys González Benítez: Methodology, Writing – review and editing.