$\newcommand\R{\mathbb R}$In general, the answer is no.
E.g., let $d=1$,
\begin{equation}
f_1(x):=\sum_{j=1}^\infty(|x|-2j)_+,\quad f_2(x):=\sum_{j=1}^\infty(|x|-2j-1)_+,
\end{equation}
where $u_+:=\max(0,u)$ --
so that $f_1$ and $f_2$ are convex functions, $f_1(x)\to\infty$ as $|x|\to\infty$, and $f_1\ge f_2$.
However, for any $b\in(0,1)$, for the function
$$h:=f_1-bf_2,$$
and for any natural $n$, we have
\begin{equation}
h(2n)-2h(2n+1)+h(2n+2)=-b<0.
\end{equation}
So, if a function $g_2$ is such that $g_2(x)-f_2(x)\to0$ as $|x|\to\infty$, then for $g:=f_1-bg_2$ we have $g(x)-h(x)\to0$ as $|x|\to\infty$, whence
\begin{equation}
g(2n)-2g(2n+1)+g(2n+2)<0
\end{equation}
for all large enough natural $n$.
So, $g$ is not convex for any $b\in(0,1)$.
In your concrete example, the answer is no as well. Indeed, let $\sigma_k=1$ for all $k$, and let $x:=(-1,u,0,\dots,0)\in\R^d$, $y:=(1,u,0,\dots,0)\in\R^d$, and $z:=(x+y)/2=(0,u,0,\dots,0)$, with $u\to\infty$. Then
\begin{equation}
f_1(x)-2f_1(z)+f_1(y)=2(1+u^2)-2u^2=2,
\end{equation}
whereas
\begin{equation}
f_2(x)-2f_2(z)+f_2(y)=2(1+u)^2-2u^2\to\infty.
\end{equation}
So, if $g_2(v)-f_2(v)\to0$ as $\|v\|\to\infty$, then for any $b\in(0,1)$ and for the function
$g:=f_1-bg_2$, we have $g(x)-2g(z)+g(y)\to-\infty$, so that $g$ is not convex for any $b\in(0,1)$.
In general, the answer is no. Moreover, the answer is no even if
\begin{equation}
\phi(t)=t\ln(1+t). \tag{1}
\end{equation}
Indeed, suppose that $P(Z_i=0)=1-2p$ and $P(Z_i=b)=p=P(Z_i=-b)$ for all $i$, where
\begin{equation*}
p:=\frac1{2\phi(b)},
\end{equation*}
$\phi$ is as given by (1),
and $b$ is a large enough positive real number so that $p\in(0,1/2)$.
Then for all $i$ we have $EZ_i=0$ and $E\phi(|Z_i|)=1$, so that $\|Z_i\|_\phi\le1$. On the other hand, for all real $c>0$ and all natural $n\ge2$
\begin{equation*}
\begin{aligned}
&E\phi\Big(\Big|\frac1n\sum_{i=1}^n Z_i\Big|/c\Big) \\
&\ge\sum_{i=1}^n \phi\Big(\frac b{cn}\Big)P(|Z_i|=b,\ Z_j=0\ \forall j\ne i) \\
&=n \phi\Big(\frac b{cn}\Big)2p(1-2p)^{n-1} \\
&=\frac{2pb}c\,\ln\Big(1+\frac b{cn}\Big)(1-2p)^{n-1}\to\frac1{2c}>1
\end{aligned}
\end{equation*}
as $n\to\infty$, if $b=n^2$ and $c\in(0,1/2)$. So, for all large enough $n$ we have $E\phi\big(\big|\frac1n\sum_{i=1}^n Z_i\big|/c\big)>1$ and hence $\|\frac1n\sum_{i=1}^n Z_i\|_\phi\ge c$ and hence
\begin{equation*}
\Big\|\frac1n\sum_{i=1}^n Z_i\Big\|_\phi\not\to0
\end{equation*}
as $n\to\infty$.
More generally, the answer will remain no if $\phi(t)=t \ell(t)$, where $\ell$ is any function such that $\ell(t)$ is slowly varying as $t\to\infty$. Yet more generally, the answer will remain no if $\phi(t)=t L(t)$, where $L$ is any function such that $\sup\limits_{K\in(0,\infty)}\limsup\limits_{t\to\infty}\dfrac{L(Kt)}{L(t)}<\infty$.
Best Answer
For each real $k>0$, \begin{equation}E\psi_\infty(|X|/k)=\infty\,P(|X|>k)+P(|X|=k) \\ =\left\{\begin{aligned}\infty\text{ if } P(|X|>k)>0,\\ P(|X|=k)\le1\text{ if } P(|X|>k)=0. \end{aligned}\right. \end{equation} So, indeed, $\|X\|_{\psi_\infty}<\infty$ iff $X$ is essentially bounded. Moreover, $\|X\|_{\psi_\infty}=\text{ess}\,\text{sup}\,|X|$.
Generally, for any non-constant nondecreasing convex function $F\colon[0,\infty)\to[0,\infty]$ such that $F(0)\le1$, the formula \begin{equation}\|X\|_F:=\inf\{t>0\colon EF(|X|/t)\le1\} \end{equation} defines a norm on the linear space, say $L_F$, of random variables (r.v.'s) $X$ on a probability space $\mathcal P$ with $\|X\|_F<\infty$. The proof of this is the same as the one in the case when $F$ is not allowed to take the value $\infty$. (If, in addition, it is assumed that $F(0+)<1$, then all bounded r.v.'s on $\mathcal P$ will be in $L_F$.)