DiuSST: a conceptual model of diurnal warm layers for idealized atmospheric simulations with interactive sea surface temperature

The diurnal variability in sea surface temperature (SST) may play an important role in cloud organization above the tropical ocean, with implications for precipitation extremes, storminess, and climate sensitivity. Recent cloud-resolving simulations demonstrate how imposed diurnal SST oscillations can strongly and delicately impact mesoscale convective organization. In spite of this nuanced interaction, many idealized modeling studies of tropical convection either assume a constant, homogeneous SST or, in the case of a responsive sea surface, represent the upper ocean by a slab with fixed thickness. Here we show that slab ocean models with constant heat capacity fail to capture the wind-dependence of observed diurnal sea surface warming. To alleviate this shortcoming, we present a simple yet explicitly depth-resolved model of upper-ocean temperature dynamics under atmospheric forcing. Our modular scheme describes turbulent mixing as diffusion with a wind-dependent diffusivity, in addition to a bulk mixing term and heat fluxes entering as sources and sinks. Using observational data, we apply Bayesian inference to calibrate the model. In contrast to a slab model, our model captures the exponential reduction in the diurnal warming amplitude with increasing wind speed. Further, our model performs comparably to a more elaborately parameterized diurnal warm-layer model. Formulated as a single partial differential equation with three key tuning parameters, the model is a suitable interactive numerical boundary condition for idealized atmospheric simulations.

1 Introduction

The role of clouds in a changing climate remains an open question that challenges predictions of climate sensitivity, regional precipitation patterns, and extreme weather events (Bony et al., 2015; Siebesma et al., 2020). Fundamental processes in cloud dynamics remain insufficiently understood, particularly how convective clouds cluster and interact with their environment. This knowledge gap calls for idealized modeling approaches simple enough to distill mechanisms yet relatable to the real world.

In the tropics, clouds organize across a wide range of spatial and temporal scales (Moncrieff, 2010). Diurnally, thunderstorms often cluster in mesoscale convective systems (MCSs), which are associated with extreme precipitation and the genesis of tropical cyclones (Tan et al., 2015; Schumacher and Rasmussen, 2020). Intraseasonally, variability is dominated by the Madden–Julian Oscillation (MJO), an eastward-propagating zone of strong deep convective activity (Madden and Julian, 1972; Zhang, 2005). Although the MJO is known to couple to the large-scale circulation and impact weather around the globe, it remains difficult to model (Jiang et al., 2020; DeMott et al., 2015).

Due to the multiscale interaction of tropical convection, small-scale processes may be key to understanding large-scale patterns (Slingo et al., 2003). Indeed, it is increasingly acknowledged that the diurnal variability in sea surface temperature (SST) can play an important role in atmospheric dynamics (Li et al., 2001; Clayson and Chen, 2002; Bernie et al., 2008; Bellenger et al., 2010; Haerter et al., 2020). Observations draw a clear link between strong diurnal SST oscillations and a diurnal cycle of marine cumulus convection (Johnson et al., 1999). Furthermore, the diurnal cycle of SST, by enhancing air–sea heat transfer, could help trigger the active phase of the MJO (Seo et al., 2014; Woolnough et al., 2007; Zhang, 2005; Zhao and Nasuno, 2020; Karlowska et al., 2023). Neglecting diurnal SST variations may lead to a bias on the order of 10 W m⁻² in monthly averaged surface heat fluxes (Weihs and Bourassa, 2014) and even stronger biases on shorter timescales. These findings emphasize the relevance of resolving diurnal air–sea interactions when modeling convective organization.

However, many idealized modeling studies of tropical convection prescribe a constant SST in space and time. A popular modeling framework is radiative–convective equilibrium (RCE), where constant solar forcing is balanced by outgoing longwave radiation above a constant, homogeneous SST (Tompkins and Craig, 1998). Under RCE, the moisture field is known to spontaneously self-organize into convective clusters separated by extended dry regions (Bretherton et al., 2005). This mechanism, termed convective self-aggregation, has been associated with real-world features such as MCS formation. Yet, self-aggregation is hampered by the realistic limit of fine horizontal model resolution, calling the realism of the mechanism into question (Yanase et al., 2020). Recent studies imposing a diurnal oscillation of SST demonstrate the emergence of diurnal self-aggregation even at a fine spatial resolution (Haerter et al., 2020; Jensen et al., 2022). Likewise, spatial variations in SST have been shown to imprint themselves on the moisture field (Müller and Hohenegger, 2020; Shamekh et al., 2020b; Skyllingstad et al., 2019). Since SST variability can substantially alter the spatiotemporal patterns of marine tropical convection, there is a need to incorporate higher-fidelity SST representations when investigating convective processes.

Over the past years, the modeling community studying convective organization has largely addressed the issue of interactive SSTs by coupling the atmosphere to a single-layer slab with fixed heat capacity, which absorbs and re-emits heat according to parameterized surface fluxes (Hohenegger and Stevens, 2016; Shamekh et al., 2020a; Coppin and Bony, 2017; Tompkins and Semie, 2021; Wing et al., 2017). These works report overall that a responsive SST slows down or even prevents the onset of convective aggregation. However, in this paper we show that slab models do not adequately produce realistic diurnal SST warming: they fail to capture its wind dependence by neglecting upper-ocean mixing. Yet wind can by no means be neglected, even in strongly idealized studies, given that cold pools, bringing gusty winds, are increasingly appreciated as an integral mechanism in convective self-organization (Tompkins, 2001; Böing, 2016; Haerter, 2019; Nissen and Haerter, 2021). To address this inconsistency, we present a simplified one-dimensional model that significantly improves the representation of diurnal SST variability compared to a slab ocean, while being adaptable and affordable for coupling to cloud-resolving atmospheric models.

SST plays a key role in governing the heat and moisture exchange between the atmosphere and ocean (Seo et al., 2023). While the observed diurnal amplitudes of surface temperature are typically largest over land, a diurnal cycle of SST is also common throughout the tropical ocean (Kawai and Wada, 2007). Under strong insolation and calm conditions, a diurnal warm layer forms during the day, raising skin SST (measured directly at the surface) by up to 3–4 K. Field studies have observed extreme diurnal warming events of 5 K or more, although these often lie in the extratropics or in coastal regions (Gentemann et al., 2008; Minnett, 2003; Ward, 2006). Importantly, the amplitude of diurnal warming is highly sensitive to wind speed, with observations suggesting an approximately exponential decay of diurnal warming with increasing wind speed (Gentemann et al., 2003; Börner, 2021).

To the first order, the diurnal warm layer results from a competition between insolation-driven thermal stratification and wind-driven mixing. On a calm, clear day after sunrise, the incident solar radiation quickly heats up the upper ocean, creating a stably stratified density profile. This stratification leads to suppressed vertical heat exchange, thus trapping further heat near the surface in a positive feedback loop that can produce strong surface warming. Later, in the early afternoon when the net surface heat flux changes sign, surface cooling causes unstable density stratification and initiates vertical mixing, resulting in a deepening of the diurnal warm layer. At night, convective mixing typically acts to “reset” the temperature profile. Furthermore, wind stress induces turbulent mixing in the water column. Known as the cool-skin effect, the molecular skin layer at the sea surface is typically a few tenths of a degree colder than the water a millimeter deeper (Donlon et al., 2002; Wong and Minnett, 2018).

Upper-ocean heat transfer has been modeled at various levels of complexity over the last 50 years, ranging from fully turbulence-resolving models to simplified bulk models and empirical models (Kawai and Wada, 2007). The class of turbulence closure models describes the dynamics of the oceanic boundary layer in terms of diffusion (Large et al., 1994; Kondo et al., 1979; Noh and Jin Kim, 1999), second moment closure (Mellor and Yamada, 1982; Kantha and Clayson, 1994), or an eddy spectrum (Stull and Kraus, 1987). Although such ocean models have been coupled to a weather forecasting model (Noh et al., 2011), they are arguably too costly and complex for the purpose of idealized atmospheric simulations at high horizontal resolution. The simplest models, derived from empirical relations, specialize in estimating the daily SST amplitude (Webster et al., 1996; Price et al., 1987; Kawai and Kawamura, 2002) or its hourly evolution (Li et al., 2001; Zeng et al., 1999; Gentemann et al., 2003) based on averaged atmospheric data. These models are not designed for the high temporal resolution of atmospheric large-eddy simulations (on the order of seconds) and are sometimes not physics informed.

Alternatively, simplified models based on boundary layer physics have been developed in the spirit of subdividing the water column into layers described by bulk dynamical equations. A popular candidate, which we compare against here, is the prognostic ZB05 scheme by Zeng and Beljaars (2005), with improvements by Takaya et al. (2010). It computes the sea skin temperature from an integrated mixed layer equation combined with a cool-skin scheme, assuming a power law temperature profile. Other bulk models include the model by Price et al. (1986) and its developments (Fairall et al., 1996, 2003; Gentemann et al., 2009; Schiller and Godfrey, 2005). Widely applied in weather and climate modeling (Brunke et al., 2008), such models could offer a suitable balance between physical rigor and computational cost. Nonetheless, they do not necessarily cater to the needs of idealized atmospheric modeling. First, some bulk models require unknown input about the oceanic background state. Second, the integral descriptions of heat and momentum transfer are typically based on elaborate parameterizations, making the models less intuitive and difficult to tune during sensitivity experiments. Third, these models do not resolve the vertical temperature profile (apart from the parameterized profile in Gentemann et al., 2009), limiting their use as an interface for sub-surface processes that might be sensitive to diurnal warming, such as microbial ecosystems (Wurl et al., 2017).

In this paper, we present an idealized model of diurnal sea surface warming that explicitly resolves temperature in the upper few meters of the ocean while being conceptually simple and efficient. Derived from first principles as a modified heat equation, the model consists of a single partial differential equation controlled by three key parameters, taking insolation and atmospheric conditions as forcing input. After describing the model (Sect. 2), we use a Bayesian approach to calibrate and evaluate the model based on cruise observations from the tropical eastern Pacific (Sect. 3). We compare the performance of our model to both a typical slab model and the ZB05 scheme (Sect. 4). While slab models fail to describe diurnal SST variability, our model accurately reproduces key features of diurnal warming, such as its wind dependence, near-surface heat trapping, and skin cooling. Our model performs comparably to the ZB05 scheme, which is more rigorously derived but does not resolve temperature vertically. In Sect. 4.4, we discuss the implications of our results for the atmospheric heat and moisture budget. We believe that our modular scheme can be useful as both a wind-responsive ocean surface for cloud-resolving modeling and an interface to study further air–sea–biosphere interactions in an idealized setting.

2 Model description

Our approach relies on the following basic assumptions. First, we assume that the oceanic, atmospherically driven diurnal temperature variability is constrained to within a few meters of the sea surface. We thus define the foundation depth z_f at which diurnal temperature changes become negligible. Conceptually, this partitions the water column into the diurnal layer above z_f and the foundation layer below. Second, we assume that we may separate the scales of diurnal warming from the ocean's slow internal temperature variability. The foundation layer may then be considered an infinite heat reservoir with constant foundation temperature T_f. Third, we neglect any horizontal inhomogeneities and flows, allowing us to treat the problem as one-dimensional. The goal is then to determine the evolution of the vertical temperature profile within the diurnal layer in response to a given time sequence of atmospheric forcing. A schematic sketch of the setup is depicted in Fig. 1.

Figure 1 Schematic of the simplified upper-ocean model, illustrating the processes considered.

2.1 Main equation

We consider the sea temperature T(z,t) as a function of time t≥0 and the vertical coordinate $z\in (z_{\mathrm{f}},\mathrm{0}]$, where z=0 defines the air–sea interface and z<0 is below the surface. Given an initial profile T(z,0), foundation temperature T_f, and atmospheric forcing ℱ(t), we propose a single prognostic partial differential equation (PDE) for the time evolution of $T\equiv T(z,t)$:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial T}{\partial t}}=\underset{\text{diffusion}}{\underbrace{{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\kappa }(z,t){\displaystyle \frac{\partial T}{\partial z}}\right)}}-\underset{\text{mixing}}{\underbrace{\mathit{\mu }{\displaystyle \frac{T-T_{\mathrm{f}}}{z-z_{\mathrm{f}}}}}}+\underset{\text{source/sink}}{\underbrace{{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{p}}}}{\displaystyle \frac{\partial Q(z,t)}{\partial z}}}}.\end{array}$$

Here κ>0 represents the (time- and depth-dependent) diffusivity, μ>0 is a constant that we term the mixing coefficient, and the constants ρ_w and c_p denote the density and specific heat capacity at the constant pressure of seawater, respectively. Finally, Q is the net vertical heat flux at depth z, defined as positive downwards (into the ocean). Explicit expressions of these quantities in relation to the forcing ℱ follow below.

The three terms on the right-hand side of Eq. (1) provide strongly idealized representations of different physical processes. Motivated by the heat equation, the first term describes the vertical diffusion of heat due to turbulent eddies within the diurnal layer. In addition, internal processes such as convection may cause additional vertical mixing of water masses between the diurnal layer and the foundation layer. This is crudely incorporated by the second term, which relaxes T to the foundation temperature with relaxation timescale $t_{\mathit{\mu }}=(z-z_{\mathrm{f}})/\mathit{\mu }$. Here the depth dependence reflects the supposition that water near the foundation depth will mix faster with foundation layer water compared to water near the surface. Lastly, the third term comprises all sources and sinks of heat, both at the air–sea interface and within the diurnal layer.

The foundation temperature $T(z_{\mathrm{f}},t)=T_{\mathrm{f}}\forall t$ acts as a Dirichlet boundary condition at z_f. At the air–sea interface (z=0), the sea temperature evolves according to the net surface heat flux $Q_{\mathrm{0}}\left(t\right)\equiv Q(\mathrm{0},t)$ (see Sect. 2.3), which closes the energy budget.

2.2 Wind-driven mixing and stratification

Wind stress at the sea surface induces shear instability, causing vertical turbulent heat transport in the upper ocean. This implies that the diffusivity κ(z,t) in our model should depend on the wind speed u(t). Since the magnitude of wind stress on a water surface is approximately proportional to the square of the wind speed (Smith, 1988; Edson et al., 2013), we propose to model the diffusivity κ as

$$\begin{array}{}\text{(2)}& \mathit{\kappa }(z,t):={\mathit{\kappa }}_{\text{mol}}+{\mathit{\kappa }}_{\mathrm{0}}\mathit{\phi }\left(z\right){\left({\displaystyle \frac{u\left(t\right)}{u_{\mathrm{0}}}}\right)}^{\mathrm{2}},\end{array}$$

where κ_mol and κ₀ are coefficients of molecular and (vertical) eddy diffusion, respectively. Typically, molecular heat conduction is negligible (κ_mol≪κ₀). Note the scaling of the eddy diffusion term with u², which we non-dimensionalize by the reference wind speed $u_{\mathrm{0}}=\mathrm{1}{\text{ m\hspace{0.17em}s}}^{-\mathrm{1}}$ simply to ensure that κ₀ has units of diffusivity.

The vertical diffusivity profile φ(z) in Eq. (2) approximates how turbulent mixing varies with depth. Physically, thermal density stratification inhibits turbulent heat transfer, since vertical mixing of a stratified fluid is energetically unfavorable. Additionally, the characteristic size of eddies decreases when approaching the air–sea boundary (Pope, 2000). Both effects motivate the suppression of the diffusivity near the surface, where stratification becomes largest during diurnal heating. For simplicity, we limit our study to a time-independent linear profile,

$$\begin{array}{}\text{(3)}& \mathit{\phi }\left(z\right)=\mathrm{1}+\mathit{\sigma }\left({\displaystyle \frac{z}{z_{\mathrm{f}}}}-\mathrm{1}\right),\end{array}$$

parameterized by the surface suppressivity $\mathit{\sigma }\in [\mathrm{0},\mathrm{1}]$. If σ>0, then φ increases linearly with depth from $\mathit{\phi }\left(\mathrm{0}\right)=\mathrm{1}-\mathit{\sigma }$ at the surface to φ(z_f)=1 at the foundation depth.

It may seem counterintuitive to decrease the diffusivity towards the surface when it is driven by surface winds. In this idealized setup, however, the diffusion term simultaneously accounts for density stratification (in a time-averaged sense). Furthermore, surface wind stress can accelerate the warm layer, leading to turbulence due to shear instability at its base (Hughes et al., 2020). Alternative choices for the diffusivity profile are discussed in Sect. 5 and in Appendix D.

2.3 Air–sea interaction

In our model, the state of the atmosphere enters as the forcing $\mathcal{F}\left(t\right)\equiv \left(R_{\mathrm{sw},\downarrow }\right(t),u(t),T_{\mathrm{a}}(t),q_{\mathrm{v}}(t\left)\right)$, comprising the downward solar, or “shortwave”, irradiance $R_{\mathrm{sw},\downarrow }$; horizontal wind speed u; air temperature T_a; and specific humidity q_v. These quantities refer to some reference height above the sea surface (usually 10 m, where they are typically measured (Friehe and Schmitt, 1976)). The diurnal layer interacts with the atmosphere through the absorption and reflection of shortwave radiation, the absorption and emission of thermal (“longwave”) radiation, and sensible and latent heat exchange. While shortwave radiation penetrates the sea surface and is absorbed at a range of depths, the other fluxes act within micrometers of the air–sea interface (Wong and Minnett, 2018).

At the air–sea interface, the net surface heat flux Q₀(t) entering the water body is given by

$$\begin{array}{}\text{(4)}& Q_{\mathrm{0}}\left(t\right)=R_{\text{sw}}\left(t\right)+R_{\text{lw}}\left(t\right)+Q_{\text{s}}\left(t\right)+Q_{\text{l}}\left(t\right),\end{array}$$

where $R_{\text{sw}}=(\mathrm{1}-\mathcal{R})R_{\text{sw,}\downarrow }$ is the penetrating shortwave irradiance after subtracting the fraction ℛ that is reflected at the sea surface from $R_{\text{sw},\downarrow }$ (see Appendix C). The net longwave radiative flux R_lw and sensible heat flux Q_s can point into or out of the ocean, whereas the latent heat flux Q_l is always negative. We use standard bulk formulae to describe the surface heat fluxes (Friehe and Schmitt, 1976; DeCosmo et al., 1996; Wells and King-Hele, 1990; Large and Yeager, 2009),

$$\begin{array}{}\text{(5a)}& {\displaystyle }{R}_{\text{lw}}\left(t\right)={\mathit{\sigma }}_{\text{SB}}\left(T_{\mathrm{a}}\right(t{)}^{\mathrm{4}}-T(\mathrm{0},t{)}^{\mathrm{4}}),\text{(5b)}& {\displaystyle }{Q}_{\text{s}}\left(t\right)={\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{p},\mathrm{a}}{C}_{\mathrm{s}}u\left(t\right)\left(T_{\mathrm{a}}\right(t)-T(\mathrm{0},t\left)\right),\text{(5c)}& {\displaystyle }{Q}_{\text{l}}\left(t\right)={\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{l}}{L}_{\mathrm{v}}u\left(t\right)\left(q_{\mathrm{v}}\right(t)-q_{\text{sat}}(T(\mathrm{0},t)\left)\right),\end{array}$$

which are given in terms of the Stefan–Boltzmann constant σ_SB; the density and specific heat capacity of air, ρ_a and c_p,a, respectively; and the Stanton number C_s, Dalton number C_l, and latent heat of vaporization L_v. We approximate these coefficients by constants based on the literature values (see Table 1), which are roughly valid for wind speeds on the order of 1 to 10 m s⁻¹ and air–sea temperature differences around 1 K, measured 10 m a.s.l. (meters above sea level) (Wells and King-Hele, 1990). Note that Eq. (5a) simply applies the Stefan–Boltzmann law for black-body radiation. When coupled to an atmospheric model, the downward longwave radiative flux can be taken from the atmospheric model output. Furthermore, the latent heat flux, Eq. (5c), involves the saturation specific humidity q_sat(T). The temperature dependence of q_sat obeys the Clausius–Clapeyron relation, which is approximated by the empirical formula

$$\begin{array}{}\text{(6)}& q_{\text{sat}}\left(T\right)\approx {\displaystyle \frac{\mathrm{611.2}}{{\mathit{\rho }}_{\mathrm{a}}{r}_{\mathrm{w}}T}}\exp \left({\displaystyle \frac{\mathrm{17.67}(T-\mathrm{273.15})}{T-\mathrm{29.65}}}\right),\end{array}$$

where r_w denotes the gas constant of water vapor, and T enters in units of Kelvin (Alduchov and Eskridge, 1996).

Table 1 Model constants and their default values used.

^a determined from data via Bayesian inference; see Table E1. ^b corresponds to mean 3 m water temperature from the fifth Marine Optical Characterization Experiment (MOCE-5) dataset.

We assume that the penetrating shortwave radiation R_sw is attenuated exponentially as it propagates downward through the water column,

$$\begin{array}{}\text{(7)}& Q(z,t)=R_{\text{sw}}\left(t\right)\exp \left({\displaystyle \frac{\mathit{\alpha }z}{\cos {\mathit{\varphi }}^{\prime }\left(t\right)}}\right)\text{for }z<\mathrm{0},\end{array}$$

where α>0 is the attenuation coefficient, and ϕ^′ denotes the refracted solar angle. It follows from the Sun's angle relative to the surface normal, $\mathit{\varphi }\in [\mathrm{0},\mathit{\pi }/\mathrm{2})$, by Snell's refraction law, that

$$\begin{array}{}\text{(8)}& {\mathit{\varphi }}^{\prime }\left(t\right)=\arcsin \left({\displaystyle \frac{n_{\mathrm{a}}}{n_{\mathrm{w}}}}\sin \mathit{\varphi }\left(t\right)\right),\end{array}$$

where n_w and n_a denote the refractive indices of seawater and air, respectively.

To summarize, our idealized model is controlled by three key parameters: (i) the eddy diffusivity κ₀ governs the magnitude of wind-driven turbulent heat diffusion, (ii) the mixing coefficient μ sets the timescale of relaxation to the foundation temperature, and (iii) the attenuation coefficient α regulates how deep shortwave radiation penetrates. Realistic values of these parameters are estimated from observational evidence (Sect. 3). We list all model constants and coefficients in Table 1.

2.4 Numerical implementation

To solve Eq. (1), we discretize the spatial coordinate z and numerically integrate the resulting system of ordinary differential equations in time (see Appendix A and Börner (2024) for a software implementation in Python). Spatial derivatives are approximated by second-order-accurate finite differences on a non-uniform vertical grid. We set the foundation depth to z_f=10 m, where diurnal temperature variability is mostly negligible (Kawai and Wada, 2007). The grid spacing is set to Δz₀=0.1 m at the sea surface and increases by a stretch factor ϵ≈1.04 with each consecutive grid point below, such that a total of N=40 grid points covers the diurnal layer $z\in (z_{\mathrm{f}},\mathrm{0}]$ (excluding the boundary point at z_f). Specifically, the depth z_n of the nth grid point is given by

$$\begin{array}{}\text{(9)}& z_n=-\mathrm{\Delta }{z}_{\mathrm{0}}\left({\displaystyle \frac{\mathrm{1}-{\mathit{ϵ}}^n}{\mathrm{1}-\mathit{ϵ}}}\right),n\in \mathit{\{}\mathrm{0},\mathrm{\dots },N\mathit{\}},\end{array}$$

with ϵ set such that z_N=z_f.

We implement the time integration as an explicit Euler scheme, which is numerically stable if the chosen time step Δt satisfies the Courant–Friedrichs–Lewy (CFL) condition,

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\mathcal{C}}_i:=\underset{n}{max}\left[\mathrm{2}{\displaystyle \frac{\mathit{\kappa }(z_n,t_i)\mathrm{\Delta }{t}_i}{(\mathrm{\Delta }{z}_n{)}^{\mathrm{2}}}}\right]\le \mathrm{1}\\ & n=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },N-\mathrm{1},\end{array}\end{array}$$

where n and i index the discrete space and time coordinates, respectively. Via the time-dependent diffusivity κ, the CFL condition depends on time, specifically on the current wind speed. To minimize computational cost, we use an adaptive time step Δt_i that maintains a maximal CFL number of 𝒞_i=0.95 at each instant of time t_i. Consequently, the time step Δt_i scales inversely with the square of the wind speed. We impose a cutoff wind speed $u_{\text{max}}=\mathrm{10}{\text{m\hspace{0.17em}s}}^{-\mathrm{1}}$, causing any wind speed u>u_max to be replaced by u_max when computing the diffusivity (Eq. 2). In the model simulations of this study, Δt_i ranged from 0.004 to 10 s, with a median time step of 5–6 s.

Comparing the explicit Euler scheme with analogous implementations of the fourth-order Runge–Kutta and implicit Euler methods, we find that the explicit Euler scheme is fastest unless the implicit method is used at a very large time step, in which case the solution lacks accuracy.

3 Model calibration and evaluation

Our model parameterizes diurnal warming in terms of three unknown constants: the diffusivity coefficient κ₀, mixing coefficient μ, and attenuation coefficient α. In principle, their values depend on the ocean properties at a given time and location. In this section, we first analyze how the parameters affect diurnal warming and then use observational data to estimate realistic values via Bayesian inference. Finally, we evaluate the performance of the calibrated model against observations.

3.1 Parameter sensitivity under idealized forcing

As an initial step, we perform a sensitivity study where we force the model with idealized atmospheric pseudo-data representing a calm and clear tropical day. Setting the foundation temperature to T_f=300 K, we let the air temperature T_a(t) oscillate harmonically around T_f, with a diurnal amplitude¹ of ΔT=2 K. Similarly, we impose a harmonically oscillating horizontal wind speed with a mean of $\stackrel{\mathrm{‾}}{u}=\mathrm{2}{\text{m\hspace{0.17em}s}}^{-\mathrm{1}}$ and amplitude $\mathrm{\Delta }u=\mathrm{2}{\text{m\hspace{0.17em}s}}^{-\mathrm{1}}$ that peak at midnight (this choice generates favorable conditions for diurnal warming and nighttime mixing). Such profiles read as follows:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& T_{\mathrm{a}}\left(t\right)=T_{\mathrm{f}}-{\displaystyle \frac{\mathrm{\Delta }T}{\mathrm{2}}}\cos (\mathrm{2}\mathit{\pi }t/t_{\mathrm{0}}),\\ & u\left(t\right)=\stackrel{\mathrm{‾}}{u}+{\displaystyle \frac{\mathrm{\Delta }u}{\mathrm{2}}}\cos (\mathrm{2}\mathit{\pi }t/t_{\mathrm{0}}),\end{array}\end{array}$$

where t is the local sun time (with respect to midnight) and t₀=1 d. The solar irradiance is given by

$$\begin{array}{}\text{(12)}& R_{\text{sw},\downarrow }\left(t\right)=\left\{\begin{array}{ll}-R_{\text{max}}\cos \left(\mathrm{2}\mathit{\pi }t/t_{\mathrm{0}}\right)& \text{if }\cos \left(\mathrm{2}\mathit{\pi }t/t_{\mathrm{0}}\right)<\mathrm{0}\\ \mathrm{0}& \text{otherwise,}\end{array}\right.\end{array}$$

where we set the peak insolation to $R_{\text{max}}=\mathrm{1000}{\text{W\hspace{0.17em}m}}^{-\mathrm{2}}$. Lastly, we fix the specific humidity at a constant value of $q_{\mathrm{v}}=\mathrm{15}{\text{ g\hspace{0.17em}kg}}^{-\mathrm{1}}$. Under these conditions, we may expect peak diurnal surface warming of around 2–3 K (Minnett, 2003).

Using the atmospheric pseudo-data described above, we now run 2 d long model simulations, consecutively varying one model parameter of the set $\mathit{\{}{\mathit{\kappa }}_{\mathrm{0}},\mathit{\mu },\mathit{\alpha }\mathit{\}}$ while fixing the other two at select default values (see Fig. 2). Increasing the eddy diffusivity κ₀ enhances the heat transport within the diurnal layer, flattening the vertical temperature gradient, diminishing surface warming, and advancing the time of peak surface warming. A decrease in the mixing coefficient μ corresponds to a slower removal of heat from the diurnal layer into the deeper ocean. This causes increased heat trapping in the upper ocean and thus stronger surface warming and a deeper warm layer. Particularly, for $\mathit{\mu }<{\mathrm{10}}^{-\mathrm{3}}{\text{m\hspace{0.17em}s}}^{-\mathrm{1}}$, excess heat remains in the diurnal layer throughout the night, accumulating heat on the next day. Setting μ to a sufficiently large value ensures that the temperature profile can reset at night, as is often observed (Kawai and Wada, 2007). Finally, the attenuation coefficient α determines the depth range at which solar radiation is absorbed. For $\mathit{\alpha }>\mathrm{1}{\text{m}}^{-\mathrm{1}}$, more than 60 % of radiation is absorbed within 1 m of the surface, leading to strong surface warming. Reducing α causes radiation to reach deeper depths, where it is more quickly transported into the foundation ocean via the mixing term.

Figure 2 Model sensitivity under variation in the parameters κ₀, μ, and α. In each sensitivity experiment, we force the model with 2 d of idealized atmospheric data (see main text), successively varying one parameter while fixing the other two. The top panels depict simulated surface warming ΔSST as a function of time for different values of the eddy diffusivity (a, d), mixing coefficient (b, e), and attenuation coefficient (c, f). The bottom panels show corresponding vertical temperature profiles at 02:00 (dashed lines, shifted by −1 K) and 14:00 local sun time (solid lines) on day 2 of the simulation. Values of the varied parameter are given in the respective legend; fixed parameter values are ${\mathit{\kappa }}_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{4}}{\text{m}}^{\mathrm{2}}{\text{s}}^{-\mathrm{1}}$, $\mathit{\mu }={\mathrm{10}}^{-\mathrm{3}}{\text{m\hspace{0.17em}s}}^{-\mathrm{1}}$, and $\mathit{\alpha }=\mathrm{3}{\text{m}}^{-\mathrm{1}}$. Further model settings are detailed in Table 1 (here T_f=300 K).

In Fig. 2, the similarity in model response between the parameters κ₀ and α suggests that they might play a similar dynamical role. However, this is not the case. The two parameters control different aspects of the air–sea coupling. The eddy diffusivity κ₀ determines the coupling to wind, whereas the attenuation coefficient α modulates the oceanic heat uptake due to solar radiation. Correlations between model parameters are discussed below (Sect. 3.3).

3.2 Observational data

To infer realistic values for the parameter set $\mathit{\{}{\mathit{\kappa }}_{\mathrm{0}},\mathit{\mu },\mathit{\alpha }\mathit{\}}$, we now conduct a case study where we force the model with real observational data and compare the modeled diurnal warming with the observed signal. Here we use cruise data (available from Börner, 2024) from the fifth Marine Optical Characterization Experiment (MOCE-5), conducted during October 1999 off the Mexican west coast (Minnett, 2003; Ward, 2006). The route went along the coast of the Baja California peninsula, both in the open Pacific Ocean and within the Gulf of California, thus including offshore as well as more coastal conditions (see Fig. 3a).

Figure 3 Observational dataset used in this study. (a) Travel route of the MOCE-5 cruise in the Pacific Ocean and Gulf of California, colored by the day since departure from near San Diego, USA. (b) Diurnal warming ΔSST by location, as recorded during the cruise (data points with higher ΔSST are enlarged). (c) Time series of observed diurnal warming, ΔSST, showing the individual data points. Gray-shaded intervals indicate the training data used for Bayesian inference. Panels (d) and (e) display the time series of the radiometric skin SST (blue), air temperature (red), horizontal wind speed (green), and downward shortwave irradiance (orange). Note that we omit data recorded after 16 October due to extended temporal gaps in the dataset.

The research vessel Melville was equipped with an infrared radiometer of type M-AERI (Marine–Atmospheric Emitted Radiance Interferometer; Minnett et al., 2001). This instrument provides precise measurements of sea skin temperature by detecting infrared radiation emitted from within micrometers of the ocean surface. Using skin SST for calibration rather than bulk SST measured at up to 1 m depth is crucial when modeling air–sea interactions because the atmosphere senses only the sea skin. Additionally, solar irradiance $R_{\text{sw},\downarrow }$, wind speed u, air temperature T_a, and water temperature at 3 m depth were recorded at time intervals of approximately 10 to 12 min (see Fig. 3). Unfortunately, the dataset does not contain air humidity. We therefore assume a constant specific humidity of $q_{\mathrm{v}}=\mathrm{15}{\text{ g\hspace{0.17em}kg}}^{-\mathrm{1}}$ throughout the time series, which may be considered typical of the tropical ocean.

Throughout the following analysis, following Minnett (2003), we define diurnal warming ΔSST as the temperature difference between the sea skin (that is, the water directly at the surface) and a reference depth d,

$$\mathrm{\Delta }\text{SST}\left(t\right):=T(\mathrm{0},t)-T(d,t),$$

where $d=-\mathrm{3}$ m for the present dataset, and T(0,t) is given by the radiometric SST measurements. Diurnal warming events exceeding 1 °C are observed on several days, with ΔSST reaching up to 5 °C on 13 October (Fig. 3c; see also Ward (2006)). The dataset also includes days without any substantial diurnal warming, such as on 2 and 9 October. These days correlate with relatively high wind speeds, while the strong warming events in the second week, such as on 10, 13, and 14 October, coincide with low winds, especially during midday. Note that the time series includes diurnal warming events at different locations (Fig. 3b) and covers a wide range of SST values from 290 to 305 K (Fig. 3d). Horizontal wind speeds (corrected for the ship's motion) did not exceed 10 m s⁻¹; downward solar irradiance peaked at around 800 to 900 W m⁻² each day.

While the MOCE-5 cruise data have a temporal resolution of around 10 min, our model requires a time step on the order of seconds to meet the CFL condition (Eq. 10). This necessitates interpolating the atmospheric data between data points; we apply linear interpolation with a wind-adaptive time, step as described in Sect. 2.4.

3.3 Bayesian parameter estimation

The MOCE-5 observations allow us to estimate the free parameters of our model using Bayesian inference. Given the data 𝒟, the conditional probability P(Θ|𝒟) that a certain value of the parameter set $\mathrm{\Theta }=\mathit{\{}{\mathit{\kappa }}_{\mathrm{0}},\mathit{\mu },\mathit{\alpha }\mathit{\}}$ represents the “true” model follows from Bayes' rule,

$$\begin{array}{}\text{(13)}& P\left(\mathrm{\Theta }\right|\mathcal{D})\propto P(\mathcal{D}\left|\mathrm{\Theta }\right)P\left(\mathrm{\Theta }\right),\end{array}$$

that is, the posterior distribution P(Θ|𝒟) is proportional to the product of the likelihood P(𝒟|Θ) and the prior distribution P(Θ). The likelihood quantifies how well the model with given parameter settings Θ describes the observed data, while the prior distribution encapsulates previous knowledge of suitable parameter ranges (Gelman et al., 2013).

To estimate the posterior distribution, we first partition the data time series into a 6 d training set and an ∼8 d validation set (see gray shading in Fig. 3c). The training data are used to evaluate the likelihood function $\mathcal{L}\left(\mathrm{\Theta }\right)\equiv P\left(\mathcal{D}\right|\mathrm{\Theta })$, which is computed in the following way. For a given Θ, the model is run using the time series of the variables $\mathcal{D}=\mathit{\{}{R}_{\text{sw},\downarrow },u,T_{\mathrm{a}},q_{\mathrm{v}}\mathit{\}}$ as atmospheric forcing. We set the foundation temperature T_f to the mean observed sea temperature at 3 m depth, T_f=298.19 K, and shift the air temperature T_a according to the current deviation of the 3 m temperature from its mean T_f. This maintains the observed air–sea temperature contrast at all times. Then, we compare the modeled time series of diurnal warming, ΔSST^model, to the observed diurnal warming signal ΔSST^obs in order to compute the likelihood,

$$\begin{array}{}\text{(14)}& \mathcal{L}\left(\mathrm{\Theta }\right)\propto \exp (-\sum _j{\displaystyle \frac{\left(\mathrm{\Delta }{\text{SST}}_j^{\text{model}}\right(\mathrm{\Theta })-\mathrm{\Delta }{\text{SST}}_j^{\text{obs}}{)}^{\mathrm{2}}}{{\mathrm{\Sigma }}_j^{\mathrm{2}}}}).\end{array}$$

Here the index j runs through all data points (over time) in the training set, and Σ_j denotes a weight (see Appendix E2).

Approximating the posterior distribution in parameter space is achieved via Markov chain Monte Carlo (MCMC) sampling using the emcee package in Python (Foreman-Mackey et al., 2013), based on an affine-invariant algorithm proposed by Goodman and Weare (2010). Further information on the Bayesian inference procedure and choice of prior is provided in Appendix E.

After convergence, the sampled posterior distribution shows a single peak in each parameter direction, indicating a well-defined optimal value for each parameter (Fig. 4). From the projected distributions, we compute the maximum, median, and mean value as parameter estimates (Table E1). The maximum value, or maximum a posteriori (MAP) estimate, corresponds to the most likely parameter value and will be used for the subsequent analysis.

Figure 4 Posterior distribution of the parameter set $\mathit{\{}{\mathit{\kappa }}_{\mathrm{0}},\mathit{\mu },\mathit{\alpha }\mathit{\}}$ obtained from Bayesian MCMC sampling. Plots on the upper diagonal are one-dimensional projections of the parameter space, showing for each parameter the prior (orange) and sampled posterior (black) distributions (vertical axis not to scale). Blue lines indicate the MAP values (Table E1). The two-dimensional projections of parameter space depict scatter points and smoothed isocontours of the posterior. Plotted using corner.py (Foreman-Mackey, 2016).

The two-dimensional projections of the posterior distribution (Fig. 4) indicate possible correlations between the three parameters. In particular, a positive correlation between κ₀ and α is visible. However, part of the correlation likely reflects the fact that the solar irradiance (modulated by α) and wind speed (modulated by κ₀) are themselves negatively correlated in the training data (not shown).

3.4 Model performance against observations

After estimating the parameters κ₀, μ, and α via Bayesian inference, we now evaluate the performance of the calibrated model against the full observational dataset, using the settings given in Table 1. We force the model with the observed 16 d time series of solar irradiance, wind speed, and air temperature, interpolated in time to match the CFL condition (Eq. 10). As in the previous section, we set the foundation temperature to the observed mean T_f=298.19 K and adjust the air temperature time series accordingly.

Figure 5 illustrates the results of this simulation in comparison with the observed diurnal warming. Similar to the observations, the model produces a wide range of diurnal warming amplitudes, from 0.2 K on day 9 to almost 4 K on day 13, which coincides with the maximal observed diurnal warming of 4.9 K (Fig. 5a). Modeled diurnal warming peaks align in time with the observations. On days 6 and 7, the model overestimates ΔSST by more than 1 K (see Sect. 5). The variation in diurnal warming amplitudes links closely to wind speed (Fig. 7b). To some extent, the model captures skin cooling at night, where ΔSST drops below 0 K in agreement with the observations (e.g., on the nights of days 3–4, 6–7, 9–10, and 12–13; see Fig. 5a). Overall, the reduction in ΔSST due to skin cooling is less by about 0.07 K in the model than in observations.

Figure 5 Simulation results for the calibrated model, forced with the observational dataset. (a) Modeled (red) and observed (blue) time series of diurnal warming. Gray-shaded intervals indicate the training data used for calibration. The subpanel below depicts wind speed (green, smoothed) for reference. (b) Modeled sea temperature $\mathrm{\Delta }T=T-T_{\mathrm{f}}$ as a function of depth and time, shown for day 13, 00:00, to day 15, 00:00 (local sun time). The green line indicates wind speed. (c) Modeled vertical temperature profiles on day 13 (black) and 14 (orange), plotted at intervals of 2 h (shifted by 0.5 K h⁻¹). The scale bar indicates a temperature difference of 1 K along the x axis.

For the whole time series, the modeled and observed ΔSSTs are correlated with a Pearson correlation coefficient of 0.74 (Fig. 6b). Considering only the training data, the correlation coefficient is 0.82, whereas the value for all data points outside of the training data is 0.67. This shows that our model has predictive skill in situations that the model has not been calibrated to. In particular, the model predicts both the absence of diurnal warming on day 9 and the strong consecutive warming event on day 10 (Fig. 5a). The model performance is further discussed in comparison with other models (Sects. 4 and 5).

Figure 6 Comparison of the calibrated DiuSST, ZB05, and slab models. (a) Time series of diurnal warming ΔSST as modeled by DiuSST (red), ZB05 (dark blue), and the slab (orange) for a select time interval. Points (light blue) indicate the observations. (b) Net surface heat loss to the atmosphere, $Q_{\text{net}}=R_{\text{lw}}+Q_{\text{l}}+Q_{\text{s}}$ (excluding shortwave radiation), for the model (red) and slab (orange). (c) Difference in Q_net between our model and the slab (black), as well as our model and a fixed SST at T_f (gray). (d) Correlations between the modeled and observed ΔSST for our model (top), as well as ZB05 and the slab (bottom). The Pearson correlation coefficients (corr.) are given. In the top panel, circles (crosses) mark training (validation) data points.

3.5 Depth-resolved temperature profile

Rather than simulating diurnal warming only at the surface, our model provides the vertical temperature profile within the upper 10 m of the ocean (Fig. 5b, c). In agreement with observations, heat trapping under calm and clear conditions occurs in the uppermost meter (Soloviev and Lukas, 1997; Ward, 2006; Gentemann et al., 2009). Following peak warming, the model exhibits a deepening of the warm layer between noon and sunset, as seen in Price et al. (1986) and Soloviev and Lukas (1997).

Analyzing the temporal evolution of the temperature profile on days 13 and 14 reveals the close connection to the wind (Fig. 5b), with increased winds causing fast diffusion of the warm layer on day 13 after noon and overall weaker temperature gradients on day 14. Nighttime skin cooling is visible as a slightly negative near-surface temperature gradient between around 22:00 and 06:00 local sun time (Fig. 5c). The fact that the model produces qualitatively realistic temperature profiles, even though its parameters were merely calibrated with respect to the temperature difference between the surface and a reference depth, supports the physical basis of the conceptual model.

4 Comparison with other models

Our model (hereafter referred to as DiuSST) is intended for use as an interactive SST boundary in idealized atmospheric simulations. This calls for a comparison with existing models that could be selected for this purpose. First, we consider a slab ocean model of the type previously used in atmospheric convection studies (Hohenegger and Stevens, 2016; Shamekh et al., 2020a; Coppin and Bony, 2017; Tompkins and Semie, 2021). Second, we test the prognostic scheme by Zeng and Beljaars (2005), which has been widely applied in weather and climate modeling but, to the best of our knowledge, has not been employed by the idealized convection community. Finally, we compare all models against the observations of the MOCE-5 dataset and discuss their effects on air–sea heat exchange.

4.1 Slab model

We use a single-layer slab with temperature dynamics described by

$$\begin{array}{}\text{(15)}& \dot{T}\left(t\right)={\displaystyle \frac{Q_{\text{0}}\left(t\right)-S}{{\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{p}}h}}-{\mathit{\xi }}_{\mathrm{1}}T\left(t\right)-{\mathit{\xi }}_{\mathrm{2}}\underset{\mathrm{0}}{\overset{t}{\int }}T\left(t^{\prime }\right)\mathrm{d}{t}^{\prime },\end{array}$$

where T denotes the slab temperature relative to T_f, h is the slab thickness, S represents a constant heat sink, and the net surface heat flux Q₀ is given by Eq. (4). The constants ρ_w, c_p, and T_f are listed in Table 1. In addition to the heat sink S, we include two correction terms that are sometimes added to control the slab temperature: a linear relaxation and an integral correction to prevent temperature drift. Their strengths are tuned via the parameters ξ₁ and ξ₂, respectively. Here t≥0 and T(0)=0.

To compare this slab model with our model, we calibrate the parameter set ${\mathrm{\Theta }}_{\text{slab}}=\mathit{\{}h,S,{\mathit{\xi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{2}}\mathit{\}}$ using Bayesian inference, taking the same data and settings as when calibrating our model (see Sect. 3.3). In the case of the slab, diurnal warming ΔSST directly corresponds to the slab temperature anomaly T. The resulting parameter estimates are given in Table E1. The Bayesian MCMC sampling converges to $S=\mathrm{92}{\text{W\hspace{0.17em}m}}^{-\mathrm{2}}$, a realistic approximation of net oceanic heat uptake in the tropics. Simulation results are discussed together with the other models (Sect. 4.3).

4.2 ZB05 model

The ZB05 model by Zeng and Beljaars (2005) consists of two components: a bulk equation describing the temperature evolution in the diurnal warm layer of depth d and a skin layer equation representing the cool-skin effect, that is, cooling due to surface heat fluxes within the upper millimeter of the ocean (Wong and Minnett, 2018; Fairall et al., 1996). The sea skin temperature relative to the foundation temperature T_d is thus given by the sum of warm-layer heating ΔT and skin cooling δT,

$$\mathrm{\Delta }\text{SST}\left(t\right)=\mathrm{\Delta }T\left(t\right)+\mathit{\delta }T\left(t\right).$$

In order to integrate over the warm layer, the model assumes a power law temperature profile parameterized by an empirical shape parameter ν=0.3. The integrated equation reads as follows (Eq. 11 in Zeng and Beljaars, 2005):

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\mathrm{d}\left(\mathrm{\Delta }T\right)}{\mathrm{d}t}}={\displaystyle \frac{Q_{\mathrm{0}}-R(-d)}{d{\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{p}}\mathit{\nu }/(\mathit{\nu }+\mathrm{1})}}-{\displaystyle \frac{(\mathit{\nu }+\mathrm{1})k{u}_{*}\mathrm{\Delta }T}{d{\mathit{\varphi }}_t(d/L)}},\end{array}$$

where $\mathrm{\Delta }T=T_{-\mathit{\delta }}-T_d$ denotes the temperature difference between the skin layer depth δ and the reference depth d=3 m, Q₀ is the net surface heat flux (Eq. 4), R(−d) is the downward-penetrating shortwave radiation at depth d, and k=0.4 is the von Kármán constant. Furthermore, u_* represents the friction velocity, which we compute in terms of the wind speed u and drag coefficient $C_{\mathrm{D}}=\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{3}}$ according to $u_{*}=\sqrt{{\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{D}}{u}^{\mathrm{2}}/{\mathit{\rho }}_{\mathrm{w}}}$ (Trenberth et al., 1989; Kara et al., 2007). The stability function ϕ_t is derived from Monin–Obukhov similarity theory and depends on the Monin–Obukhov length L. Here we use the refined formulation of Takaya et al. (2010) for ϕ_t.

In describing the cool-skin layer, we follow the scheme by Fairall et al. (1996, 2003) as implemented by the National Oceanic and Atmospheric Administration (NOAA) Physical Sciences NOAA Physical Sciences Laboratory (2023). The evolutions of the skin layer thickness δ and skin layer temperature difference T_d are computed based on each other in an iterative fashion. The reason why we do not directly use Eq. (6) in Zeng and Beljaars (2005) to determine δ is that this equation differs from the referenced version in Fairall et al. (1996) and leads to numerical instabilities.

4.3 Comparison of model performance

Overall, we find that both our DiuSST model and the ZB05 model give similar, adequate performance when tested on the MOCE-5 dataset (Figs. 6, 7, 8). On the contrary, the slab model fails to describe the observed diurnal SST variability even though it has been optimized on the training set of these data.

Figure 7 Diurnal warming in the observations and simulations with the DiuSST, ZB05, and slab models (from left to right), showing all data points of the 16 d time series. (a) ΔSST as a function of local sun time, colored by wind speed. (b) ΔSST as a function of wind speed, colored by local sun time.

Figure 8 Wind dependence of the amplitude of diurnal warming, comparing the observations against the predictions by the DiuSST, ZB05, and slab models (see figure legend). The data points represent the mean diurnal warming amplitude (see main text). Solid lines indicate least-squares exponential fits, with parameters given in Table 2. Wind speeds above 6.5 m s⁻¹ are not shown due to a small number of raw data points.

Specifically, the slab model is unable to capture the wind-dependent variation in diurnal warming amplitudes (Fig. 6a). The slab simulation exhibits a rather invariant diurnal amplitude of just above ΔSST=1 K, largely independent of wind speed (Fig. 7b). The poor agreement between the slab and the observations manifests itself in the relatively low Pearson correlation coefficient of 0.50 and a standard mean error of 0.36 K (over the whole dataset). In comparison, the DiuSST and ZB05 models yield a Pearson correlation of 0.74 and 0.72, respectively, each deviating from the observations by around 0.29 K on average. Even more than DiuSST, the ZB05 model overpredicts diurnal warming on days 6 and 7 (Fig. B1). ZB05 exhibits more high-frequency variability during diurnal warming events compared to DiuSST, which we attribute to the more sensitive skin layer model explicitly included in ZB05. Nighttime skin cooling is overall stronger in ZB05 relative to DiuSST, agreeing better with observations during the second half of the time series (days 9–16) but exaggerating skin cooling during the first, windier half (days 1–7; see Fig. B1). While capturing some of the observed nighttime dynamics of ΔSST, the slab tends to produce excess skin cooling, especially during strong winds.

Figure 7b highlights the improved wind dependence of ΔSST in DiuSST compared to the slab model. Similar to observations and ZB05, our model exhibits a roughly exponential decay of peak diurnal warming with increasing wind speed. To provide a more quantitative analysis, we bin each diurnal warming time series by wind speed, selecting a bin size of 0.5 m s⁻¹. For each bin, we calculate hourly averages of ΔSST as a function of local sun time and take the maximum of these hourly averages as an estimate of the mean diurnal warming amplitude for the given wind speed (Fig. 8). Based on exponential fits, ΔSST as a function of wind speed decays with a scaling constant of roughly 2 m s⁻¹ for the observations and the DiuSST and ZB05 models (Table 2). In the case of the slab, the exponent is almost 30 m s⁻¹, confirming the weak wind sensitivity of the diurnal amplitude. The expected diurnal warming amplitude under calm conditions, approximated by the intercept of the fit at u=0, is underestimated (overestimated) by DiuSST (ZB05) by about 0.5 K each. The slab diurnal warming amplitude is significantly too low.

Table 2 Exponential fits of the wind dependence of diurnal warming (see Fig. 8).

The fit function is $\mathrm{\Delta }\text{SST}\left(u\right)=y_{\mathrm{0}}\exp (-u/a)$

4.4 Air–sea heat fluxes

The diurnal evolution of SST impacts the atmosphere by regulating the heat, moisture, and momentum transfer at the air–sea interface. Over the course of the observational data interval (1–16 October), we compute the net surface heat flux $Q_{\text{net}}\left(t\right)=R_{\text{lw}}\left(t\right)+Q_{\text{l}}\left(t\right)+Q_{\text{s}}\left(t\right)$ (excluding shortwave radiation) for each model simulation (Fig. 6b; Q_net<0 corresponds to net heat export from ocean to atmosphere). In addition, we compute Q_net(t) for the given atmospheric forcing under the assumption that SST is fixed at $T(\mathrm{0},t)=T_{\mathrm{f}}$ for all t.

The oceanic heat loss to the atmosphere can differ by up to around 50 W m⁻² between DiuSST and the slab (Fig. 6c). Q_net is more negative for DiuSST when it produces a larger sea skin temperature compared to the slab and vice versa. During strong diurnal warming events, the difference ΔQ_net is even larger when comparing DiuSST to a fixed SST, exceeding 50 W m⁻² on 13 October. We conclude that using a fixed SST as an oceanic boundary condition leads to biases in surface heat fluxes, particularly during strong diurnal warming events, whereas using a slab model causes biases under strong insolation in both windy and calm conditions. Between DiuSST and ZB05, differences in the net surface heat flux are relatively small and mainly due to differences in the magnitude of skin cooling.

5 Discussion

This study aims to contribute a diurnal warm-layer model that offers a simple improvement over the slab ocean models previously used in idealized atmosphere–ocean studies. Upper-ocean heat transfer involves complex processes, from wave breaking and Langmuir circulations to biological productivity (Edson et al., 2007; Noh et al., 2004). To the first order, however, the diurnal variability in SST is governed by the competing effects of solar absorption and wind-driven turbulent mixing. This permits our reductionist modeling approach, which approximates the temperature dynamics in the diurnal warm layer mainly as a function of wind speed and insolation. For detailed realism of diurnal warm-layer dynamics, more comprehensive models (as listed in the introduction) would be more appropriate.

5.1 Model limitations and extensions

Due to its one-dimensional setup, DiuSST neglects horizontal flows and heat exchange. Furthermore, turbulence is dissipated instantaneously; i.e., there is no memory of momentum due to the wind stress history. These restrictions could explain the poor performance of DiuSST (and ZB05, likewise) on days 6 and 7, when the amplitude of diurnal warming is overestimated by more than 1 K (Fig. 5a). Despite strong insolation and low winds before noon on these days, observed diurnal warming does not exceed 1 K. This suggests the presence of enhanced turbulent vertical mixing not captured by the models, possibly due to non-local effects such as horizontal currents or rough seas, either advected from a windy region or remaining from a windy episode preceding the observations.

For simplicity, DiuSST represents stratification in the diurnal warm layer only in a time-averaged sense via the diffusivity profile φ(z), which suppresses turbulent diffusion near the surface. In reality, upper-ocean water column stability depends on the interplay between time- and depth-dependent density gradients and turbulence dissipation (Hughes et al., 2020). Instead of the highly simplified linear profile given in Eq. (3), a nonlinear, dynamic diffusivity profile could incorporate the nature of upper-ocean turbulence more realistically. For instance, a state-dependent profile, $\mathit{\phi }=\mathit{\phi }(z,T(z,t\left)\right)$, could reflect the temperature dependence of stratification. We discuss alternative diffusivity profiles in Appendix D.

The mixing term controlled by μ ensures that the diurnal layer temperature relaxes back to the foundation temperature, resetting the temperature profile at night. Combined with the diffusion term, it can only serve as a crude account of the complex mixing processes in the upper ocean (Hughes et al., 2020), such as convective overturning or internal waves. However, the contribution from this term is typically small compared to the other terms, particularly near the surface. In some cases, it may be realistic to further decrease μ to allow for multiday warming events (Jia et al., 2023).

With a vertical resolution set to 10 cm at the surface, our numerical implementation is too coarse to explicitly resolve the cool-skin layer, which forms within millimeters of the air–sea interface. Nonetheless, the model still captures skin cooling at night, indicated by slightly negative values of ΔSST, in agreement with the observations (Fig. 5). In fact, the model considers the cool-skin effect in a coarse-grained sense, averaging heat fluxes over the upper 10 cm of the water column. Of course, the grid resolution could be increased, albeit at the cost of a smaller integration time step to maintain numerical stability.

For future refinement, further dynamical processes could be built into DiuSST via the three key modular terms. A time-dependent diffusivity profile would allow us to parameterize the effect of precipitation, which we presently neglect in the model. Rain freshwater pools can enhance either heat trapping or mixing, depending on the competing effects of salinity and temperature (Webster et al., 1996; Bellenger et al., 2017; Witte et al., 2023). As another example, a dynamic attenuation coefficient α(t) would allow us to account for changes in the seawater optical properties, e.g., due to microbial activity (Wurl et al., 2017).

5.2 Data-driven parameter estimation

While the parameters $\mathit{\{}{\mathit{\kappa }}_{\mathrm{0}},\mathit{\alpha },\mathit{\mu }\mathit{\}}$ are fixed in our model, they generally depend on oceanic conditions that may vary in space and time. A major benefit of Bayesian inference is that the model parameters can be readily re-calibrated to additional data of interest.

The MOCE-5 cruise data used here cover measurements across approximately 10° latitude and longitude, in both the open Pacific Ocean and the coastal waters of the Gulf of California, where the dynamical and optical properties of seawater are likely to have differed. The observed foundation temperature varied by several Kelvin during the cruise, whereas in the model we fix T_f to the observed mean. Thus, spatial heterogeneity of water properties probably constitutes a main error source between the model and observations. For example, on day 13, the modeled maximum of diurnal warming is about 1 K below the observed value of ΔSST≈4.9 K. This difference may be attributed to spatial variations in optical water properties affecting the attenuation coefficient α. The cruise vessel's location on day 13 near the Midriff Islands is known for high phytoplankton concentrations, which enhance the absorption of shortwave radiation (Álvarez Borrego, 2012). Indeed, we can accurately model the maximum of ΔSST on day 13 by increasing α, but this reduces the model performance across the time series overall.

At the same time, the heterogeneity of the training data presents a benefit, as it reduces the risk of overfitting and promises a more generally applicable calibration. Even with constant parameters obtained from Bayesian inference, our model captures the observed variability in ΔSST despite the heterogeneity of the MOCE-5 dataset. It was important to train on a period including both extremes, i.e., small and large diurnal warming events, since otherwise the model performance was significantly poorer (see Börner, 2021, Sect. 6.1). Since all three parameters influence the amplitude of diurnal warming, a certain range of data is necessary to disentangle their effects (for instance, is it warming due to high insolation or low diffusivity?).

For idealized atmospheric modeling, we argue that one is usually interested in average sea properties for conditions of interest or for a parameter sensitivity experiment. While the present calibration could be specific to the background conditions of the study region, the model formulation itself is not. The parameter posterior distribution reported here can be used directly as a prior distribution for a re-calibration based on relevant additional data.

To validate the vertical profiles generated by DiuSST, depth-resolved observations are needed. Incorporating information at depth in the calibration can further constrain the parameter estimates, help resolve apparent parameter correlations, and inform the choice of diffusivity profile φ. However, the vertical temperature profiles measured during the MOCE-5 cruise (Ward, 2006) are subject to various other influences besides diurnal warming (Gentemann et al., 2009) and are currently not publicly available. Validating the profiles is thus left for future work.

5.3 Inaptitude of the slab ocean

In contrast to the DiuSST and ZB05 models, the slab model fails to capture the observed diurnal warming dynamics, including the wind dependence of ΔSST. This is because the evolution of slab temperature depends on wind speed only via the latent and sensible surface heat fluxes, whereas wind strongly influences vertical heat transport via the diffusion term in DiuSST or via the stability function in ZB05. In other words, the wind dependence of upper-ocean turbulence controls the effective heat capacity of the upper ocean, but the slab's heat capacity is fixed by the slab thickness h. This issue also pertains to ocean models with a surface layer of 1 m thickness or more, which may reproduce a mean diurnal cycle but not the full range of diurnal variability (Voldoire et al., 2022).

One might argue that the slab exhibits too little variability because the parameter ξ₁, the inverse timescale of temperature relaxation, is set too large. Indeed, decreasing ξ₁ allows for stronger diurnal warming but also leads to excessive nighttime cooling and a slow temperature decline in the afternoon, thus worsening the agreement with the observations overall (see Fig. F1). In fact, the shape of the diurnal warming curve produced by a slab with small corrector coefficients ξ₁ and ξ₂ resembles the diurnal temperature evolution observed at around 1 m depth, e.g., by moored buoys (Börner, 2021). This highlights the fact that slab oceans mimic bulk SST rather than skin SST. Also, other models of the diurnal SST cycle are sometimes calibrated with respect to bulk SST. However, skin SST is the relevant quantity for atmospheric studies since the atmosphere only senses the temperature of the sea skin.

5.4 Interactive diurnal SST in atmospheric simulations

As a more realistic alternative to a slab model, the DiuSST and ZB05 models can each be coupled to high-resolution cloud-resolving atmospheric models, e.g., to study how air–sea interactions impact atmospheric convection in idealized setups.

The appeal of the ZB05 model lies in its rigorous derivation from Monin–Obukhov similarity theory condensed into one integrated bulk equation for the warm layer, combined with a skin-layer scheme. This makes it computationally fast, although calculating the skin layer depth δ involves an implicit equation that requires iterative solving (see Eqs. 5 and 6 in Zeng and Beljaars, 2005). The empirical shape parameter ν of the vertical temperature profile could serve as a tuning parameter in sensitivity experiments. While ZB05 matches fairly well with the observations used in this study, its performance has varied in other comparative studies (Bellenger and Duvel, 2009; Jia et al., 2023). Running DiuSST requires integrating a discretized PDE, making it slightly slower than ZB05 in the Python scripts implemented here (Börner, 2024) but not significantly – especially relative to the cost of a coupled atmospheric component. Compared to ZB05 or similar existing diurnal warm-layer schemes, modelers may find DiuSST easier to tune, interpret, and adapt to their purposes due to the conceptually simple modular structure. The fact that DiuSST resolves the vertical temperature profile opens up opportunities to study, for example, how biological upper-ocean processes interact with weather and climate across the air–sea interface.

Recently, we implemented our upper-ocean model as an interactive sea surface boundary condition in the System for Atmospheric Modeling (SAM) (Khairoutdinov and Randall, 2003). At each horizontal grid point of the atmospheric model, the surface boundary condition is independently updated at each time step by numerically integrating Eq. (A1) to the next time step based on the local atmospheric forcing ℱ. This produces a responsive, spatially heterogeneous sea surface whose temperature feeds back into the atmospheric boundary layer. According to first tests, the difference in computation time between this setup and a coupled slab ocean is unnoticeable.

The influence of our model on the spatiotemporal patterns of convection is subject of future research. As a global kilometer-scale simulation study by Shevchenko et al. (2023) suggests, diurnal warm layers can have a significant impact on the atmosphere locally where diurnal warming is strong, yet their role in the large-scale tropical cloud field and convective organization remains unclear.

6 Conclusions

This paper presents a simple, one-dimensional prognostic model of diurnal sea surface temperature variability in the tropical ocean, described by Eq. (1) and formulated as a numerical scheme in Eq. (A1). Written as a single partial differential equation, the model describes upper-ocean heat transfer through three idealized terms controlled by three tuning parameters: an eddy diffusivity κ₀, a bulk mixing rate μ, and an attenuation coefficient α. κ₀ controls the strength of wind-driven turbulent heat transport, μ determines the relaxation rate towards the foundation temperature, and α specifies how deeply solar radiation penetrates into the ocean.

First, we used Bayesian inference to estimate the values of these parameters based on an observational dataset recorded on the MOCE-5 cruise in the eastern Pacific. Then, we compared the performance of our model with two other models of diurnal SST dynamics: a slab ocean model, as previously used to mimic a responsive sea surface in atmospheric simulations, and the diurnal warm-layer scheme ZB05 (Zeng and Beljaars, 2005).

Our results show that slab models with fixed heat capacity cannot capture diurnal SST variability realistically. Instead, our model reproduces an exponential dependence of the diurnal warming amplitude on wind speed, in accordance with observations and similar to ZB05. By introducing a diffusion term that scales with surface wind stress, we propose a simple solution that offers significantly improved results compared to a slab model and parameterizes the basic features of upper-ocean turbulence in a physically interpretable way. Enhancing the model in the future, e.g., by refining the diffusivity profile to include effects of precipitation, is facilitated by the concise modular code. To confirm the model's validity for diverse conditions, the Bayesian approach offers a natural way to update the parameter estimates based on additional data.

Due to its numerical and conceptual simplicity, we envision the model presented in this study serving as a generic interactive boundary condition for oceans in idealized cloud-resolving simulations of the atmosphere. This will enable a wind-responsive SST field that evolves under and feeds back into the atmospheric fluid dynamics. Given the multiscale interaction between the diurnal cycle and large-scale patterns of tropical convection, it becomes increasingly clear that we must consider diurnal warm-layer dynamics to better understand the mechanisms of marine cloud organization. We thus hope with our work to contribute to bridging the gap between idealized studies of convective aggregation and real-world process understanding.

Appendix A: DiuSST model – discretized model equation

The discretized form of Eq. (1), using an explicit Euler scheme in time and finite differences on a non-uniform grid in space, reads

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& {\displaystyle \frac{T_n^{i+\mathrm{1}}-T_n^i}{\mathrm{\Delta }{t}_i}}=\mathit{\kappa }(z_n,t_i)\left[(T_{n+\mathrm{1}}^i-\mathrm{2}{T}_n^i+T_{n-\mathrm{1}}^i)\cdot {\left({\left.{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}z}}\right|}_{z_n}\right)}^{\mathrm{2}}\right.\\ & \left.+{\displaystyle \frac{T_{n+\mathrm{1}}^i-T_{n-\mathrm{1}}^i}{\mathrm{2}}}\cdot {\left.{\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}n}{\mathrm{d}{z}^{\mathrm{2}}}}\right|}_{z_n}\right]\\ & +{\left.{\displaystyle \frac{\partial \mathit{\kappa }(z,t_i)}{\partial z}}\right|}_{z_n}\left[{\displaystyle \frac{T_{n+\mathrm{1}}^i-T_{n-\mathrm{1}}^i}{\mathrm{2}}}\cdot {\left.{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}z}}\right|}_{z_n}\right]-\mathit{\mu }{\displaystyle \frac{T_n^i-T_{\mathrm{f}}}{|z_n-z_{\mathrm{f}}|}}\\ & +{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}{c}_{\mathrm{p}}}}\left[(Q_{n+\mathrm{1}}^i-Q_n^i)\cdot {\left.{\displaystyle \frac{\mathrm{d}n}{\mathrm{d}z}}\right|}_{z_n}\right],\\ & (\text{for }n=\mathrm{0},\mathrm{\dots },N-\mathrm{1}),\end{array}\end{array}$$

where n is the depth index and i the time index (nomenclature as in the main text; see also Table 1). Here κ(z,t) is given by Eq. (2); its derivative by z is computed analytically. The derivatives $\frac{\mathrm{d}n}{\mathrm{d}z}$ and $\frac{{\mathrm{d}}^{\mathrm{2}}n}{\mathrm{d}{z}^{\mathrm{2}}}$,

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}n}{\mathrm{d}z}}={\left[\ln \mathit{ϵ}\cdot \left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{0}}}{\mathrm{1}-\mathit{ϵ}}}+z\right)\right]}^{-\mathrm{1}},\\ & {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}n}{\mathrm{d}{z}^{\mathrm{2}}}}={\left[-\ln \mathit{ϵ}\cdot {\left({\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{0}}}{\mathrm{1}-\mathit{ϵ}}}+z\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}},\end{array}\end{array}$$

map the finite differences onto the non-uniform grid spacing, where the stretch factor ϵ>1 solves the equation $\mathit{ϵ}(\mathrm{\Delta }{z}_{\mathrm{0}},N)=\mathrm{1}+\left(\mathrm{1}-{\mathit{ϵ}}^N\right)\mathrm{\Delta }{z}_{\mathrm{0}}/z_{\mathrm{f}}$ (see also Eq. 9). At the foundation depth, the Dirichlet boundary condition is $T_N^i=T_{\mathrm{f}}$ for all i. At the surface (n=0), we require a choice when dealing with the free boundary. We introduce a dummy grid point ($n=-\mathrm{1}$) at $z_{-\mathrm{1}}=\mathrm{\Delta }{z}_{\mathrm{0}}$, at which we set the temperature to $T_{-\mathrm{1}}^i=T_{\mathrm{0}}^i$ for all i (essentially assuming that the temperature gradient vanishes at the air–sea interface). The surface heat flux (at n=0) is given by

$$\begin{array}{}\text{(A3)}& Q_{\mathrm{0}}^i:=R_{\text{sw}}\left(t_i\right)+R_{\text{lw}}\left(t_i\right)+Q_{\text{s}}\left(t_i\right)+Q_{\text{l}}\left(t_i\right).\end{array}$$

For n>0, the heat flux $Q_n^i\equiv Q(z_n,t_i)$ is given by Eq. (7). Using a forward difference for $Q_n^i$ ensures that the integrated heat flux over the domain corresponds to the total heat uptake within the diurnal layer. Note that the diffusion part of Eq. (1) involves two terms in Eq. (A1), owing to the chain rule of differentiation.

Appendix B: Model comparison over the full MOCE5 time series

Figure B1 Comparison of observed diurnal warming ΔSST (light blue) against the predictions of the calibrated DiuSST (dashed red), ZB05 (dark blue), and calibrated slab (dash-dotted orange) models, showing the full 16 d time series recorded on the MOCE-5 cruise.

Simulation results over the whole time interval of the MOCE-5 observations are shown in Fig. B1 for all models. DiuSST and ZB05 produce similar results, with ZB05 occasionally overshooting the observations by more than 1 K around diurnal warming peaks, e.g., on days 6, 10, and 14. DiuSST produces a more accurate nighttime skin temperature throughout the windy first half of the time series (days 2–8), whereas ZB05 yields better nighttime skin temperatures during the second half of the time series. The slab produces a constant diurnal warming amplitude throughout, overestimating ΔSST on windy days and underestimating it on calm days.

Appendix C: Surface reflection of solar radiation in DiuSST

Based on Fresnel's equations, assuming unpolarized light, the reflected fraction ℛ of irradiance incident on the (flat and smooth) air–sea interface is given by $\mathcal{R}=({\mathcal{R}}_{⟂}+{\mathcal{R}}_{\parallel })/\mathrm{2}$, where the contributions from the two polarization directions read

$$\begin{array}{}\text{(C1)}& {\displaystyle }{\mathcal{R}}_{⟂}={\left({\displaystyle \frac{n_{\mathrm{a}}\cos \mathit{\varphi }-n_{\mathrm{w}}\sqrt{\mathrm{1}-(\frac{n_{\mathrm{a}}}{n_{\mathrm{w}}}\sin \mathit{\varphi }{)}^{\mathrm{2}}}}{n_{\mathrm{a}}\cos \mathit{\varphi }+n_{\mathrm{w}}\sqrt{\mathrm{1}-(\frac{n_{\mathrm{a}}}{n_{\mathrm{w}}}\sin \mathit{\varphi }{)}^{\mathrm{2}}}}}\right)}^{\mathrm{2}}\text{(C2)}& {\displaystyle }{\mathcal{R}}_{\parallel }={\left({\displaystyle \frac{n_{\mathrm{a}}\sqrt{\mathrm{1}-(\frac{n_{\mathrm{a}}}{n_{\mathrm{w}}}\sin \mathit{\varphi }{)}^{\mathrm{2}}}-n_{\mathrm{w}}\cos \mathit{\varphi }}{n_{\mathrm{a}}\sqrt{\mathrm{1}-(\frac{n_{\mathrm{a}}}{n_{\mathrm{w}}}\sin \mathit{\varphi }{)}^{\mathrm{2}}}+n_{\mathrm{w}}\cos \mathit{\varphi }}}\right)}^{\mathrm{2}}.\end{array}$$

Here ϕ is the solar angle with respect to the surface normal; n_a and n_w denote the refractive indices of air and water, respectively (see Table 1). Thus, we approximate the transmitted solar irradiance R_sw entering the water body as

$$\begin{array}{}\text{(C3)}& R_{\text{sw}}\left(t\right)=(\mathrm{1}-\mathcal{R}(\mathit{\varphi }\left(t\right)\left)R_{\text{sw},\downarrow }\right(t),\end{array}$$

where $R_{\text{sw},\downarrow }$ is the downward shortwave irradiance above the sea surface.

When neglecting surface reflection and the effect of surface refraction on the absorption of downward shortwave radiation, the simulated ΔSST (under the same parameter settings) can be larger by up to 0.2 K for the MOCE-5 case study (Fig. C1). However, note that if we would have calibrated the model without reflection and refraction, this difference would likely be compensated by adjusting the parameters.

Figure C1 Difference in simulated ΔSST of the DiuSST model without reflection and refraction vs. the original DiuSST model, for the parameters given in Table 1 and forcing with MOCE-5 data.

Appendix D: Alternative diffusivity profiles

Figure D1 Comparison of modeled diurnal warming ΔSST between different diffusivity profiles φ discussed in Appendix D and the main text. LIN (solid red) corresponds to the linear profile (Eq. 3) used throughout the main text. EXP (dash-dotted blue) and STAB (dashed orange) show results for the corresponding versions of φ described in Appendix D. For each case, the model parameters κ₀, μ, and α are calibrated specifically via Bayesian inference. All other model settings are kept constant, as given in Table 1.

In the main text, we introduce a highly idealized linear diffusivity profile φ(z) (Eq. 3, called LIN hereafter). Here we briefly discuss two alternatives that have been investigated in the course of this study. Consider the profile

$$\begin{array}{}\text{(D1)}& \mathit{\phi }(z,t)={\displaystyle \frac{\mathrm{1}-S\left(t\right)\mathit{\sigma }\exp (z/\mathit{\lambda })}{\mathrm{1}-S\left(t\right)\mathit{\sigma }\exp (z_{\mathrm{f}}/\mathit{\lambda })}},\end{array}$$

given in terms of the suppressivity $\mathit{\sigma }\in [\mathrm{0},\mathrm{1}]$, the trapping depth λ>0, and the (time-dependent) stability function S.

First, set S(t)=1 for all t, giving a time-independent exponential profile (EXP). If σ>0, then φ decreases exponentially towards the sea surface, where $\mathit{\phi }\left(\mathrm{0}\right)=\mathrm{1}-\mathit{\sigma }$. The parameter λ, setting the curvature of the profile, determines the depth scale at which the diffusivity is suppressed. In the limit λ→∞, we obtain the linear profile in Eq. (3).

Figure D2 Comparison of vertical temperature profiles produced by the model with diffusivity profiles EXP (dark blue) and STAB (orange). Panel (a) shows the vertical profile every 2 h, shifted by 1 K h⁻¹ along the horizontal axis. The bottom panels show the temperature difference ΔT with respect to the foundation temperature as a function of depth and time for the EXP model (b) and STAB model (c).

Second, we investigate a state-dependent stability function (STAB) that depends on the vertical temperature gradient $\mathrm{\Delta }{T}_{\mathit{\lambda }}\left(t\right):=T(\mathrm{0},t)-T(\mathit{\lambda },t)$, given by

$$\begin{array}{}\text{(D2)}& S\left(t\right)=\left\{\begin{array}{ll}min(\mathrm{1},\mathrm{\Delta }{T}_{\mathit{\lambda }}(t)/T_{\text{strat}})& \text{ if }\mathrm{\Delta }{T}_{\mathit{\lambda }}\left(t\right)>\mathrm{0}\\ \mathrm{0}& \text{otherwise,}\end{array}\right.\end{array}$$

where T_strat is a reference temperature at which thermal stratification dominates over turbulent mixing. The physical rationale behind this is that density stratification due to near-surface heat trapping locally enhances water column stability, inhibiting vertical diffusion. As soon as the sea skin temperature cools again with respect to the temperature at depth λ, diffusivity increases, mimicking enhanced afternoon mixing due to an unstable temperature profile.

Comparing simulations with the LIN and EXP profiles, we find that EXP performs slightly better than LIN for the observational dataset considered in this study, particularly with regard to nighttime skin cooling (Fig. D1). Although arguably more realistic than LIN and EXP, the STAB profile does not necessarily perform better because it kills the coarse-grained skin cooling effect in our model. Comparing the effects of EXP and STAB on the vertical temperature profile (Fig. D2), we see that STAB counteracts near-surface temperature inversions, creating a mixed warm layer of near-constant temperature that better matches observations (Soloviev and Lukas, 1997). However, STAB would have to be combined with a skin layer scheme as in ZB05 in order to capture the sea skin temperature in the presence of skin cooling. Due to the overall similarity of modeled SST warming between the LIN, EXP, and STAB versions, we choose to present the simplest option, LIN, in the main text.

Appendix E: Bayesian inference

Consider a model ℳ(Θ) controlled by the parameter set Θ. Given the data 𝒟, the probability that a certain value of Θ represents the truth follows from Bayes' rule,

$$\begin{array}{}\text{(E1)}& P\left(\mathrm{\Theta }\right|\mathcal{D})={\displaystyle \frac{P\left(\mathcal{D}\right|\mathrm{\Theta }\left)P\right(\mathrm{\Theta })}{P\left(\mathcal{D}\right)}}.\end{array}$$

Here P(A|B) denotes the conditional probability of A given B. On the right-hand side, P(𝒟|Θ) is the likelihood of observing the data under the assumption that Θ represents the true model. The prior probability P(Θ) quantifies our knowledge of the parameter set before observing the data. Lastly, the denominator states the probability of the data being true, which is independent of the choice of model and merely adds a normalization factor. Hence the posterior probability P(Θ|𝒟) is proportional to the product of the likelihood times the prior probability.

In practice, inferring the posterior distribution P(Θ|𝒟) of the parameters Θ from the data 𝒟 involves three steps: (1) construct the prior distribution P(Θ) in parameter space based on previous knowledge or belief, (2) define a suitable likelihood function $\mathcal{L}\left(\mathrm{\Theta }\right)\equiv P\left(\mathcal{D}\right|\mathrm{\Theta })$, and (3) compute the (un-normalized) posterior distribution by evaluating the product ℒ(Θ)P(Θ).

E1 Choice of prior distribution

Our sensitivity study under idealized forcing (Sect. 3.1) provides an orientation of the physical ranges of the parameters $\mathit{\{}{\mathit{\kappa }}_{\mathrm{0}},\mathit{\mu },\mathit{\alpha }\mathit{\}}$ of the DiuSST model. In the literature, common values of oceanic turbulent vertical diffusivity are on the order of $\mathit{\kappa }\sim {\mathrm{10}}^{-\mathrm{4}}{\text{m}}^{\mathrm{2}}{\text{ s}}^{-\mathrm{1}}$, varying with location and depth (Denman and Gargett, 1983). Stating a physical value for the attenuation coefficient α is difficult since, in reality, attenuation of shortwave radiation depends strongly on wavelength and the biochemical composition of the seawater. Empirical values for the diffusive attenuation coefficient of photosynthetically active radiation range from $\sim {\mathrm{10}}^{-\mathrm{2}}$ to 10 m⁻¹ (Son and Wang, 2015). However, we note that the model parameters are conceptual, particularly the mixing coefficient μ, such that they do not necessarily represent directly observable physical quantities.

Reflecting our limited knowledge, we impose uniform prior distributions for the parameters κ₀ and α but constrain their range to $\mathit{\kappa }\in [\mathrm{0},\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}]$ m² s⁻¹ and $\mathit{\alpha }\in [\mathrm{0.05},\mathrm{10}]$ m⁻¹. The parameter μ requires more subtle treatment because it acts mainly near the foundation depth, whereas the likelihood function evaluates temperature differences near the surface (see below). Figure 2 indicates that excess heat remains in the interior of the diurnal layer if μ falls below $\sim {\mathrm{10}}^{-\mathrm{4}}{\text{m\hspace{0.17em}s}}^{-\mathrm{1}}$. Based on this insight, we define a normally distributed prior for μ with a mean of $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{3}}{\text{m\hspace{0.17em}s}}^{-\mathrm{1}}$ and standard deviation of $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{3}}{\text{m\hspace{0.17em}s}}^{-\mathrm{1}}$. Thus the three-dimensional prior distribution P(Θ) is uniform and bounded along the κ₀ and α directions but Gaussian with respect to μ.

Note that the posterior distribution with respect to μ differs clearly from the prior distribution, confirming that the data added information on the parameter μ (Fig. 4).

E2 Defining the likelihood function

Our training data 𝒟 consist of a 6 d subset of the diurnal warming time series from the MOCE-5 cruise, as indicated by the gray shading in Fig. 3. Specifically, we select days 2–4 and 13–15 of October 1999 as training data, while all other days from 1 to 16 October make up the validation data. By this choice, the training data contain warming events of different amplitudes, covering the observed range from a few tenths of a Kelvin up to 5 K. Moreover, the two sub-intervals composing 𝒟 correspond to different geographical regions (open Pacific vs. Gulf of California), presumably featuring differing environmental conditions. This ensures that the resulting parameter estimates will be valid for a broader range of conditions.

We choose a likelihood function that decays exponentially with the weighted square error between model and data,

$$\begin{array}{}\text{(E2)}& \mathcal{L}\propto \exp \left(-\sum _j{\displaystyle \frac{\left(\mathrm{\Delta }{\text{SST}}_{\mathcal{M}}\right(t_j)-\mathrm{\Delta }{\text{SST}}_{\mathcal{D}}(t_j){)}^{\mathrm{2}}}{{\mathrm{\Sigma }}_j^{\mathrm{2}}}}\right),\end{array}$$

where the index j runs through all data points in 𝒟, and ΔSST_ℳ (ΔSST_𝒟) denotes the modeled (observed) diurnal warming at the corresponding time of observation. Specifically, $\mathrm{\Delta }{\text{SST}}_{\mathcal{M}}\left(t_j\right)=T(\mathrm{0},t_j)-T(z_{\mathrm{d}},t_j)$, with z_d being the depth of the vertical grid point closest to the observational reference depth $d=-\mathrm{3}$ m. Furthermore, the standard error Σ_j, determining the relative weight of the jth data point, is defined as

$$\begin{array}{}\text{(E3)}& {\mathrm{\Sigma }}_j=\mathrm{2}{\mathit{ϵ}}_j(\mathrm{1}+v_j/v_{\text{max}}),\end{array}$$

where ϵ_j denotes the uncertainty in the 3 m temperature measurement, as stated in the dataset; v_j is the current boat speed; and v_max is the maximum boat speed throughout the dataset. This choice reflects a reduced confidence in measurements taken while the ship was moving rapidly. Indeed, the MOCE-5 time series features several short, isolated warming spikes of around 1–2 °C during nighttime hours, which cannot be explained by atmospheric forcing but could have been caused by the boat's movement. Since the ship mainly moved at night while often remaining at fixed locations during the day, Σ_j emphasizes daylight data.

E3 Computing the posterior distribution

To approximate the posterior distribution, we perform Markov chain Monte Carlo (MCMC) sampling using the emcee package in Python (Foreman-Mackey et al., 2013). It is based on the affine-invariant GW10 algorithm proposed by Goodman and Weare (2010). The algorithm explores the parameter space with multiple interdependent walkers whose moves efficiently deal with correlations, yielding fast convergence.

At each step of the Markov chain, the likelihood function, Eq. (E2), is evaluated for each walker at its respective position Θ in parameter space. Each evaluation requires simulating the model for the duration of the training data, which makes the sampling computationally expensive. For each of 24 walkers we generate 4000 steps, initialized at the value $\mathit{\{}{\mathit{\kappa }}_{\mathrm{0}},\mathit{\mu },\mathit{\alpha }\mathit{\}}=\mathit{\{}{\mathrm{10}}^{-\mathrm{4}}{\text{m}}^{\mathrm{2}}{\text{ s}}^{-\mathrm{1}},\mathrm{6}\times {\mathrm{10}}^{-\mathrm{3}}{\text{m\hspace{0.17em}s}}^{-\mathrm{1}},\mathrm{4}{\text{m}}^{-\mathrm{1}}\mathit{\}}$. Finding an auto-correlation length of around 50 steps (along the Markov chain), this gives around 80 independent samples.

We summarize the resulting parameter estimates in Table E1, listing common statistical estimators as obtained from the posterior distribution of the DiuSST (Fig. 4) and slab (Fig. E1) models, respectively. Bayesian inference for the slab model was performed in analogy to the DiuSST model but in a four-dimensional parameter space $\mathit{\{}h,S,{\mathit{\xi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{2}}\mathit{\}}$ and with the slab temperature anomaly representing ΔSST_ℳ in the likelihood function. The run was initialized at values $\mathit{\{}h,S,{\mathit{\xi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{2}}\mathit{\}}=\mathit{\{}\mathrm{1},\mathrm{100},{\mathrm{10}}^{-\mathrm{4}},{\mathrm{10}}^{-\mathrm{9}}\mathit{\}}$ with uniform priors for h, ξ₁, ξ₂ and a Gaussian prior for S with mean 100 W m⁻² and standard deviation 10 W m⁻².

Table E1 Parameter estimates obtained from the sampled posterior distribution.

Figure E1 Posterior distribution for the parameters of the slab model, obtained from Bayesian MCMC sampling. Compare with Fig. 4, which shows the results for the DiuSST model.

Appendix F: Slab model parameter sensitivity

Figure F1 Parameter dependence of the slab model, comparing the calibrated slab model (slab A, see Table E1) with a version (slab B) where the parameter ξ₁ is reduced by 3 orders of magnitude. For slab B, h=1.1 m, S=70 W m⁻², ${\mathit{\xi }}_{\mathrm{1}}={\mathrm{10}}^{-\mathrm{7}}$, and ξ₂=0. For reference, the observations and DiuSST simulations (here labeled model) are shown.

As mentioned in the main text, the Bayesian optimization procedure applied to the slab model acts mainly by ramping up the relaxation parameter ξ₁, which constrains diurnal warming to 1 K. Larger diurnal amplitudes are achieved by reducing ξ₁ (see Fig. F1). However, in that case, the slab also cools excessively at night, and the diurnal warming peak shifts into the late afternoon, with a slow temperature decrease after peak warming.