I can now, with less computer calculation than I'd feared, answer Question 1. (I'll say more about Question 2 later).
Lemma:____ Let V be the vector space over Z/2 spanned by the C(r1,r2,r3) and the C(s1,s2,s3,s4). If an element of V has its power series expansion divisible by x^(l^2), it is 0.
To see this, let K be an algebraic closure of Z/2, S' be the subring of K[[x]] generated over K by the [j] and L be the field of fractions of S'. I can show that L/K is the function field of a curve, that there are exactly l(l-1)(l+1)/24 valuation rings in L/K that don't contain S', and that each of [1],...,[l-1] has a simple pole at each of them. So an element of V has at most l(l-1)(l+1)/6 poles in L/K, counted with multiplicity.
Also, the localization of S' at the maximal ideal generated by [1],...,[l-1] is dominated by exactly (l-1)/2 valuation rings in L/K. I can show that for each r prime to l there is an automorphism of L/K taking [j] to [rj] for each j, and that these automorphisms act transitively on this set of valuation rings. Now the elements of V are fixed by these automorphisms. So an element of V whose power series expansion is divisible by x^(l^2) has zeros of order at least l^2 at each of these valuation rings. Since (l^2)(l-1)/2 > l(l-1)(l+1)/6, such an element must vanish.
Suppose now that l=23. Since l=7 mod 16, my answer to my question "Existence of certain identities..." shows that C is in V. The lemma then shows that to prove 2. it's enough to show that C(3,3,1,2)+C(1,3,6)+C is divisible by x^529 in Z/2[[x]]. Now C(3,3,1,2)+C(1,3,6) is a sum of monomials in the [j]. Replacing each [j] by the sum of the first 2 terms in its power series expansion only modifies the sum by something divisible by x^529, and we're reduced to an easy computer calculation.
Establishing 3. and 4. is harder since we don't know in advance that C is in V. I'll use:
Theorem____Suppose there is an element R of V such that R+C is divisible by x^(d+1) where d=l(l+1)(l+1). Then R=C.
To see this, recall that F=x+x^9+x^25+..., that G=F(x^l), that H=G(x^l) and that C^2+C=G+H. Now there is a symmetric degree l+1 2-variable polynomial P over Z/2 with P(F,G)=0; furthermore P(z,H) is monic in z of degree l+1. This is discussed in my question "What's known about the mod 2 reduction...". Suppose a^l=1. Replacing x by ax^l in the identity P(F,G)=0, we get a root of P(z,H)=0 of the form ax^l+... This gives us l distinct roots, with G among them. By symmetry P(H,G)=0. So H(x^l) is still another root, and P(z,H) factors into linear factors over K[[x]].
Now R^2+R+G+H=(R+C)^2+ (R+C), and so is divisible by x^(d+1). Since P(G,H)=0, P(R^2+R+H,H) is divisible by x^(d+1). Also R^2+R has poles of order at most 8 at each valuation ring in L/K that doesn't contain S', while H has poles of order at most 12. It follows that P(R^2+R+H,H) has at most (l+1)l(l-1)(l+1)/2 =d(l-1)/2 poles, counted with multiplicity, in L/K. Arguing as in the proof of the lemma we see that it has more than d(l-1)/2 zeros, counted with multiplicity. So it vanishes, and R^2+R+H is a root of P(z,H)=0. Examining the roots of this equation we see that R^2+R+H can only be G. So R^2+R=G+H, and R=C.
Suppose now that l=31. To prove 3. we now see that it suffices to show that C(3,3,2,3)+C(2,3,7)+C is divisible by x^(186^2), since 186^2 >(31)(32)(32). This is carried out as in the case l=23, but now we have to use the first 12 terms in the power series expansion of each [j] rather than just the first 2. The treatment of 4. is similar, but now we show divisibility by x^(329^2) using the first 14 terms in the power series expansion of each [j].
EDIT___(UPDATE ON QUESTION 2)
Ira Gessel has now carried out computer calculations for all l<1500; here's a summary of his remarkable results. Consider the triples (r1,r2,r3) where the squares of r1,r1,r2 and r3 sum to l. To avoid duplicates, normalize each such triple so that r1,r2 and r3 are positive with r2>r3. Ira finds that for each such triple there is a unique second (normalized) triple (s1,s2,s3) such that the power series C(r1,r1,r2,r3)+C(2s1,s2+s3,s2-s3)+C is divisible by x^(l^2). Furthermore (r1,r2,r3)-->(s1,s2,s3) is an involution on the set of normalized triples with at most 1 fixed point.
When l=7 mod 16, the argument I gave above when l=23 then shows that for each normalized triple (r1,r2,r3) and the corresponding (s1,s2,s3) we have the identity C=C(r1,r1,r2,r3)+C(2s1,s2+s3,s2-s3); note that the squares of 2s1, s2+s3, and s2-s3 sum to 2l. So, for example, when l=1447 we get 17 distinct formulae for C.
When l=15 mod 16, it seems certain that once again C(r1,r1,r2,r3) and C(s1,s2,s3) sum to C. This could be proved by extending Ira's calculations to prove divisibility by x^(d+1) where d=l(l+1)(l+1), as in my treatment of l=31 and l=47. But I think this extension is unnecessary, and that one may deduce the identities simply from the divisibility by x^(l^2) established by Ira. If my idea for demonstrating this works out I'll post it as a comment.
But great mysteries remain. Why should this all be true? And can one describe the mysterious involution (r1,r2,r3)-->(s1,s2,s3) explicitly? By the way, it's known that the number of normalized triples is odd or even according as n is 7 or 15 mod 16. The proof of this goes back to Hasse; one shows that the number of triples is h/4 where h is the class number of Q(Root(-2l)) and uses results of Gauss on representations by sums of 3 squares, together with some genus theory for binary quadratic forms.