Intuition
Yes it is possible. Both, the NIG and the VG process are exponential Lévy processes, i.e. they model the stock price via $S_t=S_0e^{X_t}$, where $X_t$ is a Lévy process. Here's a recent answer to the topic. Your question boils down to the following: if $X_t$ is a general Lévy process (VG, NIG, etc.), can we find parameters of $X_t$ such that $X_t$ collapses to a standard Brownian motion (with drift). In this case, $S_t$ is log-normal and we recover the Black-Scholes model.
A Lévy process has three components: a drift, a Brownian diffusion and a jump component. The simplest Lévy process is a Brownian motion itself. If you set the jump component equal to zero, scale the Brownian part by $\sigma>0$ and choose an appropriate drift (such that the discounted (reinvested) stock price is a $\mathbb{Q}$-martingale), you obtain a geometric Brownian motion. In this sense, all exponential Lévy processes generalise a geometric Brownian motion by including different jump components. Put differently, you can always find a parameter choice which recovers the Black-Scholes solution (for exponential Lévy processes).
An exponential Lévy processes is typically not described by an SDE but instead by stating its characteristic (Lévy) triplet, i.e. mean, volatility (covariance matrix for multidimensional processes) and jump measure. This makes working with them much easier. These components also give you immediately the characteristic function of the process. I thus show you below how to reason based on those Lévy components. A fantastic source on Lévy processes is the book from Cont and Tankov (2004).
Variance Gamma Process
There are different parametrisations for the VG process. I use the first notation from the original paper from Madan, Carr, and Chang (1998) with $\theta$, $\nu$ and $\sigma$. The variance gamma process ($X_t$) is a subordinated (i.e. time changed) Brownian motion. So let $$b(t;\theta,\sigma)=\theta t+\sigma B_t,$$ where $B_t$ is a standard Brownian motion. Let $\gamma(t;1,\nu)$ be a gamma process with unit mean. Then,
\begin{align*}
X(t;\theta,\sigma,\nu)&=b(\gamma(t;1,\nu);\theta,\sigma) \\
&= \theta\gamma(t;1,\nu)+\sigma B_{\gamma(t;1;\nu)}
\end{align*}
So, $\sigma$ is a standard volatility parameter, $\theta$ corresponds to the drift and $\nu$ is the variance rate of the gamma process. This parameter controls the jumps and the kurtosis, whereas $\theta$ controls the skewness. The case $\theta=0$ is known as symmetric VG process, see Madan and Seneta (1990). The characteristic function of $X_t=X(t;\theta,\sigma,\nu)$ is given by
\begin{align*}
\varphi_{X_t}(u) &= \left(\frac{1}{1-iu\theta\nu+\frac{1}{2}\sigma^2\nu u^2}\right)^{t/\nu}.
\end{align*}
Recall that $\nu$ governs the jumps of $X_t$. The larger $\nu$, the lower the exponential decay rate of the jump measure. Thus, jumps get more likely, which in turn increases the kurtosis (tails) of $X_t$. Conversely, for $\nu=0$, jumps are impossible, $\gamma$ is deterministic and we obtain a standard Brownian motion (with drift).
I quote from Madan, Carr, and Chang (1998):
There are three option pricing formulas nested in the option pricing formula
(25). These are a) the VG model, b) the symmetric VG (obtained by restricting
$\theta$ or $\alpha$ to zero) and c) the Black Scholes model (that results on setting $\nu$ equal to zero).
Theorem 2 in their paper derives an option pricing formula specific to the VG process. It looks very similar to a ``Black-Scholes type'' formula. See below my point on Fourier pricing. The more general CGMY model from Carr et al. (2002) also nests the VG model and thus, also the BS model (like all exponential Lévy processes).
Let $S_0=1$, $t=1$, $\theta=0$, $\sigma=0.2$. I plot what happens if $\nu\to0$. We would expect that $X_t\to\sigma B_t$ and thus, $S_t$ is log-normal distributed. In fact, the target characteristic function will be $\varphi_{\sigma B_t}(u)=e^{-\frac{1}{2}\sigma^2tu^2}$.
Here we go. I plotted the logs too because you couldn't see the difference between the BS characteristic function and the $\nu=0.01$ VG characteristic function otherwise. Even then, it's hard to make out the difference between the two.
Normal Inverse Gaussian
I'm briefer here because the argument is identical to VG: identify the jump parameter, set it equal to zero and we're finished. The stock price equals $S_t=S_0e^{X_t}$ and $X_t$ is a Lévy process obtained by subordination (time change). Here,
\begin{align*}
X_t=\mu t + \beta Z^{-1}(t;\delta,\gamma)+B_{Z^{-1}(t;\delta,\gamma)},
\end{align*}
where $Z^{-1}$ is a inverse Gaussian process. One typically introduces a new parameter $\alpha$ and sets $\gamma=\sqrt{\alpha^2+\beta^2}$. Then, $X_t$ is fully described by $\alpha,\beta,\delta,\mu$. The characteristic function is given by $$\varphi_{X_t}(u)=\exp\left(\mu tiu+\delta t\left(\sqrt{\alpha^2-\beta^2}-\sqrt{\alpha^2-(\beta-iu)^2}\right) \right).$$
I cite from Barndorff-Nielsen (1997):
We also note that the normal distribution $N(\mu, \sigma^2)$ appears as a limiting case for $\beta =0$, $a\to\infty$ and $\frac{\delta}{\alpha} =\sigma^2$.
Fourier Option Pricing
There are general option pricing formulae (Lewis (2001), Carr and Madan (1999) and others). They all apply to models with known characteristic function of $\ln(S_t)$. This applies in particular to exponential Lévy models and stochastic volatility models. For example, Bakshi and Madan's (2000) formula reads as
\begin{align*}
\mathrm{Call}(S_0,K,T) &= S_0 e^{-qT} I_1 - Ke^{-rT}I_2, \\
I_1 &= \frac{1}{2}+\frac{1}{\pi}\int_0^\infty \mathrm{Re}\left(\frac{e^{-iu\ln(K)}\varphi(u-i)}{iu\varphi(-i)}\right)\mathrm{d}u, \\
I_1 &= \frac{1}{2}+\frac{1}{\pi}\int_0^\infty \mathrm{Re}\left(\frac{e^{-iu\ln(K)}\varphi(u)}{iu}\right)\mathrm{d}u.
\end{align*}
Interestingly, this formula coincides with the delta-probability decomposition from Geman, El Karoui, and Rochet (1995). So you can interpret $I_1$ as exercise probability under the stock measure (which uses $S_te^{qt}$ as numéraire) and $I_2$ as probability of $\{S_T\geq K\}$ under the standard risk-neutral measure. As you see, these formulae are very general and hold for more general models than Black-Scholes, VG or NIG. To see formally, how they relate to the exercise probabilities, use the Fourier inversion formula from Gil-Pelaez (1951):
$$\mathbb{Q}[\{X_t\leq x\}] = \frac{1}{2}-\frac{1}{\pi}\int_0^\infty \mathrm{Re}\left(\frac{e^{-iux}\varphi_{X_t}^\mathbb{Q}(u)}{iu}\right)\mathrm{d}u,$$
where $\varphi_{X_t}^\mathbb{Q}$ is the characteristic function of $X_t$ under any probability measure $\mathbb{Q}$.
In any case, as the NIG and VG process collapse to a geometric Brownian motion, their characteristic functions coincide too. Thus, the above option pricing formula collapses to the standard Black-Scholes form, that is $I_1=N(d_1)$ and $I_2=N(d_2)$. So, the answer to question 3 is yes. You can take NIG or VG, choose the right parameters and obtain a geometric Brownian motion. Then, you go to the option pricing formula above and compute the probabilities of $\{S_T\geq K\}$ under different measures. Then, you recover the BS solution.
