Error of approximating double integral with its Riemann sum

I have a double integral which I want to evaluate numerically by approximating this integral with its Riemann sum. The integral is of the following form $$\int_{0}^{a}\int_{0}^{b} f(x,y)dxdy\approx\sum_{n=1}^{K}\sum_{m=1}^{J}f(n\delta_x,m\delta_y)\delta_x\delta_y$$

There are upper-bounds on error of the approximation for one integration such as $\int_{a}^{b} f(x)dx$. However I didn't find any bounds on error of double integrals.

How can I choose $\delta_x$ and $\delta_y$ such that:

$$\left|\int_{0}^{a}\int_{0}^{b} f(x,y)dxdy-\sum_{n=1}^{K}\sum_{m=1}^{J}f(n\delta_x,m\delta_y)\delta_x\delta_y\right|\leq \epsilon$$ for a predefined $\epsilon$ knowing $f(x,y)$?

PS: In my actual problem the function $f()$ is :

$$f(I,I')=Q\left(\frac{\eta T \frac{I+I'}{v\nu}}{\sqrt{m_1'I+m_0'I'+2\sigma_n^2}}\right)\frac{1}{2I} \frac{1}{\sqrt{2\pi\sigma_X^2}} \exp\left\lbrace \frac{[\ln(I)-\ln(I_0)]^2}{8\sigma_X^2} \right\rbrace\nonumber\\ \times\frac{1}{2I'} \frac{1}{\sqrt{2\pi\sigma_X^2}} \exp\left\lbrace \frac{[\ln(I')-\ln(I_0)]^2}{8\sigma_X^2} \right\rbrace$$ where $Q(.)$ is the integration of a normalized Gaussian distribution

The easy way out would be to do a two-step approximation: Write $$\int_{0}^{b} f(x,y)\,dx\approx\sum_{n=1}^{K}f(n\delta x,y)\delta_x$$ and estimate the error in this approximation, then integrate over $y$ to get $$\int_{0}^{a}\int_{0}^{b} f(x,y)\,dx\,dy\approx \int_{0}^{a}\sum_{n=1}^{K}f(n\delta x,y)\delta_x\,dy$$ noting that your error estimate from above now gets integrated from $0$ to $a$ as well; then estimate the integral on the right with yet a Riemann sum, once more applying the one-dimensional error estimates.