20 years of trials and insights: bridging legacy and next generation in ParFlow and Land Surface Model Coupling

Groundwater plays a vital role in terrestrial water and energy cycles. Yet, it remains oversimplified in most Earth system models (ESMs), limiting their ability to represent key land-atmosphere interactions, including evapotranspiration partitioning, drought propagation, and boundary layer development. The original coupling of ParFlow with the Common Land Model (CoLM) in 2005 not only demonstrated the feasibility of integrating physically based groundwater models into ESMs, but also revealed emergent behaviors – such as lateral moisture redistribution, along with the buffering effects that emerge from enhanced subsurface connectivity – that cannot be captured by traditional land surface models (LSMs). This study reviews key findings from two decades of ParFlow–land/atmosphere coupled modeling efforts, highlighting how groundwater–land–atmosphere interactions shape surface energy balance and hydrologic connectivity across three dimensions: upward feedbacks, downward influences, and the critical zone of coupling. Given the substantial advances in LSMs such as CoLM over the past two decades, a renewed recoupling effort is warranted to enhance our understanding of groundwater's role across a broader range of Earth system processes. Preliminary efforts to recouple ParFlow with the updated water and energy modules of CoLM demonstrate improved performance when evaluated against reanalysis and observational data. To ensure long-term sustainability, we propose a modular and maintainable coupling framework addressing functional extensibility, data/code interoperability, and parallel computing needs, in which area, TerrSysMP2 has taken early steps and may be considered an initial forerunner. Finally, we summarize existing ParFlow-based coupled systems and highlight the need for a community-led model intercomparison project (PLCMIP) to benchmark performance, evaluate process coupling under varied configurations, and foster cross-community collaboration.

1 Introduction

In 2005, Maxwell and Miller published Development of a Coupled Land Surface and Groundwater Model in Journal of Hydrometeorology (Maxwell and Miller, 2005). Their work introduced the first coupling of ParFlow (Ashby and Falgout, 1996; Jones and Woodward, 2001) and the Common Land Model (CLM) (Dai et al., 2003), and validated the framework using both synthetic and real-world test cases. The study highlighted the importance of groundwater representation in land surface processes (Fan et al., 2019; Zeng et al., 2018; de Graaf and Stahl, 2022; Seuffert et al., 2002). In particular, it emphasized the role of lateral subsurface flow (Fig. 1), a component that was not explicitly represented in most land surface models (LSMs) at the time. This work represented an early step toward incorporating physically based groundwater dynamics into Earth system modeling frameworks.

Since LSMs serve as the lower boundary in ESMs, this coupling provided a practical pathway to incorporate groundwater dynamics into larger-scale Earth system frameworks. Compared to earlier coupling attempts based on tightly integrated or proprietary platforms (Yeh and Eltahir, 2005; Ivanov et al., 2004; York et al., 2002), this effort leveraged established community models and an open design philosophy, facilitating broader applicability and long-term adaptability. The resulting ParFlow-CLM model and other subsequent models coupled with ParFlow have been applied in a range of hydrological and land–atmosphere studies (Maxwell et al., 2007, 2011; Shrestha et al., 2014), contributing to improved understanding of water and energy exchanges across subsurface, land surface, and atmospheric domains (Rahman et al., 2015; Sulis et al., 2017; Keune et al., 2016; Forrester and Maxwell, 2020). Even today, groundwater–land surface coupling remains underutilized in many large-scale modeling frameworks, where groundwater models are often run offline with limited interaction with land–atmosphere processes (de Graaf et al., 2017; Reinecke et al., 2019; Verkaik et al., 2024), thereby missing dynamic feedbacks with the land–atmosphere system.

The groundwater model, ParFlow, simulates fully 3D variably saturated subsurface flow and overland flow by integrating Richards' equation with the shallow water equation in a unified numerical framework (Kollet and Maxwell, 2006; Osei-Kuffuor et al., 2014; Maxwell, 2013). Meanwhile, the Common Land Model (CLM, now CoLM) captures water and energy processes from the canopy top to the bottom of the root zone. These two models were coupled through the root zone (Fig. 1), where net fluxes from CoLM after the interactions of infiltration and evapotranspiration (ET) are treated as source/sink terms in ParFlow, while ParFlow returns soil moisture and pressure head to CoLM to close the water and energy balance. Such a coupling approach in terms of physics has been widely adopted by the following coupling works (Niu et al., 2014; Fang et al., 2022; Maina et al., 2025).

After two decades of continuous development, LSMs such as CoLM have seen substantial advancements in functionality, code architecture, data structures, I/O systems, pre-/post-processing tools, and high-performance computing capabilities. ParFlow has undergone similar progress on the hydrological modeling front (Kuffour et al., 2020). Although the original coupling between ParFlow and CoLM was once considered sustainable, it is now increasingly inadequate as both models continue to grow in complexity and scope.

At this twenty-year juncture, it is necessary to re-examine the current state of the coupled system and clarify how the next stage of development should proceed. This involves, first, synthesizing the scientific findings enabled by the coupled framework over the past two decades to fully understand its importance for groundwater–land–atmosphere interactions. It also requires an initial assessment of the new coupling – particularly its physical functionality and performance – to clarify the benefits of moving forward. Finally, it is essential to consider how future recoupling can be made sustainable. Together, these steps will support a cross-disciplinary effort that provides a robust platform for the broader community to apply coupled models efficiently, pursue advanced Earth system questions, and strengthen collaborative research. Here we take CoLM as an example to present this transitional effort.

Figure 1 Illustration of the lateral groundwater flow, the critical zone of water table depth, and the coupling strategy between ParFlow and land surface models. Modified from Yang et al. (2023).

In this paper, we begin by reviewing key insights gained from two decades of research involving ParFlow-based coupled modeling systems. Building on this foundation, we highlight how increasing model complexity and functionality are driving a shift toward a next-generation coupling paradigm. We then present a re-coupling of the latest versions of ParFlow (PF) and CoLM, focusing on core functionalities of CoLM to demonstrate feasibility and highlight improvements in overall model performance. This science-oriented integration of basic modules – built upon the original coupling interface – serves as a foundation for broader re-coupling efforts that will incorporate additional functional components under a redesigned, sustainable coupling framework. It also helps us better understand how both models have evolved since their original coupling in 2005, thereby informing the development of a next-generation framework. In recognition of the increasing number of LSMs being coupled with ParFlow, we further propose a ParFlow-Land Surface Coupled Model Intercomparison Project (PLCMIP) to promote collaboration and knowledge exchange across the community.

2 A brief review of ParFlow-Land/Atmosphere coupled modeling

The coupled model provides a more realistic representation of groundwater dynamics than traditional LSMs, while also offering more advanced ecohydrological processes at the land surface than conventional groundwater models. Over the past two decades, its major scientific contributions can be summarized in three key areas:

1. It captures the feedbacks from groundwater to land and atmospheric processes – an area often overlooked or omitted in both atmospheric and groundwater research communities.

2. It highlights the critical range of water table depth (WTD) that mediates these feedbacks.

3. It elucidates the impacts of land cover and climate change on groundwater and other complex ecohydrological interactions.

Because this manuscript is designed to synthesize the past, present, and future of the ParFlow-based coupled systems at this twenty-year juncture, this review primarily synthesizes advances within the ParFlow–LSM and ParFlow–atmosphere modeling family, while also briefly situating other groundwater–land coupling efforts in a broader community context.

2.1 Feedbacks from groundwater to land surface and atmosphere

Adding groundwater representation in ESMs reshapes the spatiotemporal distribution of soil moisture, which in turn controls surface turbulent fluxes and the evolution of the atmospheric boundary layer (Forrester and Maxwell, 2020; Rihani et al., 2015). This is primarily due to the limited simulation depth in LSMs and the absence of lateral groundwater flow. The former limits drainage in ridge areas, resulting in insufficient water release and an overestimation of soil moisture; the latter suppresses groundwater convergence in valley areas, leading to underestimation of soil moisture there.

Generally, lateral groundwater flow enhances soil moisture in topographic lows, suppresses boundary layer development, and increases the evaporative fraction, thereby weakening land–atmosphere coupling and reducing near-surface temperatures (Forrester and Maxwell, 2020; Keune et al., 2016). These responses are further modulated by the subsurface hydraulic conductivities (K), with more pronounced sensitivities to K under simplified groundwater parameterizations (Williams and Maxwell, 2011; Keune et al., 2016; Rihani et al., 2010). Notably, the impact of groundwater and subsurface properties on surface flux partitioning and boundary layer development tends to be most pronounced in the afternoon, when radiative forcing peaks and land–atmosphere interactions intensify (Rahman et al., 2015; Rihani et al., 2015; Forrester and Maxwell, 2020; Maxwell et al., 2007).

Forrester and Maxwell (2020) conducted WRF-based weather simulations over the mountainous regions of Colorado to investigate the impact of different lower boundary conditions, providing a detailed explanation of the processes mentioned above. The study included a baseline scenario and several comparative scenarios, with particular emphasis on one that used PF-WRF to explicitly represent three-dimensional groundwater flow. In the baseline scenario, conventional WRF simulation was employed, with the subsurface depth of 2 m, divided into four layers with thicknesses of 0.1, 0.3, 0.6, and 1 m from top to bottom. The bottom boundary used the native Noah model setting, which allows free drainage and further adjusts fluxes based on terrain. In the PF-WRF scenario, the subsurface depth was increased to 102 m by adding a fifth layer of 100 m in thickness, with the bottom boundary set as impermeable. The Noah model and ParFlow were coupled through the top four layers, resulting in a coupling depth of 2 m.

Simulation results showed that, in the PF-WRF scenario, enhanced drainage over ridge areas reduced soil moisture, while lateral groundwater convergence increased soil moisture in valleys. Correspondingly, the boundary layer height also exhibited increases in ridge areas and decreases in valley areas. These changes in soil moisture and boundary layer height showed significant seasonal variations. Furthermore, the results revealed that microtopography induced highly heterogeneous local variations in soil moisture. This, in turn, weakened the clear elevation-dependent trend observed in the baseline scenario.

Additionally, in the baseline scenario, the coupling strength between evaporative fraction (EF, the ratio of latent heat to the sum of latent and sensible heat) and boundary layer height was weakened or even reversed in the PF-WRF scenario. That is, the significant negative correlation between EF and boundary layer height decreased or turned positive; this may be due to the temporal variations in EF caused by lateral flow. Moreover, the PF-WRF scenario with lateral flow showed stronger morning mountain breezes (upslope) and valley breezes (downslope), which may have enhanced mountain-valley circulation. Lateral groundwater flow also modulated low-level convection in river valleys, particularly increasing convective available potential energy (CAPE) in the afternoon, thereby perturbing regional precipitation.

Keune et al. (2016) conducted simulations over the European CORDEX region using the TerrSysMP modeling system (Shrestha et al., 2014), setting up two scenarios: one with fully three-dimensional groundwater flow (3D) and the other with one-dimensional free drainage (FD). Similarly, their results revealed that different representations of groundwater led to variations in CAPE, indicating influences on the evolution of atmospheric boundary layer and free troposphere. The 3D scenario weakened land–atmosphere coupling, thereby suppressing the occurrence of extreme weather events, which is consistent with the findings of Forrester and Maxwell (2020). More specifically, the simulated 2 m air temperature was generally lower in the 3D scenario than in the FD scenario, providing useful insights for simulating European heatwaves during the study period.

The study also showed that model differences were primarily located in areas with shallow water tables (depth < 5 m), which aligns with findings of Forrester and Maxwell (2020) that humidity, potential temperature, and vertical wind exhibit more pronounced differences in mountainous valley regions. In addition, the study revealed that variations in deep soil (depth > 3 m) hydraulic conductivities led to discrepancies in simulation results. The FD scenario was more sensitive to the choice of conductivity values, suggesting that simplified physical representations may further amplify the impact of parameter uncertainty.

Williams and Maxwell (2011), using coupled PF-WRF simulations, further explored the feedbacks of geological conditions on land–atmosphere processes such as latent heat flux and wind speed. Based on idealized scenarios, they conducted ensemble simulations by perturbing the hydraulic conductivity field. The results showed that conditioning the hydraulic conductivity significantly reduced uncertainties in simulating land–atmosphere interactions compared to unconditioned cases. The ensemble mean was closer to the control scenario; for instance, the mean and distribution of simulated wind speed showed reduced uncertainty. These findings provide important implications for various wind energy applications.

Community-wide studies have also highlighted the importance of representing water-table dynamics within land-surface processes. Koirala et al. (2014) incorporated groundwater fluctuations into the MATSIRO land surface scheme and quantified the sensitivity of ET to capillary rise, showing that global mean ET increases by approximately 9 % when water-table dynamics are included. Tian et al. (2012) further examined how hydraulic conductivities regulate ET by influencing both vertical and lateral groundwater fluxes, using a coupled AquiferFlow–SiB2 modeling framework. Using ParFlow-based coupled models, Tai et al. (2018) and Fang et al. (2022) demonstrated that explicitly resolving water-table dynamics helps explain mechanisms of plant mortality, while Abbaszadeh et al. (2025) reported improved simulations of land-surface fluxes when groundwater processes are represented. Miguez-Macho and Fan (2025) incorporated lateral surface-water and groundwater subsidies simulated by the ASAP model into a humidity index, providing a more accurate depiction of the timing and magnitude of water availability in hydrologically convergent lowlands. This enhancement better explains the monthly variations of leaf area index. A more recent study (Vogelbacher et al., 2025) extends beyond physically coupled modeling frameworks by showing that integrating water-table depth into a machine-learning system yields a more robust assessment of heatwave risks.

2.2 The critical zone of WTD in groundwater–land interactions

As discussed above, numerous studies have revealed feedbacks of groundwater on land–atmosphere processes. A key scientific question thus arises: what is the quantitative relationship between land surface states/fluxes and the WTD? Maxwell and Condon (2016), in their study over the continental US, confirmed the critical role of lateral groundwater flow in modulating the partitioning between evaporation (E) and transpiration (T). This influence is most pronounced when the WTD lies between 0.5 and 5 m. Shallower WTD leads to elevated bare-soil evaporation and transpiration, while deeper WTD suppresses both fluxes. Notably, in regions where bare-soil evaporation is limited and transpiration is sustained, the $T/E$ ratio peaks.

Similarly, many studies using PF-CLM have identified a critical WTD range within which land surface variables – such as latent heat flux, sensible heat flux, and surface temperature – are highly sensitive to WTD but exhibit diminished sensitivity beyond this range (Fig. 1). For instance, Ferguson's work over the Little Washita watershed suggests a critical WTD range of approximately 1–10 m (Ferguson and Maxwell, 2012, 2011, 2010), while Yang et al. (2020, 2023) reported comparable results over the North China Plain. Rihani et al. (2015) also illustrate the coupling between WTD and planetary boundary layer depth in this transition zone from ridges to valleys along hillslopes. Fang et al. (2022) demonstrated a clear transition from hydraulic-failure-dominated mortality at shallow WTDs (< 5 m) to carbon-starvation-dominated mortality under deep water-table conditions (> 7.5–15 m). Generally, when WTD is shallower than this range, soil is nearly saturated and energy availability becomes the limiting factor, weakening the sensitivity of surface states/fluxes to WTD. Conversely, when WTD exceeds this range, gravity-driven drainage dominates, limiting moisture availability and again reducing sensitivity. The upper bound of this range is typically < 1 m, while the lower bound often aligns with the model's coupling depth (Kollet and Maxwell, 2008). However, in some cases, such as Maxwell and Condon (2016), the lower bound extends beyond the nominal 2 m coupling depth, likely due to capillary rise from the water table.

This critical WTD range varies across regions, influenced by differences in subsurface characteristics and rooting depth, though current understanding remains limited. Fan et al. (2017), through analysis of over 2200 global root depth observations and model-based inversion, showed that rooting depth is regulated by the capillary rise zone. Even within the same species and climate, rooting depth may vary with WTD conditions (Cannon, 1913). In some environments, vegetation develops both shallow fibrous roots and deep taproots to access water under varying conditions – shallow roots for near-surface moisture during wet periods, and deep roots for capillary water during droughts. On well-drained uplands, rooting depth is controlled by infiltration and may not reach significant depths (Sperry and Hacke, 2002; Schenk and Jackson, 2005; Xu and Li, 2008). In contrast, in shallow groundwater zones, oxygen stress may inhibit root growth and decouple vegetation from groundwater (Wagg, 1967; Follett et al., 1974; Martin, 1968; Armstrong et al., 1976). This adaptive rooting strategy suggests that in natural systems, the depth and intensity of groundwater–land surface coupling may exceed what models typically simulate.

Maxwell et al. (2007) further investigated how groundwater feedbacks on land surface processes change under different climate change scenarios, including hot, hot and dry, and hot and wet. It is not difficult to infer that changes in latent heat flux and recharge (precipitation minus ET) within the critical WTD range exhibit strongest spatial variability. Generally, these findings suggest that groundwater feedbacks on land surface processes are closely linked to topographic and climatic conditions. For instance, in the aforementioned mountainous regions of Colorado, spatial variability in WTD leads to diverse groundwater–land surface interactions. The transitional zones between ridges and valleys are often the key areas for such interactions. In humid regions, the water table often follows surface topography (Gleeson et al., 2011), facilitating strong groundwater–land surface coupling. However, in arid regions, WTD may exceed the lower bound of the critical range, reducing the significance of groundwater feedbacks. In natural systems, this interaction is often governed by the complex interplay between climate, topography, geology, and vegetation.

2.3 The impacts of land cover and climate changes on groundwater

Climate change has exacerbated mountain pine beetle infestations, leading to widespread tree mortality in the Rocky Mountains (Bearup et al., 2014). Mikkelson et al. (2013) studied the impacts of beetle-induced forest dieback on water and energy balances at the hillslope scale using PF-CLM simulations. An idealized hillslope model (500 m × 1000 m × 12.5 m) was used, with scenarios representing different stages of infestation – green, red, grey, and dieback – by modifying the leaf area index and stomatal conductance. Simulation results showed similar levels of ET across all scenarios in winter, but significantly higher ET in summer under the green scenario, primarily due to transpiration. In contrast, the other scenarios exhibited lower ET limited by soil moisture availability, with evaporation being the dominant process. The dieback scenario produced the highest peak in snow water equivalent (SWE), and reduced canopy cover allowed more solar radiation to penetrate, accelerating snowmelt. This earlier and more rapid melt resulted in earlier and higher streamflow peaks, as well as increased subsurface storage. A related particle tracking study (Bearup et al., 2016) further demonstrated greater groundwater contributions to streamflow during late summer.

Condon et al. (2020) conducted a continental-scale simulation across the United States using the PF-CLM CONUS 1.0 model (Maxwell et al., 2015) to examine groundwater responses to 1, 2, and 4 °C warming scenarios. Warming was found to enhance ET, with shallow groundwater providing supplementary moisture to meet the increased demand, thereby partially mitigating land surface water stress. However, prolonged warming ultimately led to continuous groundwater depletion and a decoupling of groundwater from land surface processes. The magnitude of ET increases, and groundwater storage loss varied with WTD, with the strongest responses occurring within the previously identified critical WTD range. Overall, the humid eastern U.S. exhibited greater sensitivity to warming than the arid western regions. These findings highlight the risk of underestimating groundwater–land surface feedbacks when using simplified groundwater parameterizations in ESMs.

2.4 Enhanced LSM functionality motivates recoupling

These selected representative studies have demonstrated the critical role of groundwater in Earth system processes. Over the past two decades, CoLM – like many other LSMs – has undergone substantial development, including functional extensions, improved parameterization schemes, and the introduction of multiple alternative process representations. However, our understanding of how groundwater interacts with these additional processes – including how various parameterizations respond to and influence groundwater dynamics – remains limited. These limitations underscore the pressing need to upgrade the coupling between ParFlow and LSMs. The key scientific advances of CoLM are summarized as follows (Yuan and Dai, 2025):

1. Radiation transfer. A three-dimensional vegetation shortwave (Yuan et al., 2014) and longwave radiation transfer scheme has been incorporated, the SNICAR snow radiation transfer scheme has been added to simulate snow albedo and radiation absorption within the snowpack, and an improved two-stream approximation scheme for vegetation radiation transfer has been provided (Yuan et al., 2017).

2. Turbulent fluxes. The model enhances the continuity of dynamic parameters and processes across transitions from dense to sparse vegetation. Resistance coefficients below the canopy are calculated using a profile integration method. A new turbulent exchange scheme supports multiple coexisting plant functional types (PFTs) within a three-dimensional canopy. Several soil resistance parameterizations are provided to improve surface evapotranspiration estimates. Additionally, a surface turbulent flux scheme has been introduced to account for large-eddy effects.

3. Canopy interception and plant hydraulics. The model includes multiple canopy interception schemes and a plant hydraulics module governed by Darcy's law. Different parameterizations emphasize distinct physical processes and support investigation of the evolution, drivers, and trends of interception under varying conditions. The hydraulics module replaces empirical formulations that relate plant stress to soil water potential and improves the simulation of land-atmosphere water exchange under changing environments.

4. Leaf temperature. A simplified one-dimensional two-big-leaf scheme has been implemented to improve the numerical stability of leaf temperature simulations. In addition, a new parameterization has been developed for leaf temperature in multi-PFT scenarios with a three-dimensional canopy structure.

5. Other functional extensions. Additional modules have been developed for biogeochemistry, urban systems, crop modeling, land use and land cover change, wildfire, ozone-related ecophysiological stress, and integrated hydrological processes.

3 Foundational step toward sustainable coupling

To explore the feasibility of re-coupling the latest versions of ParFlow and CoLM, we conducted a preliminary integration of the two models. This effort focuses exclusively on the basic water and energy modules of CoLM and is built upon the existing ParFlow-CLM coupling interface. That means that canopy interception, snow processes, and the surface energy balance (radiation, sensible and latent heat fluxes, and ground heat flux) were enabled in the current tests, whereas other functional extensions were kept inactive. Our goal is to understand how both models have advanced over the past two decades, in terms of functionality, code architecture (e.g., parallelism), data structures, I/O interfaces, and pre-/post-processing tools. This process also helps identify key variables and processes – along with their implementation in code – that are critical for a more comprehensive coupling effort. Although this re-coupling effort uses CoLM as an example, the experience and insights gained are also applicable to coupling ParFlow with other LSMs.

This initial re-coupling serves to evaluate model performance with respect to physical processes, particularly highlighting potential improvements gained through two decades of development. It also establishes a set of benchmarks for testing and debugging as more CoLM modules are progressively incorporated. Without this incremental approach, the complexity of multiple physical processes would make testing and debugging considerably challenging. In addition, we fully leverage the lessons learned through trial and error over the past twenty years to ensure a more stable execution of the coupled model (Ferguson et al., 2016). While this phase emphasizes gaining a deeper understanding of the physical processes, future work on sustainable coupling will likely shift toward technical aspects – such as refining the coupling interface, improving modularity, and ensuring long-term maintainability. These efforts will collectively inform our understanding of the challenges and opportunities involved in establishing a sustainable coupling framework.

3.1 Model setup, experimental design, and evaluation data

The modeling domain, selected from the CONCN domain (Yang et al., 2025), is located in the North Pearl River Basin (Fig. 2). This area was chosen as it serves as a demonstration area for the CONCN model and possesses more complete infrastructure e.g., the processed ERA5-Land reanalysis data (Muñoz-Sabater et al., 2021). To link the study domain with the coupled modeling framework, we next describe the model's vertical and horizontal discretization. Since the CONCN model and CoLM use four and ten soil layers, respectively, the CONCN model structure was adjusted to align with CoLM's vertical discretization. The ParFlow model employed in this study comprises 11 layers: the top 10 layers match CoLM in thickness, while the additional 11th layer represents the deep aquifer. ParFlow and CoLM are coupled through the top 10 layers. The coupled model maintains a horizontal resolution of ∼ 1 km, consistent with the CONCN model and includes 252 and 146 grid cells in the x and y directions, respectively. This corresponds to a spatial extent of approximately 242.35 km (x direction) × 140.41 km (y direction) × 103.43 m (z direction). Soil properties for CoLM inputs were derived from the Global Soil Dataset for Earth System Modeling (GSDE) (Shangguan et al., 2014). Soil parameters for ParFlow were reconstructed based on the sand and clay weight percentages from the same GSDE dataset, following the USDA soil classification system. Properties for the 11th layer were obtained from GLHYMPS 1.0 (Gleeson et al., 2011; Gleeson et al., 2014), as used in the CONCN model. The e-folding of aquifer hydraulic conductivity with depth was implemented using a characteristic depth of 50 m (Fan et al., 2007). Other surface input parameters – including Manning's roughness coefficients, topographic slopes, and land cover types – were adopted directly from the CONCN model configuration (Yang et al., 2025).

Figure 2 Location of the modeling domain.

We first spun up the standalone ParFlow model using potential recharge data clipped from the CONCN domain. Then we drove the coupled model using the 2018 meteorological forcing from ERA5-Land reanalysis. Year 2018 was selected to keep consistency with the evaluation year of the CONCN model. However, this is an area with limited snow. To demonstrate snow performance, we created a synthetic case by applying the water year 2003 forcing from a station (Defnet et al., 2024) located in Colorado to a single column model. To evaluate changes in model performance, we also constructed 11-layer models (with 10 coupled layers) using the old ParFlow-CLM for both real-world and synthetic cases. For the real-world case, we compared the simulated sensible heat, latent heat, skin temperature, transpiration, SWE, and the water flux exchange between the new and old models (Figs. 3 and 4). Here, water flux exchange refers to the source/sink terms in Richards' equation: positive values represent infiltration, while negative values are caused by ET. We also evaluated the simulation performance of the first four variables using ERA5-Land reanalysis. For the synthetic cases, we used data from the Snow Telemetry (SNOTEL) network maintained by the Natural Resources Conservation Service (NRCS) – specifically, the measured SWE at the same location as the meteorological forcing – to evaluate the models' overall ability to simulate the timing and magnitude of snowpack.

Figure 3 Comparison of latent heat flux, sensible heat flux, surface temperature, and transpiration between the old CLM/ParFlow and the new CoLM/ParFlow models. The corresponding values from ERA5-Land are plotted in the left column for reference. The right column shows the differences between the model simulations and ERA5-Land for each variable. Each subplot represents spatial averages over the entire modeling domain. For clarity and to prevent overlapping, the plotting order is intentionally varied across subplots.

Figure 4 Panel (a) shows the simulated snow water equivalent from the CoLM/ParFlow and CLM/ParFlow single column models, compared against SNOTEL observations. Panel (b) presents the spatially averaged net water fluxes from CoLM and CLM to ParFlow, representing the source and sink terms in the ParFlow domain; both fluxes represent spatial averages over the entire modeling domain.

3.2 Performance gains from updated CoLM support recoupling

Simulations of all variables by CoLM/PF exhibit improved performance relative to CLM/PF when evaluated against ERA5-Land reanalysis data (Fig. 3). CoLM/PF produces a more realistic partitioning of turbulent fluxes, characterized by increased latent heat and reduced sensible heat (Fig. 3a and c). Notably, transpiration simulated by CoLM/PF is substantially higher and aligns more closely with ERA5-Land data (Fig. 3g). Additionally, CoLM/PF more accurately reproduces the fluctuations in land surface temperature compared to CLM/PF (Fig. 3e). In the single column simulations (Fig. 4a), CoLM/PF also generates a higher peak SWE than CLM/PF, showing better agreement with SNOTEL observations, although both models display deviations in the timing of SWE accumulation. This discrepancy may stem from the idealized subsurface configurations used in both models. Improvements in both transpiration and SWE are further supported by previous research; for instance, O'Neill et al. (2021) reported consistently lower ET and SWE from CLM/PF in the assessment of CONUS 1.0 model. The overall advancement in model performance can likely be attributed to the more sophisticated process representations embedded in CoLM.

Figure 5 Panels (a) and (b) show the simulated water table depth and overland flow by CoLM/ParFlow and CLM/ParFlow, respectively. Each subplot represents spatial averages over the entire modeling domain.

Figure 4b illustrates the net fluxes transferred from the land surface to the subsurface, which directly influence hydrologic dynamics such as WTD and overland flow. These exchange fluxes show greater variability in the CoLM/PF simulation, suggesting more dynamic surface–subsurface interactions. Consistent with this, Fig. 5 reveals more pronounced temporal variability in both WTD and overland flow. A generally deeper water table is observed in CoLM/PF (Fig. 5a), which is likely a result of the higher plant water uptake, i.e., increased transpiration, depicted in Figs. 3a, g and 4b. Consequently, the reduction in baseflow from groundwater may explain the observed decrease in low levels of overland flow (Fig. 5b).

To improve the representation of turbulent exchanges between the vegetation canopy and the atmosphere, the model employs a profile-integrated approach to resolve key dynamical parameters (e.g., turbulent diffusivity K(z)) with explicit vertical resolution (Dai et al., 2019b). In particular, resistance-related variables – such as displacement height (d) and roughness length (z₀), which characterize canopy-atmosphere momentum exchange, as well as aerodynamic resistances for leaves (r_b) and ground surface (r_d), which govern within-canopy and near-surface heat and vapor transfer – are refined to account for structural heterogeneity (Dai et al., 2019b). Meanwhile, profile-integrated functions are dynamically computed based on vegetation structure and atmospheric stability, and directly determine resistance terms (e.g., r_ah,r_aw) (Yuan and Dai, 2025). This also includes revised roughness length formulations that explicitly account for atmospheric stability (Yuan and Dai, 2025), extending beyond the original neutral-based assumptions in schemes such as Raupach (1994, 1992). This combined approach yields a more physically consistent and vertically continuous treatment of turbulent fluxes under non-neutral stratification, enhancing realism in complex canopy conditions.

In addition, other schemes – such as those related to soil thermal parameters (including heat capacity and heat conductivity), as well as soil color and associated reflectance – also differ between the CoLM/PF and CLM/PF models (Yuan and Dai, 2025). All of these differences motivate a coupled model intercomparison project to evaluate how different schemes, either within a single model or across different models, affect the performance of the coupled model and land–hydrology process interactions. Since 2005, ParFlow has also undergone several major developments that significantly affect model performance (Kuffour et al., 2020), such as the integration of overland flow (Kollet and Maxwell, 2006) and the implementation of the terrain-following grid (Maxwell, 2013). Both test cases in this study already used these modern capabilities. Features such as GPU acceleration (Hokkanen et al., 2021) and reservoir capabilities (West et al., 2025) in ParFlow exist, but they are not relevant to the test cases here.

In this initial coupling process, we found that the main changes in CoLM were related to data structure, module organization, module names, and variable naming conventions. For example, structures were broken down into multiple arrays, modules were split and reorganized based on functionality, some modules were removed from the main program and used as preprocessing components, and module adjustments were often accompanied by renaming. A large number of variable names were also changed. Moreover, the inclusion of multiple parameterization schemes has increased code complexity to some extent, resulting in significantly larger module sizes. Nevertheless, most physical processes have retained their original core parameterizations. This means that the primary task in this initial coupling stage is to identify the key physical processes and the critical variables within the new system structure. Several new parameterization schemes have also been implemented – for example, those associated with the turbulent exchange discussed above – though the application of other schemes will require further testing in future work.

4 A sustainable recoupling framework for future development

Here, we propose a sustainable framework for future CoLM/ParFlow coupling based on our preliminary work (Fig. 6). This framework consists of four key components: a coupler-based architecture, a robust initial foundation, protocols for scalable upgrades, and a community interaction platform. The four components play different roles. The first two – coupler-based architecture and early-stage grid/process mapping – form the structural foundation that determines long-term sustainability. They address issues that cannot be solved through implementation refinements alone, such as maintaining independent model evolution and enabling robust cross-grid exchanges. The latter two components – developer protocols and community interaction – further enhance maintainability and extensibility once the structural basis is in place. These collaborative conventions support scalable development but cannot replace the structural requirements themselves, nor can they reduce the inherent complexity of multi-model coupling. Effective participation still requires sufficient domain knowledge and technical expertise.

4.1 Coupler-based architecture for long-term sustainability

While the current ParFlow-Land interface built in ParFlow supports efficient coupling with land surface and atmospheric models, demonstrates good parallel performance, and avoids the overhead associated with inter-model communication, it lacks compatibility with standardized coupling frameworks and protocols. This limits the integration of coupled models into broader Earth system modeling frameworks. In contrast, coupler-based architectures – such as ESMF/NUOPC, CESM/cpl7, and OASIS3-MCT – are now standard in modern Earth system modeling. They preserve the native data structures, domain decomposition, and parallel logic of each model, which is particularly important given the substantial structural differences between ParFlow and CoLM. For instance, this approach allows retaining ParFlow's GPU-based parallelism (Hokkanen et al., 2021) and CoLM's MPI-based structure, along with their respective domain decomposition strategies. It also enables continued use of each model's preferred data format and processing tools – for example, ParFlow's .pfb format and pftools, as well as CoLM's NetCDF-based workflow. Adopting a standardized coupler thus facilitates modular development, cross-system compatibility, and long-term maintainability. Looking forward, such coupler-based designs could also be extended to support surrogate model integration (Bennett et al., 2024; Tran et al., 2021), enabling hybrid workflows that combine physical models and AI-based components.

4.2 Strong foundations through early-stage mapping

Laying a solid foundation during the initial coupling phase is critical. First, a mapping between ParFlow's grid and CoLM's subgrids must be established to support future model extensions. The current coupling only supports integration based on the LCT (Land Cover Type) subgrid, whereas important functional extensions such as biogeochemistry and 3D vegetation canopy processes rely on plant function types (PFTs) and the associated PFT and PC (Plant Community) subgrids. The mapping between grids of different models is a key concern in the coupling community and directly affects coupling performance (Valcke, 2022). Currently, most approaches rely on physical interpolation methods; incorporating an AI-based scale transformation layer with mass conservation constraints could be a promising enhancement. Second, key variables and processes from newly introduced modules or previously untested parameterizations must be identified. For example, the plant hydraulics module requires soil hydraulic conductivity fields, which are not included in the current coupling interface. More importantly, a structured logging system should be implemented to track all exchanged variables, their associated modules, and the corresponding grid structures, thereby ensuring transparency and traceability throughout development. This logging mechanism differs from conventional version-control systems in that it records the evolution of coupling-relevant variables and grid mappings rather than general code revisions, providing scientific rather than software-level traceability.

4.3 Protocols for efficient and maintainable coupler upgrades

Two key aspects are emphasized. First, the architecture of interfaces and coupling layers should be designed for long-term clarity and ease of maintenance. Taking ESMF/NUOPC as an example, the model-side interfaces and coupler-side connector and mediator layers are implemented with a focus on modular organization, encapsulated data exchange, and well-structured control flow, ensuring that the system can be reliably extended as new model features or physical processes are introduced. For example, TerrSysMP achieves good modularity by separating tasks such as sending/receiving fields, grid mapping, and variable registration into dedicated files (Shrestha et al., 2014). TerrSysMP2 further improves on this design by organizing these components within a single Fortran module, which makes the interface structure clearer and simplifies the CMake build process (Poll et al., 2024). Second, developers introducing new modules, parallelization strategies, or grid structures, within one model must explicitly assess their potential impact on the other model, clarify any newly introduced variables or data structures to be exchanged via the coupler, and submit pull requests with corresponding explanations. Senior maintainers should review these changes and provide targeted feedback on necessary updates to the interface and coupling logic. All modifications affecting model interaction must be tracked in the logging system described above, and no update should be considered complete until it is formally registered in the log.

4.4 Community platform for collaboration and maintenance

A dedicated community platform – such as a GitHub repository, mailing list, or model portal – should be established to support developer–user interaction, technical discussion, and feedback collection. This platform will also serve to announce new model releases, coupling layer updates, or changes in the logging system. Transparent communication and community-driven collaboration are essential for the long-term sustainability and extensibility of the coupled model system. Looking forward, we envision extending this platform into a broader, community-driven environment for managing and operating ParFlow–LSM coupled models, drawing inspiration from efforts such as eWaterCycle (Hut et al., 2022). Beyond supporting collaboration and information sharing, such a platform could streamline model configuration, coupled-model execution workflows, data processing/formatting, and reproducibility – thereby accelerating adoption, improving transparency, and fostering interoperability across hydrological and Earth system science communities. This would also address a current gap, as ParFlow-based coupled systems remain largely fragmented and lack standardized tools for interfacing with upstream/downstream models, verification and benchmarking, and public release to users.

Figure 6 The coupling framework for sustainable development.

5 ParFlow-Land coupled model intercomparison project (PLCMIP)

Over the past decades, numerous land and/or atmosphere models coupled with ParFlow have been developed (Table 1). These models vary in their functional capabilities and adopt different coupling strategies, which may significantly affect computational efficiency. The two models coupled with CoLM aim to understand the fundamental interactions of water and energy between subsurface and land surface processes (Dai et al., 2003; Maxwell and Miller, 2005). In contrast, the two models coupled with ARPS and WRF (Maxwell et al., 2007; Maxwell et al., 2011; Skamarock and Klemp, 2008; Xue et al., 2000; Xue et al., 2001), along with the two generations of TerrSysMP (Shrestha et al., 2014; Oleson et al., 2008; Lawrence et al., 2019; Poll et al., 2024), provide capabilities to explore two-way feedbacks across each interface within the subsurface–land surface–atmosphere system. Furthermore, the coupling of ParFlow with TREES (Tai et al., 2018; Mackay et al., 2015), ELM-FATES (Fang et al., 2022; Caldwell et al., 2019; Fisher et al., 2015; Leung et al., 2020), and LPJ-GUESS (Jia et al., 2025) introduces advanced vegetation dynamics into land surface process representations. Finally, integration with NASA-LIS enables data assimilation within the coupled modeling framework (Maina et al., 2025; Kumar et al., 2008; Niu et al., 2011; Abbaszadeh et al., 2025), and TerrSysMP also incorporates the PDAF (Parallel Data Assimilation Framework) to support data assimilation capabilities (Kurtz et al., 2016). Overall, most of the ParFlow-based coupled systems are implemented through modular integration–embedding one model within another–whereas TSMP and ParFlow-LIS represent coupler-based architectures that mediate data exchange among components.

Table 1 ParFlow and Land/Atmosphere coupled models.

Model intercomparison provides a valuable means to assess model development and foster connections or collaborations across research communities. Several well-known intercomparison projects exist, such as the Coupled Model Intercomparison Project (CMIP) for ESM intercomparison (Eyring et al., 2016), and the Land Surface, Snow and Soil Moisture Model Intercomparison Project (LS3MIP) (van den Hurk et al., 2016), which is designed to assess the performance of land modules in current ESMs. In addition, individual intercomparison activities have also been widely conducted within the land surface modeling community (Scanlon et al., 2018; Liu et al., 2023). ParFlow has also participated in various model intercomparison projects involving hydrologic models and an increasing number of individual modeling studies, such as those by Maxwell et al. (2014), Sulis et al. (2017), Kollet et al. (2017), Sulis et al. (2010). Given the differences among the ParFlow-based coupled models mentioned above, a dedicated model intercomparison project (MIP) is needed to systematically evaluate coupled models and support the development of a community platform for benchmarking and collaboration, with the following objectives:

• To quantify the strength and spatiotemporal variability of groundwater–land–atmosphere interactions resulting from different parameterization schemes used in various land surface and atmospheric models.

• To evaluate the parameter sensitivity of each scheme, ensuring that differences attributed to model structure are not confounded with parameter choices.

• To compare computational efficiency across different coupling strategies.

• To identify the unique functionalities and strengths of each coupled model, providing users with guidance in selecting the most appropriate model for their specific research needs.

To ensure meaningful and comparable evaluations across models, the PLCMIP will encourage the use of standardized benchmark cases – either synthetic experiments or a common real-world watershed – as well as unified datasets for parameters and meteorological forcing. In addition, other groundwater-land coupled models, as well as land surface models with improved groundwater parameterizations, are likewise encouraged to participate in this intercomparison effort. Representative examples include Shen et al. (2016), Zeng et al. (2018), Tian et al. (2012), Niu et al. (2014), Sutanudjaja et al. (2014), Liao et al. (2025), Miguez-Macho and Fan (2025), Akhter et al. (2025), Dai et al. (2019a), although participation in PLCMIP is not limited to these.

6 Summary

Twenty years after the original ParFlow-CLM coupling (Maxwell and Miller, 2005), this study reaffirms the long-term scientific and technical significance of that foundational effort. Over two decades, the coupled system has made major contributions in establishing the critical role of groundwater in modulating subsurface–land–atmosphere feedbacks and identifying the existence of a critical water table depth range that governs these bidirectional interactions. Technically, this coupling demonstrated a viable approach for integrating a groundwater model with a land surface scheme – the lower boundary of Earth system models – thereby providing a template for incorporating groundwater processes into ESMs. To revisit and update this legacy, we carried out a preliminary re-coupling of the latest versions of ParFlow and CoLM, focusing on core water and energy processes. This re-coupling already reveals improved model performance and provides a functional platform for incremental expansion and benchmarking.

Looking ahead, several key steps are essential for advancing a sustainable and extensible ParFlow–LSM coupling framework. Achieving this vision will require a more comprehensive and community-oriented design. This includes adopting a lightweight coupler architecture that preserves each model's native data structures, parallel strategies, and processing tools, while supporting modular integration of new physical or surrogate components. To ensure long-term maintainability and usability, we also envision a community platform that unifies model configuration, user workflows, and benchmarking functions. Such a platform would enhance transparency, reproducibility, and ease of adoption across the hydrologic and Earth system modeling communities. In parallel, we propose launching a model intercomparison project (PLCMIP) to systematically evaluate performance, compare coupling strategies, and guide future development.