<figure><img src="image_1.png" /></figure>

HAL
open science

# A novel sea surface evaporation scheme assessed by the thermal rotating shallow water model

Masoud Rostami, Stefan Petri, Bijan Fallah, Farahnaz Fazel-Rastgar

► To cite this version:

Masoud Rostami, Stefan Petri, Bijan Fallah, Farahnaz Fazel-Rastgar. A novel sea surface evaporation scheme assessed by the thermal rotating shallow water model. Atmospheric Science Letters, 2025, 26 (1), 10.1002/asl.1287. hal-04925511

HAL Id: hal-04925511

https://enpc.hal.science/hal-04925511v1

Submitted on 2 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

Atmospheric Science Letters

<figure><img src="image_3.png" /></figure>

# A novel sea surface evaporation scheme assessed by the thermal rotating shallow water model

Masoud Rostami<fnref n="1" />,<fnref n="2" /> ID

Stefan Petri<fnref n="1" /> ID

Bijan Fallah<fnref n="3" />

Farahnaz Fazel-Rastgar<fnref n="4" />

**Correspondence**

Masoud Rostami, Potsdam Institute for Climate Impact Research (PIK), Member of the Leibniz Association, P.O. Box 60 12 03, D-14412 Potsdam, Germany.

Email: rostami@pik-potsdam.de; masoud.rostami@lmd.ipsl.fr

**Funding information**

Virgin Unite USA, Inc. through the Planetary Boundary Science Lab project; European Regional Development Fund (ERDF); German Federal Ministry of Education and Research

**Abstract**

In this study, a novel sea surface evaporation scheme, along with its corresponding bulk aerodynamic formulation, is proposed to estimate sea surface evaporation, columnar humidity, and precipitation distribution within the atmosphere. The scheme is based on three distinct functions, each dependent on a single variable: zonal wind velocity, tropospheric (potential) temperature, and free convection. It is shown that the normalized Clausius–Clapeyron formula requires an adjustable scaling factor for real-world applications, calibrated using empirical fitness curves. To validate the proposed approach, we employ a model based on the pseudo-spectral moist-convective thermal rotating shallow water model, with minimal parameterization over the entire sphere. ECMWF Reanalysis 5th Generation (ERA5) reanalysis data are used to compare the model's results with observations. The model is tested across different seasons to assess its reliability under various weather conditions. The Dedalus algorithm, which handles spin-weighted spherical harmonics, is employed to address the pseudo-spectral problem-solving tasks of the model.

**KEYWORDS**

Aeolus 2.0, bulk aerodynamic scheme, moist convection, sea surface evaporation, thermal rotating shallow water (TRSW) model

# 1 | INTRODUCTION

The bulk aerodynamic method plays a central role in modeling oceanic and atmospheric processes, particularly by enabling accurate estimations of energy fluxes into the atmosphere. The performance of General Circulation Models (GCMs), both stand-alone and coupled, is highly sensitive to the effectiveness of this parameterization. For example, in the Nucleus for European Modeling of the Ocean version 4 (NEMOv4), variations in bulk parameterizations for turbulent air-sea flux calculations can

significantly influence modeled sea surface temperatures (SSTs) (Bonino et al., 2022).

Direct measurement of surface fluxes is often limited by logistical constraints, including instrumentation and data collection challenges. Hence, parameterization remains essential for accurate modeling, even with the availability of extensive datasets. Historically, the foundation for modern parameterizations began with John Dalton's 1798 work on evaporation, followed by the formulation of bulk aerodynamic coefficients for heat and momentum by Morton et al. (1956). This

<footnote n="1">¹Potsdam Institute for Climate Impact Research (PIK), Member of the Leibniz Association, Potsdam, Germany</footnote>

<footnote n="2">²Laboratoire de Météorologie Dynamique (LMD), Sorbonne Université (SU), Ecole Normale Supérieure (ENS), Paris, France</footnote>

<footnote n="3">³Deutsches Klimarechenzentrum GmbH (DKRZ), Hamburg, Germany</footnote>

<footnote n="4">⁴School of Chemistry and Physics, University of KwaZulu Natal, Durban, South Africa</footnote>

<footnote>This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>

<footnote>© 2025 The Author(s). *Atmospheric Science Letters* published by John Wiley & Sons Ltd on behalf of Royal Meteorological Society.</footnote>
foundational work has since been refined by researchers like Roll (1965), Pond et al. (1974), and Fairall et al. (1996), forming the basis of current parameterizations for air-sea fluxes. These parameterizations estimate evaporation based on gradients of specific humidity, temperature, wind speed, and air-sea temperature differences, and are widely implemented in numerical weather and climate models to simulate Earth-system processes. Sensible heat flux ($H_s$) and moisture flux ($E$) are commonly modeled by the following equation using the bulk aerodynamic method:

$$H_s = (\rho c_p) C_T (U - U_s) \Delta\theta, \quad (1a)$$

$$E = C_q (U - U_s) \Delta q. \quad (1b)$$

Here, $H_s$ is the sensible heat flux in watts per square meter, $\rho$ is air density, $c_p$ is the specific heat of air at constant pressure, and $C_T$ and $C_q$ are the bulk aerodynamic coefficients for heat and moisture, respectively. The coefficients $C_T$ and $C_q$ depend on atmospheric stability and are determined using empirical methods or meteorological data, often through Monin-Obukhov Similarity Theory (Monin & Obukhov, 1954; Smith, 1988). In these equations, $U$ represents mean wind speed, $U_s$ is the ocean surface velocity, $\Delta\theta$ is the temperature difference between the sea surface and the air, $E$ is the moisture flux, and $\Delta q$ is the difference in humidity between the surface and the air. Although this approach has been refined to account for atmospheric dynamics, limitations persist, especially in capturing complex conditions like varying global sea surface evaporation patterns that affect humidity differences ($\Delta q$) in the bulk aerodynamic scheme. Modeling $\Delta q$ accurately, influenced by temperature and wind velocity, is crucial for realistic climate simulations.

While earlier models used empirical climatology-based parameterizations, advanced drag formulations now dominate, often with sophisticated schemes that better capture heat and moisture exchanges. Even small inaccuracies in parameterization, as highlighted by Equations (1a and 1b), can skew climate predictions and impact cloud-resolving simulations (Wang et al., 1996). In particular, the mass transfer coefficient, representing drag forces, varies across terrestrial and marine environments, influenced by factors like atmospheric stability and surface roughness, as outlined in the Penman-Monteith equation (Monteith, 1981). For non-uniform inland water surfaces, which are rarely homogeneous, rigorous assumptions are required to model evaporation and heat transfer accurately (Brutsaert & Yu, 1968). Here, $C_q$ reflects the unique transfer characteristics of specific

water bodies, influenced by their geometry, surrounding terrain, and climate, making it an essential component in refined surface flux parameterization.

In this study, deviating from conventional methodologies (as represented by Equation 1b), we propose a novel scheme for quantifying the bulk volume of atmospheric water vapor, attributable to both thermal heating and wind-induced forces. In this scheme, the bulk aerodynamic flux can be estimated as the sum of two nonlinear functions, taking into account the large-scale spatiotemporal variation of potential temperature and the velocity field of the atmosphere in the lower troposphere. The dependence of sea surface evaporation on temperature in this scheme aligns with a modified version of the normalized Clausius–Clapeyron equation. For detailed information, see Section 2.

Bridging the gap between idealized models and complex GCMs is essential for advancing climate theory. Intermediate-complexity models provide a critical framework for exploring causal mechanisms while maintaining sufficient realism. As Held (2005) emphasized: “The health of climate theory/modeling in the coming decades is threatened by a growing gap between high-end simulations and idealized theoretical work. In order to fill this gap, research with a hierarchy of models is needed.” This study employs a two-layer moist-convective thermal rotating shallow water (mcTRSW) model with an advanced evaporation scheme to address this need, offering a balance between theoretical exploration and physical applicability.

## 2 | THE NOVEL BULK AERODYNAMIC SCHEME

The change in Gibbs free energy ($dG$) for a system undergoing changes in temperature ($dT$), pressure ($dP$), and changes in chemical potential (associated with phase change) with respect to the number of moles of each component ($dN_i$) can be written as:

$$dG = -SdT + VdP + \sum_i \mu_i dN_i. \quad (2)$$

Gibbs free energy represents the maximum amount of reversible work that can be extracted from a system under constant temperature and pressure conditions. The first term on the right-hand side of Equation (1b), $-SdT$, accounts for changes in entropy due to temperature variations, and the term $VdP$ accounts for work done against external pressure. Together, they represent the contribution from changes in entropy and volume, respectively. $\sum_i \mu_i dN_i$ represents the contribution from changes in
<figure><img src="image_4.png" /></figure>

chemical potential and moles of each component. At equilibrium, the chemical potential of each component is equal between phases. This leads to a simplified expression for the change in Gibbs free energy, where only the terms related to changes in temperature (dT) and pressure (dP) are considered:

$$dG = -SdT + VdP. \qquad (3)$$

This formulation indicates that at equilibrium, the system achieves a state of maximal stability and minimal potential to perform work on its surroundings. The term “maximal stability” denotes a state of thermodynamic equilibrium in which the system achieves its minimum energy configuration, characterized by the absence or reduction of gradients in energy distribution. Similarly, “minimal potential” refers to the Gibbs free energy attaining its lowest possible value at equilibrium, signifying a condition where the system is thermodynamically stable and incapable of performing further work under the given constraints. To incorporate wind velocity into the Gibbs equation, it is important to note that in real ocean-atmosphere interaction, surface evaporation, particularly at mid-latitudes, does not occur spontaneously without wind forcing. We denote the vaporization work driven by wind energy as dW<sub>wind</sub>. This work overcomes cohesive forces at the sea surface, enabling water molecules to detach and enter the air. Additionally, wind speed influences turbulent mixing; increased wind speeds enhance turbulence and promote greater upward diffusion. Thus, the modified expression for seawater in equilibrium (dG = 0) with its atmosphere, excluding the effects of ice and clouds and accounting for wind-induced work, is given by:

$$\begin{aligned} -S_{\text{ocean}}dT + V_{\text{ocean}}dP_{\text{vap}} + dW_{\text{wind}} \\ = -S_{\text{atmosphere}}dT + V_{\text{atmosphere}}dP_{\text{vap}}. \end{aligned} \qquad (4)$$

In the case of equilibrium, where the chemical potential of water equals that of air, the equation is further simplified. This means that the chemical potential terms (μ<sub>ocean</sub> and μ<sub>atmosphere</sub>) are equal, leading to a cancelation of those terms in the equation. The simplified form accounts only for changes in temperature (dT), pressure (dP<sub>vap</sub>), and the work done by the wind (dW<sub>wind</sub>). For the system of water-atmosphere to be in equilibrium we have:

$$\frac{dP_{\text{vap}}}{dT} = \frac{\Delta S}{\Delta V} + \frac{dW_{\text{wind}}}{\Delta VdT}. \qquad (5)$$

This expression illustrates the ratio of the change in pressure (dP<sub>vap</sub>) to the change in temperature (dT)

within the system. It is articulated in terms of alterations in volume (ΔV = V<sub>atmosphere</sub> - V<sub>ocean</sub>) and entropy (ΔS = S<sub>atmosphere</sub> - S<sub>ocean</sub>). Additionally, it incorporates the influence of work done by the wind on the water, which can play a significant role in processes such as evaporation. In this context, the symbol Δ denotes the change occurring between phases. The initial term on the right-hand side of Equation (5) bears resemblance to the Clausius–Clapeyron equation. In the Clausius–Clapeyron equation, dP<sub>vap</sub>/dT = ΔH<sub>vap</sub>/(TΔV), where ΔH<sub>vap</sub> signifies the enthalpy change during vaporization. Additionally, ΔV can be presented as a function of pressure and temperature: R<sub>v</sub>T/P, where R<sub>v</sub> represents the specific gas constant for water vapor. The dependence of the Clausius–Clapeyron equation on T is given by:

$$\frac{d(\ln P_{\text{vap}})}{dT} = \frac{\Delta H_{\text{vap}}}{R_{\text{v}}T^2}. \qquad (6)$$

The application of the Clausius–Clapeyron equation in meteorology and sea surface evaporation encounters challenges stemming from idealized assumptions and simplifications. While the equation establishes a fundamental link between temperature and saturation vapor pressure, real-world complexities such as non-uniform environmental conditions, dynamic atmospheric changes (including variations in humidity, pressure, and temperature over time), feedback mechanisms, and local influences pose significant limitations. The equation's static nature and neglect of factors like phase changes, variable latent heat, and boundary layer dynamics further restrict its accuracy in capturing the intricacies of evaporation (Iribarne & Godson, 1973). Additionally, the omission of critical parameters, such as wind speed and sea surface properties, highlights the need for new schemes to account for the multifaceted nature of atmospheric and oceanic processes, particularly in the context of evolving climate conditions. To refine the sea surface evaporation scheme and the associated saturation vapor pressure, we introduce an unknown adjustable scaling factor, denoted as α, into the Clausius–Clapeyron equation (Equation 7a). This empirical factor quantifies deviations from the idealized Clausius–Clapeyron scaling, accounting for complexities inherent in real-world atmospheric and oceanic systems. These deviations arise from non-idealities such as boundary layer turbulence, vertical mixing, variability in latent heat release, and other thermodynamic constraints. The value of α, acting as the exponent of T in the Clausius–Clapeyron equation, is then determined by optimizing the empirical fit between model results and observational data. In our numerical experiments, α < 1 (≈ 0.6 – 0.7).
In our proposed methodology, the source of humidity is a function of three variables, each representing a distinct regime: vaporization due to the heating force, as a function of the lower layer temperature, $T_1$ (Equation 7a); vaporization due to the wind force, as a function of the normalized zonal wind velocity, $u_1^n$ (Equation 7b); and free convection (Equation 7c). The latter becomes predominant in scenarios where the former forces are negligible or absent. The humidity source can thus be expressed as follows:

$$P_{\text{vap}}(T_1) = P_0 \exp\left[\frac{-\Delta H_{\text{vap}}}{R_v}\left(\frac{1}{T_1^{\alpha}} - \frac{1}{T_0^{\alpha}}\right)\right], \quad (7a)$$

$$P_{\text{vap}}(|u_1^n|) = \exp\left(\frac{|u_1^n|^{\alpha_v}}{\alpha_{2v}}\right), \quad (7b)$$

$$\text{Fr}\text{.Conv.} = A_F(Q^s - q_1), \quad \text{if } (|u_1^n| \ll 1 \text{ and } b_1 \ll b_1^{\text{max}}), \quad (7c)$$

$$S_{\text{vapor}}(T_1, |u_1^n|, q_1) = \mathbb{H}(Q^s - q_1) \left( A_T \hat{P}_{\text{vap}}(T_1) \right. \quad (7d) \\ \left. + A_u \hat{P}_{\text{vap}}(|u_1^n|) + \text{Fr}\text{.Conv.} \right),$$

where $S_{\text{vapor}}$ denotes the humidity source, contingent upon the normalized value of $P_{\text{vap}}$ in relation to its maximum, while $A_T$ is an adjustable coefficient. $P_0$ is a reference vapor pressure, $\Delta H_{\text{vap}}$ symbolizes the enthalpy change during vaporization, which may be temperature-dependent. $T_0$ stands as a reference temperature. $\hat{P}_{\text{vap}}(|u_1^n|) = P_{\text{vap}}(|u_1^n|)/\text{Max}[P_{\text{vap}}(|u_1^n|)]$ represents the normalized value of water vapor pressure due to wind-induced work of zonal velocity in the lower layer, scaled

relative to its maximum value. The empirical function in Equation (7b) also includes adjustable scaling factors, $\alpha_{1v}$ and $\alpha_{2v}$, which require calibration in relation to real-world observations and the corresponding lower layer thickness in the mcTRSW model. Generally, the rate of increase in evaporation with respect to wind velocity is less pronounced compared to the rate of increase with respect to temperature (Figure 1). In this study, the values $\alpha_{1v} = 1.2$ and $\alpha_{2v} = 0.7$ have been employed a trend consistent with findings from various laboratory Wirangga et al. (2023) and field observations Penman (1948). The adjustable scaling factor $\alpha$ introduces a modification that enables the application of the equation in the context of ocean-atmosphere interactions. Figure 1 illustrates how different scaling factors cause deviations in vapor pressure concerning temperature. The term Fr.Conv., where $A_F$ is an adjustable coefficient, denotes the phenomenon of free convection occurring over the water body, independent of zonal wind and buoyancy anomaly in the lower layer, in cases where both of them are negligible with respect to the leading order magnitude.

Equation (7d) utilizes a normalized approach, ensuring vapor pressure independence from specific parameter values, like $P_0$ and $P_{\text{vap}}$. This emphasizes trends over absolute values, enhancing generalizability across diverse conditions. The normalized form facilitates consistent comparisons of vapor pressure variations concerning temperature, accounting for the scaling factor $\alpha$. The equation's versatility is augmented by the adjustable tuning parameter $A_T$, providing a means to fit the model to observational or experimental data. This approach generalizes vapor pressure behavior with respect to temperature in ocean-atmosphere interactions, mitigating sensitivity to precise values of $P_0$ and $\Delta H_{\text{vap}}$. $q_1$ denotes

<figure><img src="image_5.png" /><figcaption>FIGURE 1 Exploring the influence of various scaling factors on sea surface vapor pressure and comparing with the Clausius–Clapeyron equation.</figcaption></figure>
the actual moisture content, and $Q_s$ signifies the saturation moisture value. The Heaviside step function $H(Q_s - q_1)$ ensures that vaporization occurs only when the actual moisture content is lower than the saturation threshold. Employing a rapid escalation function, $\hat{P}_{vap}(T_1)$, enables the capture of the nonlinear reaction of sea surface evaporation to (potential-)temperature variations. As a result, the zonally-averaged distribution of total precipitation demonstrates pronounced sharpness at

latitudes adjacent to areas marked by maximum buoyancy anomalies or potential temperature, such as the equatorial zone.

The function $P_{vap}(|u_1^n|)$ acts as a source of water vapor humidity, albeit with a reduced coefficient, thereby extending its influence over land areas. The decision to solely incorporate zonal velocity in this parameterization is rooted in the pronounced impact of westerlies, which facilitate moisture transport from western water bodies to

<figure><img src="image_6.png" /><figcaption>FIGURE 2 Zonally-averaged distributions of potential temperature at various pressure levels are presented for the monthly averages of March, June, September, and December for the year 1980. Thin dashed curves represent the ERA5 data. Thick black solid and dashed lines depict vertically integrated values for the lower and upper layers, respectively, which are used as input to the model. The upper layer includes levels at 450 and 250 mbar, while the lower layer encompasses levels at 650, 800, 900, and 1000 mbar.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 3 Zonal velocity averages at different pressure levels are presented for the monthly averages of March, June, September, and December for the year 1980. Thin dashed curves represent the ERA5 data. Thick black solid and dashed lines depict vertically integrated values for the lower and upper layers, respectively, which are utilized as input to the model. The upper layer includes levels at 450 and 250 mbar, while the lower layer encompasses levels at 650, 800, 900, and 1000 mbar.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIGURE 4 The Kolmogorov–Smirnov (KS) statistic compares the empirical cumulative distribution functions of the monthly averages for March, June, September, and December of the year 1980 between the Aeolus 2.0 model and ERA5 data. Total precipitation is presented in (m month⁻¹).</figcaption></figure>

eastern regions across land masses. This rationale aligns with the prevailing philosophy underlying westerlies-dominated climate regimes in mid-latitudes (Huang et al., 2015). Hence, it is unsurprising that regions experiencing a damping of zonal velocity also exhibit a notable reduction in total precipitation. This relationship is further elucidated through a comparison of Figure 3 with Figure 5. Additionally, data utilized for evaluating the bulk aerodynamic scheme (see Section 3) indicates that only a minute fraction of cloud liquid water content (≈7%) is observed at altitudes exceeding 450 mbar. In this scheme, precipitation predominantly occurs in con-vertent zones rather than directly at the source of humid-ity. Observational studies provide further support to these dynamics, illustrating how moisture flows can indeed augment precipitation, as evidenced by their impact on the Tibetan Plateau (Jiang et al., 2023), and Arid Central Asia (ACA), where anomalous westerly winds contribute

to increased water vapor influx (Huang et al., 2013, 2015; Wei et al., 2017).

3 | ASSESSMENT AND VALIDATION OF THE SCHEME

We employ Aeolus 2.0, a pseudo-spectral mcTRSW model (Rostami, Petri, et al., 2024; Rostami, Severino, 2024, Rostami et al., 2022), to evaluate the bulk aerodynamic scheme (Section 2). The dynamical core of the two-layer version of Aeolus 2.0 is summarized in Appendix A. Simulations were conducted for key months—March, June, September, and December—corresponding to the solstices and equinoxes. The model is initialized using ERA5 data from the year 1980. ERA5 is a comprehensive atmospheric reanalysis dataset produced by the European Centre for Medium-Range Weather Forecasts
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 5 Comparison of zonally-averaged total precipitation for monthly averages in March, June, September, and December of the year 1980 between the Aeolus 2.0 model and ERA5 data. Total precipitation is presented in (m month⁻¹). Gray areas depict two distinctive regimes of mid-latitudinal and equatorial zones, which are under wind-induced and heat-induced regimes, providing sources of water vapor in the proposed bulk aerodynamic scheme.</figcaption></figure>

(ECMWF), providing high-quality global atmospheric information covering a wide range of variables (Hersbach et al., 2020). Reanalysis datasets like ERA5 are generated by combining historical observations with numerical weather prediction models to create consistent records of past weather and climate conditions. Figures 2 and 3 present the zonally-averaged distribution of initial potential temperature fields and velocity fields used as the initial condition for the corresponding months. Notably, while the model inherently provides the spatial distribution of water vapor at the initial stage, it undergoes a temporal evolution process known as ageostrophic adjustment. This adjustment period is essential for the model to attain equilibrium and accurately represent the spatial distribution of water vapor throughout the simulation period. The lack of such initialization may potentially lead to a compromised evaluation of the impact of advection within the model. This strategic initiation ensures the

expeditious capture of angular momentum advection by meridional circulations right from the model's inception. The Kolmogorov–Smirnov (KS) statistic (Massey, 1951) (Figure 4) presents that the distributions of total precipitation of Aeolus 2.0 and ERA5 data are very similar, and any observed differences can be attributed to random sampling variability rather than a systematic difference in the underlying distributions.

Figure 5 portrays a comparative analysis of zonally-averaged total precipitation for monthly averages in March, June, September, and December, juxtaposing data from the Aeolus 2.0 model with that from the ERA5 dataset.

In each of the four distinct time spans, a prominent concentration of total precipitation is discernible in the equatorial zone, accompanied by a secondary peak in zonally-averaged precipitation near the peak of the zonal wind. These discernible patterns can be attributed to
<figure><img src="image_12.png" /></figure>

distinct atmospheric regimes delineated by the equation for water vapor sources (Equation 7d), yielding two prominent zonally-distributed bands of elevated water vapor pressure: one situated within latitudinal belts characterized by high buoyancy anomalies around the equatorial zone, and the other positioned over regions marked by high zonal velocity winds. Exponential function of $P_{\text{vap}}(|u_1^n|)$ (Equation 7b) and nonlinear function of $P_{\text{vap}}(T_1)$ (Equation 7a) helps to distinct between two major drivers of heating and wind power, which is consistent with real world. The intermediate regions sandwiched between these two bands exhibit variable precipitation characteristics, characterized by localized zones of reduced precipitation where condensation processes are triggered by convergent zones induced by advected water vapors or the emergence of water vapor from free convection phenomena.

The validation procedure has primarily focused on reconciling precipitation data from the Aeolus 2.0 model with marine regions. However, it is imperative to

recognize that, unlike the Southern hemisphere, the Northern hemisphere encompasses vast terrestrial landmasses that significantly exceed the extent of oceanic surfaces. This spatial incongruity underscores the necessity for a nuanced and contextually informed approach in calibrating precipitation estimates over terrestrial landmasses. The current scheme, in its simplified version, does not incorporate factors influencing land surface evaporation, including evapotranspiration and vegetation cover. The alignment observed in the Kernel Density Estimation (KDE) outcomes of monthly precipitation averages obtained from both the Aeolus 2.0 model and ERA5 observational data serves as a benchmark for validating the model's bulk aerodynamic scheme (Figure 6).

4 | CONCLUSION

In this study, we introduced and evaluated a novel bulk aerodynamic scheme for estimating sea surface evaporation

<figure><img src="image_13.png" /><figcaption>FIGURE 6 The Kernel Density Estimation (KDE) with Scott's Rule Bandwidth compares the density functions of the monthly averages for March, June, September, and December of the year 1980 between the Aeolus 2.0 model and ERA5 data. Total precipitation is measured in (m month⁻¹).</figcaption></figure>
<figure><img src="image_14.png" /></figure>

flux implemented within the mcTRSW dynamical core of the Aeolus 2.0 model. The scheme's performance was assessed across four distinct seasons, encompassing the summer solstice, winter solstice, and two equinoxes. This approach aimed to capture various seasonal regimes and solar geometries, ensuring a robust examination of its efficiency under diverse climatological conditions.

The novel bulk aerodynamic scheme builds upon foundational principles established in the historical context of fluid dynamics and meteorology. At its core, the scheme leverages the bulk aerodynamic method to estimate air-sea fluxes, particularly focusing on turbulent heat and moisture exchange. Departing from conventional approaches, our scheme introduces a modified formulation for estimating sea surface evaporation flux, accounting for the influence of three primary regimes or forces: the heating force, the zonal wind force, and free convection. By incorporating these regimes, the proposed scheme provides a representation of the interactions between the atmosphere and the ocean. The scheme utilizes exponential functions of the potential temperature and velocity field of the atmosphere in the lower troposphere. The equilibrium vapor pressure, crucial for estimating evaporation, is determined through a modified version of the Clausius–Clapeyron equation, augmented by an adjustable scaling factor to account for ocean–atmosphere interactions. Additionally, the scheme incorporates the effects of wind velocity on evaporation, recognizing the significant role of wind forcing in facilitating moisture transport from oceanic regions to land masses. Our model initialization process utilized ERA5 reanalysis data from 1980, providing high-quality global atmospheric information to initialize the model fields. This strategic initialization ensured the accurate representation of atmospheric conditions at the model's inception, facilitating the expeditious capture of angular momentum advection by meridional circulations. The validation procedure involved comparing the model outputs with observational data, particularly focusing on reconciling precipitation data over marine regions. The alignment observed in the KDE outcomes of monthly precipitation averages between the Aeolus 2.0 model and ERA5 data served as a benchmark for validating the bulk aerodynamic scheme. The comparative analysis of zonally-averaged total precipitation revealed prominent concentrations in the equatorial zone and near the peak of zonal wind, consistent with distinct atmospheric regimes dictated by water vapor sources. However, it is essential to acknowledge the limitations of the current scheme, particularly its simplified representation of land surface processes. Unlike marine regions, terrestrial landmasses pose challenges due to factors such as evapotranspiration and vegetation cover, which are not fully accounted for in the model's formulation.

Future research will focus on refining the scheme to incorporate land surface processes and enhance accuracy by coupling Aeolus 2.0 with the Potsdam Earth Model (POEM) land (Drüke et al., 2021) and ice sheet (Kreuzer et al., 2021) components.

## AUTHOR CONTRIBUTIONS

**Masoud Rostami:** Conceptualization; data curation; formal analysis; visualization; writing – original draft; writing – review and editing; supervision; investigation; methodology; software; validation; funding acquisition; resources. **Stefan Petri:** Conceptualization; data curation; writing – original draft; methodology; investigation; software; formal analysis; validation; writing – review and editing. **Bijan Fallah:** Conceptualization; methodology; software; validation; investigation; data curation; visualization; writing – review and editing. **Farahnaz Fazel-Rastgar:** Writing – review and editing; investigation; methodology; conceptualization; data curation; visualization; software; validation.

## ACKNOWLEDGEMENTS

The authors acknowledge Stefan Rahmstorf and Georg Feulner for stimulating discussions. The support from Virgin Unite USA, Inc. through the Planetary Boundary Science Lab project, the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research, and the Land Brandenburg is gratefully acknowledged. The authors also appreciate the provision of resources on the high-performance computer system at the Potsdam Institute for Climate Impact Research, which has been instrumental in advancing their work. Open Access funding enabled and organized by Projekt DEAL.

## CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.

## DATA AVAILABILITY STATEMENT

The Aeolus 2.0 model and the corresponding parameters used in our simulations can be accessed in the Zenodo repository at https://doi.org/10.5281/zenodo.10054154.

## ORCID

Masoud Rostami https://orcid.org/0000-0003-1730-5145

Stefan Petri https://orcid.org/0000-0002-4379-4643

## REFERENCES

Bonino, G., Iovino, D., Brodeau, L. & Masina, S. (2022) The bulk parameterizations of turbulent air–sea fluxes in NEMO4: the origin of sea surface temperature differences in a global model study. *Geoscientific Model Development*, 15, 6873–6889.

Brutsaert, W. & Yu, S.L. (1968) Mass transfer aspects of pan evaporation. *Journal of Applied Meteorology and Climatology*, 7, 563–566.
<figure><img src="image_15.png" /></figure>

Dalton, J. (1798) Experimental essays, on the constitution of mixed gases; on the force of steam or vapour from water and other liquids in different temperatures, both in a {Toricellian} vacuum and in air; on evaporation; and on the expansion of gases by heat. *Memoirs and Proceedings of the Literary and Philosophical Society of Manchester*, 5, 535–602.

Drüke, M., von Bloh, W., Petri, S., Sakschewski, B., Schaphoff, S., Forkel, M. et al. (2021) CM2Mc-LPJmL v1.0: biophysical coupling of a process-based dynamic vegetation model with managed land to a general circulation model. *Geoscientific Model Development*, 14, 4117–4141.

Fairall, C.W., Bradley, E.F., Rogers, D.P., Edson, J.B. & Young, G.S. (1996) Bulk parameterization of air-sea fluxes for tropical ocean-global atmosphere coupled-ocean atmosphere response experiment. *Journal of Geophysical Research: Oceans*, 101, 3747–3764. Available from: https://doi.org/10.1029/95JC03205

Held, I.M. (2005) The gap between simulation and understanding in climate modeling. *Bulletin of the American Meteorological Society*, 86, 1609–1614.

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J. et al. (2020) The ERA5 global reanalysis. *Quarterly Journal of the Royal Meteorological Society*, 146, 1999–2049. Available from: https://doi.org/10.1002/qj.3803

Huang, W., Chen, F., Feng, S., Chen, J. & Zhang, X. (2013) Inter-annual precipitation variations in the mid-latitude Asia and their association with large-scale atmospheric circulation. *Chinese Science Bulletin*, 58, 3962–3968. Available from: https://doi.org/10.1007/s11434-013-5970-4

Huang, W., Chen, J., Zhang, X., Feng, S. & Chen, F. (2015) Definition of the core zone of the “westerlies-dominated climatic regime”, and its controlling factors during the instrumental period. *Science China Earth Sciences*, 58, 676–684. Available from: https://doi.org/10.1007/s11434-015-5057-y

Iacono, M.J., Mlawer, E.J., Clough, S.A. & Morcrette, J.-J. (2000) Impact of an improved longwave radiation model, RRTM, on the energy budget and thermodynamic properties of the NCAR community climate model, CCM3. *Journal of Geophysical Research: Atmospheres*, 105, 14873–14890. Available from: https://doi.org/10.1029/2000JD900091

Iribarne, J.V. & Godson, W.L. (1973) *Atmospheric thermodynamics*, 1st edition. Dordrecht: Springer.

Jiang, X., Cai, F., Li, Z., Wang, Z. & Zhang, T. (2023) The westerly winds control the zonal migration of rainy season over the Tibetan plateau. *Communications Earth & Environment*, 4, 363. Available from: https://doi.org/10.1038/s43247-023-01035-6

Kreuzer, M., Reese, R., Huiskamp, W.N., Petri, S., Albrecht, T., Feulner, G. et al. (2021) Coupling framework (1.0) for the PISM (1.1.4) ice sheet model and the MOM5 (5.1.0) ocean model via the PICO ice shelf cavity model in an Antarctic domain. *Geoscientific Model Development*, 14, 3697–3714.

Massey, F.J. (1951) The Kolmogorov-Smirnov test for goodness of fit. *Journal of the American Statistical Association*, 46, 68–78. Available from: https://doi.org/10.1080/01621459.1951.10500769

Mlawer, E.J., Taubman, S.J., Brown, P.D., Iacono, M.J. & Clough, S.A. (1997) Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for the longwave. *Journal of Geophysical Research: Atmospheres*, 102, 16663–16682. Available from: https://doi.org/10.1029/97JD00237

Monin, A. & Obukhov, A. (1954) Basic laws of turbulent mixing in the ground surface layer. *Trudy Geofizicheskogo Instituta Akademii Nauk SSSR*, 151, 163–187.

Monteith, J.L. (1981) Evaporation and surface temperature. *Quarterly Journal of the Royal Meteorological Society*, 107, 1–27. Available from: https://doi.org/10.1002/qj.49710745102

Morton, B.R., Taylor, G. & Turner, J.S. (1956) Turbulent gravitational convection from maintained and instantaneous sources. *Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences*, 234, 1–23.

Penman, H.L. (1948) Natural evaporation from open water, bare soil and grass. *Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences*, 193, 120–145.

Pond, S., Fissel, D.B. & Paulson, C.A. (1974) A note on bulk aerodynamic coefficients for sensible heat and moisture fluxes. * boundary-Layer Meteorology*, 6, 333–339. Available from: https://doi.org/10.1007/BF00232493

Roll, H. (1965) *Physics of the marine atmosphere*. Academic Press.

Rostami, M., Petri, S., Guimarães, S. & Fallah, B. (2024) Open-source stand-alone version of atmosphere model Aeolus 2.0 software. *Geoscience Data Journal*, 11, 1086–1093. Available from: https://doi.org/10.1002/gdj3.249

Rostami, M., Severino, L., Petri, S. & Hariri, S. (2024) Dynamics of localized extreme heatwaves in the mid-latitude atmosphere: a conceptual examination. *Atmospheric Science Letters*, 25, e1188. Available from: https://doi.org/10.1002/asl.1188

Rostami, M., Zhao, B. & Petri, S. (2022) On the genesis and dynamics of madden-Julian oscillation-like structure formed by equatorial adjustment of localized heating. *Quarterly Journal of the Royal Meteorological Society*, 148, 3788–3813.

Smith, S.D. (1988) Coefficients for sea surface wind stress, heat flux, and wind profiles as a function of wind speed and temperature. *Journal of Geophysical Research: Oceans*, 93, 15467–15472. Available from: https://doi.org/10.1029/JC093iC12p15467

Wang, Y., Tao, W.-K. & Simpson, J. (1996) The impact of ocean surface fluxes on a TOGA COARE convective system. *Monthly Weather Review*, 124, 2753–2763.

Warneford, E.S. & Dellar, P.J. (2013) The quasi-geostrophic theory of the thermal shallow water equations. *Journal of Fluid Mechanics*, 723, 374–403. Available from: https://doi.org/10.1017/jfm.2013.101

Wei, W., Zhang, R., Wen, M. & Yang, S. (2017) Relationship between the Asian westerly jet stream and summer rainfall over central Asia and North China: roles of the Indian monsoon and the South Asian high. *Journal of Climate*, 30, 537–552.

Wirangga, R., Mugisidi, D., Sayuti, A.T. & Heriyani, O. (2023) The impact of wind speed on the rate of water evaporation in a desalination chamber. *Journal of Advanced Research in Fluid Mechanics and Thermal Sciences*, 106, 39–50.

**How to cite this article:** Rostami, M., Petri, S., Fallah, B., & Fazel-Rastgar, F. (2025). A novel sea surface evaporation scheme assessed by the thermal rotating shallow water model. *Atmospheric Science Letters*, 26(1), e1287. https://doi.org/10.1002/asl.1287
APPENDIX A: MODEL DESCRIPTION

Following Rostami et al. (2022), the governing equations for the two-layer mcTRSW model are formulated as follows:

$$ (\partial_t + \mathbf{v}_1 \cdot \nabla) \mathbf{v}_1 + f \hat{\mathbf{z}} \times \mathbf{v}_1 = -\langle \nabla p_1 \rangle, \quad (8a) $$

$$ (\partial_t + \mathbf{v}_2 \cdot \nabla) \mathbf{v}_2 + f \hat{\mathbf{z}} \times \mathbf{v}_2 = -\langle \nabla p_2 \rangle - \frac{1-\gamma}{b_2 h_2} (\mathbf{v}_2 - \mathbf{v}_1)(C - D), \quad (8b) $$

$$ \partial_t h_1 + \nabla \cdot (h_1 \mathbf{v}_1) = \frac{1}{b_1} [(1-\gamma)(-C+D) - (1-\gamma^F)F_1], \quad (8c) $$

$$ \partial_t h_2 + \nabla \cdot (h_2 \mathbf{v}_2) = \frac{1}{b_2} [(1-\gamma)(+C-D) - (1-\gamma^F)F_2], \quad (8d) $$

$$ \partial_t b_1 + \mathbf{v}_1 \cdot \nabla b_1 = \frac{1}{h_1} [(+C - \mu \mathcal{E}) + F_1], \quad (8e) $$

$$ \partial_t b_2 + \mathbf{v}_2 \cdot \nabla b_2 = \frac{1}{h_2} [(-C + D) + F_2], \quad (8f) $$

$$ \partial_t q_1 + \nabla \cdot (q_1 \mathbf{v}_1) = -C + \mathcal{E}, \quad (8g) $$

$$ \partial_t q_2 + \nabla \cdot (q_2 \mathbf{v}_1) = +C - D. \quad (8h) $$

The Lagrangian derivative, denoted by $\partial/\partial t + \mathbf{v}_i \cdot \nabla$, shows the temporal evolution of quantities along the trajectory of fluid parcels. The presence of convective mass flux induces a Stokes drag, as expressed in Equation (8b), which arises from the vertical averaging of horizontal momentum equations with convective vertical velocity across the material interface. Here, $b_i$ represents the layer-averaged buoyancy, defined as $g\theta_i/\theta_s$, where $\theta_i$ is the mean horizontal potential temperature in layer $i$, $g$ is the gravitational acceleration, and $\theta_s$ is the surface-level potential temperature. The indices $i=1,2$ correspond to the lower and upper layers, respectively. $\mathbf{v}_i = (u_i, v_i)$ denotes the horizontal velocity vector, and $\langle \nabla p_i \rangle = \nabla \tilde{p}_i - \tilde{h}_i \nabla b_i$ represents the vertically-averaged pressure gradient. The Coriolis parameter is $f$, and $\hat{\mathbf{z}}$ is a unit vector in the vertical direction. The quantities $\tilde{p}_1, \tilde{p}_2, \tilde{h}_1$, and $\tilde{h}_2$ are defined as: $\tilde{p}_1 = (h_b + h_1 + h_2)b_1$, $\tilde{p}_2 = h_1 b_1 + (h_b + h_2)b_2$, $\tilde{h}_1 = h_b + (1/2)h_1 + h_2$, $\tilde{h}_2 = h_b + (1/2)h_2$, where $h_1, h_2$, and $h_b$ are the thicknesses of the lower layer, upper layer, and bottom topography, respectively.

$C$ represents the rate of condensation, $q_1$ denotes the actual moisture content, and $Q_s$ signifies the saturation moisture value. The Heaviside step function $H(q_1 - Q^s)$ ensures that condensation occurs only when the actual moisture content surpasses the saturation threshold. Mathematically, $C$ is expressed as:

$$ C = \frac{q_1 - Q_s}{\tau_c} \mathbb{H}(q_1 - Q^s). \quad (9) $$

$\mathcal{E} = S_{\text{vapor}}$ is sea surface evaporation according to Equation (7d) which acts as a cooling effect on the right-hand side of $db_1/dt$. Condensed water persists in the atmosphere as clouds, and precipitation occurs only when water droplets reach a critical size. Precipitable water, represented by Equation (10), has been incorporated into the model as an advected quantity denoted by $W$:

$$ \partial_t W_i + \nabla \cdot (W_i \mathbf{v}_i) = +C_i - V_i - \mathcal{P}. \quad (10) $$

This includes a source term due to condensation and a sink term due to vaporization ($V$) for each layer. The precipitation sink, $\mathcal{P}$, is introduced through a relaxation process with relaxation time $\tau_p$ when the bulk of precipitable water in the atmospheric column reaches a critical amount:

$$ \mathcal{P} = \frac{W_1 - W_{cr}}{\tau_p} \mathbb{H}(W_1 - W_{cr}). \quad (11) $$

The model incorporates two external forcings affecting $F_i$. First, insolation and its associated radiative transfer are accounted for using the rapid radiative transfer model (RRTM) (Iacono et al., 2000; Mlawer et al., 1997). Second, a thermal Newtonian cooling term, $-(h_i b_i - H_i B_i)/\tau_r$, is included, where $\tau_r$ denotes the relaxation time, and $H_i$ and $B_i$ signify the thickness and buoyancy at the rest state, respectively.

In the mass equations, $\mathcal{D}$ represents the downdraft from the upper layer to the lower one. The terms $\pm(C - D)$ on the right-hand side of Equation (4) account for both adiabatic convective effects (see Equations 8c and 8d) and the terms $(-C + D)$, $(+C - \mathcal{E})$, which represent diabatic nonconvective effects (see Equations 8e and 8f) acting on each layer. The coefficients $\gamma$ modulate the strength of these convective fluxes. When $\gamma = 1$, external forcing and latent heat release from condensation increase the mean potential temperature of the fluid layer, rather than venting convective plumes through the upper boundary, thus yielding a fully diabatic, nonconvective model similar to traditional approaches (Warneford & Dellar, 2013). Although the coefficients $\gamma$
<figure><img src="image_16.png" /></figure>

and $\gamma^F$ could vary between layers, we assume them to be uniform for simplicity.

The schematic representation of a spin-weighted spherical harmonic transform for shallow water equations and the nondimensionalization process are

discussed in Rostami et al. (2022); Rostami, Petri, et al. (2024). The Aeolus 2.0 model parameters used in our simulations are available in the ZENODO repository at https://doi.org/10.5281/zenodo.10054154.