In FEM, why is the stiffness matrix positive definite?

The property follows from the property of the corresponding (weak form of the) partial differential equation; this is one of the advantages of finite element methods compared to, e.g., finite difference methods.

To see that, first recall that the finite element method starts from the weak form of the Poisson equation (I'm assuming Dirichlet boundary conditions here): Find $uin H^1_0(Omega)$ such that $$ a(u,v):= int_Omega nabla ucdot nabla v ,dx = int_Omega fv,dx qquadtext{for all }vin H^1_0(Omega).$$ The important property here is that $$ a(v,v) = |nabla v|_{L^2}^2 geq c |v|_{H^1}^2 qquadtext{for all }vin H^1_0(Omega). tag{1}$$ (This follows from Poincaré's inequality.)

Now the classical finite element approach is to replace the infinite-dimensional space $H^1_0(Omega)$ by a finite-dimensional subspace $V_hsubset H^1_0(Omega)$ and find $u_hin V_h$ such that $$ a(u_h,v_h):= int_Omega nabla u_hcdot nabla v_h ,dx = int_Omega fv_h,dx qquadtext{for all }v_hin V_h.tag{2}$$ The important property here is that you are using the same $a$ and a subspace $V_hsubset H^1_0(Omega)$ (a conforming discretization); that means that you still have $$ a(v_h,v_h) geq c |v_h|_{H^1}^2 >0 qquadtext{for all }v_hin V_h. tag{3}$$

Now for the last step: To transform the variational form to a system of linear equations, you pick a basis ${varphi_1,dots,varphi_N}$ of $V_h$, write $u_h =sum_{i=1}^N u_ivarphi_i$ and insert $v_h=varphi_j$, $1leq jleq N$ into $(2)$. The stiffness matrix $K$ then has the entries $K_{ij}=a(varphi_i,varphi_j)$ (which coincides with what you wrote).

Now take an arbitrary vector $vec v=(v_1,dots,v_N)^Tin mathbb{R}^N$ and set $v_h:=sum_{i=1}^Nv_i varphi_iin V_h$. Then we have by $(3)$ and the bilinearity of $a$ (i.e., you can move scalars and sums into both arguments) $$ vec v^T K vec v = sum_{i=1}^Nsum_{j=1}^N v_iK_{ij} v_j = sum_{i=1}^Nsum_{j=1}^N a(v_ivarphi_i,v_jvarphi_j) = a(v_h,v_h) >0.$$ Since $vec v$ was arbitrary, this implies that $K$ is positive definite.

TL;DR: The stiffness matrix is positive definite because it comes from a conforming discretization of a (self-adjoint) elliptic partial differential equation.

If stiffness of element is not positive, then the system is not stable. So the model is most likely not correct. Look at the most basic equation of harmonic oscillator

$$m x''(t) + k x(t) = f(t)$$

The solution is unstable if $k$ is negative (look at the roots of the characteristic equation). It means the solution will blow up. The stiffness has to be a restoring force. At least for a physical spring. The stiffness matrix extends this to large number of elements (global stiffness matrix). That is all. But it is the same basic idea. FEM basis is in the stiffness matrix method for structural analysis where each element has a stiffness associated with it.