Thanks everyone, I think it is working correctly now. Below I am showing how I am doing it:

float height = 50.0 * sin(fs * waveVertex.x + t * _Time.y) +
125.0 * sin(fs * waveVertex.z + t * _Time.y + PI/4.0) +
135.0 * sin(-fs * waveVertex.x/4.0 - t * _Time.y * 0.75 + PI/2.0) +
165.0 * sin(-fs * waveVertex.z/8.0 - t * _Time.y * 0.5 + PI);
height *= 0.25;

float derivative = 150.0 * fs * cos(fs * waveVertex.x + t * _Time.y) +
125.0 * fs * cos(fs * waveVertex.z + t * _Time.y + PI/4.0) +
135.0 * -fs * cos(-fs * waveVertex.x/4.0 + t * _Time.y * 0.75 + PI/2.0) +
165.0 * -fs * cos(-fs * waveVertex.z/8.0 + t * _Time.y * 0.5 + PI);


float3 tan1 = float3(1.0, derivative, 0.0);
float3 bitan1 = float3(0.0, derivative, 1.0);

normal = normalize(cross(bitan1, tan1));

And this gives me the result as you see in the screenshot.

Which still doesn't look correct :S

If you have a height field y(x,z) and compute the derivatives
t(x,z) := d y(x,z) / d x
b(x,z) := d y(x,z) / d z
then t denotes the rate of change of y when stepping in direction of x, and analogously b denotes the rate of change of y when stepping in direction of z. Written as vectors, you'll choose the unit length 1 for the step in direction of x or z, resp., and t or b for the corresponding change of y, i.e.
T(x,z) := [ 1; t(x,z); 0 ]
B(x,z) := [ 0; b(x,z); 1 ]

The cross product of these 2 vectors gives you the normal. However, you have 2 possibilities T x B and B x T, but only one is valid for your purpose (depending on your choice of co-ordinate system). Above you've chosen
N(x,z) := B(x,z) x T(x,z) = [ -t(x,y); 1; -b(x,z) ]
It has to be normalized, of course.

So what does this mean w.r.t. some picked slopes?
a) For a flat height field you'll get
y(x,z) = c ==> t(x,z) = 0; b(x,z) = 0 ==> N(x,z) = [ 0; 1; 0 ]
b) For a wave in x direction only you'll get at a position with slope of 45° (i.e. increasing)
t(x,z) = 1; b(x,z) = 0 ==> N(x,z) = [ -1; 1; 0 ]
c) For a wave in z direction only you'll get at a position with slope of -45° (i.e. decreasing)
t(x,z) = 0; b(x,z) = -1 ==> N(x,z) = [ 0; 1; 1 ]

This seems me okay.


Now, there are some mistakes in your pre-previous post:

1. Notice that normals are not invariant to non-uniform scaling.

E.g. scaling the height by a factor of 0.25 is non-uniform, because it has an effect on y only but not on x and z. If you scale y(x,z) by s and compute N as above, you'll see that
N'(x,z) = [ -t(x,y) * s; 1; -b(x,z) * s ] != N(x,z) * s

2. The 1st sin() is scaled by 50, while the 1st cos() is scaled by 150.

3. The argument in the 3rd sin() has a frequency of -fs/4, but the factor of the 3rd cos() shows -fs only.

4. The argument in the 4th sin() has a frequency of -fs/8, but the factor of the 4th cos() shows -fs only.

5. The wave function depends on both x and z position, but you compute only one derivative and use it for both directions.

haegarr,

thank you for your help! I had 2 mistakes in my current code (the points from 2 - 5, I had corrected them before):
I was calculating the normal as:
normal = normalize(float3(derivativeX, 1, derivativeZ));
where it should have been with minus sign, as you said in your post.

The second mistake, was the scaling of the final accumulated height. I didn't expect it to have such a big impact!

So, thank you for your help, here is a screenshot of how the water looks now