It turns out that $K=\mathbb{Q}[i,\varphi]$ is norm-Euclidean. The proof of this fact appears as Appendix A in https://arxiv.org/abs/2205.03007. I'll explain the heft of the argument here.
Let $R=\mathbb{Z}[i,\varphi]$ and $N=N_{K/\mathbb{Q}}$. We use this formulation of norm-Euclidean: for all $\alpha\in K$ there exists $\beta\in R$ with $N(\alpha+\beta)<1$.
Model $K\cong\mathbb{Q}^4$ via $w+x\varphi+yi+zi\varphi\mapsto(w,x,y,z)$ whence $R$ identifies with $\mathbb{Z}^4$. If $K$ were "extra nice" we could simply verify that for all $\alpha\in\left[-\frac{1}{2},\frac{1}{2}\right)^4\cap K=C$, we have $N(\alpha)<1$ whence for $\alpha\in K$ we could simply let $\beta=-[\alpha]$ where $[\alpha]$ is $\alpha$ with each coordinate rounded to the nearest integer. (This is how the standard proof goes for the Gaussian number field being norm-Euclidean.) Unfortunately, this approach fails, particularly near the corners where exactly three components have the same sign.
Notice, however, that $\left(\frac{1}{2},\frac{1}{2},\frac{1}{2},\frac{1}{2}\right)$ has very small norm, so that if we shift only some points in $C$ then we may be okay, since this would resolve e.g. $\left(\frac{1}{2},\frac{1}{2},\frac{1}{2},-\frac{1}{2}\right)$, and all points nearby, since $N$ is a polynomial in the coordinates, so it is definitely continuous. The trick then is to exactly quantify continuity. We will consider $C_n=\left[-\frac{1}{2},\frac{1}{2}\right)^4\cap(\frac{1}{n}\mathbb{Z}^4)$ with $n$ to be decided later. If we can show for each $\alpha\in C_n$ that there exists $\delta\in\mathbb{Z}^4$ such that $N(\alpha+\delta+B_n)<1$ where $B_n$ is the radius-$\frac{1}{n}$ $L^\infty$ ball in $\mathbb{Q}^4$, then we are done.
We consider now finding an explicit bound for $N(\alpha+\beta)$ where $\alpha=w+x\varphi+yi+zi\varphi\in K$ is any point and $\beta=d_1+d_2\varphi+d_3i+d_4i\varphi$ has $\|\beta\|_\infty<\varepsilon$, i.e. $\lvert d_1\rvert<\varepsilon$. $N(\alpha+\beta)$ is quartic in the eight variables, so write it out and group terms by "$d$-type," that is, $N(\alpha+\beta)=\sum_{e\in\mathbb{N}^4}d^eP_e(\alpha)$ for polynomials $P_e\in\mathbb{Q}[w,x,y,z]$. (e.g., $d_1d_2^2x+d_2d_3d_4z+2d_1d_2^2y+8d_4wyz$ would have $P_{1200}=x+2y$, $P_{0111}=z$, $P_{0001}=8wyz$, and $P_e=0$ otherwise.) Then, apply the triangle inequality (and positivity of $N$): $N(\alpha+\beta)=\left\lvert\sum_{e\in\mathbb{N}^4}d^eP_e(\alpha)\right\rvert\le\sum_{e\in\mathbb{N}^4}\lvert d\rvert^e\lvert P_e(\alpha)\rvert\le\sum_{e\in\mathbb{N}^4}\varepsilon^{\sum e}\lvert P_e(\alpha)\rvert$. (For the same example, the bound becomes $8\varepsilon\lvert wyz\rvert+\varepsilon^3(\lvert x+2y\rvert+\lvert z\rvert)$.) This now depends entirely on $\alpha$ and $\varepsilon=\frac{1}{2n}$. Let that bound be $f(\alpha,n)$ (also a polynomial in $\alpha$'s coordinates).
The strategy is then to check with a computer all $n^4$ elements $\alpha\in C_n$, and for each evaluate $f(\alpha,n)$. If $f(\alpha,n)<1$ then we move on, otherwise check $f(\alpha\pm e_j,n)$ for $e_j$ one of the four Kronecker basis vectors. It turns out that when $n=6$ we only ever need to check $\alpha$ and $\alpha\pm e_j$ (i.e. 0 or 1 steps) for this number field.
Curiously, it would seem that this approach could work to show that other biquadratic number fields $\mathbb{Q}\left[i,\sqrt{m}\right]$ of relatively small discriminant are norm-Euclidean, but this approach does not yield any results for $5<m,n<100$.
