BACKGROUND: As of 2021, the National Kala-azar Elimination Programme (NKAEP) in India has achieved visceral leishmaniasis (VL) elimination (<1 case / 10,000 population/year per block) in 625 of the 633 endemic blocks (subdistricts) in four states. The programme needs to sustain this achievement and target interventions in the remaining blocks to achieve the WHO 2030 target of VL elimination as a public health problem. An effective tool to analyse programme data and predict/ forecast the spatial and temporal trends of VL incidence, elimination threshold, and risk of resurgence will be of use to the programme management at this juncture. METHODOLOGY/PRINCIPAL FINDINGS: We employed spatiotemporal models incorporating environment, climatic and demographic factors as covariates to describe monthly VL cases for 8-years (2013-2020) in 491 and 27 endemic and non-endemic blocks of Bihar and Jharkhand states. We fitted 37 models of spatial, temporal, and spatiotemporal interaction random effects with covariates to monthly VL cases for 6-years (2013-2018, training data) using Bayesian inference via Integrated Nested Laplace Approximation (INLA) approach. The best-fitting model was selected based on deviance information criterion (DIC) and Watanabe-Akaike Information Criterion (WAIC) and was validated with monthly cases for 2019-2020 (test data). The model could describe observed spatial and temporal patterns of VL incidence in the two states having widely differing incidence trajectories, with >93% and 99% coverage probability (proportion of observations falling inside 95% Bayesian credible interval for the predicted number of VL cases per month) during the training and testing periods. PIT (probability integral transform) histograms confirmed consistency between prediction and observation for the test period. Forecasting for 2021-2023 showed that the annual VL incidence is likely to exceed elimination threshold in 16-18 blocks in 4 districts of Jharkhand and 33-38 blocks in 10 districts of Bihar. The risk of VL in non-endemic neighbouring blocks of both Bihar and Jharkhand are less than 0.5 during the training and test periods, and for 2021-2023, the probability that the risk greater than 1 is negligible (P<0.1). Fitted model showed that VL occurrence was positively associated with mean temperature, minimum temperature, enhanced vegetation index, precipitation, and isothermality, and negatively with maximum temperature, land surface temperature, soil moisture and population density. CONCLUSIONS/SIGNIFICANCE: The spatiotemporal model incorporating environmental, bioclimatic, and demographic factors demonstrated that the KAMIS database of the national programmme can be used for block level predictions of long-term spatial and temporal trends in VL incidence and risk of outbreak / resurgence in endemic and non-endemic settings. The database integrated with the modelling framework and a dashboard facility can facilitate such analysis and predictions. This could aid the programme to monitor progress of VL elimination at least one-year ahead, assess risk of resurgence or outbreak in post-elimination settings, and implement timely and targeted interventions or preventive measures so that the NKAEP meet the target of achieving elimination by 2030.