We have a piecewise continuous linear function (representing topography). The shape is to be filled with a constant volume of liquid (representing an ocean). How can we find the 'sea level', and where are the boundaries of this ocean?
As an example, here is a region defined by four points: Filling problem example.
The forward problem of finding the volume $V$ of the liquid given a known liquid level $l$ is a quadratic:
$$ V = -\frac{1}{2}(l-y_1)^2\frac{x_1-x_0}{y_1-y_0} + \frac{(2l-y_1-y_2)(x_2-x_1)}{2} + \frac{1}{2}(l-y_2)^2\frac{x_3-x_2}{y_3-y_2} = 1.625 $$
This quadratic can be easily simplified and solved for $l$ if $V$ is known, and for the example problem where $l=2.5$ we obtain $V = 1.625$.
However, this method requires prior knowledge of which points the boundary of the liquid will lie between, and so does not generalise. Nonlinear optimisation tools can provide a numerical answer, but I am ideally looking for solutions that might cast the problem in terms of a linear or quadratic programming problem, or some other efficient method of solution that might ideally be able to be generalised to the 2D case.
The good news is that the volume is a monotone function of the depth. So you can simply do binary search on the depth. And if you've worked out that $$ V(a) < T $$ and $$ V(b) > T $$ where $T$ is the target volume, and $V(d)$ denotes the volume of water when the depth is $d$, then in computing $V(\frac{a+b}{2})$, it's possible that you can simply compute the volume between depth $b$ and depth $\frac{a+b}{2}$, re-using much of your prior work. But that'd take some effort to get right.
You could also do a little preprocessing, looking at each "cell" (presumably something like a quad in $xy$, divided by a well-chosen diagonal), and identifying the minimum $z$ coordinate (for depths below that, this cell can be ignored), the max $z$ coordinate (for depth $d$ above this level $s$, the cell contributes $\Delta d \cdot A$ when the depth increases by $\Delta d$, where $A$ is the area of each $xy$-rectangle), and the volume within the cell when the depth is exactly $s$. Those numbers let you rapidly estimate an upper and lower bound on the volume at some depth, so that your search-region gets small.
Anyhow, assuming an $n \times n$ grid of data, there are $n^2$ depth values, and the right answer lies between some two of these. Identifying which two involves $\log (n^2) = 2 \log n$ evaluations of volume, each of which takes time proportional to $n^2$ unless you do something clever with clustering your depths into larger blocks in some spatial data structure, so this gives about $n^2 \log n$ runtime. That's not too bad, when you realize that $n^2$ is the size of the input data rather than $n$.
By the way, you should expect this problem to be hard. Imagine a (1D) region of rough bottom topography, more or less sawtooth, and a target depth that's about half way up the "tooth". In this situation, every single cell in the domain is involved in the depth computation at the critical depth.
If, on the other hand, your ocean bottom has a gradual rise towards the shore, then at the critical level, exactly one cell will be involved in the computation, so some instances of the problem should be reasonably straightforward.
Wait...it's easier than that. Still not pretty, but...that's OK.
Let's assume the grid-squares are unit squares, i.e., the domain grid is something like $\{1, \ldots, n\} \times \{1, \ldots, n \}.$
Let's look at one "cell", the rectangular column of liquid above a unit square. To each corner of this cell is associated a height, which I'll call a $z$-value. Let's call them $z_{00}, z_{01}, z_{10}, z_{11}$. You say that the depth is piecewise linear, so lets make a triangulation: in the domain square $0 \le x,y \le 1$, we draw a diagonal from $(0,0)$ to $(1,1)$. Now instead of a rectangular column of liquid, we have two triangular columns, each with base-area $\frac12$. Let's take one, with heights I'll call $z_A, z_B, z_C$; to make things simple, let's assume $z_A \le z_B \le z_C.$ If we plot the volume of water in this triangular cylinder when the water-depth is $d$, and vary $d$, we'll get an increasing function $V(d)$.
For $d < \min(z_A, z_B, z_C)$, we'll have $V(d) = 0.$ For $z_A \le z \le z_B$, we'll have $V(d)$ is some cubic function of $d$, which can be expressed in terms of $z_A, z_B, z_C$ by doing some geometry which I'm not prepared to write out, because...it's just basic geometry. For $z_B \le z \le z_C$, we have a similar result. But for $z > z_C$, $V(d)$ has the particularly nice form $$ V(d) = K + \frac{1}{2} (z - z_C). $$
Each of the four "pieces" of this function $V$ is a cubic whose coefficients depend only on $v_A, v_B, v_C$. So you can represent it by four numbers (the coefficients $p, q, r, s$ of $d^3, d^2, d,$ and $1$).
So let's represent each of our $2n^2$ "half-cells" (triangular prisms) by a list of four quintuples, each of which looks like this: $$ (e; p, q, r, s) $$ indicating that starting at depth $e$, the next "segment" of the function is represented by the quadratic $d \mapsto pd^2 + qd + r$. So if the depth $z_A, z_B, z_C$ are $1, 2, 5.3$, then our four quadruples would look look $$ (0; 0,0,0, 0)\\ (1; *, *, *, *) \\ (2; *, *, *, *) \\ (5.3; 0, 0, \frac12, *) $$
Now suppose that we want to compute volume as a function of $d$. We have $n^2$ depth-values at the grid cells, and between any two of these, the depth is a sum of $n^2$ cubics, one contributed by each grid cell. But a sum of cubics is still a cubs, so from the $10n^2$ data items you've got, you can build a table, sorted by depth, that looks just like the list of 4 quintuples above: each item is a depth, followed by the coefficients of the cubic volume-function starting at that depth and in use until the next depth.
It might appear that at each of the $n^2$ depth values, you'd need to sum up all the $n^2$ cubics, but ... assuming all depths are distinct for the moment ... at each depth value, at most four cells change their cubics, so you need only deal with a constant amount of work.
In short, in $n^2$ time you can build a table describing the entire piecewise-cubic depth function. And then you can do binary search to find the two $z$-depths between which your goal-volume occurs, and then you just solve the cubic for the depth. (You could use Cardano's formula, but numerical methods will work just fine here.)
So...building the table takes time proportional to your input size $n^2$. Binary search on that takes time $\log (4n^2) \approx \log n$. And the you solve one cubic. All in all: linear in the input size, which is pretty good!