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