Bonnet theorem ensures that given symmetric bilinear forms $a$ and $b$ the necessary and sufficient conditions for a surface in $\mathbb{R}^3$ to have these forms as 1st and 2nd fundamental forms are the Gauss–Codazzi–Mainardi equations.
I'd like to have some geometric interpretation or intuition about such system. Why are they what they are, why is it that $a$ only enters through the Christoffel symbols, etc.
I do not know how to describe these equations purely geometrically. Below are the simplest explanations I know of. But I've omitted a lot of rigorous details.
The first fundamental form is simply the dot product on $\mathbb{R}^3$ restricted to vectors that are tangent to the surface. Note that the symmetry and positive definiteness of the dot product implies the same for the first fundamental form.
The second fundamental form can be described this way: The geometry of a surface is invariant under rigid motions. Therefore, given any point $p$ on the surface, you can move the surface so that $p$ is at the origin and the tangent plane at $p$ is the $xy$-plane. Then there is a function $f(x,y)$ such that the surface near $p$ is given by $$ z = f(x,y). $$ $f$ therefore has a critical point at the origin. The second fundamental form is then the Hessian of $f$ at the origin, $$ H(p) = \begin{bmatrix} \partial^2_{xx}f(0,0) & \partial^2_{xy}f(0,0)\\ \partial^2_{yx}f(0,0) & \partial^2_{yy}f(0) \end{bmatrix}. $$ That the second order partial derivatives of $f$ commute implies the second fundamental form is symmetric. The Codazzi-Mainardi equations now arise from the fact that the 3rd order partial derivatives of $f$ commute. In particular, if the surface is placed as described above, then the Codazzi-Mainardi equations at the origin are \begin{align*} \partial_yH_{xx} &= \partial_xH_{yx}\\ \partial_yH_{xy} &= \partial_xH_{yy}. \end{align*}
The Gauss equation is, as far as I can tell, more subtle. The idea is to find at each point on the surface a scalar or tensor invariant that depends only on the first fundamental form and is independent of the local coordinates used. If you do this, you discover both the Gauss curvature and the Gauss equation. You can find the argument for any dimension $n$ here. If you set $n=2$, the messy calculations become a lot simpler.