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