Working with time-of-day as a factor will be a challenge because I assume your sampling times are random randomly spaced. That means you might be better off fitting for a trend plus some model of diurnal variation. The problem is that latter can take many functional forms, depending on solar forcing, whether there is nighttime convection, and so on.

I would advise starting by regressing to a formula that has a trend, plus both a sine() and a cosine() with period one day. (Using both sin and cos lets you find the phase, but still in a linear regression, which is easier to deal with than a nonlinear one, unless you are experienced.) Then add a twice-daily harmonic. Then add a thrice-daily harmonic. By comparing the results you may start to learn how the diurnal variation is affecting your estimate of the trend.

I would start with artificial data, so you can test the accuracy. AIC or other methods should be used to avoid overfitting. (If you don't know what AIC means, you ought to give yourself some time to learn elementary stats, before trying anything complicated with your data.)

Whether you will see a linear trend is questionable. I've never seen such a simple pattern in ocean data. So you might want to start by fitting this sort of sequence of diurnal models, then removing the diurnal signal, and then simply plotting the results, rather than expecting a trend.

Have you tried to add the time of sampling as one of your fixed factors during the analysis?

Olga Bykova Thanks a lot.

If it is a time series, you might apply a diurnal filter to your data, that way you can extract the diurnal signal. If you only have discrete samples maybe you should use the anomaly?

Thanks Macarena Díaz Astudillo . It is not a time series but discrete stations. How would the anomaly work?

You first have to calculate the mean temperature of your data, and the substract the mean to each measure. This will give you the deviation from the mean, which is the anomaly. For more info see:

https://www.ncdc.noaa.gov/monitoring-references/dyk/anomalies-vs-temperature

Working with time-of-day as a factor will be a challenge because I assume your sampling times are random randomly spaced. That means you might be better off fitting for a trend plus some model of diurnal variation. The problem is that latter can take many functional forms, depending on solar forcing, whether there is nighttime convection, and so on.

I would advise starting by regressing to a formula that has a trend, plus both a sine() and a cosine() with period one day. (Using both sin and cos lets you find the phase, but still in a linear regression, which is easier to deal with than a nonlinear one, unless you are experienced.) Then add a twice-daily harmonic. Then add a thrice-daily harmonic. By comparing the results you may start to learn how the diurnal variation is affecting your estimate of the trend.

I would start with artificial data, so you can test the accuracy. AIC or other methods should be used to avoid overfitting. (If you don't know what AIC means, you ought to give yourself some time to learn elementary stats, before trying anything complicated with your data.)

Whether you will see a linear trend is questionable. I've never seen such a simple pattern in ocean data. So you might want to start by fitting this sort of sequence of diurnal models, then removing the diurnal signal, and then simply plotting the results, rather than expecting a trend.

It depends on your study area, but you might be able to use time series of in situ measurements. For example, you could download data from the National Data Buoy Center (https://www.ndbc.noaa.gov/). Find the buoy closest to your study area and quantify temporal variability in the data from observations at the buoy's location. You can also check Copernicus (http://marine.copernicus.eu/) to see if there is/was a glider or other in situ measurements in your area. Also check GOSUD (http://www.gosud.org/), NOAA's global drifter program (https://www.aoml.noaa.gov/phod/gdp/). Models are also an option, like NEMO, which is also available through Copernicus. You could examine temporal variability at synoptic and longer scales using AVHRR SST and/or optimally interpolated SST products. Satellite-derived SST will also quantify spatial variability of temperature in the area of your cruise track.

Dan's suggestion is excellent. From another approach, I wonder if it helps to infer decorrelation time and length scale from observations, and attempt an objective analysis.