Let $f\in L^1(\mathbb{R}^d),$ let $\widetilde{M}f$ be the unrestricted maximal function $$\widetilde{Mf}(x_0, y_0) = \sup_{Q}\frac{1}{|Q|}\int_Q |f(x,y)|\,dx\,dy. $$ where the supremum is over all $Q=[x_0-k, x_0+k]\times[y-l, y+l]$ with $k, l>0$. Let me also introduce the functions $$M_1 f(x_0, y)=\sup_{k>0}\frac{1}{2k}\int_{x_0-k}^{x_0+k}|f(x,y)|\,dx,$$ and $$M_2 f(x, y_0)=\sup_{k>0}\frac{1}{2l}\int_{y_0-l}^{y_0+l}|f(x,y)|\,dx$$
I would like to show first that $\widetilde{M}f(x_0, y_0)\le M_1M_2 f(x_0, y_0)$ for all $(x_0, y_0)\in\mathbb{R}^2$.
After that, I would like to prove the existence of $C>0$ such that, for all $f\in L^2(\mathbb{R}^2)$, $$\|\widetilde{M}f\|_{L^2(\mathbb{R})}\le C \| f\|_{L^2(\mathbb{R})}.$$ Let me first note that we may write
\begin{aligned} \widetilde{Mf}(x_0, y_0)&=\sup_{k, l>0} \frac{1}{4kl} \int_{y_0-l}^{y_0+l} \int_{x_0-k}^{x_0+k} |f(x,y)|\,dy\,dx \\ &=\sup_{k,l>0} \frac{1}{2l} \int_{y_0-l}^{y_0+l} \frac{1}{2k}\int_{x_0-k}^{x_0+k} |f(x,y)|\,dy\,dx \\ &\le \sup_{k,l>0} \frac{1}{2l} \int_{y_0-l}^{y_0+l} \sup_{k>0} \int_{x_0-k}^{x_0+k} |f(x,y)|\,dy\,dx \\ &=\sup_{l>0} \frac{1}{2l}\int_{y_0-l}^{y_0+l} M_1f(x_0, y)\,dy \\ &= M_2 M_1 f(x_0, y_0). \end{aligned} The same argument works to prove the same result with $M_1 M_2$ replacing $M_2 M_1$. I believe this should prove the first point.
For the second point, I am assuming that I should use the first result to show it. With repeated Cauchy-Schwarz, I am able to get \begin{aligned} \widetilde{M}f &\le M_1 M_2 f \\ &\le \sup_{k,l>0} \frac{1}{\sqrt{4kl}}\left(\int_{Q}|f|^2\,dx\,dy\right)^{1/2} \\ &=\sqrt{\widetilde{M}(|f|^2)}, \end{aligned} but I don't believe this can show that $Mf\in L^2(\mathbb{R})$.