Something to get you started. Define an operator:
$$
I_x(f) := \int_a^b f(x) \, dx
$$
We can see that $I_x(f + g) = I_x(f) + I_x(g)$ and if $\alpha$ does not depend on $x$, $I_x(\alpha f) = \alpha I_x(f)$. In other words the operator is linear.
If it acts on $f(x, y)$ it returns a function of $y$ only.
By analogy define another operator
$$
I_y(f) := \int_c^d f(y) \, dy
$$
Applying $I_x$ and $I_y$ consequently we get a double integral:
$$
I_y(I_x(f)) := \int_c^d \left[ \int_a^b f(x, y) \, dx \right] dy .
$$
Now define a "Simpson's rule" operator:
$$
S_x(f) := (f(a) + 4 f((a+b)/2) + f(b)) (b - a)/6
$$
It is also linear and if it acts on $f(x, y)$ it "kills" dependency on $x$. The same for $S_y$:
$$
S_y(f) := (f(c) + 4 f((c+d)/2) + f(d)) (d - c)/6
$$
By analogy with integral let compute $S_y(S_x(f(x, y)))$ and call it Simpson's rule for double integral. Apply $S_x$ to f(x, y):
$$
S_x(f(x, y)) = (f(a, y) + 4 f((a+b)/2, y) + f(b, y)) (b - a)/6
$$
and apply $S_y$ to the result using linearity
$$
S_y(S_x(f(x, y)) = (S_y(f(a, y)) + 4 S_y(f((a+b)/2, y)) + S_y(f(b, y))) (b - a)/6
$$
expand all three terms with $S_y$
$$
S_y(S_x(f(x, y)) = \\
\left(16f\left(\frac{b+a}{2},\frac{d+c}{2}\right)+4f\left(\frac{b+a}{2},d\right)+4f\left(\frac{b+a}{2},c\right)+4f\left(b,\frac{d+c}{2}\right)
+f\left(b,d\right)+f\left(b,c\right)+4f\left(a,\frac{d+c}{2}\right)+f\left(a,d\right)+f\left(a,c\right)\right)
\frac{\left(b-a\right)\left(d-c\right)}{36}
$$