This paper presents a Partial Integral Equation (PIE) representation of linear 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 L2[x, y], providing formulae for composition, adjoint, and inversion. We then extend this algebra to Rn × L2[x] × L2[y] × L2[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.