The problem that I have is essentially an NxN grid with cells (i,j) and time period T. So there are 3 indices that I am concerned about(there are actually 4 considering the index k, but I think that 3 would be enough to get an understanding). As an example, the first equation written down would be required for each unique cell (i,j) across all time intervals. I need to somehow figure out a way to write this problem formulation in matrix notation so that when I solve the model using a commercial solver from OR that I will be able to pass to it the appropriate form.
However, I am thinking that this isn't really possible given the formulation that I have because of the exponentials and non-linearities. Do you think it is possible to write this in a square matrix notation?

You're right: the exponentials make the constraints nonlinear (and in fact nonpolynomial), which rules out matrix notation. The real issue isn't matrix notation, though: it's finding a solver that handle this type of nonlinear problem. That's a question that, as Rodrigo suggested, you might want to ask on OR Stack Exchange.