I don't have anything rigorous to say, but let me share some images that may be useful or interesting to you. The Mathematica code to generate them is here. It's sparsely commented, so feel free to ask in the comments for clarification.
Below, the "out-component" is the set of vertices which are reachable by a directed path from the origin.
Here are a few example out-components at various $p$ in your model for a 161 by 161 grid:
I'm quite fond of this animated GIF file which shows the "averaged" out-component of the vertex at the origin in a 41 by 41 square grid as $p$ is tuned from 0 to 1 (in steps of 0.02). The intensity of a pixel corresponds to the frequency that that vertex was reachable from the origin in a set of 1000 pseudorandom configurations.
I'm not sure what to make of this pattern -- in particular, are they an artifact of the square boundary conditions, as they might cut off longer paths that would have made the dark regions parallel to the $x$ and $y$ axes reachable?
From the same data, here's the probability of percolation (existence of a directed path from the origin to the boundary of that 41 by 41 square grid) as a function of $p$:
And here's the mean fraction of the full grid that is reachable from the origin as a function of $p$:
Perhaps someone with more computer time can run do this with larger system sizes (my run took somewhere around an hour). I might do this for the last two graphs I showed, just to see how the transition sharpens for larger system sizes.
Edit. The last plot doesn't quite tell the full story about the distribution of out-component sizes.
Here's a plot showing the standard deviation of the out-component sizes:
Here's a sequence of plots showing histograms (from 1000 pseudorandomly generated configurations) of the fraction of total vertices reachable from the origin in a 41 by 41 grid at various $p$:
The distribution is bimodal sufficiently near $p=0.5$!
Here's a density plot of the fraction of vertices in the out-component as a function of $p$:
