Operator geostatistics from elliptic boundary value problems with Green operators interface transmission conditions and Schur complements

Many geostatistical problems live on bounded regions with complex boundaries and internal interfaces (faults, seals, stratigraphic contacts), where boundary and transmission conditions can materially alter correlation structure and uncertainty. We develop an operator-first formulation of Gaussian spatial random fields on bounded domains and manifolds with interfaces, in which the statistical model is specified by a coercive quadratic action (precision operator) together with boundary and transmission conditions. Variational theory yields a precise precision--covariance correspondence through Green operators, showing that variograms are derived quadratic functionals that depend on geometry and boundary/interface physics. Conditioning and kriging follow from exact Gaussian update identities expressed at the operator level, enabling consistent inference under mixed boundary conditions and across coupled subdomains. We provide worked micro-examples that (i) compute Dirichlet vs.\Neumann Green kernels and their induced variograms in closed form, and (ii) derive interface transmission conditions from surface energy penalties and quantify cross-interface covariance attenuation. A brief 2D synthetic experiment illustrates when boundary-condition misspecification matters most (interior-only data with long correlation range), and we give practical diagnostics and default recommendations for boundary-condition choice. The framework unifies domain reduction (Dirichlet-to-Neumann maps) and discrete Gaussian elimination (Schur complements) and aligns directly with scalable SPDE/GMRF discretizations through sparse operator representations. MSC 2020: 60G60, 62M30, 35J08, 35R60, 58J32.