Simulation model for floating gate
One of the beautiful things about floating-gate transistors, if you use them in an analog sense, is that PVT doesn't matter as much as you can actually tweak with threshold by adding or removing charge on the floating node, \$V_Q\$.  The easiest way &quot;simulate&quot; a floating-gate would be to just put a voltage source on the gate to create an effective offset.
I can guess the those papers are by Duffy, Hasler, Basu, Tor, and Krumenacher; however, I do not believe that a hot-electron injection model due to impact ionization (which is why I assume you are using a pFET) has made it into the public, but I don't keep up with those circles anymore.  Chapter 2 of Hasler's Ph.D. thesis from Caltech under Carver Mead will give you the complete modeling for the nFET.  Duffy's draft is floating around and he did pFETs but to my knowledge, he did not complete his work. The quantum effects are the same for the tunneling for both nFETs and pFETs; however, the pFET physics for injection is different.

Here's a visual description of the process from my work.  You can modify Hasler's work by figuring out the probability of impact ionization and then the gate conditions required to attract the electron to the gate.  You can use a voltage-controlled current source with an ideal BJT to give you a nice logarithmic control between the gate and the drain because the barrier \$\Phi_{DC}\$ will control what you are looking to model.
I use the EKV model for this, but it is a bit sloppy just due to my implementation.  I actually back calculate what is required by the FETs data from EKV 2.6 extraction, and then get the doping and you can go from there.
EDIT:
Based on the comments, the current through a pFET in EKV without drain dependence would be
$$I_{f,r} =I_{thp}\ln^2 \left[1+ e^{\left[{\left(\kappa \left(V_b -V_g+V_{thp}\right)\right)- \left(V_{b}-V_{s,d}\right)}\right]/\left({2 U_{T}}\right)}  \right]$$
and this gives you an equation where the surface potential is \$\kappa V_g\$ because \$\kappa\$ is the channel divider.  This is what control the surface.  When you make the device &quot;float&quot;, you end up with this capacitively coupled mess:

The surface potential with respect to the new &quot;gate&quot; terminal is
$$V_{fg} = {V_{Q}} + \frac{C_{in}}{ C_T}V_g+ \frac{C_{tun}}{ C_T}V_{tun} 
    + \frac{C_{gd}}{C_T}V_{d} + \frac{C_{gs}}{C_T}V_{s} + \frac{C_{ox}}{C_T}V_b$$
Therefore, as the floating node becomes more &quot;negative&quot;, the threshold will shift from the standpoint of the gate input.  I always reference everything from the surface when I use these devices because you then don't need to worry about different capacitor sizes and behavior.