Changes in glacier terminus position have been implicated as one of the primary drivers of the rapid changes in glacier dynamics observed across the globe in the last two decades. Iceberg calving exerts a critical control on the terminus position of the vast majority of marine-terminating glaciers, yet calving is relatively poorly understood due to the inherent difficulties in collecting observations of a stochastic process in a dangerous setting. Time-lapse camera and satellite observations suggest that the style of iceberg calving can vary tremendously in both space and time depending on the physical properties of the terminus, ranging from the detachment of giant tabular icebergs every few decades from Antarctic’s floating ice shelves to the growlers produced nearly daily from serac topples along Alaska’s coast. Here we extract quantitative metrics on the relative importance of calving driven by branching and uncorrelated fractures through application of fragmentation theory to iceberg size distributions extracted from high-resolution digital elevation models for 17 fjords around Greenland. We find that iceberg size distributions typically deviate from the widely-assumed power-law form for icebergs with surface areas >0.05 km^2, with fewer icebergs than predicted by the power-law for larger sizes. Icebergs larger than ~0.1 km^2 primarily calve as the result of full-thickness penetration of uncorrelated fractures (i.e., as tabular icebergs). Although the dataset is temporally sparse for the majority of the study sites, the data suggest that iceberg formation via branching fractures reaches a seasonal peak in summer, when icebergs up to ~0.1 km^2 follow power-law distributions. These data provide a novel means to assess the accuracy of iceberg calving models and potentially to constrain the physical characteristics of termini susceptible to the marine ice cliff instability mechanism.

Global atmospheric ‘storm-resolving’ models with horizontal grid spacing of less than 5~km resolve deep cumulus convection and flow in complex terrain. They promise to be reference models that could be used to improve computationally affordable coarse-grid global climate models across a range of climates, reducing uncertainties in regional precipitation and temperature trends. Here, machine learning of nudging tendencies as functions of column state is used to correct the physical parameterization tendencies of temperature, humidity, and optionally winds, in a real-geography coarse-grid model (FV3GFS with a 200 km grid) to be closer to those of a 40-day reference simulation using X-SHiELD, a modified version of FV3GFS with a 3 km grid. Both simulations specify the same historical sea-surface temperature fields. This methodology builds on a prior study using a global observational analysis as the reference. The coarse-grid model without machine learning corrections has too little cloud, causing too much daytime heating of land surfaces that creates excessive surface latent heat flux and rainfall. This bias is avoided by learning downwelling radiative flux from the fine-grid model. The best configuration uses learned nudging tendencies for temperature and humidity but not winds. Neural nets slightly outperform random forests. Forecasts of 850 hPa temperature gain 18 hours of skill at 3-7 day leads and time-mean precipitation patterns are improved 30% by applying the ML correction. Adding machine-learned wind tendencies improves 500 hPa height skill for the first five days of forecasts but degrades time-mean upper tropospheric temperature and zonal wind patterns thereafter. The figure shows maps of 30-day time-mean precipitation pattern difference from the fine-grid reference for prognostic simulations: (a) 200 km baseline(no machine learning correction) (b) Using random forest correction and (c) neural net correction for temperature. humidity and surface radiation corrections. RMSE is the root mean squared precipitation difference from the reference, which is 30% less for the two machine-learning corrected simulations compared to the baseline. (d) Bar charts of the land-mean, ocean-mean and global-mean precipitation biases for these three configurations, showing the machine-learning corrected simulations remove a high bias of land surface precipitation in the baseline simulation.

Coastal aquifers act as a major host of seawater-groundwater interaction and play an important role in modulating the marine elemental budget. Calcium stable isotopes (δ44/40Ca) have been used to identify mass dependent isotope fractionation processes such as carbonate dissolution and precipitation in a range of geological reservoirs that has major implications in constraining global geochemical cycles. However, there is limited Ca isotope data from coastal aquifers globally. Here we report δ44/40Ca values of groundwater collected in 2017-18 from multiple locations and depths from the Bakkhali delta front, Sundarbans, India. The sampling depth varied between 14 m below ground level (m bgl) and 333 m bgl and the salinity of the groundwater samples range from 1-25 ppt. The salinity of the water samples decreases with increasing depth indicating greater seawater incursion from the Bay of Bengal at shallower depths. Variable amounts of mixing of freshwater and seawater is also supported by Sr and Ca concentrations which vary between 1.6-62.8 μmol/l and 0.29-8.92 mmol/l, respectively, and show progressively lower concentrations with depth. The δ44/40Ca values of dissolved phase in groundwater samples (relative to NIST SRM 915a) were measured using a 43Ca-48Ca Double Spike TIMS technique at Centre for Earth Sciences, Indian Institute Science, Bangalore. The δ44/40Ca values of the groundwater samples show significant variability between 1.52-2.28‰ (2SD ~0.1‰) with several samples showing δ44/40Ca values higher than modern seawater (~1.88‰). Low δ44/40Ca values, mostly in deeper groundwater samples, is consistent with higher proportions of freshwater input. Samples showing high δ44/40Ca values are mostly from shallower depths and likely reflect carbonate precipitation which is consistent with high Sr/Ca (~8.13-12.26) in samples from 30-42 m depth. The high δ44/40Ca in groundwater samples from the Ganges-Brahmaputra delta could explain the high δ44/40Ca values reported in water samples from the Bay of Bengal1. [1] Chakrabarti et al., Goldschmidt Boston, 2018 Abstract.

The largest city in New Zealand, Auckland is home to roughly 1.5 million people -- one third of New Zealand's population. Here we assemble a bottom-up inventory of Auckland's fossil fuel carbon dioxide emissions from a variety of data sources. We use these emissions estimates in combination with the UrbanVPRM land surface model to estimate the net carbon balance of the region. This work is part of the larger CarbonWatch NZ project, which aims to produce estimates of New Zealand's net carbon balance quickly enough to assess and refine ongoing national efforts to reach carbon neutrality.

Financial institutions’ investment and lending portfolios could be affected by both physical climate risks stemming from impacts related to increasing temperatures, and from transition climate risks stemming from the economic consequences of the shift to a low-carbon economy. Here we present a consistent framework to explore near term (to 2030) transition risks and longer term (to 2050) physical risks, globally and in specific regions, for a range of plausible greenhouse gas emissions and associated temperature pathways, spanning 1.5-4oC levels of long-term warming. We draw on a technology-rich, regionally disaggregated Integrated Assessment Model representing energy system, agricultural and land-based greenhouse gas emissions, a reduced complexity climate model to simulate probabilistic global temperature changes over the 21st century, and a suite of impacts models to estimate regional climate-related physical hazards and impacts deriving from the temperature change pathways and their underlying socio-economics. We consider 11 scenarios to explore the dependence of risks on both temperature pathways, as well as socio-economic, technology and policy choices. This builds and expands on existing exercises such as the Network for Greening the Financial System (NGFS). By 2050, physical risks deriving from major heatwaves, agricultural drought, heat stress and crop duration reductions depend greatly on the temperature pathway. By 2030, transition risks most sensitive to temperature pathways stem from economy-wide mitigation costs, carbon price increases, fossil fuel demand reductions and potential stranding of carbon-intensive assets such as coal-fired power stations. The more stringent the mitigation action, the higher the abatement costs and sector-specific transition risks. However, such scenarios result in lower physical climate hazards throughout the century. Our study also explores multiple 2 deg C pathways which demonstrate that scenarios with similar longer-term physical risks could have very different near-term transition risks depending on technological, policy and socio-economic factors. As such, “a single scenario will not answer all questions”.

Climate warming is projected to increase the frequency of climate extremes. The physiological response of tree species to these changes has not been well characterized. Particularly, how hydraulic parameters may adjust with temperature and what are the impacts during a severe drought in mesic forest tree species have limited study. Further, the pattern of responses to drought and recovery has been linked with isohydry-anisohydry strategies, but we do not know if these strategies are still relevant in the hotter and drier future. To fill this knowledge gap, we built a factorial experiment with ambient (AT) and high temperature (HT; ambient +4.5 0C) for two growing seasons, with six weeks of the initial acclimation period, four weeks of drought, and six weeks of recovery during the first year on seven tree species. Throughout, we followed water relations, leaf gas exchange, stem water potential, and hydraulic conductivity. Seedlings were acclimated to prolonged HT by adjusting key hydraulic traits and increasing anisohydry strategies, but these traits were not coordinated, as some isohydric species had surprisingly high stem water potentials at 50% loss of xylem conductivity (P50) and some anisohydric species had negative safety margins. In addition to this, anisohydric species reduced stem water potential during drought and delayed recovery, while isohydric species experienced earlier reductions in photosynthesis and increased dieback at HT. Our findings highlight that hydraulic traits may acclimate to HT and reduce some of the additive effects during drought in anisohydric species, but hotter drought intensified soil and atmospheric drought and increased mortality in isohydric species. Thus, acclimation of hydraulic traits to HT may help anisohydric species to avoid some of the consequences but can’t compensate for the negative effects.

Constraining the moment release of injection-induced earthquakes is of paramount importance to reduce the seismic hazard in the geo-energy industry. Recent studies suggest that a significant part of the moment release during fluid injections can be due to aseismic motion, namely, aseismic moment M0. Current models of injection-induced aseismic moment do not incorporate fault rupture mechanics. Here, we present a theoretical and numerical analysis that highlights a possible scaling relation between the aseismic moment and a key operational parameter, the injected volume of fluid V. The scaling relation emerges from the model of a stable frictional shear crack that propagates in mixed mode (II+III) on a planar fault interface. The interface is characterized by a constant hydraulic transmissivity and a shear strength that is equal to the product of a constant friction coefficient and the local effective normal stress. Fluid is injected right into the fault interface at a constant flow rate. The resulting relation between the aseismic moment and the injected volume is M0=A⋅ V^(3/2). The prefactor A accounts for the dependence of the aseismic moment on the pre-injection stress state, the parameters of the injection (notably, the injection flow rate), and the fault elasto-frictional and hydraulic properties. Unlike previous studies, our model accounts for the possibility that ruptures can propagate beyond the fluid-pressurized fault patch, a condition that is expected to occur in critically stressed and/or highly-pressurized fractures/faults. We test the scaling relation against estimates of moment release due to aseismic motion during fluid injections that vary in size from laboratory experiments to industrial applications. Our predictions are in good agreement with these observations. These results provide a simple means to quantify the size of aseismic ruptures in response to fluid injections related to both natural and human sources.

Radiance measurements from several types of passive microwave (PMW) sensors are combined in the Global Precipitation Measurement mission (GPM) to increase the temporal and spatial coverage of precipitation observations. The measurements of these sensors are converted to precipitation estimates by the GPM Profiling Algorithm (GPROF). High frequency PMW-channels are used to retrieve precipitation estimates over land, as these frequencies can measure the radiance scattered by ice particles in rain clouds. Scattering related to shallow and low-intensity precipitation events, however, is limited. Hence, radiometric signals associated with these events are hard to distinguish from the naturally emitted radiation from the Earth’s surface, especially since this so-called background radiation is dependent on the surface type. A better understanding of the physical processes that occur during precipitation events can help to identify possible weaknesses in the GPROF algorithm. Hence, this study couples overpasses of GPM radiometers over the Netherlands to two dual-polarization radars from the Royal Netherlands Meteorological Institute (KNMI) in 2019. All rainy overpasses (>0.1 mm/hr) within a 75 km radius around one of the radars are selected. This coupling provides the opportunity to relate GPROFs performance to physical characteristics of precipitation events, such as the vertical reflectivity profile and dual-polarization information on the melting layer. Additionally, simultaneous observations from both the PMW sensor and the dual-frequency precipitation radar (DPR, used as a-priori database in GPROF) aboard the GPM core satellite are available. Hence, space-based and ground-based reflectivity profiles can be compared and coupled to discrepancies of the GPROF algorithm.

Wet Prairies are unique communities occurring in isolated, saturated depressions that typically have standing water from late fall to late spring and then dry up and are burned by lightning or human-set fires in summer. Wet prairies are highly sensitive to precipitation and evapotranspiration patterns as they are highly dependent on shallow groundwater levels. Three years (May 2015-May 2018) of hourly watertable data recorded by a datalogger in a piezometer, and hourly precipitation, temperature, relative humidity, and wind speed data from the Toledo Express Airport (Ohio) weather station were used as inputs for a nonlinear autoregressive neural network with exogenous inputs (NARX) model. The data was used to train the model and to forecast watertable levels at an hourly and a daily time step. The NARX model used a Levenberg-Marquardt learning algorithm along with a combination of input delay, feedback delay, and hidden layers. The training set encompassed 70% of the total dataset, with cross-validation and testing covering the remaining 30% of the total dataset split equally among them. The model was trained with the hourly data and then the hourly data was aggregated to a daily record, and a new NARX model trained using this new time-step. The intention was to see if a much smaller data record, but typically all that is needed by a land manager, was capable of producing a satisfactory watertable forecast. The NARX model’s performance was assessed using R2, RMSE, and the Nash-Sutcliffe index. The NARX model provided successful short-term forecasts (6 months) for hourly and daily temporal resolutions. The R², RMSE, and Nash-Sutcliffe for the hourly testing period were 0.85, 0.08, and 0.85 respectively. For the daily testing period the R², RMSE, and Nash-Sutcliffe were 0.91, 0.07, and 90 respectively. The NARX model was not able to predict a sudden increase in water table levels due to a large snowmelt event in 2018, but this is not surprising as the model was not trained using snow melt events or a snow depth variable. Regardless of its current limitations, land managers could use this NARX model to better understand watertable patterns in wet prairies, one of the main drivers of this natural community.

Climate change will amplify erosion rates as erosive power of rainfall will likely increase due to change in rainfall characteristics (e.g., energy, intensity, duration as well as frequency). According to Intergovernmental Panel on Climate Change Fourth Assessment report (IPCC AR4), it is projected that by mid-21st Century, the rainfall across southeastern US will both increase and decrease in intensity, which will substantially affect rainfall erosivity. Few studies have estimated the impact of climate change (e.g., rainfall intensity) on rainfall erosivity across US and around the world. However, previously published erosion indices have discrepancies due to differences in methodologies (e.g., primarily omission of small and low rainfall intensity) adopted in those studies. Therefore, the objective of this study was to estimate change in erosion indices for the period 2030-2059 using the benchmark rainfall indices established for southeastern region of US. Hourly precipitation data were retrieved from NA-CORDEX under Representative concentration pathways (RCP) 8.5. Results on change in erosion indices as a result of climate change will be presented.

How can we use our current wealth of terrestrial data, encompassing biogenic and abiogenic systems, to determine the distinguishing properties of life? SCOBI (Statistical Classification of Biosignature Information) uses machine learning techniques to algorithmically identify combinations of measurements that are “indicative of life”. A set of ~1000 observations, comprising elemental abundance, isotopic fractionation, VNIR reflectance, and (in progress) Raman spectra, have been assembled from existing literature and databases. The observations cover systems classified as “indicative alive” (e.g., cells, vegetation), “indicative non-alive” (e.g., fossils, teeth), “mixed indicative” (e.g., soil, pond water), or “non-indicative” (e.g., rocks, meteorites). VNIR data was preprocessed by linear interpolation from 400-2100 nm and smoothed with a Savitzky-Golay filter. To limit the amount of Earth-biochemistry-specific (non-agnostic) information included, the first five spectral features extracted were number of peaks, number of troughs, mean reflectance, mean peak width, and broadest peak width. To help further emphasize agnostic biosignatures, Earth-specific features such as chlorophylls have been manually flagged so that feature importance with and without them can be compared. Classifiers including k-nearest neighbors (KNN), Gaussian Naïve Bayes (GNB), logistic regression (LR), random forest (RF), and support vector machine (SVM) were implemented, as was a combination voting classifier. Performance metrics included false positive rates, false negative rates, and AUC with 50-50 test/train splits (Monte Carlo simulations). Key takeaways from this stage, prior to the inclusion of Raman spectra, are (1) the overall success rate of 0.933 AUC was most heavily influenced by the elemental abundance data; and (2) VNIR reflectance had the lowest classification performance with 0.52 AUC (58% of objects correctly classified). The next steps are to complete integration of Raman spectral data and to improve the approach to pre-processing and feature extraction for both types of spectral data, such as automated baseline removal, whole spectrum matching, and dimensionality reduction.

​​We believe that the future of Earth Observation (EO) is in fusion, harmonization, and interoperability of satellite imagery. Intensified monitoring leads to better understanding of land use and reduction of maintenance costs for all Land Cover products. RapidAI4EO is an initiative that aims to establish the foundations for the next generation of Copernicus Land Monitoring Service (CLMS) products. The goal is to provide intensified monitoring of Land Use (LU), Land Cover (LC) changes at a much higher spatial resolution and temporal cadence than is possible today. Key objectives are to explore, evaluate, and quantify state of the art deep learning algorithms and methodologies that leverage three meter, daily time series, in conjunction with higher spectral resolution Sentinel-2 imagery.  Consortium Partners: The RapidAI4EO projects brings together Planet Labs PBC, the operator of the world’s largest fleet of Earth-imaging satellites and the recognized leader of the CubeSat revolution, VITO, the main production center of the Copernicus Global Land Service, Vision Impulse, a recent spin-off of German Research Center for Artificial Intelligence (DFKI, the largest research center for Artificial Intelligence in the world and one of the two European NVIDIA AI Labs), the International Institute for Applied Systems Analysis (IIASA) whose Center for Earth Observation and Citizen Science (EOCS) devises new approaches and technologies to collect data on land cover and land use, and Serco Italia, a worldwide service provider to governments, international agencies and industries, and operator of the ONDA DIAS platform. The objectives of RapidAI4EO are: 1) the creation and release of the most comprehensive spatiotemporal EO training sets ever produced for machine learning;  2) the development and implementation of novel AI solutions for continuous change detection that leverage these data sets; 3)the ability to drive frequent temporal updates of the Corine Land Cover (CLC) product; and 4)to demonstrate improved LULC mapping using harmonized Sentinel-2 and very high resolution, high cadence data streams.

The Dijkstra’s algorithm is employed to find the shortest paths between the epicenter and all the wet nodes on a global triangular unstructured mesh using the phase speed of surface gravity waves and progressive acoustic modes. The phase speed estimator takes into account the simultaneous effects of the slight compressibility of water, sea-bed elasticity, and static compression of the ocean under gravity, leading to precise calculation of the arrival time [1]. Coupled to a mathematical solution for multi-fault rupture along the transect extracted from the shortest path between the source and destination is then employed to characterize the genesis of the tsunami [2]. The effectiveness of the aforementioned algorithms and distinct arrival time of different the tsunami and its precursor components (acoustic modes) at DART and deep ocean acoustic observations can be used for tsunami early warning systems. [1] Abdolali, A., Kadri, U. & J.T. Kirby, 2019, Effect of Water Compressibility, Sea-floor Elasticity, and Field Gravitational Potential on Tsunami Phase Speed, Scientific Reports, Nature, doi:10.1038/s41598-019-52475-0 [2] Williams, B., Kadri, U., & Abdolali, A. (2021). Acoustic–gravity waves from multi-fault rupture. Journal of Fluid Mechanics, 915, A108. doi:10.1017/jfm.2021.101

Coastal wetlands provide a livelihood for local communities and simultaneously provide services that include coastal protection, enhanced carbon storage and habitats provision for both terrestrial- and marine life. Yet, the development of climate resilient management strategies that need to account for changes that might occur in the coming decades is challenging due to the variable response of coastal ecosystems. This is because coastal ecosystems, including mangroves, wetlands and tidal flats, are directly linked to coastline dynamics, as coastlines can move seaward, when waves are low and sediment availability is high, or retreat when the opposite is true. Particularly the long-term fate of the mud dominated coastline of Suriname, part of the Guianas coastal system stretching from the Amazon to the Orinoco delta, is determined by migrating subtidal mudbanks that cause a cyclic instability of alternating erosion and progradation phases. We present a semi-automatic remote sensing approach to quantify the influence of mudbank migration on coastline dynamics along the entire coast of Suriname. We validated our approach with high resolution drone imagery collected at contrasting locations, indicating average accuracy of changes in coastline position to be within 50 meters. This suggest we can apply our methodology on all available Landsat images between 1985 and 2020, acquired over the heterogonous and complex coastal landscapes. The results show that regional forcing mechanisms, related to migration of six to eight subtidal mudbanks in front of the Suriname coast, have a strong imprint on local coastline dynamics with an average 25 m/yr-1 expansion during mudbank presence, and 8 m/yr-1 retreat of the coastline during mudbank absence between 1986 and 2020. More importantly, we found that not all spatial and temporal variations in the magnitude and timing of local changes can be explained by migrating mudbanks. This demonstrates the importance of incorporating changes that cannot be explained by regional forcing mechanisms in management frameworks that aim at explaining current variability and predicting future coastline changes.

Floods are among the most destructive natural hazards on Earth. In paleohydrology, sediments are generally considered as one of the best archives to extend flood records to pre-historical timescales. Doing so requires being able to identify flood deposits from sediment archives and decipher between flood types. The latter is particularly important in glacierized regions, where meteorological floods frequently co-occur with Glacial Lake Outburst Floods (GLOFs). In Patagonia, GLOFs are particularly pronounced in the Baker River watershed (48°S), where 23 events occurred between 2008 and 2020. Since 1976, the same region experienced three intense rain-on-snow events. To identify the sedimentary signature of these flood events, ten sediment cores collected in the fjord immediately downstream of the Baker River (Martínez Channel) were investigated and compared to the recent flood history of the river. Results show that sediments accumulate on the fjord head delta at 2.0 to 3.4 cm yr-1 and that GLOF deposits can be distinguished from background sediments by their finer grain size (5.98 ± 0.82 μm) and lower organic carbon content (0.31 ± 0.06%), reflecting the release and transport in suspension of high amounts of glacial rock flour. Our results also show that the rain-on-snow events that occur in summer, and therefore primarily affect the glacierized part of the watershed, have the same sedimentary signature as GLOFs. In contrast, rain-on-snow events occurring in winter have a distinct coarse and organic-rich signature, reflecting sediment input from the non-glacierized part of the watershed. In summary, this study shows that (a) GLOF deposits in fjord sediments are distinct from typical flood turbidites and are best identified by their low grain size and total organic carbon content, and (b) the sedimentary signature of rain-on-snow floods in partially glacierized watersheds depends on the season during which they occur. We anticipate that our findings will contribute to a better interpretation of flood records from partially glacierized watersheds.

Assessing and understanding the extent and trajectory of change in inland waters is a great challenge, due in part to both differing methods — and cultures — of agencies that provide synoptic observations of Earth’s systems as well as the community of lake scientists whose research generates heterogeneous and distributed in situ data. Advancements require socio-technological initiatives that harness the resources of the highly diverse and distributed community of ecologists, as well as the products and expertise of the satellite remote sensing community. Here we describe a prototype for linking in situ and remotely sensed data for lakes through the collaborative efforts of the Global Lake Ecological Observatory Network (GLEON), the Environmental Data Initiative (EDI), and NASA. GLEON provides a community of lake scientists and data from lake observatories. EDI curates and publishes data and ensures conformity to rigorous FAIR principles. NASA provides the expertise and workflows to deliver remotely sensed data products on demand. The integration of the data and the communities provides a foundation for a new generation of lake science.

The Named Storm Event Model (NSEM) which includes the atmospheric suite (HWRF 1, HRRR 2, RTMA-URMA 3 and WRF-LES 4) [1], the coupled wave-surge-riverine modeling system (WW3 5-ADCIRC 6-NWM 7) [2,3] and a validation suite has been developed to provide definitive estimates of wind and water variables of major landfalling hurricanes, in compliance with the United States COASTAL Act of 2012. Within this framework, the performance of atmospheric products is evaluated and blended fields are generated on a high-resolution grid, using a master blend recipe. The atmospheric data forces the downstream hydro-coupled component which includes the wave (WAVEWATCH III [4,5]), ocean circulation (ADCIRC) and hydrological models (NWM). These have been coupled using the community-based National Unified Operational Prediction Capability (NUOPC) layer based on the Earth Systems Modeling Framework (ESMF). The wind and hydro products are evaluated against offshore and coastal wave buoys, nearshore water level stations, radars and satellite altimeters, as well as USGS rapid-deployment water level and wave gauges placed in the nearshore regions and overland. The NSEM workflow will be presented at the conference, highlighting the advantages of ESMF in model performance improvement and challenges for upstream atmospheric model blending, model coupling, and validation against observations. 1- HWRF (Hurricane Weather Research and Forecasting) 2- HRRR (High-Resolution Rapid Refresh) 3- RTMA-URMA (Real Time Mesoscale Analysis-UnRestricted Mesoscale Analysis) 4- WRF-LES (Weather Research and Forecasting-Large Eddy Simulation) 5- WAVEWATCH III (WAVE-height, WATer depth and Current Hindcasting) 6- ADCIRC (The ADvanced CIRCulation model) 7- NWM (National Water Model)

India is the largest user of groundwater in the world, accounting for a quarter of the world’s annual groundwater extraction. In Kerala, a state in Southern India, groundwater is the drinking water source for 62% of the population. Groundwater management is vital to ensure the sustainable availability of water. However, it is complex due to uncertainties in subsurface parameters. Groundwater modeling can be used as an efficient tool to understand the regional flow as well as the hydrogeological condition of the groundwater system. The modular finite-difference groundwater flow model (MODFLOW) developed by the U.S. Geological Survey (USGS) is a widely used modular three-dimensional block-centered finite difference code used for simulating dynamic features and scenarios in groundwater systems. In this study, the Visual MODFLOW Flex 7.0 is used to simulate the flow of groundwater through aquifers in the Thuthapuzha sub-basin, a sixth-order sub-basin of the Bharathapuzha river basin Kerala’s second-largest river basin. A four-layer model covering an area of 940 km2 has been prepared using Visual MODFLOW Flex 7.0. Groundwater level data of 37 wells spread over the study area has been used for the model simulation. Model design, calibration, and validation have been carried out. Boundary conditions including recharge and evapotranspiration are derived from the Soil and Water Assessment Tool (SWAT) hydrologic model. The model calibration is performed in two stages, which include steady-state calibration and transient state calibration using observed groundwater levels from 2000 - 2016. The validation is done by using observed groundwater levels from 2017 - 2019. The spatial distribution of hydraulic conductivity and storage properties are optimized using a combination of the trial-and-error method. Model performance evaluation is done using Standard Error of the Estimate, Root Mean Squared and Normalized RMS, found to be in the optimum range. The results of the model are found in good agreement with the observed groundwater level data. Therefore, the model is suitable for studying the groundwater level changes in the study area. In addition, the model can be extended to study the groundwater flow dynamics of the whole Bharathapuzha river basin and similar subbasins for the sustainable planning and management of groundwater resources.