An empirical model for the statistics of sea surface diurnal warming

A statistical model is derived relating the diurnal variation of sea surface temperature (SST) to the net surface heat flux and surface wind speed from a numerical weather prediction (NWP) model. The model is derived using fluxes and winds from the European Centre for Medium-Range Weather Forecasting (ECMWF) NWP model and SSTs from the Spinning Enhanced Visible and Infrared Imager (SE- VIRI). In the model, diurnal warming has a linear depen- dence on the net surface heat flux integrated since (approx- imately) dawn and an inverse quadratic dependence on the maximum of the surface wind speed in the same period. The model coefficients are found by matching, for a given inte- grated heat flux, the frequency distributions of the maximum wind speed and the observed warming. Diurnal cooling, where it occurs, is modelled as proportional to the integrated heat flux divided by the heat capacity of the seasonal mixed layer. The model reproduces the statistics (mean, standard deviation, and 95-percentile) of the diurnal variation of SST seen by SEVIRI and reproduces the geographical pattern of mean warming seen by the Advanced Microwave Scanning Radiometer (AMSR-E). We use the functional dependencies in the statistical model to test the behaviour of two physical model of diurnal warming that display contrasting systematic errors. and its geophysical impacts, together with a review of exist- ing diurnal warming models, emphasizing how they differ in purpose and nature from the model developed in this paper. Section 2 describes the statistical model, with the details of data and fitting procedures given in Sect. 3. Section 4 de- scribes the performance and some validation of the model. The penultimate section of the paper gives an illustration of the use of the statistical model to assess the distributions of diurnal warming predicted by two physically based models, neither of which seem accurately to match the observed dis- tributions over the full functional range tested. We conclude the paper with a final section which puts the results presented in the paper in a wider context. Diurnal warming of the sea surface is the sub-daily varia- tion in sea surface temperature (SST) associated principally with the daily cycle in solar heating (although other influen- tial factors may also have a daily cycle). During the day, the upper few metres of the ocean are heated by short-wave solar radiation. This heating is usually partially offset by cooling via net outgoing long-wave radiation and sensible and latent heat fluxes. The top 5 m of the ocean absorbs 60 % of the in- coming solar radiation (Fairall et al., 1996), and thus there is a tendency for the near-surface to warm more than the photic zone of the water column as a whole. During the night, the water column will usually cool from the surface due to the other heat fluxes. This daily progression in heating and cool- ing gives a diurnal cycle in the sea-surface temperature. Generally, the diurnal cycle in SST is modest. Kennedy et al. (2007) observed a peak-to-peak mean amplitude in drift- ing buoy SSTs of 0.4 K for observations within 20 of lati- tude of the equator. The mean amplitude for the ocean as a whole (based on satellite observations to be discussed below)


Introduction
The purpose of this paper is to present an empirical model that properly captures the statistical distribution of diurnal warming of the sea surface. We begin, in this introduction, by summarizing the nature of sea surface diurnal warming and its geophysical impacts, together with a review of existing diurnal warming models, emphasizing how they differ in purpose and nature from the model developed in this paper. Section 2 describes the statistical model, with the details of data and fitting procedures given in Sect. 3. Section 4 describes the performance and some validation of the model. The penultimate section of the paper gives an illustration of the use of the statistical model to assess the distributions of diurnal warming predicted by two physically based models, neither of which seem accurately to match the observed distributions over the full functional range tested. We conclude the paper with a final section which puts the results presented in the paper in a wider context. Diurnal warming of the sea surface is the sub-daily variation in sea surface temperature (SST) associated principally with the daily cycle in solar heating (although other influential factors may also have a daily cycle). During the day, the upper few metres of the ocean are heated by short-wave solar radiation. This heating is usually partially offset by cooling via net outgoing long-wave radiation and sensible and latent heat fluxes. The top 5 m of the ocean absorbs 60 % of the incoming solar radiation (Fairall et al., 1996), and thus there is a tendency for the near-surface to warm more than the photic zone of the water column as a whole. During the night, the water column will usually cool from the surface due to the other heat fluxes. This daily progression in heating and cooling gives a diurnal cycle in the sea-surface temperature.
Generally, the diurnal cycle in SST is modest. Kennedy et al. (2007) observed a peak-to-peak mean amplitude in drifting buoy SSTs of 0.4 K for observations within 20 • of latitude of the equator. The mean amplitude for the ocean as a whole (based on satellite observations to be discussed below) is 0.25 K. Low amplitudes arise because surface wind stress causes mixing in the upper ocean, which mixes the energy from the Sun downwards.
However, in conditions with light winds and strong sunlight, a shallow (0.3 m to 3 m, e.g. Fairall et al., 1996, Fig. 7) and stably stratified layer can form in which the temperature increases by one kelvin or more. Indeed, the diurnal excursion of the surface temperature can exceed 3.5 K in extreme cases. Excursions exceeding 3.5 K were observed in day-night differences from polar orbiting satellites relatively early in the satellite era (e.g. Saunders et al., 1982;Stramma et al., 1986). The ability to track the hour-by-hour progression of large-amplitude events from geostationary satellites (e.g.  and to confirm them in multiple satellite data sets  has demonstrated that extreme cycles of up to 6 K in amplitude are routinely observable from space (although they have not, to our knowledge, been documented in situ in the open ocean). Diurnal variability in satellite imagery varies over horizontal scales of 10 to 1000 km (Stramma et al., 1986) with spatial variability linked to orographic effects and turbidity in coastal areas .
Prediction of the diurnal variation in SST is important for atmospheric and ocean forecasting. SST affects air-sea fluxes and convection: Haman and Clayson (2007) found in modelling studies of the tropical Pacific that the use of diurnally-varying SST rather than daily-averaged SST can change modelled mid-troposphere cloudiness. Atmosphereocean coupling including diurnal timescales is required to maintain in weather forecasts the phase relationship between SSTs and rainfall in the tropics (N. P. Klingaman, personal communication, 2009). Ocean forecasting models need SST to constrain the upper ocean thermal structure and dynamics (e.g. Stark et al., 2007). The "surface" temperature for ocean forecasting models is the temperature of their top model layer (which has typically been of order 10 m thick) and has been equated to the "foundation temperature" concept (Donlon et al., 2009). Satellite observations are sensitive to sea temperature at or near the air-sea interface, and are therefore decoupled from this foundation temperature whenever a warm-layer event occurs. Therefore, assimilation of diurnally warmed observations represents a danger of significant over-estimation of upper mixed-layer heat content. Some assimilation and analysis systems therefore discard SST data that are at risk of significant diurnal warming (Stark et al., 2007).
Aspects of climate are responsive to the diurnal cycle in SST. Bernie et al. (2005) found that the increase of the daily mean SST by the diurnal cycle of SST accounts for about one-third of the magnitude of intraseasonal variability of SST associated with the Madden-Julian oscillation in the western Pacific warm pool, and that diurnally resolved oceanatmosphere coupling improved the strength and coherency of representation of the Madden Julian Oscillation in a coupled general circulation model (Bernie et al., 2008). Both heat fluxes (Haman and Clayson, 2007) and gas exchange (e.g. Kettle et al., 2009) are non-linearly responsive to the air-sea interface temperature, and therefore may be affected in the mean by any significant diurnal cycle.
The recognition of the importance of sub-daily SST variability to larger scale phenomena has led to increased interest in characterising (e.g. Kennedy et al., 2007) and modelling diurnal warming, using both physical and statistical approaches (see the review by Kawai and Wada, 2007).
The statistical approach to modelling the diurnal cycle that we present here is arises from , which documented the frequency of diurnal cycle amplitudes of different magnitudes in the summer western Mediterranean Sea and European shelf seas. The peak of the diurnal cycle typically occurs around 14:00 LT. It was pointed out that the largest events (>4 K, say) happen much less frequently than dead calms (i.e. a wind speed of zero, the winds being from numerical weather prediction, NWP). This statistical observation is consistent with what we understand about the dynamics of diurnal warming from mooring observations and turbulence modelling: the warmest warm layers occur when the wind speed has been persistently close to zero throughout a period of sufficiently large solar heating. This is a less common event than the wind speed being close to zero at a given time. A consequence is that any statistical model based on regression against the instantaneous wind at time of observation will fail to correctly capture the distribution of diurnal warming amplitudes.
Given this, an obvious approach is to use NWP or other temporally resolved wind fields matched to the observations of the diurnal cycle. This is the approach we adopt in this paper, using NWP fields from the European Centre for Medium-range Weather Forecasting (ECMWF) and observations from a sensor on a geostationary satellite. However, we do not use conventional regression, and it is useful to comment at this stage on why we do not do so. The limitation that would arise using conventional regression is that errors in temporally resolved wind fields result in a distortion of the distribution of predicted diurnal warming, because the wind field errors are asymmetric (contravening the assumptions behind conventional regression). There are two aspects to this asymmetry. First, diurnal warming occurs when the magnitude of wind stress is small, and so in the relevant regime errors are inevitably asymmetric. (To see this, consider the limiting case of the true wind speed being zero; any error in the NWP wind speed can then only be positive.) Second, the error in a wind field includes spatial displacement of features such as wind minima relative to their true position. Large diurnal warming events may be observed in locations where the NWP winds are not at their weakest, therefore; and diurnal warming may be modest where the NWP is minimal. Thus, spatial mismatch errors lead the regression model to flatten out the dependence of diurnal warming on wind speed, even if using an appropriate measure of persistent wind speed.
In this paper, therefore, we take an approach aimed at ensuring that the correct statistical distribution of diurnal warming amplitudes is predicted by the model. More precisely, it is designed such that amplitudes of a certain magnitude will arise at the correct frequency (when driven by NWP fields consistent with those used to derive the model coefficients). Moreover, to the extent that the observations and NWP fields used are realistic, the model we derive captures the dependence of diurnal variability on persistent wind stress at different times of day (although with some limitations to be discussed).
The above discussion focuses on the wind dependence of diurnal warming. As will be seen below, the model also represents the effect of cumulative net heating of the sea surface, which is also available from NWP.
In this introduction we have attempted to give the motivation for this study, both in terms of the importance of diurnal variability in various contexts, and in terms of the need for a statistical model that encapsulates aspects that variability whose distribution is not captured in other statistical models. The next sections describe the development and properties of the model in detail.

A statistical model of diurnal variability
We define the diurnal warming (D) at a particular time of day as the difference between the SST at that time and a reference temperature. The reference temperature is the SST soon after dawn at the end of the period through which the ocean (usually) has cooled overnight. In practice, our observations of SST are available hourly (see below), and we find it adequate to take the reference SST as the last hourly observation before the net heat flux into the ocean becomes positive (as insolation increases during the morning). We consider that this reference temperature approximates closely the foundation temperature concept (the SST at a time and depth showing no influence of any diurnal warming), thus the model will not be applicable in cases where there is a warm layer that persists from day to day.
The SST observations are at infra-red wavelengths and respond to the ocean skin temperature, and observations of D are differences between these SSTs. We do not separately consider any diurnal variation in the cool skin effect (the fact that the air-sea interface temperature is ∼0.2 K cooler than the water a few millimetres below). Fairall et al. (1996) and others have argued that there is an effect of direct solar heating on the skin effect, although this seems to be contradicted by recent observations (Minnett et al., 2011). It is clear that the average cool skin is of a greater magnitude (e.g. ∼0.5 K) at very low wind speeds, tending to offset partially the temperature increase associated with any underlying warm layer. Systematic cool skin variability is thus part of the signal described by this model. The diurnal variation of SST (D) depends on many variables, the most important being the net heat flux during the day and the wind speed. The model expresses the diurnal warming (D) as a function of hour of day (t), integrated net heat flux (Q) since net heat flux (q) became positive (near dawn) and maximum wind speed since net heat flux became positive (W ).
The measure used for the heating term, Q, is the integrated net heat flux (solar, long-wave, sensible and latent), starting the integration when the net heating becomes positive soon after dawn. For a constant depth diurnal layer, Q minus the portion of the solar heating absorbed below this constant depth would be proportional to the SST change. The measure of wind mixing used, W , is the maximum wind speed (from 6-hourly data) since the heat integration started. Using W introduces some sensitivity to the wind history, particularly to the degree of persistence of low wind speeds, which is not captured by using the instantaneous wind.
The functional form of the model was developed empirically from the data. The method is detailed in Sect. 3, below. The functional form derived for parametrising warm-layer formation is This is applied for heating periods when Q > 0, and the coefficients a, b, and c depend on the hour of day (see Appendix A). This form was chosen in the light of inspection of the observations D against W stratified by Q and time of day. For cooling periods (Q < 0) the decrease in SST is represented as where ρ and c p are the density and specific heat of water and d is the climatological mixed layer depth from de Boyer Montegut et al. (2004). The coefficient f can be interpreted as the fraction of the climatological mixed layer that appears to cool on sub-daily time-scales in our observations. Thus, cumulated cooling is assumed to spread throughout a fraction f of the seasonal mixed layer by mechanical and convective mixing (and thus is small), whereas cumulated heating can create a shallow diurnal layer with larger rate of change of SST.
In practice, the diurnal warm layer profile evolves in depth during the day partly in response to wind stress, and the complex effect of this on the relationship between Q and D is fitted on average by the time-dependence of the coefficients a, b, and c. This average time dependence will be less accurate for day lengths outside the range (8 to 16 h) sampled in the observations used to fit the model.
A factor not accounted for (although presumably present in the observations of D) is the dependence on temperature and salinity dependencies of sea-water's thermal expansivity.
The range of latitudes of the observations fitted in the model is 50 • S to 50 • N, over which the annual mean expansivity varies by 30 %; the model will be less applicable at higher latitudes with colder waters. Nor does the model account for the variability of ocean optical properties affecting the depth of solar absorption.
The model is derived by matching, for each hour of day independently, the cumulative distribution functions for D (from satellite observations) and W (from an NWP model), binned by Q. For cases with Q > 0, the resulting empirical relationship between the parameters is then fitted to Eq. (1) for each hour of day independently. For cases with Q < 0, the empirical relationship between D and Q is fitted to Eq. (2).
The model has been derived for relatively large scale (0.25 • ) NWP and observational inputs, so is only likely to applicable to the same time and space scales. The model has no seasonal or regional (e.g. latitude) dependence, being derived from one year's data over the whole of the Atlantic and Mediterranean. The model coefficients would need to be re-fitted if different NWP fields were used (for example, the winds used here have a 6 h time resolution; using winds with a 3 h time resolution will give different values for W , closer to the true maximum wind speed) but the functional form would likely remain the same.
3 Details of the model and fitting procedure

General comments
The model is formulated in local time. The model is constructed using NWP fields and observational data from June 2007 to June 2008 and tested (Sect. 4) using fields and observations from May 2006 to May 2007.

Diurnal SST observations
The SST data are the hourly observations by the Spinning Enhanced Visible and Infrared Imager (SEVIRI) on the geostationary Meteosat Second Generation spacecraft (MSG) (Le Borgne et al., 2006). These are nominally observations of sub-skin SST, but variability in skin effect is not accounted for -i.e. the skin SST to which the sensor is actually sensitive differs from the reported sub-skin SST by a constant that cancels when an SST difference is formed. We use the observations mapped to a 0.05 • resolution grid, covering the area 50 • S-50 • N, 55 • W-45 • E with SST values every hour (UT). The retrieval of the SSTs from the SEVIRI radiances at 11 and 12 µm uses the same algorithm by day and by night, allowing the calculation of the diurnal variation. Only the SST at satellite zenith angles less than 60 • are used, since beyond this limit the noisiness of the SST retrieval increases.
The quality of the SEVIRI SSTs is indicated using a quality flag (cfd) which ranges from 0 to 5. The SEVIRI SSTs with "good" confidence (cfd ≥ 3) are used. This gives a fair balance between coverage and accuracy. Observations can be contaminated by Saharan dust in the atmosphere (Merchant et al., 2006), and (since April 2006) a correction has been applied to dust-affected pixels. This correction is strictly valid for night-time only (although a new technique will deliver day-time correction in a future reprocessing). The last night-time correction is simply preserved for 6 to 10 h, after which the pixel is given a lower confidence level (cfd = 2) because the Saharan dust will have moved significantly by advection. Thus, the afternoon SST values, when the peak diurnal warming would occur, are often given "bad" confidence in summer over the Mediterranean Sea, Red Sea, and tropical north Atlantic. Rejecting these pixels gives an unrealistic decrease in the average diurnal SST during the afternoon. Inspection of several days of data showed that these cfd = 2 pixels can be used in the diurnal SST calculation, if they are not consistently "bad". Therefore, we use afternoon observations for pixels which have fewer than 3 occurrences of cfd = 2 during the period of the day with valid Saharan dust correction.
The model is constructed to work in local time. The SST values at 00:00-23:00 UT are converted to local time by binning to the nearest local time (LT) hour. The reference temperature is the SST at the LT hour before the NWP heat flux becomes positive. The diurnal signal, D, at a later time is the later SST minus the SST at the reference time.
The resolution on which D is calculated for model construction and testing is 0.25 • . The diurnal warming D is calculated at full 0.05 • resolution and the median D in each 5 × 5 pixel block with more than 12 cloud-free pixels gives the 0.25 • resolution diurnal warming. This is the observational data used for the model construction and testing. The distribution of diurnal warming amplitude is dependent on the resolution of observation, but this resolution is adequate to resolve the length scales typically associated with events exceeding 6 K . This spatial averaging reduces the radiometric noise in the SST observations. It conveniently matches the spatial scale of the global Advanced Microwave Scanning Radiometer-E observations that we use in Sect. 4 to assess our model results for "all-sky" conditions. The estimated radiometric noise in the observational D at 0.25 • resolution is 0.1 K (compare this with the average peak warming of 0.5 K in the SEVIRI data used).

NWP fluxes and winds
The integrated net heat flux Q and the maximum wind speed W since net heat flux (q) became positive (which is generally just after dawn) are derived from operational forecast and analysis fields from the European Centre for Medium-Range Weather Forecasting (online). The ECMWF fields available are at approximately 1.125 • resolution (on a reduced Gaussian grid). The wind fields are the instantaneous 6-hourly analyses, the heat flux fields are the 24-h accumulated forecasts. Land and ice grid points were masked (rigorously excluding coastal points, which mix the very different land and For a given frequency, curves for larger Q tend to give larger D. Note that there is covariance of W and Q (from varying atmospheric stability and varying sensible and latent heat fluxes) and this is included in the resulting model. The concept of frequency matching is that (following the thin grey arrow), for a particular Q bin (bold curves), persistent wind speeds less than 1.4 m s −1 occur with frequency 10 −2 in the NWP data, which corresponds to the frequency of warming greater than 2.5 K. sea heat fluxes) and the fields interpolated to 0.25 • using the Climate Data Operators software (online). The wind speed was calculated and linearly interpolated to each LT hour. The heat fluxes were separated into solar flux and non-solar flux (long-wave radiation, sensible heat, and latent heat). The non-solar heat flux is the 24-h average derived from the 24-h accumulation. The solar flux is derived from the 24-hour average solar flux (derived from the 24-h accumulated flux) by modulating with the cosine of the solar zenith angle to give a realistic daily variation.
The integrated net heat flux Q is then calculated at each LT hour from when the net heat flux q (solar plus non-solar) first becomes positive (just after dawn). Until Q first becomes negative (late in the evening) a warm layer can form and then decay; after Q becomes negative there will be no warm layer (and the model simply assumes the cooling occurs throughout a fraction of the climatological mixed layer). Similarly, the maximum wind speed W since q became positive is calculated at each LT hour using the linearly interpolated 6-hourly instantaneous wind speed.

Frequency matching for warming model
We discussed above the problem presented to traditional regression by spatial mismatch between the true locations of wind speed minima (where diurnal warming is observed) and the locations of corresponding minima in the NWP fields. To avoid biased model coefficients, we therefore derive the coefficients by matching the overall distributions of W and D, at a particular value of (positive) integrated heat flux. We call this frequency matching. This method is the same as the probability matching method used to correlate rain rate measurements with radar reflectivity (Calheiros and Zawadzki, 1987;Rosenfeld et al., 1994). The hypothesis is that, at a given time of day and for equal integrated heat flux, D is a function of W (low W will give large D) and that this function can be found by matching the probability of wind speed less than W to the probability of warming greater than D. This is valid if NWP fields have a realistic distribution of wind speed in the vicinity of diurnal warming events over all, without requiring that the precise location of wind minima be correct.
The sampling in time and space is identical between the satellite observations and the NWP fields so that the probability distributions can validly be compared -i.e. the W is sampled only where we have valid D. The interpolation of 1.125 • NWP fields to 0.25 • to match D means that there are spatial correlations in W at the finer resolution output grid. The integrated heat flux is smoother spatially than the wind speed, and simple binning of the results by Q appears to be satisfactory.
The  The first step of directly frequency-matching the W and D distributions is illustrated for 14:00 LT in Fig. 1. Matching the curves point-by-point gives a functional relationship between W and D for each Q, see Fig. 2, from which an empirical relationship between D, Q and W is derived. For low wind speeds, D is proportional to Q, which is physically reasonable. Several functional forms for the W dependence were tested and 1/D was found to be well fitted by a quadratic function of W . Since D ∼ Q/L, where L is the warm layer depth, this implies that L ∼ W 2 . Although one might expect that the depth of the warm layer should be proportional to the Monin-Obukhov length (l ∼ w 3 /q), Noh (1996) showed that, because the turbulent dissipation increases with l, the warm layer depth is not a linear function of l and is in fact close to l 2/3 ∼ w 2 for moderate wind speeds (see Fig. 6 of Noh, 1996).
Extrapolating to zero Q in Fig. 2a we see that D is not zero at zero Q. Likewise, D becomes negative at large W . Both these limiting behaviours are unphysical and arise as a result of errors in the observed D. This is because the simple frequency matching above is one-sided: small values of W are matched to large values of D. Uncertainty in observations of D will therefore give rise to an apparent W dependence even for the sub-set of D observations whose true value of D is zero.
Our best information on the error distribution on D is the distribution of D for near-zero heat flux (lowest Q bin). Figure 3 shows the probability density of this distribution for 14:00 LT. Below, we first make some comments on the sources of error in D, and then describe how we adjust the frequency matching procedure for this effect.
The central portion of the distribution in Fig. 3 is closely fitted by Gaussian noise with standard deviation varying from 0.1 K to 0.2 K over the day. This is broadly consis--1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 D / K 10 -4 10 -3 10 -2 10 -1 10 0 probability density Fig. 3. Normalised probability density of observed warming D for near-zero integrated heat flux Q at 14:00 LT (thin curve) and fitted Gaussian-exponential density function (thick curve). The median warming is 0.08 K and the median integrated heat flux is 1.7 MJ m −2 . tent with estimated radiometric noise in the SST observations and the degree of averaging used in calculating D at 0.25 • . Other factors may also contribute, however. The atmospheric correction implicit in the retrieval of SST from the SEVIRI radiances decorrelates on synoptic timescales as atmospheric systems move. The seasonal trend in foundation temperature is of order 0.05 K over one day. These provide additional variability that will increase with time from the reference time.  Fig. 4. Fit between the probability distributions of modelled warming > D and observed warming > D at 18:00 LT for small integrated heat flux (4.3 MJ m −2 ). The modelled warming ranges from zero to 0.8 K. When the background noise is added, the modelled distribution matches the observed distribution. The fit is poor at the very largest values of warming but does include non-zero probability of negative warming.
In addition there are types of error capable of causing the non-Gaussian tails of the D distribution in Fig. 3. An example which increases steadily with time since the reference time is advection of horizontal variations in SST. Mostly, SST gradients are modest and this effect will be small. But for a fraction of cases, advection of strong fronts could cause a variability in SST up to 0.3 K in 24 h, easily affecting the tails of the frequency distribution. Likewise, cloud screening and the treatment of Saharan dust are not perfect and occasionally will give rise to outliers in D. There may be a contribution from errors in the fluxes and winds in the NWP (the variation in D from the finite bin size for Q is estimated to be minor). Finally, there will be high-frequency temporal variations in wind speed and heat flux that are not resolved in the low time resolution NWP fields used; these variations will also give rise to uncertainty in D estimated by the statistical model.
We fit this distribution of D for the near-zero heat flux bin with a Gaussian plus exponential error model, the exponential part accounting for the observed tails in Fig. 3. The distribution is slightly shifted to positive values of D because the heat flux is slightly larger than zero. We assume this distribution (shifted to have zero mean) adequately describes the error distribution for the full range of Q, that is, that this error distribution is convolved with the real D distribution to give the observed D distribution. We find it is not possible because of numerical instability to deconvolve the fitted error distribution from the observed D distribution, and therefore adopt a "forward" approach to deriving the coefficients for our model, as follows. (1) with an assumed set of coefficients and apply this model to the observed W distribution to get a modelled "true" D distribution, for a given LT hour. This is then convolved with the fitted noise distribution to get a modelled "observed" D distribution. This is then repeated many times for different assumed sets of coefficients sampling the plausible parameter space. The selected model coefficients are those that give the least squared deviation between the modelled and actual D distributions. The least squares deviation is weighted by the distribution itself, equivalent to using number weighting. The fitting is performed on the range D > median(D), because these are the data where there is a sufficient variation in D for successful fitting. Examples of the resulting fit between the modelled and observed distributions for small and large Q are shown in Figs. 4 and 5 (at 18:00 LT to show the effect of the background noise when the modelled warming is small).
The fit between the modelled warming (with background noise added) and the observed warming is shown for 14:00 LT in Fig. 6; the results are similar for other times of day. The poorer fit at low wind speeds seen in Fig. 6b is probably because the number of observations becomes very small and the results are then dominated by noise. At high wind speeds and middling integrated heat flux the model also diverges from the observations, see Fig. 6a. This is likely due to the fit being restricted to D > median(D) as noted above, although it is not clear why the results are poorer in this particular range. Similar curves for the model without background noise (not shown), show similar behaviour but D is zero at zero Q and does not become negative with increasing W (at high wind speed the model curves match those produced by binning D against Q and W ).

Cooling
By construction, the frequency matching model is only applicable in cases where the integrated heat flux since dawn is positive. When the integrated net heat flux drops below zero (typically late in the evening) we assume that any warm layer formed during the day has been completely removed by cooling and convective mixing. At this point we assume that wind-driven and convective mixing will distribute the cooling throughout the seasonal mixed layer and that D will be given by Eq. (2). It was found that the constant of proportionality f was very close to one.  with the variance of the observations, the background variability, approximated by the variance of the "zero heat flux" D described above, is added to the variance of the model. With this correction, the model and observation standard deviations show fair agreement. The main reason to develop the statistical model was to produce a realistic distribution of large warming events. From the agreement of the model and observation 95-percentile curves shown in Fig. 7, this appears to have been achieved. Note that, as for the variance, the model distribution needs to have background variability added (by convolving the distributions) to match the observed distribution. Due to the background variability we do not expect the model to perform well beyond one day and so D is reset to zero each day at dawn. The average geographical pattern of all-sky diurnal warming from the model was compared with AMSR-E during the period June 2006-May 2007. AMSR-E is a microwave imager on the polar orbiting NASA Aqua spacecraft and observes SST globally, in clear and cloudy conditions, with one observation at night and one in day-time at any one point of the Earth. We use the AMSR-E observations of SST which are available from Remote Sensing Systems (online). The observations are on a 0.25 • longitude by 0.25 • latitude grid. The day-night difference at each grid point was calculated and used as measure of the diurnal warming. The times of observations are approximately 01:30 and 13:30 LT for night and day respectively. For each day the diurnal model was used to estimate the warming from dawn at the integral hour closest to the day-time observation. The cooling between the night-time observation time and dawn was calculated using the model for the previous day. The magnitude of this cooling is ∼0.1 K on average. Since the model does not extend beyond 06:00 LT of the next day, the night-time cooling is approximate when dawn is later than 06:00 LT (i.e. high latitudes in winter). The averages over the period June 2006-May 2007 of day-night differences at each AMSR-E grid point are compared in Fig. 8. The global mean, standard deviation, and 95-percentiles of the modelled and observed day-night differences are given in Table 1. Note that the model values do not have any estimate for AMSR-E uncertainty (random noise and discretisation error) added, so the modelled standard deviation and 95-percentile will be less than the observed AMSR-E values.

Validation of the model
In the tropics the magnitude and patterns of diurnal warming are similar in the model and AMSR-E, although there are regional biases in both directions. At high latitudes the model shows little diurnal warming, while AMSR-E shows consistent warming near 0.2 K and sporadically higher values. In the Southern Ocean, where high wind speeds generally prevail, the model warming near zero seems more reasonable. The large warmings observed in parts of the northern hemisphere may be due to enhanced absorption in turbid water, which is not accounted for in the model. The other main dif- ference between model and observations is in the southern sub-tropics, where the model estimates larger warming than observations in the sub-tropical stratus regions. This appears to be caused by regional biases in the NWP model fluxes. The versions of the ECMWF model used predict 15-20 % less cloud than the ISCCP D2 climatology in these regions, as discussed by Köhler et al. (2011, Fig . 7), which would lead to an over-estimate of the solar flux in these regions. Since the statistical model is derived from global statistics, regional biases in the NWP model fluxes will lead to regional biases in the predicted warming, in this case a larger warming than observations.

Applications
Applications of a diurnal warming model include: adjustment of satellite observations of SST to a standard local observation time for climate studies, calculation of a background SST for assimilation of SST observations and for the retrieval of SST from satellite radiance observations, and parametrisation of diurnal warming in free-running climate www.ocean-sci.net/8/197/2012/ Ocean Sci., 8, 197-209, 2012 models. A further application of the statistical model described above is as a convenient summary of SEVIRI observations for testing the performance of other models over a wide range of conditions. When comparing a diurnal model driven by NWP fields with observations, the differences can arise from errors in the NWP fields used (predominantly coincidence errors in the wind fields) and not only from deficiencies in the model. It can be hard to attribute the errors between these sources. Using the statistical model driven by the same NWP inputs as the tested model will produce a set of synthetic observations that can be used to help interpret model-observation differences.
Here, we have compared the results of the statistical model for an idealised test situation against the Profiles Of Surface Heating (POSH) model  and the model of Zeng and Beljaars (2005) (hereafter, the "ZB model"). The POSH model is based on the bulk model of Fairall et al. (1996) with varying profiles of temperature in the warm layer. The ZB model is based on Monin-Obukhov theory with a fixed profile and warm layer depth. The test situation is as follows. The wind speed is constant throughout the day, set to a value ranging from 1 to 6 m s −1 . The net solar heat flux is proportional to cos(π/2·(t − 12)/6) between 06:00 and 18:00 LT, zero otherwise, with the mean net solar heat flux ranging from 100 to 320 W m −2 . The net non-solar heat flux (long-wave IR radiation, sensible heat flux and latent heat flux) is constant throughout the day and for all cases is set to −100 W m −2 . Note that the net non-solar heat flux used as input to the statistical model is the flux which would be measured if there were no diurnal warming, since this is implicit in the NWP fields used. The flux input to the other models should be that which would be measured including any effects on the fluxes due to the diurnal warming; the importance of this is discussed further below.
In Fig. 9 we compare the peak warming (approximated as the SST difference between 06:00 and 14:00 LT) for each model. The POSH model was developed using IR measurements of the skin SST reduced to sub-skin SST by subtracting the skin effect parametrised by the method of Donlon et al. (2002). Since this parametrisation depends only on the wind speed, the SST differences for the POSH model are equal to skin SST differences in the test situation of constant wind speed throughout the day.
The results shown in Fig. 9 indicate that the ZB model has similar overall magnitude to the statistical model, but predicts larger warming at moderate wind speeds (∼5 m s −1 ), as noted by Bellenger and Duvel (2009), while predicting smaller warming at low wind speeds (<3 m s −1 ). This can be interpreted as a result of the use of a fixed depth scale for the stratified layer. At moderate wind speeds, a somewhat deeper diurnal mixed layer is formed than assumed by ZB, and thus the magnitude of warming is over-estimated by the model (for a given amount of heating). At low wind speeds, diurnal warm layers can be shallower, and therefore warmer, than assumed in the ZB model. The POSH model slightly under-estimates the warming at larger wind speeds compared to the statistical model. But most notably, POSH predicts very large warming of up to 8 K at the highest insolation and lowest wind speed. One possible explanation for this difference is spatial scale. The POSH model was developed specifically with reference to "point" ship-borne measurements of diurnal warming, while the 0.25 • resolution of the observations used to construct the statistical model will average over the highest, most localized peak warmings . However, there were no observations available to constrain POSH in the regime of warming exceeding 5.2 K . Thus, it is reasonable that POSH may have systematic errors (apparently an over-estimation) in this regime. Before firmly concluding this, however, it is necessary to consider the effect of driving the POSH with idealized fluxes that do not respond to the diurnal variation itself.
For these large values of warming, the perturbation of the surface fluxes by the diurnal warming itself is large enough to modify the results. To show this we used the POSH model together with the Donlon et al. (2002) skin effect and the COARE 2.5 model for the atmospheric surface layer (Fairall et al., 2010) to adjust the flux as the model stepped through the day. We assumed that the air-sea temperature difference remained constant at −1 K (Chen and Houze, 1997) and that the relative humidity remained fixed at 80 % -i.e. that air temperature and absolute humidity are modified in concert with the diurnal warming of the surface water. The resulting perturbation in the downward radiative flux was estimated using the method used for SEVIRI radiative flux retrieval (OSI SAF, 2005). The sign is such that the effect of the surface warming is to reduce the net flux into the warm layer. Including this negative feedback via the fluxes reduces the peak warming modelled by POSH from 8 K to 6.5 K, but has little effect on peak warmings less than 3 K, see Fig. 9d. Thus, this feedback brings POSH closer to the statistical model, but does not change the ZB model results significantly. The assumption of constant air-sea temperature difference will result in a smaller perturbation to the fluxes than if the amplitude of the cycle in air temperature is less than that in SST. So it is possible that negative feedback is somewhat larger, which could further improve the agreement between the POSH model and the statistical model.

Discussion and conclusions
The statistical model presented here has two distinctive features: the fitting of the model to observations is done on the basis of frequency distributions of events of different magnitudes, and the effect of wind speed is parametrised using the maximum 6-hourly wind speed during the period between the time when the net heat flux into the ocean became positive (generally, this is shortly after dawn) and the time of the model prediction. Using the maximum in wind speed  Zeng and Beljaars (2005), (c) the Profiles of Surface Heating (POSH) model , and (d) the change in the warming predicted by the POSH model when interaction with the atmosphere is included (the values range from −1.5 K to 0 K). captures the need for persistently low wind speeds to explain the incidence of large diurnal warming events, but has the disadvantage that the model predictions are less precise beyond the late afternoon. For example, if there is high wind speed in the morning and low wind speed in the afternoon, there will be some residual warming by early evening that the statistical model does not capture; conversely, if there is high wind speed in the afternoon, any warming will be mixed through a deeper layer and will be smaller than that predicted by the model. A possible future development of the model is to parametrise the model to account for the history of wind speed during the day, by fitting the model for the maximum wind speed in several periods (e.g. 06:00-12:00, 12:00-18:00, 18:00-24:00, and 24:00-06:00 LT), with a dependence on the warming at the beginning of the period. The functional form would likely be different but the frequency matching method can still be used. There will still remain uncertainty from seasonal effects and advection.
Matching the frequency distributions bypasses the problem that NWP model winds are imperfectly spatially coincident with the observed warming. However, the method is sensitive to noise and the estimate of the noise distribution is assumed to be valid at all values of integrated heat flux Q. There are indications, see Fig. 6a, that for medium values of Q and high wind speeds, the noise (or variability) in observations may be larger. The model does reproduce the statistics of the observed warming, including the large values (>3 K) regularly seen in the SEVIRI observations and thus a possible use of the model is to invert the SEVIRI observations to give estimates of the surface wind speed for NWP analyses or for assimilation into NWP models. The frequency matching method is of general applicability and could be with other models to "tune" them to large scale observations.
The model is rapidly calculable (being a statistical fit) and differentiable with respect to persistent wind and cumulated heat flux (potentially useful in variational analysis of diurnal variability combining model and observations). late afternoon onwards. It is also less applicable for latitudes polewards of ∼50 • north and south, outside the latitude range of the observations of diurnal warming used in its formulation. The model will not meet every need for diurnal warming estimation, but by design, has two particular strengths. First, since it captures the distribution of diurnal warming and the dependence on wind and insolation of this distribution, the model is a useful comparison for exploring the parameter space of physical models and examining the distribution of diurnal warming they predict. Second, the model is relevant for calculating the mean effect of diurnal variability on processes that are non-linear in SST, since in such cases the full distribution of diurnal variability determines the mean effect.
We have shown that the statistical model can be used as a surrogate for observations when testing other, more physically based models. A complementary way to test these models has been developed by Bellenger and Duvel (2009) who calculated the frequency distribution of warming for two models using NWP fields from the ECMWF ERA-40 re-analysis and compared the distributions with those from surface drifting buoy measurements. Their results, and those of Bernie et al. (2008), show that the inclusion of diurnal warming in climate models is required to reproduce mean and intra-seasonal effects. The model described in this paper is very simple to implement and reproduces the observed statistics of diurnal warming and is a candidate parametrisation for use in climate models. Table A1 gives the coefficients for the warming model, Eq. (1). The fit for hours before 09:00 LT and after 24:00+05:00 LT were very poor or did not converge, so the coefficients for these times are set to the values for 09:00 LT and 24:00+05:00 LT, respectively. The coefficient f in Eq. (2) was very close to 1 for all times of day, so we take f (t) = 1.