Ok, so I know that
$\ln(1+e^x)\approx x$ when $x$ is large.
But what about
$\ln(1+e^x+e^y)$ when both $x$ and $y$ are large?
I can figure out cases when $x\gg y$ or $y\gg x$ since that simplifies to the earlier case. But, not the case when both are large. Are there any approximations in this example?
One can certainly write for large $x, y$ that $$\ln(1 + e^x + e^y) \approx \ln(e^x + e^y) = \ln(e^x) \ln(1 + e^{y - x}) = x + \ln(1 + e^{y - x}),$$ which quickly leads to the given formulas when $x$ and $y$ are far apart.
In the limit where they are close, that is when $\epsilon := y - x$ is close to zero, we have that $$\ln(1 + e^{\epsilon}) = \ln \left[1 + \left(1 + \epsilon + \frac{1}{2} \epsilon^2 + O(\epsilon^3)\right)\right] = \ln 2 + \frac{1}{2} \epsilon + \frac{1}{8} \epsilon^2 + O(\epsilon^3),$$ and so for $x, y \gg 0$ we have $$\ln(1 + e^x + e^y) \approx \ln2 + \frac{1}{2} (x + y) + \frac{1}{8} (x - y)^2 + O((x - y)^3).$$ (It turns out that the part of the approximation that is polynomial in $(x - y)$ is even, and in particular we can replace $O((x - y)^3)$ with $O((x - y)^4)$.) In the limit $y \to x$ (i.e., $\epsilon \to 0$) this is just $\log 2 + x$, recovering vadim123's observation in his comment.
We will break this into 3 parts.
Part $1$ provides derivation of an series expansion for $\log(1+e^x+e^y)$ in terms of powers of $e^{-\left(x-y\right)}$, where it is assumed that $e^x>1+e^y$. This is equivalent to assuming that $x-y>\log\left(1+e^{-y}\right)$
In Part $2$, we develop a series converges rapidly when $x-y$ is adequately small. This development is motivated since the series in Part 1 converges slowly when $x-y$ is small.
Part $3$ provides upper and lower bounds that are global (i.e., We impose no conditions on $x$ or $y$.).
PART $1$:
Without loss of generality, assume $x>y$. Let $M=e^x$ and $m=e^y$ and further assume that $m+1<M$ to write
$$\begin{align} 1+e^x+e^y&=1+m+M\\\\ &=M\left(1+\frac{m+1}{M}\right) \end{align}$$
Then$$\log(1+e^x+e^y)=\log(M)+\log\left(1+\frac{m+1}{M}\right)$$
Using the Taylor series for $\log(1+x)$ yields
$$\log(1+e^x+e^y)=\log(M)+\sum_{n=1}^{\infty} (-1)^{n+1}\frac{1}{n}\left(\frac{m+1}{M}\right)^n$$
Next, using the binomial theorem to expand $\left(\frac{m+1}{M}\right)^n$ yields
$$\begin{align} \log(1+e^x+e^y)&=\log(M)+\sum_{n=1}^{\infty} \sum_{p=0}^{n}(-1)^{n+1}{n \choose p} \frac{1}{n}\left(\frac{1}{M}\right)^nm^{p}\\\\ &=x+\sum_{n=1}^{\infty} \sum_{p=0}^{n}(-1)^{n+1}{n \choose p} \frac{1}{n}e^{-\left(nx-py\right)} \end{align}$$
Examining the first few terms reveals
$$\begin{align} \log(1+e^x+e^y)&=x+e^{-(x-y)}-\frac12e^{-2(x-y)}+\frac13e^{-3(x-y)}+e^{-x}\left(1-e^{-(x-y)}+e^{-2(x-y)}\right)\\\\ &+O\left(e^{-2x}\right) \end{align}$$
This provides a good series expansion in powers of $e^{-(x-y)}$ when $e^{-(x-y)}<1-\frac{1}{1+e^y}<1$. However, the speed of series convergence decreases with decreasing $x-y>0$.
Part $2$ considers the case in which $x-y$ is "small."
PART $2$:
To address the case for which $0<x-y<1$, let's reexamine the term $\log\left(1+\frac{m+1}{M}\right)$. Set $s=-(x-y)$, and $w=M^{-1}$, and consider the function
$$f(s)=\log\left(1+w+e^s\right)$$
We proceed to expand $f$ in a Taylor series. Evaluating the derivatives up to third order at $s=0$ yields
$$\begin{align} f^{\left(0\right)}(s)&=\log(2+w)\\\\ f^{\left(1\right)}(s)&=\frac{1}{2+w}\\\\ f^{\left(2\right)}(s)&=\frac{1}{2+w}-\frac{1}{\left(2+w\right)^2}\\\\ f^{\left(3\right)}(s)&=\frac{1}{2+w}-\frac{3}{\left(2+w\right)^2}+\frac{2}{\left(2+w\right)^3} \end{align}$$
Next, note that the terms $(2+w)^{-k}=\frac{1}{2^k}\left(1+\frac{w}{2}\right)^k$ can be expanded as $$(2+w)^{-k}=\frac{1}{2^k}\left(1+O\left(w\right)\right)$$ where only the leading term was retained under the tacit assumption that $x>>1$ so that $w=e^{-x}<<1$. This permits simplifying the derivatives of $f$ as
$$\begin{align} f^{\left(0\right)}(0)&=\log(2)\\\\ f^{\left(1\right)}(0)&=\frac12\\\\ f^{\left(2\right)}(0)&=\frac14\\\\ f^{\left(3\right)}(0)&=0 \end{align}$$
Finally, we can express $\log(1+e^x+e^y)$ as
$$\begin{align} \log(1+e^x+e^y)&=x + f(0) + f^{(1)}(0) \left(-(x-y)\right) + \frac12 f^{(2)}(0) \left(-(x-y)\right) + O\left(\left(x-y\right)^4\right) + O\left(e^{-x}\right)\\\\ &=\log(2)+\frac12(x-y)+\frac18(x-y)^2+O\left(\left(x-y\right)^4\right)+O\left(e^{-x}\right) \end{align}$$
since $f^{\left(3\right)}(0)=0$. This validates the answer posted by @Travis for small $x-y$ and adds the last term to provide a strict equality.
PART $3$: We may also obtain crude bounds as follows. For the upper bound we have
$$\begin{align} 1+e^x+e^y&=(e^x+e^y)\left(1+\frac{1}{e^x+e^y}\right)\\\\ &<2\max(e^x,e^y)\left(1+\frac{1}{2\min(e^x,e^y)}\right)\\\\ \end{align}$$
Now use (i) the series for $\log(1+x)$ and (ii) the fact that the $\log$ is a monotonically increasing function, which implies $\log\left(\max(a,b)\right)=\max\left(\log(a),\log(b)\right)$ to reveal
$$\log\left(1+e^x+e^y\right)<\log(2)+\max(x,y)+\frac12\max\left(e^{-x},e^{-y}\right)$$
Under the assumption that $x>y$, this upper bound becomes
$$\log\left(1+e^x+e^y\right)\le\log(2)+x+\frac12e^{-y}$$
For the lower bound, we have
$$\begin{align} 1+e^x+e^y&=(e^x+e^y)\left(1+\frac{1}{e^x+e^y}\right)\\\\ &\ge2\sqrt(e^xe^y)\left(1+\frac{1}{e^x+e^y}\right) \end{align}$$
Now use (i) the series for $\log(1+x)$ and (ii) the inequality for the geometric and arithmetic means to reveal
$$\begin{align} \log(1+e^x+e^y)&\ge \log(2)+\frac12\left(x+y\right)+\frac{1}{e^x+e^y}-\frac12\left(\frac{1}{e^x+e^y}\right)^2\\\\ &\ge\log\left(2\right)+\frac12\left(x+y\right)+\frac12\min\left(e^{-x},e^{-y}\right)-\frac18\ \max\left(e^{-2x},e^{-2y}\right) \end{align}$$
Under the assumption that $x>y$, this lower bound becomes
$$\log(1+e^x+e^y)\ge\log\left(2\right)+\frac12\left(x+y\right)+\frac12e^{-x}-\frac18\ e^{-2y}$$
As @TravisWillse noted in his answer, for positive $\{x,\ y\}\gg 1$ you will have: $$\begin{array}{r c l} \ln(1+e^x+e^y)&\approx &\ln(e^x+e^y)\\ & = & x+\ln(1+e^{y-x}) \end{array}$$
from here you could do something different to better represent the function where $x\approx y$:
Working in other question I found this other approximation: $$\ln(1+e^x) \approx \begin{cases} \frac{x}{1-e^{-\frac{x}{\ln(2)}}},\quad x\neq 0\\ \ln(2),\quad x=0\end{cases}$$
For the last one the maximum error is below $0.0092$, and the plot is quite close:
Results that applying this to your problem leads to the following approximation for splitting the logarithm of the sum of exponentials LogSumExp(x,y): $$\ln(e^x+e^y) \approx \begin{cases} \displaystyle{\frac{xe^{\frac{x}{\ln(2)}}-ye^{\frac{y}{\ln(2)}}}{e^{\frac{x}{\ln(2)}}-e^{\frac{y}{\ln(2)}}}},\quad x\neq y\\ x+\ln(2),\quad x\equiv y\end{cases}$$
it works quite good if $e^x+e^y$ is not too close to zero.
