I figured out the mistake that I was making and confirmed the hypothesis in my comment above. I was rendering each quad with its distinct set of vertex normals calculated from the local parametrisation. With lighting enabled, the vertices were therefore shaded differently for adjacent quads owing to the different dot products between the local normal and the light direction vector.
What I did was average out the normals to create a unique normal at each vertex. Now when I fly into the clouds, I do not see the quad boundaries.
As Dan suggested, I could avoid the problem altogether by disabling lighting, but it would eliminate the ability to simulate different times of the day and in any case, the fix above worked and only involved about 30 lines of code!
Problem solved.
In case anyone else finds this useful, here's the code (I wrote it in VBA) to calculate averaged vertex normals given a set of vertices on a surface that is topologically rectangular.
Sub Calculate_Averaged_Vertex_Normals(Vertx#(), M%, N%, VnormX#())
' INPUT:
'L,M: Array Dimensions
' Vertx - 3D array of vertex coordinates
' OUTPUT: VnormX: 3D array of Averaged Vertex Normals
Dim i%, j%, k%, p%, q%, w%, NormX#(3), Ncount%(), Vx#(), Fn#(3),Gn#(3),Hn#(3),A#
ReDim Ncount%(0 To L, NF + 1), Vx(0 To L, NF + 1, 4, 3)
For k = 0 To M - 1
For i = 0 To N - 1
Call Quad_Normal_Vector_Parametrisation(Vertx,k, i,Fn,Gn,Hn) ' Local Normals at each Vertex
p = k: q = i: Ncount(p, q) = Ncount(p, q) + 1: j = Ncount(p, q)
Call Quad_Vertex_Normal(0, 0, NormX,Fn,Gn,Hn)
For w = 1 To 3: Vx(p, q, j, w) = NormX(w): Next w
p = k: q = i + 1: Ncount(p, q) = Ncount(p, q) + 1: j = Ncount(p, q)
Call Quad_Vertex_Normal(1, 0, NormX,Fn,Gn,Hn)
For w = 1 To 3: Vx(p, q, j, w) = NormX(w): Next w
p = k + 1: q = i + 1: Ncount(p, q) = Ncount(p, q) + 1: j = Ncount(p, q)
Call Quad_Vertex_Normal(1, 1, NormX,Fn,Gn,Hn)
For w = 1 To 3: Vx(p, q, j, w) = NormX(w): Next w
p = k + 1: q = i: Ncount(p, q) = Ncount(p, q) + 1: j = Ncount(p, q)
Call Quad_Vertex_Normal(0, 1, NormX,Fn,Gn,Hn)
For w = 1 To 3: Vx(p, q, j, w) = NormX(w): Next w
Next i
Next k
For k = 0 To M
For i = 0 To N
For w = 1 To 3
A = 0#
For j = 1 To Ncount(k, i)
A = A + Vx(k, i, j, w)
Next j
VnormX(k, i, w) = A / Ncount(k, i) ' Averaged Vertex Normal
Next w
Next i
Next k
End Sub
Sub Quad_Normal_Vector_Parametrisation(Vertx#(),k%, i%,Fn#(),Gn#(),Hn#())
Dim j%, Bv#(3), CV#(3), Dv#(3)
For j = 1 To 3
Bv(j) = Vertx(k, i + 1, j) - Vertx(k, i, j) ' Tangent Vectors from Bilinear Parametrisation
CV(j) = Vertx(k + 1, i, j) - Vertx(k, i, j)
Dv(j) = Vertx(k + 1, i + 1, j) - Vertx(k, i + 1, j) - Vertx(k + 1, i, j) + Vertx(k, i, j)
Next j
Call CrossProd(Bv, CV, Fn) ' Normal Vector from Cross Product of Tangent Vectors.
Call CrossProd(Bv, Dv, Gn)
Call CrossProd(Dv, CV, Hn)
End Sub
Sub Quad_Vertex_Normal(x#, y#, NormX#(),Fn#(),Gn#(),Hn#())
Dim j%, Norm#
Norm = 0#
For j = 1 To 3
NormX(j) = Fn(j) + Gn(j) * x + Hn(j) * y
Norm = Norm + NormX(j) * NormX(j)
Next j
Norm = Sqr(Norm)
For j = 1 To 3
NormX(j) = NormX(j) / Norm
Next j
End Sub
Sub CrossProd(x#(), y#(), Z#())
Z(1) = x(2) * y(3) - y(2) * x(3)
Z(2) = y(1) * x(3) - x(1) * y(3)
Z(3) = x(1) * y(2) - y(1) * x(2)
End Sub
Any suggestions on how to get around this artifact would be appreciated. I'm using 50 grid segments in both the azimuthal and polar directions of the ellipsoid which is actually truncated 10 degrees below the poles. The texture map uses GL_MODULATE with GL_LINEAR and GL_LINEAR_MIPMAP_LINEAR for the Min and Mag filters. The ellipsoid base colour is white, so that there are no shading issues involved. I have also verified that the texture coordinates are correctly applied and the normals are oriented correctly.