CO2 fertilization enhances vegetation productivity and reduces ecological drought in India

Higher warming will affect more regions globally with intensified agricultural and ecological droughts. Higher CO2 concentration improves vegetation’s water use efficiency (WUE), but its potential to alleviate extreme agricultural and ecological droughts is unclear. India is the second-highest contributor to global greening, having two of the eight global hottest biodiversity hotspots. Here, for the first time, using the CMIP6 earth system models (ESMs), we found an increase in the net vegetation productivity in India at the rate of 10.552 TgC year−1 with 1% per year increase in atmospheric CO2 concentration from 285 ppm to 1140 ppm, contrary to global trends. The improved WUE resulting from carbon fertilization and higher rain under warming will supersede the increased evapotranspiration water loss due to radiative effects. We found that the substantial increase in vegetation productivity in India attributes to plant physiology, and such factor needs to be considered in the drought projections.


Introduction
Human-induced climate change has resulted in increased warming, a rise in net radiation, and low relative humidity that, in turn, has increased the atmospheric evaporative demand (AED) and evapotranspiration (ET) globally (IPCC 2021). Increased AED causes more intense and frequent agricultural and ecological drought (Park Williams et al 2013, Allen et al 2015, Anderegg et al 2013, McDowell et al 2016, Douville et al 2021, Grossiord et al 2020. IPCC assessment report 6, WG1, chapter 5 (Canadell et al 2021), and chapter 11 (Seneviratne et al 2021) have assessed with high confidence that the global land carbon sink will be less efficient in a warmer world. Several studies have suggested that the enhanced AED will increase plant water stress under dry conditions (Brodribb et al 2020), resulting in plant mortality (Anderegg et al 2013, Park Williams et al 2013, McDowell et al 2016. Multiple regions around the globe have already experienced increased vegetation stress due to increased AED , Sanginés de Cárcer et al 2018, Yuan et al 2019. However, literature reports that a major driving factor controlling vegetation productivity is global water availability (Christian et al 2010, Jung et al 2017, Humphrey et al 2018 in terms of soil moisture (SM) and total water storage. Lower water storage results in a reduced land carbon uptake and carbon growth rate (Humphrey et al 2018). A few studies also suggest the increasing importance of vapor pressure deficit (VPD) in stomatal regulation (Novick et al 2016, Grossiord et al 2020).
Observational evidence suggests that high atmospheric CO 2 concentration benefits the plant by increasing vegetation productivity and water use efficiency (WUE) due CO 2 fertilization effect (Keenan et al 2013, Cheng et al 2017, El Masri et al 2019. Higher atmospheric CO 2 concentration reduces stomatal conductance with the increase of concentration gradient between the atmosphere and plant leaves (Swann et al 2016). Hence, with constant leaf area, the physiological responses of plants reduce water stress (Swann et al 2016), which is opposite to the impacts of increasing AED. However, most drought studies (Dai 2013, Prudhomme et al 2014, Cook et al 2015, Huang et al 2016, Pokhrel et al 2021, Satoh et al 2022 do not consider the increasing impact of atmospheric CO 2 on plant water intake. The scientific question is-can the reduced stomatal conductance balance the water stress arising from increasing warming and AED?-Remains unaddressed globally and regionally. A water deficit that impacts the ecosystem function and services and triggers feedback in natural or human climate systems is called ecological drought (Crausbay et al 2017). The changes in agricultural or ecological aridity/drought can be measured by changes in vegetation productivity and WUE in changing climate (Roderick et al 2015, Greve et al 2017). WUE, the ratio of gross primary productivity (GPP) to transpiration (TRAN), is a critical physiological indicator to quantify the rate of carbon uptake per unit of water loss by plants in terms of TRAN (Keenan et al 2013, Tang et al 2014 Betts et al 2007). On the other hand, leaf biomass increases with rising CO 2 , generating a larger evaporative surface that can partly offset the soil water saving caused by stomata closure (Lemordant et al 2016, Lemordant andGentine 2019). The trade-off between this increased WUE and augmented water uptake caused by increased biomass plays a fundamental role in structuring the plant response to ecological drought, emphasizing the need for understanding the impacts of rising CO 2 levels on ecological drought (Brodribb et al 2020).
The ratio of net carbon uptake (net primary productivity (NPP)) to the gross carbon uptake (GPP) is referred to as carbon use efficiency (CUE; de Lucia et al 2007). The CUE plays a vital role in maintaining the carbon balance on the land (de Lucia et al 2007) and determines the carbon allocated to biomass (El Masri et al 2019). The changes in plant productivity depend on multiple interconnected and conflicting factors, the rising CO 2 , the AED, the increased leaf biomass, and the water storage variability. However, understanding their individual contributions to plant productivity variations is unclear, which makes the projections of plant productivity uncertain in the near future under low warming (below 2 • C) (Canadell et al 2021, Douville et al 2021, Seneviratne et al 2021. India is the second-highest contributor to global greening (Chen et al 2019), and hence, the vegetation productivity in India will likely play a significant role in CO 2 uptake. However, the number of analyses is limited in this direction. A recent study shows that the climate drivers of Indian vegetation productivity differ from those of global tropical regions (Verma et al 2022). The projected changes in future drought in India have low confidence in the IPCC assessments due to limited analyses and uncertain monsoon simulations by the models (Seneviratne et al 2021 a). To the best of our knowledge, no study for India has recognized the climate implications on plant physiological effects of CO 2 .
In summary, there is a significant gap in understanding the trajectory of vegetation productivity and ecological drought in India, affected by covarying and conflicting factors such as CO 2 concentration, plant physiology, warming, increased AED, and uncertain monsoon. Here, we used the Coupled Model Intercomparison Project, phase 6 (CMIP6) (Eyring et al 2016) simulations, including the biogeochemical (BGC), radiative (RAD), and combined (FULL) experiments. Based on multimodel simulations, we projected the plant productivity in India under increased CO 2 concentration. Using different simulations, we have also separated the BGC effects from the atmospheric RAD effects of rising CO 2 on Indian vegetation productivity.

Data and methods
We evaluated the model performance against reanalysis data for the historical period (a historical experiment in CMIP6 models) for the period 1982-2014. The reanalysis data for hydro-meteorological variables is obtained from ERA5 (Hersbach et al 2019; 10.24381/cds.f17050d7). The GPP data is obtained from Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC) (Madani and Parazoo 2020). The reanalysis GPP is based on the Monteith light use efficiency (LUE) equation and was improved with optimized LUE values derived by FLUXNET tower sites. The gridded reanalysis GPP is then derived using the optimized LUE, Global Inventory Modeling and Mapping Studies (GIMMS3g), canopy fraction of photosynthetically active radiation, and Modern-Era Retrospective analysis for Research and Applications, Version 2, (MERRA-2) meteorological information (Madani and Parazoo 2020).
We selected six earth system model (ESM) outputs from the CMIP6 (Eyring et al 2016) archive, all selected variables with CO 2 forcing experiments. In the FULL experiments, there is a 1% increase in CO 2 concentration per year from a pre-industrial value of about 285 ppm (1st year; 1850). The CO 2 concentration in all experiments is increased from the global annual mean 1850 value until quadrupling (140th year; 1989). RAD experiments have the same experimental design, but only the RAD model of the ESMs has an increase in CO 2 concentration, and the BGC component has pre-industrial CO 2 concentration. In contrast, in the BGC experiment, only the BGC component of the ESMs has an increase in CO 2 concentration, and the RAD part has a constant preindustrial CO 2 concentration. In the abrupt-4 × CO 2 , an instantaneous quadrupling of the atmospheric CO 2 concentration from the pre-industrial mean value is imposed, then held fixed. The other forcing (e.g. anthropogenic aerosol forcings and land-use land-cover changes) have fixed pre-industrial ). Representation of the terrestrial and vegetation process varies among the models. ACCESS-ESM1-5, CMCC-ESM2, and MPI-ESM1-2-LR models explicitly use the nitrogen cycle to simulate the terrestrial process, which is essential for vegetation processes because of the potential limitation of nutrient availability during photosynthesis (Boé 2021).
All the standardization of variables is done with respect to piControl. We calculated the multimodel mean maps after regridding all the model outputs to a consistent 1.5 • × 1.5 • grid. The ecosystem WUE is calculated as the ratio of GPP to TRAN. The CUE is calculated as the ratio of GPP to NPP. The VPD is calculated as the difference between saturation and actual vapor pressure. Changes in the climatology for all variables are calculated by dividing the data of 120 years (11th−130th; 1860-1979) into equal 40 years span (i.e. 11th-50th, 51st-90th, and 91st-130th). The spatio-temporal changes in the variable are calculated as where X fut is the mean of X over the years 89th-118th while X hist is the mean of X over the years 1st-20th. Periods of averaging are defined so that the mean CO 2 concentrations of the hist and fut are close to 313 ppm and 800 ppm, respectively (Lemordant et al 2018).

Physiological and RAD effects on the eco-hydrological variables over India
We have assessed the changes in the critical ecohydrological variables over Indian land as simulated by C4MIP (Jones et al 2016) experiments, FULL, BGC, and RAD (details of experiments and ESMs in methods). Before understanding the trajectory of eco-hydrological variables as simulated by the ESMs, we first evaluated supplementary figure 1) their historical simulations against the reanalysis data (details in methods). We found that all the models performed well in simulating the TAS on the Indian landmass (supplementary figures 1-3). Though biases exist, the models are doing well in simulating VPD, PR, ET, and SM climatologies (supplementary figure 2). Most models, except for CanESM5 and MPI-ESM2-LR, simulate the climatology of GPP but with biases. The CanESM5 and MPI-ESM2-LR failed to simulate the climatology of GPP with a flat line across months. The other models show high GPP as expected during the retreat and end of the monsoon season. The multimodel average captured the climatology of GPP. We have also plotted the Taylor diagram (supplementary figure 3) to compare the model simulations with reanalysis data in simulating the mean annual values across each grid point. It shows how the model can capture the spatial variability of different processes. In general, models can capture the spatial variability of temperature and correlate highly with reanalysis data. However, the model performance in simulating PR, GPP, ET, and SM is comparatively low, which induces uncertainty, as observed in the results. However, there is still much scope to improve the land surface models (LSMs) used in ESMs, such as improving the acclimation parameter in models using photosynthesis parameters (Mengoli et al 2022). Figure 1 shows the trajectory of the eco-hydroclimatological variables as the function of CO 2 concentration obtained from the FULL simulations. We also used the piControl (pre-industrial CO 2 concentration) and abrupt-4 × CO 2 simulations. The trajectories were obtained from the FULL simulations standardized with respect to the piControl simulated values. The circles at the beginning and end of each trajectory show the mean of piControl and abrupt-4 × CO 2 simulations, respectively. It should be noted that outputs for all the variables were unavailable for all GCMs. Figure 1(a) shows the trajectories of TAS that show consistent and apparent increases with the rise in CO 2 concentration. However, the rise in temperature will not necessarily result in increasing VPD ( figure 1(b)). For example, CanESM5 shows decreasing VPD, despite the increasing temperature. This could be associated with the increased moisture supply from the ocean to the Indian landmass. Model simulations project a consistent increase in PR over the Indian landmass (figure 1(c)), getting reflected in the rise in SM (figure 1(e)) and Q (figure 1(f)), despite increases in ET ( figure 1(d)).
The trajectory of the LAI in most models shows the increase of vegetation biomass in India (figure 1(g)), which is consistent with the recently observed greening trend (Chen et al 2019). Model  figure 1). All the models show an increase in the GPP, with CMCC-ESM2 and CanESM5 being the two outliers on either side ( figure 1(i)). The NPP also shows a similar positive trend (figure 1(j)), with a multimodel average of 10.552 TgC year −1 (table 1). The CanESM5 and CMCC-ESM2 project has the highest (24.242 TgC year −1 ) and the lowest increase (1.936 TgC year −1 ), respectively. An increase in the TAS and LAI will result in a rise in TRAN at lower increasing levels of CO 2 . In contrast, TRAN stabilizes at a high CO 2 concentration (>600 ppm), probably because of the changes in the stomatal conductivity. This result is consistent across the models ( figure 1(h)). A high increase in GPP and NPP, with a lower increase (stable at a high CO 2 level) in TRAN, will increase the WUE (figure 1(l)). GPP and NPP show a high increase, but the increase in NPP is lower than that of GPP for a few models due to an increase in respiration. Hence, the projected changing patterns of CUE show high uncertainty from decreasing to increasing patterns.
The FULL simulations indicate that the increased productivity and WUE attribute to the changing vegetation physiological factors. To ascertain the same, we presented the trajectories of all the variables shown in figure 1 for the two simulations, RAD (supplementary figure 4) and BGC (supplementary figure 5). As expected, in the RAD simulations, the changes are observed in the hydro-meteorological variables (left column). In contrast, the vegetation variables remain constant with a few exceptions (right column). Temperature and PR in India show an increasing trend (supplementary figures 4(a) and (c)). Increasing PR improves SM, with an insignificant impact from increasing temperature and ET. The vegetation variables remain almost constant with a slight reduction in NPP and CUE (supplementary figures 4(j) and (k)). As expected, the WUE remains steady in the RAD simulations (supplementary figure 4(l)). The BGC simulations show no trend in the hydro-meteorological variables (supplementary figure 5, left column) but a similar trend in the vegetation variables (supplementary figure 5, right column) as observed in FULL simulations. The vegetation physiology increases the LAI, GPP, WUE, and NPP at a high atmospheric CO 2 concentration. The same is observed in FULL simulations and indicates that vegetation physiology governs productivity in India compared to the increasing temperature. This is in contrast to the global results and IPCC AR6 global assessments. Interestingly, all the models consistently show an increase in CUE in the BGC experiment, whereas FULL simulations show high uncertainty. As the temperature is not increasing in BGC simulations, respiration does not increase. However, the GPP increases due to carbon fertilization, which increases the CUE.
Supplementary figure 7 shows the CO 2 fertilization effect for different agroecological zones of India. The 20 agroecological zones are defined by the national bureau of soil survey and land use planning (Sehgal et al 1990). The agroecological zones are classified based on soil, climate, physiographic, and conducive moisture availability periods or length of the growing period. These agroecological zones are grouped into five ecosystems based on the climate zones (arid, semi-arid, sub-humid, humid, and coastal; supplementary figure 6). The results for all agroecological zones (supplementary figure 7) are similar to India's average (figures 1(i) and (l)). GPP and WUE overcome the negative impact of warming (supplementary figure 8) and increase in all the regions due to changes in the biogeochemistry of plants (supplementary figure 9), attributing to CO 2 fertilization effects.
The trend values for all the variables in each experiment are given in supplementary table 2. To understand the changes in eco-hydro-meteorological variables with respect to baseline (piControl). We have also plotted climatology of each variable for FULL, RAD, BGC and piControl for multimodel mean (supplementary figures 10-14). All the ecological variables (GPP, NPP) increase with respect to baseline in BGC and thus increase in FULL (supplementary figure 10, whereas hydro-meteorological variables (PR, SM, VPD, TAS, ET) increase in RAD and remain unchanged in BGC, thus increase in FULL (supplementary figures 11-13). However, TRAN increases in RAD but reduces in BGC after around 600 ppm (supplementary figure 13). This may be because of the water saving by plants at high CO 2 concentrations. The multimodel mean shows an overall increase in CUE and WUE for the FULL experiment The spatial distribution of changes (see details in methods) in GPP and NPP with changing CO 2 shows a positive increase in the FULL simulations (figures 2(a) and (d)) because of dominating BGC effects over vegetation for the whole of India (figures 2(c) and (f)). The RAD simulations show a slight decrease in NPP and GPP over India due to warming. The plant's physiological responses are not visible in these simulations. Hence, the net GPP and NPP are driven by plant physiology in India. The projected LAI changes with increasing CO 2 concentration, as simulated by FULL, RAD, and BGC simulations, follow the same pattern as observed for GPP and NPP.
The spatial patterns of changes in TRAN are different (figures 2 (j), (k), and (l)). There is a very clear west-east asymmetry. The TRAN shows an increase in the dry western regions even in the FULL simulations and a decrease in the wet eastern side. The northeast is a dominant forest region resulting in an all-India average stable TRAN. However, an increase in TRAN over the dry west region and a decrease in TRAN over the wet eastern region is not desirable for overall drought scenarios. The RAD effect and temperature impacts dominate the high TRAN over the western region. It is noteworthy that increasing TRAN does  Figure 3 presents the monthly changes in WUE and CUE with the increase in atmospheric CO 2 concentration as obtained from FULL simulations by the GCMs. We present the results divided into three periods 11th-50th, 51st-90th, and 91st-120th year of the simulation, with a uniform increase in the CO 2 concentration. WUE is a representative variable for water-carbon cycle interactions. Figures 3(a)-(d) shows the changing WUE. Except for the CanESM5, all other models show an increase in the WUE, Uncertainties exist in the projected CUE under increasing CO 2 (figures 3(e)-(j)) in the FULL simulations. ACCESS-ESM1-5, CMCC-ESM2, and MPI-ESM2-LR show a decrease in CUE, the first one uniformly over the entire year, while the latter only in the monsoon months. IPSL-CM6A-LR shows a slight increase, and CanESM5 shows a very high increase in CUE in the non-monsoon months. Interestingly, the RAD simulations for all the models project decreases in CUE, with the highest changes for CanESM5 (supplementary figures 16(e)-(j)). The sharp increase in the TAS causes the reduction in CUE in the RAD experiment. It may be because of the lower nutrient deposition at higher CO 2 concentrations, higher temperature, and high respiratory cost. The BGC simulations show exactly the opposite signs of changes (supplementary figures 17(e)-(j)) due to CO 2 fertilization. The uncertainties observed (figures 3(e)-(j)) in the FULL simulations result from the differences in magnitudes of positive and negative changes from RAD and BGC simulations by individual GCMs. Supplementary figure 18 shows the spatial distribution of the changes in the CUE and WUE over India. CUE change is similar in all regions except for a small region in the west. However, the increase in WUE in the western region is less than that of the eastern and southern peninsula because the RAD effects counteract the physiological effects on the dry regions of west India.

Resilience of vegetation productivity to meteorological drought
Literature suggests that PR, SM, and total water storage control the interannual variability of vegetation productivity (Humphrey et al 2018(Humphrey et al , 2021. Analysis of Indian vegetation shows PR to be the driving factor behind the variability of vegetation productivity (Verma et al 2022). Here, we investigated if carbon fertilization can reduce the adverse impacts of water stress on vegetation productivity. The possible pathway to increase ecosystem resilience could be increasing the WUE even under water-limited conditions (Sharma 2017). Figure 4 presents the variability of standardized GPP with the meteorological drought indicator, standardized precipitation index −12 (SPI-12), and standardized SM over India, as obtained from the FULL simulations. The scatter plots are presented for piControl, abrupt 4 × CO 2, and a 1% per year increase in CO 2 concentration. For the latter case, the results are separately presented for the periods 11th-50th, 51st-90th, and 91st-130th year of simulation.
We found a consistent increase in the standardized GPP with increased atmospheric CO 2 concentration. The impacts of carbon fertilization offset the adverse effect of water stress and improve GPP multiple times. Looking at scatter plots for individual cases, such as pi-control, we found that for all the GCMs, except for CMCC-ESM2, the total variability of GPP is low with the changing SPI-12, or standardized SM. For most cases, where the drought indicators vary from +2 to −2, signifying extreme wet or dry conditions, the standardized GPP values are near zero. These results indicate a high resilience of Indian vegetation to meteorological stresses. However, vegetation reacts significantly to carbon fertilization, with more than ten-fold increases in the 4 × CO 2 simulations. Hence, all the GCMs show distinctly identifiable scatter plots for different levels of atmospheric CO 2 concentrations. Supplementary figure 9 shows similar scatter plots for the RAD simulations. The GCMs either show decreasing GPP or almost unchanged GPP with increasing CO 2 , most likely because of temperature effects. For the GCMs ACCESS-ESM1-5, BCC-CSM2-MR, and CanESM5, the spreads of GPP across periods are more than SPI-12 (or standardized SM) within a single period. We can conclude that the GPP is more resilient to drought than to high temperatures when physiological effects are not considered. Supplementary figure 19 shows no change to decrease in SM drought expressed as standardized SM under increased CO 2 concentration. Such a change is attributed to the increased monsoon PR under warming in India, though there is low confidence in such projections due to poor model performance (Seneviratne et al 2021). We found exactly opposite results for BGC simulations (supplementary figure 20), which are consistent with FULL simulations. This clearly shows that the physiological vegetation factors are more influential than the radiation impacts resulting in net future productivity improvements under increased CO 2 concentration.
However, uncertainty exists among model simulations, for example, CanESM5 and BCC-CSM2-MR (figure 1). The CanESM5 shows a very low increase in WUE despite a high increase in GPP. This is because of a very high TRAN resulting from high PR and SM. On the contrary, BCC-CSM2-MR simulates a very high rise in WUE, leading from a moderate increase in GPP but a decrease in TRAN above a CO 2 concentration of 600 ppm. The differences across ESM simulations of plant variables can be attributed to different parameterization schemes across LSMs and the uncertainty across the meteorological variables. Supplementary table 1 shows that the land component of CanESM5 uses only nine PFTs, which differs from most models. CanESM5 also does not consider the Nitrogen cycle. These could be the reason behind a distinctly different simulation by CanESM5 from other ESMs.
Similarly, in figure 4, the impact of carbon fertilization on the GPP during drought is highest in CanESM5 and lowest in CMCC-ESM2. This is consistent with the results presented in figure 1(i). The latter uses 15 PFTs, whereas the former uses only nine PFTs. CMCC-ESM2 also considers fire and Nitrogen cycle, not considered in CanESM5. Such differences in the model contribute to higher differences in the changes in GPP under increased atmospheric CO 2 concentration. Despite high uncertainty, all the models agree to project an increase in GPP because of carbon fertilization under drought, which increases the confidence in the conclusion derived from the present study.
To understand and summarize the results, we presented the probability distribution functions (pdfs) of standardized GPP for drought scenarios (when SPI-12 < −1) for the 11th-50th, 51st-90th, and 91st-130th year of simulation (figure 5). For the initial period 11th-50th year, the GPP also shows ecological drought (standardized GPP < 0) scenarios with different probabilities, varying across GCMs. However, the ecological drought scenarios vanish for all the GCMs (BGC simulations) with increased CO 2, even under meteorological drought conditions. The periods 51st-90th, and 91st-130th year have higher atmospheric CO 2 concentrations than 11th-50th year, with the reduction of ecological drought and improved vegetation productivity.

Conclusions
There is high confidence (Canadell et al 2021) that the global land carbon sink will become less efficient at a higher emission scenario due to SM limitations. We found a very high increase in the carbon sink potential of Indian vegetation. The improved vegetation productivity in India under an increased CO 2 concentration attributes to the physiological changes in vegetation due to carbon fertilization. The other reason could be that CMIP6 models project increased monsoon PR and improved SM. Hence, there is a lack of increased water stress or SM limitation in the increased CO 2 concentration in India, leading to an increase in productivity. However, the IPCC has assessed low confidence in the future drought projections in India primarily due to poor model performance. Hence, directly from the simulations, the question remains unaddressed, if the future rainfall declines, will it impact the productivity of the 2nd highest contributor to global greening? We addressed the same with the results presented in figures 4 and 5. The vegetation productivity in India is not very sensitive to variations of SPI-12 or standardized SM. Instead, they react significantly to the changing atmospheric CO 2 concentrations, showing high productivity due to CO 2 fertilization. The comparison between the RAD and BGC simulations shows that the vegetation's physiological factors and biogeochemistry will play a dominant role compared to the RAD effects resulting in very high net improvements in productivity.
There are two conflicting factors related to vegetation affecting SM or land water availability. First, reduction of stomatal conductance with increasing WUE resulting in low water intake and, second, an increase in the LAI resulting in high water intake. We found that the physiological factor resulting in reduced stomatal conductance dominates the process in Indian vegetation, with an overall improvement in the drought scenarios. Notably, the future increase in LAI will also depend on human land management practices and the associated government policies. The CMIP6 models do not explicitly consider such region-specific information, which is a limitation of the present study. For example, LSMs used in the ESMs of CMIP6 have different irrigation modules. The methodology of irrigation input in the models is empirical and does not change temporally nor reflect varying human irrigation practices. These factors can also impact the GPP and WUE by affecting land-atmosphere feedback. Studying these impacts will require regional climate models coupled with LSMs. Also, a decrease in ecological drought in the future under warming due to increased carbon fertilization and monsoon rainfall may have significant implications on climate adaptation as well. For example, increased plant productivity and reduced water use may have the potential to ensure water and food security associated with sustainable development goals 2 and 6. However, projecting the changes in land management practices needs a more detailed holistic regional analysis, which may include changing fertilizer use in agriculture and subsequently reducing pollution, switching of crops, and intensification of agriculture. These may be considered as the potential future research area of the present analysis.
Another limitation of this study could be not considering the response of C3 and C4 plants separately. C3 plants are more responsive to the CO 2 fertilization effect than the C4 plants because of the difference in photosynthesis mechanism (Boretti andFlorentine 2019, Wang et al 2020). However, few studies have also shown that these responses could be reversed in long-term elevated CO 2 concentrations (Reich et al 2018). ESMs consider C3 and C4 plant mechanisms and their response to changing forcings separately (Boysen et al 2020). Depending on the fractions of different types of plants within a grid, the models produce the gridded average GPP. The available simulations do not provide the differential responses of C3 and C4 plants to CO 2 concentrations; hence, we could not separately produce the carbon fertilization impacts on C3 and C4 plants for India (supplementary text S2). Further, the analysis of daily simulations could have provided more robust insights into the effects of elevated CO 2 on the WUE and other environmental variables. However, we have used monthly temporal resolution due to the non-availability of daily data for BGC and RAD simulations.
Few studies also have pointed out that the vegetation productivity increase due to CO 2 fertilization may get saturated in future warming (Shi et al 2021). Extreme warming and low water availability can offset the positive impact of CO 2 fertilization by drying out the plants and increasing plant mortality due to xylem embolism, forest fire, insect and invasive plant species attacks due to climate change (Choat et al 2018, Brodribb et al 2020, Grossiord et al 2020. However, we found that the GPP continuously increases throughout this period, and we do not get a saturation point. This may be one of the limitations of the experiments used in the present study. Overall, there is high consistency across the models that India will experience a reduced ecological drought with increased productivity under increased CO 2 concentration. We found that the plant can sustain better even in drought scenarios in the future, compared to the present, due to carbon fertilization. India has committed to Net Zero by 2070, and planning for investments is in process for the same. Such planning involves afforestation, carbon capture and sequestration technology development, renewable energy development, etc. India is the 2nd highest contributor to global greening (Chen et al 2019). Priorities on afforestation depend on the trajectory of vegetation carbon uptake potential in the future under increased atmospheric CO 2 concentration and warming. The trajectory depends on meteorological factors and carbon fertilization. While the former has been addressed in the literature, the latter has not been considered for analyzing the ecological impacts of climate change. The present study concludes the higher potential of vegetation in India under increased GHG concentration and recommends the continuation of the existing high priority of afforestation by the Government of India. Our results indicate that the increased vegetation productivity will contribute to India's climate mitigation pledge and net zero plan by 2070.

Data availability statement
No new data were created or analysed in this study.