These are just several thoughts. but it seems that they show in particular that the answr is $0$ for $k< p^p-1$.

$\def\FF{{\mathbb F}}$
**1.** Firstly, denoting $d=\mathop{\rm lcm}(p-1,p^2-1,\dots,p^p-1)$ we see that $A^{pd+p}=A^p$ for every matrix $A$ (the order of a semisimple component divides $d$, and we need $p$ in order for a nilpotent component to vanish). Thus we may assume that $k < pd+p$.

**2.** Your sum is equal to
$$
S=\sum_{a_{11}\in\FF_p} \dots\sum_{a_{pp}\in\FF_p}
\left(\sum_{i=1}^p\sum_{j=1}^p a_{ij}E_{ij}\right)^k,
$$
Now expand the inner brackets; we obtain
$$
S=\sum_J \left(\sum_{a_{11}\in\FF_p} \dots\sum_{a_{pp}\in\FF_p} a^J\right)M_J,
$$
where the outer summation is taken over some multiindices $J$ with $|J|=k$, and $M_J$ are some matrices. The summation in brackets over $a_{ij}$ gives zero unless the exponent of $a_{ij}$ is nonzero and divisible by $p-1$. Thus the sum in the brackets vanishes unless all the coordinates of $J$ are nonzero and divisible by $p-1$. So, if $k < p^2(p-1)$ then $S=0$, as is in the case $k=80$ and $p=17$.

**3.** For the remaining case, we need to calculate $M_J$. Assume that $J=(j_{11},\dots,j_{pp})$ (all $j$'s are multiples of $p-1$). Consider a digraph $G_J$ with $\FF_p$ as the set of vertices and $j_{k\ell}$ edges from $k$ to $\ell$. Now, if $M_J=[m_{k\ell}]$ then $m_{kk}$ is the number of Eulerian paths starting at $k$ (multiplies by the sum in brackets which is $-1$), and all other entries are zeroes.

Now let us show that the number of such cycles is divisible by $p$. We will assume that $k=1$, so that the cycles start and end at $1$. Split each cycle into subcycles starting at $1$, ending at $1$ and not passing through $1$ any more. Correspond to each cycle all the cycles obtained by permutatons of subcycles; thus the set of all cycles is partitioned into such equivalence classes, and the number of elements in each class is a corresponding multinomial coefficient $\binom{s}{s_1,\dots,s_t}$ where $s_1,\dots,s_t$ are the numbers of occurences of different subcycles. This coefficient is divisible by $p$ unless in the $p$-adic notation there are no transitions in the addition $s_1+\dots+s_t=k$.

Notice that there are at least $p$ distinct subcycles --- at least one starting from each of the edges $1\to 1$, $1\to 2$, \dots, $1\to p$. Moreover, we may partition all the subcycles into such classes --- the number in each will be divisible by $p-1$. THus, we have $p$ nonzero numbers divisble by $p-1$, and there shpuld be no transitions in addition of these $p$ numbers; this may happen only if $k\geq (p-1)+p(p-1)+\dots+p^{p-1}(p-1)=p^p-1$. So, for $k< p^p-1$ we definitely have $S=0$.

**EDIT.** Below in the comments, several improvements of this bound are shown.

`$\left( \begin{array}{cc} a & b \\ c & d \end{array}\right)^k $`

has the form`$\sum_{i=0}^k \sum_{j=0}^{\frac{i-k}{2}} \left(\begin{array}{c} i+j \\ j \end{array} \right) \left(\begin{array}{c} k-i-j-1 \\ j-1 \end{array} \right) a^i b^j c^j d^{k-i-2j}$`

If we sum over all matrices, the monomial terms become $0$ unless $i,k,k-i-2j$ are all positive multiples of $p-1$ and $1$ if they are. So we just get a combinatoric sum mod $p$. But the combinatorics seem really hard in general already for the $2 \times 2$ case. $\endgroup$