I need to calculate $f(k)=\frac{dg(k)}{dk}$, where
$g(k)=\iiint_{\{(x,y,z)\in\mathbb{R}^3:h(x,y,z)\leq k\}}e^{-(x+y+z)}dxdydz$
$h(x,y,z)\leq k$ is the domain of integration (and quite nasty, by the way). Obviously, integrating first and later differentiating does not work (too complicated). I have to solve it numerically and later look for a sound fit. If I get a fit for $g(k)$ and then differentiate to obtain an approximation for $f(k)$, it does not work very well.
Compute $g(k+\mathrm dk)-g(k)$ to the first order in $\mathrm dk$. This gives $$g(k+\mathrm dk)-g(k)=\iiint_{(x,y,z)\in\mathbb{R}^3,\; k\leq h(x,y,z)\leq k+\mathrm dk}\mathrm e^{-x-y-z}\;\mathrm dx\mathrm dy\mathrm dz$$ The domain $\{(x,y,z)\in\mathbb{R}^3,\; k\leq h(x,y,z)\leq k+\mathrm dk\}$ is a domain in space delimited by two surfaces $h(x,y,z)=k$ and $h(x,y,z)=k+\mathrm dk$. This volume is a "thickened" surface of thickness equal to $\epsilon(x,y,z)=\mathrm dk/\|\nabla h(x,y,z)\|$. You can rewrite $$g(k+\mathrm dk)-g(k)=\iiint_{(x,y,z)\in\mathbb{R}^3,\; h(x,y,z)=k} \epsilon(x,y,z)\mathrm e^{-x-y-z}\;\mathrm dx\mathrm dy\mathrm dz.$$ Keeping the first order in $\mathrm dk$ you get $$\frac{\mathrm dg}{\mathrm dk}=\iiint_{(x,y,z)\in\mathbb{R}^3,\; h(x,y,z)=k} \mathrm e^{-x-y-z}\frac{\mathrm dx\mathrm dy\mathrm dz}{\|\nabla h(x,y,z)\|}.$$
Beware that the integral is performed now on a surface...