Reporting phase differences in percent is usually not a good idea, because phase either wraps around (0 to 360 degrees to 720 degrees etc...) or has a branch line (where 360 degrees equals 0 degrees for instance).
In either case, for instance, is 360 degrees 0%, 100 %, or infinite % error relative to 0 degrees?... you get the point I think.
The fact is that angles (phase) are ultimately dimensionless quantities. That is, 2 Pi radians is the ratio of the circumference of a circle divided by its radius. And angles furthermore wrap around or have branch lines. Both of these argue for reporting differences in phase simply as the difference in phase, in degrees or radians, and not as a percentage.
Having said that, you might take the complex Fourier Transform (FFT) of the time series or each within a time window spanning many tide cycles, and then look at the dominant peak at the tide frequency for each, which will have a magnitude and phase from the Fourier Transform. Simply take the difference between the two phases given by the FFT. If you slide the time window along in time, then you estimate the phase difference as a function of time, with time identified as the time at the center of the sliding time window.
If the peaks are broad (spread across a narrow band of frequencies near the tide frequency in the frequency domain), then take the weighted mean phase difference over that frequency band, where the weights are proportional to the spectral energy at each frequency in one of your signals. This way you focus the measure of the phase difference at the frequency where the tidal energy tends to be.
I am assuming that the tide level was digitally sampled at a fixed sampling frequency, and that a low-pass filter was expressly applied prior to sampling. If an irregular sampling period was used, then you will have to resample the signals at a fixed sampling period using interpolation between data points. If no low-pass filter was applied prior to sampling, then there is nothing you can do, your data is ruined by aliasing.
Doing a DFT is essentially least-squared fitting a set of sinusoids to a time sequence. So performing an FFT operation is already a "statistical" operation. So the complex coefficients that it returns are already averaged quantities.