Let $S\subseteq\mathbb R^3$ be some smooth surface (smooth 2-submanifold of $\mathbb R^3$) with global parameterization $\psi\colon\mathbb R^2\supseteq U\to S,(s,t)\mapsto\psi(s,t)$.
The surface area $\int_SdA$ of $S$ is given by
$$ \int_SdA % =\int_U\lvert\psi_s\times\psi_t\rvert ds\,dt,\qquad\psi_s:=\frac{\partial\psi}{\partial s},\quad\psi_t:=\frac{\partial\psi}{\partial t}. $$
Most resources I can find online justify this equation as follows, where I highlighted in $\color{red}{\text{red}}$ the one detail no explanation goes into more detail and which is the subject of my question. Keep in mind that the following explanation is supposed to be more of a detailed idea than a rigorous proof.
If we look at a sufficiently small surface patch $\psi(Q)\subseteq S$ with negligible curvature, where $Q=[s,s+\delta s]\times[t,t+\delta t]\subseteq U$ is an axis-aligned rectangle, then the surface area $\delta A$ of $\psi(Q)$ is $\color{red}{\text{approximately equal to}}$ the surface area $\delta\widetilde A$ of a parallelogram spanned by the vectors $\psi(s+\delta s,t)-\psi(s,t)$ and $\psi(s,t+\delta t)-\psi(s,t)$. According to Taylor's theorem these two vectors can be written as $\psi_s(s,t)\delta s+O(\delta s^2)$ and $\psi_t(s,t)\delta t+O(\delta t^2)$, respectively, as $\delta s,\delta t\to0$. Therefore, we have
$$ \begin{aligned} \delta A=\delta\widetilde A+(\delta A-\delta\widetilde A) % &=\bigl\lvert\bigl(\psi(s+\delta s,t)-\psi(s,t)\bigr)\times\bigl(\psi(s,t+\delta t)-\psi(s,t)\bigr)\bigr\rvert+(\delta A-\delta\widetilde A) % \\[2mm] % &=\bigl\lvert\bigl(\psi_s(s,t)\delta s+O(\delta s^2)\bigr)\times\bigl(\psi_t(s,t)\delta t+O(\delta t^2)\bigr)\bigr\rvert+(\delta A-\delta\widetilde A) % \\[2mm] % &=\lvert\psi_s(s,t)\times\psi_t(s,t)\rvert\delta s\,\delta t+O(\delta s\,\delta t^2)+O(\delta s^2\delta t)+(\delta A-\delta\widetilde A)\quad\text{as }\delta s,\delta t\to0, \end{aligned} $$ It follows for every integrable function $f\colon S\to\mathbb R$ that $$ \begin{aligned} \int_Sf\,dA &=\lim_{\delta s,\delta t\to0}\sum_{\substack{(s,t)\in U\cap\dots\\\dots(\delta s\mathbb Z\times\delta t\mathbb Z)}}f(\psi(s,t))\cdot\delta A_{(s,t)} \\[2mm] &=\lim_{\delta s,\delta t\to0}\Biggl[\sum_{(s,t)}f(\psi(s,t))\lvert\psi_s(s,t)\times\psi_t(s,t)\rvert\delta s\,\delta t+\underbrace{\sum_{(s,t)}O(\delta s\,\delta t^2)}_{=O(\delta t)~~(**)}+\underbrace{\sum_{(s,t)}O(\delta s^2\delta t)}_{=O(\delta s)~~(**)}+\sum_{(s,t)}\bigl(\delta A_{(s,t)}-\delta\widetilde A_{(s,t)}\bigr)\Biggr] \\[2mm] &=\int_Uf\lvert\psi_s\times\psi_t\rvert ds\,dt+\underbrace{\lim_{\delta s,\delta t\to0}\sum_{(s,t)}\bigl(\delta A_{(s,t)}-\delta\widetilde A_{(s,t)}\bigr)}_{\color{red}{\text{vanishes}}}, \end{aligned} $$ where $(**)$ follows from $$ \sum_{(s,t)}O(\delta s\,\delta t^2) =O\Biggl(\frac{\text{width of }U}{\delta s}\cdot\frac{\text{height of }U}{\delta t}\cdot\delta s\,\delta t^2\Biggr) =O(\delta t) $$ and, analogously, $\sum_{(s,t)}O(\delta s^2\delta t)=O(\delta s)$.
Now, to make this last term $\sum_{(s,t)}\bigl(\delta A_{(s,t)}-\delta\widetilde A_{(s,t)}\bigr)$ vanish as $\delta s,\delta t\to0$, we need $\delta A-\delta\widetilde A$ to go to zero strictly faster than $\delta s\,\delta t$, but is that actually the case?
If $P$ denotes the parallelogram mentioned before, then the surface patch $\psi(Q)$ is just a smooth, planar (curvature is negligible) deformation of $P$ keeping the three corners $\psi(s,t)$, $\psi(s+\delta s,t)$, and $\psi(s,t+\delta t)$ fixed in place.
The lengths of the lower edge $e:=\psi(s+\delta s,t)-\psi(s,t)$ of $P$ and the lower boundary $b:=\psi([s,s+\delta s]\times\{t\})$ of $\psi(Q)$ are bounded by
$$ \ell:=\delta s\underbrace{\sup_{\sigma\in[s,s+\delta s]}\lvert\psi_s(\sigma,t)\rvert}_{=O(1)\text{ as }\delta s,\delta t\to0}. $$
The area between $e$ and $b$ lies inside a rectangle with side lengths $\ell$, so the size of the area between $e$ and $b$ is bounded by $\ell^2=O(\delta s^2)$ as $\delta s\to0$. This means that the lower edge/boundary contributes a term $O(\delta s^2)$ to $\delta A-\delta\widetilde A$ which under summation $\sum_{(s,t)}$ becomes $O(\frac{\delta s}{\delta t})$. Clearly, this term does not vanish as $\delta s,\delta t\to0$. Luckily, this term is just an upper bound for the area between $e$ and $b$, so there is still room for improvement, I just don't see how.
Maybe someone else can enlighten me as to why $\delta A-\delta\widetilde A$ goes to zero faster than $\delta s\,\delta t$?