As advised by The Photon, I redid my calculation to compensate for the effect of an inhomogeneous microstrip surrounded by air. To do this, rather than using the \$ \epsilon_{r} \$ of the substrate, the effective relative permittivity \$ \epsilon_{r_{eff}} \$ should be used.
Rule-of-Thumb
First, I used a quick estimation I found online.
$$ \epsilon_{r_{eff}} = 0.64 \epsilon_{r} + 0.36 $$
It gives a resonant frequency of \$ f = \frac{c}{4l\sqrt{\epsilon_{reff}}} = 14.6 \text{GHz} \$, not too bad for a ballpark figure!
Hammerstad and Jensen
Next, I used the estimation by Hammerstad and Jensen, et al., found online from the textbook Microwave and RF design, a Systems Apporach, page 223 (Later I noticed Qucs also included this model, and actually the Schneider model is much easier to use than Hammerstad and Jensen, the formula is significantly shorter, but I already finished the original calculation at this point, sigh...)
Given \$ \epsilon_{r} \$, \$ w \$, and \$ h \$, Hammerstad and Jensen says the effective relative permittivity is:
$$ \epsilon_{r_{eff}} = \frac{\epsilon_r + 1}{2} + \frac{\epsilon_r - 1 }{2} (1 + \frac{10h}{w})^{-ab} $$
where
$$ u = \frac{w}{h} $$
$$ a(u) = 1 + \frac{1}{49} \ln[{\frac{u^4 + (u/52)^2}{u^4 + 0.432}}] + \frac{1}{18.7} \ln{[1+(\frac{u}{18.1})^3]} $$
$$ b(\epsilon_r) = 0.564 (\frac{\epsilon_r - 0.9}{\epsilon_r + 3})^{0.053} $$
Note that \$ u = \frac{w}{h} \$ is a dimensionless number, so absolute unit of measurement is irrelevant. The formula shows \$ \epsilon_{r_{eff}} = 3.05 \$ for my microstrip. Thus the resonant frequency \$ f = \frac{c}{4l\sqrt{\epsilon_{reff}}} = 14.3 \text{GHz} \$. It's exactly what the simulation shows.
Modeling for Fringing Effects
I also redid the simulation after connecting the open end of the "microstrip" to a "microstrip open" following The Photon's suggestions.
So after considering the fringing effects, a more realistic resonant frequency is 14.0 GHz.
Python code
For me to copy-paste in the future...
import math
def er_eff(er, w, h):
u = w / h
a = 1 + (1 / 49) * math.log((u ** 4 + (u / 52) ** 2) / (u ** 4 + 0.432)) + (1 / 18.7) * math.log(1 + (u / 18.1) ** 3)
b = 0.564 * ((er - 0.9) / (er + 3)) ** 0.053
return (er + 1) / 2 + ((er - 1) / 2) * (1 + (10 * h) / w) ** (-a * b)
>>> er_eff(4, 0.188, 0.1)
3.0544624555012434