We introduce a Partial Integral Equation (PIE) representation of Partial
Differential Equations (PDEs) in two spatial variables. PIEs are an algebraic
state-space representation of infinite-dimensional systems and have been used
to model 1D PDEs and time-delay systems without continuity constraints or
boundary conditions -- making these PIE representations amenable to stability
analysis using convex optimization. To extend the PIE framework to 2D PDEs, we
first construct an algebra of Partial Integral (PI) operators on the function
space L_2[x,y], providing formulae for composition, adjoint, and inversion. We
then extend this algebra to R^n x L_2[x] x L_2[y] x L_2[x,y] and demonstrate
that, for any suitable coupled, linear PDE in 2 spatial variables, there exists
an associated PIE whose solutions bijectively map to solutions of the original
PDE -- providing conversion formulae between these representations. Next, we
use positive matrices to parameterize the convex cone of 2D PI operators --
allowing us to optimize PI operators and solve Linear PI Inequality (LPI)
feasibility problems. Finally, we use the 2D LPI framework to provide
conditions for stability of 2D linear PDEs. We test these conditions on 2D heat
and wave equations and demonstrate that the stability condition has little to
no conservatism