Process‐Based Intercomparison of Water Isotope‐Enabled Models and Reanalysis Nudging Effects

The products from the Stable Water Isotope Intercomparison Group, Phase 2, are currently used for numerous studies, allowing water isotope model‐data comparisons with various isotope‐enabled atmospheric general circulation model (AGCMs) outputs. However, the simulations under this framework were performed using different parameterizations and forcings. Therefore, a uniform experimental design with state‐of‐the‐art AGCMs is required to interpret isotope observations rigorously. Here, we evaluate the outputs from three isotope‐enabled numerical models nudged by three different reanalysis products and investigate the ability of the isotope‐enabled AGCMs to reproduce the spatial and temporal patterns of water isotopic composition observed at the surface and in the atmospheric airborne water. Through correlation analyses at various spatial and temporal scales, we found that the model's performance depends on the model or reanalysis we use, the observations we compare, and the vertical levels we select. Moreover, we employed the stable isotope mass balance method to conduct decomposition analyses on the ratio of isotopic changes in the atmosphere. Our goal was to elucidate the spread in simulated atmospheric column δ18O, which is influenced by factors such as evaporation, precipitation, and horizontal moisture flux. Satisfying the law of conservation of water isotopes, this budget method is expected to explain various fractionation phenomena in atmospheric meteorological and climatic events. It also aims to highlight the spreads in modeled isotope results among different experiments using multiple models and reanalyses, which are primarily dominated by uncertainties in moisture flux and precipitation, respectively.

The traditional notation for isotopes is denoted by a delta (δ), and it is used to express a small amount of heavy water isotopes in a sample, specifically H 2 18 O and HD 16 O (where deuterium D is 2 H), which constitute 0.20% and 0.02% of total water, respectively.This notation is expressed in ‰ (per mil) and is calculated using the × 1, 000 .In this equation, R represents the ratio of heavy to light isotopes; R sample is the isotope ratio in the water sample, and R standard is the isotope ratio in the reference standard, known as Vienna Standard Mean Ocean Water (D/ 1 H = 155.76± 0.1 ppm, 18 O/ 16 O = 2005.20 ± 0.43 ppm).
Early comparisons of modeled water isotopes examined the performance of present-day and Last Glacial Maximum simulations using GISS and ECHAM models (Hoffmann et al., 2000;Jouzel et al., 2000).Subsequently, a larger comparison effort was undertaken by the Stable Water Isotope Intercomparison Group (SWING) to assess the reliability of simulation results amongst multiple models (Werner et al., 2004).The latest and most commonly used intercomparison simulation data set is the SWING phase 2 (Risi et al., 2012).This data set has supported numerous studies that interpret measured isotope observations or compare those findings with other models (Conroy et al., 2013;S. Dee et al., 2015;Hu et al., 2018;Shi et al., 2022).However, simulations under this framework have not been nudged using consensus reanalysis data sets, and some results were free-running simulations that followed the experiment design for the comparison of Atmospheric Model Intercomparison Project simulations (Gates, 1992).Therefore, it is unclear whether the differences in water isotope values are caused by the choice of model or the choice of nudged reanalysis.This lack of clarity has motivated our research, so we propose a new uniform experimental design with multiple models for a more rigorous intercomparison and precise uncertainty quantification.
One challenge in comparing water isotope models lies in determining the causes of variations in simulation spreads.Different reanalyses for nudging and variations in model structures-such as tracer transport and moist convection schemes-can significantly affect the spread in modeled isotope values.To achieve simulations and analyses for comparison, we first reconstructed the isotopic fractionation processes from 1979 to 2020 using the nudging technique with representative reanalysis fields in the three models.These were run under consistent present-day boundary conditions, such as orbital parameters and greenhouse gases.The model's performance was verified by comparing it with observed water isotope data, and it was evaluated based on regions, seasons, BONG ET AL. 10.1029/2023JD038719 3 of 28 and variables such as temperature or precipitation that are highly correlated with water isotope values.However, the model's performance in representing water isotopes can be highly dependent on the observation locations used for evaluation.Additionally, it is difficult to quantify the uncertainties associated with the spreads in model simulations for each fractionation process, from evaporation to precipitation.
In light of this, we need to employ quantitative methods for diagnosing differences caused by model structures, water isotope-related physics, and atmospheric circulations.The necessary methodology involves utilizing a mass budget analysis for the heavy water isotopes related to their change rates in the atmosphere.This rate is expressed as "‰/day" and accounts for various decomposed processes in the water cycle, including evaporation, water transport, and precipitation (while considering fractionation effects).Since the numerical models are designed to solve the governing equations and thereby conserve mass, momentum, and energy, we expect the masses of the water isotopes to be conserved.This means that any changes in the distribution of water isotopes can be accounted for within a closed system.Therefore, by applying this methodology to a range of nudging simulations using different models and reanalyses, we not only endeavor to pinpoint the factors leading to a variety of isotope values but also capitalize on the opportunity to quantify uncertainties in fractionation throughout hydrological processes, at the surface and in the atmosphere.
For the first time in our study, we estimated the decomposed isotopic change rates.The global estimation over the past 42 years explains the fractionation imbalance between the evaporation-derived enriching rate (about 2‰/ day) and the precipitation-derived depleting rate (about −3‰/day), highlighting the important role of alleviating this imbalance by the enrichment of horizontal moisture flux-derived fractionation (about 1‰/day) in the atmosphere.Through the incorporation of variations in these estimates, a process-based decomposition approach can be employed to quantify uncertainties in water isotope models.Furthermore, the analysis can indicate how far away the water vapor was transported after being evaporated, as determined by the surrounding temperature and humidity conditions in meteorological phenomena on both the synoptic and mesoscale levels.From a long-range perspective, trend analysis can be used to identify which hydrological processes have influenced the isotope ratios in the atmosphere, and by how much some of decomposed processes will likely change in the future.
This study is the first comparison of water isotope models since the SWING2 project, conducted over 10 years ago.In order to advance to the next phase of the modeled isotope comparison, we used updated observation data sets and performed historical, present-day simulations with newer versions of three of the seven models from SWING2.The remainder of this paper is structured as follows: Section 2 describes the experimental design, models adopted, data sets for evaluation, and surface-atmosphere integrated analysis method.Section 3 evaluates the overall performance of the modeled isotope values against various observations, focusing on the surface and vertical distribution.Section 4 quantifies model spreads and their underlying causes in different regions.Lastly, Section 5 summarizes and discusses our findings and insights gained from applying water vapor budgets using isotopic information.
In essence, this study addresses two crucial questions for a precise comparison that has not been previously attempted: 1. How can we validate modeled water isotopes and interpret the discrepancies against observed water isotopes?2. How do specific hydrological processes affect atmospheric isotopic ratios in simulations where different models or reanalyses are used for nudging?
To explore the answers to these questions, we conducted targeted experiments to test our hypothesis that the choice of model exerts a more significant influence on the uncertainty in water isotope simulations than the use of different reanalyses for forcing, particularly in key water cycle processes such as evaporation, horizontal moisture flux, and precipitation.Our hypothesis, combined with the methods introduced in this study, allows us to systematically assess both the uncertainty of isotope-enabled models and the associated spread in modeled water isotopes.

Model Description and Data Set
The model characteristics, reanalysis data sets utilized, and observation databases used for evaluating our simulations are summarized in Table 1.We provide a detailed description of these products below.
MIROC5-iso (Okazaki & Yoshimura, 2019) is an isotope-enabled atmospheric component of the Model for Interdisciplinary Research on Climate, Earth System Model (MIROC-ESM) (Watanabe et al., 2010) developed by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) in collaboration with the University of Tokyo and the National Institute for Environmental Studies.The model utilizes the Chikira-Sugiyama cumulus scheme (Chikira, 2004) and the entraining plume model (Chikira & Sugiyama, 2010) for convective parameterization.It also employs a semi-Lagrangian scheme for the tracer advection.Additionally, the land component, known as Minimal Advanced Treatments of Surface Interaction and RunOff (MATSIRO) (Nitta et al., 2014;Takata et al., 2003), from MIROC-ESM is used as the land surface model, and isotope tracers are implemented by Yoshimura et al. (2006).Sea surface and lake water isotopes were constant at 0‰, whereas sea ice δ 18 O and δD were constant at 3 and 20‰, respectively (Joussaume & Jouzel, 1993).The model has a resolution of the Note.Additionally, JRA55 is utilized for nudging across all three models: IsoGSM, MIROC5-iso, and ECHAM6-wiso.To validate the models' performance, water isotope observations in monthly precipitation (GNIP), hourly surface vapor (SWVID), multi-year Database of Antarctic snow isotopic composition, and daily vertical vapor (TES) are utilized.(Kanamitsu, Ebisuzaki, et al., 2002;Kanamitsu, Kumar, et al., 2002).

Water Isotope Observational Products and Climate Data Sets: GNIP, SWVID, Antarctic Snow, TES, CRU TS, and GPCP
The observation data sets used to evaluate the model results consisted of precipitation, snow, and water vapor measurements near the surface and atmosphere.The Global Network of Isotopes in Precipitation (GNIP) measures water isotopes globally every month (IAEA/ WMO, 2006).For this study, we utilized a part of 1,017 stations by selecting the sites where measurements had been taken for at least 12 consecutive months from 1979 to 2020.Seasonality was removed from the regional correlation analysis.The Stable Water Vapor Isotope Database (SWVID) was also used for near-surface isotope observations (Wei et al., 2019).This data set gathers isotope measurements made via infrared isotopic spectroscopy over the land surface and sea surface via ship and in the free troposphere via aircraft.The time resolution of the raw data is hourly.We converted these data to daily timesteps using an amount-weighted average for our analyses.To evaluate our simulations for the Antarctic area, we used the Database of Antarctic snow isotopic composition (Masson-Delmotte et al., 2008).The data set includes 1,279 data points: 1,125 locations for δ 18 O, 938 sites for δD, and 794 sites with data for both isotopes, enabling the calculation of the deuterium excess.
For the atmospheric variables, we compared our simulations with water isotope satellite data for water vapor isotopes in the atmosphere, as explained in Worden et al. (2006).This data is included in the Tropospheric Emission Spectrometer (TES-NASA)/Aura Level 3 Water Vapor Daily Gridded Version 5 data set.It has 15 BONG ET AL.
10.1029/2023JD038719 6 of 28 vertical levels, and this study considers three layers at 825.4, 681.29, and 464.16 hPa, known as the most sensitive levels for estimating the TES δD profile (Yoshimura et al., 2011).The modeled target variables (water vapor and its δD) were interpolated to the TES grids for model-data comparisons (three layers on the vertical, 4° × 2° grid on the horizontal).To remove the impacts of cloudiness (clear-sky bias) induced by poor estimation from the space-borne infrared sensor, all simulations and TES data sets were vertically filtered to remove rainy days according to the Global Precipitation Climatology Project (GPCP) v1.3 Daily Precipitation Analysis Climate Data Record (Adler et al., 2017).The TES data covers the period from 2004 to 2018.However, we note that 2008 was excluded from the analysis due to weak spatial data coverage.Additionally, because of the "Sun Synchronous Orbit" of the TES satellite, there are deficiencies in the water vapor isotope data set at 2-day intervals.
To evaluate modeled temperature and precipitation, which significantly affect isotope values, we used the climate data sets: the Climatic Research Unit's gridded Time Series (CRU TS) v4.06 for near-surface temperature (Harris et al., 2020), supported by an extensive network of weather station observations, and the GPCP v2.3 monthly satellite-gauge combined precipitation (Adler et al., 2016), which integrates satellite precipitation estimates from passive microwave and infrared, along with other low Earth orbit data and in-situ observations.

Experimental Design and Nudging
Simulated water isotopes differ depending on the models or the forcing conditions employed.Thus, two sets of model simulations were designed to make a rigorous comparison of the water isotope-enabled models.One uses a single reanalysis (JRA55) to nudge the different models (IsoGSM, MIROC5-iso, ECHAM6-wiso), and the other uses a single model (IsoGSM) nudged with different reanalyses (JRA55, NCEP-R2, ERA5).Using these two sets of experiments-Group 1 (IsoGSM.JRA55, MIROC5-iso.JRA55, ECHAM6-wiso.JRA55), which uses a single reanalysis to nudge three models, and Group 2 (IsoGSM.JRA55, IsoGSM.NCEP-R2, IsoGSM.ERA5), which uses multiple reanalyses to nudge one model-it is possible to separately diagnose (a) differences in model physics, and (b) differences in how the reanalysis product applied affects atmospheric dynamics.For nudging simulations, our water isotope-enabled models were nudged toward the chosen reanalysis data considering 3D-fields of temperature (for GSM and ECHAM, except MIROC) and atmospheric circulation (GSM and MIROC: zonal and meridional wind fields; ECHAM: vorticity and divergence), every 6 hr for our period of interest.However, humidity reanalysis data were not utilized for nudging during the simulation of water isotopes.Consequently, while the water isotopes were modeled based on reanalysis circulation information, the humidity in the simulations may vary due to the distinct schemes of each model, even though all three models were nudged by the same reanalysis circulation for the Group 1 experiment.Typically, it takes up to a 1-month model period to reconstruct a reliable atmospheric distribution of water vapor isotopes, regardless of the initial conditions.The sea surface temperature and sea ice fields (GSM: daily mean, MIROC: 6-hourly, ECHAM: monthly mean) from the corresponding reanalysis data set were applied as sea surface conditions for all our experiments.

Isotope Mass Budget Analysis
This study introduces a new methodology, water isotopic change rates, to quantitatively decompose all transformations of the water isotope ratios in meteoric water.This method quantifies isotope value changes in the atmosphere by considering the gradient of isotope values between the atmospheric column and the surroundings with a weighting amount of water exchange.First, the atmospheric water balance equation describes the change in the water amount (Oki et al., 1995;Väisänen, 1961).(2) Here where F in and F out are the horizontal incoming and outgoing moisture flux components in the convergence, respectively.Second, the water balance equation for the gradient of isotope values is multiplied by the surrounding isotope ratios (δe t , δp t ,  in  , and  out  ) accompanying the water exchange to approximate the isotopic mass balance form in Equation 7 (Hayes, 2004;Yoshimura et al., 2003), and the isotope ratio in the column in Equation 8.In this study, δw refers to the oxygen isotope value in the vertically mass-weighted sum of the water vapor in the atmosphere.

+∆
where δe t , δp t ,  in  ,  out  are isotope values in evaporation, precipitation, and the ambient isotope ratio of oxygen entering and exiting through horizontal moisture flux, respectively.By re-organizing Equations 7 and 8 (Krabbenhoft et al., 1990), we defined ∆δw t , isotopic change rate in a column as follows: where ∆δw e , ∆δw p , ∆δw q are the isotopic change rates derived from the evaporation, precipitation, and horizontal moisture flux, respectively.
Please note that water amounts ( 10.1029/2023JD038719 8 of 28 tionation factors for each process are constant during 6 hr.However, during our preliminary testing prior to fully adopting the method, we confirmed that the conservation of mass for water isotopes becomes less reliable as the time interval expands, due to the inability to account for changes in the fractionation factor over extended periods.This limitation presents challenges for validation against low-time-resolution satellite water vapor isotope data sets, as the isotopic change rate (‰/s) is currently only calculable through isotope-enabled models.
As we assumed in this study that water and heavy water isotope masses are conserved, we can define contributions to ∆δw in the following manner: where C δe , C δp , and C δq (instead of C δX ) are the relative contributions of evaporation (∆δw e ), precipitation (∆δw p ), and horizontal moisture flux (∆δw q ), respectively, to the total isotopic change rate in a column, ∆δw.This method takes absolute values to consider all processes, regardless of enrichment or depletion of water isotopes (see Text S1: Method details in Supporting Information S1).

Mean Conditions and Global Distribution of Water Isotopes
The annual mean of model ensembles of Group 1-three models nudged toward JRA55-presents temperature (Figure 1a), precipitation (Figure 1b), evaporation (Figure 1c), and near-surface humidity (Figure 1d).Figures 1e-1h illustrate their variables' spread in the models.When constraining sea surface temperature and sea ice concentration to a single reanalysis, simulations over the ocean exhibit a significantly smaller surface temperature spread (as indicated by the standard deviation of Group 1) compared to the land (Figure 1e).On the land surface, even when the same atmospheric temperature and circulation are applied to the models (except for MIROC5, which differs with respect to nudging vertical temperature), the different land components of the three models can cause variations in radiative energy and moisture exchange.This can enlarge the spread of modeled skin temperature.Particularly in Antarctica, where the spread of modeled temperature is most prominent, various model schemes for cloud physics and snow and resolutions (Tewari et al., 2022), which affect the topography in the models, could induce different temperatures due to the different reflections of incoming solar radiation and outgoing ocean thermal radiation.
Precipitation spread is most noticeable in the tropics (Figure 1f), where warm surface conditions trigger uplifting air and maintain low-level convergence.One possible cause may be the different moist convection schemes influencing cloud formation.Another factor could be the different strengths of nudging among the models.These differences result in different low-level moisture convergences for precipitation, especially in the Intertropical Convergence Zone (ITCZ).Such differences in model characteristics amplify the spread of rainfall and the depth of convective clouds in the subtropical North and South Pacific and Atlantic.Moreover, land models and monsoons can impact the spread of precipitation in Amazonia and even more significantly in Central Africa.This is because the Atlantic and Indian oceans supply differing amounts of water vapor depending on the season.
In the Amazonia and Central Africa regions, the amount of evaporation varies depending on the model considered (Figure 1g).Meanwhile, processes like evapotranspiration from plants and re-evaporation from the soil, which influence moisture flux release into the atmosphere, are defined differently across our models.IsoGSM and ECHAM6-wiso do not account for fractionation processes during evaporation over the land surface, whereas MIROC5-iso incorporates kinetic fractionation for surface water evaporation.Over the ocean surface, the modeled surface moisture and drag coefficient may contribute to the spread of evaporation, but they are much smaller contributors compared to the spread of precipitation.
The spread of surface layer humidity, or the humidity in the bottom-most layer of the atmosphere in a hybrid sigma-pressure vertical coordinate (Figure 1h), can be explained by multiple spreads of other variables (as shown BONG ET AL. 10.1029/2023JD038719 9 of 28 in Figures 1e-1g).Since precipitation moistens soil, the spread of modeled surface humidity in Central Africa correlates with the spread of rainfall.However, the spreads of modeled temperature, which arise from the different land schemes in the grasslands east of the Andes Mountains, can increase the range of water vapor the near-surface air can contain.As a result, these spread patterns in four key variables (temperature, precipitation, evaporation, and near-surface humidity) play a role in explaining the differences in modeled surface water When investigating the disparities in the modeled water cycle, examining arid regions such as deserts or high-latitude areas has limitations in measuring these spreads of modeled water amounts at the surface (Figures 1f-1h) and in the atmosphere.In this context, the relative amount of heavy water, δ 18 O in precipitation, serves as a more versatile and effective measure for evaluating the modeled water cycle.This is because δ 18 O, determined by fractionation processes, is more sensitive to phase change processes than changes in typical water amount, especially in extreme climate conditions like very dry or cold environments.
As shown in Figure 2a (from Group 1), lighter water (lower value of δ 18 O) travels further from the ocean to inland because heavy water isotopes preferentially condense early in the atmosphere.In the subtropical high-pressure regions that span both oceanic and terrestrial zones, descending air motions result in less precipitation (Figure 1b), offering fewer opportunities for the depletion of atmospheric vapor isotopes by rainfall.However, the isotopic composition of precipitation may differ depending on surface conditions.Over the oceans, heavy isotopes are supplied to the atmosphere through evaporation; hence, when precipitation occurs, it contains a higher concentration of heavy isotopes without previous significant depletion.In contrast, in the desert regions of Africa and Australia, where heavy isotopes are not supplied to the atmosphere via surface evaporation, a pronounced continent effect emerges, leading to precipitation with reduced isotope values.
In the Maritime Continents, the amount effect, which is sensitively tied to deep convection (Tharammal et al., 2017), depletes heavy water.The well-known latitude effect also occurs as moisture moves to cold, high latitudes; its isotope ratio decreases progressively due to Rayleigh fractionation processes (Clark & Fritz, 1997).
In regions where these fractionation effects (Dansgaard, 1964) are strong or overlap, such as in high mountains or dry areas like the Sahara Desert, Antarctica, and Greenland (Figure 2e), the spread of modeled δ 18 O in precipitation is significant.This is because different models have different cloud physics or convection schemes.However, patterns of spread in δ 18 O in precipitation and precipitation amounts (Figure 1f) differ.This is because water vapor δ 18 O in the atmosphere can be enriched by re-evaporation or affected by different vapor sources while being modeled in the atmospheric circulation system.Notably, the spread of δ 18 O in the atmosphere nudged by different reanalyses (Group 2, Figure 2f) is much smaller than the impact from models (Group 1, Figure 2e).
The second-order parameter in isotopic analysis, known as deuterium excess (abbreviated as d-excess), reflects the characteristics of water source environments (Figures 2c and 2d).It is determined by the local relationship between δD and δ 18 O.Cold, dry, and land-sourced water is known to increase the kinetic fractionation (in humid-deficit non-equilibrium conditions) of evaporation and d-excess values (Craig, 1961;Harmon & Schwarcz, 1981).The modeled d-excess value is notably diverse in the Sahara Desert, irrespective of the model or reanalysis nudging choice, mainly due to deficient precipitation (Figures 2g and 2h).Antarctica is the region with the highest spread of modeled isotope values.Intrinsically, uncertainties in modeled surface temperatures contribute to this spread in δ 18 O and d-excess.Moreover, Antarctica's vast, cold, and dry land conditions may make water isotope simulations more sensitive and variable.Conversely, the spreads of δ 18 O and d-excess in the precipitation according to the reanalysis choice are much less than the spreads among models (standard deviation for models δ 18 O: 3.8‰ and d-excess: 11.7‰, and for reanalyses δ 18 O: 0.9‰ and d-excess: 1.6‰), as shown in Figures 2e and 2g (Group 1) and Figures 2f and 2h (Group 2).Despite the small amount of precipitation in Antarctica, different horizontal resolutions and a lack of polar-specific physics result in precipitation biases in GCMs (Genthon et al., 2009;Tewari et al., 2022), leading to a greater spread of isotope values than the reanalysis forcing difference.

Evaluation of Water Isotopes in Precipitation and Surface Water Vapor From Monthly to Daily Time Scales
The model's capability to simulate water isotopes was evaluated using both in situ and satellite data.The Taylor diagram (Figure 3) presents the model-data correlation for monthly variations in temperature, precipitation, δ 18 O, and d-excess at the GNIP sites.Points closer to the reference (REF) data indicate higher model-data correlations and more similar standard deviations between the modeled and observed signals.
The nudged temperature showed a high correlation of about 0.9, and the variance of the modeled precipitation was closest to the observations (standardized deviation ∼1.0 in the Taylor diagram).For water isotope values in precipitation, δ 18 O tended to follow the precipitation performance; however, the correlation was lower by about 0.1, and its variance was slightly larger (standardized deviation above 1.0) than precipitation.The d-excess parameter presented challenges due to its significantly smaller correlation and larger variance than other variables.This disparity arises from the complexities in accounting for kinetic fractionation during evaporation processes.Moreover, excluding this kinetic process might result in a more unrealistic d-excess during water evaporation from land or when it re-evaporates in the atmosphere.To improve the model-data correlation for isotope values in precipitation, one approach is to create a model ensemble (δ 18 O, δD, d-excess in Table 2 and Table S1 in Supporting Information S1) using the same reanalysis for nudging (JRA55), which corresponds to Group 1 (IsoGSM.JRA55, MIROC5-iso.JRA55, ECHAM6-wiso.JRA55).In the model ensemble (Mo ENS), precipitation (GNIP) and surface water vapor (SWVID) isotope observations were always better correlated than with the individual models.In the case of δ 18 O, the global mean correlation value with GNIP data (0.601) was similar to that of the SWVID data set (0.570).One aspect to note is that the GNIP sites are densely located in the European region (Figure 4), so correlations in SWVID data cannot be directly compared due to the differences in time resolution and seasonality across various areas.
When considering the same model (Group 2 using IsoGSM), the correlation between observation and simulation nudged to ERA5 was higher than that achieved using other reanalyses for nudging.On a regional scale, selecting the appropriate reanalysis data set for nudging becomes crucial and depends on the targeted area (Table S2 in Supporting Information S1).For instance, between 60° north and the pole, the dynamics of JRA55, such as temperature and wind fields, correlate best with isotopic observations.Specifically, between 30° and 60° north, JRA55 for nudging is the most suitable for explaining temporal variation in Eastern North America (0.588) and Western Eurasia (0.396).However, its performance declines in Eastern Eurasia (0.275) due to the influence of a strong, wavering jet stream, which leads to highly variable precipitation in the mid-latitudes of the Northern Hemisphere.In Antarctica, the worst reanalysis choice (NCEP-R2) for nudging IsoGSM, in terms of the evaluation of modeled water isotopes at the global scale, resulted in the highest model-data correlations (0.930) and regression slopes (0.74) closest to 1 with the Antarctic snow data set.Thus, a better regional-scale analysis for water isotope studies requires careful selection of the nudging reanalysis data set.

Mean Conditions and Atmospheric Moisture Cycling Between Ocean and Land
In the water cycle, water vapor absorbs surplus solar heat and is transported from the tropics toward the poles through meridional eddy moisture transport, and from the ocean to land, releasing heat into energy-deficient regions with variations depending on the hemisphere (Figure 5a).Processes that result in fractionation, such as the amount effect in the Maritime Continents, the continental effect from seacoasts to inland areas, and the altitude effect in mountains and plateaus, influence the isotopic composition of water vapor (Figure 5b).
In all oceans, a subtropical high-pressure system increases evaporation minus precipitation by surface winds and inhibits the uplift of air mass.This keeps heavy water in the atmosphere without being removed by rainfall.Also, western boundary currents play a crucial role in shaping the hemispheric contrast of the water isotope distribution.These narrow and fast meridional currents transports heat energy poleward and release energy through latent heat and water vapor into the atmosphere that becomes dry from mid to high latitudes.In a dry atmosphere above the ocean, particularly where the air sinks in the mid-latitudes, the Northern Hemisphere ends up being wetter than the Southern Hemisphere, highlighting different roles in water transportation toward the poles.
As the water moves from the humid equator toward drier high latitudes, the observed water vapor concentration and its isotopic composition vary between the Northern and Southern Hemispheres due to differences in the spatial coverage of land (Figure 5c).In latitudes between 20° and 60°, the Northern Hemisphere has more land coverage, which places the source of water vapor further away.In contrast, the expansive Southern Ocean supplies abundant water to the land and drives a more extensive water cycle compared to the Northern Hemisphere.As a result, the atmosphere over the northern continents is significantly drier than that over the southern continents.BONG ET AL. 10.1029/2023JD038719 13 of 28 At the same time, depletion of δD becomes more sensitive to humidity decreases between land and ocean in the higher latitudes because the remaining vapor in the cloud becomes colder, thereby enhancing fractionation through the Rayleigh process (Clark & Fritz, 1997).The models (Figures 5d and 5e) underestimate the amount of atmospheric water vapor compared to the observations (Figure 5c), being 0.69-1.45mmol/mol lower between 20° and 60° latitudes, and 2.45-3.05mmol/mol lower in the tropics.However, the models perform well in simulating the hemispheric contrasts of the water vapor concentrations and their isotope ratios.We find significant differences in model spread (measured as the standard deviation of scattered points) in the land-ocean moistening and dehydrating processes on water vapor and δD: 0.16 mmol/mol and 5.64‰ standard deviation over land and 0.22 mmol/mol and 3.75‰ over the ocean.This implies that the modeled heavy water isotope values also exhibit non-negligible uncertainty, particularly when considering the atmosphere over land.

Evaluation of the Atmospheric Isotopic Composition of Water Vapor
We considered the mid-troposphere, which includes the most sensitive levels (825, 681, and 464 hPa) in estimating the TES δD profile (Yoshimura et al., 2011) containing a substantial portion of water vapor in the vertical profile.
The temporal correlation of vertically integrated daily δD between observations and simulations is positive from 60°S to 60°N latitude (Figure 6), and simulations nudged to JRA55 and ERA5 perform better (Table S3 in Supporting Information S1 and Figure 7) over the ocean and land, respectively.Overall, these patterns become increasingly apparent as the altitude rises to 464 hPa.The model ensembles of δD come closer to observation (Figures 7b-7d), and the root mean square deviation (RMSD) decreases significantly from 54.94‰ at 825 hPa to 24.36‰ at 464 hPa.
Also, the overall vertically integrated δD values, when considering the three levels' amount of water vapor, reproduce a better time variation of the TES δD than when considering only single levels in the daily climatology of the spatial mean from 60°S to 60°N latitude (correlation = 0.894 in Figure S2a in Supporting Information S1,  BONG ET AL.

10.1029/2023JD038719
14 of 28 model-data slope = 0.99 in Table S3 in Supporting Information S1).Integrating multiple vertical levels not only offsets the bias but also reduces the negative correlation value in high-latitude regions.
Low correlations are closely related to the temporal variations in humidity.In the tropics, the modeled vertical humidity along the ITCZ seasonal movement range was drier than that in the TES observations (Figure S3b in Supporting Information S1).It becomes even drier from mid to high latitudes due to the cold bias in   BONG ET AL.

10.1029/2023JD038719
16 of 28 the winter hemisphere (Figure S3c in Supporting Information S1).The dryness in the tropics significantly affects the most humid level, 825 hPa, and results in a very large RMSD (70.4‰) and a very low correlation (0.24) in Figure S2c in Supporting Information S1, as seen in Figure 6.The strong cold bias, which causes substantial depletion of δD by decreasing cloud temperature and reducing vapor amount in Figure S3a in Supporting Information S1 in the high latitudes of the winter hemisphere, reduces water vapor by approximately 20%-60% in Figure S3d in Supporting Information S1, especially in Antarctica, and leads to negative correlations.
The difference in the heat capacity between land and ocean influences the correlation of the vertical daily δD with the observations.As humidity over land varies more because of the more extensive changes in daily and seasonal surface temperatures, the temporal correlation is higher in the Northern Hemisphere, which has a more significant proportion of land-covered area.As shown in Figure 6 and Table S3 in Supporting Information S1, the correlation coefficients increased in the order of tropics, Southern Hemisphere, and Northern Hemisphere; correlation over land was higher than over the ocean.As a result, vertically integrated daily δD shows the highest correlation with observations (0.425) in the northern continent and the least in the tropical ocean (0.135).While the land-rich Northern Hemisphere exhibits clearer seasonality and correlates more strongly with observed δD (0.95) than the ocean-rich Southern Hemisphere (0.78), equatorial regions in the model ensemble were found To explain the series of heavy water isotope transports in the water cycle, we utilized an atmospheric water budget analysis to investigate and decompose process-based fractionations.Unlike traditional moisture transport within the water cycle, our method for measuring water isotopic change rate (∆δw = ∆δ 18 O‰/day in the atmosphere) can help quantify fractionation processes and their proportional contributions.Furthermore, it helps establish the extent to which atmospheric heavy water isotopes are enriched or depleted after evaporation, horizontal moisture flux, and precipitation.

Fractionations Derived by Evaporation, Precipitation, and Horizontal Moisture Flux
Evaporation-derived fractionation enriches the atmosphere by supplying heavy water through the surface-air interface (Figure 8a).This is more pronounced in high-latitude oceans (south: 3.49‰/day, north: 2.81‰/day,  S4 in Supporting Information S1), indicating that the gradient of isotope values between the ocean and atmosphere is more prominent, even though the evaporation amount is smaller than in low-latitude oceans.Moreover, the dehydrated marine atmospheric boundary layer in the cold, high-latitude atmospheres is more strongly influenced by warm, moist western boundary currents (Figure 5a).
Other atmospheric fractionation processes related to precipitation, such as the continental effect globally and the altitude effect in high mountains and the Tibetan Plateau, have also been simulated.In tropical oceans, rainfall generates a large amount of precipitable water at the surface (Figure 1b).However, the amount effect is less significant than the other effects mentioned above (average of four tropical regions in India and Pacific Oceans: −1.82‰/day, compared to North Atlantic Current: −5.06‰/day and Antarctic Circumpolar Current: −6.57‰/ day).A lower  ∆  implies a smaller gradient of isotope values between surface precipitation (δp t in Figure 2a) and atmosphere (δw t in Figure 5b) in the tropics.This is explained by precipitation-derived fractionation term,  ∆  = −( − )∕( + ∆→+∆) , as seen in Equation 9of the isotope mass budget analysis method.
In other words, moisture recycling in the tropics occurs more rapidly than in high-latitude regions due to shorter atmospheric water vapor residence time (Gimeno et al., 2021).Consequently, there are fewer chances for the water isotopes to be depleted by precipitation in the tropics.Meanwhile, precipitation-derived highly depleted regions broadly correspond to where moisture converges, such as in Western North America, Tibet, Greenland, and East Antarctica (Figure 8c).
The horizontal gradient of the atmospheric isotope value affects the isotopic change rate induced by fractionation resulting from horizontal moisture transport.The global annual mean of isotopic change rate, which is derived from horizontal moisture flux, is positive, with an average value of  ∆ = 1.01‰∕day .This suggests that atmospheric circulation globally enriches water vapor and alleviates the imbalance between evaporation-derived  ∆ = 2.19‰∕day and precipitation-derived  ∆ = −3.13‰∕dayisotope change rates (Figure 9 and Table S5 in Supporting Information S1).In other words, the impact of circulation on changes in atmospheric water isotope values is more closely aligned with precipitation's depleting effects rather than evaporation's enriching effects.Indeed, regions that are enriched by circulation tend to correspond with areas that are depleted by precipitation.

Understanding Atmospheric Water Isotope Residuals
Taken together, considering that the atmospheric oxygen isotope residual is small (  ∆ = ∆ + ∆ + ∆ = 0.07 ‰∕day , which is equivalent to 1.1% of  |∆| + |∆| + |∆| ), we conclude that isotopic change rates introduced in this study explain a significant feature of large-scale heavy water isotope cycle.As a result, this study has successfully demonstrated that the mass of heavy water remains globally balanced, as indicated by no residual (white bar) in Figure 9a, while providing global estimates for changes in isotopic water rates.
Another approach for understanding the heavy water isotope cycle is to calculate the contributions of the isotopic change rate (Figures 8d-8f).Taking the absolute values of positive enrichment and negative depletion rates and expressing them proportionally in percentages can clarify the contributions of changes in isotope influx and outflux derived by the three different processes.In turn, we can determine which process predominantly determines water isotope values in various regions.The horizontal moisture flux contributes   = 50.0% of changes in atmospheric isotope value globally, and the contribution of evaporation  = 27.1 % is 4% higher than that of precipitation   = 22.9% (Table S5 in Supporting Information S1).This is because the evaporation impact occurs over a broad subtropical region (Figure 8d), whereas precipitation mainly controls the isotopic change in the narrow tropical convergence zone (Figure 8e).High-latitude continental water isotopes are more sensitive to horizontal moisture flux (Figure 8f).This accounts for 90% of the isotopic change signal in the vast Sahara Desert, where both precipitation and evaporation are nearly zero.
Strong positive residuals in Greenland and East Antarctica are attributed to similar contributions from horizontal moisture flux (∼68%) and precipitation (∼22%), which play roles in enrichment and depletion, respectively (Figures 9a and 9b; details in Table S5 in Supporting Information S1).In other words, the positive residuals (white bar) suggest that atmospheric condensation for rain (with its depleting role) does not sufficiently offset the enrichment caused by horizontal moisture flux.However, the mechanisms influencing the residuals in these BONG ET AL.

10.1029/2023JD038719
19 of 28 regions differ.The isotopic depletion and enrichment become more pronounced under cold and dry conditions, which are prevalent in these two regions.Moreover, environmental factors such as temperature or relative humidity during evaporation, as well as the distance traveled through dry air, affect precipitation and result in variations in the extent of fractionation changes.In Greenland, the warm North Atlantic Drift transports heat energy to high latitudes, and the newly evaporated oceanic water establishes a more significant horizontal gradient in isotope values.In contrast, in the vast expanse of East Antarctica, fractionations due to phase changes diminish signif- 10.1029/2023JD038719 20 of 28 icantly from the seashore, as water vapor is more likely to condense into precipitation before being transported further inland.Consequently, while the contributions from three isotopic change rates are similar, the rates in Greenland exceed those in Antarctica (details in Table S4 in Supporting Information S1).The positive residuals in the atmosphere also imply that a surplus of isotope residuals can be effectively removed by incorporating more heavy water isotopes once the atmospheric vapor condenses and precipitates.These characteristics have implications for isotope signals recorded in ice cores for past climate studies, especially in Antarctica and Greenland, where annual rainfall is scarce.
In the mid-latitudes, we found substantial negative residuals (white bar in Figure 9a) in Central and Eastern North America and Far East Asia.In these regions, the precipitation-derived isotopic change rate (red bar in Figure 9a) is approximately −3.3‰/day, contributing to about 19.1% of the isotope signal (details in Table S4 in Supporting Information S1).While evaporation depletes (blue bar in Figure 9a), horizontal moisture flux can dampen or enhance atmospheric water isotope differently (yellow bar in Figure 9a).In these regions, it is noteworthy that horizontal moisture flux and precipitation both contribute to decreasing water isotope values in concert (red and yellow bars are all negative), unlike in scenarios where the mechanisms are opposite, such as when increasing moisture due to air convergence leads to decreasing the values due to subsequent rainfall (negative red bar while positive yellow bar).Also, the more significant the contribution of evaporation to the isotopic change rate (Figure 9b), the more horizontal transport depletes heavy water isotopes (Figure 9a).To account for the negative residuals and characteristics observed in these mid-latitude regions, potential mechanisms may involve oceanic vapor sources: westerlies transporting vapor from the ocean to western continents, the polar vortex carrying cold-origin vapor from the Arctic Ocean, and monsoons on the eastern coasts of continents.
In Central North America, the prevailing westerly carries a stream of Pacific Oceanic vapor eastward, enriching water vapor over the western continents and windward coastal mountains (Figures 8b and 8c).Following rainfall, the depleted water vapor isotopes move further east.In Eastern North America and Far East Asia, in winter, a polar vortex allows low-isotope value water originating from the Arctic Ocean to flow into the mid-latitudes.In summer, vapor evaporated from hot and humid western boundary currents, such as the Gulf Stream and Kuroshio Current (Figure 8a), is carried to the east coast under subtropical high pressure (Figure 8c).As a result, the summer monsoon advances northward, bringing increased water vapor isotopes from the evaporation-derived enriched regions to the continents and depleting the atmosphere through horizontal moisture flux-derived fractionation.Similarly, additional water isotope studies can be conducted in various areas to explore different weather and climate change scenarios using this decomposed budget analysis beyond simply focusing on the annual mean analysis.

Exploring Water Source Conditions Through D-Excess Analysis
As a further step, changes in the water source conditions can also be explored by applying our method to d-excess = δD − 8 × δ 18 O (Figure S5 in Supporting Information S1).Unlike the isotopic change rates for δ 18 O, evaporation lower d-excess in the atmosphere (global d-excess ∼10‰) by supplying vapor with lower d-excess values originating from the ocean (which are close to zero ‰), as shown in Figure S5a in Supporting Information S1.Similarly, the transport of ocean-oriented water vapor tends to decrease the d-excess over continents, especially in regions where moisture is blocked and converged, such as in the windward coastal Rockies, the Andes Mountains, and Tibet (Figure S5c in Supporting Information S1).All decomposition processes intensify at higher latitudes, owing to increased kinetic fractionation caused by the slower diffusion of H 2 18 O compared to HD 16 O at lower relative humidity.This results in an excess of D relative to 18 O.As a consequence, in dry atmospheric conditions where the slope is lower than 8 (the global meteoric water line) in xy scatter plot (with x representing δ 18 O and y representing δD), precipitation-derived depletion causes an increase in the change rate of d-excess at higher latitudes (Figure S5b in Supporting Information S1).Also, it's noteworthy that the contributions from precipitation-derived change rates for d-excess (12.6%) are reduced to half of those for δ 18 O (22.9%).This phenomenon arises because equilibrium fractionation becomes dominant during condensation under 100% relative humidity conditions.Specifically, the average contributions are as follows:  -excess  = 23.2% ,  -excess  = 12.6% , and  -excess  = 64.2%, while for δ 18 O, they are

Process-Based Quantification of Uncertainty in Water Isotope Models
Quantifying the uncertainty, which serves as the cause for the spread of modeled variables, provides information about how different water isotopes can be simulated depending on the AGCM used.The modeled spreads of isotopic change rates (represented in brown in Figure 10) are decomposed according to evaporation, precipitation, and horizontal moisture flux, as shown in Figures 11a-11c, respectively.
The moisture flux at high latitudes and mountains accounts for 57% of the total model spread, while precipitation, predominantly on land, accounts for 25% (Figure 10a).Although the same reanalysis is nudged to the models (Group 1), the modeled wind fields carrying moisture flux can vary slightly due to model biases.As such, it is difficult to strictly identify the causes of the spread of moisture flux-and precipitation-derived fractionations.However, they may be influenced by model structural errors caused by schemes for atmospheric transport, convection, or cloud physics.The spreads due to the choice of reanalysis are much smaller than that caused by the model choice (with a ratio of 1:5.6).Although precipitation accounts for 60% (global area average: 0.18‰/ day) of the total spread of the three processes (0.30‰/day) resulting from different reanalyses (Group 2), it is smaller than the minor cause, which is evaporation-derived fractionation (0.29‰/day) accounting for 17% of the total spread from the models (1.68‰/day).Upon closer inspection of this global average locally, by averaging the three spreads of processes (Figure 10b), the influence of the model choice is more pronounced in most regions.One exceptional area is the eastern tropical Pacific, where the choice of reanalysis has a more significant impact.Generally, spreads are more prominent where the magnitudes of isotopic change rates are high, such as in Greenland and East Antarctica, at high latitudes.Similar to Tibet and Western North America, some regions in the mid-latitudes also exhibit more significant uncertainty.In such cases, our method can help identify which isotopic fractionation process contributes to higher uncertainty in those regions through regional decomposition analysis (Figures 11d-11f).
In the eastern tropical Pacific region, where air-sea coupling influences El Niño-Southern Oscillation (ENSO), the cause of a larger spread due to reanalysis choice, compared to the models in Figure 10b, is the spread of precipitation-derived fractionation displayed in Figure 11e.This spread (specific values are in Table S4 in Supporting Information S1) is nearly three times larger (0.47‰/day) than that caused by different models (0.16‰/day).Because the physics of the model remains consistent within a single model, reanalyses constrain the large-scale atmospheric circulations differently, affecting factors such as the positions of the ITCZ.This results in differences in atmospheric depletion by precipitation.Therefore, careful attention is required when conducting isotope studies using nudging simulations in the ENSO region, especially when determining precipitation conditions associated with stratus clouds over this area.
Additionally, the Antarctic Circumpolar Current exhibits a notable spread in evaporation-derived fractionation (Figure 11d).This spread is more influenced by variations in sea ice concentration nudging over open seawater than by the choice of model.Specifically, the spread measures 0.67‰/day in simulations from different models in Group 1 and 0.81‰/day in simulations from different reanalyses in Group 2. In the case of Tibet and Western North America, horizontal moisture flux is the most crucial factor contributing to model spread for simulations involving atmospheric water vapor isotopes.Before water isotopes were transported and converged due to atmospheric convergence over Tibet, and before they were transported and converged through atmospheric rivers into Western North America, variations in simulated and transported water vapor within the models (Group 1) can amplify both the spread of precipitable water and the uncertainty in isotopic change rates.

Summary and Conclusion
Although numerous comparative studies have been conducted to date, this study presents a unique approach that compared and validated the impact of both models and nudging data on modeled surface and atmospheric water isotopes.Furthermore, it quantified simulation uncertainties in terms of isotopic mass conservation.
With validation against global surface observations, our findings indicate that correlations of modeled isotopes in precipitation and near-surface water vapor are 0.60 for δ 18 O, while correlations for d-excess are considerably smaller (0.26).Notably, a discernible trend emerged: better performance in simulating temperature, precipitation, and humidity corresponds to higher correlations in modeled water isotope values with observations.However, it is imperative to acknowledge that optimal performance does not always align with specific models or reanalyses for nudging.Consequently, in order to enhance the performance of water cycle simulation, selecting an opti-  Another approach proposed to improve simulation performance in this multi-model study involves the calculation of averages.When applied to surface water isotopes such as precipitation and near-surface vapor, this method reduces susceptibility to errors by addressing outliers that arise when one model's simulation significantly deviates from the others.Hence, the ensemble of models provides a more comprehensive explanation for surface water isotopes, underscoring the importance of considering multiple water isotope-enabled models.Additionally, by vertically integrating water vapor isotopes, a more robust correlation emerges in the seasonal cycle of the modeled isotope values compared to considering individual altitude layers.Furthermore, in addition to averaging results or selecting a combination of a certain model and forcing data, gaining a profound understanding of the uncertainty associated with modeling water isotopes in a specific region of interest remains crucial.In light of this, we propose a method that involves examining the change rates of isotope values (expressed in ‰ per day) for δ 18 O and d-excess in the atmosphere.This method facilitates the quantitative breakdown of fractionations during the water cycle, spanning from evaporation to precipitation, thereby enabling the identification of processes contributing to the intermodel spread of modeled water isotopes.Consequently, we conclude that the spread in global simulations of modeled water isotopes, which results from model and reanalysis choices, is primarily driven by modeled moisture flux and precipitation uncertainty.These uncertainties originate from distinct model structures and dynamics, specifically atmospheric circulation.
Importantly, this method enables the decomposition of changes in atmospheric water isotopes into hydrological processes, a task that has been challenging due to the limitations of satellite observations.Through this method, the global estimates presented in this study reveal that surface evaporation (which enriches the atmosphere) and atmospheric precipitation (which depletes the atmosphere) do not offset each other's isotopic changes (specifically, depletion by precipitation outweighs enrichment by evaporation) in the atmosphere.This offers a unique perspective on water isotopic fractionation, deviating from the conventional understanding that evaporation and precipitation are balanced in the natural water cycle.Notably, our budget method explains a novel mechanism, viewed from a water isotope perspective, that leads to an enrichment in the isotopic value of atmospheric water through atmospheric circulation, addressing the disparity between evaporation and precipitation.Our findings have the potential not only to estimate how specific processes have influenced atmospheric water isotopes in response to climate changes but also to analyze isotopic variations in phase change in water.From this perspective, our budget analysis utilizes a process-based isotopic mass balance approach to provide a helpful tool for investigating isotopic signals from various regional to global scales, ranging from short-term weather to longterm climate timescales, and for measuring their uncertainty.
Historically, nearly four decades have passed since the inception of the numerical model for globally simulating atmospheric water isotopes (Joussaume et al., 1984).Recent efforts by various modeling groups aim to elucidate novel and intricate water isotope processes across different time periods, achieved by coupling existing atmospheric models with the ocean or by incorporating isotopic processes on land (Brady et al., 2019;Bühler et al., 2021;Cauquoin et al., 2019;Parker et al., 2021).These endeavors could be further enhanced by considering factors such as soil moisture (Zhou et al., 2021), plants (Beyer et al., 2020), surface water (Marković et al., 2020), and anthropogenic combustion-derived vapor (Xing et al., 2020).Also, water isotope simulations play a practical role in improving model performance by aiding in the determination of physics-related parameters for cloud and convection processes (Ramos et al., 2022).In this regard, the proposition put forth by the US Climate Variability and Predictability Program (CLIVAR) report (Bailey et al., 2021), as well as the perspective shared by S. G. Dee et al. (2023), underscores the significant benefits of incorporating water isotope physics into the physics of Intergovernmental Panel on Climate Change (IPCC)-class GCMs.This integration stands to substantially advance the model, particularly by introducing an observable tracer, water isotopes, that can effectively constrain model physics.Thus, in the context of modeling collaborations, such as the Coupled Model Intercomparison Project or Paleoclimate Modeling Intercomparison Project, the project SWING2, which provides results from water isotope models, has been very useful and important in terms of explaining a large number of observations, enabling the selection of the most appropriate modeling results, and presenting a wide range of modeling outputs to overcome model uncertainties.
Nevertheless, it's important to acknowledge that comprehensive discussions with diverse modeling groups are lacking, and a consensus has yet to be established for a collaborative framework design.Such a design, tentatively termed iso-CMIP or iso-PMIP, can pave the way for future collaborative efforts to improve climate models and analyze the water cycle using the water isotope, a tracer sensitive to changes such as temperature and humidity.Although this study lays the groundwork necessary for a full-scale SWING3, it does not encompass all the existing water isotope-enabled climate models, which include ECHAM (Hoffmann et al., 1998), GENESIS (Mathieu et al., 2002), GISS ModelE (Schmidt et al., 2005), REMOiso (Sturm et al., 2005), DHARMA (Smith et al., 2006), IsoGSM (Yoshimura et al., 2008) HadCM (Tindall et al., 2009), LMDZ-iso (Risi et al., 2010), Isotope-enabled SAM (Blossey et al., 2010), MIROC (Kurita et al., 2011), COSMOiso (Pfahl et al., 2012), ISOLES (X. Lee et al., 2012), UVic (Brennan et al., 2012), iLOVECLIM (Roche, 2013), SPEEDY-IER (S.Dee et al., 2015), Isotope-enabled WRF (Moore et al., 2016), MPI-ESM-wiso (Cauquoin et al., 2019), iCESM (Brady et al., 2019).BONG ET AL.Moving forward, it would be advantageous for future full-scale SWING3 projects to prioritize uniform experimental designs, similar to those introduced in this study, by using the same input and experimental framework.Such standardization could enhance both the rigor and interpretability of inter-model comparisons, especially when addressing the challenges of historical climate reconstructions and forward-looking climate simulations.A focused examination of the detailed differences in the models, particularly in their physical components, is also crucial.This would refine our understanding of the uncertainties among models but would also contribute significantly to the robustness of SWING3 project outcomes.Additionally, the budget analysis introduced in this study will help clarify the uncertainty associated with modeled water isotopes.
Particularly, amidst the growing impact of global warming, water isotope modeling emerges as a crucial tool in addressing responses to extreme weather events in the climate change crisis.Moving beyond traditional water isotope modeling for past reconstructions and proxies, studies that focus on understanding water isotopic behaviors under high-moisture conditions, such as during flood events or atmospheric extremes, offer vital insights for disaster preparedness.Through a collaborative, multi-model approach, nudging simulations can tightly constrain atmospheric circulation by leveraging the same present-day reanalysis data, reflecting recent temperature increases.This approach facilitates a more rigorous comparison of hydrological and physical processes simulated in each model, leading to a better grasp of their strengths and weaknesses in simulating water isotopic signals, which are sensitively responded to by changes in climate and weather.The enhanced comprehension of the simulated water cycle not only drives significant advancements in model accuracy but also plays a pivotal role in improving the precision of current-time extreme event predictions through a collaborative ensemble approach involving these refined models.
, and   refers to the precipitable flux.The horizontal moisture flux divergence can be decomposed into the following equations:

Figure 2 .
Figure 2. The annual means of ensemble models (a-d) are as follows: (a, c) represent the average of three models from Group 1, which includes IsoGSM.JRA55, MIROC5-iso.JRA55, ECHAM6-wiso.JRA55.Panels (b, d) are based on nudging to three reanalyses from Group 2: IsoGSM.JRA55, IsoGSM.NCEPR2, and IsoGSM.ERA5.The standard deviations for these ensembles (a-d) are represented as (e-h): (e) and (f) for δ 18 O (‰ or per mil) and (g, h) for d-excess in precipitation.The colored circles depict the annual mean of the respective Global Network of Isotopes in Precipitation observations.

Figure 3 .
Figure 3. Taylor diagram demonstrates model performance at Global Network of Isotopes in Precipitation (GNIP) sites for monthly mean temperature, precipitation, and isotopic ratio (δ 18 O and d-excess) in precipitation, benchmarked against observations from the respective data sets from Climatic Research Unit gridded Time Series (CRU TS) v4.06 near-surface temperature, Global Precipitation Climatology Project v2.3 monthly satellite-gauge combined precipitation, and GNIP, as references.The metrics for assessing performance skill include the averaged Pearson's correlation coefficients and the standardized deviations (represented as the ratio of simulation variance to reference observation variance).

Figure 4 .
Figure 4. (a-d) Spatial distribution of temporal Pearson's correlation between model ensemble of Group 1 and observations at Global Network of Isotopes in Precipitation (GNIP) station locations for temperature (CRU TS), precipitation (Global Precipitation Climatology Project), δ 18 O and d-excess, respectively.GNIP observation data are interpolated to T42 horizontal resolution for model comparison by averaging the values of the points closest to T42 grids.Seasonality is removed in modeled and observed isotope values in GNIP sites to alleviate more substantial seasonal variation in the tropics.Panels (e-h) as for (a-d) but with daily water vapor isotopes data set Stable Water Vapor Isotope Database.The circles are the site locations.

Figure 5 .
Figure 5.The annual mean of (a) integrated column of water vapor (colored background, kg/m 2 ) and vertically integrated water vapor transport (kg/m/s) (b) δ 18 O in the model ensemble of Group 1.Vapor amounts of each pressure thickness are weighted and integrated from 1,000 hPa to 10 hPa.q − δD pair diagram for moistening and dehydrating process over ocean and land in the collocated annual mean of (c)TES and (d, e)  simulations between vertical integrated (825 to 464 hPa) layer in which estimated HDO is primarily sensitive.Black lines connect the average of all simulations' ocean and land in different latitude areas.

Figure 6 .
Figure 6.Temporal Pearson's correlation of vertical daily δD between collocated TES and the model ensemble of Group 1 (a) in the vapor amount-weighted integrated column from 825 to 464 hPa, (b) at 464 hPa, (c) at 681 hPa, and (d) at 825 hPa.

Figure 7 .
Figure 7.Comparison of global (60°S to 60°N) daily climatology of δD from TES with model ensemble of Group 1 (based on Figure S2a in Supporting Information S1).Plots (a) to (d) are for the integrated column (825-464 hPa) and the 464, 681, and 825 hPa pressure levels, respectively.Yellow scatter points represent the ensemble means of our three models.Multipleyear daily data TES covering from 2004 to 2018 in a non-constant time frame is used to make the annual cycle, and the year 2008 was removed due to a lack of spatial coverage of TES.

Figure 8 .
Figure 8. Global time means of decomposed oxygen isotopic change rates in the model ensemble from Group 1, derived from (a) evaporation, (b) precipitation, and (c) horizontal moisture flux.The contributions of decomposed processes are shown in (d) evaporation flux-derived, (e) precipitation flux-derived, and (f) horizontal moisture flux-derived changing δ 18 O in the atmosphere.

Figure 9 .
Figure 9. Regional isotope mass budget analysis considering (a) isotopic change rate (‰/day) and (b) relative importance for evaporation, precipitation, and horizontal moisture flux in the model ensemble of Group 1.The white bar (∆δw) represents the residual after considering all isotopic processes (∆δw e + ∆δw p + ∆δw q ), indicating whether the water vapor oxygen isotope (δ 18 O) is depleting or enriching in the atmosphere over the long-term period from 1979 to 2020.The regions are shown in Figure 8a and Figure S4 in Supporting Information S1, and the values are averaged by the areas.

Figure 10 .
Figure10.The standard deviation of decomposed isotopic change rate is presented for multiple models (Group 1, represented in brown) and multiple reanalyses (Group 2, represented in lighter color).Global averages of spreads in (a) are from three models and three reanalyses for nudging.Model spread (Group 1) is 5.6 times larger than simulations conducted with different reanalyses (Group 2), and percentages represent proportions to total spread.The bars for uncertainty in (b) present the regional averages of spread for evaporation, precipitation, and horizontal moisture flux, as shown in Figures11d-11f.
comparing the effects of different models and reanalyses, depending on the specific phenomena and water phases under investigation.

Figure 11 .
Figure 11.Standard deviation map of decomposed isotopic change rate from multiple models (a) for evaporation, (b) for precipitation, and (c) for horizontal moisture flux.The regional standard deviation of decomposed isotopic change rate from three models (brown for Group 1) and three reanalyses (light color for Group 2): (d) for evaporation from shown in (a, e) for precipitation from (b), and (f) for horizontal moisture flux from (c)

Table 1
Model Descriptions, Reanalysis Data Sets, and Observations Used inThis Study Are Detailed as Follows: IsoGSM is Nudged by the Japanese 55-Year Reanalysis (JRA55), the National Centers for Environmental Prediction Reanalysis Version 2 (NCEP-R2), and the Fifth Generation European Center for Medium-Range Weather Forecasts Atmospheric Reanalysis (ERA5) The National Centers for Environmental Prediction-Department of Energy (NCEP-DOE) Atmospheric Model Intercomparison Project (AMIP-II) reanalysis II (R2) was conducted by the mission of the National Energy Research Scientific Computing Center (NERSC).It is an improved version of NCEP-NCAR Reanalysis I ( , where  cos  = 1 at the equator) components of the moisture flux vector.
,  ⃖⃖ ⃗  = (, ) represents the zonal (longitude,   ) and meridional (latitude,  with respect to longitude and latitude, respectively.  is the equatorial radius.Changes in precipitable water in the atmosphere (   ) are determined by the inflow and outflow of moisture flux over a specific time interval (  ∆ = 6 hr ,  6 × 60 × 60 s) :

Table 2 The
Global Average of Pearson's Correlation Coefficient Is Calculated From Model Simulations, Which Include Ensembles for Models (Mo ENS, Group 1) and Different Reanalyses (Re ENS, Group 2) 18O  = 27.1% ,    18 O  = 22.9% , and    18 O  = 50.0% .Despite the fact that simulating d-excess poses one of the most formidable challenges in water isotope modeling, our new approach holds the promise of enhancing our comprehension of kinetic fractionation.