Playing around a bit, it's clear that finding some solutions symbolically is easier if you accept some (hopefully satisfiable) constraints on the values of b1, b2.. and c1, c2... In effect, if we allow there to be only two unknowns parameters for "b" and "c" we can simplify things significantly.
Consider replacing the left hand side with b1 (y + z)^2 + (-b2) (y - 1) (z - 1) + b2. Expanding it, we that it is a restriction on the original system:
Expand[b1 (y + z)^2 + (-b2) (y - 1) (z - 1) + b2]
b2 y + b1 y^2 + b2 z + 2 b1 y z - b2 y z + b1 z^2
So our equations are:
b1 (y + z)^2 - b2 (y - 1) (z - 1) == a1-b2
c1 (y + z)^2 - c2 (y - 1) (z - 1) == a2-c2
We have a linear system with the unknown values of (y + z)^2 and (y - 1) (z - 1)
{res1, res2} = LinearSolve[{{b1, -b2}, {c1, -c2}}, {a1 - b2, a2 - c2}]
We know have a simpler system to solve:
Solve[{(y + z)^2 == res1, (y - 1) (z - 1) == res2}, {y, z}, Reals]
The results are enormous, heavily conditional, incomplete, and probably of no practical use if I did everything right. But maybe this gives you an idea of what some of the solutions you are asking for look like. It's also possible to change up the constraints and simplify the problem similarly.
