WEFNC Nexus Analysis for Agricultural Systems

Lalit Kumar Rolaniya

Introduction

The wefnexus package provides a comprehensive toolkit for analyzing Water-Energy-Food-Nutrient-Carbon (WEFNC) nexus interactions in agricultural production systems. This vignette demonstrates a complete workflow using a conservation agriculture dataset from western Rajasthan, India.

If you have replicated field data and want the whole analysis in one step, jump to Section 7. Sections 1 to 6 show the individual modules.

Load Sample Data

data(arid_pulse_nexus)
d <- arid_pulse_nexus
str(d)
#> 'data.frame':    6 obs. of  26 variables:
#>  $ treatment          : chr  "CT+Flood" "CT+Sprinkler" "ZT+Drip" "ZT+Drip+Residue" ...
#>  $ tillage            : chr  "CT" "CT" "ZT" "ZT" ...
#>  $ irrigation         : chr  "Flood" "Sprinkler" "Drip" "Drip" ...
#>  $ residue            : logi  FALSE FALSE FALSE TRUE FALSE TRUE
#>  $ grain_yield        : num  980 1120 1280 1380 1450 1580
#>  $ straw_yield        : num  1450 1680 1850 2050 2180 2350
#>  $ irrigation_applied : num  280 220 160 155 135 130
#>  $ effective_rainfall : num  185 185 185 185 185 185
#>  $ total_water        : num  465 405 345 340 320 315
#>  $ crop_et            : num  390 355 310 305 290 285
#>  $ energy_input       : num  8950 7850 6480 6950 5850 6380
#>  $ energy_output_grain: num  14406 16464 18816 20286 21315 ...
#>  $ energy_output_straw: num  18125 21000 23125 25625 27250 ...
#>  $ n_applied          : num  20 20 20 20 20 20
#>  $ p_applied          : num  40 40 40 40 40 40
#>  $ k_applied          : num  0 0 0 0 0 0
#>  $ n_uptake           : num  32.5 37.8 43.2 48.5 52.1 56.8
#>  $ p_uptake           : num  4.8 5.5 6.3 7.1 7.6 8.2
#>  $ grain_n_uptake     : num  24.5 28.6 33 37.2 40.1 43.8
#>  $ diesel_use         : num  65 55 38 40 32 35
#>  $ electricity_kwh    : num  180 120 75 70 55 50
#>  $ soc_pct            : num  0.32 0.34 0.38 0.42 0.4 0.45
#>  $ bulk_density       : num  1.58 1.56 1.52 1.48 1.5 1.46
#>  $ ghg_emission       : num  1850 1520 1180 1250 1020 1100
#>  $ cost_cultivation   : num  22500 24800 21500 23200 20800 22500
#>  $ gross_return       : num  53900 61600 70400 75900 79750 ...

1. Water Module

# Water Use Efficiency
wue <- water_use_efficiency(d$grain_yield, d$total_water)

# Crop Water Productivity
cwp <- crop_water_productivity(d$grain_yield, d$crop_et)

# Water Footprint (green + blue)
wf <- water_footprint(
  green_water = d$effective_rainfall,
  blue_water = d$irrigation_applied,
  yield = d$grain_yield
)
wf
#>   green_wf blue_wf grey_wf total_wf
#> 1     1.89    2.86       0     4.75
#> 2     1.65    1.96       0     3.61
#> 3     1.45    1.25       0     2.70
#> 4     1.34    1.12       0     2.46
#> 5     1.28    0.93       0     2.21
#> 6     1.17    0.82       0     1.99

2. Energy Module

# Total energy output
e_out <- d$energy_output_grain + d$energy_output_straw

# Energy Use Efficiency
eue <- energy_use_efficiency(e_out, d$energy_input)

# Energy Return on Investment (EROI)
eroi_values <- eroi(e_out, d$energy_input)

# Net Energy Balance
ne <- net_energy(e_out, d$energy_input)

# Compare treatments
data.frame(Treatment = d$treatment, EUE = eue, EROI = eroi_values,
           Net_Energy = ne)
#>         Treatment  EUE EROI Net_Energy
#> 1        CT+Flood 3.63 3.63      23581
#> 2    CT+Sprinkler 4.77 4.77      29614
#> 3         ZT+Drip 6.47 6.47      35461
#> 4 ZT+Drip+Residue 6.61 6.61      38961
#> 5         PB+SSDI 8.30 8.30      42715
#> 6 PB+SSDI+Residue 8.24 8.24      46221

Energy budget from physical quantities

Most field records hold litres of diesel, kilograms of fertilizer and labour hours rather than megajoules. energy_budget() converts them with the energy_coefficients table and splits input energy into direct, indirect, renewable and non-renewable shares.

head(energy_coefficients[, c("input", "unit", "mj_per_unit", "form",
                             "source")])
#>                input    unit mj_per_unit     form        source
#> 1   human_labour_man   man-h        1.96   direct     renewable
#> 2 human_labour_woman woman-h        1.57   direct     renewable
#> 3       bullock_pair  pair-h       10.10   direct     renewable
#> 4             diesel       L       56.31   direct non-renewable
#> 5        electricity     kWh       11.93   direct non-renewable
#> 6          machinery      kg       62.70 indirect non-renewable

inputs <- data.frame(
  treatment = c("CT+Flood", "ZT+Drip"),
  human_labour_man = c(320, 240),
  diesel = c(65, 38),
  electricity = c(180, 75),
  nitrogen = c(20, 20),
  phosphorus = c(40, 40),
  seed = c(20, 20),
  herbicide = c(0, 0.75),
  grain = c(980, 1280),
  straw = c(1450, 1850)
)
eb <- energy_budget(inputs)
eb$summary[, c("treatment", "total_input_MJ", "direct_pct",
               "renewable_pct", "EUE", "energy_intensity_MJ_kg")]
#>   treatment total_input_MJ direct_pct renewable_pct  EUE energy_intensity_MJ_kg
#> 1  CT+Flood        8384.75      76.74         10.99 3.88                   8.56
#> 2   ZT+Drip        5633.43      62.22         13.57 7.45                   4.40

Energy equivalents vary between sources. To use your own, edit a copy of the table and pass it as coefficients.

plot(eb)

plot(eb, type = "source", percent = TRUE)

3. Food Module

# Harvest Index
bio_yield <- d$grain_yield + d$straw_yield
hi <- harvest_index(d$grain_yield, bio_yield)

# Protein Yield (assuming 22% protein in pulses)
py <- protein_yield(d$grain_yield, protein_content = 22)

data.frame(Treatment = d$treatment, HI = hi, Protein_kg_ha = py)
#>         Treatment     HI Protein_kg_ha
#> 1        CT+Flood 0.4033         215.6
#> 2    CT+Sprinkler 0.4000         246.4
#> 3         ZT+Drip 0.4089         281.6
#> 4 ZT+Drip+Residue 0.4023         303.6
#> 5         PB+SSDI 0.3994         319.0
#> 6 PB+SSDI+Residue 0.4020         347.6

4. Nutrient Module

# Partial Factor Productivity of N
pfp <- partial_factor_productivity(d$grain_yield, d$n_applied)

# Nutrient Harvest Index for N
nhi <- nutrient_harvest_index(d$grain_n_uptake, d$n_uptake)

data.frame(Treatment = d$treatment, PFP_N = pfp, NHI_N = nhi)
#>         Treatment PFP_N  NHI_N
#> 1        CT+Flood  49.0 0.7538
#> 2    CT+Sprinkler  56.0 0.7566
#> 3         ZT+Drip  64.0 0.7639
#> 4 ZT+Drip+Residue  69.0 0.7670
#> 5         PB+SSDI  72.5 0.7697
#> 6 PB+SSDI+Residue  79.0 0.7711

5. Carbon Module

The default direct N2O emission factor is the IPCC Tier 1 value of 0.01. For dry climates the 2019 IPCC Refinement gives 0.005, which can be set with ef_n2o.

# Carbon Footprint (IPCC AR6 GWP: CH4=27, N2O=273)
cf <- carbon_footprint(
  diesel_use = d$diesel_use[1],
  electricity_use = d$electricity_kwh[1],
  n_fertilizer = d$n_applied[1],
  p_fertilizer = d$p_applied[1],
  yield = d$grain_yield[1],
  ef_n2o = 0.005
)
cf$breakdown
#>         source emission_kg_CO2eq share_pct
#> 1       Diesel            174.20      34.1
#> 2  Electricity            130.86      25.6
#> 3 N_fertilizer             99.20      19.4
#> 4 P_fertilizer             64.40      12.6
#> 5 K_fertilizer              0.00       0.0
#> 6    Pesticide              0.00       0.0
#> 7         Seed              0.00       0.0
#> 8   N2O_direct             42.90       8.4
#> 9          CH4              0.00       0.0

# Soil Carbon Stock
soc <- soil_carbon_stock(d$soc_pct, d$bulk_density, depth = 30)

# Global Warming Potential (example)
gwp100 <- global_warming_potential(co2 = 500, ch4 = 10, n2o = 2)
gwp20 <- global_warming_potential(co2 = 500, ch4 = 10, n2o = 2,
                                   time_horizon = "20yr")

6. Nexus Integration

Every treatment in this trial received the same starter dose of 20 kg N/ha, so partial factor productivity would only mirror yield. The nutrient dimension therefore uses internal N utilization efficiency (nutrient_metric = "iue"), and carbon is scored per kg of grain.

# Full Nexus Summary
ns <- nexus_summary(
  yield = d$grain_yield,
  water_consumed = d$total_water,
  energy_input = d$energy_input,
  energy_output = e_out,
  n_applied = d$n_applied,
  n_uptake = d$n_uptake,
  carbon_emission = d$ghg_emission,
  treatment_names = d$treatment,
  nutrient_metric = "iue",
  carbon_basis = "intensity"
)
ns[, c("treatment", "EROI", "nexus_index")]
#> WEFNC Nexus Analysis Results
#> ---------------------------------------- 
#> Treatments: 6 
#> Nexus Index range: 0.2 - 0.797 
#> Best performer: PB+SSDI+Residue 
#> 
#>        treatment EROI nexus_index
#>         CT+Flood 3.63      0.2000
#>     CT+Sprinkler 4.77      0.3852
#>          ZT+Drip 6.47      0.6493
#>  ZT+Drip+Residue 6.61      0.6140
#>          PB+SSDI 8.30      0.7230
#>  PB+SSDI+Residue 8.24      0.7974

Radar Plot

scores <- as.matrix(ns[, c("W_score", "E_score", "F_score",
                             "N_score", "C_score")])
nexus_radar(scores, treatment_names = d$treatment)

Sustainability Scoring

nss <- nexus_sustainability_score(
  water_score = ns$W_score,
  energy_score = ns$E_score,
  food_score = ns$F_score,
  nutrient_score = ns$N_score,
  carbon_score = ns$C_score
)
nss[, c("nexus_score", "category")]
#> WEFNC Sustainability Assessment
#> ---------------------------------------- 
#>   Score: 0.2000 -> Unsustainable
#>   Score: 0.3852 -> Unsustainable
#>   Score: 0.6493 -> Moderately Sustainable
#>   Score: 0.6140 -> Moderately Sustainable
#>   Score: 0.7230 -> Moderately Sustainable
#>   Score: 0.7974 -> Moderately Sustainable
#> 
#>  nexus_score               category
#>       0.2000          Unsustainable
#>       0.3852          Unsustainable
#>       0.6493 Moderately Sustainable
#>       0.6140 Moderately Sustainable
#>       0.7230 Moderately Sustainable
#>       0.7974 Moderately Sustainable

Objective weights and TOPSIS ranking

Equal weights are a value judgement. nexus_weights() derives weights from the data (entropy, CRITIC or PCA), and nexus_topsis() ranks treatments on raw indicators in their own units.

ind <- data.frame(
  WUE = ns$WUE_kg_mm,
  EUE = ns$EUE_ratio,
  Yield = ns$yield_kg_ha,
  IUE = ns$IUE_kg_kgN,
  CI = ns$carbon_intensity
)
benefit <- c(TRUE, TRUE, TRUE, TRUE, FALSE)  # lower carbon intensity is better

rbind(
  entropy = nexus_weights(ind, "entropy", benefit),
  critic  = nexus_weights(ind, "critic", benefit),
  pca     = nexus_weights(ind, "pca", benefit)
)
#>               WUE       EUE     Yield       IUE        CI
#> entropy 0.1800941 0.1803205 0.1802616 0.3166415 0.1426824
#> critic  0.1169555 0.1252187 0.1161999 0.5209748 0.1206510
#> pca     0.2036357 0.2024941 0.2029193 0.1936396 0.1973113

nexus_topsis(ind, weights = "entropy", benefit = benefit,
             treatment_names = d$treatment)
#>         treatment d_ideal d_anti_ideal closeness rank
#> 1        CT+Flood  0.1029       0.0104    0.0919    6
#> 2    CT+Sprinkler  0.0722       0.0341    0.3208    5
#> 3         ZT+Drip  0.0383       0.0683    0.6408    4
#> 4 ZT+Drip+Residue  0.0315       0.0739    0.7014    3
#> 5         PB+SSDI  0.0157       0.0956    0.8588    2
#> 6 PB+SSDI+Residue  0.0104       0.1025    0.9076    1

Is the ranking robust?

A recommendation is more convincing when it does not depend on one choice of weights. nexus_rank_compare() ranks the treatments eight ways (composite index and TOPSIS, each with equal, entropy, CRITIC and PCA weights) and measures agreement with Kendall’s W (0 = none, 1 = identical rankings).

rc <- nexus_rank_compare(ind, benefit = benefit,
                         treatment_names = d$treatment)
rc
#> Ranking robustness across weighting methods
#> -------------------------------------------------- 
#>        treatment index_equal topsis_equal index_entropy topsis_entropy
#>  PB+SSDI+Residue           1            1             1              1
#>          PB+SSDI           2            2             3              2
#>          ZT+Drip           3            4             2              4
#>  ZT+Drip+Residue           4            3             4              3
#>     CT+Sprinkler           5            5             5              5
#>         CT+Flood           6            6             6              6
#>  index_critic topsis_critic index_pca topsis_pca mean_rank consensus
#>             4             1         1          1     1.375         1
#>             6             2         2          2     2.625         2
#>             1             4         3          4     3.125         3
#>             5             3         4          3     3.625         4
#>             2             5         5          5     4.625         5
#>             3             6         6          6     5.625         6
#> 
#> Kendall's W = 0.641 (moderate agreement across 8 rankings).
#> The rankings share the same data, so read W descriptively.
plot(rc)

Synergies and trade-offs

nexus_tradeoff() returns the correlation matrix with p-values attached; plot_tradeoff() draws it. Carbon intensity is lower-is-better, so pass benefit to orient the indicators: blue then always means the two improve together (synergy) and red means one improves as the other gets worse (trade-off).

plot_tradeoff(nexus_tradeoff(ind), benefit = benefit)

7. One-call workflow for replicated trials

arid_pulse_trial holds plot-level data from a randomized complete block design with three replications. nexus_analysis() computes plot-level indicators, runs ANOVA, builds the nexus index from treatment means and ranks treatments by TOPSIS.

data(arid_pulse_trial)
res <- nexus_analysis(
  arid_pulse_trial,
  treatment = "treatment", block = "replication",
  yield = "grain_yield", water = "total_water",
  energy_in = "energy_input", energy_out = "energy_output",
  emission = "ghg_emission", n_uptake = "n_uptake",
  weights = "entropy", carbon_basis = "intensity"
)
res
#> WEFNC Nexus Analysis
#> ================================================== 
#> Weights: Water = 0.181, Energy = 0.178, Food = 0.176, Nutrient = 0.308, Carbon = 0.157 
#> 
#>        treatment nexus_index nexus_rank topsis_closeness topsis_rank
#>          PB+SSDI      0.6852          1           0.8301           1
#>  PB+SSDI+Residue      0.6827          2           0.8252           2
#>          ZT+Drip      0.5758          3           0.6751           4
#>  ZT+Drip+Residue      0.5273          4           0.6883           3
#>     CT+Sprinkler      0.5055          5           0.3452           5
#>         CT+Flood      0.2117          6           0.1243           6
#> 
#> Components: $nexus, $anova, $topsis, $plot_data, $treatment_means

The ANOVA table is laid out as in most agronomy journals: means followed by letter groups (Fisher’s LSD), then SEm, CD at 5 percent and CV.

res$anova
#> Treatment comparison (means followed by the same letter do not differ significantly)
#> 
#>  treatment       grain_yield WUE    EUE    NUE      CI    
#>  CT+Flood        979.67 e    2.11 f 3.50 d 30.64 ab 1.84 a
#>  CT+Sprinkler    1121.00 d   2.77 e 4.77 c 32.24 a  1.42 b
#>  ZT+Drip         1288.33 c   3.73 d 6.73 b 29.20 bc 0.93 c
#>  ZT+Drip+Residue 1351.00 bc  3.97 c 6.62 b 27.86 c  0.91 c
#>  PB+SSDI         1420.33 b   4.44 b 8.22 a 28.15 bc 0.70 d
#>  PB+SSDI+Residue 1512.00 a   4.80 a 8.09 a 27.35 c  0.74 d
#>  SEm +/-         22.02       0.06   0.10   0.82     0.03  
#>  CD (p = 0.05)   69.37       0.20   0.32   2.60     0.09  
#>  CV (%)          2.98        3.03   2.81   4.88     4.50

nexus_anova() can also be used on its own for any set of response variables, in a CRD (no block) or RCBD.

Figures for the manuscript

Every result has a plot() method that returns a ggplot2 object, so you can add your own theme or labels and save it with ggplot2::ggsave().

plot(res$anova, response = c("grain_yield", "WUE", "EUE", "CI"),
     labels = c(grain_yield = "Grain yield (kg/ha)",
                WUE = "WUE (kg/ha/mm)",
                EUE = "Energy use efficiency",
                CI = "Carbon intensity (kg CO2-eq/kg)"),
     ncol = 2)

plot(res)                    # one radar per treatment

plot(res, which = "topsis")

References