Q:

Build a sandpile

An abelian sandpile, for our purposes, is an infinite grid with integer coordinates, initially empty of sand. After each second, a grain of sand is placed at (0,0). Whenever a grid cell has 4 or more grains of sand, it spills one grain of sand to each of its four neighbors simultaneously. The neighbors of (x,y) are (x-1,y), (x+1,y), (x,y-1), and (x,y+1).
When a cell spills, it may cause its neighbors to spill. Some facts:

This cascade will eventually stop.
The order in which the cells spill is irrelevant; the result will be the same.

Example
After 3 seconds, the grid looks like
.....
.....
..3..
.....
.....

After 4 seconds:
.....
..1..
.1.1.
..1..
.....

After 15 seconds:
.....
..3..
.333.
..3..
.....

And after 16 seconds:
..1..
.212.
11.11
.212.
..1..

The challenge
In as few bytes as possible, write a function that takes a single positive
integer t and outputs a picture of the sandpile after t seconds.
Input
A single positive integer t, in any format you choose.
Output
A picture of the sandpile after t seconds, using the characters
 . 1 2 3

Edit: Use any four distinct characters you like, or draw a picture. If you're not using ".123" or "0123", specify in your answer what the characters signify.
Unlike in the examples, your output should contain the minimal number of rows and columns necessary to show the nonzero part of the sandpile.
That is, for input 3, the output should be
 3

For 4, the output should be
 .1.
 1.1
 .1.

Scoring
Standard golf scoring applies.
Rules
No language functions or libraries that already know what a sandpile is are allowed.
Edit: The output section has been edited, the character set restriction has been completely lifted. Use any four distinct characters or colors you like.

A:

R, 378 343 297 291 bytes
As usually, the user supplies his/her input via scan() (I already used the variable t, so let us take z instead), so the second line should be launched separately, and then the rest:
e=numeric
a=1%*%scan()
x=1
o=a>3
n=1
while(any(o)){
v=which(o,T)
if(any(v==1)){a=rbind(e(n+2),cbind(e(n),a,e(n)),e(n+2));x=x+1;n=n+2;v=which(a>3,T)}
q=nrow(v)
u=cbind(e(q),1)
l=v-u[,1:2];r=v+u[,1:2];t=v-u[,2:1];b=v+u[,2:1]
a[l]=a[l]+1;a[r]=a[r]+1;a[t]=a[t]+1;a[b]=a[b]+1
a[v]=a[v]-4
o=a>3}
a

Outputs an array that contains the values of a at tth generation (0, 1, 2 or 3).
Test cases:
z=3
     [,1]
[1,]    3
z=4
     [,1] [,2] [,3]
[1,]    0    1    0
[2,]    1    0    1
[3,]    0    1    0
z=16
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    1    0    0
[2,]    0    2    1    2    0
[3,]    1    1    0    1    1
[4,]    0    2    1    2    0
[5,]    0    0    1    0    0

It helps us that this thing is symmetrical both vertcially and horizontally, which means that is the leftmost point gets a height of 4, this means that the topmost, rightmost and lowest points are also 4.
Oh, and did I say that you can make beautiful visualisations?
After 1000 drops:

After 50000 drops (≈4 seconds):

After 333333 drops (≈15 minutes):

You can draw it, too!
image(1:n,1:n,a,col=colorRampPalette(c("#FFFFFF","#000000"))(4), axes=F, xlab="", ylab="")

This thing took 4 seconds for 10000 iterations but slows down considerably for larger array sizes (e.g. a couple of minutes for 100000 iterations). This is why it gets so slow (I estimated the growth rate as in  and obtained  τ(i)≈689·i^1.08, so the average time per one additional grain until the whole sandpile settles after ith step is slightly larger than one), and the total time as a function of the number of grains grows a little bit slower than quadratically (T(i)≈0.028*i^1.74):

And now with a complete explanation:
e=numeric # Convenient abbreviation for further repeated use
a=1%*%scan() # Creates a 1×1 array with a user-supplied number
x=1 # The coordinate of the centre
o=a>3 # Remember which cells were overflown
n=1 # Array height that is going to change over time
while(any(o)){ # If there is still any overflow
  v=which(o,T) # Get overflown cells' indices
  if(any(v==1)){ # If overflow occurred at the border, grow the array
    a=rbind(e(n+2),cbind(e(n),a,e(n)),e(n+2)) # Growing
    x=x+1 # Move the centre
    n=n+2 # Change the height
    v=which(a>3,T) # Re-index the overflowed cells
    }
  q=nrow(v) # See how many indices are overflown
  u=cbind(e(q),1) # Building block for neighbours' indices
  l=v-u[,1:2];r=v+u[,1:2];t=v-u[,2:1];b=v+u[,2:1] # L, R, T, B neighbours
  a[l]=a[l]+1;a[r]=a[r]+1;a[t]=a[t]+1;a[b]=a[b]+1 # Increment neighbours
  a[v]=a[v]-4 # Remove 4 grains from the overflown indices
  o=a>3} # See if still overflown indices remain
a # Output the matrix

This is the first time in my life when growing objects (like a <- c(a, 1)) works so much faster than pre-allocating a big empty matrix for values and filling it gradually with a ton of zeros unused.
Update. Golfed 18 bytes by removing arr.ind in which owing to Billywob and replacing rep(0,n) with e=numeric;e(n) in 5 instances owing to JDL, and 17 more bytes owing to JDL.
Update 2. Since the sandpile is Abelian, it may start with a stack of desired height, so I removed the redundant loop and gained a huge boost in productivity!

A:

MATL, 55 53 48 43 42 bytes
Inspired by @flawr's answer.
Graphical output:
0i:"Gto~+XytP*+t"t4=t1Y6Z+b+w~*]]tat3$)1YG

Try it at MATL Online!. It takes about 10 seconds for input 30. You may need to refresh the page and press "Run" again if it doesn't work.
Here's an example result for input 100:

ASCII output (43 bytes):
0i:"Gto~+XytP*+t"t4=t1Y6Z+b+w~*]]tat3$)48+c

Try it online!
Explanation
0          % Push a 0. This is the initial array. Will be resized in first iteration
i:         % Take input n. Generate range [1 2 ... n]
"          % For each, i.e. repeat n times
  Gto~+    %   Push input and add negate parity. This "rounds up" n to odd number
           %   m = n or n+1
  Xy       %   Identity matrix with that size
  tP*      %   Multiply element-wise by vertically flipped copy. This produces a
           %   matrix with a 1 in the center and the rest entries equal to 0
  +        %   Add to previous array. This updates the sandpile array
  t        %   Duplicate
  "        %   For each column (i.e. repeat m times)
    t4=    %     Duplicate. Compare with 4 element-wise. This gives a 2D mask that
           %     contains 1 for entries of the sandpile array that equal 4, and 0
           %     for the rest
    t      %     Duplicate
    1Y6    %     Predefined literal: [0 1 0; 1 0 1; 0 1 0]
    Z+     %     2D convolution, maintaining size
    b      %     Bubble up to bring sandpile array to top
    +      %     Element-wise addition. This adds 1 to the neighbours of a 4
    w      %     Swap to bring copy of mask to top
    ~*     %     Multiply bu negated mask. This removes all previous 4
  ]        %  End
]          % End
t          % Duplicate the updated sandpile array
a          % 1D mask that contains 1 for columns that contain a 1. This will be
           % used as a logical index to select columns
t          % Duplicate. This will be used as logical index to select rows (this
           % can be done because of symmetry)
3$)        % Keep only those rows and columns. This trims the outer zeros in the
           % sandpile array
1YG        % Display as scaled image

A:

Matlab, 160 156 148 bytes
n=input('');z=zeros(3*n);z(n+1,n+1)=n;for k=1:n;x=z>3;z=z+conv2(+x,1-mod(spiral(3),2),'s');z(x)=z(x)-4;end;v=find(sum(z));z=z(v,v);[z+48-(z<1)*2,'']

First a way too big matrix is created, with the n in the middle somewhere. Then the cascading is calculated with a very convenient 2d convolution. At the end the excess is trimmed off and the whole thing is converted to a string.
Example output for t=100
...121...
..32.23..
.3.323.3.
123.3.321
2.23.32.2
123.3.321
.3.323.3.
..32.23..
...121...

As always:

Convolution is the key to success.

