On the Wikipedia page for Adaptive Simpson's method, the criterion for stopping the bisection of the integration interval $[a,b]$ when approximating an integral $\int_a^b f(x)\ dx$ is given:
$|S(a,e) + S(e,b) - S(a,b)|/15 < \epsilon$
where $S(x,y)$ is the result of applying Simpson's rule to the interval $[x,y]$, $\epsilon$ is the error tolerance, and $e = (a+b)/2$. If the criterion is not met, one applies the adaptive method to the intervals $[a,e]$ and $[e,b]$, with error tolerances $\epsilon/2$ in each, so the tolerances still add to $\epsilon$.
I want to perform adaptive Simpson quadrature for a function of two variables. I know how to apply Simpson's rule to a rectangle $[a,b]\times [c,d]$, and if some criterion is not met, I'll subdivide into four smaller rectangles $[a,e]\times [c,f], [a,e]\times [f,d], [e,b]\times [c,f], [e,b]\times [f,d]$ where $e=(a+b)/2$ and $f=(c+d)/2$.
What is the suitable criterion for stopping this subdivision, corresponding to the criterion for 1D above? And if this criterion is not met, what error tolerances should I have for the four new regions? I'm guessing $\epsilon/4$, is that reasonable?