Abstract Nonlinear wave interactions describe the resonant energy transfer between wave components, playing a fundamental role in the evolution of ocean wave spectra. Nonlinear wave interactions significantly influence wave growth and development, making them essential for accurate wave modeling. However, resolving the full six‐dimensional Boltzmann integral of the exact nonlinear wave interactions (Webb‐Resio‐Tracy method, WRT) is computationally expensive, limiting its application in real‐time operational wave forecasting and for research purposes. Current approximations, such as the Discrete Interaction Approximation (DIA), prioritize computational speed over accuracy, resulting in significant errors in wave mean parameters. Here, we introduce NLML, a machine learning (ML) emulator designed to approximate the exact nonlinear wave interactions within WAVEWATCH III (WW3), with the goal of achieving the accuracy of WRT while maintaining the stability and computational speed of DIA. By leveraging GPU capabilities such as half precision inference, we achieved substantial speedups, up to 136× faster than the WRT and only a modest 1.04× slowdown relative to DIA, while achieving 2× the accuracy of DIA in global wave spectral energy and mean wave parameters, with up to 7× higher accuracy in some regions. Unlike previous ML approaches, NLML maintained inherent stability throughout model integration in a standalone, year‐long WW3 simulation, without requiring additional constraints. Our new ML parameterization bridges the gap between accuracy and efficiency, offering a promising alternative for improving wave modeling in operational settings and research purposes.

Abstract Many Earth System Models (ESMs) approximate surface emissivity as constant over all wavelengths (broadband). This approximation misrepresents the spectral distribution of longwave (LW) flux and introduces error in the LW atmospheric heating rate. Many ESMs employ the “Effective Greybody” (EG) method, which creates bias in the effective temperature used by the atmosphere component. We implement a novel “Full Greybody” (FG) method, modifying the treatment of LW flux in the atmosphere to match that of the ocean surface component in the Energy Exascale Earth System model. We find that, relative to greybody emissivity, the ocean blackbody approximation over the ocean overestimates instantaneous mean upwelling LW surface flux by 1.38±0.05 W/m2. The common EG method introduces a mean bias of −0.24 K in the effective temperature. The FG method reduces this bias to 0.0039 K, comparable to the default blackbody bias (0.0070 K) while making ocean emission in ESMs more physically realistic.

The Energy Exascale Earth System Model version 3 (E3SMv3) represents the latest advancement in Earth system modeling developed by the U.S. Department of Energy (DOE). Building upon previous versions, E3SMv3 introduces significant updates across its coupled components to enhance capability and improve fidelity. The atmosphere component incorporates advancements in chemistry, aerosol-cloud interactions, convection, and microphysics. The ocean features a new time-stepping scheme and a higher-resolution unstructured mesh with sub-ice-shelf cavities, while the sea ice model integrates advanced snow and ice physics for more realistic cryospheric simulations. The land model introduces prognostic vegetation dynamics and a new sub-grid topographic treatment of solar radiation. A new tri-grid configuration harmonizes the horizontal grids of the land and river components for improved process coupling. It is enabled by a new non-linear remapping between the atmosphere and land. E3SMv3 underwent extensive testing through a comprehensive simulation campaign, including pre-industrial control, idealized experiments, and historical simulations spanning 1850-2024. The model demonstrates significant improvements in simulating the evolution of the historical surface temperature, particularly addressing the "pothole cooling" bias in earlier versions. Reduced aerosol-related forcing contributes to more realistic radiative forcing and better alignment with the observational record. Ocean heat content (OHC) and sea ice trends are also improved as a result.

The solar radiative processes that contribute to Earth's surface and atmospheric energy budgets are strongly dependent on wavelength. For example, snow and water vapor become more absorptive as incident radiation shifts from visible (VIS) to near-infrared (NIR) wavelengths. Some earth system models (ESMs) aggregate solar radiation into just two bands (VIS and NIR) partitioned at 0.7 m for transmission between the atmospheric and surface components. In the widely used radiative transfer model RRTMG_SW this partition is near the center of the overlap spectral band spanning 0.625-0.778 m, whose flux is often approximated as being evenly divided between the VIS and NIR surface bands. Using a hyperspectral radiative transfer model, we show that the fractional downwelling surface flux within the overlap band is usually distributed about 55.5%:44.5% VIS:NIR. This improved approximation shifts as much as 4.27 W m-2 from the NIR to the VIS band, leading to an instantaneous decrease in surface absorption of up to 0.91 W m-2 over snow-covered surfaces. Century-long fully coupled ESM simulations show surface absorption over snow decreases by over 1.6 W m-2. The coupled response in sea ice regions amplifies the initial forcing ten-fold, and increases seasonal sea ice area by up to 4.9%. These results highlight the importance of accurately representing the spectral distribution of solar radiation in the cryosphere.

This study investigates the sensitivity of the Madden-Julian oscillation (MJO) to changes to the bulk flux parameterization and the role of ocean surface waves in air-sea coupling using a fully coupled ocean-atmosphere-wave model. The atmospheric and ocean model components of the Energy Exascale Earth System Model (E3SM) are coupled to a spectral wave model, WAVEWATCH III (WW3). Two experiments with wind speed-dependent bulk algorithms (NCAR and COARE3.0a) and one experiment with wave-state-dependent flux (COR3.0a-WAV) were conducted. We modify COARE3.0a to include surface roughness calculated within WW3 and also account for the buffering effect of waves on the relative difference between air-side and ocean-side momentum flux. Differences in surface fluxes, primarily caused by discrepancies in drag coefficients, result in significant differences in MJO's properties. While COARE3.0a has better convection-circulation coupling than NCAR, it exhibits anomalous MJO convection east of the date line. The wavestate-dependent flux (COR3.0-WAV) improves the MJO representation over the default COARE3.0 algorithm. Strong easterlies over the Pacific Ocean in COARE3.0a enhance the latent heat flux (LHFLX). This is responsible for the anomalous MJO propagation after the date line. In COR3.0a-WAV, waves reduce the anomalous easterlies, leading to a decrease in LHFLX and MJO dissipation after the date line. These findings highlight the role of surface fluxes in MJO simulation fidelity. Most importantly, we show that the proper treatment of wave-induced effects in bulk flux parameterization improves the simulation of coupled climate variability.

The Arctic is undergoing increased warming compared to the global mean, which has major implications for freshwater runoff into the oceans from seasonal snow and glaciers. Here, we present high-resolution (2.5 km) simulations of glacier mass balance, runoff, and snow conditions on Svalbard from 1991-2022, one of the fastest warming regions in the world. The simulations are created using the CryoGrid community model forced by Copernicus Arctic Regional ReAnalysis (CARRA) (1991-2021) and AROME-ARCTIC forecasts (2016-2022). Updates to the water percolation and runoff schemes are implemented in the CryoGrid model for the simulations. In situ observations available for Svalbard, including automatic weather station data, stake measurements, and discharge observations, are used to carefully evaluate the quality of the simulations and model forcing.We find a slightly negative climatic mass balance (CMB) over the simulation period of -0.08 mw.e.yr-1 but with no statistically significant trend. The most negative annual CMB is found for Nordenskioldland (-0.73 mw.e.yr-1), with a significant negative trend of -0.27 mw.e. per decade for the region. Although there is no trend in the annual CMB, we do find a significant increasing trend in the runoff from glaciers of 0.14 mw.e. per decade. The average runoff was found to be 0.8 mw.e.yr-1. We also find a significant negative trend in the refreezing of -0.13 mw.e. per decade.Using AROME-ARCTIC forcing, we find that 2021/22 has the most negative CMB and highest runoff over the 1991-2022 simulation period investigated in this study. We find the simulated climatic mass balance and runoff using CARRA and AROME-ARCTIC forcing are similar and differ by only 0.1 mw.e.yr-1 in climatic mass balance and by 0.2 mw.e.yr-1 in glacier runoff when averaged over all of Svalbard. There is, however, a clear difference over Nordenskioldland, where AROME-ARCTIC simulates significantly higher mass balance and significantly lower runoff. This indicates that AROME-ARCTIC may provide similar high-quality predictions of the total mass balance of Svalbard as CARRA, but regional uncertainties should be taken into consideration.The simulations produced for this study are made publicly available at a daily and monthly resolution, and these high-resolution simulations may be re-used in a wide range of applications including studies on glacial runoff, ocean currents, and ecosystems.

Surface gravity waves play a critical role in several processes, including mixing, coastal inundation, and surface fluxes. Despite the growing literature on the importance of ocean surface waves, wind–wave processes have traditionally been excluded from Earth system models (ESMs) due to the high computational costs of running spectral wave models. The development of the Next Generation Ocean Model for the DOE’s (Department of Energy) E3SM (Energy Exascale Earth System Model) Project partly focuses on the inclusion of a wave model, WAVEWATCH III (WW3), into E3SM. WW3, which was originally developed for operational wave forecasting, needs to be computationally less expensive before it can be integrated into ESMs. To accomplish this, we take advantage of heterogeneous architectures at DOE leadership computing facilities and the increasing computing power of general-purpose graphics processing units (GPUs). This paper identifies the wave action source terms, W3SRCEMD, as the most computationally intensive module in WW3 and then accelerates them via GPU. Our experiments on two computing platforms, Kodiak (P100 GPU and Intel(R) Xeon(R) central processing unit, CPU, E5-2695 v4) and Summit (V100 GPU and IBM POWER9 CPU) show respective average speedups of 2× and 4× when mapping one Message Passing Interface (MPI) per GPU. An average speedup of 1.4× was achieved using all 42 CPU cores and 6 GPUs on a Summit node (with 7 MPI ranks per GPU). However, the GPU speedup over the 42 CPU cores remains relatively unchanged (∼ 1.3×) even when using 4 MPI ranks per GPU (24 ranks in total) and 3 MPI ranks per GPU (18 ranks in total). This corresponds to a 35 %–40 % decrease in both simulation time and usage of resources. Due to too many local scalars and arrays in the W3SRCEMD subroutine and the huge WW3 memory requirement, GPU performance is currently limited by the data transfer bandwidth between the CPU and the GPU. Ideally, OpenACC routine directives could be used to further improve performance. However, W3SRCEMD would require significant code refactoring to make this possible. We also discuss how the trade-off between the occupancy, register, and latency affects the GPU performance of WW3.

In-situ observations of subsurface ocean temperatures are, in many regions, inconsistently distributed in time and space. These spatio-temporal inconsistencies in the observational network lead to difficulties in utilizing those observations effectively for ocean model evaluation or understanding larger-scale ocean characteristics. Model accuracy of subsurface ocean characteristics is especially important within regions that contain complex ocean structures. One such region is the European Arctic which not only contains several types of water masses with unique characteristics, but also wintertime sea ice coverage and complex bathymetry. This study presents an unsupervised neural networking technique that can be used in combination with traditional ocean model evaluation techniques to provide additional information on the accuracy of modeled vertical ocean temperature profiles. Self-organizing maps is an unsupervised machine learning technique that we apply to approximately twenty thousand Argo and CTD temperature profiles from 2012 to 2020 in the European Arctic to categorize the observed vertical ocean temperature structures in the top 150 m. The observed ocean profile categories, or neurons, defined by the self-organizing map show strong spatial and temporal dependencies. We then use the neuron weights, or the learned temperature profile structure of each neuron, to validate the spatial and temporal variability of modeled vertical temperature structures. This analysis gives us new insights about the model's capabilities to reproduce specific vertical structures of the top-most ocean layer within different regions and seasons. Mapping modeled ocean temperature profiles onto the neuron-space of the observationally-defined self organized map highlights the potential of this method to advance our understanding of model deficiencies in that region.

Accurately simulating the interactions between the components of a coupled Earth modelling system (atmosphere, sea-ice, and wave) on a kilometer-scale resolution is a new challenge in operational numerical weather prediction. It is difficult due to the complexity of interactive mechanisms, the limited accuracy of model components and scarcity of observations available for assessing relevant coupled processes. This study presents a newly developed convective-scale atmosphere-wave coupled forecasting system for the European Arctic. The HARMONIE-AROME configuration of the ALADIN-HIRLAM numerical weather prediction system is coupled to the spectral wave model WAVEWATCH III using the OASIS3 model coupling toolkit. We analyze the impact of representing the kilometer-scale atmosphere-wave interactions through coupled and uncoupled forecasts on a model domain with 2.5 km spatial resolution. In order to assess the coupled model’s accuracy and uncertainties we compare 48-hour model forecasts against satellite observational products such as Advanced Scatterometer 10 m wind speed, and altimeter based significant wave height. The fully coupled atmosphere-wave model results closely match both satellite-based wind speed and significant wave height observations as well as surface pressure and wind speed measurements from selected coastal station observation sites. Furthermore, the coupled model contains smaller standard deviation of errors in both 10m wind speed and significant wave height parameters when compared to the uncoupled model forecasts. Atmosphere and wave coupling reduces the short term forecast error variability of 10 m wind speed and significant wave height with the greatest benefit occurring for high wind and wave conditions.

Increased economic activity and research interest in the Arctic raise the need for better wave forecasts in the marginal ice zone (MIZ). Mathematical and numerical models of wave propagation in sea ice would benefit from more in situ data for validation. This study presents shipborne wave measurements from the MIZ where altimeter readings are corrected for ship motion to obtain estimated single point ocean surface elevation. From the combined measurements, we obtain significant wave height and zero up-crossing period, as well as one-dimensional wave spectra. In addition, we provide spectra and integrated parameters obtained from inertial motion units (IMU) placed on ice floes inside the MIZ. The results are compared with integrated parameters from the WAM-4 spectral wave model over a period of three days in the open ocean. We also compare our measurements outside and inside the MIZ with hindcast data from the new pan-Arctic WAM-3 model and the Wave Watch III model for the European Arctic, which both model wave attenuation in sea ice. A good agreement is found with WAM-4 and WW3 in zero up-crossing period and significant wave height outside the MIZ, where deviations are less than 23 %. WAM-3 is on the other hand up to 60 % higher than observations. WW3 and WAM-3 are able to estimate the trends for significant wave height and zero up-crossing period inside the MIZ, although the discrepancies with respect to the observations were larger than in the open ocean. Wave damping by sea ice is investigated by looking at the spatial attenuation coefficients. Predicted attenuation coefficients are found to be 72-78 % smaller for WW3 and 3-131 % larger for WAM-3 compared to observations. Hence, further model tuning is necessary to better estimate wave parameters in the ice.

Numerous oceanic and atmospheric phenomena influence El Niño–Southern Oscillation (ENSO) variability, complicating both prediction and analysis of the mechanisms responsible for generating ENSO diversity. Predictability of ENSO events depends on the characteristics of both the forecast initial conditions and the stochastic forcing that occurs subsequent to forecast initialization. Within a linear inverse model framework, stochastic forcing reduces ENSO predictability when it excites unpredictable growth or interference after the forecast is initialized, but also enhances ENSO predictability when it excites optimal initial conditions that maximize deterministic ENSO growth. Linear inverse modeling (LIM) allows for straightforward separation between predictable signal and unpredictable noise and so can diagnose its own skill. While previous LIM studies of ENSO focused on deterministic dynamics, here we explore how noise forcing influences ENSO diversity and predictability. This study identifies stochastic forcing details potentially contributing to the development of central Pacific (CP) or eastern Pacific (EP) ENSO characteristics. The technique is then used to diagnose the relative roles of initial conditions and noise forcing throughout the evolution of several ENSO events. LIM results show varying roles of noise forcing for any given event, highlighting its utility in separating deterministic from noise-forced contributions to the evolution of individual ENSO events. For example, the strong 1982 event was considerably more influenced by noise forcing late in its evolution than the strong 1997 event, which was more predictable with long lead times due to its deterministic growth. Furthermore, the 2014 deterministic trajectory suggests that a strong event in 2014 was unlikely.

Interactions between the Pacific meridional mode (PMM) and El Nino-Southern Oscillation (ENSO) are investigated using the National Center for Atmospheric Research (NCAR) Community Earth System Model (CESM) and an intermediate coupled model (ICM). The two models are configured so that the CESM simulates the PMM but not ENSO, and the ICM simulates ENSO but not the PMM, allowing for a clean separation between the PMM evolution and the subsequent ENSO response. An ensemble of CESM simulations is run with an imposed surface heat flux associated with the North Pacific Oscillation (NPO) generating a sea surface temperature (SST) and wind response representative of the PMM. The PMM wind is then applied as a forcing to the ICM to simulate the ENSO response. The positive (negative) ensemble-mean PMM wind forcing results in a warm (cold) ENSO event although the responses are not symmetric (warm ENSO events are larger in amplitude than cold ENSO events), and large variability between ensemble members suggests that any individual ENSO event is strongly influenced by natural variability contained within the CESM simulations. Sensitivity experiments show that 1) direct forcing of Kelvin waves by PMM winds dominates the ENSO response, 2) seasonality of PMM forcing and ENSO growth rates influences the resulting ENSO amplitude, 3) ocean dynamics within the ICM dominate the ENSO asymmetry, and 4) the nonlinear relationship between PMM wind anomalies and surface wind stress may enhance the La Nina response to negative PMM variations. Implications for ENSO variability are discussed.